Merge pull request #7290 from KDMOE/nonhomogeneous_euler_eq

ODE: Added non homogeneoous linear euler eq fixes #7172

sympy/solvers/ode.py

@@ -292,7 +292,9 @@
     "nth_linear_constant_coeff_homogeneous",
     "nth_linear_euler_eq_homogeneous",
     "nth_linear_constant_coeff_undetermined_coefficients",
+    "nth_linear_euler_eq_nonhomogeneous_undetermined_coefficients",
     "nth_linear_constant_coeff_variation_of_parameters",
+    "nth_linear_euler_eq_nonhomogeneous_variation_of_parameters",
     "Liouville",
     "2nd_power_series_ordinary",
     "2nd_power_series_regular",
@@ -306,6 +308,7 @@
     "linear_coefficients_Integral",
     "separable_reduced_Integral",
     "nth_linear_constant_coeff_variation_of_parameters_Integral",
+    "nth_linear_euler_eq_nonhomogeneous_variation_of_parameters_Integral",
     "Liouville_Integral",
     )
 
@@ -1121,10 +1124,11 @@ class in it.  Note that a hint may do this anyway if
             else:
                 matching_hints["nth_linear_constant_coeff_homogeneous"] = r
 
-        # Euler equation case (a_i * x**i for all i)
+        # nth order Euler equation a_n*x**n*y^(n) + ... + a_1*x*y' + a_0*y = F(x)
+        #In case of Homogeneous euler equation F(x) = 0
         def _test_term(coeff, order):
             r"""
-            Linear Euler ODEs have the form  K*x**order*diff(y(x),x,order),
+            Linear Euler ODEs have the form  K*x**order*diff(y(x),x,order) = F(x),
             where K is independent of x and y(x), order>= 0.
             So we need to check that for each term, coeff == K*x**order from
             some K.  We have a few cases, since coeff may have several
@@ -1137,7 +1141,7 @@ def _test_term(coeff, order):
             if order == 0:
                 if x in coeff.free_symbols:
                     return False
-                return f(x) not in coeff.atoms()
+                return True
             if coeff.is_Mul:
                 if coeff.has(f(x)):
                     return False
@@ -1150,6 +1154,15 @@ def _test_term(coeff, order):
         if r and not any(not _test_term(r[i], i) for i in r if i >= 0):
             if not r[-1]:
                 matching_hints["nth_linear_euler_eq_homogeneous"] = r
+            else:
+                matching_hints["nth_linear_euler_eq_nonhomogeneous_variation_of_parameters"] = r
+                matching_hints["nth_linear_euler_eq_nonhomogeneous_variation_of_parameters_Integral"] = r
+                e, re = posify(r[-1].subs(x, exp(x)))
+                undetcoeff = _undetermined_coefficients_match(e.subs(re), x)
+                if undetcoeff['test']:
+                    r['trialset'] = undetcoeff['trialset']
+                    matching_hints["nth_linear_euler_eq_nonhomogeneous_undetermined_coefficients"] = r
+
 
     # Order keys based on allhints.
     retlist = [i for i in allhints if i in matching_hints]
@@ -1227,6 +1240,8 @@ def odesimp(eq, func, order, hint):
 
     # First, integrate if the hint allows it.
     eq = _handle_Integral(eq, func, order, hint)
+    if hint.startswith("nth_linear_euler_eq_nonhomogeneous"):
+        eq = simplify(eq)
     if not isinstance(eq, Equality):
         raise TypeError("eq should be an instance of Equality")
 
@@ -3235,6 +3250,141 @@ def ode_nth_linear_euler_eq_homogeneous(eq, func, order, match, returns='sol'):
     else:
         raise ValueError('Unknown value for key "returns".')
 
+
+def ode_nth_linear_euler_eq_nonhomogeneous_undetermined_coefficients(eq, func, order, match, returns='sol'):
+    r"""
+    Solves an `n`\th order linear non homogeneous Cauchy-Euler equidimensional
+    ordinary differential equation using undetermined coefficients.
+
+    This is an equation with form `g(x) = a_0 f(x) + a_1 x f'(x) + a_2 x^2 f''(x)
+    \cdots`.
+
+    These equations can be solved in a general manner, by substituting
+    solutions of the form `x = exp(t)`, and deriving a characteristic equation
+    of form `g(exp(t)) = b_0 f(t) + b_1 f'(t) + b_2 f''(t) \cdots` which can
+    be then solved by nth_linear_constant_coeff_undetermined_coefficients if
+    g(exp(t)) has finite number of lineary independent derivatives.
+
+    Functions that fit this requirement are finite sums functions of the form
+    `a x^i e^{b x} \sin(c x + d)` or `a x^i e^{b x} \cos(c x + d)`, where `i`
+    is a non-negative integer and `a`, `b`, `c`, and `d` are constants.  For
+    example any polynomial in `x`, functions like `x^2 e^{2 x}`, `x \sin(x)`,
+    and `e^x \cos(x)` can all be used.  Products of `\sin`'s and `\cos`'s have
+    a finite number of derivatives, because they can be expanded into `\sin(a
+    x)` and `\cos(b x)` terms.  However, SymPy currently cannot do that
+    expansion, so you will need to manually rewrite the expression in terms of
+    the above to use this method.  So, for example, you will need to manually
+    convert `\sin^2(x)` into `(1 + \cos(2 x))/2` to properly apply the method
+    of undetermined coefficients on it.
+
+    After replacement of x by exp(t), this method works by creating a trial function
+    from the expression and all of its linear independent derivatives and
+    substituting them into the original ODE.  The coefficients for each term
+    will be a system of linear equations, which are be solved for and
+    substituted, giving the solution. If any of the trial functions are linearly
+    dependent on the solution to the homogeneous equation, they are multiplied
+    by sufficient `x` to make them linearly independent.
+
+    Examples
+    ========
+
+    >>> from sympy import dsolve, Function, Derivative, log
+    >>> from sympy.abc import x
+    >>> f = Function('f')
+    >>> eq = x**2*Derivative(f(x), x, x) - 2*x*Derivative(f(x), x) + 2*f(x) - log(x)
+    >>> dsolve(eq, f(x),
+    ... hint='nth_linear_euler_eq_nonhomogeneous_undetermined_coefficients').expand()
+    f(x) == C1*x + C2*x**2 + log(x)/2 + 3/4
+
+    """
+    x = func.args[0]
+    f = func.func
+    r = match
+
+    chareq, eq, symbol = S.Zero, S.Zero, Dummy('x')
+
+    for i in r.keys():
+        if not isinstance(i, str) and i >= 0:
+            chareq += (r[i]*diff(x**symbol, x, i)*x**-symbol).expand()
+
+    for i in range(1,degree(Poly(chareq, symbol))+1):
+        eq += chareq.coeff(symbol**i)*diff(f(x), x, i)
+
+    if chareq.as_coeff_add(symbol)[0]:
+        eq += chareq.as_coeff_add(symbol)[0]*f(x)
+    e, re = posify(r[-1].subs(x, exp(x)))
+    eq += e.subs(re)
+
+    match = _nth_linear_match(eq, f(x), ode_order(eq, f(x)))
+    match['trialset'] = r['trialset']
+    return ode_nth_linear_constant_coeff_undetermined_coefficients(eq, func, order, match).subs(x, log(x)).subs(f(log(x)), f(x)).expand()
+
+
+def ode_nth_linear_euler_eq_nonhomogeneous_variation_of_parameters(eq, func, order, match, returns='sol'):
+    r"""
+    Solves an `n`\th order linear non homogeneous Cauchy-Euler equidimensional
+    ordinary differential equation using variation of parameters.
+
+    This is an equation with form `g(x) = a_0 f(x) + a_1 x f'(x) + a_2 x^2 f''(x)
+    \cdots`.
+
+    This method works by assuming that the particular solution takes the form
+
+    .. math:: \sum_{x=1}^{n} c_i(x) y_i(x) {a_n} {x^n} \text{,}
+
+    where `y_i` is the `i`\th solution to the homogeneous equation.  The
+    solution is then solved using Wronskian's and Cramer's Rule.  The
+    particular solution is given by multiplying eq given below with `a_n x^{n}`
+
+    .. math:: \sum_{x=1}^n \left( \int \frac{W_i(x)}{W(x)} \,dx
+                \right) y_i(x) \text{,}
+
+    where `W(x)` is the Wronskian of the fundamental system (the system of `n`
+    linearly independent solutions to the homogeneous equation), and `W_i(x)`
+    is the Wronskian of the fundamental system with the `i`\th column replaced
+    with `[0, 0, \cdots, 0, \frac{x^{- n}}{a_n} g{\left (x \right )}]`.
+
+    This method is general enough to solve any `n`\th order inhomogeneous
+    linear differential equation, but sometimes SymPy cannot simplify the
+    Wronskian well enough to integrate it.  If this method hangs, try using the
+    ``nth_linear_constant_coeff_variation_of_parameters_Integral`` hint and
+    simplifying the integrals manually.  Also, prefer using
+    ``nth_linear_constant_coeff_undetermined_coefficients`` when it
+    applies, because it doesn't use integration, making it faster and more
+    reliable.
+
+    Warning, using simplify=False with
+    'nth_linear_constant_coeff_variation_of_parameters' in
+    :py:meth:`~sympy.solvers.ode.dsolve` may cause it to hang, because it will
+    not attempt to simplify the Wronskian before integrating.  It is
+    recommended that you only use simplify=False with
+    'nth_linear_constant_coeff_variation_of_parameters_Integral' for this
+    method, especially if the solution to the homogeneous equation has
+    trigonometric functions in it.
+
+    Examples
+    ========
+
+    >>> from sympy import Function, dsolve, Derivative
+    >>> from sympy.abc import x
+    >>> f = Function('f')
+    >>> eq = x**2*Derivative(f(x), x, x) - 2*x*Derivative(f(x), x) + 2*f(x) - x**4
+    >>> dsolve(eq, f(x),
+    ... hint='nth_linear_euler_eq_nonhomogeneous_variation_of_parameters').expand()
+    f(x) == C1*x + C2*x**2 + x**4/6
+
+    """
+    x = func.args[0]
+    f = func.func
+    r = match
+
+    gensol = ode_nth_linear_euler_eq_homogeneous(eq, func, order, match, returns='both')
+    match.update(gensol)
+    r[-1] = r[-1]/r[ode_order(eq, f(x))]
+    sol = _solve_variation_of_parameters(eq, func, order, match)
+    return Eq(f(x), r['sol'].rhs + (sol.rhs - r['sol'].rhs)*r[ode_order(eq, f(x))])
+
+
 def ode_almost_linear(eq, func, order, match):
     r"""
     Solves an almost-linear differential equation.
