ODE

User Functions

These are functions that are imported into the global namespace with from sympy import *. These functions (unlike Hint Functions, below) are intended for use by ordinary users of SymPy.

dsolve()

sympy.solvers.ode.dsolve(eq, func=None, hint='default', simplify=True, ics=None, xi=None, eta=None, x0=0, n=6, **kwargs)[source]

Solves any (supported) kind of ordinary differential equation and system of ordinary differential equations.

Examples

>>> from sympy import Function, dsolve, Eq, Derivative, sin, cos, symbols
>>> from sympy.abc import x
>>> f = Function('f')
>>> dsolve(Derivative(f(x), x, x) + 9*f(x), f(x))
Eq(f(x), C1*sin(3*x) + C2*cos(3*x))
>>> eq = sin(x)*cos(f(x)) + cos(x)*sin(f(x))*f(x).diff(x)
>>> dsolve(eq, hint='1st_exact')
[Eq(f(x), -acos(C1/cos(x)) + 2*pi), Eq(f(x), acos(C1/cos(x)))]
>>> dsolve(eq, hint='almost_linear')
[Eq(f(x), -acos(C1/cos(x)) + 2*pi), Eq(f(x), acos(C1/cos(x)))]
>>> t = symbols('t')
>>> x, y = symbols('x, y', function=True)
>>> eq = (Eq(Derivative(x(t),t), 12*t*x(t) + 8*y(t)), Eq(Derivative(y(t),t), 21*x(t) + 7*t*y(t)))
>>> dsolve(eq)
[Eq(x(t), C1*x0 + C2*x0*Integral(8*exp(Integral(7*t, t))*exp(Integral(12*t, t))/x0**2, t)),
Eq(y(t), C1*y0 + C2(y0*Integral(8*exp(Integral(7*t, t))*exp(Integral(12*t, t))/x0**2, t) +
exp(Integral(7*t, t))*exp(Integral(12*t, t))/x0))]
>>> eq = (Eq(Derivative(x(t),t),x(t)*y(t)*sin(t)), Eq(Derivative(y(t),t),y(t)**2*sin(t)))
>>> dsolve(eq)
{Eq(x(t), -exp(C1)/(C2*exp(C1) - cos(t))), Eq(y(t), -1/(C1 - cos(t)))}

For Single Ordinary Differential Equation

It is classified under this when number of equation in eq is one. Usage

dsolve(eq, f(x), hint) -> Solve ordinary differential equation eq for function f(x), using method hint.

Details

eq can be any supported ordinary differential equation (see the

ode docstring for supported methods). This can either be an Equality, or an expression, which is assumed to be equal to 0.

f(x) is a function of one variable whose derivatives in that

variable make up the ordinary differential equation eq. In many cases it is not necessary to provide this; it will be autodetected (and an error raised if it couldn’t be detected).

hint is the solving method that you want dsolve to use. Use

classify_ode(eq, f(x)) to get all of the possible hints for an ODE. The default hint, default, will use whatever hint is returned first by classify_ode(). See Hints below for more options that you can use for hint.

simplify enables simplification by

odesimp(). See its docstring for more information. Turn this off, for example, to disable solving of solutions for func or simplification of arbitrary constants. It will still integrate with this hint. Note that the solution may contain more arbitrary constants than the order of the ODE with this option enabled.

xi and eta are the infinitesimal functions of an ordinary

differential equation. They are the infinitesimals of the Lie group of point transformations for which the differential equation is invariant. The user can specify values for the infinitesimals. If nothing is specified, xi and eta are calculated using infinitesimals() with the help of various heuristics.

ics is the set of initial/boundary conditions for the differential equation.

It should be given in the form of {f(x0): x1, f(x).diff(x).subs(x, x2): x3} and so on. For power series solutions, if no initial conditions are specified f(0) is assumed to be C0 and the power series solution is calculated about 0.

x0 is the point about which the power series solution of a differential

equation is to be evaluated.

n gives the exponent of the dependent variable up to which the power series

solution of a differential equation is to be evaluated.

Hints

Aside from the various solving methods, there are also some meta-hints that you can pass to dsolve():

default:

This uses whatever hint is returned first by classify_ode(). This is the default argument to dsolve().

all:

To make dsolve() apply all relevant classification hints, use dsolve(ODE, func, hint="all"). This will return a dictionary of hint:solution terms. If a hint causes dsolve to raise the NotImplementedError, value of that hint’s key will be the exception object raised. The dictionary will also include some special keys:

  • order: The order of the ODE. See also ode_order() in deutils.py.

  • best: The simplest hint; what would be returned by best below.

  • best_hint: The hint that would produce the solution given by best. If more than one hint produces the best solution, the first one in the tuple returned by classify_ode() is chosen.

  • default: The solution that would be returned by default. This is the one produced by the hint that appears first in the tuple returned by classify_ode().

all_Integral:

This is the same as all, except if a hint also has a corresponding _Integral hint, it only returns the _Integral hint. This is useful if all causes dsolve() to hang because of a difficult or impossible integral. This meta-hint will also be much faster than all, because integrate() is an expensive routine.

best:

To have dsolve() try all methods and return the simplest one. This takes into account whether the solution is solvable in the function, whether it contains any Integral classes (i.e. unevaluatable integrals), and which one is the shortest in size.

See also the classify_ode() docstring for more info on hints, and the ode docstring for a list of all supported hints.

Tips

  • You can declare the derivative of an unknown function this way:

    >>> from sympy import Function, Derivative
    >>> from sympy.abc import x # x is the independent variable
    >>> f = Function("f")(x) # f is a function of x
    >>> # f_ will be the derivative of f with respect to x
    >>> f_ = Derivative(f, x)
    
  • See test_ode.py for many tests, which serves also as a set of examples for how to use dsolve().

  • dsolve() always returns an Equality class (except for the case when the hint is all or all_Integral). If possible, it solves the solution explicitly for the function being solved for. Otherwise, it returns an implicit solution.

  • Arbitrary constants are symbols named C1, C2, and so on.

  • Because all solutions should be mathematically equivalent, some hints may return the exact same result for an ODE. Often, though, two different hints will return the same solution formatted differently. The two should be equivalent. Also note that sometimes the values of the arbitrary constants in two different solutions may not be the same, because one constant may have “absorbed” other constants into it.

  • Do help(ode.ode_<hintname>) to get help more information on a specific hint, where <hintname> is the name of a hint without _Integral.

For System Of Ordinary Differential Equations

Usage

dsolve(eq, func) -> Solve a system of ordinary differential equations eq for func being list of functions including \(x(t)\), \(y(t)\), \(z(t)\) where number of functions in the list depends upon the number of equations provided in eq.

Details

eq can be any supported system of ordinary differential equations This can either be an Equality, or an expression, which is assumed to be equal to 0.

func holds x(t) and y(t) being functions of one variable which together with some of their derivatives make up the system of ordinary differential equation eq. It is not necessary to provide this; it will be autodetected (and an error raised if it couldn’t be detected).

Hints

The hints are formed by parameters returned by classify_sysode, combining them give hints name used later for forming method name.

classify_ode()

sympy.solvers.ode.classify_ode(eq, func=None, dict=False, ics=None, **kwargs)[source]

Returns a tuple of possible dsolve() classifications for an ODE.

The tuple is ordered so that first item is the classification that dsolve() uses to solve the ODE by default. In general, classifications at the near the beginning of the list will produce better solutions faster than those near the end, thought there are always exceptions. To make dsolve() use a different classification, use dsolve(ODE, func, hint=<classification>). See also the dsolve() docstring for different meta-hints you can use.

If dict is true, classify_ode() will return a dictionary of hint:match expression terms. This is intended for internal use by dsolve(). Note that because dictionaries are ordered arbitrarily, this will most likely not be in the same order as the tuple.

You can get help on different hints by executing help(ode.ode_hintname), where hintname is the name of the hint without _Integral.

See allhints or the ode docstring for a list of all supported hints that can be returned from classify_ode().

Notes

These are remarks on hint names.

_Integral

If a classification has _Integral at the end, it will return the expression with an unevaluated Integral class in it. Note that a hint may do this anyway if integrate() cannot do the integral, though just using an _Integral will do so much faster. Indeed, an _Integral hint will always be faster than its corresponding hint without _Integral because integrate() is an expensive routine. If dsolve() hangs, it is probably because integrate() is hanging on a tough or impossible integral. Try using an _Integral hint or all_Integral to get it return something.

Note that some hints do not have _Integral counterparts. This is because integrate() is not used in solving the ODE for those method. For example, \(n\)th order linear homogeneous ODEs with constant coefficients do not require integration to solve, so there is no nth_linear_homogeneous_constant_coeff_Integrate hint. You can easily evaluate any unevaluated Integrals in an expression by doing expr.doit().

Ordinals

Some hints contain an ordinal such as 1st_linear. This is to help differentiate them from other hints, as well as from other methods that may not be implemented yet. If a hint has nth in it, such as the nth_linear hints, this means that the method used to applies to ODEs of any order.

indep and dep

Some hints contain the words indep or dep. These reference the independent variable and the dependent function, respectively. For example, if an ODE is in terms of \(f(x)\), then indep will refer to \(x\) and dep will refer to \(f\).

subs

If a hints has the word subs in it, it means the the ODE is solved by substituting the expression given after the word subs for a single dummy variable. This is usually in terms of indep and dep as above. The substituted expression will be written only in characters allowed for names of Python objects, meaning operators will be spelled out. For example, indep/dep will be written as indep_div_dep.

coeff

The word coeff in a hint refers to the coefficients of something in the ODE, usually of the derivative terms. See the docstring for the individual methods for more info (help(ode)). This is contrast to coefficients, as in undetermined_coefficients, which refers to the common name of a method.

_best

Methods that have more than one fundamental way to solve will have a hint for each sub-method and a _best meta-classification. This will evaluate all hints and return the best, using the same considerations as the normal best meta-hint.

Examples

>>> from sympy import Function, classify_ode, Eq
>>> from sympy.abc import x
>>> f = Function('f')
>>> classify_ode(Eq(f(x).diff(x), 0), f(x))
('separable', '1st_linear', '1st_homogeneous_coeff_best',
'1st_homogeneous_coeff_subs_indep_div_dep',
'1st_homogeneous_coeff_subs_dep_div_indep',
'1st_power_series', 'lie_group',
'nth_linear_constant_coeff_homogeneous',
'separable_Integral', '1st_linear_Integral',
'1st_homogeneous_coeff_subs_indep_div_dep_Integral',
'1st_homogeneous_coeff_subs_dep_div_indep_Integral')
>>> classify_ode(f(x).diff(x, 2) + 3*f(x).diff(x) + 2*f(x) - 4)
('nth_linear_constant_coeff_undetermined_coefficients',
'nth_linear_constant_coeff_variation_of_parameters',
'nth_linear_constant_coeff_variation_of_parameters_Integral')

checkodesol()

sympy.solvers.ode.checkodesol(ode, sol, func=None, order='auto', solve_for_func=True)[source]

Substitutes sol into ode and checks that the result is 0.

This only works when func is one function, like \(f(x)\). sol can be a single solution or a list of solutions. Each solution may be an Equality that the solution satisfies, e.g. Eq(f(x), C1), Eq(f(x) + C1, 0); or simply an Expr, e.g. f(x) - C1. In most cases it will not be necessary to explicitly identify the function, but if the function cannot be inferred from the original equation it can be supplied through the func argument.

If a sequence of solutions is passed, the same sort of container will be used to return the result for each solution.

It tries the following methods, in order, until it finds zero equivalence:

  1. Substitute the solution for \(f\) in the original equation. This only works if ode is solved for \(f\). It will attempt to solve it first unless solve_for_func == False.

  2. Take \(n\) derivatives of the solution, where \(n\) is the order of ode, and check to see if that is equal to the solution. This only works on exact ODEs.

  3. Take the 1st, 2nd, …, \(n\)th derivatives of the solution, each time solving for the derivative of \(f\) of that order (this will always be possible because \(f\) is a linear operator). Then back substitute each derivative into ode in reverse order.

This function returns a tuple. The first item in the tuple is True if the substitution results in 0, and False otherwise. The second item in the tuple is what the substitution results in. It should always be 0 if the first item is True. Sometimes this function will return False even when an expression is identically equal to 0. This happens when simplify() does not reduce the expression to 0. If an expression returned by this function vanishes identically, then sol really is a solution to the ode.

If this function seems to hang, it is probably because of a hard simplification.

To use this function to test, test the first item of the tuple.

Examples

>>> from sympy import Eq, Function, checkodesol, symbols
>>> x, C1 = symbols('x,C1')
>>> f = Function('f')
>>> checkodesol(f(x).diff(x), Eq(f(x), C1))
(True, 0)
>>> assert checkodesol(f(x).diff(x), C1)[0]
>>> assert not checkodesol(f(x).diff(x), x)[0]
>>> checkodesol(f(x).diff(x, 2), x**2)
(False, 2)

homogeneous_order()

sympy.solvers.ode.homogeneous_order(eq, *symbols)[source]

Returns the order \(n\) if \(g\) is homogeneous and None if it is not homogeneous.

Determines if a function is homogeneous and if so of what order. A function \(f(x, y, \cdots)\) is homogeneous of order \(n\) if \(f(t x, t y, \cdots) = t^n f(x, y, \cdots)\).

If the function is of two variables, \(F(x, y)\), then \(f\) being homogeneous of any order is equivalent to being able to rewrite \(F(x, y)\) as \(G(x/y)\) or \(H(y/x)\). This fact is used to solve 1st order ordinary differential equations whose coefficients are homogeneous of the same order (see the docstrings of ode_1st_homogeneous_coeff_subs_dep_div_indep() and ode_1st_homogeneous_coeff_subs_indep_div_dep()).

Symbols can be functions, but every argument of the function must be a symbol, and the arguments of the function that appear in the expression must match those given in the list of symbols. If a declared function appears with different arguments than given in the list of symbols, None is returned.

Examples

>>> from sympy import Function, homogeneous_order, sqrt
>>> from sympy.abc import x, y
>>> f = Function('f')
>>> homogeneous_order(f(x), f(x)) is None
True
>>> homogeneous_order(f(x,y), f(y, x), x, y) is None
True
>>> homogeneous_order(f(x), f(x), x)
1
>>> homogeneous_order(x**2*f(x)/sqrt(x**2+f(x)**2), x, f(x))
2
>>> homogeneous_order(x**2+f(x), x, f(x)) is None
True

