diff --git a/doc/src/modules/solvers/ode.rst b/doc/src/modules/solvers/ode.rst index 8f9c74fe856c..1c5e07e76944 100644 --- a/doc/src/modules/solvers/ode.rst +++ b/doc/src/modules/solvers/ode.rst @@ -102,12 +102,12 @@ factorable 2nd_linear_airy ^^^^^^^^^^^^^^^ -.. autoclass:: sympy.solvers.ode.single::SecondLinearAiry +.. autoclass:: sympy.solvers.ode.single::LinearAiry2nd :members: 2nd_linear_bessel ^^^^^^^^^^^^^^^^^ -.. autoclass:: sympy.solvers.ode.single::SecondLinearBessel +.. autoclass:: sympy.solvers.ode.single::LinearBessel2nd :members: Bernoulli @@ -187,11 +187,12 @@ separable_reduced lie_group ^^^^^^^^^ -.. autofunction:: sympy.solvers.ode.ode::ode_lie_group +.. autoclass:: sympy.solvers.ode.ode::LieGroup + :members: 2nd_hypergeometric ^^^^^^^^^^^^^^^^^^ -.. autoclass:: sympy.solvers.ode.single::SecondHypergeometric +.. autoclass:: sympy.solvers.ode.single::Hypergeometric2nd :members: 1st_power_series @@ -214,39 +215,39 @@ implemented for the various heuristics. abaco1_simple ^^^^^^^^^^^^^ -.. autofunction:: sympy.solvers.ode.ode::lie_heuristic_abaco1_simple +.. autofunction:: sympy.solvers.ode.lie_group::lie_heuristic_abaco1_simple abaco1_product ^^^^^^^^^^^^^^ -.. autofunction:: sympy.solvers.ode.ode::lie_heuristic_abaco1_product +.. autofunction:: sympy.solvers.ode.lie_group::lie_heuristic_abaco1_product bivariate ^^^^^^^^^ -.. autofunction:: sympy.solvers.ode.ode::lie_heuristic_bivariate +.. autofunction:: sympy.solvers.ode.lie_group::lie_heuristic_bivariate chi ^^^ -.. autofunction:: sympy.solvers.ode.ode::lie_heuristic_chi +.. autofunction:: sympy.solvers.ode.lie_group::lie_heuristic_chi abaco2_similar ^^^^^^^^^^^^^^ -.. autofunction:: sympy.solvers.ode.ode::lie_heuristic_abaco2_similar +.. autofunction:: sympy.solvers.ode.lie_group::lie_heuristic_abaco2_similar function_sum ^^^^^^^^^^^^ -.. autofunction:: sympy.solvers.ode.ode::lie_heuristic_function_sum +.. autofunction:: sympy.solvers.ode.lie_group::lie_heuristic_function_sum abaco2_unique_unknown ^^^^^^^^^^^^^^^^^^^^^ -.. autofunction:: sympy.solvers.ode.ode::lie_heuristic_abaco2_unique_unknown +.. autofunction:: sympy.solvers.ode.lie_group::lie_heuristic_abaco2_unique_unknown abaco2_unique_general ^^^^^^^^^^^^^^^^^^^^^ -.. autofunction:: sympy.solvers.ode.ode::lie_heuristic_abaco2_unique_general +.. autofunction:: sympy.solvers.ode.lie_group::lie_heuristic_abaco2_unique_general linear ^^^^^^ -.. autofunction:: sympy.solvers.ode.ode::lie_heuristic_linear +.. autofunction:: sympy.solvers.ode.lie_group::lie_heuristic_linear System of ODEs -------------- diff --git a/sympy/solvers/ode/__init__.py b/sympy/solvers/ode/__init__.py index 127ac789146d..2b543425251d 100644 --- a/sympy/solvers/ode/__init__.py +++ b/sympy/solvers/ode/__init__.py @@ -1,5 +1,7 @@ from .ode import (allhints, checkinfsol, classify_ode, - constantsimp, dsolve, homogeneous_order, infinitesimals) + constantsimp, dsolve, homogeneous_order) + +from .lie_group import infinitesimals from .subscheck import checkodesol diff --git a/sympy/solvers/ode/lie_group.py b/sympy/solvers/ode/lie_group.py new file mode 100644 index 000000000000..57d1d4772523 --- /dev/null +++ b/sympy/solvers/ode/lie_group.py @@ -0,0 +1,1082 @@ +from itertools import islice + +from sympy.core import Add, S, Mul, Pow +from sympy.core.exprtools import factor_terms +from sympy.core.function import Function, AppliedUndef, expand +from sympy.core.relational import Equality, Eq +from sympy.core.symbol import Symbol, Wild, Dummy, symbols +from sympy.functions import exp, log +from sympy.integrals.integrals import integrate +from sympy.polys import Poly +from sympy.polys.polytools import cancel, div +from sympy.simplify import (collect, powsimp, # type: ignore + separatevars, simplify) +from sympy.solvers import solve +from sympy.solvers.pde import pdsolve + +from sympy.utilities import numbered_symbols +from sympy.solvers.deutils import _preprocess, ode_order +from .ode import checkinfsol + + +lie_heuristics = ( + "abaco1_simple", + "abaco1_product", + "abaco2_similar", + "abaco2_unique_unknown", + "abaco2_unique_general", + "linear", + "function_sum", + "bivariate", + "chi" + ) + + +def _ode_lie_group_try_heuristic(eq, heuristic, func, match, inf): + + xi = Function("xi") + eta = Function("eta") + f = func.func + x = func.args[0] + y = match['y'] + h = match['h'] + tempsol = [] + if not inf: + try: + inf = infinitesimals(eq, hint=heuristic, func=func, order=1, match=match) + except ValueError: + return None + for infsim in inf: + xiinf = (infsim[xi(x, func)]).subs(func, y) + etainf = (infsim[eta(x, func)]).subs(func, y) + # This condition creates recursion while using pdsolve. + # Since the first step while solving a PDE of form + # a*(f(x, y).diff(x)) + b*(f(x, y).diff(y)) + c = 0 + # is to solve the ODE dy/dx = b/a + if simplify(etainf/xiinf) == h: + continue + rpde = f(x, y).diff(x)*xiinf + f(x, y).diff(y)*etainf + r = pdsolve(rpde, func=f(x, y)).rhs + s = pdsolve(rpde - 1, func=f(x, y)).rhs + newcoord = [_lie_group_remove(coord) for coord in [r, s]] + r = Dummy("r") + s = Dummy("s") + C1 = Symbol("C1") + rcoord = newcoord[0] + scoord = newcoord[-1] + try: + sol = solve([r - rcoord, s - scoord], x, y, dict=True) + if sol == []: + continue + except NotImplementedError: + continue + else: + sol = sol[0] + xsub = sol[x] + ysub = sol[y] + num = simplify(scoord.diff(x) + scoord.diff(y)*h) + denom = simplify(rcoord.diff(x) + rcoord.diff(y)*h) + if num and denom: + diffeq = simplify((num/denom).subs([(x, xsub), (y, ysub)])) + sep = separatevars(diffeq, symbols=[r, s], dict=True) + if sep: + # Trying to separate, r and s coordinates + deq = integrate((1/sep[s]), s) + C1 - integrate(sep['coeff']*sep[r], r) + # Substituting and reverting back to original coordinates + deq = deq.subs([(r, rcoord), (s, scoord)]) + try: + sdeq = solve(deq, y) + except NotImplementedError: + tempsol.append(deq) + else: + return [Eq(f(x), sol) for sol in sdeq] + + + elif denom: # (ds/dr) is zero which means s is constant + return [Eq(f(x), solve(scoord - C1, y)[0])] + + elif num: # (dr/ds) is zero which means r is constant + return [Eq(f(x), solve(rcoord - C1, y)[0])] + + # If nothing works, return solution as it is, without solving for y + if tempsol: + return [Eq(sol.subs(y, f(x)), 0) for sol in tempsol] + return None + + +def _ode_lie_group( s, func, order, match): + + heuristics = lie_heuristics + inf = {} + f = func.func + x = func.args[0] + df = func.diff(x) + xi = Function("xi") + eta = Function("eta") + xis = match['xi'] + etas = match['eta'] + y = match.pop('y', None) + if y: + h = -simplify(match[match['d']]/match[match['e']]) + y = y + else: + y = Dummy("y") + h = s.subs(func, y) + + if xis is not None and etas is not None: + inf = [{xi(x, f(x)): S(xis), eta(x, f(x)): S(etas)}] + + if checkinfsol(Eq(df, s), inf, func=f(x), order=1)[0][0]: + heuristics = ["user_defined"] + list(heuristics) + + match = {'h': h, 'y': y} + + # This is done so that if any heuristic raises a ValueError + # another heuristic can be used. + sol = None + for heuristic in heuristics: + sol = _ode_lie_group_try_heuristic(Eq(df, s), heuristic, func, match, inf) + if sol: + return sol + return sol + + +def infinitesimals(eq, func=None, order=None, hint='default', match=None): + r""" + The infinitesimal functions of an ordinary differential equation, `\xi(x,y)` + and `\eta(x,y)`, are the infinitesimals of the Lie group of point transformations + for which the differential equation is invariant. So, the ODE `y'=f(x,y)` + would admit a Lie group `x^*=X(x,y;\varepsilon)=x+\varepsilon\xi(x,y)`, + `y^*=Y(x,y;\varepsilon)=y+\varepsilon\eta(x,y)` such that `(y^*)'=f(x^*, y^*)`. + A change of coordinates, to `r(x,y)` and `s(x,y)`, can be performed so this Lie group + becomes the translation group, `r^*=r` and `s^*=s+\varepsilon`. + They are tangents to the coordinate curves of the new system. + + Consider the transformation `(x, y) \to (X, Y)` such that the + differential equation remains invariant. `\xi` and `\eta` are the tangents to + the transformed coordinates `X` and `Y`, at `\varepsilon=0`. + + .. math:: \left(\frac{\partial X(x,y;\varepsilon)}{\partial\varepsilon + }\right)|_{\varepsilon=0} = \xi, + \left(\frac{\partial Y(x,y;\varepsilon)}{\partial\varepsilon + }\right)|_{\varepsilon=0} = \eta, + + The infinitesimals can be found by solving the following PDE: + + >>> from sympy import Function, Eq, pprint + >>> from sympy.abc import x, y + >>> xi, eta, h = map(Function, ['xi', 'eta', 'h']) + >>> h = h(x, y) # dy/dx = h + >>> eta = eta(x, y) + >>> xi = xi(x, y) + >>> genform = Eq(eta.diff(x) + (eta.diff(y) - xi.diff(x))*h + ... - (xi.diff(y))*h**2 - xi*(h.diff(x)) - eta*(h.diff(y)), 0) + >>> pprint(genform) + /d d \ d 2 d + |--(eta(x, y)) - --(xi(x, y))|*h(x, y) - eta(x, y)*--(h(x, y)) - h (x, y)*--(x + \dy dx / dy dy + + d d + i(x, y)) - xi(x, y)*--(h(x, y)) + --(eta(x, y)) = 0 + dx dx + + Solving the above mentioned PDE is not trivial, and can be solved only by + making intelligent assumptions for `\xi` and `\eta` (heuristics). Once an + infinitesimal is found, the attempt to find more heuristics stops. This is done to + optimise the speed of solving the differential equation. If a list of all the + infinitesimals is needed, ``hint`` should be flagged as ``all``, which gives + the complete list of infinitesimals. If the infinitesimals for a particular + heuristic needs to be found, it can be passed as a flag to ``hint``. + + Examples + ======== + + >>> from sympy import Function + >>> from sympy.solvers.ode.lie_group import infinitesimals + >>> from sympy.abc import x + >>> f = Function('f') + >>> eq = f(x).diff(x) - x**2*f(x) + >>> infinitesimals(eq) + [{eta(x, f(x)): exp(x**3/3), xi(x, f(x)): 0}] + + References + ========== + + - Solving differential equations by Symmetry Groups, + John Starrett, pp. 1 - pp. 14 + + """ + + if isinstance(eq, Equality): + eq = eq.lhs - eq.rhs + if not func: + eq, func = _preprocess(eq) + variables = func.args + if len(variables) != 1: + raise ValueError("ODE's have only one independent variable") + else: + x = variables[0] + if not order: + order = ode_order(eq, func) + if order != 1: + raise NotImplementedError("Infinitesimals for only " + "first order ODE's have been implemented") + else: + df = func.diff(x) + # Matching differential equation of the form a*df + b + a = Wild('a', exclude = [df]) + b = Wild('b', exclude = [df]) + if match: # Used by lie_group hint + h = match['h'] + y = match['y'] + else: + match = collect(expand(eq), df).match(a*df + b) + if match: + h = -simplify(match[b]/match[a]) + else: + try: + sol = solve(eq, df) + except NotImplementedError: + raise NotImplementedError("Infinitesimals for the " + "first order ODE could not be found") + else: + h = sol[0] # Find infinitesimals for one solution + y = Dummy("y") + h = h.subs(func, y) + + u = Dummy("u") + hx = h.diff(x) + hy = h.diff(y) + hinv = ((1/h).subs([(x, u), (y, x)])).subs(u, y) # Inverse ODE + match = {'h': h, 'func': func, 'hx': hx, 'hy': hy, 'y': y, 'hinv': hinv} + if hint == 'all': + xieta = [] + for heuristic in lie_heuristics: + function = globals()['lie_heuristic_' + heuristic] + inflist = function(match, comp=True) + if inflist: + xieta.extend([inf for inf in inflist if inf not in xieta]) + if xieta: + return xieta + else: + raise NotImplementedError("Infinitesimals could not be found for " + "the given ODE") + + elif hint == 'default': + for heuristic in lie_heuristics: + function = globals()['lie_heuristic_' + heuristic] + xieta = function(match, comp=False) + if xieta: + return xieta + + raise NotImplementedError("Infinitesimals could not be found for" + " the given ODE") + + elif hint not in lie_heuristics: + raise ValueError("Heuristic not recognized: " + hint) + + else: + function = globals()['lie_heuristic_' + hint] + xieta = function(match, comp=True) + if xieta: + return xieta + else: + raise ValueError("Infinitesimals could not be found using the" + " given heuristic") + + +def lie_heuristic_abaco1_simple(match, comp=False): + r""" + The first heuristic uses the following four sets of + assumptions on `\xi` and `\eta` + + .. math:: \xi = 0, \eta = f(x) + + .. math:: \xi = 0, \eta = f(y) + + .. math:: \xi = f(x), \eta = 0 + + .. math:: \xi = f(y), \eta = 0 + + The success of this heuristic is determined by algebraic factorisation. + For the first assumption `\xi = 0` and `\eta` to be a function of `x`, the PDE + + .. math:: \frac{\partial \eta}{\partial x} + (\frac{\partial \eta}{\partial y} + - \frac{\partial \xi}{\partial x})*h + - \frac{\partial \xi}{\partial y}*h^{2} + - \xi*\frac{\partial h}{\partial x} - \eta*\frac{\partial h}{\partial y} = 0 + + reduces to `f'(x) - f\frac{\partial h}{\partial y} = 0` + If `\frac{\partial h}{\partial y}` is a function of `x`, then this can usually + be integrated easily. A similar idea is applied to the other 3 assumptions as well. + + + References + ========== + + - E.S Cheb-Terrab, L.G.S Duarte and L.A,C.P da Mota, Computer Algebra + Solving of First Order ODEs Using Symmetry Methods, pp. 8 + + + """ + + xieta = [] + y = match['y'] + h = match['h'] + func = match['func'] + x = func.args[0] + hx = match['hx'] + hy = match['hy'] + xi = Function('xi')(x, func) + eta = Function('eta')(x, func) + + hysym = hy.free_symbols + if y not in hysym: + try: + fx = exp(integrate(hy, x)) + except NotImplementedError: + pass + else: + inf = {xi: S.Zero, eta: fx} + if not comp: + return [inf] + if comp and inf not in xieta: + xieta.append(inf) + + factor = hy/h + facsym = factor.free_symbols + if x not in facsym: + try: + fy = exp(integrate(factor, y)) + except NotImplementedError: + pass + else: + inf = {xi: S.Zero, eta: fy.subs(y, func)} + if not comp: + return [inf] + if comp and inf not in xieta: + xieta.append(inf) + + factor = -hx/h + facsym = factor.free_symbols + if y not in facsym: + try: + fx = exp(integrate(factor, x)) + except NotImplementedError: + pass + else: + inf = {xi: fx, eta: S.Zero} + if not comp: + return [inf] + if comp and inf not in xieta: + xieta.append(inf) + + factor = -hx/(h**2) + facsym = factor.free_symbols + if x not in facsym: + try: + fy = exp(integrate(factor, y)) + except NotImplementedError: + pass + else: + inf = {xi: fy.subs(y, func), eta: S.Zero} + if not comp: + return [inf] + if comp and inf not in xieta: + xieta.append(inf) + + if xieta: + return xieta + +def lie_heuristic_abaco1_product(match, comp=False): + r""" + The second heuristic uses the following two assumptions on `\xi` and `\eta` + + .. math:: \eta = 0, \xi = f(x)*g(y) + + .. math:: \eta = f(x)*g(y), \xi = 0 + + The first assumption of this heuristic holds good if + `\frac{1}{h^{2}}\frac{\partial^2}{\partial x \partial y}\log(h)` is + separable in `x` and `y`, then the separated factors containing `x` + is `f(x)`, and `g(y)` is obtained by + + .. math:: e^{\int f\frac{\partial}{\partial x}\left(\frac{1}{f*h}\right)\,dy} + + provided `f\frac{\partial}{\partial x}\left(\frac{1}{f*h}\right)` is a function + of `y` only. + + The second assumption holds good if `\frac{dy}{dx} = h(x, y)` is rewritten as + `\frac{dy}{dx} = \frac{1}{h(y, x)}` and the same properties of the first assumption + satisfies. After obtaining `f(x)` and `g(y)`, the coordinates are again + interchanged, to get `\eta` as `f(x)*g(y)` + + + References + ========== + - E.S. Cheb-Terrab, A.D. Roche, Symmetries and First Order + ODE Patterns, pp. 7 - pp. 8 + + """ + + xieta = [] + y = match['y'] + h = match['h'] + hinv = match['hinv'] + func = match['func'] + x = func.args[0] + xi = Function('xi')(x, func) + eta = Function('eta')(x, func) + + + inf = separatevars(((log(h).diff(y)).diff(x))/h**2, dict=True, symbols=[x, y]) + if inf and inf['coeff']: + fx = inf[x] + gy = simplify(fx*((1/(fx*h)).diff(x))) + gysyms = gy.free_symbols + if x not in gysyms: + gy = exp(integrate(gy, y)) + inf = {eta: S.Zero, xi: (fx*gy).subs(y, func)} + if not comp: + return [inf] + if comp and inf not in xieta: + xieta.append(inf) + + u1 = Dummy("u1") + inf = separatevars(((log(hinv).diff(y)).diff(x))/hinv**2, dict=True, symbols=[x, y]) + if inf and inf['coeff']: + fx = inf[x] + gy = simplify(fx*((1/(fx*hinv)).diff(x))) + gysyms = gy.free_symbols + if x not in gysyms: + gy = exp(integrate(gy, y)) + etaval = fx*gy + etaval = (etaval.subs([(x, u1), (y, x)])).subs(u1, y) + inf = {eta: etaval.subs(y, func), xi: S.Zero} + if not comp: + return [inf] + if comp and inf not in xieta: + xieta.append(inf) + + if xieta: + return xieta + +def lie_heuristic_bivariate(match, comp=False): + r""" + The third heuristic assumes the infinitesimals `\xi` and `\eta` + to be bi-variate polynomials in `x` and `y`. The assumption made here + for the logic below is that `h` is a rational function in `x` and `y` + though that may not be necessary for the infinitesimals to be + bivariate polynomials. The coefficients of the infinitesimals + are found out by substituting them in the PDE and grouping similar terms + that are polynomials and since they form a linear system, solve and check + for non trivial solutions. The degree of the assumed bivariates + are increased till a certain maximum value. + + References + ========== + - Lie Groups and Differential Equations + pp. 327 - pp. 329 + + """ + + h = match['h'] + hx = match['hx'] + hy = match['hy'] + func = match['func'] + x = func.args[0] + y = match['y'] + xi = Function('xi')(x, func) + eta = Function('eta')(x, func) + + if h.is_rational_function(): + # The maximum degree that the infinitesimals can take is + # calculated by this technique. + etax, etay, etad, xix, xiy, xid = symbols("etax etay etad xix xiy xid") + ipde = etax + (etay - xix)*h - xiy*h**2 - xid*hx - etad*hy + num, denom = cancel(ipde).as_numer_denom() + deg = Poly(num, x, y).total_degree() + deta = Function('deta')(x, y) + dxi = Function('dxi')(x, y) + ipde = (deta.diff(x) + (deta.diff(y) - dxi.diff(x))*h - (dxi.diff(y))*h**2 + - dxi*hx - deta*hy) + xieq = Symbol("xi0") + etaeq = Symbol("eta0") + + for i in range(deg + 1): + if i: + xieq += Add(*[ + Symbol("xi_" + str(power) + "_" + str(i - power))*x**power*y**(i - power) + for power in range(i + 1)]) + etaeq += Add(*[ + Symbol("eta_" + str(power) + "_" + str(i - power))*x**power*y**(i - power) + for power in range(i + 1)]) + pden, denom = (ipde.subs({dxi: xieq, deta: etaeq}).doit()).as_numer_denom() + pden = expand(pden) + + # If the individual terms are monomials, the coefficients + # are grouped + if pden.is_polynomial(x, y) and pden.is_Add: + polyy = Poly(pden, x, y).as_dict() + if polyy: + symset = xieq.free_symbols.union(etaeq.free_symbols) - {x, y} + soldict = solve(polyy.values(), *symset) + if isinstance(soldict, list): + soldict = soldict[0] + if any(soldict.values()): + xired = xieq.subs(soldict) + etared = etaeq.subs(soldict) + # Scaling is done by substituting one for the parameters + # This can be any number except zero. + dict_ = {sym: 1 for sym in symset} + inf = {eta: etared.subs(dict_).subs(y, func), + xi: xired.subs(dict_).subs(y, func)} + return [inf] + +def lie_heuristic_chi(match, comp=False): + r""" + The aim of the fourth heuristic is to find the function `\chi(x, y)` + that satisfies the PDE `\frac{d\chi}{dx} + h\frac{d\chi}{dx} + - \frac{\partial h}{\partial y}\chi = 0`. + + This assumes `\chi` to be a bivariate polynomial in `x` and `y`. By intuition, + `h` should be a rational function in `x` and `y`. The method used here is + to substitute a general binomial for `\chi` up to a certain maximum degree + is reached. The coefficients of the polynomials, are calculated by by collecting + terms of the same order in `x` and `y`. + + After finding `\chi`, the next step is to use `\eta = \xi*h + \chi`, to + determine `\xi` and `\eta`. This can be done by dividing `\chi` by `h` + which would give `-\xi` as the quotient and `\eta` as the remainder. + + + References + ========== + - E.S Cheb-Terrab, L.G.S Duarte and L.A,C.P da Mota, Computer Algebra + Solving of First Order ODEs Using Symmetry Methods, pp. 8 + + """ + + h = match['h'] + hy = match['hy'] + func = match['func'] + x = func.args[0] + y = match['y'] + xi = Function('xi')(x, func) + eta = Function('eta')(x, func) + + if h.is_rational_function(): + schi, schix, schiy = symbols("schi, schix, schiy") + cpde = schix + h*schiy - hy*schi + num, denom = cancel(cpde).as_numer_denom() + deg = Poly(num, x, y).total_degree() + + chi = Function('chi')(x, y) + chix = chi.diff(x) + chiy = chi.diff(y) + cpde = chix + h*chiy - hy*chi + chieq = Symbol("chi") + for i in range(1, deg + 1): + chieq += Add(*[ + Symbol("chi_" + str(power) + "_" + str(i - power))*x**power*y**(i - power) + for power in range(i + 1)]) + cnum, cden = cancel(cpde.subs({chi : chieq}).doit()).as_numer_denom() + cnum = expand(cnum) + if cnum.is_polynomial(x, y) and cnum.is_Add: + cpoly = Poly(cnum, x, y).as_dict() + if cpoly: + solsyms = chieq.free_symbols - {x, y} + soldict = solve(cpoly.values(), *solsyms) + if isinstance(soldict, list): + soldict = soldict[0] + if any(soldict.values()): + chieq = chieq.subs(soldict) + dict_ = {sym: 1 for sym in solsyms} + chieq = chieq.subs(dict_) + # After finding chi, the main aim is to find out + # eta, xi by the equation eta = xi*h + chi + # One method to set xi, would be rearranging it to + # (eta/h) - xi = (chi/h). This would mean dividing + # chi by h would give -xi as the quotient and eta + # as the remainder. Thanks to Sean Vig for suggesting + # this method. + xic, etac = div(chieq, h) + inf = {eta: etac.subs(y, func), xi: -xic.subs(y, func)} + return [inf] + +def lie_heuristic_function_sum(match, comp=False): + r""" + This heuristic uses the following two assumptions on `\xi` and `\eta` + + .. math:: \eta = 0, \xi = f(x) + g(y) + + .. math:: \eta = f(x) + g(y), \xi = 0 + + The first assumption of this heuristic holds good if + + .. math:: \frac{\partial}{\partial y}[(h\frac{\partial^{2}}{ + \partial x^{2}}(h^{-1}))^{-1}] + + is separable in `x` and `y`, + + 1. The separated factors containing `y` is `\frac{\partial g}{\partial y}`. + From this `g(y)` can be determined. + 2. The separated factors containing `x` is `f''(x)`. + 3. `h\frac{\partial^{2}}{\partial x^{2}}(h^{-1})` equals + `\frac{f''(x)}{f(x) + g(y)}`. From this `f(x)` can be determined. + + The second assumption holds good if `\frac{dy}{dx} = h(x, y)` is rewritten as + `\frac{dy}{dx} = \frac{1}{h(y, x)}` and the same properties of the first + assumption satisfies. After obtaining `f(x)` and `g(y)`, the coordinates + are again interchanged, to get `\eta` as `f(x) + g(y)`. + + For both assumptions, the constant factors are separated among `g(y)` + and `f''(x)`, such that `f''(x)` obtained from 3] is the same as that + obtained from 2]. If not possible, then this heuristic fails. + + + References + ========== + - E.S. Cheb-Terrab, A.D. Roche, Symmetries and First Order + ODE Patterns, pp. 7 - pp. 8 + + """ + + xieta = [] + h = match['h'] + func = match['func'] + hinv = match['hinv'] + x = func.args[0] + y = match['y'] + xi = Function('xi')(x, func) + eta = Function('eta')(x, func) + + for odefac in [h, hinv]: + factor = odefac*((1/odefac).diff(x, 2)) + sep = separatevars((1/factor).diff(y), dict=True, symbols=[x, y]) + if sep and sep['coeff'] and sep[x].has(x) and sep[y].has(y): + k = Dummy("k") + try: + gy = k*integrate(sep[y], y) + except NotImplementedError: + pass + else: + fdd = 1/(k*sep[x]*sep['coeff']) + fx = simplify(fdd/factor - gy) + check = simplify(fx.diff(x, 2) - fdd) + if fx: + if not check: + fx = fx.subs(k, 1) + gy = (gy/k) + else: + sol = solve(check, k) + if sol: + sol = sol[0] + fx = fx.subs(k, sol) + gy = (gy/k)*sol + else: + continue + if odefac == hinv: # Inverse ODE + fx = fx.subs(x, y) + gy = gy.subs(y, x) + etaval = factor_terms(fx + gy) + if etaval.is_Mul: + etaval = Mul(*[arg for arg in etaval.args if arg.has(x, y)]) + if odefac == hinv: # Inverse ODE + inf = {eta: etaval.subs(y, func), xi : S.Zero} + else: + inf = {xi: etaval.subs(y, func), eta : S.Zero} + if not comp: + return [inf] + else: + xieta.append(inf) + + if xieta: + return xieta + +def lie_heuristic_abaco2_similar(match, comp=False): + r""" + This heuristic uses the following two assumptions on `\xi` and `\eta` + + .. math:: \eta = g(x), \xi = f(x) + + .. math:: \eta = f(y), \xi = g(y) + + For the first assumption, + + 1. First `\frac{\frac{\partial h}{\partial y}}{\frac{\partial^{2} h}{ + \partial yy}}` is calculated. Let us say this value is A + + 2. If this is constant, then `h` is matched to the form `A(x) + B(x)e^{ + \frac{y}{C}}` then, `\frac{e^{\int \frac{A(x)}{C} \,dx}}{B(x)}` gives `f(x)` + and `A(x)*f(x)` gives `g(x)` + + 3. Otherwise `\frac{\frac{\partial A}{\partial X}}{\frac{\partial A}{ + \partial Y}} = \gamma` is calculated. If + + a] `\gamma` is a function of `x` alone + + b] `\frac{\gamma\frac{\partial h}{\partial y} - \gamma'(x) - \frac{ + \partial h}{\partial x}}{h + \gamma} = G` is a function of `x` alone. + then, `e^{\int G \,dx}` gives `f(x)` and `-\gamma*f(x)` gives `g(x)` + + The second assumption holds good if `\frac{dy}{dx} = h(x, y)` is rewritten as + `\frac{dy}{dx} = \frac{1}{h(y, x)}` and the same properties of the first assumption + satisfies. After obtaining `f(x)` and `g(x)`, the coordinates are again + interchanged, to get `\xi` as `f(x^*)` and `\eta` as `g(y^*)` + + References + ========== + - E.S. Cheb-Terrab, A.D. Roche, Symmetries and First Order + ODE Patterns, pp. 10 - pp. 12 + + """ + + h = match['h'] + hx = match['hx'] + hy = match['hy'] + func = match['func'] + hinv = match['hinv'] + x = func.args[0] + y = match['y'] + xi = Function('xi')(x, func) + eta = Function('eta')(x, func) + + factor = cancel(h.diff(y)/h.diff(y, 2)) + factorx = factor.diff(x) + factory = factor.diff(y) + if not factor.has(x) and not factor.has(y): + A = Wild('A', exclude=[y]) + B = Wild('B', exclude=[y]) + C = Wild('C', exclude=[x, y]) + match = h.match(A + B*exp(y/C)) + try: + tau = exp(-integrate(match[A]/match[C]), x)/match[B] + except NotImplementedError: + pass + else: + gx = match[A]*tau + return [{xi: tau, eta: gx}] + + else: + gamma = cancel(factorx/factory) + if not gamma.has(y): + tauint = cancel((gamma*hy - gamma.diff(x) - hx)/(h + gamma)) + if not tauint.has(y): + try: + tau = exp(integrate(tauint, x)) + except NotImplementedError: + pass + else: + gx = -tau*gamma + return [{xi: tau, eta: gx}] + + factor = cancel(hinv.diff(y)/hinv.diff(y, 2)) + factorx = factor.diff(x) + factory = factor.diff(y) + if not factor.has(x) and not factor.has(y): + A = Wild('A', exclude=[y]) + B = Wild('B', exclude=[y]) + C = Wild('C', exclude=[x, y]) + match = h.match(A + B*exp(y/C)) + try: + tau = exp(-integrate(match[A]/match[C]), x)/match[B] + except NotImplementedError: + pass + else: + gx = match[A]*tau + return [{eta: tau.subs(x, func), xi: gx.subs(x, func)}] + + else: + gamma = cancel(factorx/factory) + if not gamma.has(y): + tauint = cancel((gamma*hinv.diff(y) - gamma.diff(x) - hinv.diff(x))/( + hinv + gamma)) + if not tauint.has(y): + try: + tau = exp(integrate(tauint, x)) + except NotImplementedError: + pass + else: + gx = -tau*gamma + return [{eta: tau.subs(x, func), xi: gx.subs(x, func)}] + + +def lie_heuristic_abaco2_unique_unknown(match, comp=False): + r""" + This heuristic assumes the presence of unknown functions or known functions + with non-integer powers. + + 1. A list of all functions and non-integer powers containing x and y + 2. Loop over each element `f` in the list, find `\frac{\frac{\partial f}{\partial x}}{ + \frac{\partial f}{\partial x}} = R` + + If it is separable in `x` and `y`, let `X` be the factors containing `x`. Then + + a] Check if `\xi = X` and `\eta = -\frac{X}{R}` satisfy the PDE. If yes, then return + `\xi` and `\eta` + b] Check if `\xi = \frac{-R}{X}` and `\eta = -\frac{1}{X}` satisfy the PDE. + If yes, then return `\xi` and `\eta` + + If not, then check if + + a] :math:`\xi = -R,\eta = 1` + + b] :math:`\xi = 1, \eta = -\frac{1}{R}` + + are solutions. + + References + ========== + - E.S. Cheb-Terrab, A.D. Roche, Symmetries and First Order + ODE Patterns, pp. 10 - pp. 12 + + """ + + h = match['h'] + hx = match['hx'] + hy = match['hy'] + func = match['func'] + x = func.args[0] + y = match['y'] + xi = Function('xi')(x, func) + eta = Function('eta')(x, func) + + funclist = [] + for atom in h.atoms(Pow): + base, exp = atom.as_base_exp() + if base.has(x) and base.has(y): + if not exp.is_Integer: + funclist.append(atom) + + for function in h.atoms(AppliedUndef): + syms = function.free_symbols + if x in syms and y in syms: + funclist.append(function) + + for f in funclist: + frac = cancel(f.diff(y)/f.diff(x)) + sep = separatevars(frac, dict=True, symbols=[x, y]) + if sep and sep['coeff']: + xitry1 = sep[x] + etatry1 = -1/(sep[y]*sep['coeff']) + pde1 = etatry1.diff(y)*h - xitry1.diff(x)*h - xitry1*hx - etatry1*hy + if not simplify(pde1): + return [{xi: xitry1, eta: etatry1.subs(y, func)}] + xitry2 = 1/etatry1 + etatry2 = 1/xitry1 + pde2 = etatry2.diff(x) - (xitry2.diff(y))*h**2 - xitry2*hx - etatry2*hy + if not simplify(expand(pde2)): + return [{xi: xitry2.subs(y, func), eta: etatry2}] + + else: + etatry = -1/frac + pde = etatry.diff(x) + etatry.diff(y)*h - hx - etatry*hy + if not simplify(pde): + return [{xi: S.One, eta: etatry.subs(y, func)}] + xitry = -frac + pde = -xitry.diff(x)*h -xitry.diff(y)*h**2 - xitry*hx -hy + if not simplify(expand(pde)): + return [{xi: xitry.subs(y, func), eta: S.One}] + + +def lie_heuristic_abaco2_unique_general(match, comp=False): + r""" + This heuristic finds if infinitesimals of the form `\eta = f(x)`, `\xi = g(y)` + without making any assumptions on `h`. + + The complete sequence of steps is given in the paper mentioned below. + + References + ========== + - E.S. Cheb-Terrab, A.D. Roche, Symmetries and First Order + ODE Patterns, pp. 10 - pp. 12 + + """ + hx = match['hx'] + hy = match['hy'] + func = match['func'] + x = func.args[0] + y = match['y'] + xi = Function('xi')(x, func) + eta = Function('eta')(x, func) + + A = hx.diff(y) + B = hy.diff(y) + hy**2 + C = hx.diff(x) - hx**2 + + if not (A and B and C): + return + + Ax = A.diff(x) + Ay = A.diff(y) + Axy = Ax.diff(y) + Axx = Ax.diff(x) + Ayy = Ay.diff(y) + D = simplify(2*Axy + hx*Ay - Ax*hy + (hx*hy + 2*A)*A)*A - 3*Ax*Ay + if not D: + E1 = simplify(3*Ax**2 + ((hx**2 + 2*C)*A - 2*Axx)*A) + if E1: + E2 = simplify((2*Ayy + (2*B - hy**2)*A)*A - 3*Ay**2) + if not E2: + E3 = simplify( + E1*((28*Ax + 4*hx*A)*A**3 - E1*(hy*A + Ay)) - E1.diff(x)*8*A**4) + if not E3: + etaval = cancel((4*A**3*(Ax - hx*A) + E1*(hy*A - Ay))/(S(2)*A*E1)) + if x not in etaval: + try: + etaval = exp(integrate(etaval, y)) + except NotImplementedError: + pass + else: + xival = -4*A**3*etaval/E1 + if y not in xival: + return [{xi: xival, eta: etaval.subs(y, func)}] + + else: + E1 = simplify((2*Ayy + (2*B - hy**2)*A)*A - 3*Ay**2) + if E1: + E2 = simplify( + 4*A**3*D - D**2 + E1*((2*Axx - (hx**2 + 2*C)*A)*A - 3*Ax**2)) + if not E2: + E3 = simplify( + -(A*D)*E1.diff(y) + ((E1.diff(x) - hy*D)*A + 3*Ay*D + + (A*hx - 3*Ax)*E1)*E1) + if not E3: + etaval = cancel(((A*hx - Ax)*E1 - (Ay + A*hy)*D)/(S(2)*A*D)) + if x not in etaval: + try: + etaval = exp(integrate(etaval, y)) + except NotImplementedError: + pass + else: + xival = -E1*etaval/D + if y not in xival: + return [{xi: xival, eta: etaval.subs(y, func)}] + + +def lie_heuristic_linear(match, comp=False): + r""" + This heuristic assumes + + 1. `\xi = ax + by + c` and + 2. `\eta = fx + gy + h` + + After substituting the following assumptions in the determining PDE, it + reduces to + + .. math:: f + (g - a)h - bh^{2} - (ax + by + c)\frac{\partial h}{\partial x} + - (fx + gy + c)\frac{\partial h}{\partial y} + + Solving the reduced PDE obtained, using the method of characteristics, becomes + impractical. The method followed is grouping similar terms and solving the system + of linear equations obtained. The difference between the bivariate heuristic is that + `h` need not be a rational function in this case. + + References + ========== + - E.S. Cheb-Terrab, A.D. Roche, Symmetries and First Order + ODE Patterns, pp. 10 - pp. 12 + + """ + h = match['h'] + hx = match['hx'] + hy = match['hy'] + func = match['func'] + x = func.args[0] + y = match['y'] + xi = Function('xi')(x, func) + eta = Function('eta')(x, func) + + coeffdict = {} + symbols = numbered_symbols("c", cls=Dummy) + symlist = [next(symbols) for _ in islice(symbols, 6)] + C0, C1, C2, C3, C4, C5 = symlist + pde = C3 + (C4 - C0)*h - (C0*x + C1*y + C2)*hx - (C3*x + C4*y + C5)*hy - C1*h**2 + pde, denom = pde.as_numer_denom() + pde = powsimp(expand(pde)) + if pde.is_Add: + terms = pde.args + for term in terms: + if term.is_Mul: + rem = Mul(*[m for m in term.args if not m.has(x, y)]) + xypart = term/rem + if xypart not in coeffdict: + coeffdict[xypart] = rem + else: + coeffdict[xypart] += rem + else: + if term not in coeffdict: + coeffdict[term] = S.One + else: + coeffdict[term] += S.One + + sollist = coeffdict.values() + soldict = solve(sollist, symlist) + if soldict: + if isinstance(soldict, list): + soldict = soldict[0] + subval = soldict.values() + if any(t for t in subval): + onedict = dict(zip(symlist, [1]*6)) + xival = C0*x + C1*func + C2 + etaval = C3*x + C4*func + C5 + xival = xival.subs(soldict) + etaval = etaval.subs(soldict) + xival = xival.subs(onedict) + etaval = etaval.subs(onedict) + return [{xi: xival, eta: etaval}] + + +def _lie_group_remove(coords): + r""" + This function is strictly meant for internal use by the Lie group ODE solving + method. It replaces arbitrary functions returned by pdsolve as follows: + + 1] If coords is an arbitrary function, then its argument is returned. + 2] An arbitrary function in an Add object is replaced by zero. + 3] An arbitrary function in a Mul object is replaced by one. + 4] If there is no arbitrary function coords is returned unchanged. + + Examples + ======== + + >>> from sympy.solvers.ode.lie_group import _lie_group_remove + >>> from sympy import Function + >>> from sympy.abc import x, y + >>> F = Function("F") + >>> eq = x**2*y + >>> _lie_group_remove(eq) + x**2*y + >>> eq = F(x**2*y) + >>> _lie_group_remove(eq) + x**2*y + >>> eq = x*y**2 + F(x**3) + >>> _lie_group_remove(eq) + x*y**2 + >>> eq = (F(x**3) + y)*x**4 + >>> _lie_group_remove(eq) + x**4*y + + """ + if isinstance(coords, AppliedUndef): + return coords.args[0] + elif coords.is_Add: + subfunc = coords.atoms(AppliedUndef) + if subfunc: + for func in subfunc: + coords = coords.subs(func, 0) + return coords + elif coords.is_Pow: + base, expr = coords.as_base_exp() + base = _lie_group_remove(base) + expr = _lie_group_remove(expr) + return base**expr + elif coords.is_Mul: + mulargs = [] + coordargs = coords.args + for arg in coordargs: + if not isinstance(coords, AppliedUndef): + mulargs.append(_lie_group_remove(arg)) + return Mul(*mulargs) + return coords diff --git a/sympy/solvers/ode/ode.py b/sympy/solvers/ode/ode.py index 438336c3263f..7611783202a5 100644 --- a/sympy/solvers/ode/ode.py +++ b/sympy/solvers/ode/ode.py @@ -227,12 +227,10 @@ of those tests will surely fail. """ -from itertools import islice from sympy.core import Add, S, Mul, Pow, oo from sympy.core.compatibility import ordered, iterable from sympy.core.containers import Tuple -from sympy.core.exprtools import factor_terms from sympy.core.expr import AtomicExpr, Expr from sympy.core.function import (Function, Derivative, AppliedUndef, diff, expand, expand_mul, Subs) @@ -246,16 +244,15 @@ BooleanFalse) from sympy.functions import exp, log, sqrt from sympy.functions.combinatorial.factorials import factorial -from sympy.integrals.integrals import Integral, integrate +from sympy.integrals.integrals import Integral from sympy.polys import (Poly, terms_gcd, PolynomialError, lcm) -from sympy.polys.polytools import cancel, div +from sympy.polys.polytools import cancel from sympy.series import Order from sympy.series.series import series from sympy.simplify import (collect, logcombine, powsimp, # type: ignore separatevars, simplify, cse) from sympy.simplify.radsimp import collect_const from sympy.solvers import checksol, solve -from sympy.solvers.pde import pdsolve from sympy.utilities import numbered_symbols, default_sort_key, sift from sympy.utilities.iterables import uniq @@ -323,18 +320,6 @@ "2nd_nonlinear_autonomous_conserved_Integral", ) -lie_heuristics = ( - "abaco1_simple", - "abaco1_product", - "abaco2_similar", - "abaco2_unique_unknown", - "abaco2_unique_general", - "linear", - "function_sum", - "bivariate", - "chi" - ) - def get_numbered_constants(eq, num=1, start=1, prefix='C'): @@ -985,7 +970,6 @@ class in it. Note that a hint may do this anyway if a3 = Wild('a3', exclude=[f(x), df, f(x).diff(x, 2)]) b3 = Wild('b3', exclude=[f(x), df, f(x).diff(x, 2)]) c3 = Wild('c3', exclude=[f(x), df, f(x).diff(x, 2)]) - r3 = {'xi': xi, 'eta': eta} # Used for the lie_group hint boundary = {} # Used to extract initial conditions C1 = Symbol("C1") @@ -1033,7 +1017,7 @@ class in it. Note that a hint may do this anyway if # Any ODE that can be solved with a combination of algebra and # integrals e.g.: # d^3/dx^3(x y) = F(x) - ode = SingleODEProblem(eq_orig, func, x, prep=prep) + ode = SingleODEProblem(eq_orig, func, x, prep=prep, xi=xi, eta=eta) solvers = { NthAlgebraic: ('nth_algebraic',), FirstExact:('1st_exact',), @@ -1060,6 +1044,7 @@ class in it. Note that a hint may do this anyway if NthLinearEulerEqNonhomogeneousUndeterminedCoefficients: ('nth_linear_euler_eq_nonhomogeneous_undetermined_coefficients',), SecondLinearBessel: ('2nd_linear_bessel',), SecondLinearAiry: ('2nd_linear_airy',), + LieGroup: ('lie_group',), } for solvercls in solvers: solver = solvercls(ode) @@ -1110,12 +1095,6 @@ class in it. Note that a hint may do this anyway if rseries.update({'terms': terms, 'f0': point, 'f0val': value}) matching_hints["1st_power_series"] = rseries - r3.update(r) - - # Any first order ODE can be ideally solved by the Lie Group - # method - matching_hints["lie_group"] = r3 - elif order == 2: # Homogeneous second order differential equation of the form # a3*f(x).diff(x, 2) + b3*f(x).diff(x) + c3 @@ -2752,54 +2731,6 @@ def _is_special_case_of(soln1, soln2, eq, order, var): return False -def _nth_linear_match(eq, func, order): - r""" - Matches a differential equation to the linear form: - - .. math:: a_n(x) y^{(n)} + \cdots + a_1(x)y' + a_0(x) y + B(x) = 0 - - Returns a dict of order:coeff terms, where order is the order of the - derivative on each term, and coeff is the coefficient of that derivative. - The key ``-1`` holds the function `B(x)`. Returns ``None`` if the ODE is - not linear. This function assumes that ``func`` has already been checked - to be good. - - Examples - ======== - - >>> from sympy import Function, cos, sin - >>> from sympy.abc import x - >>> from sympy.solvers.ode.ode import _nth_linear_match - >>> f = Function('f') - >>> _nth_linear_match(f(x).diff(x, 3) + 2*f(x).diff(x) + - ... x*f(x).diff(x, 2) + cos(x)*f(x).diff(x) + x - f(x) - - ... sin(x), f(x), 3) - {-1: x - sin(x), 0: -1, 1: cos(x) + 2, 2: x, 3: 1} - >>> _nth_linear_match(f(x).diff(x, 3) + 2*f(x).diff(x) + - ... x*f(x).diff(x, 2) + cos(x)*f(x).diff(x) + x - f(x) - - ... sin(f(x)), f(x), 3) == None - True - - """ - x = func.args[0] - one_x = {x} - terms = {i: S.Zero for i in range(-1, order + 1)} - for i in Add.make_args(eq): - if not i.has(func): - terms[-1] += i - else: - c, f = i.as_independent(func) - if (isinstance(f, Derivative) - and set(f.variables) == one_x - and f.args[0] == func): - terms[f.derivative_count] += c - elif f == func: - terms[len(f.args[1:])] += c - else: - return None - return terms - - def ode_1st_power_series(eq, func, order, match): r""" The power series solution is a method which gives the Taylor series expansion @@ -2951,1121 +2882,6 @@ def checkinfsol(eq, infinitesimals, func=None, order=None): soltup.append((True, 0)) return soltup -def _ode_lie_group_try_heuristic(eq, heuristic, func, match, inf): - - xi = Function("xi") - eta = Function("eta") - f = func.func - x = func.args[0] - y = match['y'] - h = match['h'] - tempsol = [] - if not inf: - try: - inf = infinitesimals(eq, hint=heuristic, func=func, order=1, match=match) - except ValueError: - return None - for infsim in inf: - xiinf = (infsim[xi(x, func)]).subs(func, y) - etainf = (infsim[eta(x, func)]).subs(func, y) - # This condition creates recursion while using pdsolve. - # Since the first step while solving a PDE of form - # a*(f(x, y).diff(x)) + b*(f(x, y).diff(y)) + c = 0 - # is to solve the ODE dy/dx = b/a - if simplify(etainf/xiinf) == h: - continue - rpde = f(x, y).diff(x)*xiinf + f(x, y).diff(y)*etainf - r = pdsolve(rpde, func=f(x, y)).rhs - s = pdsolve(rpde - 1, func=f(x, y)).rhs - newcoord = [_lie_group_remove(coord) for coord in [r, s]] - r = Dummy("r") - s = Dummy("s") - C1 = Symbol("C1") - rcoord = newcoord[0] - scoord = newcoord[-1] - try: - sol = solve([r - rcoord, s - scoord], x, y, dict=True) - if sol == []: - continue - except NotImplementedError: - continue - else: - sol = sol[0] - xsub = sol[x] - ysub = sol[y] - num = simplify(scoord.diff(x) + scoord.diff(y)*h) - denom = simplify(rcoord.diff(x) + rcoord.diff(y)*h) - if num and denom: - diffeq = simplify((num/denom).subs([(x, xsub), (y, ysub)])) - sep = separatevars(diffeq, symbols=[r, s], dict=True) - if sep: - # Trying to separate, r and s coordinates - deq = integrate((1/sep[s]), s) + C1 - integrate(sep['coeff']*sep[r], r) - # Substituting and reverting back to original coordinates - deq = deq.subs([(r, rcoord), (s, scoord)]) - try: - sdeq = solve(deq, y) - except NotImplementedError: - tempsol.append(deq) - else: - return [Eq(f(x), sol) for sol in sdeq] - - - elif denom: # (ds/dr) is zero which means s is constant - return [Eq(f(x), solve(scoord - C1, y)[0])] - - elif num: # (dr/ds) is zero which means r is constant - return [Eq(f(x), solve(rcoord - C1, y)[0])] - - # If nothing works, return solution as it is, without solving for y - if tempsol: - return [Eq(sol.subs(y, f(x)), 0) for sol in tempsol] - return None - -def _ode_lie_group( s, func, order, match): - - heuristics = lie_heuristics - inf = {} - f = func.func - x = func.args[0] - df = func.diff(x) - xi = Function("xi") - eta = Function("eta") - xis = match['xi'] - etas = match['eta'] - y = match.pop('y', None) - if y: - h = -simplify(match[match['d']]/match[match['e']]) - y = y - else: - y = Dummy("y") - h = s.subs(func, y) - - if xis is not None and etas is not None: - inf = [{xi(x, f(x)): S(xis), eta(x, f(x)): S(etas)}] - - if checkinfsol(Eq(df, s), inf, func=f(x), order=1)[0][0]: - heuristics = ["user_defined"] + list(heuristics) - - match = {'h': h, 'y': y} - - # This is done so that if any heuristic raises a ValueError - # another heuristic can be used. - sol = None - for heuristic in heuristics: - sol = _ode_lie_group_try_heuristic(Eq(df, s), heuristic, func, match, inf) - if sol: - return sol - return sol - -def ode_lie_group(eq, func, order, match): - r""" - This hint implements the Lie group method of solving first order differential - equations. The aim is to convert the given differential equation from the - given coordinate system into another coordinate system where it becomes - invariant under the one-parameter Lie group of translations. The converted - ODE can be easily solved by quadrature. It makes use of the - :py:meth:`sympy.solvers.ode.infinitesimals` function which returns the - infinitesimals of the transformation. - - The coordinates `r` and `s` can be found by solving the following Partial - Differential Equations. - - .. math :: \xi\frac{\partial r}{\partial x} + \eta\frac{\partial r}{\partial y} - = 0 - - .. math :: \xi\frac{\partial s}{\partial x} + \eta\frac{\partial s}{\partial y} - = 1 - - The differential equation becomes separable in the new coordinate system - - .. math :: \frac{ds}{dr} = \frac{\frac{\partial s}{\partial x} + - h(x, y)\frac{\partial s}{\partial y}}{ - \frac{\partial r}{\partial x} + h(x, y)\frac{\partial r}{\partial y}} - - After finding the solution by integration, it is then converted back to the original - coordinate system by substituting `r` and `s` in terms of `x` and `y` again. - - Examples - ======== - - >>> from sympy import Function, dsolve, exp, pprint - >>> from sympy.abc import x - >>> f = Function('f') - >>> pprint(dsolve(f(x).diff(x) + 2*x*f(x) - x*exp(-x**2), f(x), - ... hint='lie_group')) - / 2\ 2 - | x | -x - f(x) = |C1 + --|*e - \ 2 / - - - References - ========== - - - Solving differential equations by Symmetry Groups, - John Starrett, pp. 1 - pp. 14 - - """ - - x = func.args[0] - df = func.diff(x) - - try: - eqsol = solve(eq, df) - except NotImplementedError: - eqsol = [] - - desols = [] - for s in eqsol: - sol = _ode_lie_group(s, func, order, match=match) - if sol: - desols.extend(sol) - - if desols == []: - raise NotImplementedError("The given ODE " + str(eq) + " cannot be solved by" - + " the lie group method") - return desols - -def _lie_group_remove(coords): - r""" - This function is strictly meant for internal use by the Lie group ODE solving - method. It replaces arbitrary functions returned by pdsolve as follows: - - 1] If coords is an arbitrary function, then its argument is returned. - 2] An arbitrary function in an Add object is replaced by zero. - 3] An arbitrary function in a Mul object is replaced by one. - 4] If there is no arbitrary function coords is returned unchanged. - - Examples - ======== - - >>> from sympy.solvers.ode.ode import _lie_group_remove - >>> from sympy import Function - >>> from sympy.abc import x, y - >>> F = Function("F") - >>> eq = x**2*y - >>> _lie_group_remove(eq) - x**2*y - >>> eq = F(x**2*y) - >>> _lie_group_remove(eq) - x**2*y - >>> eq = x*y**2 + F(x**3) - >>> _lie_group_remove(eq) - x*y**2 - >>> eq = (F(x**3) + y)*x**4 - >>> _lie_group_remove(eq) - x**4*y - - """ - if isinstance(coords, AppliedUndef): - return coords.args[0] - elif coords.is_Add: - subfunc = coords.atoms(AppliedUndef) - if subfunc: - for func in subfunc: - coords = coords.subs(func, 0) - return coords - elif coords.is_Pow: - base, expr = coords.as_base_exp() - base = _lie_group_remove(base) - expr = _lie_group_remove(expr) - return base**expr - elif coords.is_Mul: - mulargs = [] - coordargs = coords.args - for arg in coordargs: - if not isinstance(coords, AppliedUndef): - mulargs.append(_lie_group_remove(arg)) - return Mul(*mulargs) - return coords - -def infinitesimals(eq, func=None, order=None, hint='default', match=None): - r""" - The infinitesimal functions of an ordinary differential equation, `\xi(x,y)` - and `\eta(x,y)`, are the infinitesimals of the Lie group of point transformations - for which the differential equation is invariant. So, the ODE `y'=f(x,y)` - would admit a Lie group `x^*=X(x,y;\varepsilon)=x+\varepsilon\xi(x,y)`, - `y^*=Y(x,y;\varepsilon)=y+\varepsilon\eta(x,y)` such that `(y^*)'=f(x^*, y^*)`. - A change of coordinates, to `r(x,y)` and `s(x,y)`, can be performed so this Lie group - becomes the translation group, `r^*=r` and `s^*=s+\varepsilon`. - They are tangents to the coordinate curves of the new system. - - Consider the transformation `(x, y) \to (X, Y)` such that the - differential equation remains invariant. `\xi` and `\eta` are the tangents to - the transformed coordinates `X` and `Y`, at `\varepsilon=0`. - - .. math:: \left(\frac{\partial X(x,y;\varepsilon)}{\partial\varepsilon - }\right)|_{\varepsilon=0} = \xi, - \left(\frac{\partial Y(x,y;\varepsilon)}{\partial\varepsilon - }\right)|_{\varepsilon=0} = \eta, - - The infinitesimals can be found by solving the following PDE: - - >>> from sympy import Function, Eq, pprint - >>> from sympy.abc import x, y - >>> xi, eta, h = map(Function, ['xi', 'eta', 'h']) - >>> h = h(x, y) # dy/dx = h - >>> eta = eta(x, y) - >>> xi = xi(x, y) - >>> genform = Eq(eta.diff(x) + (eta.diff(y) - xi.diff(x))*h - ... - (xi.diff(y))*h**2 - xi*(h.diff(x)) - eta*(h.diff(y)), 0) - >>> pprint(genform) - /d d \ d 2 d - |--(eta(x, y)) - --(xi(x, y))|*h(x, y) - eta(x, y)*--(h(x, y)) - h (x, y)*--(x - \dy dx / dy dy - - d d - i(x, y)) - xi(x, y)*--(h(x, y)) + --(eta(x, y)) = 0 - dx dx - - Solving the above mentioned PDE is not trivial, and can be solved only by - making intelligent assumptions for `\xi` and `\eta` (heuristics). Once an - infinitesimal is found, the attempt to find more heuristics stops. This is done to - optimise the speed of solving the differential equation. If a list of all the - infinitesimals is needed, ``hint`` should be flagged as ``all``, which gives - the complete list of infinitesimals. If the infinitesimals for a particular - heuristic needs to be found, it can be passed as a flag to ``hint``. - - Examples - ======== - - >>> from sympy import Function - >>> from sympy.solvers.ode.ode import infinitesimals - >>> from sympy.abc import x - >>> f = Function('f') - >>> eq = f(x).diff(x) - x**2*f(x) - >>> infinitesimals(eq) - [{eta(x, f(x)): exp(x**3/3), xi(x, f(x)): 0}] - - References - ========== - - - Solving differential equations by Symmetry Groups, - John Starrett, pp. 1 - pp. 14 - - """ - - if isinstance(eq, Equality): - eq = eq.lhs - eq.rhs - if not func: - eq, func = _preprocess(eq) - variables = func.args - if len(variables) != 1: - raise ValueError("ODE's have only one independent variable") - else: - x = variables[0] - if not order: - order = ode_order(eq, func) - if order != 1: - raise NotImplementedError("Infinitesimals for only " - "first order ODE's have been implemented") - else: - df = func.diff(x) - # Matching differential equation of the form a*df + b - a = Wild('a', exclude = [df]) - b = Wild('b', exclude = [df]) - if match: # Used by lie_group hint - h = match['h'] - y = match['y'] - else: - match = collect(expand(eq), df).match(a*df + b) - if match: - h = -simplify(match[b]/match[a]) - else: - try: - sol = solve(eq, df) - except NotImplementedError: - raise NotImplementedError("Infinitesimals for the " - "first order ODE could not be found") - else: - h = sol[0] # Find infinitesimals for one solution - y = Dummy("y") - h = h.subs(func, y) - - u = Dummy("u") - hx = h.diff(x) - hy = h.diff(y) - hinv = ((1/h).subs([(x, u), (y, x)])).subs(u, y) # Inverse ODE - match = {'h': h, 'func': func, 'hx': hx, 'hy': hy, 'y': y, 'hinv': hinv} - if hint == 'all': - xieta = [] - for heuristic in lie_heuristics: - function = globals()['lie_heuristic_' + heuristic] - inflist = function(match, comp=True) - if inflist: - xieta.extend([inf for inf in inflist if inf not in xieta]) - if xieta: - return xieta - else: - raise NotImplementedError("Infinitesimals could not be found for " - "the given ODE") - - elif hint == 'default': - for heuristic in lie_heuristics: - function = globals()['lie_heuristic_' + heuristic] - xieta = function(match, comp=False) - if xieta: - return xieta - - raise NotImplementedError("Infinitesimals could not be found for" - " the given ODE") - - elif hint not in lie_heuristics: - raise ValueError("Heuristic not recognized: " + hint) - - else: - function = globals()['lie_heuristic_' + hint] - xieta = function(match, comp=True) - if xieta: - return xieta - else: - raise ValueError("Infinitesimals could not be found using the" - " given heuristic") - - -def lie_heuristic_abaco1_simple(match, comp=False): - r""" - The first heuristic uses the following four sets of - assumptions on `\xi` and `\eta` - - .. math:: \xi = 0, \eta = f(x) - - .. math:: \xi = 0, \eta = f(y) - - .. math:: \xi = f(x), \eta = 0 - - .. math:: \xi = f(y), \eta = 0 - - The success of this heuristic is determined by algebraic factorisation. - For the first assumption `\xi = 0` and `\eta` to be a function of `x`, the PDE - - .. math:: \frac{\partial \eta}{\partial x} + (\frac{\partial \eta}{\partial y} - - \frac{\partial \xi}{\partial x})*h - - \frac{\partial \xi}{\partial y}*h^{2} - - \xi*\frac{\partial h}{\partial x} - \eta*\frac{\partial h}{\partial y} = 0 - - reduces to `f'(x) - f\frac{\partial h}{\partial y} = 0` - If `\frac{\partial h}{\partial y}` is a function of `x`, then this can usually - be integrated easily. A similar idea is applied to the other 3 assumptions as well. - - - References - ========== - - - E.S Cheb-Terrab, L.G.S Duarte and L.A,C.P da Mota, Computer Algebra - Solving of First Order ODEs Using Symmetry Methods, pp. 8 - - - """ - - xieta = [] - y = match['y'] - h = match['h'] - func = match['func'] - x = func.args[0] - hx = match['hx'] - hy = match['hy'] - xi = Function('xi')(x, func) - eta = Function('eta')(x, func) - - hysym = hy.free_symbols - if y not in hysym: - try: - fx = exp(integrate(hy, x)) - except NotImplementedError: - pass - else: - inf = {xi: S.Zero, eta: fx} - if not comp: - return [inf] - if comp and inf not in xieta: - xieta.append(inf) - - factor = hy/h - facsym = factor.free_symbols - if x not in facsym: - try: - fy = exp(integrate(factor, y)) - except NotImplementedError: - pass - else: - inf = {xi: S.Zero, eta: fy.subs(y, func)} - if not comp: - return [inf] - if comp and inf not in xieta: - xieta.append(inf) - - factor = -hx/h - facsym = factor.free_symbols - if y not in facsym: - try: - fx = exp(integrate(factor, x)) - except NotImplementedError: - pass - else: - inf = {xi: fx, eta: S.Zero} - if not comp: - return [inf] - if comp and inf not in xieta: - xieta.append(inf) - - factor = -hx/(h**2) - facsym = factor.free_symbols - if x not in facsym: - try: - fy = exp(integrate(factor, y)) - except NotImplementedError: - pass - else: - inf = {xi: fy.subs(y, func), eta: S.Zero} - if not comp: - return [inf] - if comp and inf not in xieta: - xieta.append(inf) - - if xieta: - return xieta - -def lie_heuristic_abaco1_product(match, comp=False): - r""" - The second heuristic uses the following two assumptions on `\xi` and `\eta` - - .. math:: \eta = 0, \xi = f(x)*g(y) - - .. math:: \eta = f(x)*g(y), \xi = 0 - - The first assumption of this heuristic holds good if - `\frac{1}{h^{2}}\frac{\partial^2}{\partial x \partial y}\log(h)` is - separable in `x` and `y`, then the separated factors containing `x` - is `f(x)`, and `g(y)` is obtained by - - .. math:: e^{\int f\frac{\partial}{\partial x}\left(\frac{1}{f*h}\right)\,dy} - - provided `f\frac{\partial}{\partial x}\left(\frac{1}{f*h}\right)` is a function - of `y` only. - - The second assumption holds good if `\frac{dy}{dx} = h(x, y)` is rewritten as - `\frac{dy}{dx} = \frac{1}{h(y, x)}` and the same properties of the first assumption - satisfies. After obtaining `f(x)` and `g(y)`, the coordinates are again - interchanged, to get `\eta` as `f(x)*g(y)` - - - References - ========== - - E.S. Cheb-Terrab, A.D. Roche, Symmetries and First Order - ODE Patterns, pp. 7 - pp. 8 - - """ - - xieta = [] - y = match['y'] - h = match['h'] - hinv = match['hinv'] - func = match['func'] - x = func.args[0] - xi = Function('xi')(x, func) - eta = Function('eta')(x, func) - - - inf = separatevars(((log(h).diff(y)).diff(x))/h**2, dict=True, symbols=[x, y]) - if inf and inf['coeff']: - fx = inf[x] - gy = simplify(fx*((1/(fx*h)).diff(x))) - gysyms = gy.free_symbols - if x not in gysyms: - gy = exp(integrate(gy, y)) - inf = {eta: S.Zero, xi: (fx*gy).subs(y, func)} - if not comp: - return [inf] - if comp and inf not in xieta: - xieta.append(inf) - - u1 = Dummy("u1") - inf = separatevars(((log(hinv).diff(y)).diff(x))/hinv**2, dict=True, symbols=[x, y]) - if inf and inf['coeff']: - fx = inf[x] - gy = simplify(fx*((1/(fx*hinv)).diff(x))) - gysyms = gy.free_symbols - if x not in gysyms: - gy = exp(integrate(gy, y)) - etaval = fx*gy - etaval = (etaval.subs([(x, u1), (y, x)])).subs(u1, y) - inf = {eta: etaval.subs(y, func), xi: S.Zero} - if not comp: - return [inf] - if comp and inf not in xieta: - xieta.append(inf) - - if xieta: - return xieta - -def lie_heuristic_bivariate(match, comp=False): - r""" - The third heuristic assumes the infinitesimals `\xi` and `\eta` - to be bi-variate polynomials in `x` and `y`. The assumption made here - for the logic below is that `h` is a rational function in `x` and `y` - though that may not be necessary for the infinitesimals to be - bivariate polynomials. The coefficients of the infinitesimals - are found out by substituting them in the PDE and grouping similar terms - that are polynomials and since they form a linear system, solve and check - for non trivial solutions. The degree of the assumed bivariates - are increased till a certain maximum value. - - References - ========== - - Lie Groups and Differential Equations - pp. 327 - pp. 329 - - """ - - h = match['h'] - hx = match['hx'] - hy = match['hy'] - func = match['func'] - x = func.args[0] - y = match['y'] - xi = Function('xi')(x, func) - eta = Function('eta')(x, func) - - if h.is_rational_function(): - # The maximum degree that the infinitesimals can take is - # calculated by this technique. - etax, etay, etad, xix, xiy, xid = symbols("etax etay etad xix xiy xid") - ipde = etax + (etay - xix)*h - xiy*h**2 - xid*hx - etad*hy - num, denom = cancel(ipde).as_numer_denom() - deg = Poly(num, x, y).total_degree() - deta = Function('deta')(x, y) - dxi = Function('dxi')(x, y) - ipde = (deta.diff(x) + (deta.diff(y) - dxi.diff(x))*h - (dxi.diff(y))*h**2 - - dxi*hx - deta*hy) - xieq = Symbol("xi0") - etaeq = Symbol("eta0") - - for i in range(deg + 1): - if i: - xieq += Add(*[ - Symbol("xi_" + str(power) + "_" + str(i - power))*x**power*y**(i - power) - for power in range(i + 1)]) - etaeq += Add(*[ - Symbol("eta_" + str(power) + "_" + str(i - power))*x**power*y**(i - power) - for power in range(i + 1)]) - pden, denom = (ipde.subs({dxi: xieq, deta: etaeq}).doit()).as_numer_denom() - pden = expand(pden) - - # If the individual terms are monomials, the coefficients - # are grouped - if pden.is_polynomial(x, y) and pden.is_Add: - polyy = Poly(pden, x, y).as_dict() - if polyy: - symset = xieq.free_symbols.union(etaeq.free_symbols) - {x, y} - soldict = solve(polyy.values(), *symset) - if isinstance(soldict, list): - soldict = soldict[0] - if any(soldict.values()): - xired = xieq.subs(soldict) - etared = etaeq.subs(soldict) - # Scaling is done by substituting one for the parameters - # This can be any number except zero. - dict_ = {sym: 1 for sym in symset} - inf = {eta: etared.subs(dict_).subs(y, func), - xi: xired.subs(dict_).subs(y, func)} - return [inf] - -def lie_heuristic_chi(match, comp=False): - r""" - The aim of the fourth heuristic is to find the function `\chi(x, y)` - that satisfies the PDE `\frac{d\chi}{dx} + h\frac{d\chi}{dx} - - \frac{\partial h}{\partial y}\chi = 0`. - - This assumes `\chi` to be a bivariate polynomial in `x` and `y`. By intuition, - `h` should be a rational function in `x` and `y`. The method used here is - to substitute a general binomial for `\chi` up to a certain maximum degree - is reached. The coefficients of the polynomials, are calculated by by collecting - terms of the same order in `x` and `y`. - - After finding `\chi`, the next step is to use `\eta = \xi*h + \chi`, to - determine `\xi` and `\eta`. This can be done by dividing `\chi` by `h` - which would give `-\xi` as the quotient and `\eta` as the remainder. - - - References - ========== - - E.S Cheb-Terrab, L.G.S Duarte and L.A,C.P da Mota, Computer Algebra - Solving of First Order ODEs Using Symmetry Methods, pp. 8 - - """ - - h = match['h'] - hy = match['hy'] - func = match['func'] - x = func.args[0] - y = match['y'] - xi = Function('xi')(x, func) - eta = Function('eta')(x, func) - - if h.is_rational_function(): - schi, schix, schiy = symbols("schi, schix, schiy") - cpde = schix + h*schiy - hy*schi - num, denom = cancel(cpde).as_numer_denom() - deg = Poly(num, x, y).total_degree() - - chi = Function('chi')(x, y) - chix = chi.diff(x) - chiy = chi.diff(y) - cpde = chix + h*chiy - hy*chi - chieq = Symbol("chi") - for i in range(1, deg + 1): - chieq += Add(*[ - Symbol("chi_" + str(power) + "_" + str(i - power))*x**power*y**(i - power) - for power in range(i + 1)]) - cnum, cden = cancel(cpde.subs({chi : chieq}).doit()).as_numer_denom() - cnum = expand(cnum) - if cnum.is_polynomial(x, y) and cnum.is_Add: - cpoly = Poly(cnum, x, y).as_dict() - if cpoly: - solsyms = chieq.free_symbols - {x, y} - soldict = solve(cpoly.values(), *solsyms) - if isinstance(soldict, list): - soldict = soldict[0] - if any(soldict.values()): - chieq = chieq.subs(soldict) - dict_ = {sym: 1 for sym in solsyms} - chieq = chieq.subs(dict_) - # After finding chi, the main aim is to find out - # eta, xi by the equation eta = xi*h + chi - # One method to set xi, would be rearranging it to - # (eta/h) - xi = (chi/h). This would mean dividing - # chi by h would give -xi as the quotient and eta - # as the remainder. Thanks to Sean Vig for suggesting - # this method. - xic, etac = div(chieq, h) - inf = {eta: etac.subs(y, func), xi: -xic.subs(y, func)} - return [inf] - -def lie_heuristic_function_sum(match, comp=False): - r""" - This heuristic uses the following two assumptions on `\xi` and `\eta` - - .. math:: \eta = 0, \xi = f(x) + g(y) - - .. math:: \eta = f(x) + g(y), \xi = 0 - - The first assumption of this heuristic holds good if - - .. math:: \frac{\partial}{\partial y}[(h\frac{\partial^{2}}{ - \partial x^{2}}(h^{-1}))^{-1}] - - is separable in `x` and `y`, - - 1. The separated factors containing `y` is `\frac{\partial g}{\partial y}`. - From this `g(y)` can be determined. - 2. The separated factors containing `x` is `f''(x)`. - 3. `h\frac{\partial^{2}}{\partial x^{2}}(h^{-1})` equals - `\frac{f''(x)}{f(x) + g(y)}`. From this `f(x)` can be determined. - - The second assumption holds good if `\frac{dy}{dx} = h(x, y)` is rewritten as - `\frac{dy}{dx} = \frac{1}{h(y, x)}` and the same properties of the first - assumption satisfies. After obtaining `f(x)` and `g(y)`, the coordinates - are again interchanged, to get `\eta` as `f(x) + g(y)`. - - For both assumptions, the constant factors are separated among `g(y)` - and `f''(x)`, such that `f''(x)` obtained from 3] is the same as that - obtained from 2]. If not possible, then this heuristic fails. - - - References - ========== - - E.S. Cheb-Terrab, A.D. Roche, Symmetries and First Order - ODE Patterns, pp. 7 - pp. 8 - - """ - - xieta = [] - h = match['h'] - func = match['func'] - hinv = match['hinv'] - x = func.args[0] - y = match['y'] - xi = Function('xi')(x, func) - eta = Function('eta')(x, func) - - for odefac in [h, hinv]: - factor = odefac*((1/odefac).diff(x, 2)) - sep = separatevars((1/factor).diff(y), dict=True, symbols=[x, y]) - if sep and sep['coeff'] and sep[x].has(x) and sep[y].has(y): - k = Dummy("k") - try: - gy = k*integrate(sep[y], y) - except NotImplementedError: - pass - else: - fdd = 1/(k*sep[x]*sep['coeff']) - fx = simplify(fdd/factor - gy) - check = simplify(fx.diff(x, 2) - fdd) - if fx: - if not check: - fx = fx.subs(k, 1) - gy = (gy/k) - else: - sol = solve(check, k) - if sol: - sol = sol[0] - fx = fx.subs(k, sol) - gy = (gy/k)*sol - else: - continue - if odefac == hinv: # Inverse ODE - fx = fx.subs(x, y) - gy = gy.subs(y, x) - etaval = factor_terms(fx + gy) - if etaval.is_Mul: - etaval = Mul(*[arg for arg in etaval.args if arg.has(x, y)]) - if odefac == hinv: # Inverse ODE - inf = {eta: etaval.subs(y, func), xi : S.Zero} - else: - inf = {xi: etaval.subs(y, func), eta : S.Zero} - if not comp: - return [inf] - else: - xieta.append(inf) - - if xieta: - return xieta - -def lie_heuristic_abaco2_similar(match, comp=False): - r""" - This heuristic uses the following two assumptions on `\xi` and `\eta` - - .. math:: \eta = g(x), \xi = f(x) - - .. math:: \eta = f(y), \xi = g(y) - - For the first assumption, - - 1. First `\frac{\frac{\partial h}{\partial y}}{\frac{\partial^{2} h}{ - \partial yy}}` is calculated. Let us say this value is A - - 2. If this is constant, then `h` is matched to the form `A(x) + B(x)e^{ - \frac{y}{C}}` then, `\frac{e^{\int \frac{A(x)}{C} \,dx}}{B(x)}` gives `f(x)` - and `A(x)*f(x)` gives `g(x)` - - 3. Otherwise `\frac{\frac{\partial A}{\partial X}}{\frac{\partial A}{ - \partial Y}} = \gamma` is calculated. If - - a] `\gamma` is a function of `x` alone - - b] `\frac{\gamma\frac{\partial h}{\partial y} - \gamma'(x) - \frac{ - \partial h}{\partial x}}{h + \gamma} = G` is a function of `x` alone. - then, `e^{\int G \,dx}` gives `f(x)` and `-\gamma*f(x)` gives `g(x)` - - The second assumption holds good if `\frac{dy}{dx} = h(x, y)` is rewritten as - `\frac{dy}{dx} = \frac{1}{h(y, x)}` and the same properties of the first assumption - satisfies. After obtaining `f(x)` and `g(x)`, the coordinates are again - interchanged, to get `\xi` as `f(x^*)` and `\eta` as `g(y^*)` - - References - ========== - - E.S. Cheb-Terrab, A.D. Roche, Symmetries and First Order - ODE Patterns, pp. 10 - pp. 12 - - """ - - h = match['h'] - hx = match['hx'] - hy = match['hy'] - func = match['func'] - hinv = match['hinv'] - x = func.args[0] - y = match['y'] - xi = Function('xi')(x, func) - eta = Function('eta')(x, func) - - factor = cancel(h.diff(y)/h.diff(y, 2)) - factorx = factor.diff(x) - factory = factor.diff(y) - if not factor.has(x) and not factor.has(y): - A = Wild('A', exclude=[y]) - B = Wild('B', exclude=[y]) - C = Wild('C', exclude=[x, y]) - match = h.match(A + B*exp(y/C)) - try: - tau = exp(-integrate(match[A]/match[C]), x)/match[B] - except NotImplementedError: - pass - else: - gx = match[A]*tau - return [{xi: tau, eta: gx}] - - else: - gamma = cancel(factorx/factory) - if not gamma.has(y): - tauint = cancel((gamma*hy - gamma.diff(x) - hx)/(h + gamma)) - if not tauint.has(y): - try: - tau = exp(integrate(tauint, x)) - except NotImplementedError: - pass - else: - gx = -tau*gamma - return [{xi: tau, eta: gx}] - - factor = cancel(hinv.diff(y)/hinv.diff(y, 2)) - factorx = factor.diff(x) - factory = factor.diff(y) - if not factor.has(x) and not factor.has(y): - A = Wild('A', exclude=[y]) - B = Wild('B', exclude=[y]) - C = Wild('C', exclude=[x, y]) - match = h.match(A + B*exp(y/C)) - try: - tau = exp(-integrate(match[A]/match[C]), x)/match[B] - except NotImplementedError: - pass - else: - gx = match[A]*tau - return [{eta: tau.subs(x, func), xi: gx.subs(x, func)}] - - else: - gamma = cancel(factorx/factory) - if not gamma.has(y): - tauint = cancel((gamma*hinv.diff(y) - gamma.diff(x) - hinv.diff(x))/( - hinv + gamma)) - if not tauint.has(y): - try: - tau = exp(integrate(tauint, x)) - except NotImplementedError: - pass - else: - gx = -tau*gamma - return [{eta: tau.subs(x, func), xi: gx.subs(x, func)}] - - -def lie_heuristic_abaco2_unique_unknown(match, comp=False): - r""" - This heuristic assumes the presence of unknown functions or known functions - with non-integer powers. - - 1. A list of all functions and non-integer powers containing x and y - 2. Loop over each element `f` in the list, find `\frac{\frac{\partial f}{\partial x}}{ - \frac{\partial f}{\partial x}} = R` - - If it is separable in `x` and `y`, let `X` be the factors containing `x`. Then - - a] Check if `\xi = X` and `\eta = -\frac{X}{R}` satisfy the PDE. If yes, then return - `\xi` and `\eta` - b] Check if `\xi = \frac{-R}{X}` and `\eta = -\frac{1}{X}` satisfy the PDE. - If yes, then return `\xi` and `\eta` - - If not, then check if - - a] :math:`\xi = -R,\eta = 1` - - b] :math:`\xi = 1, \eta = -\frac{1}{R}` - - are solutions. - - References - ========== - - E.S. Cheb-Terrab, A.D. Roche, Symmetries and First Order - ODE Patterns, pp. 10 - pp. 12 - - """ - - h = match['h'] - hx = match['hx'] - hy = match['hy'] - func = match['func'] - x = func.args[0] - y = match['y'] - xi = Function('xi')(x, func) - eta = Function('eta')(x, func) - - funclist = [] - for atom in h.atoms(Pow): - base, exp = atom.as_base_exp() - if base.has(x) and base.has(y): - if not exp.is_Integer: - funclist.append(atom) - - for function in h.atoms(AppliedUndef): - syms = function.free_symbols - if x in syms and y in syms: - funclist.append(function) - - for f in funclist: - frac = cancel(f.diff(y)/f.diff(x)) - sep = separatevars(frac, dict=True, symbols=[x, y]) - if sep and sep['coeff']: - xitry1 = sep[x] - etatry1 = -1/(sep[y]*sep['coeff']) - pde1 = etatry1.diff(y)*h - xitry1.diff(x)*h - xitry1*hx - etatry1*hy - if not simplify(pde1): - return [{xi: xitry1, eta: etatry1.subs(y, func)}] - xitry2 = 1/etatry1 - etatry2 = 1/xitry1 - pde2 = etatry2.diff(x) - (xitry2.diff(y))*h**2 - xitry2*hx - etatry2*hy - if not simplify(expand(pde2)): - return [{xi: xitry2.subs(y, func), eta: etatry2}] - - else: - etatry = -1/frac - pde = etatry.diff(x) + etatry.diff(y)*h - hx - etatry*hy - if not simplify(pde): - return [{xi: S.One, eta: etatry.subs(y, func)}] - xitry = -frac - pde = -xitry.diff(x)*h -xitry.diff(y)*h**2 - xitry*hx -hy - if not simplify(expand(pde)): - return [{xi: xitry.subs(y, func), eta: S.One}] - - -def lie_heuristic_abaco2_unique_general(match, comp=False): - r""" - This heuristic finds if infinitesimals of the form `\eta = f(x)`, `\xi = g(y)` - without making any assumptions on `h`. - - The complete sequence of steps is given in the paper mentioned below. - - References - ========== - - E.S. Cheb-Terrab, A.D. Roche, Symmetries and First Order - ODE Patterns, pp. 10 - pp. 12 - - """ - hx = match['hx'] - hy = match['hy'] - func = match['func'] - x = func.args[0] - y = match['y'] - xi = Function('xi')(x, func) - eta = Function('eta')(x, func) - - A = hx.diff(y) - B = hy.diff(y) + hy**2 - C = hx.diff(x) - hx**2 - - if not (A and B and C): - return - - Ax = A.diff(x) - Ay = A.diff(y) - Axy = Ax.diff(y) - Axx = Ax.diff(x) - Ayy = Ay.diff(y) - D = simplify(2*Axy + hx*Ay - Ax*hy + (hx*hy + 2*A)*A)*A - 3*Ax*Ay - if not D: - E1 = simplify(3*Ax**2 + ((hx**2 + 2*C)*A - 2*Axx)*A) - if E1: - E2 = simplify((2*Ayy + (2*B - hy**2)*A)*A - 3*Ay**2) - if not E2: - E3 = simplify( - E1*((28*Ax + 4*hx*A)*A**3 - E1*(hy*A + Ay)) - E1.diff(x)*8*A**4) - if not E3: - etaval = cancel((4*A**3*(Ax - hx*A) + E1*(hy*A - Ay))/(S(2)*A*E1)) - if x not in etaval: - try: - etaval = exp(integrate(etaval, y)) - except NotImplementedError: - pass - else: - xival = -4*A**3*etaval/E1 - if y not in xival: - return [{xi: xival, eta: etaval.subs(y, func)}] - - else: - E1 = simplify((2*Ayy + (2*B - hy**2)*A)*A - 3*Ay**2) - if E1: - E2 = simplify( - 4*A**3*D - D**2 + E1*((2*Axx - (hx**2 + 2*C)*A)*A - 3*Ax**2)) - if not E2: - E3 = simplify( - -(A*D)*E1.diff(y) + ((E1.diff(x) - hy*D)*A + 3*Ay*D + - (A*hx - 3*Ax)*E1)*E1) - if not E3: - etaval = cancel(((A*hx - Ax)*E1 - (Ay + A*hy)*D)/(S(2)*A*D)) - if x not in etaval: - try: - etaval = exp(integrate(etaval, y)) - except NotImplementedError: - pass - else: - xival = -E1*etaval/D - if y not in xival: - return [{xi: xival, eta: etaval.subs(y, func)}] - - -def lie_heuristic_linear(match, comp=False): - r""" - This heuristic assumes - - 1. `\xi = ax + by + c` and - 2. `\eta = fx + gy + h` - - After substituting the following assumptions in the determining PDE, it - reduces to - - .. math:: f + (g - a)h - bh^{2} - (ax + by + c)\frac{\partial h}{\partial x} - - (fx + gy + c)\frac{\partial h}{\partial y} - - Solving the reduced PDE obtained, using the method of characteristics, becomes - impractical. The method followed is grouping similar terms and solving the system - of linear equations obtained. The difference between the bivariate heuristic is that - `h` need not be a rational function in this case. - - References - ========== - - E.S. Cheb-Terrab, A.D. Roche, Symmetries and First Order - ODE Patterns, pp. 10 - pp. 12 - - """ - h = match['h'] - hx = match['hx'] - hy = match['hy'] - func = match['func'] - x = func.args[0] - y = match['y'] - xi = Function('xi')(x, func) - eta = Function('eta')(x, func) - - coeffdict = {} - symbols = numbered_symbols("c", cls=Dummy) - symlist = [next(symbols) for _ in islice(symbols, 6)] - C0, C1, C2, C3, C4, C5 = symlist - pde = C3 + (C4 - C0)*h - (C0*x + C1*y + C2)*hx - (C3*x + C4*y + C5)*hy - C1*h**2 - pde, denom = pde.as_numer_denom() - pde = powsimp(expand(pde)) - if pde.is_Add: - terms = pde.args - for term in terms: - if term.is_Mul: - rem = Mul(*[m for m in term.args if not m.has(x, y)]) - xypart = term/rem - if xypart not in coeffdict: - coeffdict[xypart] = rem - else: - coeffdict[xypart] += rem - else: - if term not in coeffdict: - coeffdict[term] = S.One - else: - coeffdict[term] += S.One - - sollist = coeffdict.values() - soldict = solve(sollist, symlist) - if soldict: - if isinstance(soldict, list): - soldict = soldict[0] - subval = soldict.values() - if any(t for t in subval): - onedict = dict(zip(symlist, [1]*6)) - xival = C0*x + C1*func + C2 - etaval = C3*x + C4*func + C5 - xival = xival.subs(soldict) - etaval = etaval.subs(soldict) - xival = xival.subs(onedict) - etaval = etaval.subs(onedict) - return [{xi: xival, eta: etaval}] - def sysode_linear_2eq_order1(match_): x = match_['func'][0].func @@ -4772,5 +3588,5 @@ def _nonlinear_3eq_order1_type5(x, y, z, t, eq): HomogeneousCoeffBest, LinearCoefficients, NthOrderReducible, SecondHypergeometric, NthLinearConstantCoeffHomogeneous, NthLinearConstantCoeffVariationOfParameters, NthLinearConstantCoeffUndeterminedCoefficients, NthLinearEulerEqHomogeneous, - NthLinearEulerEqNonhomogeneousVariationOfParameters, + NthLinearEulerEqNonhomogeneousVariationOfParameters, LieGroup, NthLinearEulerEqNonhomogeneousUndeterminedCoefficients, SecondLinearBessel, SecondLinearAiry) diff --git a/sympy/solvers/ode/single.py b/sympy/solvers/ode/single.py index efdc382e530c..45e0421600c7 100644 --- a/sympy/solvers/ode/single.py +++ b/sympy/solvers/ode/single.py @@ -32,6 +32,7 @@ from .nonhomogeneous import _get_euler_characteristic_eq_sols, _get_const_characteristic_eq_sols, \ _solve_undetermined_coefficients, _solve_variation_of_parameters, _test_term, _undetermined_coefficients_match, \ _get_simplified_sol +from .lie_group import _ode_lie_group class ODEMatchError(NotImplementedError): @@ -90,7 +91,7 @@ class SingleODEProblem: _eq_preprocessed = None # type: Expr _eq_high_order_free = None - def __init__(self, eq, func, sym, prep=True): + def __init__(self, eq, func, sym, prep=True, **kwargs): assert isinstance(eq, Expr) assert isinstance(func, AppliedUndef) assert isinstance(sym, Symbol) @@ -99,6 +100,7 @@ def __init__(self, eq, func, sym, prep=True): self.func = func self.sym = sym self.prep = prep + self.params = kwargs @cached_property def order(self) -> int: @@ -2793,6 +2795,109 @@ def _get_general_solution(self, *, simplify_flag: bool = True): return [Eq(f(x), C1*airyai(arg) + C2*airybi(arg))] +class LieGroup(SingleODESolver): + r""" + This hint implements the Lie group method of solving first order differential + equations. The aim is to convert the given differential equation from the + given coordinate system into another coordinate system where it becomes + invariant under the one-parameter Lie group of translations. The converted + ODE can be easily solved by quadrature. It makes use of the + :py:meth:`sympy.solvers.ode.infinitesimals` function which returns the + infinitesimals of the transformation. + + The coordinates `r` and `s` can be found by solving the following Partial + Differential Equations. + + .. math :: \xi\frac{\partial r}{\partial x} + \eta\frac{\partial r}{\partial y} + = 0 + + .. math :: \xi\frac{\partial s}{\partial x} + \eta\frac{\partial s}{\partial y} + = 1 + + The differential equation becomes separable in the new coordinate system + + .. math :: \frac{ds}{dr} = \frac{\frac{\partial s}{\partial x} + + h(x, y)\frac{\partial s}{\partial y}}{ + \frac{\partial r}{\partial x} + h(x, y)\frac{\partial r}{\partial y}} + + After finding the solution by integration, it is then converted back to the original + coordinate system by substituting `r` and `s` in terms of `x` and `y` again. + + Examples + ======== + + >>> from sympy import Function, dsolve, exp, pprint + >>> from sympy.abc import x + >>> f = Function('f') + >>> pprint(dsolve(f(x).diff(x) + 2*x*f(x) - x*exp(-x**2), f(x), + ... hint='lie_group')) + / 2\ 2 + | x | -x + f(x) = |C1 + --|*e + \ 2 / + + + References + ========== + + - Solving differential equations by Symmetry Groups, + John Starrett, pp. 1 - pp. 14 + + """ + hint = "lie_group" + has_integral = False + + def _has_additional_params(self): + return 'xi' in self.ode_problem.params and 'eta' in self.ode_problem.params + + def _matches(self): + eq = self.ode_problem.eq + f = self.ode_problem.func.func + order = self.ode_problem.order + x = self.ode_problem.sym + df = f(x).diff(x) + y = Dummy('y') + d = Wild('d', exclude=[df, f(x).diff(x, 2)]) + e = Wild('e', exclude=[df]) + does_match = False + if self._has_additional_params() and order == 1: + xi = self.ode_problem.params['xi'] + eta = self.ode_problem.params['eta'] + self.r3 = {'xi': xi, 'eta': eta} + r = collect(eq, df, exact=True).match(d + e * df) + if r: + r['d'] = d + r['e'] = e + r['y'] = y + r[d] = r[d].subs(f(x), y) + r[e] = r[e].subs(f(x), y) + self.r3.update(r) + does_match = True + return does_match + + def _get_general_solution(self, *, simplify_flag: bool = True): + eq = self.ode_problem.eq + x = self.ode_problem.sym + func = self.ode_problem.func + order = self.ode_problem.order + df = func.diff(x) + + try: + eqsol = solve(eq, df) + except NotImplementedError: + eqsol = [] + + desols = [] + for s in eqsol: + sol = _ode_lie_group(s, func, order, match=self.r3) + if sol: + desols.extend(sol) + + if desols == []: + raise NotImplementedError("The given ODE " + str(eq) + " cannot be solved by" + + " the lie group method") + return desols + # Avoid circular import: from .ode import dsolve, ode_sol_simplicity, odesimp, homogeneous_order diff --git a/sympy/solvers/ode/tests/test_lie_group.py b/sympy/solvers/ode/tests/test_lie_group.py index 07c2cd651173..0f6c37e1210a 100644 --- a/sympy/solvers/ode/tests/test_lie_group.py +++ b/sympy/solvers/ode/tests/test_lie_group.py @@ -1,7 +1,7 @@ from sympy import (atan, Eq, exp, Function, log, Rational, sin, sqrt, Symbol, tan, symbols) -from sympy.solvers.ode import (classify_ode, infinitesimals, checkinfsol, dsolve) +from sympy.solvers.ode import (classify_ode, checkinfsol, dsolve, infinitesimals) from sympy.solvers.ode.subscheck import checkodesol