@@ -4130,7 +4280,7 @@ def ode_nth_linear_constant_coeff_variation_of_parameters(eq, func, order, match
 
 def _solve_variation_of_parameters(eq, func, order, match):
     r"""
-    Helper function for the method of variation of parameters.
+    Helper function for the method of variation of parameters and nonhomogeneous euler eq.
 
     See the
     :py:meth:`~sympy.solvers.ode.ode_nth_linear_constant_coeff_variation_of_parameters`

sympy/solvers/tests/test_ode.py

@@ -1396,6 +1396,93 @@ def test_nth_order_linear_euler_eq_homogeneous():
     assert checkodesol(eq, sol, order=2, solve_for_func=False)[0]
 
 
+def test_nth_order_linear_euler_eq_nonhomogeneous_undetermined_coefficients():
+    x, t = symbols('x t')
+    a, b, c, d = symbols('a b c d', integer=True)
+    our_hint = "nth_linear_euler_eq_nonhomogeneous_undetermined_coefficients"
+
+    eq = x**4*diff(f(x), x, 4) - 13*x**2*diff(f(x), x, 2) + 36*f(x)	+ x
+    assert our_hint in classify_ode(eq, f(x))
+
+    eq = a*x**2*diff(f(x), x, 2) + b*x*diff(f(x), x) + c*f(x) + d*log(x)
+    assert our_hint in classify_ode(eq, f(x))
+
+    eq = Eq(x**2*diff(f(x), x, x) + x*diff(f(x), x), 1)
+    sol =  C1 + C2*log(x) + log(x)**2/2
+    sols = constant_renumber(sol, 'C', 1, 2)
+    assert our_hint in classify_ode(eq, f(x))
+    assert dsolve(eq, f(x), hint=our_hint).rhs in (sol, sols)
+    assert checkodesol(eq, sol, order=2, solve_for_func=False)[0]
+
+    eq = Eq(x**2*diff(f(x), x, x) - 2*x*diff(f(x), x) + 2*f(x), x**3)
+    sol = x*(C1 + C2*x + Rational(1, 2)*x**2)
+    sols = constant_renumber(sol, 'C', 1, 2)
+    assert our_hint in classify_ode(eq, f(x))
+    assert dsolve(eq, f(x), hint=our_hint).rhs in (sol, sols)
+    assert checkodesol(eq, sol, order=2, solve_for_func=False)[0]
+
+    eq = Eq(x**2*diff(f(x), x, x) - x*diff(f(x), x) - 3*f(x), log(x)/x)
+    sol =  C1/x + C2*x**3 - Rational(1, 16)*log(x)/x - Rational(1, 8)*log(x)**2/x
+    sols = constant_renumber(sol, 'C', 1, 2)
+    assert our_hint in classify_ode(eq, f(x))
+    assert dsolve(eq, f(x), hint=our_hint).rhs.expand() in (sol, sols)
+    assert checkodesol(eq, sol, order=2, solve_for_func=False)[0]
+
+    eq = Eq(x**2*diff(f(x), x, x) + 3*x*diff(f(x), x) - 8*f(x), log(x)**3 - log(x))
+    sol = C1/x**4 + C2*x**2 - Rational(1,8)*log(x)**3 - Rational(3,32)*log(x)**2 - Rational(1,64)*log(x) - Rational(7, 256)
+    sols = constant_renumber(sol, 'C', 1, 2)
+    assert our_hint in classify_ode(eq)
+    assert dsolve(eq, f(x), hint=our_hint).rhs.expand() in (sol, sols)
+    assert checkodesol(eq, sol, order=2, solve_for_func=False)[0]
+
+    eq = Eq(x**3*diff(f(x), x, x, x) - 3*x**2*diff(f(x), x, x) + 6*x*diff(f(x), x) - 6*f(x), log(x))
+    sol = C1*x + C2*x**2 + C3*x**3 - Rational(1, 6)*log(x) - Rational(11, 36)
+    sols = constant_renumber(sol, 'C', 1, 3)
+    assert our_hint in classify_ode(eq)
+    assert dsolve(eq, f(x), hint=our_hint).rhs.expand() in (sol, sols)
+    assert checkodesol(eq, sol, order=2, solve_for_func=False)[0]
+
+
+def test_nth_order_linear_euler_eq_nonhomogeneous_variation_of_parameters():
+    x, t = symbols('x, t')
+    a, b, c, d = symbols('a, b, c, d', integer=True)
+    our_hint = "nth_linear_euler_eq_nonhomogeneous_variation_of_parameters"
+
+    eq = Eq(x**2*diff(f(x),x,2) - 8*x*diff(f(x),x) + 12*f(x), x**2)
+    assert our_hint in classify_ode(eq, f(x))
+
+    eq = Eq(a*x**3*diff(f(x),x,3) + b*x**2*diff(f(x),x,2) + c*x*diff(f(x),x) + d*f(x), x*log(x))
+    assert our_hint in classify_ode(eq, f(x))
+
+    eq = Eq(x**2*Derivative(f(x), x, x) - 2*x*Derivative(f(x), x) + 2*f(x), x**4)
+    sol = C1*x + C2*x**2 + x**4/6
+    sols = constant_renumber(sol, 'C', 1, 2)
+    assert our_hint in classify_ode(eq)
+    assert dsolve(eq, f(x), hint=our_hint).rhs.expand() in (sol, sols)
+    assert checkodesol(eq, sol, order=2, solve_for_func=False)[0]
+
+    eq = Eq(3*x**2*diff(f(x), x, x) + 6*x*diff(f(x), x) - 6*f(x), x**3*exp(x))
+    sol = C1/x**2 + C2*x + x*exp(x)/3 - 4*exp(x)/3 + 8*exp(x)/(3*x) - 8*exp(x)/(3*x**2)
+    sols = constant_renumber(sol, 'C', 1, 2)
+    assert our_hint in classify_ode(eq)
+    assert dsolve(eq, f(x), hint=our_hint).rhs.expand() in (sol, sols)
+    assert checkodesol(eq, sol, order=2, solve_for_func=False)[0]
+
+    eq = Eq(x**2*Derivative(f(x), x, x) - 2*x*Derivative(f(x), x) + 2*f(x), x**4*exp(x))
+    sol = C1*x + C2*x**2 + x**2*exp(x) - 2*x*exp(x)
+    sols = constant_renumber(sol, 'C', 1, 2)
+    assert our_hint in classify_ode(eq)
+    assert dsolve(eq, f(x), hint=our_hint).rhs.expand() in (sol, sols)
+    assert checkodesol(eq, sol, order=2, solve_for_func=False)[0]
+
+    eq = x**2*Derivative(f(x), x, x) - 2*x*Derivative(f(x), x) + 2*f(x) - log(x)
+    sol = C1*x + C2*x**2 + log(x)/2 + 3/4
+    sols = constant_renumber(sol, 'C', 1, 2)
+    assert our_hint in classify_ode(eq)
+    assert dsolve(eq, f(x), hint=our_hint).rhs in (sol, sols)
+    assert checkodesol(eq, sol, order=2, solve_for_func=False)[0]
+
+
 def test_issue_5095():
     f = Function('f')
     raises(ValueError, lambda: dsolve(f(x).diff(x)**2, f(x), 'separable'))