infinitesimals()

sympy.solvers.ode.infinitesimals(eq, func=None, order=None, hint='default', match=None)[source]

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\).

\[\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, diff, 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.

References

  • Solving differential equations by Symmetry Groups, John Starrett, pp. 1 - pp. 14

Examples

>>> from sympy import Function, diff
>>> from sympy.solvers.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}]

checkinfsol()

sympy.solvers.ode.checkinfsol(eq, infinitesimals, func=None, order=None)[source]

This function is used to check if the given infinitesimals are the actual infinitesimals of the given first order differential equation. This method is specific to the Lie Group Solver of ODEs.

As of now, it simply checks, by substituting the infinitesimals in the partial differential equation.

\[\frac{\partial \eta}{\partial x} + \left(\frac{\partial \eta}{\partial y} - \frac{\partial \xi}{\partial x}\right)*h - \frac{\partial \xi}{\partial y}*h^{2} - \xi\frac{\partial h}{\partial x} - \eta\frac{\partial h}{\partial y} = 0\]

where \(\eta\), and \(\xi\) are the infinitesimals and \(h(x,y) = \frac{dy}{dx}\)

The infinitesimals should be given in the form of a list of dicts [{xi(x, y): inf, eta(x, y): inf}], corresponding to the output of the function infinitesimals. It returns a list of values of the form [(True/False, sol)] where sol is the value obtained after substituting the infinitesimals in the PDE. If it is True, then sol would be 0.

Hint Functions

These functions are intended for internal use by dsolve() and others. Unlike User Functions, above, these are not intended for every-day use by ordinary SymPy users. Instead, functions such as dsolve() should be used. Nonetheless, these functions contain useful information in their docstrings on the various ODE solving methods. For this reason, they are documented here.

allhints

sympy.solvers.ode.allhints = ('separable', '1st_exact', '1st_linear', 'Bernoulli', 'Riccati_special_minus2', '1st_homogeneous_coeff_best', '1st_homogeneous_coeff_subs_indep_div_dep', '1st_homogeneous_coeff_subs_dep_div_indep', 'almost_linear', 'linear_coefficients', 'separable_reduced', '1st_power_series', 'lie_group', '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', 'separable_Integral', '1st_exact_Integral', '1st_linear_Integral', 'Bernoulli_Integral', '1st_homogeneous_coeff_subs_indep_div_dep_Integral', '1st_homogeneous_coeff_subs_dep_div_indep_Integral', 'almost_linear_Integral', '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')

This is a list of hints in the order that they should be preferred by classify_ode(). In general, hints earlier in the list should produce simpler solutions than those later in the list (for ODEs that fit both). For now, the order of this list is based on empirical observations by the developers of SymPy.

The hint used by dsolve() for a specific ODE can be overridden (see the docstring).

In general, _Integral hints are grouped at the end of the list, unless there is a method that returns an unevaluable integral most of the time (which go near the end of the list anyway). default, all, best, and all_Integral meta-hints should not be included in this list, but _best and _Integral hints should be included.

odesimp

sympy.solvers.ode.odesimp(eq, func, order, constants, hint)[source]

Simplifies ODEs, including trying to solve for func and running constantsimp().

It may use knowledge of the type of solution that the hint returns to apply additional simplifications.

It also attempts to integrate any Integrals in the expression, if the hint is not an _Integral hint.

This function should have no effect on expressions returned by dsolve(), as dsolve() already calls odesimp(), but the individual hint functions do not call odesimp() (because the dsolve() wrapper does). Therefore, this function is designed for mainly internal use.

Examples

>>> from sympy import sin, symbols, dsolve, pprint, Function
>>> from sympy.solvers.ode import odesimp
>>> x , u2, C1= symbols('x,u2,C1')
>>> f = Function('f')
>>> eq = dsolve(x*f(x).diff(x) - f(x) - x*sin(f(x)/x), f(x),
... hint='1st_homogeneous_coeff_subs_indep_div_dep_Integral',
... simplify=False)
>>> pprint(eq, wrap_line=False)
                        x
                       ----
                       f(x)
                         /
                        |
                        |   /        1   \
                        |  -|u2 + -------|
                        |   |        /1 \|
                        |   |     sin|--||
                        |   \        \u2//
log(f(x)) = log(C1) +   |  ---------------- d(u2)
                        |          2
                        |        u2
                        |
                       /
>>> pprint(odesimp(eq, f(x), 1, {C1},
... hint='1st_homogeneous_coeff_subs_indep_div_dep'
... )) 
    x
--------- = C1
   /f(x)\
tan|----|
   \2*x /

constant_renumber

sympy.solvers.ode.constant_renumber(expr, symbolname, startnumber, endnumber)[source]

Renumber arbitrary constants in expr to have numbers 1 through \(N\) where \(N\) is endnumber - startnumber + 1 at most. In the process, this reorders expression terms in a standard way.

This is a simple function that goes through and renumbers any Symbol with a name in the form symbolname + num where num is in the range from startnumber to endnumber.

Symbols are renumbered based on .sort_key(), so they should be numbered roughly in the order that they appear in the final, printed expression. Note that this ordering is based in part on hashes, so it can produce different results on different machines.

The structure of this function is very similar to that of constantsimp().

Examples

>>> from sympy import symbols, Eq, pprint
>>> from sympy.solvers.ode import constant_renumber
>>> x, C0, C1, C2, C3, C4 = symbols('x,C:5')

Only constants in the given range (inclusive) are renumbered; the renumbering always starts from 1:

>>> constant_renumber(C1 + C3 + C4, 'C', 1, 3)
C1 + C2 + C4
>>> constant_renumber(C0 + C1 + C3 + C4, 'C', 2, 4)
C0 + 2*C1 + C2
>>> constant_renumber(C0 + 2*C1 + C2, 'C', 0, 1)
C1 + 3*C2
>>> pprint(C2 + C1*x + C3*x**2)
                2
C1*x + C2 + C3*x
>>> pprint(constant_renumber(C2 + C1*x + C3*x**2, 'C', 1, 3))
                2
C1 + C2*x + C3*x

constantsimp

sympy.solvers.ode.constantsimp(expr, constants)[source]

Simplifies an expression with arbitrary constants in it.

This function is written specifically to work with dsolve(), and is not intended for general use.

Simplification is done by “absorbing” the arbitrary constants into other arbitrary constants, numbers, and symbols that they are not independent of.

The symbols must all have the same name with numbers after it, for example, C1, C2, C3. The symbolname here would be ‘C’, the startnumber would be 1, and the endnumber would be 3. If the arbitrary constants are independent of the variable x, then the independent symbol would be x. There is no need to specify the dependent function, such as f(x), because it already has the independent symbol, x, in it.

Because terms are “absorbed” into arbitrary constants and because constants are renumbered after simplifying, the arbitrary constants in expr are not necessarily equal to the ones of the same name in the returned result.

If two or more arbitrary constants are added, multiplied, or raised to the power of each other, they are first absorbed together into a single arbitrary constant. Then the new constant is combined into other terms if necessary.

Absorption of constants is done with limited assistance:

  1. terms of Adds are collected to try join constants so \(e^x (C_1 \cos(x) + C_2 \cos(x))\) will simplify to \(e^x C_1 \cos(x)\);

  2. powers with exponents that are Adds are expanded so \(e^{C_1 + x}\) will be simplified to \(C_1 e^x\).

Use constant_renumber() to renumber constants after simplification or else arbitrary numbers on constants may appear, e.g. \(C_1 + C_3 x\).

In rare cases, a single constant can be “simplified” into two constants. Every differential equation solution should have as many arbitrary constants as the order of the differential equation. The result here will be technically correct, but it may, for example, have \(C_1\) and \(C_2\) in an expression, when \(C_1\) is actually equal to \(C_2\). Use your discretion in such situations, and also take advantage of the ability to use hints in dsolve().

Examples

>>> from sympy import symbols
>>> from sympy.solvers.ode import constantsimp
>>> C1, C2, C3, x, y = symbols('C1, C2, C3, x, y')
>>> constantsimp(2*C1*x, {C1, C2, C3})
C1*x
>>> constantsimp(C1 + 2 + x, {C1, C2, C3})
C1 + x
>>> constantsimp(C1*C2 + 2 + C2 + C3*x, {C1, C2, C3})
C1 + C3*x

sol_simplicity

sympy.solvers.ode.ode_sol_simplicity(sol, func, trysolving=True)[source]

Returns an extended integer representing how simple a solution to an ODE is.

