sympy/sympy

Switch branches/tags
Nothing to show
Fetching contributors…
Cannot retrieve contributors at this time
8444 lines (7089 sloc) 328 KB
 r""" This module contains :py:meth:~sympy.solvers.ode.dsolve and different helper functions that it uses. :py:meth:~sympy.solvers.ode.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 :py:meth:~sympy.solvers.ode.dsolve. **Functions in this module** These are the user functions in this module: - :py:meth:~sympy.solvers.ode.dsolve - Solves ODEs. - :py:meth:~sympy.solvers.ode.classify_ode - Classifies ODEs into possible hints for :py:meth:~sympy.solvers.ode.dsolve. - :py:meth:~sympy.solvers.ode.checkodesol - Checks if an equation is the solution to an ODE. - :py:meth:~sympy.solvers.ode.homogeneous_order - Returns the homogeneous order of an expression. - :py:meth:~sympy.solvers.ode.infinitesimals - Returns the infinitesimals of the Lie group of point transformations of an ODE, such that it is invariant. - :py:meth:~sympy.solvers.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 :py:meth:~sympy.solvers.ode.dsolve to obtain the functionality provided by these functions: - :py:meth:~sympy.solvers.ode.odesimp - Does all forms of ODE simplification. - :py:meth:~sympy.solvers.ode.ode_sol_simplicity - A key function for comparing solutions by simplicity. - :py:meth:~sympy.solvers.ode.constantsimp - Simplifies arbitrary constants. - :py:meth:~sympy.solvers.ode.constant_renumber - Renumber arbitrary constants. - :py:meth:~sympy.solvers.ode._handle_Integral - Evaluate unevaluated Integrals. 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 :py:meth:~sympy.solvers.ode.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_. That function takes in the ODE and any match expression gathered by :py:meth:~sympy.solvers.ode.classify_ode and returns a solved result. If this result has any integrals in it, the hint function will return an unevaluated :py:class:~sympy.integrals.Integral class. :py:meth:~sympy.solvers.ode.dsolve, which is the user wrapper function around all of this, will then call :py:meth:~sympy.solvers.ode.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 :py:meth:~sympy.solvers.ode.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 :py:mod:~sympy.solvers.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 _best hint. Your ode__best() function should choose the best using min with ode_sol_simplicity as the key argument. See :py:meth:~sympy.solvers.ode.ode_1st_homogeneous_coeff_best, for example. The function that uses your method will be called ode_(), 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 (:py:meth:~sympy.solvers.ode.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 :py:meth:~sympy.solvers.ode.dsolve code. See also the :py:meth:~sympy.solvers.ode.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 :py:data:~sympy.solvers.ode.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 :py:meth:~sympy.solvers.ode.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 :py:meth:~sympy.solvers.ode.classify_ode (if the match function is more than just a few lines, like :py:meth:~sympy.solvers.ode._undetermined_coefficients_match, it should go outside of :py:meth:~sympy.solvers.ode.classify_ode). It should match the ODE without solving for it as much as possible, so that :py:meth:~sympy.solvers.ode.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 :py:meth:~sympy.solvers.ode.classify_ode. :py:meth:~sympy.solvers.ode.classify_ode will then send this to :py:meth:~sympy.solvers.ode.dsolve, which will send it to your function as the match argument. Your function should be named ode_(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 :py:meth:~sympy.solvers.ode.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. :py:meth:~sympy.solvers.ode.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 :py:meth:Integral()  instead of :py:meth:~sympy.core.expr.Expr.integrate. This allows the user to bypass hard/slow integration by using the _Integral variant of your hint. In most cases, calling :py:meth:sympy.core.basic.Basic.doit will integrate your solution. If this is not the case, you will need to write special code in :py:meth:~sympy.solvers.ode._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_() function. :py:meth:~sympy.solvers.ode.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 :py:meth:~sympy.solvers.ode.odesimp for it. For example, solutions returned from the 1st_homogeneous_coeff hints often have many :py:meth:~sympy.functions.log terms, so :py:meth:~sympy.solvers.ode.odesimp calls :py:meth:~sympy.simplify.simplify.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 :py:meth:~sympy.solvers.ode.constantsimp take better advantage of it. It is better to put simplification in :py:meth:~sympy.solvers.ode.odesimp than in your method, because it can then be turned off with the simplify flag in :py:meth:~sympy.solvers.ode.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 :py:meth:~sympy.solvers.ode.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 :py:meth:~sympy.solvers.ode.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. """ from __future__ import print_function, division from collections import defaultdict from itertools import islice from sympy.core import Add, S, Mul, Pow, oo from sympy.core.compatibility import ordered, iterable, is_sequence, range from sympy.core.containers import Tuple from sympy.core.exprtools import factor_terms from sympy.core.expr import AtomicExpr, Expr from sympy.core.function import (Function, Derivative, AppliedUndef, diff, expand, expand_mul, Subs, _mexpand) from sympy.core.multidimensional import vectorize from sympy.core.numbers import NaN, zoo, I, Number from sympy.core.relational import Equality, Eq from sympy.core.symbol import Symbol, Wild, Dummy, symbols from sympy.core.sympify import sympify from sympy.logic.boolalg import BooleanAtom, And, Or, Not from sympy.functions import cos, exp, im, log, re, sin, tan, sqrt, \ atan2, conjugate, Piecewise from sympy.functions.combinatorial.factorials import factorial from sympy.integrals.integrals import Integral, integrate from sympy.matrices import wronskian, Matrix, eye, zeros from sympy.polys import (Poly, RootOf, rootof, terms_gcd, PolynomialError, lcm) from sympy.polys.polyroots import roots_quartic from sympy.polys.polytools import cancel, degree, div from sympy.series import Order from sympy.series.series import series from sympy.simplify import collect, logcombine, powsimp, separatevars, \ simplify, trigsimp, denom, posify, cse from sympy.simplify.powsimp import powdenest from sympy.simplify.radsimp import collect_const from sympy.solvers import solve from sympy.solvers.pde import pdsolve from sympy.utilities import numbered_symbols, default_sort_key, sift from sympy.solvers.deutils import _preprocess, ode_order, _desolve #: This is a list of hints in the order that they should be preferred by #: :py:meth:~sympy.solvers.ode.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 :py:meth:~sympy.solvers.ode.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. 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", ) lie_heuristics = ( "abaco1_simple", "abaco1_product", "abaco2_similar", "abaco2_unique_unknown", "abaco2_unique_general", "linear", "function_sum", "bivariate", "chi" ) def sub_func_doit(eq, func, new): r""" When replacing the func with something else, we usually want the derivative evaluated, so this function helps in making that happen. To keep subs from having to look through all derivatives, we mask them off with dummy variables, do the func sub, and then replace masked-off derivatives with their doit values. Examples ======== >>> from sympy import Derivative, symbols, Function >>> from sympy.solvers.ode import sub_func_doit >>> x, z = symbols('x, z') >>> y = Function('y') >>> sub_func_doit(3*Derivative(y(x), x) - 1, y(x), x) 2 >>> sub_func_doit(x*Derivative(y(x), x) - y(x)**2 + y(x), y(x), ... 1/(x*(z + 1/x))) x*(-1/(x**2*(z + 1/x)) + 1/(x**3*(z + 1/x)**2)) + 1/(x*(z + 1/x)) ...- 1/(x**2*(z + 1/x)**2) """ reps = {} repu = {} for d in eq.atoms(Derivative): u = Dummy('u') repu[u] = d.subs(func, new).doit() reps[d] = u # Make sure that expressions such as Derivative(f(x), (x, 2)) get # replaced before Derivative(f(x), x): reps = sorted(reps.items(), key=lambda x: -x[0].derivative_count) return eq.subs(reps).subs(func, new).subs(repu) def get_numbered_constants(eq, num=1, start=1, prefix='C'): """ Returns a list of constants that do not occur in eq already. """ if isinstance(eq, Expr): eq = [eq] elif not iterable(eq): raise ValueError("Expected Expr or iterable but got %s" % eq) atom_set = set().union(*[i.free_symbols for i in eq]) ncs = numbered_symbols(start=start, prefix=prefix, exclude=atom_set) Cs = [next(ncs) for i in range(num)] return (Cs[0] if num == 1 else tuple(Cs)) def dsolve(eq, func=None, hint="default", simplify=True, ics= None, xi=None, eta=None, x0=0, n=6, **kwargs): r""" Solves any (supported) kind of ordinary differential equation and system of ordinary differential equations. 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 :py:mod:~sympy.solvers.ode docstring for supported methods). This can either be an :py:class:~sympy.core.relational.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 :py:meth:~sympy.solvers.ode.classify_ode. See Hints below for more options that you can use for hint. simplify enables simplification by :py:meth:~sympy.solvers.ode.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 :py:meth:~sympy.solvers.ode.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 :py:meth:~sympy.solvers.ode.dsolve: default: This uses whatever hint is returned first by :py:meth:~sympy.solvers.ode.classify_ode. This is the default argument to :py:meth:~sympy.solvers.ode.dsolve. all: To make :py:meth:~sympy.solvers.ode.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 :py:meth:~sympy.solvers.deutils.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 :py:meth:~sympy.solvers.ode.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 :py:meth:~sympy.solvers.ode.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 :py:meth:~sympy.solvers.ode.dsolve to hang because of a difficult or impossible integral. This meta-hint will also be much faster than all, because :py:meth:~sympy.core.expr.Expr.integrate is an expensive routine. best: To have :py:meth:~sympy.solvers.ode.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 :py:meth:~sympy.solvers.ode.classify_ode docstring for more info on hints, and the :py:mod:~sympy.solvers.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 :py:meth:~sympy.solvers.ode.dsolve. - :py:meth:~sympy.solvers.ode.dsolve always returns an :py:class:~sympy.core.relational.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_) to get help more information on a specific hint, where  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 :py:class:~sympy.core.relational.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. 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)))} """ if iterable(eq): match = classify_sysode(eq, func) eq = match['eq'] order = match['order'] func = match['func'] t = list(list(eq[0].atoms(Derivative))[0].atoms(Symbol))[0] # keep highest order term coefficient positive for i in range(len(eq)): for func_ in func: if isinstance(func_, list): pass else: if eq[i].coeff(diff(func[i],t,ode_order(eq[i], func[i]))).is_negative: eq[i] = -eq[i] match['eq'] = eq if len(set(order.values()))!=1: raise ValueError("It solves only those systems of equations whose orders are equal") match['order'] = list(order.values())[0] def recur_len(l): return sum(recur_len(item) if isinstance(item,list) else 1 for item in l) if recur_len(func) != len(eq): raise ValueError("dsolve() and classify_sysode() work with " "number of functions being equal to number of equations") if match['type_of_equation'] is None: raise NotImplementedError else: if match['is_linear'] == True: if match['no_of_equation'] > 3: solvefunc = globals()['sysode_linear_neq_order%(order)s' % match] else: solvefunc = globals()['sysode_linear_%(no_of_equation)seq_order%(order)s' % match] else: solvefunc = globals()['sysode_nonlinear_%(no_of_equation)seq_order%(order)s' % match] sols = solvefunc(match) if ics: constants = Tuple(*sols).free_symbols - Tuple(*eq).free_symbols solved_constants = solve_ics(sols, func, constants, ics) return [sol.subs(solved_constants) for sol in sols] return sols else: given_hint = hint # hint given by the user # See the docstring of _desolve for more details. hints = _desolve(eq, func=func, hint=hint, simplify=True, xi=xi, eta=eta, type='ode', ics=ics, x0=x0, n=n, **kwargs) eq = hints.pop('eq', eq) all_ = hints.pop('all', False) if all_: retdict = {} failed_hints = {} gethints = classify_ode(eq, dict=True) orderedhints = gethints['ordered_hints'] for hint in hints: try: rv = _helper_simplify(eq, hint, hints[hint], simplify) except NotImplementedError as detail: failed_hints[hint] = detail else: retdict[hint] = rv func = hints[hint]['func'] retdict['best'] = min(list(retdict.values()), key=lambda x: ode_sol_simplicity(x, func, trysolving=not simplify)) if given_hint == 'best': return retdict['best'] for i in orderedhints: if retdict['best'] == retdict.get(i, None): retdict['best_hint'] = i break retdict['default'] = gethints['default'] retdict['order'] = gethints['order'] retdict.update(failed_hints) return retdict else: # The key 'hint' stores the hint needed to be solved for. hint = hints['hint'] return _helper_simplify(eq, hint, hints, simplify, ics=ics) def _helper_simplify(eq, hint, match, simplify=True, ics=None, **kwargs): r""" Helper function of dsolve that calls the respective :py:mod:~sympy.solvers.ode functions to solve for the ordinary differential equations. This minimizes the computation in calling :py:meth:~sympy.solvers.deutils._desolve multiple times. """ r = match if hint.endswith('_Integral'): solvefunc = globals()['ode_' + hint[:-len('_Integral')]] else: solvefunc = globals()['ode_' + hint] func = r['func'] order = r['order'] match = r[hint] free = eq.free_symbols cons = lambda s: s.free_symbols.difference(free) if simplify: # odesimp() will attempt to integrate, if necessary, apply constantsimp(), # attempt to solve for func, and apply any other hint specific # simplifications sols = solvefunc(eq, func, order, match) if isinstance(sols, Expr): rv = odesimp(sols, func, order, cons(sols), hint) else: rv = [odesimp(s, func, order, cons(s), hint) for s in sols] else: # We still want to integrate (you can disable it separately with the hint) match['simplify'] = False # Some hints can take advantage of this option rv = _handle_Integral(solvefunc(eq, func, order, match), func, order, hint) if ics and not 'power_series' in hint: if isinstance(rv, Expr): solved_constants = solve_ics([rv], [r['func']], cons(rv), ics) rv = rv.subs(solved_constants) else: rv1 = [] for s in rv: solved_constants = solve_ics([s], [r['func']], cons(s), ics) rv1.append(s.subs(solved_constants)) rv = rv1 return rv def solve_ics(sols, funcs, constants, ics): """ Solve for the constants given initial conditions sols is a list of solutions. funcs is a list of functions. constants is a list of constants. 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. Returns a dictionary mapping constants to values. solution.subs(constants) will replace the constants in solution. Example ======= >>> # From dsolve(f(x).diff(x) - f(x), f(x)) >>> from sympy import symbols, Eq, exp, Function >>> from sympy.solvers.ode import solve_ics >>> f = Function('f') >>> x, C1 = symbols('x C1') >>> sols = [Eq(f(x), C1*exp(x))] >>> funcs = [f(x)] >>> constants = [C1] >>> ics = {f(0): 2} >>> solved_constants = solve_ics(sols, funcs, constants, ics) >>> solved_constants {C1: 2} >>> sols[0].subs(solved_constants) Eq(f(x), 2*exp(x)) """ # Assume ics are of the form f(x0): value or Subs(diff(f(x), x, n), (x, # x0)): value (currently checked by classify_ode). To solve, replace x # with x0, f(x0) with value, then solve for constants. For f^(n)(x0), # differentiate the solution n times, so that f^(n)(x) appears. x = funcs[0].args[0] diff_sols = [] subs_sols = [] diff_variables = set() for funcarg, value in ics.items(): if isinstance(funcarg, AppliedUndef): x0 = funcarg.args[0] matching_func = [f for f in funcs if f.func == funcarg.func][0] S = sols elif isinstance(funcarg, (Subs, Derivative)): if isinstance(funcarg, Subs): # Make sure it stays a subs. Otherwise subs below will produce # a different looking term. funcarg = funcarg.doit() if isinstance(funcarg, Subs): deriv = funcarg.expr x0 = funcarg.point[0] variables = funcarg.expr.variables matching_func = deriv elif isinstance(funcarg, Derivative): deriv = funcarg x0 = funcarg.variables[0] variables = (x,)*len(funcarg.variables) matching_func = deriv.subs(x0, x) if variables not in diff_variables: for sol in sols: if sol.has(deriv.expr.func): diff_sols.append(Eq(sol.lhs.diff(*variables), sol.rhs.diff(*variables))) diff_variables.add(variables) S = diff_sols else: raise NotImplementedError("Unrecognized initial condition") for sol in S: if sol.has(matching_func): sol2 = sol sol2 = sol2.subs(x, x0) sol2 = sol2.subs(funcarg, value) subs_sols.append(sol2) # TODO: Use solveset here try: solved_constants = solve(subs_sols, constants, dict=True) except NotImplementedError: solved_constants = [] # XXX: We can't differentiate between the solution not existing because of # invalid initial conditions, and not existing because solve is not smart # enough. If we could use solveset, this might be improvable, but for now, # we use NotImplementedError in this case. if not solved_constants: raise NotImplementedError("Couldn't solve for initial conditions") if solved_constants == True: raise ValueError("Initial conditions did not produce any solutions for constants. Perhaps they are degenerate.") if len(solved_constants) > 1: raise NotImplementedError("Initial conditions produced too many solutions for constants") if len(solved_constants[0]) != len(constants): raise ValueError("Initial conditions did not produce a solution for all constants. Perhaps they are under-specified.") return solved_constants[0] def classify_ode(eq, func=None, dict=False, ics=None, **kwargs): r""" Returns a tuple of possible :py:meth:~sympy.solvers.ode.dsolve classifications for an ODE. The tuple is ordered so that first item is the classification that :py:meth:~sympy.solvers.ode.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 :py:meth:~sympy.solvers.ode.dsolve use a different classification, use dsolve(ODE, func, hint=). See also the :py:meth:~sympy.solvers.ode.dsolve docstring for different meta-hints you can use. If dict is true, :py:meth:~sympy.solvers.ode.classify_ode will return a dictionary of hint:match expression terms. This is intended for internal use by :py:meth:~sympy.solvers.ode.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 :py:data:~sympy.solvers.ode.allhints or the :py:mod:~sympy.solvers.ode docstring for a list of all supported hints that can be returned from :py:meth:~sympy.solvers.ode.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 :py:class:~sympy.integrals.Integral class in it. Note that a hint may do this anyway if :py:meth:~sympy.core.expr.Expr.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 :py:meth:~sympy.core.expr.Expr.integrate is an expensive routine. If :py:meth:~sympy.solvers.ode.dsolve hangs, it is probably because :py:meth:~sympy.core.expr.Expr.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 :py:meth:~sympy.solvers.ode.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 :py:class:~sympy.integrals.Integral\s 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') """ ics = sympify(ics) prep = kwargs.pop('prep', True) if func and len(func.args) != 1: raise ValueError("dsolve() and classify_ode() only " "work with functions of one variable, not %s" % func) if prep or func is None: eq, func_ = _preprocess(eq, func) if func is None: func = func_ x = func.args[0] f = func.func y = Dummy('y') xi = kwargs.get('xi') eta = kwargs.get('eta') terms = kwargs.get('n') if isinstance(eq, Equality): if eq.rhs != 0: return classify_ode(eq.lhs - eq.rhs, func, dict=dict, ics=ics, xi=xi, n=terms, eta=eta, prep=False) eq = eq.lhs order = ode_order(eq, f(x)) # hint:matchdict or hint:(tuple of matchdicts) # Also will contain "default": and "order":order items. matching_hints = {"order": order} if not order: if dict: matching_hints["default"] = None return matching_hints else: return () df = f(x).diff(x) a = Wild('a', exclude=[f(x)]) b = Wild('b', exclude=[f(x)]) c = Wild('c', exclude=[f(x)]) d = Wild('d', exclude=[df, f(x).diff(x, 2)]) e = Wild('e', exclude=[df]) k = Wild('k', exclude=[df]) n = Wild('n', exclude=[x, f(x), df]) c1 = Wild('c1', exclude=[x]) a2 = Wild('a2', exclude=[x, f(x), df]) b2 = Wild('b2', exclude=[x, f(x), df]) c2 = Wild('c2', exclude=[x, f(x), df]) d2 = Wild('d2', exclude=[x, f(x), df]) a3 = Wild('a3', exclude=[f(x), df, f(x).diff(x, 2)]) b3 = Wild('b3', exclude=[f(x), df, f(x).diff(x, 2)]) c3 = Wild('c3', exclude=[f(x), df, f(x).diff(x, 2)]) r3 = {'xi': xi, 'eta': eta} # Used for the lie_group hint boundary = {} # Used to extract initial conditions C1 = Symbol("C1") eq = expand(eq) # Preprocessing to get the initial conditions out if ics is not None: for funcarg in ics: # Separating derivatives if isinstance(funcarg, (Subs, Derivative)): # f(x).diff(x).subs(x, 0) is a Subs, but f(x).diff(x).subs(x, # y) is a Derivative if isinstance(funcarg, Subs): deriv = funcarg.expr old = funcarg.variables[0] new = funcarg.point[0] elif isinstance(funcarg, Derivative): deriv = funcarg # No information on this. Just assume it was x old = x new = funcarg.variables[0] if (isinstance(deriv, Derivative) and isinstance(deriv.args[0], AppliedUndef) and deriv.args[0].func == f and len(deriv.args[0].args) == 1 and old == x and not new.has(x) and all(i == deriv.variables[0] for i in deriv.variables) and not ics[funcarg].has(f)): dorder = ode_order(deriv, x) temp = 'f' + str(dorder) boundary.update({temp: new, temp + 'val': ics[funcarg]}) else: raise ValueError("Enter valid boundary conditions for Derivatives") # Separating functions elif isinstance(funcarg, AppliedUndef): if (funcarg.func == f and len(funcarg.args) == 1 and not funcarg.args[0].has(x) and not ics[funcarg].has(f)): boundary.update({'f0': funcarg.args[0], 'f0val': ics[funcarg]}) else: raise ValueError("Enter valid boundary conditions for Function") else: raise ValueError("Enter boundary conditions of the form ics={f(point}: value, f(x).diff(x, order).subs(x, point): value}") # Precondition to try remove f(x) from highest order derivative reduced_eq = None if eq.is_Add: deriv_coef = eq.coeff(f(x).diff(x, order)) if deriv_coef not in (1, 0): r = deriv_coef.match(a*f(x)**c1) if r and r[c1]: den = f(x)**r[c1] reduced_eq = Add(*[arg/den for arg in eq.args]) if not reduced_eq: reduced_eq = eq if order == 1: ## Linear case: a(x)*y'+b(x)*y+c(x) == 0 if eq.is_Add: ind, dep = reduced_eq.as_independent(f) else: u = Dummy('u') ind, dep = (reduced_eq + u).as_independent(f) ind, dep = [tmp.subs(u, 0) for tmp in [ind, dep]] r = {a: dep.coeff(df), b: dep.coeff(f(x)), c: ind} # double check f[a] since the preconditioning may have failed if not r[a].has(f) and not r[b].has(f) and ( r[a]*df + r[b]*f(x) + r[c]).expand() - reduced_eq == 0: r['a'] = a r['b'] = b r['c'] = c matching_hints["1st_linear"] = r matching_hints["1st_linear_Integral"] = r ## Bernoulli case: a(x)*y'+b(x)*y+c(x)*y**n == 0 r = collect( reduced_eq, f(x), exact=True).match(a*df + b*f(x) + c*f(x)**n) if r and r[c] != 0 and r[n] != 1: # See issue 4676 r['a'] = a r['b'] = b r['c'] = c r['n'] = n matching_hints["Bernoulli"] = r matching_hints["Bernoulli_Integral"] = r ## Riccati special n == -2 case: a2*y'+b2*y**2+c2*y/x+d2/x**2 == 0 r = collect(reduced_eq, f(x), exact=True).match(a2*df + b2*f(x)**2 + c2*f(x)/x + d2/x**2) if r and r[b2] != 0 and (r[c2] != 0 or r[d2] != 0): r['a2'] = a2 r['b2'] = b2 r['c2'] = c2 r['d2'] = d2 matching_hints["Riccati_special_minus2"] = r # NON-REDUCED FORM OF EQUATION matches r = collect(eq, df, exact=True).match(d + e * df) if r: r['d'] = d r['e'] = e r['y'] = y r[d] = r[d].subs(f(x), y) r[e] = r[e].subs(f(x), y) # FIRST ORDER POWER SERIES WHICH NEEDS INITIAL CONDITIONS # TODO: Hint first order series should match only if d/e is analytic. # For now, only d/e and (d/e).diff(arg) is checked for existence at # at a given point. # This is currently done internally in ode_1st_power_series. point = boundary.get('f0', 0) value = boundary.get('f0val', C1) check = cancel(r[d]/r[e]) check1 = check.subs({x: point, y: value}) if not check1.has(oo) and not check1.has(zoo) and \ not check1.has(NaN) and not check1.has(-oo): check2 = (check1.diff(x)).subs({x: point, y: value}) if not check2.has(oo) and not check2.has(zoo) and \ not check2.has(NaN) and not check2.has(-oo): rseries = r.copy() rseries.update({'terms': terms, 'f0': point, 'f0val': value}) matching_hints["1st_power_series"] = rseries r3.update(r) ## Exact Differential Equation: P(x, y) + Q(x, y)*y' = 0 where # dP/dy == dQ/dx try: if r[d] != 0: numerator = simplify(r[d].diff(y) - r[e].diff(x)) # The following few conditions try to convert a non-exact # differential equation into an exact one. # References : Differential equations with applications # and historical notes - George E. Simmons if numerator: # If (dP/dy - dQ/dx) / Q = f(x) # then exp(integral(f(x))*equation becomes exact factor = simplify(numerator/r[e]) variables = factor.free_symbols if len(variables) == 1 and x == variables.pop(): factor = exp(Integral(factor).doit()) r[d] *= factor r[e] *= factor matching_hints["1st_exact"] = r matching_hints["1st_exact_Integral"] = r else: # If (dP/dy - dQ/dx) / -P = f(y) # then exp(integral(f(y))*equation becomes exact factor = simplify(-numerator/r[d]) variables = factor.free_symbols if len(variables) == 1 and y == variables.pop(): factor = exp(Integral(factor).doit()) r[d] *= factor r[e] *= factor matching_hints["1st_exact"] = r matching_hints["1st_exact_Integral"] = r else: matching_hints["1st_exact"] = r matching_hints["1st_exact_Integral"] = r except NotImplementedError: # Differentiating the coefficients might fail because of things # like f(2*x).diff(x). See issue 4624 and issue 4719. pass # Any first order ODE can be ideally solved by the Lie Group # method matching_hints["lie_group"] = r3 # This match is used for several cases below; we now collect on # f(x) so the matching works. r = collect(reduced_eq, df, exact=True).match(d + e*df) if r: # Using r[d] and r[e] without any modification for hints # linear-coefficients and separable-reduced. num, den = r[d], r[e] # ODE = d/e + df r['d'] = d r['e'] = e r['y'] = y r[d] = num.subs(f(x), y) r[e] = den.subs(f(x), y) ## Separable Case: y' == P(y)*Q(x) r[d] = separatevars(r[d]) r[e] = separatevars(r[e]) # m1[coeff]*m1[x]*m1[y] + m2[coeff]*m2[x]*m2[y]*y' m1 = separatevars(r[d], dict=True, symbols=(x, y)) m2 = separatevars(r[e], dict=True, symbols=(x, y)) if m1 and m2: r1 = {'m1': m1, 'm2': m2, 'y': y} matching_hints["separable"] = r1 matching_hints["separable_Integral"] = r1 ## First order equation with homogeneous coefficients: # dy/dx == F(y/x) or dy/dx == F(x/y) ordera = homogeneous_order(r[d], x, y) if ordera is not None: orderb = homogeneous_order(r[e], x, y) if ordera == orderb: # u1=y/x and u2=x/y u1 = Dummy('u1') u2 = Dummy('u2') s = "1st_homogeneous_coeff_subs" s1 = s + "_dep_div_indep" s2 = s + "_indep_div_dep" if simplify((r[d] + u1*r[e]).subs({x: 1, y: u1})) != 0: matching_hints[s1] = r matching_hints[s1 + "_Integral"] = r if simplify((r[e] + u2*r[d]).subs({x: u2, y: 1})) != 0: matching_hints[s2] = r matching_hints[s2 + "_Integral"] = r if s1 in matching_hints and s2 in matching_hints: matching_hints["1st_homogeneous_coeff_best"] = r ## Linear coefficients of the form # y'+ F((a*x + b*y + c)/(a'*x + b'y + c')) = 0 # that can be reduced to homogeneous form. F = num/den params = _linear_coeff_match(F, func) if params: xarg, yarg = params u = Dummy('u') t = Dummy('t') # Dummy substitution for df and f(x). dummy_eq = reduced_eq.subs(((df, t), (f(x), u))) reps = ((x, x + xarg), (u, u + yarg), (t, df), (u, f(x))) dummy_eq = simplify(dummy_eq.subs(reps)) # get the re-cast values for e and d r2 = collect(expand(dummy_eq), [df, f(x)]).match(e*df + d) if r2: orderd = homogeneous_order(r2[d], x, f(x)) if orderd is not None: ordere = homogeneous_order(r2[e], x, f(x)) if orderd == ordere: # Match arguments are passed in such a way that it # is coherent with the already existing homogeneous # functions. r2[d] = r2[d].subs(f(x), y) r2[e] = r2[e].subs(f(x), y) r2.update({'xarg': xarg, 'yarg': yarg, 'd': d, 'e': e, 'y': y}) matching_hints["linear_coefficients"] = r2 matching_hints["linear_coefficients_Integral"] = r2 ## Equation of the form y' + (y/x)*H(x^n*y) = 0 # that can be reduced to separable form factor = simplify(x/f(x)*num/den) # Try representing factor in terms of x^n*y # where n is lowest power of x in factor; # first remove terms like sqrt(2)*3 from factor.atoms(Mul) u = None for mul in ordered(factor.atoms(Mul)): if mul.has(x): _, u = mul.as_independent(x, f(x)) break if u and u.has(f(x)): h = x**(degree(Poly(u.subs(f(x), y), gen=x)))*f(x) p = Wild('p') if (u/h == 1) or ((u/h).simplify().match(x**p)): t = Dummy('t') r2 = {'t': t} xpart, ypart = u.as_independent(f(x)) test = factor.subs(((u, t), (1/u, 1/t))) free = test.free_symbols if len(free) == 1 and free.pop() == t: r2.update({'power': xpart.as_base_exp()[1], 'u': test}) matching_hints["separable_reduced"] = r2 matching_hints["separable_reduced_Integral"] = r2 ## Almost-linear equation of the form f(x)*g(y)*y' + k(x)*l(y) + m(x) = 0 r = collect(eq, [df, f(x)]).match(e*df + d) if r: r2 = r.copy() r2[c] = S.Zero if r2[d].is_Add: # Separate the terms having f(x) to r[d] and # remaining to r[c] no_f, r2[d] = r2[d].as_independent(f(x)) r2[c] += no_f factor = simplify(r2[d].diff(f(x))/r[e]) if factor and not factor.has(f(x)): r2[d] = factor_terms(r2[d]) u = r2[d].as_independent(f(x), as_Add=False)[1] r2.update({'a': e, 'b': d, 'c': c, 'u': u}) r2[d] /= u r2[e] /= u.diff(f(x)) matching_hints["almost_linear"] = r2 matching_hints["almost_linear_Integral"] = r2 elif order == 2: # Liouville ODE in the form # f(x).diff(x, 2) + g(f(x))*(f(x).diff(x))**2 + h(x)*f(x).diff(x) # See Goldstein and Braun, "Advanced Methods for the Solution of # Differential Equations", pg. 98 s = d*f(x).diff(x, 2) + e*df**2 + k*df r = reduced_eq.match(s) if r and r[d] != 0: y = Dummy('y') g = simplify(r[e]/r[d]).subs(f(x), y) h = simplify(r[k]/r[d]) if h.has(f(x)) or g.has(x): pass else: r = {'g': g, 'h': h, 'y': y} matching_hints["Liouville"] = r matching_hints["Liouville_Integral"] = r # Homogeneous second order differential equation of the form # a3*f(x).diff(x, 2) + b3*f(x).diff(x) + c3, where # for simplicity, a3, b3 and c3 are assumed to be polynomials. # It has a definite power series solution at point x0 if, b3/a3 and c3/a3 # are analytic at x0. deq = a3*(f(x).diff(x, 2)) + b3*df + c3*f(x) r = collect(reduced_eq, [f(x).diff(x, 2), f(x).diff(x), f(x)]).match(deq) ordinary = False if r and r[a3] != 0: if all([r[key].is_polynomial() for key in r]): p = cancel(r[b3]/r[a3]) # Used below q = cancel(r[c3]/r[a3]) # Used below point = kwargs.get('x0', 0) check = p.subs(x, point) if not check.has(oo) and not check.has(NaN) and \ not check.has(zoo) and not check.has(-oo): check = q.subs(x, point) if not check.has(oo) and not check.has(NaN) and \ not check.has(zoo) and not check.has(-oo): ordinary = True r.update({'a3': a3, 'b3': b3, 'c3': c3, 'x0': point, 'terms': terms}) matching_hints["2nd_power_series_ordinary"] = r # Checking if the differential equation has a regular singular point # at x0. It has a regular singular point at x0, if (b3/a3)*(x - x0) # and (c3/a3)*((x - x0)**2) are analytic at x0. if not ordinary: p = cancel((x - point)*p) check = p.subs(x, point) if not check.has(oo) and not check.has(NaN) and \ not check.has(zoo) and not check.has(-oo): q = cancel(((x - point)**2)*q) check = q.subs(x, point) if not check.has(oo) and not check.has(NaN) and \ not check.has(zoo) and not check.has(-oo): coeff_dict = {'p': p, 'q': q, 'x0': point, 'terms': terms} matching_hints["2nd_power_series_regular"] = coeff_dict if order > 0: # nth order linear ODE # a_n(x)y^(n) + ... + a_1(x)y' + a_0(x)y = F(x) = b r = _nth_linear_match(reduced_eq, func, order) # Constant coefficient case (a_i is constant for all i) if r and not any(r[i].has(x) for i in r if i >= 0): # Inhomogeneous case: F(x) is not identically 0 if r[-1]: undetcoeff = _undetermined_coefficients_match(r[-1], x) s = "nth_linear_constant_coeff_variation_of_parameters" matching_hints[s] = r matching_hints[s + "_Integral"] = r if undetcoeff['test']: r['trialset'] = undetcoeff['trialset'] matching_hints[ "nth_linear_constant_coeff_undetermined_coefficients" ] = r # Homogeneous case: F(x) is identically 0 else: matching_hints["nth_linear_constant_coeff_homogeneous"] = r # nth order Euler equation a_n*x**n*y^(n) + ... + a_1*x*y' + a_0*y = F(x) #In case of Homogeneous euler equation F(x) = 0 def _test_term(coeff, order): r""" Linear Euler ODEs have the form K*x**order*diff(y(x),x,order) = F(x), where K is independent of x and y(x), order>= 0. So we need to check that for each term, coeff == K*x**order from some K. We have a few cases, since coeff may have several different types. """ if order < 0: raise ValueError("order should be greater than 0") if coeff == 0: return True if order == 0: if x in coeff.free_symbols: return False return True if coeff.is_Mul: if coeff.has(f(x)): return False return x**order in coeff.args elif coeff.is_Pow: return coeff.as_base_exp() == (x, order) elif order == 1: return x == coeff return False if r and not any(not _test_term(r[i], i) for i in r if i >= 0): if not r[-1]: matching_hints["nth_linear_euler_eq_homogeneous"] = r else: matching_hints["nth_linear_euler_eq_nonhomogeneous_variation_of_parameters"] = r matching_hints["nth_linear_euler_eq_nonhomogeneous_variation_of_parameters_Integral"] = r e, re = posify(r[-1].subs(x, exp(x))) undetcoeff = _undetermined_coefficients_match(e.subs(re), x) if undetcoeff['test']: r['trialset'] = undetcoeff['trialset'] matching_hints["nth_linear_euler_eq_nonhomogeneous_undetermined_coefficients"] = r # Order keys based on allhints. retlist = [i for i in allhints if i in matching_hints] if dict: # Dictionaries are ordered arbitrarily, so make note of which # hint would come first for dsolve(). Use an ordered dict in Py 3. matching_hints["default"] = retlist[0] if retlist else None matching_hints["ordered_hints"] = tuple(retlist) return matching_hints else: return tuple(retlist) def classify_sysode(eq, funcs=None, **kwargs): r""" Returns a dictionary of parameter names and values that define the system of ordinary differential equations in eq. The parameters are further used in :py:meth:~sympy.solvers.ode.dsolve for solving that system. The parameter names and values are: 'is_linear' (boolean), which tells whether the given system is linear. Note that "linear" here refers to the operator: terms such as x*diff(x,t) are nonlinear, whereas terms like sin(t)*diff(x,t) are still linear operators. 'func' (list) contains the :py:class:~sympy.core.function.Functions that appear with a derivative in the ODE, i.e. those that we are trying to solve the ODE for. 'order' (dict) with the maximum derivative for each element of the 'func' parameter. 'func_coeff' (dict) with the coefficient for each triple (equation number, function, order). The coefficients are those subexpressions that do not appear in 'func', and hence can be considered constant for purposes of ODE solving. 'eq' (list) with the equations from eq, sympified and transformed into expressions (we are solving for these expressions to be zero). 'no_of_equations' (int) is the number of equations (same as len(eq)). 'type_of_equation' (string) is an internal classification of the type of ODE. References ========== -http://eqworld.ipmnet.ru/en/solutions/sysode/sode-toc1.htm -A. D. Polyanin and A. V. Manzhirov, Handbook of Mathematics for Engineers and Scientists Examples ======== >>> from sympy import Function, Eq, symbols, diff >>> from sympy.solvers.ode import classify_sysode >>> from sympy.abc import t >>> f, x, y = symbols('f, x, y', function=True) >>> k, l, m, n = symbols('k, l, m, n', Integer=True) >>> x1 = diff(x(t), t) ; y1 = diff(y(t), t) >>> x2 = diff(x(t), t, t) ; y2 = diff(y(t), t, t) >>> eq = (Eq(5*x1, 12*x(t) - 6*y(t)), Eq(2*y1, 11*x(t) + 3*y(t))) >>> classify_sysode(eq) {'eq': [-12*x(t) + 6*y(t) + 5*Derivative(x(t), t), -11*x(t) - 3*y(t) + 2*Derivative(y(t), t)], 'func': [x(t), y(t)], 'func_coeff': {(0, x(t), 0): -12, (0, x(t), 1): 5, (0, y(t), 0): 6, (0, y(t), 1): 0, (1, x(t), 0): -11, (1, x(t), 1): 0, (1, y(t), 0): -3, (1, y(t), 1): 2}, 'is_linear': True, 'no_of_equation': 2, 'order': {x(t): 1, y(t): 1}, 'type_of_equation': 'type1'} >>> eq = (Eq(diff(x(t),t), 5*t*x(t) + t**2*y(t)), Eq(diff(y(t),t), -t**2*x(t) + 5*t*y(t))) >>> classify_sysode(eq) {'eq': [-t**2*y(t) - 5*t*x(t) + Derivative(x(t), t), t**2*x(t) - 5*t*y(t) + Derivative(y(t), t)], 'func': [x(t), y(t)], 'func_coeff': {(0, x(t), 0): -5*t, (0, x(t), 1): 1, (0, y(t), 0): -t**2, (0, y(t), 1): 0, (1, x(t), 0): t**2, (1, x(t), 1): 0, (1, y(t), 0): -5*t, (1, y(t), 1): 1}, 'is_linear': True, 'no_of_equation': 2, 'order': {x(t): 1, y(t): 1}, 'type_of_equation': 'type4'} """ # Sympify equations and convert iterables of equations into # a list of equations def _sympify(eq): return list(map(sympify, eq if iterable(eq) else [eq])) eq, funcs = (_sympify(w) for w in [eq, funcs]) for i, fi in enumerate(eq): if isinstance(fi, Equality): eq[i] = fi.lhs - fi.rhs matching_hints = {"no_of_equation":i+1} matching_hints['eq'] = eq if i==0: raise ValueError("classify_sysode() works for systems of ODEs. " "For scalar ODEs, classify_ode should be used") t = list(list(eq[0].atoms(Derivative))[0].atoms(Symbol))[0] # find all the functions if not given order = dict() if funcs==[None]: funcs = [] for eqs in eq: derivs = eqs.atoms(Derivative) func = set().union(*[d.atoms(AppliedUndef) for d in derivs]) for func_ in func: funcs.append(func_) funcs = list(set(funcs)) if len(funcs) < len(eq): raise ValueError("Number of functions given is less than number of equations %s" % funcs) func_dict = dict() for func in funcs: if not order.get(func, False): max_order = 0 for i, eqs_ in enumerate(eq): order_ = ode_order(eqs_,func) if max_order < order_: max_order = order_ eq_no = i if eq_no in func_dict: list_func = [] list_func.append(func_dict[eq_no]) list_func.append(func) func_dict[eq_no] = list_func else: func_dict[eq_no] = func order[func] = max_order funcs = [func_dict[i] for i in range(len(func_dict))] matching_hints['func'] = funcs for func in funcs: if isinstance(func, list): for func_elem in func: if len(func_elem.args) != 1: raise ValueError("dsolve() and classify_sysode() work with " "functions of one variable only, not %s" % func) else: if func and len(func.args) != 1: raise ValueError("dsolve() and classify_sysode() work with " "functions of one variable only, not %s" % func) # find the order of all equation in system of odes matching_hints["order"] = order # find coefficients of terms f(t), diff(f(t),t) and higher derivatives # and similarly for other functions g(t), diff(g(t),t) in all equations. # Here j denotes the equation number, funcs[l] denotes the function about # which we are talking about and k denotes the order of function funcs[l] # whose coefficient we are calculating. def linearity_check(eqs, j, func, is_linear_): for k in range(order[func] + 1): func_coef[j, func, k] = collect(eqs.expand(), [diff(func, t, k)]).coeff(diff(func, t, k)) if is_linear_ == True: if func_coef[j, func, k] == 0: if k == 0: coef = eqs.as_independent(func, as_Add=True)[1] for xr in range(1, ode_order(eqs,func) + 1): coef -= eqs.as_independent(diff(func, t, xr), as_Add=True)[1] if coef != 0: is_linear_ = False else: if eqs.as_independent(diff(func, t, k), as_Add=True)[1]: is_linear_ = False else: for func_ in funcs: if isinstance(func_, list): for elem_func_ in func_: dep = func_coef[j, func, k].as_independent(elem_func_, as_Add=True)[1] if dep != 0: is_linear_ = False else: dep = func_coef[j, func, k].as_independent(func_, as_Add=True)[1] if dep != 0: is_linear_ = False return is_linear_ func_coef = {} is_linear = True for j, eqs in enumerate(eq): for func in funcs: if isinstance(func, list): for func_elem in func: is_linear = linearity_check(eqs, j, func_elem, is_linear) else: is_linear = linearity_check(eqs, j, func, is_linear) matching_hints['func_coeff'] = func_coef matching_hints['is_linear'] = is_linear if len(set(order.values()))==1: order_eq = list(matching_hints['order'].values())[0] if matching_hints['is_linear'] == True: if matching_hints['no_of_equation'] == 2: if order_eq == 1: type_of_equation = check_linear_2eq_order1(eq, funcs, func_coef) elif order_eq == 2: type_of_equation = check_linear_2eq_order2(eq, funcs, func_coef) else: type_of_equation = None elif matching_hints['no_of_equation'] == 3: if order_eq == 1: type_of_equation = check_linear_3eq_order1(eq, funcs, func_coef) if type_of_equation==None: type_of_equation = check_linear_neq_order1(eq, funcs, func_coef) else: type_of_equation = None else: if order_eq == 1: type_of_equation = check_linear_neq_order1(eq, funcs, func_coef) else: type_of_equation = None else: if matching_hints['no_of_equation'] == 2: if order_eq == 1: type_of_equation = check_nonlinear_2eq_order1(eq, funcs, func_coef) else: type_of_equation = None elif matching_hints['no_of_equation'] == 3: if order_eq == 1: type_of_equation = check_nonlinear_3eq_order1(eq, funcs, func_coef) else: type_of_equation = None else: type_of_equation = None else: type_of_equation = None matching_hints['type_of_equation'] = type_of_equation return matching_hints def check_linear_2eq_order1(eq, func, func_coef): x = func[0].func y = func[1].func fc = func_coef t = list(list(eq[0].atoms(Derivative))[0].atoms(Symbol))[0] r = dict() # for equations Eq(a1*diff(x(t),t), b1*x(t) + c1*y(t) + d1) # and Eq(a2*diff(y(t),t), b2*x(t) + c2*y(t) + d2) r['a1'] = fc[0,x(t),1] ; r['a2'] = fc[1,y(t),1] r['b1'] = -fc[0,x(t),0]/fc[0,x(t),1] ; r['b2'] = -fc[1,x(t),0]/fc[1,y(t),1] r['c1'] = -fc[0,y(t),0]/fc[0,x(t),1] ; r['c2'] = -fc[1,y(t),0]/fc[1,y(t),1] forcing = [S(0),S(0)] for i in range(2): for j in Add.make_args(eq[i]): if not j.has(x(t), y(t)): forcing[i] += j if not (forcing[0].has(t) or forcing[1].has(t)): # We can handle homogeneous case and simple constant forcings r['d1'] = forcing[0] r['d2'] = forcing[1] else: # Issue #9244: nonhomogeneous linear systems are not supported return None # Conditions to check for type 6 whose equations are Eq(diff(x(t),t), f(t)*x(t) + g(t)*y(t)) and # Eq(diff(y(t),t), a*[f(t) + a*h(t)]x(t) + a*[g(t) - h(t)]*y(t)) p = 0 q = 0 p1 = cancel(r['b2']/(cancel(r['b2']/r['c2']).as_numer_denom()[0])) p2 = cancel(r['b1']/(cancel(r['b1']/r['c1']).as_numer_denom()[0])) for n, i in enumerate([p1, p2]): for j in Mul.make_args(collect_const(i)): if not j.has(t): q = j if q and n==0: if ((r['b2']/j - r['b1'])/(r['c1'] - r['c2']/j)) == j: p = 1 elif q and n==1: if ((r['b1']/j - r['b2'])/(r['c2'] - r['c1']/j)) == j: p = 2 # End of condition for type 6 if r['d1']!=0 or r['d2']!=0: if not r['d1'].has(t) and not r['d2'].has(t): if all(not r[k].has(t) for k in 'a1 a2 b1 b2 c1 c2'.split()): # Equations for type 2 are Eq(a1*diff(x(t),t),b1*x(t)+c1*y(t)+d1) and Eq(a2*diff(y(t),t),b2*x(t)+c2*y(t)+d2) return "type2" else: return None else: if all(not r[k].has(t) for k in 'a1 a2 b1 b2 c1 c2'.split()): # Equations for type 1 are Eq(a1*diff(x(t),t),b1*x(t)+c1*y(t)) and Eq(a2*diff(y(t),t),b2*x(t)+c2*y(t)) return "type1" else: r['b1'] = r['b1']/r['a1'] ; r['b2'] = r['b2']/r['a2'] r['c1'] = r['c1']/r['a1'] ; r['c2'] = r['c2']/r['a2'] if (r['b1'] == r['c2']) and (r['c1'] == r['b2']): # Equation for type 3 are Eq(diff(x(t),t), f(t)*x(t) + g(t)*y(t)) and Eq(diff(y(t),t), g(t)*x(t) + f(t)*y(t)) return "type3" elif (r['b1'] == r['c2']) and (r['c1'] == -r['b2']) or (r['b1'] == -r['c2']) and (r['c1'] == r['b2']): # Equation for type 4 are Eq(diff(x(t),t), f(t)*x(t) + g(t)*y(t)) and Eq(diff(y(t),t), -g(t)*x(t) + f(t)*y(t)) return "type4" elif (not cancel(r['b2']/r['c1']).has(t) and not cancel((r['c2']-r['b1'])/r['c1']).has(t)) \ or (not cancel(r['b1']/r['c2']).has(t) and not cancel((r['c1']-r['b2'])/r['c2']).has(t)): # Equations for type 5 are Eq(diff(x(t),t), f(t)*x(t) + g(t)*y(t)) and Eq(diff(y(t),t), a*g(t)*x(t) + [f(t) + b*g(t)]*y(t) return "type5" elif p: return "type6" else: # Equations for type 7 are Eq(diff(x(t),t), f(t)*x(t) + g(t)*y(t)) and Eq(diff(y(t),t), h(t)*x(t) + p(t)*y(t)) return "type7" def check_linear_2eq_order2(eq, func, func_coef): x = func[0].func y = func[1].func fc = func_coef t = list(list(eq[0].atoms(Derivative))[0].atoms(Symbol))[0] r = dict() a = Wild('a', exclude=[1/t]) b = Wild('b', exclude=[1/t**2]) u = Wild('u', exclude=[t, t**2]) v = Wild('v', exclude=[t, t**2]) w = Wild('w', exclude=[t, t**2]) p = Wild('p', exclude=[t, t**2]) r['a1'] = fc[0,x(t),2] ; r['a2'] = fc[1,y(t),2] r['b1'] = fc[0,x(t),1] ; r['b2'] = fc[1,x(t),1] r['c1'] = fc[0,y(t),1] ; r['c2'] = fc[1,y(t),1] r['d1'] = fc[0,x(t),0] ; r['d2'] = fc[1,x(t),0] r['e1'] = fc[0,y(t),0] ; r['e2'] = fc[1,y(t),0] const = [S(0), S(0)] for i in range(2): for j in Add.make_args(eq[i]): if not (j.has(x(t)) or j.has(y(t))): const[i] += j r['f1'] = const[0] r['f2'] = const[1] if r['f1']!=0 or r['f2']!=0: if all(not r[k].has(t) for k in 'a1 a2 d1 d2 e1 e2 f1 f2'.split()) \ and r['b1']==r['c1']==r['b2']==r['c2']==0: return "type2" elif all(not r[k].has(t) for k in 'a1 a2 b1 b2 c1 c2 d1 d2 e1 e1'.split()): p = [S(0), S(0)] ; q = [S(0), S(0)] for n, e in enumerate([r['f1'], r['f2']]): if e.has(t): tpart = e.as_independent(t, Mul)[1] for i in Mul.make_args(tpart): if i.has(exp): b, e = i.as_base_exp() co = e.coeff(t) if co and not co.has(t) and co.has(I): p[n] = 1 else: q[n] = 1 else: q[n] = 1 else: q[n] = 1 if p[0]==1 and p[1]==1 and q[0]==0 and q[1]==0: return "type4" else: return None else: return None else: if r['b1']==r['b2']==r['c1']==r['c2']==0 and all(not r[k].has(t) \ for k in 'a1 a2 d1 d2 e1 e2'.split()): return "type1" elif r['b1']==r['e1']==r['c2']==r['d2']==0 and all(not r[k].has(t) \ for k in 'a1 a2 b2 c1 d1 e2'.split()) and r['c1'] == -r['b2'] and \ r['d1'] == r['e2']: return "type3" elif cancel(-r['b2']/r['d2'])==t and cancel(-r['c1']/r['e1'])==t and not \ (r['d2']/r['a2']).has(t) and not (r['e1']/r['a1']).has(t) and \ r['b1']==r['d1']==r['c2']==r['e2']==0: return "type5" elif ((r['a1']/r['d1']).expand()).match((p*(u*t**2+v*t+w)**2).expand()) and not \ (cancel(r['a1']*r['d2']/(r['a2']*r['d1']))).has(t) and not (r['d1']/r['e1']).has(t) and not \ (r['d2']/r['e2']).has(t) and r['b1'] == r['b2'] == r['c1'] == r['c2'] == 0: return "type10" elif not cancel(r['d1']/r['e1']).has(t) and not cancel(r['d2']/r['e2']).has(t) and not \ cancel(r['d1']*r['a2']/(r['d2']*r['a1'])).has(t) and r['b1']==r['b2']==r['c1']==r['c2']==0: return "type6" elif not cancel(r['b1']/r['c1']).has(t) and not cancel(r['b2']/r['c2']).has(t) and not \ cancel(r['b1']*r['a2']/(r['b2']*r['a1'])).has(t) and r['d1']==r['d2']==r['e1']==r['e2']==0: return "type7" elif cancel(-r['b2']/r['d2'])==t and cancel(-r['c1']/r['e1'])==t and not \ cancel(r['e1']*r['a2']/(r['d2']*r['a1'])).has(t) and r['e1'].has(t) \ and r['b1']==r['d1']==r['c2']==r['e2']==0: return "type8" elif (r['b1']/r['a1']).match(a/t) and (r['b2']/r['a2']).match(a/t) and not \ (r['b1']/r['c1']).has(t) and not (r['b2']/r['c2']).has(t) and \ (r['d1']/r['a1']).match(b/t**2) and (r['d2']/r['a2']).match(b/t**2) \ and not (r['d1']/r['e1']).has(t) and not (r['d2']/r['e2']).has(t): return "type9" elif -r['b1']/r['d1']==-r['c1']/r['e1']==-r['b2']/r['d2']==-r['c2']/r['e2']==t: return "type11" else: return None def check_linear_3eq_order1(eq, func, func_coef): x = func[0].func y = func[1].func z = func[2].func fc = func_coef t = list(list(eq[0].atoms(Derivative))[0].atoms(Symbol))[0] r = dict() r['a1'] = fc[0,x(t),1]; r['a2'] = fc[1,y(t),1]; r['a3'] = fc[2,z(t),1] r['b1'] = fc[0,x(t),0]; r['b2'] = fc[1,x(t),0]; r['b3'] = fc[2,x(t),0] r['c1'] = fc[0,y(t),0]; r['c2'] = fc[1,y(t),0]; r['c3'] = fc[2,y(t),0] r['d1'] = fc[0,z(t),0]; r['d2'] = fc[1,z(t),0]; r['d3'] = fc[2,z(t),0] forcing = [S(0), S(0), S(0)] for i in range(3): for j in Add.make_args(eq[i]): if not j.has(x(t), y(t), z(t)): forcing[i] += j if forcing[0].has(t) or forcing[1].has(t) or forcing[2].has(t): # We can handle homogeneous case and simple constant forcings. # Issue #9244: nonhomogeneous linear systems are not supported return None if all(not r[k].has(t) for k in 'a1 a2 a3 b1 b2 b3 c1 c2 c3 d1 d2 d3'.split()): if r['c1']==r['d1']==r['d2']==0: return 'type1' elif r['c1'] == -r['b2'] and r['d1'] == -r['b3'] and r['d2'] == -r['c3'] \ and r['b1'] == r['c2'] == r['d3'] == 0: return 'type2' elif r['b1'] == r['c2'] == r['d3'] == 0 and r['c1']/r['a1'] == -r['d1']/r['a1'] \ and r['d2']/r['a2'] == -r['b2']/r['a2'] and r['b3']/r['a3'] == -r['c3']/r['a3']: return 'type3' else: return None else: for k1 in 'c1 d1 b2 d2 b3 c3'.split(): if r[k1] == 0: continue else: if all(not cancel(r[k1]/r[k]).has(t) for k in 'd1 b2 d2 b3 c3'.split() if r[k]!=0) \ and all(not cancel(r[k1]/(r['b1'] - r[k])).has(t) for k in 'b1 c2 d3'.split() if r['b1']!=r[k]): return 'type4' else: break return None def check_linear_neq_order1(eq, func, func_coef): x = func[0].func y = func[1].func z = func[2].func fc = func_coef t = list(list(eq[0].atoms(Derivative))[0].atoms(Symbol))[0] r = dict() n = len(eq) for i in range(n): for j in range(n): if (fc[i,func[j],0]/fc[i,func[i],1]).has(t): return None if len(eq)==3: return 'type6' return 'type1' def check_nonlinear_2eq_order1(eq, func, func_coef): t = list(list(eq[0].atoms(Derivative))[0].atoms(Symbol))[0] f = Wild('f') g = Wild('g') u, v = symbols('u, v', cls=Dummy) def check_type(x, y): r1 = eq[0].match(t*diff(x(t),t) - x(t) + f) r2 = eq[1].match(t*diff(y(t),t) - y(t) + g) if not (r1 and r2): r1 = eq[0].match(diff(x(t),t) - x(t)/t + f/t) r2 = eq[1].match(diff(y(t),t) - y(t)/t + g/t) if not (r1 and r2): r1 = (-eq[0]).match(t*diff(x(t),t) - x(t) + f) r2 = (-eq[1]).match(t*diff(y(t),t) - y(t) + g) if not (r1 and r2): r1 = (-eq[0]).match(diff(x(t),t) - x(t)/t + f/t) r2 = (-eq[1]).match(diff(y(t),t) - y(t)/t + g/t) if r1 and r2 and not (r1[f].subs(diff(x(t),t),u).subs(diff(y(t),t),v).has(t) \ or r2[g].subs(diff(x(t),t),u).subs(diff(y(t),t),v).has(t)): return 'type5' else: return None for func_ in func: if isinstance(func_, list): x = func[0][0].func y = func[0][1].func eq_type = check_type(x, y) if not eq_type: eq_type = check_type(y, x) return eq_type x = func[0].func y = func[1].func fc = func_coef n = Wild('n', exclude=[x(t),y(t)]) f1 = Wild('f1', exclude=[v,t]) f2 = Wild('f2', exclude=[v,t]) g1 = Wild('g1', exclude=[u,t]) g2 = Wild('g2', exclude=[u,t]) for i in range(2): eqs = 0 for terms in Add.make_args(eq[i]): eqs += terms/fc[i,func[i],1] eq[i] = eqs r = eq[0].match(diff(x(t),t) - x(t)**n*f) if r: g = (diff(y(t),t) - eq[1])/r[f] if r and not (g.has(x(t)) or g.subs(y(t),v).has(t) or r[f].subs(x(t),u).subs(y(t),v).has(t)): return 'type1' r = eq[0].match(diff(x(t),t) - exp(n*x(t))*f) if r: g = (diff(y(t),t) - eq[1])/r[f] if r and not (g.has(x(t)) or g.subs(y(t),v).has(t) or r[f].subs(x(t),u).subs(y(t),v).has(t)): return 'type2' g = Wild('g') r1 = eq[0].match(diff(x(t),t) - f) r2 = eq[1].match(diff(y(t),t) - g) if r1 and r2 and not (r1[f].subs(x(t),u).subs(y(t),v).has(t) or \ r2[g].subs(x(t),u).subs(y(t),v).has(t)): return 'type3' r1 = eq[0].match(diff(x(t),t) - f) r2 = eq[1].match(diff(y(t),t) - g) num, den = ( (r1[f].subs(x(t),u).subs(y(t),v))/ (r2[g].subs(x(t),u).subs(y(t),v))).as_numer_denom() R1 = num.match(f1*g1) R2 = den.match(f2*g2) phi = (r1[f].subs(x(t),u).subs(y(t),v))/num if R1 and R2: return 'type4' return None def check_nonlinear_2eq_order2(eq, func, func_coef): return None def check_nonlinear_3eq_order1(eq, func, func_coef): x = func[0].func y = func[1].func z = func[2].func fc = func_coef t = list(list(eq[0].atoms(Derivative))[0].atoms(Symbol))[0] u, v, w = symbols('u, v, w', cls=Dummy) a = Wild('a', exclude=[x(t), y(t), z(t), t]) b = Wild('b', exclude=[x(t), y(t), z(t), t]) c = Wild('c', exclude=[x(t), y(t), z(t), t]) f = Wild('f') F1 = Wild('F1') F2 = Wild('F2') F3 = Wild('F3') for i in range(3): eqs = 0 for terms in Add.make_args(eq[i]): eqs += terms/fc[i,func[i],1] eq[i] = eqs r1 = eq[0].match(diff(x(t),t) - a*y(t)*z(t)) r2 = eq[1].match(diff(y(t),t) - b*z(t)*x(t)) r3 = eq[2].match(diff(z(t),t) - c*x(t)*y(t)) if r1 and r2 and r3: num1, den1 = r1[a].as_numer_denom() num2, den2 = r2[b].as_numer_denom() num3, den3 = r3[c].as_numer_denom() if solve([num1*u-den1*(v-w), num2*v-den2*(w-u), num3*w-den3*(u-v)],[u, v]): return 'type1' r = eq[0].match(diff(x(t),t) - y(t)*z(t)*f) if r: r1 = collect_const(r[f]).match(a*f) r2 = ((diff(y(t),t) - eq[1])/r1[f]).match(b*z(t)*x(t)) r3 = ((diff(z(t),t) - eq[2])/r1[f]).match(c*x(t)*y(t)) if r1 and r2 and r3: num1, den1 = r1[a].as_numer_denom() num2, den2 = r2[b].as_numer_denom() num3, den3 = r3[c].as_numer_denom() if solve([num1*u-den1*(v-w), num2*v-den2*(w-u), num3*w-den3*(u-v)],[u, v]): return 'type2' r = eq[0].match(diff(x(t),t) - (F2-F3)) if r: r1 = collect_const(r[F2]).match(c*F2) r1.update(collect_const(r[F3]).match(b*F3)) if r1: if eq[1].has(r1[F2]) and not eq[1].has(r1[F3]): r1[F2], r1[F3] = r1[F3], r1[F2] r1[c], r1[b] = -r1[b], -r1[c] r2 = eq[1].match(diff(y(t),t) - a*r1[F3] + r1[c]*F1) if r2: r3 = (eq[2] == diff(z(t),t) - r1[b]*r2[F1] + r2[a]*r1[F2]) if r1 and r2 and r3: return 'type3' r = eq[0].match(diff(x(t),t) - z(t)*F2 + y(t)*F3) if r: r1 = collect_const(r[F2]).match(c*F2) r1.update(collect_const(r[F3]).match(b*F3)) if r1: if eq[1].has(r1[F2]) and not eq[1].has(r1[F3]): r1[F2], r1[F3] = r1[F3], r1[F2] r1[c], r1[b] = -r1[b], -r1[c] r2 = (diff(y(t),t) - eq[1]).match(a*x(t)*r1[F3] - r1[c]*z(t)*F1) if r2: r3 = (diff(z(t),t) - eq[2] == r1[b]*y(t)*r2[F1] - r2[a]*x(t)*r1[F2]) if r1 and r2 and r3: return 'type4' r = (diff(x(t),t) - eq[0]).match(x(t)*(F2 - F3)) if r: r1 = collect_const(r[F2]).match(c*F2) r1.update(collect_const(r[F3]).match(b*F3)) if r1: if eq[1].has(r1[F2]) and not eq[1].has(r1[F3]): r1[F2], r1[F3] = r1[F3], r1[F2] r1[c], r1[b] = -r1[b], -r1[c] r2 = (diff(y(t),t) - eq[1]).match(y(t)*(a*r1[F3] - r1[c]*F1)) if r2: r3 = (diff(z(t),t) - eq[2] == z(t)*(r1[b]*r2[F1] - r2[a]*r1[F2])) if r1 and r2 and r3: return 'type5' return None def check_nonlinear_3eq_order2(eq, func, func_coef): return None def checksysodesol(eqs, sols, func=None): r""" Substitutes corresponding sols for each functions into each eqs and checks that the result of substitutions for each equation is 0. The equations and solutions passed can be any iterable. This only works when each sols have one function only, like x(t) or y(t). For each function, sols can have a single solution or a list of solutions. 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. When a sequence of equations is passed, the same sequence is used to return the result for each equation with each function substituted with corresponding solutions. It tries the following method to find zero equivalence for each equation: Substitute the solutions for functions, like x(t) and y(t) into the original equations containing those functions. This function returns a tuple. The first item in the tuple is True if the substitution results for each equation is 0, and False otherwise. The second item in the tuple is what the substitution results in. Each element of the list should always be 0 corresponding to each equation if the first item is True. Note that sometimes this function may return False, but with an expression that is identically equal to 0, instead of returning True. This is because :py:meth:~sympy.simplify.simplify.simplify cannot reduce the expression to 0. If an expression returned by each function vanishes identically, then sols really is a solution to eqs. If this function seems to hang, it is probably because of a difficult simplification. Examples ======== >>> from sympy import Eq, diff, symbols, sin, cos, exp, sqrt, S >>> from sympy.solvers.ode import checksysodesol >>> C1, C2 = symbols('C1:3') >>> t = symbols('t') >>> x, y = symbols('x, y', function=True) >>> eq = (Eq(diff(x(t),t), x(t) + y(t) + 17), Eq(diff(y(t),t), -2*x(t) + y(t) + 12)) >>> sol = [Eq(x(t), (C1*sin(sqrt(2)*t) + C2*cos(sqrt(2)*t))*exp(t) - S(5)/3), ... Eq(y(t), (sqrt(2)*C1*cos(sqrt(2)*t) - sqrt(2)*C2*sin(sqrt(2)*t))*exp(t) - S(46)/3)] >>> checksysodesol(eq, sol) (True, [0, 0]) >>> eq = (Eq(diff(x(t),t),x(t)*y(t)**4), Eq(diff(y(t),t),y(t)**3)) >>> sol = [Eq(x(t), C1*exp(-1/(4*(C2 + t)))), Eq(y(t), -sqrt(2)*sqrt(-1/(C2 + t))/2), ... Eq(x(t), C1*exp(-1/(4*(C2 + t)))), Eq(y(t), sqrt(2)*sqrt(-1/(C2 + t))/2)] >>> checksysodesol(eq, sol) (True, [0, 0]) """ def _sympify(eq): return list(map(sympify, eq if iterable(eq) else [eq])) eqs = _sympify(eqs) for i in range(len(eqs)): if isinstance(eqs[i], Equality): eqs[i] = eqs[i].lhs - eqs[i].rhs if func is None: funcs = [] for eq in eqs: derivs = eq.atoms(Derivative) func = set().union(*[d.atoms(AppliedUndef) for d in derivs]) for func_ in func: funcs.append(func_) funcs = list(set(funcs)) if not all(isinstance(func, AppliedUndef) and len(func.args) == 1 for func in funcs)\ and len({func.args for func in funcs})!=1: raise ValueError("func must be a function of one variable, not %s" % func) for sol in sols: if len(sol.atoms(AppliedUndef)) != 1: raise ValueError("solutions should have one function only") if len(funcs) != len({sol.lhs for sol in sols}): raise ValueError("number of solutions provided does not match the number of equations") t = funcs[0].args[0] dictsol = dict() for sol in sols: func = list(sol.atoms(AppliedUndef))[0] if sol.rhs == func: sol = sol.reversed solved = sol.lhs == func and not sol.rhs.has(func) if not solved: rhs = solve(sol, func) if not rhs: raise NotImplementedError else: rhs = sol.rhs dictsol[func] = rhs checkeq = [] for eq in eqs: for func in funcs: eq = sub_func_doit(eq, func, dictsol[func]) ss = simplify(eq) if ss != 0: eq = ss.expand(force=True) else: eq = 0 checkeq.append(eq) if len(set(checkeq)) == 1 and list(set(checkeq))[0] == 0: return (True, checkeq) else: return (False, checkeq) @vectorize(0) def odesimp(eq, func, order, constants, hint): r""" Simplifies ODEs, including trying to solve for func and running :py:meth:~sympy.solvers.ode.constantsimp. It may use knowledge of the type of solution that the hint returns to apply additional simplifications. It also attempts to integrate any :py:class:~sympy.integrals.Integral\s in the expression, if the hint is not an _Integral hint. This function should have no effect on expressions returned by :py:meth:~sympy.solvers.ode.dsolve, as :py:meth:~sympy.solvers.ode.dsolve already calls :py:meth:~sympy.solvers.ode.odesimp, but the individual hint functions do not call :py:meth:~sympy.solvers.ode.odesimp (because the :py:meth:~sympy.solvers.ode.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' ... )) #doctest: +SKIP x --------- = C1 /f(x)\ tan|----| \2*x / """ x = func.args[0] f = func.func C1 = get_numbered_constants(eq, num=1) # First, integrate if the hint allows it. eq = _handle_Integral(eq, func, order, hint) if hint.startswith("nth_linear_euler_eq_nonhomogeneous"): eq = simplify(eq) if not isinstance(eq, Equality): raise TypeError("eq should be an instance of Equality") # Second, clean up the arbitrary constants. # Right now, nth linear hints can put as many as 2*order constants in an # expression. If that number grows with another hint, the third argument # here should be raised accordingly, or constantsimp() rewritten to handle # an arbitrary number of constants. eq = constantsimp(eq, constants) # Lastly, now that we have cleaned up the expression, try solving for func. # When CRootOf is implemented in solve(), we will want to return a CRootOf # every time instead of an Equality. # Get the f(x) on the left if possible. if eq.rhs == func and not eq.lhs.has(func): eq = [Eq(eq.rhs, eq.lhs)] # make sure we are working with lists of solutions in simplified form. if eq.lhs == func and not eq.rhs.has(func): # The solution is already solved eq = [eq] # special simplification of the rhs if hint.startswith("nth_linear_constant_coeff"): # Collect terms to make the solution look nice. # This is also necessary for constantsimp to remove unnecessary # terms from the particular solution from variation of parameters # # Collect is not behaving reliably here. The results for # some linear constant-coefficient equations with repeated # roots do not properly simplify all constants sometimes. # 'collectterms' gives different orders sometimes, and results # differ in collect based on that order. The # sort-reverse trick fixes things, but may fail in the # future. In addition, collect is splitting exponentials with # rational powers for no reason. We have to do a match # to fix this using Wilds. global collectterms try: collectterms.sort(key=default_sort_key) collectterms.reverse() except Exception: pass assert len(eq) == 1 and eq[0].lhs == f(x) sol = eq[0].rhs sol = expand_mul(sol) for i, reroot, imroot in collectterms: sol = collect(sol, x**i*exp(reroot*x)*sin(abs(imroot)*x)) sol = collect(sol, x**i*exp(reroot*x)*cos(imroot*x)) for i, reroot, imroot in collectterms: sol = collect(sol, x**i*exp(reroot*x)) del collectterms # Collect is splitting exponentials with rational powers for # no reason. We call powsimp to fix. sol = powsimp(sol) eq[0] = Eq(f(x), sol) else: # The solution is not solved, so try to solve it try: floats = any(i.is_Float for i in eq.atoms(Number)) eqsol = solve(eq, func, force=True, rational=False if floats else None) if not eqsol: raise NotImplementedError except (NotImplementedError, PolynomialError): eq = [eq] else: def _expand(expr): numer, denom = expr.as_numer_denom() if denom.is_Add: return expr else: return powsimp(expr.expand(), combine='exp', deep=True) # XXX: the rest of odesimp() expects each t to be in a # specific normal form: rational expression with numerator # expanded, but with combined exponential functions (at # least in this setup all tests pass). eq = [Eq(f(x), _expand(t)) for t in eqsol] # special simplification of the lhs. if hint.startswith("1st_homogeneous_coeff"): for j, eqi in enumerate(eq): newi = logcombine(eqi, force=True) if isinstance(newi.lhs, log) and newi.rhs == 0: newi = Eq(newi.lhs.args[0]/C1, C1) eq[j] = newi # We cleaned up the constants before solving to help the solve engine with # a simpler expression, but the solved expression could have introduced # things like -C1, so rerun constantsimp() one last time before returning. for i, eqi in enumerate(eq): eq[i] = constantsimp(eqi, constants) eq[i] = constant_renumber(eq[i], 'C', 1, 2*order) # If there is only 1 solution, return it; # otherwise return the list of solutions. if len(eq) == 1: eq = eq[0] return eq def checkodesol(ode, sol, func=None, order='auto', solve_for_func=True): r""" 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 :py:class:~sympy.core.relational.Equality that the solution satisfies, e.g. Eq(f(x), C1), Eq(f(x) + C1, 0); or simply an :py:class:~sympy.core.expr.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 :py:meth:~sympy.simplify.simplify.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) """ if not isinstance(ode, Equality): ode = Eq(ode, 0) if func is None: try: _, func = _preprocess(ode.lhs) except ValueError: funcs = [s.atoms(AppliedUndef) for s in ( sol if is_sequence(sol, set) else [sol])] funcs = set().union(*funcs) if len(funcs) != 1: raise ValueError( 'must pass func arg to checkodesol for this case.') func = funcs.pop() if not isinstance(func, AppliedUndef) or len(func.args) != 1: raise ValueError( "func must be a function of one variable, not %s" % func) if is_sequence(sol, set): return type(sol)([checkodesol(ode, i, order=order, solve_for_func=solve_for_func) for i in sol]) if not isinstance(sol, Equality): sol = Eq(func, sol) elif sol.rhs == func: sol = sol.reversed if order == 'auto': order = ode_order(ode, func) solved = sol.lhs == func and not sol.rhs.has(func) if solve_for_func and not solved: rhs = solve(sol, func) if rhs: eqs = [Eq(func, t) for t in rhs] if len(rhs) == 1: eqs = eqs[0] return checkodesol(ode, eqs, order=order, solve_for_func=False) s = True testnum = 0 x = func.args[0] while s: if testnum == 0: # First pass, try substituting a solved solution directly into the # ODE. This has the highest chance of succeeding. ode_diff = ode.lhs - ode.rhs if sol.lhs == func: s = sub_func_doit(ode_diff, func, sol.rhs) else: testnum += 1 continue ss = simplify(s) if ss: # with the new numer_denom in power.py, if we do a simple # expansion then testnum == 0 verifies all solutions. s = ss.expand(force=True) else: s = 0 testnum += 1 elif testnum == 1: # Second pass. If we cannot substitute f, try seeing if the nth # derivative is equal, this will only work for odes that are exact, # by definition. s = simplify( trigsimp(diff(sol.lhs, x, order) - diff(sol.rhs, x, order)) - trigsimp(ode.lhs) + trigsimp(ode.rhs)) # s2 = simplify( # diff(sol.lhs, x, order) - diff(sol.rhs, x, order) - \ # ode.lhs + ode.rhs) testnum += 1 elif testnum == 2: # Third pass. Try solving for df/dx and substituting that into the # ODE. Thanks to Chris Smith for suggesting this method. Many of # the comments below are his, too. # The method: # - Take each of 1..n derivatives of the solution. # - Solve each nth derivative for d^(n)f/dx^(n) # (the differential of that order) # - Back substitute into the ODE in decreasing order # (i.e., n, n-1, ...) # - Check the result for zero equivalence if sol.lhs == func and not sol.rhs.has(func): diffsols = {0: sol.rhs} elif sol.rhs == func and not sol.lhs.has(func): diffsols = {0: sol.lhs} else: diffsols = {} sol = sol.lhs - sol.rhs for i in range(1, order + 1): # Differentiation is a linear operator, so there should always # be 1 solution. Nonetheless, we test just to make sure. # We only need to solve once. After that, we automatically # have the solution to the differential in the order we want. if i == 1: ds = sol.diff(x) try: sdf = solve(ds, func.diff(x, i)) if not sdf: raise NotImplementedError except NotImplementedError: testnum += 1 break else: diffsols[i] = sdf[0] else: # This is what the solution says df/dx should be. diffsols[i] = diffsols[i - 1].diff(x) # Make sure the above didn't fail. if testnum > 2: continue else: # Substitute it into ODE to check for self consistency. lhs, rhs = ode.lhs, ode.rhs for i in range(order, -1, -1): if i == 0 and 0 not in diffsols: # We can only substitute f(x) if the solution was # solved for f(x). break lhs = sub_func_doit(lhs, func.diff(x, i), diffsols[i]) rhs = sub_func_doit(rhs, func.diff(x, i), diffsols[i]) ode_or_bool = Eq(lhs, rhs) ode_or_bool = simplify(ode_or_bool) if isinstance(ode_or_bool, (bool, BooleanAtom)): if ode_or_bool: lhs = rhs = S.Zero else: lhs = ode_or_bool.lhs rhs = ode_or_bool.rhs # No sense in overworking simplify -- just prove that the # numerator goes to zero num = trigsimp((lhs - rhs).as_numer_denom()[0]) # since solutions are obtained using force=True we test # using the same level of assumptions ## replace function with dummy so assumptions will work _func = Dummy('func') num = num.subs(func, _func) ## posify the expression num, reps = posify(num) s = simplify(num).xreplace(reps).xreplace({_func: func}) testnum += 1 else: break if not s: return (True, s) elif s is True: # The code above never was able to change s raise NotImplementedError("Unable to test if " + str(sol) + " is a solution to " + str(ode) + ".") else: return (False, s) def ode_sol_simplicity(sol, func, trysolving=True): r""" 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 :py:class:~sympy.integrals.Integral\s, 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 | len(str(sol)) | | func | | +----------------------------------------------+-------------------+ | sol contains an | oo | | :py:class:~sympy.integrals.Integral | | +----------------------------------------------+-------------------+ oo here means the SymPy infinity, which should compare greater than any integer. If you already know :py:meth:~sympy.solvers.solvers.solve cannot solve sol, you can use trysolving=False to skip that step, which is the only potentially slow step. For example, :py:meth:~sympy.solvers.ode.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) """ # TODO: if two solutions are solved for f(x), we still want to be # able to get the simpler of the two # See the docstring for the coercion rules. We check easier (faster) # things here first, to save time. if iterable(sol): # See if there are Integrals for i in sol: if ode_sol_simplicity(i, func, trysolving=trysolving) == oo: return oo return len(str(sol)) if sol.has(Integral): return oo # Next, try to solve for func. This code will change slightly when CRootOf # is implemented in solve(). Probably a CRootOf solution should fall # somewhere between a normal solution and an unsolvable expression. # First, see if they are already solved if sol.lhs == func and not sol.rhs.has(func) or \ sol.rhs == func and not sol.lhs.has(func): return -2 # We are not so lucky, try solving manually if trysolving: try: sols = solve(sol, func) if not sols: raise NotImplementedError except NotImplementedError: pass else: return -1 # Finally, a naive computation based on the length of the string version # of the expression. This may favor combined fractions because they # will not have duplicate denominators, and may slightly favor expressions # with fewer additions and subtractions, as those are separated by spaces # by the printer. # Additional ideas for simplicity heuristics are welcome, like maybe # checking if a equation has a larger domain, or if constantsimp has # introduced arbitrary constants numbered higher than the order of a # given ODE that sol is a solution of. return len(str(sol)) def _get_constant_subexpressions(expr, Cs): Cs = set(Cs) Ces = [] def _recursive_walk(expr): expr_syms = expr.free_symbols if len(expr_syms) > 0 and expr_syms.issubset(Cs): Ces.append(expr) else: if expr.func == exp: expr = expr.expand(mul=True) if expr.func in (Add, Mul): d = sift(expr.args, lambda i : i.free_symbols.issubset(Cs)) if len(d[True]) > 1: x = expr.func(*d[True]) if not x.is_number: Ces.append(x) elif isinstance(expr, Integral): if expr.free_symbols.issubset(Cs) and \ all(len(x) == 3 for x in expr.limits): Ces.append(expr) for i in expr.args: _recursive_walk(i) return _recursive_walk(expr) return Ces def __remove_linear_redundancies(expr, Cs): cnts = {i: expr.count(i) for i in Cs} Cs = [i for i in Cs if cnts[i] > 0] def _linear(expr): if isinstance(expr, Add): xs = [i for i in Cs if expr.count(i)==cnts[i] \ and 0 == expr.diff(i, 2)] d = {} for x in xs: y = expr.diff(x) if y not in d: d[y]=[] d[y].append(x) for y in d: if len(d[y]) > 1: d[y].sort(key=str) for x in d[y][1:]: expr = expr.subs(x, 0) return expr def _recursive_walk(expr): if len(expr.args) != 0: expr = expr.func(*[_recursive_walk(i) for i in expr.args]) expr = _linear(expr) return expr if isinstance(expr, Equality): lhs, rhs = [_recursive_walk(i) for i in expr.args] f = lambda i: isinstance(i, Number) or i in Cs if isinstance(lhs, Symbol) and lhs in Cs: rhs, lhs = lhs, rhs if lhs.func in (Add, Symbol) and rhs.func in (Add, Symbol): dlhs = sift([lhs] if isinstance(lhs, AtomicExpr) else lhs.args, f) drhs = sift([rhs] if isinstance(rhs, AtomicExpr) else rhs.args, f) for i in [True, False]: for hs in [dlhs, drhs]: if i not in hs: hs[i] = [0] # this calculation can be simplified lhs = Add(*dlhs[False]) - Add(*drhs[False]) rhs = Add(*drhs[True]) - Add(*dlhs[True]) elif lhs.func in (Mul, Symbol) and rhs.func in (Mul, Symbol): dlhs = sift([lhs] if isinstance(lhs, AtomicExpr) else lhs.args, f) if True in dlhs: if False not in dlhs: dlhs[False] = [1] lhs = Mul(*dlhs[False]) rhs = rhs/Mul(*dlhs[True]) return Eq(lhs, rhs) else: return _recursive_walk(expr) @vectorize(0) def constantsimp(expr, constants): r""" Simplifies an expression with arbitrary constants in it. This function is written specifically to work with :py:meth:~sympy.solvers.ode.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 :py:class:~sympy.core.add.Add\s 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 :py:class:~sympy.core.add.Add\s are expanded so e^{C_1 + x} will be simplified to C_1 e^x. Use :py:meth:~sympy.solvers.ode.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 :py:meth:~sympy.solvers.ode.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 """ # This function works recursively. The idea is that, for Mul, # Add, Pow, and Function, if the class has a constant in it, then # we can simplify it, which we do by recursing down and # simplifying up. Otherwise, we can skip that part of the # expression. Cs = constants orig_expr = expr constant_subexprs = _get_constant_subexpressions(expr, Cs) for xe in constant_subexprs: xes = list(xe.free_symbols) if not xes: continue if all([expr.count(c) == xe.count(c) for c in xes]): xes.sort(key=str) expr = expr.subs(xe, xes[0]) # try to perform common sub-expression elimination of constant terms try: commons, rexpr = cse(expr) commons.reverse() rexpr = rexpr[0] for s in commons: cs = list(s[1].atoms(Symbol)) if len(cs) == 1 and cs[0] in Cs and \ cs[0] not in rexpr.atoms(Symbol) and \ not any(cs[0] in ex for ex in commons if ex != s): rexpr = rexpr.subs(s[0], cs[0]) else: rexpr = rexpr.subs(*s) expr = rexpr except Exception: pass expr = __remove_linear_redundancies(expr, Cs) def _conditional_term_factoring(expr): new_expr = terms_gcd(expr, clear=False, deep=True, expand=False) # we do not want to factor exponentials, so handle this separately if new_expr.is_Mul: infac = False asfac = False for m in new_expr.args: if isinstance(m, exp): asfac = True elif m.is_Add: infac = any(isinstance(fi, exp) for t in m.args for fi in Mul.make_args(t)) if asfac and infac: new_expr = expr break return new_expr expr = _conditional_term_factoring(expr) # call recursively if more simplification is possible if orig_expr != expr: return constantsimp(expr, Cs) return expr def constant_renumber(expr, symbolname, startnumber, endnumber): r""" 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 :py:class:~sympy.core.symbol.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 :py:meth:~sympy.solvers.ode.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 """ if type(expr) in (set, list, tuple): return type(expr)( [constant_renumber(i, symbolname=symbolname, startnumber=startnumber, endnumber=endnumber) for i in expr] ) global newstartnumber newstartnumber = 1 constants_found = [None]*(endnumber + 2) constantsymbols = [Symbol( symbolname + "%d" % t) for t in range(startnumber, endnumber + 1)] # make a mapping to send all constantsymbols to S.One and use # that to make sure that term ordering is not dependent on # the indexed value of C C_1 = [(ci, S.One) for ci in constantsymbols] sort_key=lambda arg: default_sort_key(arg.subs(C_1)) def _constant_renumber(expr): r""" We need to have an internal recursive function so that newstartnumber maintains its values throughout recursive calls. """ global newstartnumber if isinstance(expr, Equality): return Eq( _constant_renumber(expr.lhs), _constant_renumber(expr.rhs)) if type(expr) not in (Mul, Add, Pow) and not expr.is_Function and \ not expr.has(*constantsymbols): # Base case, as above. Hope there aren't constants inside # of some other class, because they won't be renumbered. return expr elif expr.is_Piecewise: return expr elif expr in constantsymbols: if expr not in constants_found: constants_found[newstartnumber] = expr newstartnumber += 1 return expr elif expr.is_Function or expr.is_Pow or isinstance(expr, Tuple): return expr.func( *[_constant_renumber(x) for x in expr.args]) else: sortedargs = list(expr.args) sortedargs.sort(key=sort_key) return expr.func(*[_constant_renumber(x) for x in sortedargs]) expr = _constant_renumber(expr) # Renumbering happens here newconsts = symbols('C1:%d' % newstartnumber) expr = expr.subs(zip(constants_found[1:], newconsts), simultaneous=True) return expr def _handle_Integral(expr, func, order, hint): r""" Converts a solution with Integrals in it into an actual solution. For most hints, this simply runs expr.doit(). """ global y x = func.args[0] f = func.func if hint == "1st_exact": sol = (expr.doit()).subs(y, f(x)) del y elif hint == "1st_exact_Integral": sol = Eq(Subs(expr.lhs, y, f(x)), expr.rhs) del y elif hint == "nth_linear_constant_coeff_homogeneous": sol = expr elif not hint.endswith("_Integral"): sol = expr.doit() else: sol = expr return sol # FIXME: replace the general solution in the docstring with # dsolve(equation, hint='1st_exact_Integral'). You will need to be able # to have assumptions on P and Q that dP/dy = dQ/dx. def ode_1st_exact(eq, func, order, match): r""" 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 .. math:: 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. 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) References ========== - http://en.wikipedia.org/wiki/Exact_differential_equation - M. Tenenbaum & H. Pollard, "Ordinary Differential Equations", Dover 1963, pp. 73 # indirect doctest """ x = func.args[0] f = func.func r = match # d+e*diff(f(x),x) e = r[r['e']] d = r[r['d']] global y # This is the only way to pass dummy y to _handle_Integral y = r['y'] C1 = get_numbered_constants(eq, num=1) # Refer Joel Moses, "Symbolic Integration - The Stormy Decade", # Communications of the ACM, Volume 14, Number 8, August 1971, pp. 558 # which gives the method to solve an exact differential equation. sol = Integral(d, x) + Integral((e - (Integral(d, x).diff(y))), y) return Eq(sol, C1) def ode_1st_homogeneous_coeff_best(eq, func, order, match): r""" 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 :py:meth:~sympy.solvers.ode.ode_sol_simplicity. See the :py:meth:~sympy.solvers.ode.ode_1st_homogeneous_coeff_subs_indep_div_dep and :py:meth:~sympy.solvers.ode.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. 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 References ========== - http://en.wikipedia.org/wiki/Homogeneous_differential_equation - M. Tenenbaum & H. Pollard, "Ordinary Differential Equations", Dover 1963, pp. 59 # indirect doctest """ # There are two substitutions that solve the equation, u1=y/x and u2=x/y # They produce different integrals, so try them both and see which # one is easier. sol1 = ode_1st_homogeneous_coeff_subs_indep_div_dep(eq, func, order, match) sol2 = ode_1st_homogeneous_coeff_subs_dep_div_indep(eq, func, order, match) simplify = match.get('simplify', True) if simplify: # why is odesimp called here? Should it be at the usual spot? constants = sol1.free_symbols.difference(eq.free_symbols) sol1 = odesimp( sol1, func, order, constants, "1st_homogeneous_coeff_subs_indep_div_dep") constants = sol2.free_symbols.difference(eq.free_symbols) sol2 = odesimp( sol2, func, order, constants, "1st_homogeneous_coeff_subs_dep_div_indep") return min([sol1, sol2], key=lambda x: ode_sol_simplicity(x, func, trysolving=not simplify)) def ode_1st_homogeneous_coeff_subs_dep_div_indep(eq, func, order, match): r""" Solves a 1st order differential equation with homogeneous coefficients using the substitution u_1 = \frac{\text{}}{\text{}}. This is a differential equation .. math:: 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 :py:meth:~sympy.solvers.ode.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 :py:meth:~sympy.solvers.ode.ode_1st_homogeneous_coeff_best and :py:meth:~sympy.solvers.ode.ode_1st_homogeneous_coeff_subs_indep_div_dep. 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 References ========== - http://en.wikipedia.org/wiki/Homogeneous_differential_equation - M. Tenenbaum & H. Pollard, "Ordinary Differential Equations", Dover 1963, pp. 59 # indirect doctest """ x = func.args[0] f = func.func u = Dummy('u') u1 = Dummy('u1') # u1 == f(x)/x r = match # d+e*diff(f(x),x) C1 = get_numbered_constants(eq, num=1) xarg = match.get('xarg', 0) yarg = match.get('yarg', 0) int = Integral( (-r[r['e']]/(r[r['d']] + u1*r[r['e']])).subs({x: 1, r['y']: u1}), (u1, None, f(x)/x)) sol = logcombine(Eq(log(x), int + log(C1)), force=True) sol = sol.subs(f(x), u).subs(((u, u - yarg), (x, x - xarg), (u, f(x)))) return sol def ode_1st_homogeneous_coeff_subs_indep_div_dep(eq, func, order, match): r""" Solves a 1st order differential equation with homogeneous coefficients using the substitution u_2 = \frac{\text{}}{\text{}}. This is a differential equation .. math:: 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 :py:meth:~sympy.solvers.ode.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 :py:meth:~sympy.solvers.ode.ode_1st_homogeneous_coeff_best and :py:meth:~sympy.solvers.ode.ode_1st_homogeneous_coeff_subs_dep_div_indep. 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 References ========== - http://en.wikipedia.org/wiki/Homogeneous_differential_equation - M. Tenenbaum & H. Pollard, "Ordinary Differential Equations", Dover 1963, pp. 59 # indirect doctest """ x = func.args[0] f = func.func u = Dummy('u') u2 = Dummy('u2') # u2 == x/f(x) r = match # d+e*diff(f(x),x) C1 = get_numbered_constants(eq, num=1) xarg = match.get('xarg', 0) # If xarg present take xarg, else zero yarg = match.get('yarg', 0) # If yarg present take yarg, else zero int = Integral( simplify( (-r[r['d']]/(r[r['e']] + u2*r[r['d']])).subs({x: u2, r['y']: 1})), (u2, None, x/f(x))) sol = logcombine(Eq(log(f(x)), int + log(C1)), force=True) sol = sol.subs(f(x), u).subs(((u, u - yarg), (x, x - xarg), (u, f(x)))) return sol # XXX: Should this function maybe go somewhere else? def homogeneous_order(eq, *symbols): r""" 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 :py:meth:~solvers.ode.ode_1st_homogeneous_coeff_subs_dep_div_indep and :py:meth:~solvers.ode.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 """ if not symbols: raise ValueError("homogeneous_order: no symbols were given.") symset = set(symbols) eq = sympify(eq) # The following are not supported if eq.has(Order, Derivative): return None # These are all constants if (eq.is_Number or eq.is_NumberSymbol or eq.is_number ): return S.Zero # Replace all functions with dummy variables dum = numbered_symbols(prefix='d', cls=Dummy) newsyms = set() for i in [j for j in symset if getattr(j, 'is_Function')]: iargs = set(i.args) if iargs.difference(symset): return None else: dummyvar = next(dum) eq = eq.subs(i, dummyvar) symset.remove(i) newsyms.add(dummyvar) symset.update(newsyms) if not eq.free_symbols & symset: return None # assuming order of a nested function can only be equal to zero if isinstance(eq, Function): return None if homogeneous_order( eq.args[0], *tuple(symset)) != 0 else S.Zero # make the replacement of x with x*t and see if t can be factored out t = Dummy('t', positive=True) # It is sufficient that t > 0 eqs = separatevars(eq.subs([(i, t*i) for i in symset]), [t], dict=True)[t] if eqs is S.One: return S.Zero # there was no term with only t i, d = eqs.as_independent(t, as_Add=False) b, e = d.as_base_exp() if b == t: return e def ode_1st_linear(eq, func, order, match): r""" Solves 1st order linear differential equations. These are differential equations of the form .. math:: 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::