The following things are considered, in order from most simple to least:

  • sol is solved for func.

  • sol is not solved for func, but can be if passed to solve (e.g., a solution returned by dsolve(ode, func, simplify=False).

  • If sol is not solved for func, then base the result on the length of sol, as computed by len(str(sol)).

  • If sol has any unevaluated Integrals, this will automatically be considered less simple than any of the above.

This function returns an integer such that if solution A is simpler than solution B by above metric, then ode_sol_simplicity(sola, func) < ode_sol_simplicity(solb, func).

Currently, the following are the numbers returned, but if the heuristic is ever improved, this may change. Only the ordering is guaranteed.

Simplicity

Return

sol solved for func

-2

sol not solved for func but can be

-1

sol is not solved nor solvable for func

len(str(sol))

sol contains an Integral

oo

oo here means the SymPy infinity, which should compare greater than any integer.

If you already know solve() cannot solve sol, you can use trysolving=False to skip that step, which is the only potentially slow step. For example, dsolve() with the simplify=False flag should do this.

If sol is a list of solutions, if the worst solution in the list returns oo it returns that, otherwise it returns len(str(sol)), that is, the length of the string representation of the whole list.

Examples

This function is designed to be passed to min as the key argument, such as min(listofsolutions, key=lambda i: ode_sol_simplicity(i, f(x))).

>>> from sympy import symbols, Function, Eq, tan, cos, sqrt, Integral
>>> from sympy.solvers.ode import ode_sol_simplicity
>>> x, C1, C2 = symbols('x, C1, C2')
>>> f = Function('f')
>>> ode_sol_simplicity(Eq(f(x), C1*x**2), f(x))
-2
>>> ode_sol_simplicity(Eq(x**2 + f(x), C1), f(x))
-1
>>> ode_sol_simplicity(Eq(f(x), C1*Integral(2*x, x)), f(x))
oo
>>> eq1 = Eq(f(x)/tan(f(x)/(2*x)), C1)
>>> eq2 = Eq(f(x)/tan(f(x)/(2*x) + f(x)), C2)
>>> [ode_sol_simplicity(eq, f(x)) for eq in [eq1, eq2]]
[28, 35]
>>> min([eq1, eq2], key=lambda i: ode_sol_simplicity(i, f(x)))
Eq(f(x)/tan(f(x)/(2*x)), C1)

1st_exact

sympy.solvers.ode.ode_1st_exact(eq, func, order, match)[source]

Solves 1st order exact ordinary differential equations.

A 1st order differential equation is called exact if it is the total differential of a function. That is, the differential equation

\[P(x, y) \,\partial{}x + Q(x, y) \,\partial{}y = 0\]

is exact if there is some function \(F(x, y)\) such that \(P(x, y) = \partial{}F/\partial{}x\) and \(Q(x, y) = \partial{}F/\partial{}y\). It can be shown that a necessary and sufficient condition for a first order ODE to be exact is that \(\partial{}P/\partial{}y = \partial{}Q/\partial{}x\). Then, the solution will be as given below:

>>> from sympy import Function, Eq, Integral, symbols, pprint
>>> x, y, t, x0, y0, C1= symbols('x,y,t,x0,y0,C1')
>>> P, Q, F= map(Function, ['P', 'Q', 'F'])
>>> pprint(Eq(Eq(F(x, y), Integral(P(t, y), (t, x0, x)) +
... Integral(Q(x0, t), (t, y0, y))), C1))
            x                y
            /                /
           |                |
F(x, y) =  |  P(t, y) dt +  |  Q(x0, t) dt = C1
           |                |
          /                /
          x0               y0

Where the first partials of \(P\) and \(Q\) exist and are continuous in a simply connected region.

A note: SymPy currently has no way to represent inert substitution on an expression, so the hint 1st_exact_Integral will return an integral with \(dy\). This is supposed to represent the function that you are solving for.

References

# indirect doctest

Examples

>>> from sympy import Function, dsolve, cos, sin
>>> from sympy.abc import x
>>> f = Function('f')
>>> dsolve(cos(f(x)) - (x*sin(f(x)) - f(x)**2)*f(x).diff(x),
... f(x), hint='1st_exact')
Eq(x*cos(f(x)) + f(x)**3/3, C1)

1st_homogeneous_coeff_best

sympy.solvers.ode.ode_1st_homogeneous_coeff_best(eq, func, order, match)[source]

Returns the best solution to an ODE from the two hints 1st_homogeneous_coeff_subs_dep_div_indep and 1st_homogeneous_coeff_subs_indep_div_dep.

This is as determined by ode_sol_simplicity().

See the ode_1st_homogeneous_coeff_subs_indep_div_dep() and ode_1st_homogeneous_coeff_subs_dep_div_indep() docstrings for more information on these hints. Note that there is no ode_1st_homogeneous_coeff_best_Integral hint.

References

# indirect doctest

Examples

>>> from sympy import Function, dsolve, pprint
>>> from sympy.abc import x
>>> f = Function('f')
>>> pprint(dsolve(2*x*f(x) + (x**2 + f(x)**2)*f(x).diff(x), f(x),
... hint='1st_homogeneous_coeff_best', simplify=False))
                         /    2    \
                         | 3*x     |
                      log|----- + 1|
                         | 2       |
                         \f (x)    /
log(f(x)) = log(C1) - --------------
                            3

1st_homogeneous_coeff_subs_dep_div_indep

sympy.solvers.ode.ode_1st_homogeneous_coeff_subs_dep_div_indep(eq, func, order, match)[source]

Solves a 1st order differential equation with homogeneous coefficients using the substitution \(u_1 = \frac{\text{<dependent variable>}}{\text{<independent variable>}}\).

This is a differential equation

\[P(x, y) + Q(x, y) dy/dx = 0\]

such that \(P\) and \(Q\) are homogeneous and of the same order. A function \(F(x, y)\) is homogeneous of order \(n\) if \(F(x t, y t) = t^n F(x, y)\). Equivalently, \(F(x, y)\) can be rewritten as \(G(y/x)\) or \(H(x/y)\). See also the docstring of homogeneous_order().

If the coefficients \(P\) and \(Q\) in the differential equation above are homogeneous functions of the same order, then it can be shown that the substitution \(y = u_1 x\) (i.e. \(u_1 = y/x\)) will turn the differential equation into an equation separable in the variables \(x\) and \(u\). If \(h(u_1)\) is the function that results from making the substitution \(u_1 = f(x)/x\) on \(P(x, f(x))\) and \(g(u_2)\) is the function that results from the substitution on \(Q(x, f(x))\) in the differential equation \(P(x, f(x)) + Q(x, f(x)) f'(x) = 0\), then the general solution is:

>>> from sympy import Function, dsolve, pprint
>>> from sympy.abc import x
>>> f, g, h = map(Function, ['f', 'g', 'h'])
>>> genform = g(f(x)/x) + h(f(x)/x)*f(x).diff(x)
>>> pprint(genform)
 /f(x)\    /f(x)\ d
g|----| + h|----|*--(f(x))
 \ x  /    \ x  / dx
>>> pprint(dsolve(genform, f(x),
... hint='1st_homogeneous_coeff_subs_dep_div_indep_Integral'))
               f(x)
               ----
                x
                 /
                |
                |       -h(u1)
log(x) = C1 +   |  ---------------- d(u1)
                |  u1*h(u1) + g(u1)
                |
               /

Where \(u_1 h(u_1) + g(u_1) \ne 0\) and \(x \ne 0\).

See also the docstrings of ode_1st_homogeneous_coeff_best() and ode_1st_homogeneous_coeff_subs_indep_div_dep().

References

# indirect doctest

Examples

>>> from sympy import Function, dsolve
>>> from sympy.abc import x
>>> f = Function('f')
>>> pprint(dsolve(2*x*f(x) + (x**2 + f(x)**2)*f(x).diff(x), f(x),
... hint='1st_homogeneous_coeff_subs_dep_div_indep', simplify=False))
                      /          3   \
                      |3*f(x)   f (x)|
                   log|------ + -----|
                      |  x         3 |
                      \           x  /
log(x) = log(C1) - -------------------
                            3

1st_homogeneous_coeff_subs_indep_div_dep

sympy.solvers.ode.ode_1st_homogeneous_coeff_subs_indep_div_dep(eq, func, order, match)[source]

Solves a 1st order differential equation with homogeneous coefficients using the substitution \(u_2 = \frac{\text{<independent variable>}}{\text{<dependent variable>}}\).

This is a differential equation

\[P(x, y) + Q(x, y) dy/dx = 0\]

such that \(P\) and \(Q\) are homogeneous and of the same order. A function \(F(x, y)\) is homogeneous of order \(n\) if \(F(x t, y t) = t^n F(x, y)\). Equivalently, \(F(x, y)\) can be rewritten as \(G(y/x)\) or \(H(x/y)\). See also the docstring of homogeneous_order().

If the coefficients \(P\) and \(Q\) in the differential equation above are homogeneous functions of the same order, then it can be shown that the substitution \(x = u_2 y\) (i.e. \(u_2 = x/y\)) will turn the differential equation into an equation separable in the variables \(y\) and \(u_2\). If \(h(u_2)\) is the function that results from making the substitution \(u_2 = x/f(x)\) on \(P(x, f(x))\) and \(g(u_2)\) is the function that results from the substitution on \(Q(x, f(x))\) in the differential equation \(P(x, f(x)) + Q(x, f(x)) f'(x) = 0\), then the general solution is:

>>> from sympy import Function, dsolve, pprint
>>> from sympy.abc import x
>>> f, g, h = map(Function, ['f', 'g', 'h'])
>>> genform = g(x/f(x)) + h(x/f(x))*f(x).diff(x)
>>> pprint(genform)
 / x  \    / x  \ d
g|----| + h|----|*--(f(x))
 \f(x)/    \f(x)/ dx
>>> pprint(dsolve(genform, f(x),
... hint='1st_homogeneous_coeff_subs_indep_div_dep_Integral'))
             x
            ----
            f(x)
              /
             |
             |       -g(u2)
             |  ---------------- d(u2)
             |  u2*g(u2) + h(u2)
             |
            /

f(x) = C1*e

Where \(u_2 g(u_2) + h(u_2) \ne 0\) and \(f(x) \ne 0\).

See also the docstrings of ode_1st_homogeneous_coeff_best() and ode_1st_homogeneous_coeff_subs_dep_div_indep().

References

# indirect doctest

Examples

>>> from sympy import Function, pprint, dsolve
>>> from sympy.abc import x
>>> f = Function('f')
>>> pprint(dsolve(2*x*f(x) + (x**2 + f(x)**2)*f(x).diff(x), f(x),
... hint='1st_homogeneous_coeff_subs_indep_div_dep',
... simplify=False))
                         /    2    \
                         | 3*x     |
                      log|----- + 1|
                         | 2       |
                         \f (x)    /
log(f(x)) = log(C1) - --------------
                            3

1st_linear

sympy.solvers.ode.ode_1st_linear(eq, func, order, match)[source]

Solves 1st order linear differential equations.

These are differential equations of the form

\[dy/dx + P(x) y = Q(x)\text{.}\]

These kinds of differential equations can be solved in a general way. The integrating factor \(e^{\int P(x) \,dx}\) will turn the equation into a separable equation. The general solution is:

>>> from sympy import Function, dsolve, Eq, pprint, diff, sin
>>> from sympy.abc import x
>>> f, P, Q = map(Function, ['f', 'P', 'Q'])
>>> genform = Eq(f(x).diff(x) + P(x)*f(x), Q(x))
>>> pprint(genform)
            d
P(x)*f(x) + --(f(x)) = Q(x)
            dx
>>> pprint(dsolve(genform, f(x), hint='1st_linear_Integral'))
       /       /                   \
       |      |                    |
       |      |         /          |     /
       |      |        |           |    |
       |      |        | P(x) dx   |  - | P(x) dx
       |      |        |           |    |
       |      |       /            |   /
f(x) = |C1 +  | Q(x)*e           dx|*e
       |      |                    |
       \     /                     /

References

# indirect doctest

Examples

>>> f = Function('f')
>>> pprint(dsolve(Eq(x*diff(f(x), x) - f(x), x**2*sin(x)),
... f(x), '1st_linear'))
f(x) = x*(C1 - cos(x))

Bernoulli

sympy.solvers.ode.ode_Bernoulli(eq, func, order, match)[source]

Solves Bernoulli differential equations.

These are equations of the form

\[dy/dx + P(x) y = Q(x) y^n\text{, }n \ne 1`\text{.}\]

The substitution \(w = 1/y^{1-n}\) will transform an equation of this form into one that is linear (see the docstring of ode_1st_linear()). The general solution is:

>>> from sympy import Function, dsolve, Eq, pprint
>>> from sympy.abc import x, n
>>> f, P, Q = map(Function, ['f', 'P', 'Q'])
>>> genform = Eq(f(x).diff(x) + P(x)*f(x), Q(x)*f(x)**n)
>>> pprint(genform)
            d                n
P(x)*f(x) + --(f(x)) = Q(x)*f (x)
            dx
>>> pprint(dsolve(genform, f(x), hint='Bernoulli_Integral')) 
                                                                               1
                                                                              ----
                                                                             1 - n
       //                /                            \                     \
       ||               |                             |                     |
       ||               |                  /          |             /       |
       ||               |                 |           |            |        |
       ||               |        (1 - n)* | P(x) dx   |  (-1 + n)* | P(x) dx|
       ||               |                 |           |            |        |
       ||               |                /            |           /         |
f(x) = ||C1 + (-1 + n)* | -Q(x)*e                   dx|*e                   |
       ||               |                             |                     |
       \\               /                            /                     /

Note that the equation is separable when \(n = 1\) (see the docstring of ode_separable()).

>>> pprint(dsolve(Eq(f(x).diff(x) + P(x)*f(x), Q(x)*f(x)), f(x),
... hint='separable_Integral'))
 f(x)
   /
  |                /
  |  1            |
  |  - dy = C1 +  | (-P(x) + Q(x)) dx
  |  y            |
  |              /
 /

References

# indirect doctest

Examples

>>> from sympy import Function, dsolve, Eq, pprint, log
>>> from sympy.abc import x
>>> f = Function('f')
>>> pprint(dsolve(Eq(x*f(x).diff(x) + f(x), log(x)*f(x)**2),
... f(x), hint='Bernoulli'))
                1
f(x) = -------------------
         /     log(x)   1\
       x*|C1 + ------ + -|
         \       x      x/

Liouville

sympy.solvers.ode.ode_Liouville(eq, func, order, match)[source]

Solves 2nd order Liouville differential equations.

The general form of a Liouville ODE is

\[\frac{d^2 y}{dx^2} + g(y) \left(\! \frac{dy}{dx}\!\right)^2 + h(x) \frac{dy}{dx}\text{.}\]

The general solution is:

>>> from sympy import Function, dsolve, Eq, pprint, diff
>>> from sympy.abc import x
>>> f, g, h = map(Function, ['f', 'g', 'h'])
>>> genform = Eq(diff(f(x),x,x) + g(f(x))*diff(f(x),x)**2 +
... h(x)*diff(f(x),x), 0)
>>> pprint(genform)
                  2                    2
        /d       \         d          d
g(f(x))*|--(f(x))|  + h(x)*--(f(x)) + ---(f(x)) = 0
        \dx      /         dx           2
                                      dx
>>> pprint(dsolve(genform, f(x), hint='Liouville_Integral'))
                                  f(x)
          /                     /
         |                     |
         |     /               |     /
         |    |                |    |
         |  - | h(x) dx        |    | g(y) dy
         |    |                |    |
         |   /                 |   /
C1 + C2* | e            dx +   |  e           dy = 0
         |                     |
        /                     /

References

# indirect doctest

Examples

>>> from sympy import Function, dsolve, Eq, pprint
>>> from sympy.abc import x
>>> f = Function('f')
>>> pprint(dsolve(diff(f(x), x, x) + diff(f(x), x)**2/f(x) +
... diff(f(x), x)/x, f(x), hint='Liouville'))
           ________________           ________________
[f(x) = -\/ C1 + C2*log(x) , f(x) = \/ C1 + C2*log(x) ]

Riccati_special_minus2

sympy.solvers.ode.ode_Riccati_special_minus2(eq, func, order, match)[source]

The general Riccati equation has the form

\[dy/dx = f(x) y^2 + g(x) y + h(x)\text{.}\]

While it does not have a general solution [1], the “special” form, \(dy/dx = a y^2 - b x^c\), does have solutions in many cases [2]. This routine returns a solution for \(a(dy/dx) = b y^2 + c y/x + d/x^2\) that is obtained by using a suitable change of variables to reduce it to the special form and is valid when neither \(a\) nor \(b\) are zero and either \(c\) or \(d\) is zero.

>>> from sympy.abc import x, y, a, b, c, d
>>> from sympy.solvers.ode import dsolve, checkodesol
>>> from sympy import pprint, Function
>>> f = Function('f')
>>> y = f(x)
>>> genform = a*y.diff(x) - (b*y**2 + c*y/x + d/x**2)
>>> sol = dsolve(genform, y)
>>> pprint(sol, wrap_line=False)
        /                                 /        __________________       \\
        |           __________________    |       /                2        ||
        |          /                2     |     \/  4*b*d - (a + c)  *log(x)||
       -|a + c - \/  4*b*d - (a + c)  *tan|C1 + ----------------------------||
        \                                 \                 2*a             //
f(x) = ------------------------------------------------------------------------
                                        2*b*x
>>> checkodesol(genform, sol, order=1)[0]
True

References

  1. http://www.maplesoft.com/support/help/Maple/view.aspx?path=odeadvisor/Riccati

  2. http://eqworld.ipmnet.ru/en/solutions/ode/ode0106.pdf - http://eqworld.ipmnet.ru/en/solutions/ode/ode0123.pdf

nth_linear_constant_coeff_homogeneous

sympy.solvers.ode.ode_nth_linear_constant_coeff_homogeneous(eq, func, order, match, returns='sol')[source]

Solves an \(n\)th order linear homogeneous differential equation with constant coefficients.

This is an equation of the form

\[a_n f^{(n)}(x) + a_{n-1} f^{(n-1)}(x) + \cdots + a_1 f'(x) + a_0 f(x) = 0\text{.}\]

These equations can be solved in a general manner, by taking the roots of the characteristic equation \(a_n m^n + a_{n-1} m^{n-1} + \cdots + a_1 m + a_0 = 0\). The solution will then be the sum of \(C_n x^i e^{r x}\) terms, for each where \(C_n\) is an arbitrary constant, \(r\) is a root of the characteristic equation and \(i\) is one of each from 0 to the multiplicity of the root - 1 (for example, a root 3 of multiplicity 2 would create the terms \(C_1 e^{3 x} + C_2 x e^{3 x}\)). The exponential is usually expanded for complex roots using Euler’s equation \(e^{I x} = \cos(x) + I \sin(x)\). Complex roots always come in conjugate pairs in polynomials with real coefficients, so the two roots will be represented (after simplifying the constants) as \(e^{a x} \left(C_1 \cos(b x) + C_2 \sin(b x)\right)\).

If SymPy cannot find exact roots to the characteristic equation, a CRootOf instance will be return instead.

>>> from sympy import Function, dsolve, Eq
>>> from sympy.abc import x
>>> f = Function('f')
>>> dsolve(f(x).diff(x, 5) + 10*f(x).diff(x) - 2*f(x), f(x),
... hint='nth_linear_constant_coeff_homogeneous')
... 
Eq(f(x), C1*exp(x*CRootOf(_x**5 + 10*_x - 2, 0)) +
C2*exp(x*CRootOf(_x**5 + 10*_x - 2, 1)) +
C3*exp(x*CRootOf(_x**5 + 10*_x - 2, 2)) +
C4*exp(x*CRootOf(_x**5 + 10*_x - 2, 3)) +
C5*exp(x*CRootOf(_x**5 + 10*_x - 2, 4)))

Note that because this method does not involve integration, there is no nth_linear_constant_coeff_homogeneous_Integral hint.

The following is for internal use:

  • returns = 'sol' returns the solution to the ODE.

  • returns = 'list' returns a list of linearly independent solutions, for use with non homogeneous solution methods like variation of parameters and undetermined coefficients. Note that, though the solutions should be linearly independent, this function does not explicitly check that. You can do assert simplify(wronskian(sollist)) != 0 to check for linear independence. Also, assert len(sollist) == order will need to pass.

  • returns = 'both', return a dictionary {'sol': <solution to ODE>, 'list': <list of linearly independent solutions>}.

References

# indirect doctest

Examples

>>> from sympy import Function, dsolve, pprint
>>> from sympy.abc import x
>>> f = Function('f')
>>> pprint(dsolve(f(x).diff(x, 4) + 2*f(x).diff(x, 3) -
... 2*f(x).diff(x, 2) - 6*f(x).diff(x) + 5*f(x), f(x),
... hint='nth_linear_constant_coeff_homogeneous'))
                    x                            -2*x
f(x) = (C1 + C2*x)*e  + (C3*sin(x) + C4*cos(x))*e

nth_linear_constant_coeff_undetermined_coefficients

sympy.solvers.ode.ode_nth_linear_constant_coeff_undetermined_coefficients(eq, func, order, match)[source]

Solves an \(n\)th order linear differential equation with constant coefficients using the method of undetermined coefficients.

This method works on differential equations of the form

\[a_n f^{(n)}(x) + a_{n-1} f^{(n-1)}(x) + \cdots + a_1 f'(x) + a_0 f(x) = P(x)\text{,}\]

where \(P(x)\) is a function that has a finite number of linearly 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.

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.

References

# indirect doctest

Examples

>>> from sympy import Function, dsolve, pprint, exp, cos
>>> from sympy.abc import x
>>> f = Function('f')
>>> pprint(dsolve(f(x).diff(x, 2) + 2*f(x).diff(x) + f(x) -
... 4*exp(-x)*x**2 + cos(2*x), f(x),
... hint='nth_linear_constant_coeff_undetermined_coefficients'))
       /             4\
       |            x |  -x   4*sin(2*x)   3*cos(2*x)
f(x) = |C1 + C2*x + --|*e   - ---------- + ----------
       \            3 /           25           25

nth_linear_constant_coeff_variation_of_parameters

sympy.solvers.ode.ode_nth_linear_constant_coeff_variation_of_parameters(eq, func, order, match)[source]

Solves an \(n\)th order linear differential equation with constant coefficients using the method of variation of parameters.

This method works on any differential equations of the form

\[f^{(n)}(x) + a_{n-1} f^{(n-1)}(x) + \cdots + a_1 f'(x) + a_0 f(x) = P(x)\text{.}\]

This method works by assuming that the particular solution takes the form

\[\sum_{x=1}^{n} c_i(x) y_i(x)\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

\[\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, P(x)]\).

This method is general enough to solve any \(n\)th order inhomogeneous linear differential equation with constant coefficients, 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 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.

References

# indirect doctest

Examples

>>> from sympy import Function, dsolve, pprint, exp, log
>>> from sympy.abc import x
>>> f = Function('f')
>>> pprint(dsolve(f(x).diff(x, 3) - 3*f(x).diff(x, 2) +
... 3*f(x).diff(x) - f(x) - exp(x)*log(x), f(x),
... hint='nth_linear_constant_coeff_variation_of_parameters'))
       /                     3                \
       |                2   x *(6*log(x) - 11)|  x
f(x) = |C1 + C2*x + C3*x  + ------------------|*e
       \                            36        /

separable

sympy.solvers.ode.ode_separable(eq, func, order, match)[source]

Solves separable 1st order differential equations.

This is any differential equation that can be written as \(P(y) \tfrac{dy}{dx} = Q(x)\). The solution can then just be found by rearranging terms and integrating: \(\int P(y) \,dy = \int Q(x) \,dx\). This hint uses sympy.simplify.simplify.separatevars() as its back end, so if a separable equation is not caught by this solver, it is most likely the fault of that function. separatevars() is smart enough to do most expansion and factoring necessary to convert a separable equation \(F(x, y)\) into the proper form \(P(x)\cdot{}Q(y)\). The general solution is:

>>> from sympy import Function, dsolve, Eq, pprint
>>> from sympy.abc import x
>>> a, b, c, d, f = map(Function, ['a', 'b', 'c', 'd', 'f'])
>>> genform = Eq(a(x)*b(f(x))*f(x).diff(x), c(x)*d(f(x)))
>>> pprint(genform)
             d
a(x)*b(f(x))*--(f(x)) = c(x)*d(f(x))
             dx
>>> pprint(dsolve(genform, f(x), hint='separable_Integral'))
     f(x)
   /                  /
  |                  |
  |  b(y)            | c(x)
  |  ---- dy = C1 +  | ---- dx
  |  d(y)            | a(x)
  |                  |
 /                  /

References

  • M. Tenenbaum & H. Pollard, “Ordinary Differential Equations”, Dover 1963, pp. 52

# indirect doctest

Examples

>>> from sympy import Function, dsolve, Eq
>>> from sympy.abc import x
>>> f = Function('f')
>>> pprint(dsolve(Eq(f(x)*f(x).diff(x) + x, 3*x*f(x)**2), f(x),
... hint='separable', simplify=False))
   /   2       \         2
log\3*f (x) - 1/        x
---------------- = C1 + --
       6                2

almost_linear

sympy.solvers.ode.ode_almost_linear(eq, func, order, match)[source]

Solves an almost-linear differential equation.

The general form of an almost linear differential equation is

\[f(x) g(y) y + k(x) l(y) + m(x) = 0 \text{where} l'(y) = g(y)\text{.}\]

This can be solved by substituting \(l(y) = u(y)\). Making the given substitution reduces it to a linear differential equation of the form \(u' + P(x) u + Q(x) = 0\).

The general solution is

>>> from sympy import Function, dsolve, Eq, pprint
>>> from sympy.abc import x, y, n
>>> f, g, k, l = map(Function, ['f', 'g', 'k', 'l'])
>>> genform = Eq(f(x)*(l(y).diff(y)) + k(x)*l(y) + g(x))
>>> pprint(genform)
     d
f(x)*--(l(y)) + g(x) + k(x)*l(y) = 0
     dy
>>> pprint(dsolve(genform, hint = 'almost_linear'))
       /     //       y*k(x)                \\
       |     ||       ------                ||
       |     ||        f(x)                 ||  -y*k(x)
       |     ||-g(x)*e                      ||  --------
       |     ||--------------  for k(x) != 0||    f(x)
l(y) = |C1 + |<     k(x)                    ||*e
       |     ||                             ||
       |     ||   -y*g(x)                   ||
       |     ||   --------       otherwise  ||
       |     ||     f(x)                    ||
       \     \\                             //

References

  • Joel Moses, “Symbolic Integration - The Stormy Decade”, Communications of the ACM, Volume 14, Number 8, August 1971, pp. 558

Examples

>>> from sympy import Function, Derivative, pprint
>>> from sympy.solvers.ode import dsolve, classify_ode
>>> from sympy.abc import x
>>> f = Function('f')
>>> d = f(x).diff(x)
>>> eq = x*d + x*f(x) + 1
>>> dsolve(eq, f(x), hint='almost_linear')
Eq(f(x), (C1 - Ei(x))*exp(-x))
>>> pprint(dsolve(eq, f(x), hint='almost_linear'))
                     -x
f(x) = (C1 - Ei(x))*e

linear_coefficients

sympy.solvers.ode.ode_linear_coefficients(eq, func, order, match)[source]

Solves a differential equation with linear coefficients.

The general form of a differential equation with linear coefficients is

\[y' + F\left(\!\frac{a_1 x + b_1 y + c_1}{a_2 x + b_2 y + c_2}\!\right) = 0\text{,}\]

where \(a_1\), \(b_1\), \(c_1\), \(a_2\), \(b_2\), \(c_2\) are constants and \(a_1 b_2 - a_2 b_1 \ne 0\).

This can be solved by substituting:

\[ \begin{align}\begin{aligned}x = x' + \frac{b_2 c_1 - b_1 c_2}{a_2 b_1 - a_1 b_2}\\y = y' + \frac{a_1 c_2 - a_2 c_1}{a_2 b_1 - a_1 b_2}\text{.}\end{aligned}\end{align} \]

This substitution reduces the equation to a homogeneous differential equation.

References

  • Joel Moses, “Symbolic Integration - The Stormy Decade”, Communications of the ACM, Volume 14, Number 8, August 1971, pp. 558

Examples

>>> from sympy import Function, Derivative, pprint
>>> from sympy.solvers.ode import dsolve, classify_ode
>>> from sympy.abc import x
>>> f = Function('f')
>>> df = f(x).diff(x)
>>> eq = (x + f(x) + 1)*df + (f(x) - 6*x + 1)
>>> dsolve(eq, hint='linear_coefficients')
[Eq(f(x), -x - sqrt(C1 + 7*x**2) - 1), Eq(f(x), -x + sqrt(C1 + 7*x**2) - 1)]
>>> pprint(dsolve(eq, hint='linear_coefficients'))
                  ___________                     ___________
               /         2                     /         2
[f(x) = -x - \/  C1 + 7*x   - 1, f(x) = -x + \/  C1 + 7*x   - 1]

separable_reduced

sympy.solvers.ode.ode_separable_reduced(eq, func, order, match)[source]

Solves a differential equation that can be reduced to the separable form.

The general form of this equation is

\[y' + (y/x) H(x^n y) = 0\text{}.\]

This can be solved by substituting \(u(y) = x^n y\). The equation then reduces to the separable form \(\frac{u'}{u (\mathrm{power} - H(u))} - \frac{1}{x} = 0\).

The general solution is:

>>> from sympy import Function, dsolve, Eq, pprint
>>> from sympy.abc import x, n
>>> f, g = map(Function, ['f', 'g'])
>>> genform = f(x).diff(x) + (f(x)/x)*g(x**n*f(x))
>>> pprint(genform)
                 / n     \
d          f(x)*g\x *f(x)/
--(f(x)) + ---------------
dx                x
>>> pprint(dsolve(genform, hint='separable_reduced'))
 n
x *f(x)
  /
 |
 |         1
 |    ------------ dy = C1 + log(x)
 |    y*(n - g(y))
 |
 /

References

  • Joel Moses, “Symbolic Integration - The Stormy Decade”, Communications of the ACM, Volume 14, Number 8, August 1971, pp. 558

Examples

>>> from sympy import Function, Derivative, pprint
>>> from sympy.solvers.ode import dsolve, classify_ode
>>> from sympy.abc import x
>>> f = Function('f')
>>> d = f(x).diff(x)
>>> eq = (x - x**2*f(x))*d - f(x)
>>> dsolve(eq, hint='separable_reduced')
[Eq(f(x), (-sqrt(C1*x**2 + 1) + 1)/x), Eq(f(x), (sqrt(C1*x**2 + 1) + 1)/x)]
>>> pprint(dsolve(eq, hint='separable_reduced'))
             ___________                ___________
            /     2                    /     2
        - \/  C1*x  + 1  + 1         \/  C1*x  + 1  + 1
[f(x) = --------------------, f(x) = ------------------]
                 x                           x

lie_group

sympy.solvers.ode.ode_lie_group(eq, func, order, match)[source]

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 given system into another coordinate system where it becomes invariant under the one-parameter Lie group of translations. The converted ODE is quadrature and can be solved easily. It makes use of the 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.

\[\xi\frac{\partial r}{\partial x} + \eta\frac{\partial r}{\partial y} = 0\]
\[\xi\frac{\partial s}{\partial x} + \eta\frac{\partial s}{\partial y} = 1\]

The differential equation becomes separable in the new coordinate system

\[\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.

References

  • Solving differential equations by Symmetry Groups, John Starrett, pp. 1 - pp. 14

Examples

>>> from sympy import Function, dsolve, Eq, 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 /

1st_power_series

sympy.solvers.ode.ode_1st_power_series(eq, func, order, match)[source]

The power series solution is a method which gives the Taylor series expansion to the solution of a differential equation.

For a first order differential equation \(\frac{dy}{dx} = h(x, y)\), a power series solution exists at a point \(x = x_{0}\) if \(h(x, y)\) is analytic at \(x_{0}\). The solution is given by

\[y(x) = y(x_{0}) + \sum_{n = 1}^{\infty} \frac{F_{n}(x_{0},b)(x - x_{0})^n}{n!},\]

where \(y(x_{0}) = b\) is the value of y at the initial value of \(x_{0}\). To compute the values of the \(F_{n}(x_{0},b)\) the following algorithm is followed, until the required number of terms are generated.

  1. \(F_1 = h(x_{0}, b)\)

  2. \(F_{n+1} = \frac{\partial F_{n}}{\partial x} + \frac{\partial F_{n}}{\partial y}F_{1}\)

References

  • Travis W. Walker, Analytic power series technique for solving first-order differential equations, p.p 17, 18

Examples

>>> from sympy import Function, Derivative, pprint, exp
>>> from sympy.solvers.ode import dsolve
>>> from sympy.abc import x
>>> f = Function('f')
>>> eq = exp(x)*(f(x).diff(x)) - f(x)
>>> pprint(dsolve(eq, hint='1st_power_series'))
                       3       4       5
                   C1*x    C1*x    C1*x     / 6\
f(x) = C1 + C1*x - ----- + ----- + ----- + O\x /
                     6       24      60

2nd_power_series_ordinary

sympy.solvers.ode.ode_2nd_power_series_ordinary(eq, func, order, match)[source]

Gives a power series solution to a second order homogeneous differential equation with polynomial coefficients at an ordinary point. A homogenous differential equation is of the form

\[P(x)\frac{d^2y}{dx^2} + Q(x)\frac{dy}{dx} + R(x) = 0\]

For simplicity it is assumed that \(P(x)\), \(Q(x)\) and \(R(x)\) are polynomials, it is sufficient that \(\frac{Q(x)}{P(x)}\) and \(\frac{R(x)}{P(x)}\) exists at \(x_{0}\). A recurrence relation is obtained by substituting \(y\) as \(\sum_{n=0}^\infty a_{n}x^{n}\), in the differential equation, and equating the nth term. Using this relation various terms can be generated.

References

Examples

>>> from sympy import dsolve, Function, pprint
>>> from sympy.abc import x, y
>>> f = Function("f")
>>> eq = f(x).diff(x, 2) + f(x)
>>> pprint(dsolve(eq, hint='2nd_power_series_ordinary'))
          / 4    2    \        /   2    \
          |x    x     |        |  x     |    / 6\
f(x) = C2*|-- - -- + 1| + C1*x*|- -- + 1| + O\x /
          \24   2     /        \  6     /

2nd_power_series_regular

sympy.solvers.ode.ode_2nd_power_series_regular(eq, func, order, match)[source]

Gives a power series solution to a second order homogeneous differential equation with polynomial coefficients at a regular point. A second order homogenous differential equation is of the form

\[P(x)\frac{d^2y}{dx^2} + Q(x)\frac{dy}{dx} + R(x) = 0\]

A point is said to regular singular at \(x0\) if \(x - x0\frac{Q(x)}{P(x)}\) and \((x - x0)^{2}\frac{R(x)}{P(x)}\) are analytic at \(x0\). For simplicity \(P(x)\), \(Q(x)\) and \(R(x)\) are assumed to be polynomials. The algorithm for finding the power series solutions is:

  1. Try expressing \((x - x0)P(x)\) and \(((x - x0)^{2})Q(x)\) as power series solutions about x0. Find \(p0\) and \(q0\) which are the constants of the power series expansions.

  2. Solve the indicial equation \(f(m) = m(m - 1) + m*p0 + q0\), to obtain the roots \(m1\) and \(m2\) of the indicial equation.

  3. If \(m1 - m2\) is a non integer there exists two series solutions. If \(m1 = m2\), there exists only one solution. If \(m1 - m2\) is an integer, then the existence of one solution is confirmed. The other solution may or may not exist.

The power series solution is of the form \(x^{m}\sum_{n=0}^\infty a_{n}x^{n}\). The coefficients are determined by the following recurrence relation. \(a_{n} = -\frac{\sum_{k=0}^{n-1} q_{n-k} + (m + k)p_{n-k}}{f(m + n)}\). For the case in which \(m1 - m2\) is an integer, it can be seen from the recurrence relation that for the lower root \(m\), when \(n\) equals the difference of both the roots, the denominator becomes zero. So if the numerator is not equal to zero, a second series solution exists.

References

  • George E. Simmons, “Differential Equations with Applications and Historical Notes”, p.p 176 - 184

Examples

>>> from sympy import dsolve, Function, pprint
>>> from sympy.abc import x, y
>>> f = Function("f")
>>> eq = x*(f(x).diff(x, 2)) + 2*(f(x).diff(x)) + x*f(x)
>>> pprint(dsolve(eq))
                              /    6    4    2    \
                              |   x    x    x     |
          /  4    2    \   C1*|- --- + -- - -- + 1|
          | x    x     |      \  720   24   2     /    / 6\
f(x) = C2*|--- - -- + 1| + ------------------------ + O\x /
          \120   6     /              x

Lie heuristics

These functions are intended for internal use of the Lie Group Solver. Nonetheless, they contain useful information in their docstrings on the algorithms implemented for the various heuristics.

abaco1_simple

sympy.solvers.ode.lie_heuristic_abaco1_simple(match, comp=False)[source]

The first heuristic uses the following four sets of assumptions on \(\xi\) and \(\eta\)

\[\xi = 0, \eta = f(x)\]
\[\xi = 0, \eta = f(y)\]
\[\xi = f(x), \eta = 0\]
\[\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

\[\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

abaco1_product

sympy.solvers.ode.lie_heuristic_abaco1_product(match, comp=False)[source]

The second heuristic uses the following two assumptions on \(\xi\) and \(\eta\)

\[\eta = 0, \xi = f(x)*g(y)\]
\[\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

\[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

bivariate

sympy.solvers.ode.lie_heuristic_bivariate(match, comp=False)[source]

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

chi

sympy.solvers.ode.lie_heuristic_chi(match, comp=False)[source]

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

abaco2_similar

sympy.solvers.ode.lie_heuristic_abaco2_similar(match, comp=False)[source]

This heuristic uses the following two assumptions on \(\xi\) and \(\eta\)

\[\eta = g(x), \xi = f(x)\]
\[\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

function_sum

sympy.solvers.ode.lie_heuristic_function_sum(match, comp=False)[source]

This heuristic uses the following two assumptions on \(\xi\) and \(\eta\)

\[\eta = 0, \xi = f(x) + g(y)\]
\[\eta = f(x) + g(y), \xi = 0\]

The first assumption of this heuristic holds good if

\[\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

abaco2_unique_unknown

sympy.solvers.ode.lie_heuristic_abaco2_unique_unknown(match, comp=False)[source]

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] \(\xi = -R,\eta = 1\)

    b] \(\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

abaco2_unique_general

sympy.solvers.ode.lie_heuristic_abaco2_unique_general(match, comp=False)[source]

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

linear

sympy.solvers.ode.lie_heuristic_linear(match, comp=False)[source]

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

\[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

System of ODEs

These functions are intended for internal use by dsolve() for system of differential equations.

system_of_odes_linear_2eq_order1_type1

sympy.solvers.ode._linear_2eq_order1_type1(x, y, t, r, eq)[source]

It is classified under system of two linear homogeneous first-order constant-coefficient ordinary differential equations.

The equations which come under this type are

\[x' = ax + by,\]
\[y' = cx + dy\]

The characteristics equation is written as

\[\lambda^{2} + (a+d) \lambda + ad - bc = 0\]

and its discriminant is \(D = (a-d)^{2} + 4bc\). There are several cases

1. Case when \(ad - bc \neq 0\). The origin of coordinates, \(x = y = 0\), is the only stationary point; it is - a node if \(D = 0\) - a node if \(D > 0\) and \(ad - bc > 0\) - a saddle if \(D > 0\) and \(ad - bc < 0\) - a focus if \(D < 0\) and \(a + d \neq 0\) - a centre if \(D < 0\) and \(a + d \neq 0\).

1.1. If \(D > 0\). The characteristic equation has two distinct real roots \(\lambda_1\) and \(\lambda_ 2\) . The general solution of the system in question is expressed as

\[x = C_1 b e^{\lambda_1 t} + C_2 b e^{\lambda_2 t}\]
\[y = C_1 (\lambda_1 - a) e^{\lambda_1 t} + C_2 (\lambda_2 - a) e^{\lambda_2 t}\]

where \(C_1\) and \(C_2\) being arbitrary constants

1.2. If \(D < 0\). The characteristics equation has two conjugate roots, \(\lambda_1 = \sigma + i \beta\) and \(\lambda_2 = \sigma - i \beta\). The general solution of the system is given by

\[x = b e^{\sigma t} (C_1 \sin(\beta t) + C_2 \cos(\beta t))\]
\[y = e^{\sigma t} ([(\sigma - a) C_1 - \beta C_2] \sin(\beta t) + [\beta C_1 + (\sigma - a) C_2 \cos(\beta t)])\]

1.3. If \(D = 0\) and \(a \neq d\). The characteristic equation has two equal roots, \(\lambda_1 = \lambda_2\). The general solution of the system is written as

\[x = 2b (C_1 + \frac{C_2}{a-d} + C_2 t) e^{\frac{a+d}{2} t}\]
\[y = [(d - a) C_1 + C_2 + (d - a) C_2 t] e^{\frac{a+d}{2} t}\]

1.4. If \(D = 0\) and \(a = d \neq 0\) and \(b = 0\)

\[x = C_1 e^{a t} , y = (c C_1 t + C_2) e^{a t}\]

1.5. If \(D = 0\) and \(a = d \neq 0\) and \(c = 0\)

\[x = (b C_1 t + C_2) e^{a t} , y = C_1 e^{a t}\]

2. Case when \(ad - bc = 0\) and \(a^{2} + b^{2} > 0\). The whole straight line \(ax + by = 0\) consists of singular points. The original system of differential equations can be rewritten as

\[x' = ax + by , y' = k (ax + by)\]

2.1 If \(a + bk \neq 0\), solution will be

\[x = b C_1 + C_2 e^{(a + bk) t} , y = -a C_1 + k C_2 e^{(a + bk) t}\]

2.2 If \(a + bk = 0\), solution will be

\[x = C_1 (bk t - 1) + b C_2 t , y = k^{2} b C_1 t + (b k^{2} t + 1) C_2\]

system_of_odes_linear_2eq_order1_type2

sympy.solvers.ode._linear_2eq_order1_type2(x, y, t, r, eq)[source]

The equations of this type are

\[x' = ax + by + k1 , y' = cx + dy + k2\]

The general solution of this system is given by sum of its particular solution and the general solution of the corresponding homogeneous system is obtained from type1.

1. When \(ad - bc \neq 0\). The particular solution will be \(x = x_0\) and \(y = y_0\) where \(x_0\) and \(y_0\) are determined by solving linear system of equations

\[a x_0 + b y_0 + k1 = 0 , c x_0 + d y_0 + k2 = 0\]
  1. When \(ad - bc = 0\) and \(a^{2} + b^{2} > 0\). In this case, the system of equation becomes

\[x' = ax + by + k_1 , y' = k (ax + by) + k_2\]

2.1 If \(\sigma = a + bk \neq 0\), particular solution is given by

\[x = b \sigma^{-1} (c_1 k - c_2) t - \sigma^{-2} (a c_1 + b c_2)\]
\[y = kx + (c_2 - c_1 k) t\]

2.2 If \(\sigma = a + bk = 0\), particular solution is given by

\[x = \frac{1}{2} b (c_2 - c_1 k) t^{2} + c_1 t\]
\[y = kx + (c_2 - c_1 k) t\]

system_of_odes_linear_2eq_order1_type3

sympy.solvers.ode._linear_2eq_order1_type3(x, y, t, r, eq)[source]

The equations of this type of ode are

\[x' = f(t) x + g(t) y\]
\[y' = g(t) x + f(t) y\]

The solution of such equations is given by

\[x = e^{F} (C_1 e^{G} + C_2 e^{-G}) , y = e^{F} (C_1 e^{G} - C_2 e^{-G})\]

where \(C_1\) and \(C_2\) are arbitrary constants, and

\[F = \int f(t) \,dt , G = \int g(t) \,dt\]

system_of_odes_linear_2eq_order1_type4

sympy.solvers.ode._linear_2eq_order1_type4(x, y, t, r, eq)[source]

The equations of this type of ode are .

\[x' = f(t) x + g(t) y\]
\[y' = -g(t) x + f(t) y\]

The solution is given by

\[x = F (C_1 \cos(G) + C_2 \sin(G)), y = F (-C_1 \sin(G) + C_2 \cos(G))\]

where \(C_1\) and \(C_2\) are arbitrary constants, and

\[F = \int f(t) \,dt , G = \int g(t) \,dt\]

system_of_odes_linear_2eq_order1_type5

sympy.solvers.ode._linear_2eq_order1_type5(x, y, t, r, eq)[source]

The equations of this type of ode are .

\[x' = f(t) x + g(t) y\]
\[y' = a g(t) x + [f(t) + b g(t)] y\]

The transformation of

\[x = e^{\int f(t) \,dt} u , y = e^{\int f(t) \,dt} v , T = \int g(t) \,dt\]

leads to a system of constant coefficient linear differential equations

\[u'(T) = v , v'(T) = au + bv\]

system_of_odes_linear_2eq_order1_type6

sympy.solvers.ode._linear_2eq_order1_type6(x, y, t, r, eq)[source]

The equations of this type of ode are .

\[x' = f(t) x + g(t) y\]
\[y' = a [f(t) + a h(t)] x + a [g(t) - h(t)] y\]

This is solved by first multiplying the first equation by \(-a\) and adding it to the second equation to obtain

\[y' - a x' = -a h(t) (y - a x)\]

Setting \(U = y - ax\) and integrating the equation we arrive at

\[y - ax = C_1 e^{-a \int h(t) \,dt}\]

and on substituting the value of y in first equation give rise to first order ODEs. After solving for \(x\), we can obtain \(y\) by substituting the value of \(x\) in second equation.

system_of_odes_linear_2eq_order1_type7

sympy.solvers.ode._linear_2eq_order1_type7(x, y, t, r, eq)[source]

The equations of this type of ode are .

\[x' = f(t) x + g(t) y\]
\[y' = h(t) x + p(t) y\]

Differentiating the first equation and substituting the value of \(y\) from second equation will give a second-order linear equation

\[g x'' - (fg + gp + g') x' + (fgp - g^{2} h + f g' - f' g) x = 0\]

This above equation can be easily integrated if following conditions are satisfied.

  1. \(fgp - g^{2} h + f g' - f' g = 0\)

  2. \(fgp - g^{2} h + f g' - f' g = ag, fg + gp + g' = bg\)

If first condition is satisfied then it is solved by current dsolve solver and in second case it becomes a constant coefficient differential equation which is also solved by current solver.

Otherwise if the above condition fails then, a particular solution is assumed as \(x = x_0(t)\) and \(y = y_0(t)\) Then the general solution is expressed as

\[x = C_1 x_0(t) + C_2 x_0(t) \int \frac{g(t) F(t) P(t)}{x_0^{2}(t)} \,dt\]
\[y = C_1 y_0(t) + C_2 [\frac{F(t) P(t)}{x_0(t)} + y_0(t) \int \frac{g(t) F(t) P(t)}{x_0^{2}(t)} \,dt]\]

where C1 and C2 are arbitrary constants and

\[F(t) = e^{\int f(t) \,dt} , P(t) = e^{\int p(t) \,dt}\]

system_of_odes_linear_2eq_order2_type1

sympy.solvers.ode._linear_2eq_order2_type1(x, y, t, r, eq)[source]

System of two constant-coefficient second-order linear homogeneous differential equations

\[x'' = ax + by\]
\[y'' = cx + dy\]

The characteristic equation for above equations

\[\lambda^4 - (a + d) \lambda^2 + ad - bc = 0\]

whose discriminant is \(D = (a - d)^2 + 4bc \neq 0\)

  1. When \(ad - bc \neq 0\)

1.1. If \(D \neq 0\). The characteristic equation has four distinct roots, \(\lambda_1, \lambda_2, \lambda_3, \lambda_4\). The general solution of the system is

\[x = C_1 b e^{\lambda_1 t} + C_2 b e^{\lambda_2 t} + C_3 b e^{\lambda_3 t} + C_4 b e^{\lambda_4 t}\]
\[y = C_1 (\lambda_1^{2} - a) e^{\lambda_1 t} + C_2 (\lambda_2^{2} - a) e^{\lambda_2 t} + C_3 (\lambda_3^{2} - a) e^{\lambda_3 t} + C_4 (\lambda_4^{2} - a) e^{\lambda_4 t}\]

where \(C_1,..., C_4\) are arbitrary constants.

1.2. If \(D = 0\) and \(a \neq d\):

\[x = 2 C_1 (bt + \frac{2bk}{a - d}) e^{\frac{kt}{2}} + 2 C_2 (bt + \frac{2bk}{a - d}) e^{\frac{-kt}{2}} + 2b C_3 t e^{\frac{kt}{2}} + 2b C_4 t e^{\frac{-kt}{2}}\]
\[y = C_1 (d - a) t e^{\frac{kt}{2}} + C_2 (d - a) t e^{\frac{-kt}{2}} + C_3 [(d - a) t + 2k] e^{\frac{kt}{2}} + C_4 [(d - a) t - 2k] e^{\frac{-kt}{2}}\]

where \(C_1,..., C_4\) are arbitrary constants and \(k = \sqrt{2 (a + d)}\)

1.3. If \(D = 0\) and \(a = d \neq 0\) and \(b = 0\):

\[x = 2 \sqrt{a} C_1 e^{\sqrt{a} t} + 2 \sqrt{a} C_2 e^{-\sqrt{a} t}\]
\[y = c C_1 t e^{\sqrt{a} t} - c C_2 t e^{-\sqrt{a} t} + C_3 e^{\sqrt{a} t} + C_4 e^{-\sqrt{a} t}\]

1.4. If \(D = 0\) and \(a = d \neq 0\) and \(c = 0\):

\[x = b C_1 t e^{\sqrt{a} t} - b C_2 t e^{-\sqrt{a} t} + C_3 e^{\sqrt{a} t} + C_4 e^{-\sqrt{a} t}\]
\[y = 2 \sqrt{a} C_1 e^{\sqrt{a} t} + 2 \sqrt{a} C_2 e^{-\sqrt{a} t}\]
  1. When \(ad - bc = 0\) and \(a^2 + b^2 > 0\). Then the original system becomes

\[x'' = ax + by\]
\[y'' = k (ax + by)\]

2.1. If \(a + bk \neq 0\):

\[x = C_1 e^{t \sqrt{a + bk}} + C_2 e^{-t \sqrt{a + bk}} + C_3 bt + C_4 b\]
\[y = C_1 k e^{t \sqrt{a + bk}} + C_2 k e^{-t \sqrt{a + bk}} - C_3 at - C_4 a\]

2.2. If \(a + bk = 0\):

\[x = C_1 b t^3 + C_2 b t^2 + C_3 t + C_4\]
\[y = kx + 6 C_1 t + 2 C_2\]

system_of_odes_linear_2eq_order2_type2

sympy.solvers.ode._linear_2eq_order2_type2(x, y, t, r, eq)[source]

The equations in this type are

\[x'' = a_1 x + b_1 y + c_1\]
\[y'' = a_2 x + b_2 y + c_2\]

The general solution of this system is given by the sum of its particular solution and the general solution of the homogeneous system. The general solution is given by the linear system of 2 equation of order 2 and type 1

1. If \(a_1 b_2 - a_2 b_1 \neq 0\). A particular solution will be \(x = x_0\) and \(y = y_0\) where the constants \(x_0\) and \(y_0\) are determined by solving the linear algebraic system

\[a_1 x_0 + b_1 y_0 + c_1 = 0, a_2 x_0 + b_2 y_0 + c_2 = 0\]
  1. If \(a_1 b_2 - a_2 b_1 = 0\) and \(a_1^2 + b_1^2 > 0\). In this case, the system in question becomes

\[x'' = ax + by + c_1, y'' = k (ax + by) + c_2\]

2.1. If \(\sigma = a + bk \neq 0\), the particular solution will be

\[x = \frac{1}{2} b \sigma^{-1} (c_1 k - c_2) t^2 - \sigma^{-2} (a c_1 + b c_2)\]
\[y = kx + \frac{1}{2} (c_2 - c_1 k) t^2\]

2.2. If \(\sigma = a + bk = 0\), the particular solution will be

\[x = \frac{1}{24} b (c_2 - c_1 k) t^4 + \frac{1}{2} c_1 t^2\]
\[y = kx + \frac{1}{2} (c_2 - c_1 k) t^2\]

system_of_odes_linear_2eq_order2_type3

sympy.solvers.ode._linear_2eq_order2_type3(x, y, t, r, eq)[source]

These type of equation is used for describing the horizontal motion of a pendulum taking into account the Earth rotation. The solution is given with \(a^2 + 4b > 0\):

\[x = C_1 \cos(\alpha t) + C_2 \sin(\alpha t) + C_3 \cos(\beta t) + C_4 \sin(\beta t)\]
\[y = -C_1 \sin(\alpha t) + C_2 \cos(\alpha t) - C_3 \sin(\beta t) + C_4 \cos(\beta t)\]

where \(C_1,...,C_4\) and

\[\alpha = \frac{1}{2} a + \frac{1}{2} \sqrt{a^2 + 4b}, \beta = \frac{1}{2} a - \frac{1}{2} \sqrt{a^2 + 4b}\]

system_of_odes_linear_2eq_order2_type4

sympy.solvers.ode._linear_2eq_order2_type4(x, y, t, r, eq)[source]

These equations are found in the theory of oscillations

\[x'' + a_1 x' + b_1 y' + c_1 x + d_1 y = k_1 e^{i \omega t}\]
\[y'' + a_2 x' + b_2 y' + c_2 x + d_2 y = k_2 e^{i \omega t}\]

The general solution of this linear nonhomogeneous system of constant-coefficient differential equations is given by the sum of its particular solution and the general solution of the corresponding homogeneous system (with \(k_1 = k_2 = 0\))

  1. A particular solution is obtained by the method of undetermined coefficients:

\[x = A_* e^{i \omega t}, y = B_* e^{i \omega t}\]

On substituting these expressions into the original system of differential equations, one arrive at a linear nonhomogeneous system of algebraic equations for the coefficients \(A\) and \(B\).

2. The general solution of the homogeneous system of differential equations is determined by a linear combination of linearly independent particular solutions determined by the method of undetermined coefficients in the form of exponentials:

\[x = A e^{\lambda t}, y = B e^{\lambda t}\]

On substituting these expressions into the original system and collecting the coefficients of the unknown \(A\) and \(B\), one obtains

\[(\lambda^{2} + a_1 \lambda + c_1) A + (b_1 \lambda + d_1) B = 0\]
\[(a_2 \lambda + c_2) A + (\lambda^{2} + b_2 \lambda + d_2) B = 0\]

The determinant of this system must vanish for nontrivial solutions A, B to exist. This requirement results in the following characteristic equation for \(\lambda\)

\[(\lambda^2 + a_1 \lambda + c_1) (\lambda^2 + b_2 \lambda + d_2) - (b_1 \lambda + d_1) (a_2 \lambda + c_2) = 0\]

If all roots \(k_1,...,k_4\) of this equation are distinct, the general solution of the original system of the differential equations has the form

\[x = C_1 (b_1 \lambda_1 + d_1) e^{\lambda_1 t} - C_2 (b_1 \lambda_2 + d_1) e^{\lambda_2 t} - C_3 (b_1 \lambda_3 + d_1) e^{\lambda_3 t} - C_4 (b_1 \lambda_4 + d_1) e^{\lambda_4 t}\]
\[y = C_1 (\lambda_1^{2} + a_1 \lambda_1 + c_1) e^{\lambda_1 t} + C_2 (\lambda_2^{2} + a_1 \lambda_2 + c_1) e^{\lambda_2 t} + C_3 (\lambda_3^{2} + a_1 \lambda_3 + c_1) e^{\lambda_3 t} + C_4 (\lambda_4^{2} + a_1 \lambda_4 + c_1) e^{\lambda_4 t}\]

system_of_odes_linear_2eq_order2_type5

sympy.solvers.ode._linear_2eq_order2_type5(x, y, t, r, eq)[source]

The equation which come under this category are

\[x'' = a (t y' - y)\]
\[y'' = b (t x' - x)\]

The transformation

\[u = t x' - x, b = t y' - y\]

leads to the first-order system

\[u' = atv, v' = btu\]

The general solution of this system is given by

If \(ab > 0\):

\[u = C_1 a e^{\frac{1}{2} \sqrt{ab} t^2} + C_2 a e^{-\frac{1}{2} \sqrt{ab} t^2}\]
\[v = C_1 \sqrt{ab} e^{\frac{1}{2} \sqrt{ab} t^2} - C_2 \sqrt{ab} e^{-\frac{1}{2} \sqrt{ab} t^2}\]

If \(ab < 0\):

\[u = C_1 a \cos(\frac{1}{2} \sqrt{\left|ab\right|} t^2) + C_2 a \sin(-\frac{1}{2} \sqrt{\left|ab\right|} t^2)\]
\[v = C_1 \sqrt{\left|ab\right|} \sin(\frac{1}{2} \sqrt{\left|ab\right|} t^2) + C_2 \sqrt{\left|ab\right|} \cos(-\frac{1}{2} \sqrt{\left|ab\right|} t^2)\]

where \(C_1\) and \(C_2\) are arbitrary constants. On substituting the value of \(u\) and \(v\) in above equations and integrating the resulting expressions, the general solution will become

\[x = C_3 t + t \int \frac{u}{t^2} \,dt, y = C_4 t + t \int \frac{u}{t^2} \,dt\]

where \(C_3\) and \(C_4\) are arbitrary constants.

system_of_odes_linear_2eq_order2_type6

sympy.solvers.ode._linear_2eq_order2_type6(x, y, t, r, eq)[source]

The equations are

\[x'' = f(t) (a_1 x + b_1 y)\]
\[y'' = f(t) (a_2 x + b_2 y)\]

If \(k_1\) and \(k_2\) are roots of the quadratic equation

\[k^2 - (a_1 + b_2) k + a_1 b_2 - a_2 b_1 = 0\]

Then by multiplying appropriate constants and adding together original equations we obtain two independent equations:

\[z_1'' = k_1 f(t) z_1, z_1 = a_2 x + (k_1 - a_1) y\]
\[z_2'' = k_2 f(t) z_2, z_2 = a_2 x + (k_2 - a_1) y\]

Solving the equations will give the values of \(x\) and \(y\) after obtaining the value of \(z_1\) and \(z_2\) by solving the differential equation and substituting the result.

system_of_odes_linear_2eq_order2_type7

sympy.solvers.ode._linear_2eq_order2_type7(x, y, t, r, eq)[source]

The equations are given as

\[x'' = f(t) (a_1 x' + b_1 y')\]
\[y'' = f(t) (a_2 x' + b_2 y')\]

If \(k_1\) and ‘k_2` are roots of the quadratic equation

\[k^2 - (a_1 + b_2) k + a_1 b_2 - a_2 b_1 = 0\]

Then the system can be reduced by adding together the two equations multiplied by appropriate constants give following two independent equations:

\[z_1'' = k_1 f(t) z_1', z_1 = a_2 x + (k_1 - a_1) y\]
\[z_2'' = k_2 f(t) z_2', z_2 = a_2 x + (k_2 - a_1) y\]

Integrating these and returning to the original variables, one arrives at a linear algebraic system for the unknowns \(x\) and \(y\):

\[a_2 x + (k_1 - a_1) y = C_1 \int e^{k_1 F(t)} \,dt + C_2\]
\[a_2 x + (k_2 - a_1) y = C_3 \int e^{k_2 F(t)} \,dt + C_4\]

where \(C_1,...,C_4\) are arbitrary constants and \(F(t) = \int f(t) \,dt\)

system_of_odes_linear_2eq_order2_type8

sympy.solvers.ode._linear_2eq_order2_type8(x, y, t, r, eq)[source]

The equation of this category are

\[x'' = a f(t) (t y' - y)\]
\[y'' = b f(t) (t x' - x)\]

The transformation

\[u = t x' - x, v = t y' - y\]

leads to the system of first-order equations

\[u' = a t f(t) v, v' = b t f(t) u\]

The general solution of this system has the form

If \(ab > 0\):

\[u = C_1 a e^{\sqrt{ab} \int t f(t) \,dt} + C_2 a e^{-\sqrt{ab} \int t f(t) \,dt}\]
\[v = C_1 \sqrt{ab} e^{\sqrt{ab} \int t f(t) \,dt} - C_2 \sqrt{ab} e^{-\sqrt{ab} \int t f(t) \,dt}\]

If \(ab < 0\):

\[u = C_1 a \cos(\sqrt{\left|ab\right|} \int t f(t) \,dt) + C_2 a \sin(-\sqrt{\left|ab\right|} \int t f(t) \,dt)\]
\[v = C_1 \sqrt{\left|ab\right|} \sin(\sqrt{\left|ab\right|} \int t f(t) \,dt) + C_2 \sqrt{\left|ab\right|} \cos(-\sqrt{\left|ab\right|} \int t f(t) \,dt)\]

where \(C_1\) and \(C_2\) are arbitrary constants. On substituting the value of \(u\) and \(v\) in above equations and integrating the resulting expressions, the general solution will become

\[x = C_3 t + t \int \frac{u}{t^2} \,dt, y = C_4 t + t \int \frac{u}{t^2} \,dt\]

where \(C_3\) and \(C_4\) are arbitrary constants.

system_of_odes_linear_2eq_order2_type9

sympy.solvers.ode._linear_2eq_order2_type9(x, y, t, r, eq)[source]
\[t^2 x'' + a_1 t x' + b_1 t y' + c_1 x + d_1 y = 0\]
\[t^2 y'' + a_2 t x' + b_2 t y' + c_2 x + d_2 y = 0\]

These system of equations are euler type.

The substitution of \(t = \sigma e^{\tau} (\sigma \neq 0)\) leads to the system of constant coefficient linear differential equations

\[x'' + (a_1 - 1) x' + b_1 y' + c_1 x + d_1 y = 0\]
\[y'' + a_2 x' + (b_2 - 1) y' + c_2 x + d_2 y = 0\]

The general solution of the homogeneous system of differential equations is determined by a linear combination of linearly independent particular solutions determined by the method of undetermined coefficients in the form of exponentials

\[x = A e^{\lambda t}, y = B e^{\lambda t}\]

On substituting these expressions into the original system and collecting the coefficients of the unknown \(A\) and \(B\), one obtains

\[(\lambda^{2} + (a_1 - 1) \lambda + c_1) A + (b_1 \lambda + d_1) B = 0\]
\[(a_2 \lambda + c_2) A + (\lambda^{2} + (b_2 - 1) \lambda + d_2) B = 0\]

The determinant of this system must vanish for nontrivial solutions A, B to exist. This requirement results in the following characteristic equation for \(\lambda\)

\[(\lambda^2 + (a_1 - 1) \lambda + c_1) (\lambda^2 + (b_2 - 1) \lambda + d_2) - (b_1 \lambda + d_1) (a_2 \lambda + c_2) = 0\]

If all roots \(k_1,...,k_4\) of this equation are distinct, the general solution of the original system of the differential equations has the form

\[x = C_1 (b_1 \lambda_1 + d_1) e^{\lambda_1 t} - C_2 (b_1 \lambda_2 + d_1) e^{\lambda_2 t} - C_3 (b_1 \lambda_3 + d_1) e^{\lambda_3 t} - C_4 (b_1 \lambda_4 + d_1) e^{\lambda_4 t}\]
\[y = C_1 (\lambda_1^{2} + (a_1 - 1) \lambda_1 + c_1) e^{\lambda_1 t} + C_2 (\lambda_2^{2} + (a_1 - 1) \lambda_2 + c_1) e^{\lambda_2 t} + C_3 (\lambda_3^{2} + (a_1 - 1) \lambda_3 + c_1) e^{\lambda_3 t} + C_4 (\lambda_4^{2} + (a_1 - 1) \lambda_4 + c_1) e^{\lambda_4 t}\]

system_of_odes_linear_2eq_order2_type10

sympy.solvers.ode._linear_2eq_order2_type10(x, y, t, r, eq)[source]

The equation of this category are

\[(\alpha t^2 + \beta t + \gamma)^{2} x'' = ax + by\]
\[(\alpha t^2 + \beta t + \gamma)^{2} y'' = cx + dy\]

The transformation

\[\tau = \int \frac{1}{\alpha t^2 + \beta t + \gamma} \,dt , u = \frac{x}{\sqrt{\left|\alpha t^2 + \beta t + \gamma\right|}} , v = \frac{y}{\sqrt{\left|\alpha t^2 + \beta t + \gamma\right|}}\]

leads to a constant coefficient linear system of equations

\[u'' = (a - \alpha \gamma + \frac{1}{4} \beta^{2}) u + b v\]
\[v'' = c u + (d - \alpha \gamma + \frac{1}{4} \beta^{2}) v\]

These system of equations obtained can be solved by type1 of System of two constant-coefficient second-order linear homogeneous differential equations.

system_of_odes_linear_2eq_order2_type11

sympy.solvers.ode._linear_2eq_order2_type11(x, y, t, r, eq)[source]

The equations which comes under this type are

\[x'' = f(t) (t x' - x) + g(t) (t y' - y)\]
\[y'' = h(t) (t x' - x) + p(t) (t y' - y)\]

The transformation

\[u = t x' - x, v = t y' - y\]

leads to the linear system of first-order equations

\[u' = t f(t) u + t g(t) v, v' = t h(t) u + t p(t) v\]

On substituting the value of \(u\) and \(v\) in transformed equation gives value of \(x\) and \(y\) as

\[x = C_3 t + t \int \frac{u}{t^2} \,dt , y = C_4 t + t \int \frac{v}{t^2} \,dt.\]

where \(C_3\) and \(C_4\) are arbitrary constants.

system_of_odes_linear_3eq_order1_type1

sympy.solvers.ode._linear_3eq_order1_type1(x, y, z, t, r, eq)[source]
\[x' = ax\]
\[y' = bx + cy\]
\[z' = dx + ky + pz\]

Solution of such equations are forward substitution. Solving first equations gives the value of \(x\), substituting it in second and third equation and solving second equation gives \(y\) and similarly substituting \(y\) in third equation give \(z\).

\[x = C_1 e^{at}\]
\[y = \frac{b C_1}{a - c} e^{at} + C_2 e^{ct}\]
\[z = \frac{C_1}{a - p} (d + \frac{bk}{a - c}) e^{at} + \frac{k C_2}{c - p} e^{ct} + C_3 e^{pt}\]

where \(C_1, C_2\) and \(C_3\) are arbitrary constants.

system_of_odes_linear_3eq_order1_type2

sympy.solvers.ode._linear_3eq_order1_type2(x, y, z, t, r, eq)[source]

The equations of this type are

\[x' = cy - bz\]
\[y' = az - cx\]
\[z' = bx - ay\]
  1. First integral:

\[ax + by + cz = A \qquad - (1)\]
\[x^2 + y^2 + z^2 = B^2 \qquad - (2)\]

where \(A\) and \(B\) are arbitrary constants. It follows from these integrals that the integral lines are circles formed by the intersection of the planes \((1)\) and sphere \((2)\)

  1. Solution:

\[x = a C_0 + k C_1 \cos(kt) + (c C_2 - b C_3) \sin(kt)\]
\[y = b C_0 + k C_2 \cos(kt) + (a C_2 - c C_3) \sin(kt)\]
\[z = c C_0 + k C_3 \cos(kt) + (b C_2 - a C_3) \sin(kt)\]

where \(k = \sqrt{a^2 + b^2 + c^2}\) and the four constants of integration, \(C_1,...,C_4\) are constrained by a single relation,

\[a C_1 + b C_2 + c C_3 = 0\]

system_of_odes_linear_3eq_order1_type3

sympy.solvers.ode._linear_3eq_order1_type3(x, y, z, t, r, eq)[source]

Equations of this system of ODEs

\[a x' = bc (y - z)\]
\[b y' = ac (z - x)\]
\[c z' = ab (x - y)\]
  1. First integral:

\[a^2 x + b^2 y + c^2 z = A\]

where A is an arbitrary constant. It follows that the integral lines are plane curves.

  1. Solution:

\[x = C_0 + k C_1 \cos(kt) + a^{-1} bc (C_2 - C_3) \sin(kt)\]
\[y = C_0 + k C_2 \cos(kt) + a b^{-1} c (C_3 - C_1) \sin(kt)\]
\[z = C_0 + k C_3 \cos(kt) + ab c^{-1} (C_1 - C_2) \sin(kt)\]

where \(k = \sqrt{a^2 + b^2 + c^2}\) and the four constants of integration, \(C_1,...,C_4\) are constrained by a single relation

\[a^2 C_1 + b^2 C_2 + c^2 C_3 = 0\]

system_of_odes_linear_3eq_order1_type4

sympy.solvers.ode._linear_3eq_order1_type4(x, y, z, t, r, eq)[source]

Equations:

\[x' = (a_1 f(t) + g(t)) x + a_2 f(t) y + a_3 f(t) z\]
\[y' = b_1 f(t) x + (b_2 f(t) + g(t)) y + b_3 f(t) z\]
\[z' = c_1 f(t) x + c_2 f(t) y + (c_3 f(t) + g(t)) z\]

The transformation

\[x = e^{\int g(t) \,dt} u, y = e^{\int g(t) \,dt} v, z = e^{\int g(t) \,dt} w, \tau = \int f(t) \,dt\]

leads to the system of constant coefficient linear differential equations

\[u' = a_1 u + a_2 v + a_3 w\]
\[v' = b_1 u + b_2 v + b_3 w\]
\[w' = c_1 u + c_2 v + c_3 w\]

These system of equations are solved by homogeneous linear system of constant coefficients of \(n\) equations of first order. Then substituting the value of \(u, v\) and \(w\) in transformed equation gives value of \(x, y\) and \(z\).

system_of_odes_linear_neq_order1_type1

sympy.solvers.ode._linear_neq_order1_type1(match_)[source]

System of n first-order constant-coefficient linear nonhomogeneous differential equation

\[y'_k = a_{k1} y_1 + a_{k2} y_2 +...+ a_{kn} y_n; k = 1,2,...,n\]

or that can be written as \(\vec{y'} = A . \vec{y}\) where \(\vec{y}\) is matrix of \(y_k\) for \(k = 1,2,...n\) and \(A\) is a \(n \times n\) matrix.

Since these equations are equivalent to a first order homogeneous linear differential equation. So the general solution will contain \(n\) linearly independent parts and solution will consist some type of exponential functions. Assuming \(y = \vec{v} e^{rt}\) is a solution of the system where \(\vec{v}\) is a vector of coefficients of \(y_1,...,y_n\). Substituting \(y\) and \(y' = r v e^{r t}\) into the equation \(\vec{y'} = A . \vec{y}\), we get

\[r \vec{v} e^{rt} = A \vec{v} e^{rt}\]
\[r \vec{v} = A \vec{v}\]

where \(r\) comes out to be eigenvalue of \(A\) and vector \(\vec{v}\) is the eigenvector of \(A\) corresponding to \(r\). There are three possibilities of eigenvalues of \(A\)

  • \(n\) distinct real eigenvalues

  • complex conjugate eigenvalues

  • eigenvalues with multiplicity \(k\)

1. When all eigenvalues \(r_1,..,r_n\) are distinct with \(n\) different eigenvectors \(v_1,...v_n\) then the solution is given by

\[\vec{y} = C_1 e^{r_1 t} \vec{v_1} + C_2 e^{r_2 t} \vec{v_2} +...+ C_n e^{r_n t} \vec{v_n}\]

where \(C_1,C_2,...,C_n\) are arbitrary constants.

2. When some eigenvalues are complex then in order to make the solution real, we take a linear combination: if \(r = a + bi\) has an eigenvector \(\vec{v} = \vec{w_1} + i \vec{w_2}\) then to obtain real-valued solutions to the system, replace the complex-valued solutions \(e^{rx} \vec{v}\) with real-valued solution \(e^{ax} (\vec{w_1} \cos(bx) - \vec{w_2} \sin(bx))\) and for \(r = a - bi\) replace the solution \(e^{-r x} \vec{v}\) with \(e^{ax} (\vec{w_1} \sin(bx) + \vec{w_2} \cos(bx))\)

3. If some eigenvalues are repeated. Then we get fewer than \(n\) linearly independent eigenvectors, we miss some of the solutions and need to construct the missing ones. We do this via generalized eigenvectors, vectors which are not eigenvectors but are close enough that we can use to write down the remaining solutions. For a eigenvalue \(r\) with eigenvector \(\vec{w}\) we obtain \(\vec{w_2},...,\vec{w_k}\) using

\[(A - r I) . \vec{w_2} = \vec{w}\]
\[(A - r I) . \vec{w_3} = \vec{w_2}\]
\[\vdots\]
\[(A - r I) . \vec{w_k} = \vec{w_{k-1}}\]

Then the solutions to the system for the eigenspace are \(e^{rt} [\vec{w}], e^{rt} [t \vec{w} + \vec{w_2}], e^{rt} [\frac{t^2}{2} \vec{w} + t \vec{w_2} + \vec{w_3}], ...,e^{rt} [\frac{t^{k-1}}{(k-1)!} \vec{w} + \frac{t^{k-2}}{(k-2)!} \vec{w_2} +...+ t \vec{w_{k-1}} + \vec{w_k}]\)

So, If \(\vec{y_1},...,\vec{y_n}\) are \(n\) solution of obtained from three categories of \(A\), then general solution to the system \(\vec{y'} = A . \vec{y}\)

\[\vec{y} = C_1 \vec{y_1} + C_2 \vec{y_2} + \cdots + C_n \vec{y_n}\]

system_of_odes_nonlinear_2eq_order1_type1

sympy.solvers.ode._nonlinear_2eq_order1_type1(x, y, t, eq)[source]

Equations:

\[x' = x^n F(x,y)\]
\[y' = g(y) F(x,y)\]

Solution:

\[x = \varphi(y), \int \frac{1}{g(y) F(\varphi(y),y)} \,dy = t + C_2\]

where

if \(n \neq 1\)

\[\varphi = [C_1 + (1-n) \int \frac{1}{g(y)} \,dy]^{\frac{1}{1-n}}\]

if \(n = 1\)

\[\varphi = C_1 e^{\int \frac{1}{g(y)} \,dy}\]

where \(C_1\) and \(C_2\) are arbitrary constants.

system_of_odes_nonlinear_2eq_order1_type2

sympy.solvers.ode._nonlinear_2eq_order1_type2(x, y, t, eq)[source]

Equations:

\[x' = e^{\lambda x} F(x,y)\]
\[y' = g(y) F(x,y)\]

Solution:

\[x = \varphi(y), \int \frac{1}{g(y) F(\varphi(y),y)} \,dy = t + C_2\]

where

if \(\lambda \neq 0\)

\[\varphi = -\frac{1}{\lambda} log(C_1 - \lambda \int \frac{1}{g(y)} \,dy)\]

if \(\lambda = 0\)

\[\varphi = C_1 + \int \frac{1}{g(y)} \,dy\]

where \(C_1\) and \(C_2\) are arbitrary constants.

system_of_odes_nonlinear_2eq_order1_type3

sympy.solvers.ode._nonlinear_2eq_order1_type3(x, y, t, eq)[source]

Autonomous system of general form

\[x' = F(x,y)\]
\[y' = G(x,y)\]

Assuming \(y = y(x, C_1)\) where \(C_1\) is an arbitrary constant is the general solution of the first-order equation

\[F(x,y) y'_x = G(x,y)\]

Then the general solution of the original system of equations has the form

\[\int \frac{1}{F(x,y(x,C_1))} \,dx = t + C_1\]

system_of_odes_nonlinear_2eq_order1_type4

sympy.solvers.ode._nonlinear_2eq_order1_type4(x, y, t, eq)[source]

Equation:

\[x' = f_1(x) g_1(y) \phi(x,y,t)\]
\[y' = f_2(x) g_2(y) \phi(x,y,t)\]

First integral:

\[\int \frac{f_2(x)}{f_1(x)} \,dx - \int \frac{g_1(y)}{g_2(y)} \,dy = C\]

where \(C\) is an arbitrary constant.

On solving the first integral for \(x\) (resp., \(y\) ) and on substituting the resulting expression into either equation of the original solution, one arrives at a first-order equation for determining \(y\) (resp., \(x\) ).

system_of_odes_nonlinear_2eq_order1_type5

sympy.solvers.ode._nonlinear_2eq_order1_type5(func, t, eq)[source]

Clairaut system of ODEs

\[x = t x' + F(x',y')\]
\[y = t y' + G(x',y')\]

The following are solutions of the system

\((i)\) straight lines:

\[x = C_1 t + F(C_1, C_2), y = C_2 t + G(C_1, C_2)\]

where \(C_1\) and \(C_2\) are arbitrary constants;

\((ii)\) envelopes of the above lines;

\((iii)\) continuously differentiable lines made up from segments of the lines \((i)\) and \((ii)\).

system_of_odes_nonlinear_3eq_order1_type1

sympy.solvers.ode._nonlinear_3eq_order1_type1(x, y, z, t, eq)[source]

Equations:

\[a x' = (b - c) y z, \enspace b y' = (c - a) z x, \enspace c z' = (a - b) x y\]

First Integrals:

\[a x^{2} + b y^{2} + c z^{2} = C_1\]
\[a^{2} x^{2} + b^{2} y^{2} + c^{2} z^{2} = C_2\]

where \(C_1\) and \(C_2\) are arbitrary constants. On solving the integrals for \(y\) and \(z\) and on substituting the resulting expressions into the first equation of the system, we arrives at a separable first-order equation on \(x\). Similarly doing that for other two equations, we will arrive at first order equation on \(y\) and \(z\) too.

References

-http://eqworld.ipmnet.ru/en/solutions/sysode/sode0401.pdf

system_of_odes_nonlinear_3eq_order1_type2

sympy.solvers.ode._nonlinear_3eq_order1_type2(x, y, z, t, eq)[source]

Equations:

\[a x' = (b - c) y z f(x, y, z, t)\]
\[b y' = (c - a) z x f(x, y, z, t)\]
\[c z' = (a - b) x y f(x, y, z, t)\]

First Integrals:

\[a x^{2} + b y^{2} + c z^{2} = C_1\]
\[a^{2} x^{2} + b^{2} y^{2} + c^{2} z^{2} = C_2\]

where \(C_1\) and \(C_2\) are arbitrary constants. On solving the integrals for \(y\) and \(z\) and on substituting the resulting expressions into the first equation of the system, we arrives at a first-order differential equations on \(x\). Similarly doing that for other two equations we will arrive at first order equation on \(y\) and \(z\).

References

-http://eqworld.ipmnet.ru/en/solutions/sysode/sode0402.pdf

system_of_odes_nonlinear_3eq_order1_type3

sympy.solvers.ode._nonlinear_3eq_order1_type3(x, y, z, t, eq)[source]

Equations:

\[x' = c F_2 - b F_3, \enspace y' = a F_3 - c F_1, \enspace z' = b F_1 - a F_2\]

where \(F_n = F_n(x, y, z, t)\).

  1. First Integral:

\[a x + b y + c z = C_1,\]

where C is an arbitrary constant.

2. If we assume function \(F_n\) to be independent of \(t\),i.e, \(F_n\) = \(F_n (x, y, z)\) Then, on eliminating \(t\) and \(z\) from the first two equation of the system, one arrives at the first-order equation

\[\frac{dy}{dx} = \frac{a F_3 (x, y, z) - c F_1 (x, y, z)}{c F_2 (x, y, z) - b F_3 (x, y, z)}\]

where \(z = \frac{1}{c} (C_1 - a x - b y)\)

References

-http://eqworld.ipmnet.ru/en/solutions/sysode/sode0404.pdf

system_of_odes_nonlinear_3eq_order1_type4

sympy.solvers.ode._nonlinear_3eq_order1_type4(x, y, z, t, eq)[source]

Equations:

\[x' = c z F_2 - b y F_3, \enspace y' = a x F_3 - c z F_1, \enspace z' = b y F_1 - a x F_2\]

where \(F_n = F_n (x, y, z, t)\)

  1. First integral:

\[a x^{2} + b y^{2} + c z^{2} = C_1\]

where \(C\) is an arbitrary constant.

2. Assuming the function \(F_n\) is independent of \(t\): \(F_n = F_n (x, y, z)\). Then on eliminating \(t\) and \(z\) from the first two equations of the system, one arrives at the first-order equation

\[\frac{dy}{dx} = \frac{a x F_3 (x, y, z) - c z F_1 (x, y, z)} {c z F_2 (x, y, z) - b y F_3 (x, y, z)}\]

where \(z = \pm \sqrt{\frac{1}{c} (C_1 - a x^{2} - b y^{2})}\)

References

-http://eqworld.ipmnet.ru/en/solutions/sysode/sode0405.pdf

system_of_odes_nonlinear_3eq_order1_type5

sympy.solvers.ode._nonlinear_3eq_order1_type5(x, y, t, eq)[source]
\[x' = x (c F_2 - b F_3), \enspace y' = y (a F_3 - c F_1), \enspace z' = z (b F_1 - a F_2)\]

where \(F_n = F_n (x, y, z, t)\) and are arbitrary functions.

First Integral:

\[\left|x\right|^{a} \left|y\right|^{b} \left|z\right|^{c} = C_1\]

where \(C\) is an arbitrary constant. If the function \(F_n\) is independent of \(t\), then, by eliminating \(t\) and \(z\) from the first two equations of the system, one arrives at a first-order equation.

References

-http://eqworld.ipmnet.ru/en/solutions/sysode/sode0406.pdf

Information on the ode module

This module contains dsolve() and different helper functions that it uses.

dsolve() solves ordinary differential equations. See the docstring on the various functions for their uses. Note that partial differential equations support is in pde.py. Note that hint functions have docstrings describing their various methods, but they are intended for internal use. Use dsolve(ode, func, hint=hint) to solve an ODE using a specific hint. See also the docstring on dsolve().

Functions in this module

These are the user functions in this module:

  • dsolve() - Solves ODEs.

  • classify_ode() - Classifies ODEs into possible hints for dsolve().

  • checkodesol() - Checks if an equation is the solution to an ODE.

  • homogeneous_order() - Returns the homogeneous order of an expression.

  • infinitesimals() - Returns the infinitesimals of the Lie group of point transformations of an ODE, such that it is invariant.

  • ode_checkinfsol() - Checks if the given infinitesimals are the actual infinitesimals of a first order ODE.

These are the non-solver helper functions that are for internal use. The user should use the various options to dsolve() to obtain the functionality provided by these functions:

See also the docstrings of these functions.

Currently implemented solver methods

The following methods are implemented for solving ordinary differential equations. See the docstrings of the various hint functions for more information on each (run help(ode)):

  • 1st order separable differential equations.

  • 1st order differential equations whose coefficients or \(dx\) and \(dy\) are functions homogeneous of the same order.

  • 1st order exact differential equations.

  • 1st order linear differential equations.

  • 1st order Bernoulli differential equations.

  • Power series solutions for first order differential equations.

  • Lie Group method of solving first order differential equations.

  • 2nd order Liouville differential equations.

  • Power series solutions for second order differential equations at ordinary and regular singular points.

  • \(n\)th order linear homogeneous differential equation with constant coefficients.

  • \(n\)th order linear inhomogeneous differential equation with constant coefficients using the method of undetermined coefficients.

  • \(n\)th order linear inhomogeneous differential equation with constant coefficients using the method of variation of parameters.

Philosophy behind this module

This module is designed to make it easy to add new ODE solving methods without having to mess with the solving code for other methods. The idea is that there is a classify_ode() function, which takes in an ODE and tells you what hints, if any, will solve the ODE. It does this without attempting to solve the ODE, so it is fast. Each solving method is a hint, and it has its own function, named ode_<hint>. That function takes in the ODE and any match expression gathered by classify_ode() and returns a solved result. If this result has any integrals in it, the hint function will return an unevaluated Integral class. dsolve(), which is the user wrapper function around all of this, will then call odesimp() on the result, which, among other things, will attempt to solve the equation for the dependent variable (the function we are solving for), simplify the arbitrary constants in the expression, and evaluate any integrals, if the hint allows it.

How to add new solution methods

If you have an ODE that you want dsolve() to be able to solve, try to avoid adding special case code here. Instead, try finding a general method that will solve your ODE, as well as others. This way, the ode module will become more robust, and unhindered by special case hacks. WolphramAlpha and Maple’s DETools[odeadvisor] function are two resources you can use to classify a specific ODE. It is also better for a method to work with an \(n\)th order ODE instead of only with specific orders, if possible.

To add a new method, there are a few things that you need to do. First, you need a hint name for your method. Try to name your hint so that it is unambiguous with all other methods, including ones that may not be implemented yet. If your method uses integrals, also include a hint_Integral hint. If there is more than one way to solve ODEs with your method, include a hint for each one, as well as a <hint>_best hint. Your ode_<hint>_best() function should choose the best using min with ode_sol_simplicity as the key argument. See ode_1st_homogeneous_coeff_best(), for example. The function that uses your method will be called ode_<hint>(), so the hint must only use characters that are allowed in a Python function name (alphanumeric characters and the underscore ‘_’ character). Include a function for every hint, except for _Integral hints (dsolve() takes care of those automatically). Hint names should be all lowercase, unless a word is commonly capitalized (such as Integral or Bernoulli). If you have a hint that you do not want to run with all_Integral that doesn’t have an _Integral counterpart (such as a best hint that would defeat the purpose of all_Integral), you will need to remove it manually in the dsolve() code. See also the classify_ode() docstring for guidelines on writing a hint name.

Determine in general how the solutions returned by your method compare with other methods that can potentially solve the same ODEs. Then, put your hints in the allhints tuple in the order that they should be called. The ordering of this tuple determines which hints are default. Note that exceptions are ok, because it is easy for the user to choose individual hints with dsolve(). In general, _Integral variants should go at the end of the list, and _best variants should go before the various hints they apply to. For example, the undetermined_coefficients hint comes before the variation_of_parameters hint because, even though variation of parameters is more general than undetermined coefficients, undetermined coefficients generally returns cleaner results for the ODEs that it can solve than variation of parameters does, and it does not require integration, so it is much faster.

Next, you need to have a match expression or a function that matches the type of the ODE, which you should put in classify_ode() (if the match function is more than just a few lines, like _undetermined_coefficients_match(), it should go outside of classify_ode()). It should match the ODE without solving for it as much as possible, so that classify_ode() remains fast and is not hindered by bugs in solving code. Be sure to consider corner cases. For example, if your solution method involves dividing by something, make sure you exclude the case where that division will be 0.

In most cases, the matching of the ODE will also give you the various parts that you need to solve it. You should put that in a dictionary (.match() will do this for you), and add that as matching_hints['hint'] = matchdict in the relevant part of classify_ode(). classify_ode() will then send this to dsolve(), which will send it to your function as the match argument. Your function should be named ode_<hint>(eq, func, order, match)`.  If you need to send more information, put it in the ``match dictionary. For example, if you had to substitute in a dummy variable in classify_ode() to match the ODE, you will need to pass it to your function using the \(match\) dict to access it. You can access the independent variable using func.args[0], and the dependent variable (the function you are trying to solve for) as func.func. If, while trying to solve the ODE, you find that you cannot, raise NotImplementedError. dsolve() will catch this error with the all meta-hint, rather than causing the whole routine to fail.

Add a docstring to your function that describes the method employed. Like with anything else in SymPy, you will need to add a doctest to the docstring, in addition to real tests in test_ode.py. Try to maintain consistency with the other hint functions’ docstrings. Add your method to the list at the top of this docstring. Also, add your method to ode.rst in the docs/src directory, so that the Sphinx docs will pull its docstring into the main SymPy documentation. Be sure to make the Sphinx documentation by running make html from within the doc directory to verify that the docstring formats correctly.

If your solution method involves integrating, use Integral() instead of integrate(). This allows the user to bypass hard/slow integration by using the _Integral variant of your hint. In most cases, calling sympy.core.basic.Basic.doit() will integrate your solution. If this is not the case, you will need to write special code in _handle_Integral(). Arbitrary constants should be symbols named C1, C2, and so on. All solution methods should return an equality instance. If you need an arbitrary number of arbitrary constants, you can use constants = numbered_symbols(prefix='C', cls=Symbol, start=1). If it is possible to solve for the dependent function in a general way, do so. Otherwise, do as best as you can, but do not call solve in your ode_<hint>() function. odesimp() will attempt to solve the solution for you, so you do not need to do that. Lastly, if your ODE has a common simplification that can be applied to your solutions, you can add a special case in odesimp() for it. For example, solutions returned from the 1st_homogeneous_coeff hints often have many log() terms, so odesimp() calls logcombine() on them (it also helps to write the arbitrary constant as log(C1) instead of C1 in this case). Also consider common ways that you can rearrange your solution to have constantsimp() take better advantage of it. It is better to put simplification in odesimp() than in your method, because it can then be turned off with the simplify flag in dsolve(). If you have any extraneous simplification in your function, be sure to only run it using if match.get('simplify', True):, especially if it can be slow or if it can reduce the domain of the solution.

Finally, as with every contribution to SymPy, your method will need to be tested. Add a test for each method in test_ode.py. Follow the conventions there, i.e., test the solver using dsolve(eq, f(x), hint=your_hint), and also test the solution using checkodesol() (you can put these in a separate tests and skip/XFAIL if it runs too slow/doesn’t work). Be sure to call your hint specifically in dsolve(), that way the test won’t be broken simply by the introduction of another matching hint. If your method works for higher order (>1) ODEs, you will need to run sol = constant_renumber(sol, 'C', 1, order) for each solution, where order is the order of the ODE. This is because constant_renumber renumbers the arbitrary constants by printing order, which is platform dependent. Try to test every corner case of your solver, including a range of orders if it is a \(n\)th order solver, but if your solver is slow, such as if it involves hard integration, try to keep the test run time down.

Feel free to refactor existing hints to avoid duplicating code or creating inconsistencies. If you can show that your method exactly duplicates an existing method, including in the simplicity and speed of obtaining the solutions, then you can remove the old, less general method. The existing code is tested extensively in test_ode.py, so if anything is broken, one of those tests will surely fail.