Coverage for /usr/lib/python3/dist-packages/sympy/solvers/ode/ode.py: 5%
1548 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-14 15:55 +0200
1r"""
2This module contains :py:meth:`~sympy.solvers.ode.dsolve` and different helper
3functions that it uses.
5:py:meth:`~sympy.solvers.ode.dsolve` solves ordinary differential equations.
6See the docstring on the various functions for their uses. Note that partial
7differential equations support is in ``pde.py``. Note that hint functions
8have docstrings describing their various methods, but they are intended for
9internal use. Use ``dsolve(ode, func, hint=hint)`` to solve an ODE using a
10specific hint. See also the docstring on
11:py:meth:`~sympy.solvers.ode.dsolve`.
13**Functions in this module**
15 These are the user functions in this module:
17 - :py:meth:`~sympy.solvers.ode.dsolve` - Solves ODEs.
18 - :py:meth:`~sympy.solvers.ode.classify_ode` - Classifies ODEs into
19 possible hints for :py:meth:`~sympy.solvers.ode.dsolve`.
20 - :py:meth:`~sympy.solvers.ode.checkodesol` - Checks if an equation is the
21 solution to an ODE.
22 - :py:meth:`~sympy.solvers.ode.homogeneous_order` - Returns the
23 homogeneous order of an expression.
24 - :py:meth:`~sympy.solvers.ode.infinitesimals` - Returns the infinitesimals
25 of the Lie group of point transformations of an ODE, such that it is
26 invariant.
27 - :py:meth:`~sympy.solvers.ode.checkinfsol` - Checks if the given infinitesimals
28 are the actual infinitesimals of a first order ODE.
30 These are the non-solver helper functions that are for internal use. The
31 user should use the various options to
32 :py:meth:`~sympy.solvers.ode.dsolve` to obtain the functionality provided
33 by these functions:
35 - :py:meth:`~sympy.solvers.ode.ode.odesimp` - Does all forms of ODE
36 simplification.
37 - :py:meth:`~sympy.solvers.ode.ode.ode_sol_simplicity` - A key function for
38 comparing solutions by simplicity.
39 - :py:meth:`~sympy.solvers.ode.constantsimp` - Simplifies arbitrary
40 constants.
41 - :py:meth:`~sympy.solvers.ode.ode.constant_renumber` - Renumber arbitrary
42 constants.
43 - :py:meth:`~sympy.solvers.ode.ode._handle_Integral` - Evaluate unevaluated
44 Integrals.
46 See also the docstrings of these functions.
48**Currently implemented solver methods**
50The following methods are implemented for solving ordinary differential
51equations. See the docstrings of the various hint functions for more
52information on each (run ``help(ode)``):
54 - 1st order separable differential equations.
55 - 1st order differential equations whose coefficients or `dx` and `dy` are
56 functions homogeneous of the same order.
57 - 1st order exact differential equations.
58 - 1st order linear differential equations.
59 - 1st order Bernoulli differential equations.
60 - Power series solutions for first order differential equations.
61 - Lie Group method of solving first order differential equations.
62 - 2nd order Liouville differential equations.
63 - Power series solutions for second order differential equations
64 at ordinary and regular singular points.
65 - `n`\th order differential equation that can be solved with algebraic
66 rearrangement and integration.
67 - `n`\th order linear homogeneous differential equation with constant
68 coefficients.
69 - `n`\th order linear inhomogeneous differential equation with constant
70 coefficients using the method of undetermined coefficients.
71 - `n`\th order linear inhomogeneous differential equation with constant
72 coefficients using the method of variation of parameters.
74**Philosophy behind this module**
76This module is designed to make it easy to add new ODE solving methods without
77having to mess with the solving code for other methods. The idea is that
78there is a :py:meth:`~sympy.solvers.ode.classify_ode` function, which takes in
79an ODE and tells you what hints, if any, will solve the ODE. It does this
80without attempting to solve the ODE, so it is fast. Each solving method is a
81hint, and it has its own function, named ``ode_<hint>``. That function takes
82in the ODE and any match expression gathered by
83:py:meth:`~sympy.solvers.ode.classify_ode` and returns a solved result. If
84this result has any integrals in it, the hint function will return an
85unevaluated :py:class:`~sympy.integrals.integrals.Integral` class.
86:py:meth:`~sympy.solvers.ode.dsolve`, which is the user wrapper function
87around all of this, will then call :py:meth:`~sympy.solvers.ode.ode.odesimp` on
88the result, which, among other things, will attempt to solve the equation for
89the dependent variable (the function we are solving for), simplify the
90arbitrary constants in the expression, and evaluate any integrals, if the hint
91allows it.
93**How to add new solution methods**
95If you have an ODE that you want :py:meth:`~sympy.solvers.ode.dsolve` to be
96able to solve, try to avoid adding special case code here. Instead, try
97finding a general method that will solve your ODE, as well as others. This
98way, the :py:mod:`~sympy.solvers.ode` module will become more robust, and
99unhindered by special case hacks. WolphramAlpha and Maple's
100DETools[odeadvisor] function are two resources you can use to classify a
101specific ODE. It is also better for a method to work with an `n`\th order ODE
102instead of only with specific orders, if possible.
104To add a new method, there are a few things that you need to do. First, you
105need a hint name for your method. Try to name your hint so that it is
106unambiguous with all other methods, including ones that may not be implemented
107yet. If your method uses integrals, also include a ``hint_Integral`` hint.
108If there is more than one way to solve ODEs with your method, include a hint
109for each one, as well as a ``<hint>_best`` hint. Your ``ode_<hint>_best()``
110function should choose the best using min with ``ode_sol_simplicity`` as the
111key argument. See
112:obj:`~sympy.solvers.ode.single.HomogeneousCoeffBest`, for example.
113The function that uses your method will be called ``ode_<hint>()``, so the
114hint must only use characters that are allowed in a Python function name
115(alphanumeric characters and the underscore '``_``' character). Include a
116function for every hint, except for ``_Integral`` hints
117(:py:meth:`~sympy.solvers.ode.dsolve` takes care of those automatically).
118Hint names should be all lowercase, unless a word is commonly capitalized
119(such as Integral or Bernoulli). If you have a hint that you do not want to
120run with ``all_Integral`` that does not have an ``_Integral`` counterpart (such
121as a best hint that would defeat the purpose of ``all_Integral``), you will
122need to remove it manually in the :py:meth:`~sympy.solvers.ode.dsolve` code.
123See also the :py:meth:`~sympy.solvers.ode.classify_ode` docstring for
124guidelines on writing a hint name.
126Determine *in general* how the solutions returned by your method compare with
127other methods that can potentially solve the same ODEs. Then, put your hints
128in the :py:data:`~sympy.solvers.ode.allhints` tuple in the order that they
129should be called. The ordering of this tuple determines which hints are
130default. Note that exceptions are ok, because it is easy for the user to
131choose individual hints with :py:meth:`~sympy.solvers.ode.dsolve`. In
132general, ``_Integral`` variants should go at the end of the list, and
133``_best`` variants should go before the various hints they apply to. For
134example, the ``undetermined_coefficients`` hint comes before the
135``variation_of_parameters`` hint because, even though variation of parameters
136is more general than undetermined coefficients, undetermined coefficients
137generally returns cleaner results for the ODEs that it can solve than
138variation of parameters does, and it does not require integration, so it is
139much faster.
141Next, you need to have a match expression or a function that matches the type
142of the ODE, which you should put in :py:meth:`~sympy.solvers.ode.classify_ode`
143(if the match function is more than just a few lines. It should match the
144ODE without solving for it as much as possible, so that
145:py:meth:`~sympy.solvers.ode.classify_ode` remains fast and is not hindered by
146bugs in solving code. Be sure to consider corner cases. For example, if your
147solution method involves dividing by something, make sure you exclude the case
148where that division will be 0.
150In most cases, the matching of the ODE will also give you the various parts
151that you need to solve it. You should put that in a dictionary (``.match()``
152will do this for you), and add that as ``matching_hints['hint'] = matchdict``
153in the relevant part of :py:meth:`~sympy.solvers.ode.classify_ode`.
154:py:meth:`~sympy.solvers.ode.classify_ode` will then send this to
155:py:meth:`~sympy.solvers.ode.dsolve`, which will send it to your function as
156the ``match`` argument. Your function should be named ``ode_<hint>(eq, func,
157order, match)`. If you need to send more information, put it in the ``match``
158dictionary. For example, if you had to substitute in a dummy variable in
159:py:meth:`~sympy.solvers.ode.classify_ode` to match the ODE, you will need to
160pass it to your function using the `match` dict to access it. You can access
161the independent variable using ``func.args[0]``, and the dependent variable
162(the function you are trying to solve for) as ``func.func``. If, while trying
163to solve the ODE, you find that you cannot, raise ``NotImplementedError``.
164:py:meth:`~sympy.solvers.ode.dsolve` will catch this error with the ``all``
165meta-hint, rather than causing the whole routine to fail.
167Add a docstring to your function that describes the method employed. Like
168with anything else in SymPy, you will need to add a doctest to the docstring,
169in addition to real tests in ``test_ode.py``. Try to maintain consistency
170with the other hint functions' docstrings. Add your method to the list at the
171top of this docstring. Also, add your method to ``ode.rst`` in the
172``docs/src`` directory, so that the Sphinx docs will pull its docstring into
173the main SymPy documentation. Be sure to make the Sphinx documentation by
174running ``make html`` from within the doc directory to verify that the
175docstring formats correctly.
177If your solution method involves integrating, use :py:obj:`~.Integral` instead of
178:py:meth:`~sympy.core.expr.Expr.integrate`. This allows the user to bypass
179hard/slow integration by using the ``_Integral`` variant of your hint. In
180most cases, calling :py:meth:`sympy.core.basic.Basic.doit` will integrate your
181solution. If this is not the case, you will need to write special code in
182:py:meth:`~sympy.solvers.ode.ode._handle_Integral`. Arbitrary constants should be
183symbols named ``C1``, ``C2``, and so on. All solution methods should return
184an equality instance. If you need an arbitrary number of arbitrary constants,
185you can use ``constants = numbered_symbols(prefix='C', cls=Symbol, start=1)``.
186If it is possible to solve for the dependent function in a general way, do so.
187Otherwise, do as best as you can, but do not call solve in your
188``ode_<hint>()`` function. :py:meth:`~sympy.solvers.ode.ode.odesimp` will attempt
189to solve the solution for you, so you do not need to do that. Lastly, if your
190ODE has a common simplification that can be applied to your solutions, you can
191add a special case in :py:meth:`~sympy.solvers.ode.ode.odesimp` for it. For
192example, solutions returned from the ``1st_homogeneous_coeff`` hints often
193have many :obj:`~sympy.functions.elementary.exponential.log` terms, so
194:py:meth:`~sympy.solvers.ode.ode.odesimp` calls
195:py:meth:`~sympy.simplify.simplify.logcombine` on them (it also helps to write
196the arbitrary constant as ``log(C1)`` instead of ``C1`` in this case). Also
197consider common ways that you can rearrange your solution to have
198:py:meth:`~sympy.solvers.ode.constantsimp` take better advantage of it. It is
199better to put simplification in :py:meth:`~sympy.solvers.ode.ode.odesimp` than in
200your method, because it can then be turned off with the simplify flag in
201:py:meth:`~sympy.solvers.ode.dsolve`. If you have any extraneous
202simplification in your function, be sure to only run it using ``if
203match.get('simplify', True):``, especially if it can be slow or if it can
204reduce the domain of the solution.
206Finally, as with every contribution to SymPy, your method will need to be
207tested. Add a test for each method in ``test_ode.py``. Follow the
208conventions there, i.e., test the solver using ``dsolve(eq, f(x),
209hint=your_hint)``, and also test the solution using
210:py:meth:`~sympy.solvers.ode.checkodesol` (you can put these in a separate
211tests and skip/XFAIL if it runs too slow/does not work). Be sure to call your
212hint specifically in :py:meth:`~sympy.solvers.ode.dsolve`, that way the test
213will not be broken simply by the introduction of another matching hint. If your
214method works for higher order (>1) ODEs, you will need to run ``sol =
215constant_renumber(sol, 'C', 1, order)`` for each solution, where ``order`` is
216the order of the ODE. This is because ``constant_renumber`` renumbers the
217arbitrary constants by printing order, which is platform dependent. Try to
218test every corner case of your solver, including a range of orders if it is a
219`n`\th order solver, but if your solver is slow, such as if it involves hard
220integration, try to keep the test run time down.
222Feel free to refactor existing hints to avoid duplicating code or creating
223inconsistencies. If you can show that your method exactly duplicates an
224existing method, including in the simplicity and speed of obtaining the
225solutions, then you can remove the old, less general method. The existing
226code is tested extensively in ``test_ode.py``, so if anything is broken, one
227of those tests will surely fail.
229"""
231from sympy.core import Add, S, Mul, Pow, oo
232from sympy.core.containers import Tuple
233from sympy.core.expr import AtomicExpr, Expr
234from sympy.core.function import (Function, Derivative, AppliedUndef, diff,
235 expand, expand_mul, Subs)
236from sympy.core.multidimensional import vectorize
237from sympy.core.numbers import nan, zoo, Number
238from sympy.core.relational import Equality, Eq
239from sympy.core.sorting import default_sort_key, ordered
240from sympy.core.symbol import Symbol, Wild, Dummy, symbols
241from sympy.core.sympify import sympify
242from sympy.core.traversal import preorder_traversal
244from sympy.logic.boolalg import (BooleanAtom, BooleanTrue,
245 BooleanFalse)
246from sympy.functions import exp, log, sqrt
247from sympy.functions.combinatorial.factorials import factorial
248from sympy.integrals.integrals import Integral
249from sympy.polys import (Poly, terms_gcd, PolynomialError, lcm)
250from sympy.polys.polytools import cancel
251from sympy.series import Order
252from sympy.series.series import series
253from sympy.simplify import (collect, logcombine, powsimp, # type: ignore
254 separatevars, simplify, cse)
255from sympy.simplify.radsimp import collect_const
256from sympy.solvers import checksol, solve
258from sympy.utilities import numbered_symbols
259from sympy.utilities.iterables import uniq, sift, iterable
260from sympy.solvers.deutils import _preprocess, ode_order, _desolve
263#: This is a list of hints in the order that they should be preferred by
264#: :py:meth:`~sympy.solvers.ode.classify_ode`. In general, hints earlier in the
265#: list should produce simpler solutions than those later in the list (for
266#: ODEs that fit both). For now, the order of this list is based on empirical
267#: observations by the developers of SymPy.
268#:
269#: The hint used by :py:meth:`~sympy.solvers.ode.dsolve` for a specific ODE
270#: can be overridden (see the docstring).
271#:
272#: In general, ``_Integral`` hints are grouped at the end of the list, unless
273#: there is a method that returns an unevaluable integral most of the time
274#: (which go near the end of the list anyway). ``default``, ``all``,
275#: ``best``, and ``all_Integral`` meta-hints should not be included in this
276#: list, but ``_best`` and ``_Integral`` hints should be included.
277allhints = (
278 "factorable",
279 "nth_algebraic",
280 "separable",
281 "1st_exact",
282 "1st_linear",
283 "Bernoulli",
284 "1st_rational_riccati",
285 "Riccati_special_minus2",
286 "1st_homogeneous_coeff_best",
287 "1st_homogeneous_coeff_subs_indep_div_dep",
288 "1st_homogeneous_coeff_subs_dep_div_indep",
289 "almost_linear",
290 "linear_coefficients",
291 "separable_reduced",
292 "1st_power_series",
293 "lie_group",
294 "nth_linear_constant_coeff_homogeneous",
295 "nth_linear_euler_eq_homogeneous",
296 "nth_linear_constant_coeff_undetermined_coefficients",
297 "nth_linear_euler_eq_nonhomogeneous_undetermined_coefficients",
298 "nth_linear_constant_coeff_variation_of_parameters",
299 "nth_linear_euler_eq_nonhomogeneous_variation_of_parameters",
300 "Liouville",
301 "2nd_linear_airy",
302 "2nd_linear_bessel",
303 "2nd_hypergeometric",
304 "2nd_hypergeometric_Integral",
305 "nth_order_reducible",
306 "2nd_power_series_ordinary",
307 "2nd_power_series_regular",
308 "nth_algebraic_Integral",
309 "separable_Integral",
310 "1st_exact_Integral",
311 "1st_linear_Integral",
312 "Bernoulli_Integral",
313 "1st_homogeneous_coeff_subs_indep_div_dep_Integral",
314 "1st_homogeneous_coeff_subs_dep_div_indep_Integral",
315 "almost_linear_Integral",
316 "linear_coefficients_Integral",
317 "separable_reduced_Integral",
318 "nth_linear_constant_coeff_variation_of_parameters_Integral",
319 "nth_linear_euler_eq_nonhomogeneous_variation_of_parameters_Integral",
320 "Liouville_Integral",
321 "2nd_nonlinear_autonomous_conserved",
322 "2nd_nonlinear_autonomous_conserved_Integral",
323 )
327def get_numbered_constants(eq, num=1, start=1, prefix='C'):
328 """
329 Returns a list of constants that do not occur
330 in eq already.
331 """
333 ncs = iter_numbered_constants(eq, start, prefix)
334 Cs = [next(ncs) for i in range(num)]
335 return (Cs[0] if num == 1 else tuple(Cs))
338def iter_numbered_constants(eq, start=1, prefix='C'):
339 """
340 Returns an iterator of constants that do not occur
341 in eq already.
342 """
344 if isinstance(eq, (Expr, Eq)):
345 eq = [eq]
346 elif not iterable(eq):
347 raise ValueError("Expected Expr or iterable but got %s" % eq)
349 atom_set = set().union(*[i.free_symbols for i in eq])
350 func_set = set().union(*[i.atoms(Function) for i in eq])
351 if func_set:
352 atom_set |= {Symbol(str(f.func)) for f in func_set}
353 return numbered_symbols(start=start, prefix=prefix, exclude=atom_set)
356def dsolve(eq, func=None, hint="default", simplify=True,
357 ics= None, xi=None, eta=None, x0=0, n=6, **kwargs):
358 r"""
359 Solves any (supported) kind of ordinary differential equation and
360 system of ordinary differential equations.
362 For single ordinary differential equation
363 =========================================
365 It is classified under this when number of equation in ``eq`` is one.
366 **Usage**
368 ``dsolve(eq, f(x), hint)`` -> Solve ordinary differential equation
369 ``eq`` for function ``f(x)``, using method ``hint``.
371 **Details**
373 ``eq`` can be any supported ordinary differential equation (see the
374 :py:mod:`~sympy.solvers.ode` docstring for supported methods).
375 This can either be an :py:class:`~sympy.core.relational.Equality`,
376 or an expression, which is assumed to be equal to ``0``.
378 ``f(x)`` is a function of one variable whose derivatives in that
379 variable make up the ordinary differential equation ``eq``. In
380 many cases it is not necessary to provide this; it will be
381 autodetected (and an error raised if it could not be detected).
383 ``hint`` is the solving method that you want dsolve to use. Use
384 ``classify_ode(eq, f(x))`` to get all of the possible hints for an
385 ODE. The default hint, ``default``, will use whatever hint is
386 returned first by :py:meth:`~sympy.solvers.ode.classify_ode`. See
387 Hints below for more options that you can use for hint.
389 ``simplify`` enables simplification by
390 :py:meth:`~sympy.solvers.ode.ode.odesimp`. See its docstring for more
391 information. Turn this off, for example, to disable solving of
392 solutions for ``func`` or simplification of arbitrary constants.
393 It will still integrate with this hint. Note that the solution may
394 contain more arbitrary constants than the order of the ODE with
395 this option enabled.
397 ``xi`` and ``eta`` are the infinitesimal functions of an ordinary
398 differential equation. They are the infinitesimals of the Lie group
399 of point transformations for which the differential equation is
400 invariant. The user can specify values for the infinitesimals. If
401 nothing is specified, ``xi`` and ``eta`` are calculated using
402 :py:meth:`~sympy.solvers.ode.infinitesimals` with the help of various
403 heuristics.
405 ``ics`` is the set of initial/boundary conditions for the differential equation.
406 It should be given in the form of ``{f(x0): x1, f(x).diff(x).subs(x, x2):
407 x3}`` and so on. For power series solutions, if no initial
408 conditions are specified ``f(0)`` is assumed to be ``C0`` and the power
409 series solution is calculated about 0.
411 ``x0`` is the point about which the power series solution of a differential
412 equation is to be evaluated.
414 ``n`` gives the exponent of the dependent variable up to which the power series
415 solution of a differential equation is to be evaluated.
417 **Hints**
419 Aside from the various solving methods, there are also some meta-hints
420 that you can pass to :py:meth:`~sympy.solvers.ode.dsolve`:
422 ``default``:
423 This uses whatever hint is returned first by
424 :py:meth:`~sympy.solvers.ode.classify_ode`. This is the
425 default argument to :py:meth:`~sympy.solvers.ode.dsolve`.
427 ``all``:
428 To make :py:meth:`~sympy.solvers.ode.dsolve` apply all
429 relevant classification hints, use ``dsolve(ODE, func,
430 hint="all")``. This will return a dictionary of
431 ``hint:solution`` terms. If a hint causes dsolve to raise the
432 ``NotImplementedError``, value of that hint's key will be the
433 exception object raised. The dictionary will also include
434 some special keys:
436 - ``order``: The order of the ODE. See also
437 :py:meth:`~sympy.solvers.deutils.ode_order` in
438 ``deutils.py``.
439 - ``best``: The simplest hint; what would be returned by
440 ``best`` below.
441 - ``best_hint``: The hint that would produce the solution
442 given by ``best``. If more than one hint produces the best
443 solution, the first one in the tuple returned by
444 :py:meth:`~sympy.solvers.ode.classify_ode` is chosen.
445 - ``default``: The solution that would be returned by default.
446 This is the one produced by the hint that appears first in
447 the tuple returned by
448 :py:meth:`~sympy.solvers.ode.classify_ode`.
450 ``all_Integral``:
451 This is the same as ``all``, except if a hint also has a
452 corresponding ``_Integral`` hint, it only returns the
453 ``_Integral`` hint. This is useful if ``all`` causes
454 :py:meth:`~sympy.solvers.ode.dsolve` to hang because of a
455 difficult or impossible integral. This meta-hint will also be
456 much faster than ``all``, because
457 :py:meth:`~sympy.core.expr.Expr.integrate` is an expensive
458 routine.
460 ``best``:
461 To have :py:meth:`~sympy.solvers.ode.dsolve` try all methods
462 and return the simplest one. This takes into account whether
463 the solution is solvable in the function, whether it contains
464 any Integral classes (i.e. unevaluatable integrals), and
465 which one is the shortest in size.
467 See also the :py:meth:`~sympy.solvers.ode.classify_ode` docstring for
468 more info on hints, and the :py:mod:`~sympy.solvers.ode` docstring for
469 a list of all supported hints.
471 **Tips**
473 - You can declare the derivative of an unknown function this way:
475 >>> from sympy import Function, Derivative
476 >>> from sympy.abc import x # x is the independent variable
477 >>> f = Function("f")(x) # f is a function of x
478 >>> # f_ will be the derivative of f with respect to x
479 >>> f_ = Derivative(f, x)
481 - See ``test_ode.py`` for many tests, which serves also as a set of
482 examples for how to use :py:meth:`~sympy.solvers.ode.dsolve`.
483 - :py:meth:`~sympy.solvers.ode.dsolve` always returns an
484 :py:class:`~sympy.core.relational.Equality` class (except for the
485 case when the hint is ``all`` or ``all_Integral``). If possible, it
486 solves the solution explicitly for the function being solved for.
487 Otherwise, it returns an implicit solution.
488 - Arbitrary constants are symbols named ``C1``, ``C2``, and so on.
489 - Because all solutions should be mathematically equivalent, some
490 hints may return the exact same result for an ODE. Often, though,
491 two different hints will return the same solution formatted
492 differently. The two should be equivalent. Also note that sometimes
493 the values of the arbitrary constants in two different solutions may
494 not be the same, because one constant may have "absorbed" other
495 constants into it.
496 - Do ``help(ode.ode_<hintname>)`` to get help more information on a
497 specific hint, where ``<hintname>`` is the name of a hint without
498 ``_Integral``.
500 For system of ordinary differential equations
501 =============================================
503 **Usage**
504 ``dsolve(eq, func)`` -> Solve a system of ordinary differential
505 equations ``eq`` for ``func`` being list of functions including
506 `x(t)`, `y(t)`, `z(t)` where number of functions in the list depends
507 upon the number of equations provided in ``eq``.
509 **Details**
511 ``eq`` can be any supported system of ordinary differential equations
512 This can either be an :py:class:`~sympy.core.relational.Equality`,
513 or an expression, which is assumed to be equal to ``0``.
515 ``func`` holds ``x(t)`` and ``y(t)`` being functions of one variable which
516 together with some of their derivatives make up the system of ordinary
517 differential equation ``eq``. It is not necessary to provide this; it
518 will be autodetected (and an error raised if it could not be detected).
520 **Hints**
522 The hints are formed by parameters returned by classify_sysode, combining
523 them give hints name used later for forming method name.
525 Examples
526 ========
528 >>> from sympy import Function, dsolve, Eq, Derivative, sin, cos, symbols
529 >>> from sympy.abc import x
530 >>> f = Function('f')
531 >>> dsolve(Derivative(f(x), x, x) + 9*f(x), f(x))
532 Eq(f(x), C1*sin(3*x) + C2*cos(3*x))
534 >>> eq = sin(x)*cos(f(x)) + cos(x)*sin(f(x))*f(x).diff(x)
535 >>> dsolve(eq, hint='1st_exact')
536 [Eq(f(x), -acos(C1/cos(x)) + 2*pi), Eq(f(x), acos(C1/cos(x)))]
537 >>> dsolve(eq, hint='almost_linear')
538 [Eq(f(x), -acos(C1/cos(x)) + 2*pi), Eq(f(x), acos(C1/cos(x)))]
539 >>> t = symbols('t')
540 >>> x, y = symbols('x, y', cls=Function)
541 >>> 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)))
542 >>> dsolve(eq)
543 [Eq(x(t), C1*x0(t) + C2*x0(t)*Integral(8*exp(Integral(7*t, t))*exp(Integral(12*t, t))/x0(t)**2, t)),
544 Eq(y(t), C1*y0(t) + C2*(y0(t)*Integral(8*exp(Integral(7*t, t))*exp(Integral(12*t, t))/x0(t)**2, t) +
545 exp(Integral(7*t, t))*exp(Integral(12*t, t))/x0(t)))]
546 >>> eq = (Eq(Derivative(x(t),t),x(t)*y(t)*sin(t)), Eq(Derivative(y(t),t),y(t)**2*sin(t)))
547 >>> dsolve(eq)
548 {Eq(x(t), -exp(C1)/(C2*exp(C1) - cos(t))), Eq(y(t), -1/(C1 - cos(t)))}
549 """
550 if iterable(eq):
551 from sympy.solvers.ode.systems import dsolve_system
553 # This may have to be changed in future
554 # when we have weakly and strongly
555 # connected components. This have to
556 # changed to show the systems that haven't
557 # been solved.
558 try:
559 sol = dsolve_system(eq, funcs=func, ics=ics, doit=True)
560 return sol[0] if len(sol) == 1 else sol
561 except NotImplementedError:
562 pass
564 match = classify_sysode(eq, func)
566 eq = match['eq']
567 order = match['order']
568 func = match['func']
569 t = list(list(eq[0].atoms(Derivative))[0].atoms(Symbol))[0]
571 # keep highest order term coefficient positive
572 for i in range(len(eq)):
573 for func_ in func:
574 if isinstance(func_, list):
575 pass
576 else:
577 if eq[i].coeff(diff(func[i],t,ode_order(eq[i], func[i]))).is_negative:
578 eq[i] = -eq[i]
579 match['eq'] = eq
580 if len(set(order.values()))!=1:
581 raise ValueError("It solves only those systems of equations whose orders are equal")
582 match['order'] = list(order.values())[0]
583 def recur_len(l):
584 return sum(recur_len(item) if isinstance(item,list) else 1 for item in l)
585 if recur_len(func) != len(eq):
586 raise ValueError("dsolve() and classify_sysode() work with "
587 "number of functions being equal to number of equations")
588 if match['type_of_equation'] is None:
589 raise NotImplementedError
590 else:
591 if match['is_linear'] == True:
592 solvefunc = globals()['sysode_linear_%(no_of_equation)seq_order%(order)s' % match]
593 else:
594 solvefunc = globals()['sysode_nonlinear_%(no_of_equation)seq_order%(order)s' % match]
595 sols = solvefunc(match)
596 if ics:
597 constants = Tuple(*sols).free_symbols - Tuple(*eq).free_symbols
598 solved_constants = solve_ics(sols, func, constants, ics)
599 return [sol.subs(solved_constants) for sol in sols]
600 return sols
601 else:
602 given_hint = hint # hint given by the user
604 # See the docstring of _desolve for more details.
605 hints = _desolve(eq, func=func,
606 hint=hint, simplify=True, xi=xi, eta=eta, type='ode', ics=ics,
607 x0=x0, n=n, **kwargs)
608 eq = hints.pop('eq', eq)
609 all_ = hints.pop('all', False)
610 if all_:
611 retdict = {}
612 failed_hints = {}
613 gethints = classify_ode(eq, dict=True, hint='all')
614 orderedhints = gethints['ordered_hints']
615 for hint in hints:
616 try:
617 rv = _helper_simplify(eq, hint, hints[hint], simplify)
618 except NotImplementedError as detail:
619 failed_hints[hint] = detail
620 else:
621 retdict[hint] = rv
622 func = hints[hint]['func']
624 retdict['best'] = min(list(retdict.values()), key=lambda x:
625 ode_sol_simplicity(x, func, trysolving=not simplify))
626 if given_hint == 'best':
627 return retdict['best']
628 for i in orderedhints:
629 if retdict['best'] == retdict.get(i, None):
630 retdict['best_hint'] = i
631 break
632 retdict['default'] = gethints['default']
633 retdict['order'] = gethints['order']
634 retdict.update(failed_hints)
635 return retdict
637 else:
638 # The key 'hint' stores the hint needed to be solved for.
639 hint = hints['hint']
640 return _helper_simplify(eq, hint, hints, simplify, ics=ics)
642def _helper_simplify(eq, hint, match, simplify=True, ics=None, **kwargs):
643 r"""
644 Helper function of dsolve that calls the respective
645 :py:mod:`~sympy.solvers.ode` functions to solve for the ordinary
646 differential equations. This minimizes the computation in calling
647 :py:meth:`~sympy.solvers.deutils._desolve` multiple times.
648 """
649 r = match
650 func = r['func']
651 order = r['order']
652 match = r[hint]
654 if isinstance(match, SingleODESolver):
655 solvefunc = match
656 elif hint.endswith('_Integral'):
657 solvefunc = globals()['ode_' + hint[:-len('_Integral')]]
658 else:
659 solvefunc = globals()['ode_' + hint]
661 free = eq.free_symbols
662 cons = lambda s: s.free_symbols.difference(free)
664 if simplify:
665 # odesimp() will attempt to integrate, if necessary, apply constantsimp(),
666 # attempt to solve for func, and apply any other hint specific
667 # simplifications
668 if isinstance(solvefunc, SingleODESolver):
669 sols = solvefunc.get_general_solution()
670 else:
671 sols = solvefunc(eq, func, order, match)
672 if iterable(sols):
673 rv = [odesimp(eq, s, func, hint) for s in sols]
674 else:
675 rv = odesimp(eq, sols, func, hint)
676 else:
677 # We still want to integrate (you can disable it separately with the hint)
678 if isinstance(solvefunc, SingleODESolver):
679 exprs = solvefunc.get_general_solution(simplify=False)
680 else:
681 match['simplify'] = False # Some hints can take advantage of this option
682 exprs = solvefunc(eq, func, order, match)
683 if isinstance(exprs, list):
684 rv = [_handle_Integral(expr, func, hint) for expr in exprs]
685 else:
686 rv = _handle_Integral(exprs, func, hint)
688 if isinstance(rv, list):
689 if simplify:
690 rv = _remove_redundant_solutions(eq, rv, order, func.args[0])
691 if len(rv) == 1:
692 rv = rv[0]
693 if ics and 'power_series' not in hint:
694 if isinstance(rv, (Expr, Eq)):
695 solved_constants = solve_ics([rv], [r['func']], cons(rv), ics)
696 rv = rv.subs(solved_constants)
697 else:
698 rv1 = []
699 for s in rv:
700 try:
701 solved_constants = solve_ics([s], [r['func']], cons(s), ics)
702 except ValueError:
703 continue
704 rv1.append(s.subs(solved_constants))
705 if len(rv1) == 1:
706 return rv1[0]
707 rv = rv1
708 return rv
710def solve_ics(sols, funcs, constants, ics):
711 """
712 Solve for the constants given initial conditions
714 ``sols`` is a list of solutions.
716 ``funcs`` is a list of functions.
718 ``constants`` is a list of constants.
720 ``ics`` is the set of initial/boundary conditions for the differential
721 equation. It should be given in the form of ``{f(x0): x1,
722 f(x).diff(x).subs(x, x2): x3}`` and so on.
724 Returns a dictionary mapping constants to values.
725 ``solution.subs(constants)`` will replace the constants in ``solution``.
727 Example
728 =======
729 >>> # From dsolve(f(x).diff(x) - f(x), f(x))
730 >>> from sympy import symbols, Eq, exp, Function
731 >>> from sympy.solvers.ode.ode import solve_ics
732 >>> f = Function('f')
733 >>> x, C1 = symbols('x C1')
734 >>> sols = [Eq(f(x), C1*exp(x))]
735 >>> funcs = [f(x)]
736 >>> constants = [C1]
737 >>> ics = {f(0): 2}
738 >>> solved_constants = solve_ics(sols, funcs, constants, ics)
739 >>> solved_constants
740 {C1: 2}
741 >>> sols[0].subs(solved_constants)
742 Eq(f(x), 2*exp(x))
744 """
745 # Assume ics are of the form f(x0): value or Subs(diff(f(x), x, n), (x,
746 # x0)): value (currently checked by classify_ode). To solve, replace x
747 # with x0, f(x0) with value, then solve for constants. For f^(n)(x0),
748 # differentiate the solution n times, so that f^(n)(x) appears.
749 x = funcs[0].args[0]
750 diff_sols = []
751 subs_sols = []
752 diff_variables = set()
753 for funcarg, value in ics.items():
754 if isinstance(funcarg, AppliedUndef):
755 x0 = funcarg.args[0]
756 matching_func = [f for f in funcs if f.func == funcarg.func][0]
757 S = sols
758 elif isinstance(funcarg, (Subs, Derivative)):
759 if isinstance(funcarg, Subs):
760 # Make sure it stays a subs. Otherwise subs below will produce
761 # a different looking term.
762 funcarg = funcarg.doit()
763 if isinstance(funcarg, Subs):
764 deriv = funcarg.expr
765 x0 = funcarg.point[0]
766 variables = funcarg.expr.variables
767 matching_func = deriv
768 elif isinstance(funcarg, Derivative):
769 deriv = funcarg
770 x0 = funcarg.variables[0]
771 variables = (x,)*len(funcarg.variables)
772 matching_func = deriv.subs(x0, x)
773 for sol in sols:
774 if sol.has(deriv.expr.func):
775 diff_sols.append(Eq(sol.lhs.diff(*variables), sol.rhs.diff(*variables)))
776 diff_variables.add(variables)
777 S = diff_sols
778 else:
779 raise NotImplementedError("Unrecognized initial condition")
781 for sol in S:
782 if sol.has(matching_func):
783 sol2 = sol
784 sol2 = sol2.subs(x, x0)
785 sol2 = sol2.subs(funcarg, value)
786 # This check is necessary because of issue #15724
787 if not isinstance(sol2, BooleanAtom) or not subs_sols:
788 subs_sols = [s for s in subs_sols if not isinstance(s, BooleanAtom)]
789 subs_sols.append(sol2)
791 # TODO: Use solveset here
792 try:
793 solved_constants = solve(subs_sols, constants, dict=True)
794 except NotImplementedError:
795 solved_constants = []
797 # XXX: We can't differentiate between the solution not existing because of
798 # invalid initial conditions, and not existing because solve is not smart
799 # enough. If we could use solveset, this might be improvable, but for now,
800 # we use NotImplementedError in this case.
801 if not solved_constants:
802 raise ValueError("Couldn't solve for initial conditions")
804 if solved_constants == True:
805 raise ValueError("Initial conditions did not produce any solutions for constants. Perhaps they are degenerate.")
807 if len(solved_constants) > 1:
808 raise NotImplementedError("Initial conditions produced too many solutions for constants")
810 return solved_constants[0]
812def classify_ode(eq, func=None, dict=False, ics=None, *, prep=True, xi=None, eta=None, n=None, **kwargs):
813 r"""
814 Returns a tuple of possible :py:meth:`~sympy.solvers.ode.dsolve`
815 classifications for an ODE.
817 The tuple is ordered so that first item is the classification that
818 :py:meth:`~sympy.solvers.ode.dsolve` uses to solve the ODE by default. In
819 general, classifications at the near the beginning of the list will
820 produce better solutions faster than those near the end, thought there are
821 always exceptions. To make :py:meth:`~sympy.solvers.ode.dsolve` use a
822 different classification, use ``dsolve(ODE, func,
823 hint=<classification>)``. See also the
824 :py:meth:`~sympy.solvers.ode.dsolve` docstring for different meta-hints
825 you can use.
827 If ``dict`` is true, :py:meth:`~sympy.solvers.ode.classify_ode` will
828 return a dictionary of ``hint:match`` expression terms. This is intended
829 for internal use by :py:meth:`~sympy.solvers.ode.dsolve`. Note that
830 because dictionaries are ordered arbitrarily, this will most likely not be
831 in the same order as the tuple.
833 You can get help on different hints by executing
834 ``help(ode.ode_hintname)``, where ``hintname`` is the name of the hint
835 without ``_Integral``.
837 See :py:data:`~sympy.solvers.ode.allhints` or the
838 :py:mod:`~sympy.solvers.ode` docstring for a list of all supported hints
839 that can be returned from :py:meth:`~sympy.solvers.ode.classify_ode`.
841 Notes
842 =====
844 These are remarks on hint names.
846 ``_Integral``
848 If a classification has ``_Integral`` at the end, it will return the
849 expression with an unevaluated :py:class:`~.Integral`
850 class in it. Note that a hint may do this anyway if
851 :py:meth:`~sympy.core.expr.Expr.integrate` cannot do the integral,
852 though just using an ``_Integral`` will do so much faster. Indeed, an
853 ``_Integral`` hint will always be faster than its corresponding hint
854 without ``_Integral`` because
855 :py:meth:`~sympy.core.expr.Expr.integrate` is an expensive routine.
856 If :py:meth:`~sympy.solvers.ode.dsolve` hangs, it is probably because
857 :py:meth:`~sympy.core.expr.Expr.integrate` is hanging on a tough or
858 impossible integral. Try using an ``_Integral`` hint or
859 ``all_Integral`` to get it return something.
861 Note that some hints do not have ``_Integral`` counterparts. This is
862 because :py:func:`~sympy.integrals.integrals.integrate` is not used in
863 solving the ODE for those method. For example, `n`\th order linear
864 homogeneous ODEs with constant coefficients do not require integration
865 to solve, so there is no
866 ``nth_linear_homogeneous_constant_coeff_Integrate`` hint. You can
867 easily evaluate any unevaluated
868 :py:class:`~sympy.integrals.integrals.Integral`\s in an expression by
869 doing ``expr.doit()``.
871 Ordinals
873 Some hints contain an ordinal such as ``1st_linear``. This is to help
874 differentiate them from other hints, as well as from other methods
875 that may not be implemented yet. If a hint has ``nth`` in it, such as
876 the ``nth_linear`` hints, this means that the method used to applies
877 to ODEs of any order.
879 ``indep`` and ``dep``
881 Some hints contain the words ``indep`` or ``dep``. These reference
882 the independent variable and the dependent function, respectively. For
883 example, if an ODE is in terms of `f(x)`, then ``indep`` will refer to
884 `x` and ``dep`` will refer to `f`.
886 ``subs``
888 If a hints has the word ``subs`` in it, it means that the ODE is solved
889 by substituting the expression given after the word ``subs`` for a
890 single dummy variable. This is usually in terms of ``indep`` and
891 ``dep`` as above. The substituted expression will be written only in
892 characters allowed for names of Python objects, meaning operators will
893 be spelled out. For example, ``indep``/``dep`` will be written as
894 ``indep_div_dep``.
896 ``coeff``
898 The word ``coeff`` in a hint refers to the coefficients of something
899 in the ODE, usually of the derivative terms. See the docstring for
900 the individual methods for more info (``help(ode)``). This is
901 contrast to ``coefficients``, as in ``undetermined_coefficients``,
902 which refers to the common name of a method.
904 ``_best``
906 Methods that have more than one fundamental way to solve will have a
907 hint for each sub-method and a ``_best`` meta-classification. This
908 will evaluate all hints and return the best, using the same
909 considerations as the normal ``best`` meta-hint.
912 Examples
913 ========
915 >>> from sympy import Function, classify_ode, Eq
916 >>> from sympy.abc import x
917 >>> f = Function('f')
918 >>> classify_ode(Eq(f(x).diff(x), 0), f(x))
919 ('nth_algebraic',
920 'separable',
921 '1st_exact',
922 '1st_linear',
923 'Bernoulli',
924 '1st_homogeneous_coeff_best',
925 '1st_homogeneous_coeff_subs_indep_div_dep',
926 '1st_homogeneous_coeff_subs_dep_div_indep',
927 '1st_power_series', 'lie_group', 'nth_linear_constant_coeff_homogeneous',
928 'nth_linear_euler_eq_homogeneous',
929 'nth_algebraic_Integral', 'separable_Integral', '1st_exact_Integral',
930 '1st_linear_Integral', 'Bernoulli_Integral',
931 '1st_homogeneous_coeff_subs_indep_div_dep_Integral',
932 '1st_homogeneous_coeff_subs_dep_div_indep_Integral')
933 >>> classify_ode(f(x).diff(x, 2) + 3*f(x).diff(x) + 2*f(x) - 4)
934 ('factorable', 'nth_linear_constant_coeff_undetermined_coefficients',
935 'nth_linear_constant_coeff_variation_of_parameters',
936 'nth_linear_constant_coeff_variation_of_parameters_Integral')
938 """
939 ics = sympify(ics)
941 if func and len(func.args) != 1:
942 raise ValueError("dsolve() and classify_ode() only "
943 "work with functions of one variable, not %s" % func)
945 if isinstance(eq, Equality):
946 eq = eq.lhs - eq.rhs
948 # Some methods want the unprocessed equation
949 eq_orig = eq
951 if prep or func is None:
952 eq, func_ = _preprocess(eq, func)
953 if func is None:
954 func = func_
955 x = func.args[0]
956 f = func.func
957 y = Dummy('y')
958 terms = 5 if n is None else n
960 order = ode_order(eq, f(x))
961 # hint:matchdict or hint:(tuple of matchdicts)
962 # Also will contain "default":<default hint> and "order":order items.
963 matching_hints = {"order": order}
965 df = f(x).diff(x)
966 a = Wild('a', exclude=[f(x)])
967 d = Wild('d', exclude=[df, f(x).diff(x, 2)])
968 e = Wild('e', exclude=[df])
969 n = Wild('n', exclude=[x, f(x), df])
970 c1 = Wild('c1', exclude=[x])
971 a3 = Wild('a3', exclude=[f(x), df, f(x).diff(x, 2)])
972 b3 = Wild('b3', exclude=[f(x), df, f(x).diff(x, 2)])
973 c3 = Wild('c3', exclude=[f(x), df, f(x).diff(x, 2)])
974 boundary = {} # Used to extract initial conditions
975 C1 = Symbol("C1")
977 # Preprocessing to get the initial conditions out
978 if ics is not None:
979 for funcarg in ics:
980 # Separating derivatives
981 if isinstance(funcarg, (Subs, Derivative)):
982 # f(x).diff(x).subs(x, 0) is a Subs, but f(x).diff(x).subs(x,
983 # y) is a Derivative
984 if isinstance(funcarg, Subs):
985 deriv = funcarg.expr
986 old = funcarg.variables[0]
987 new = funcarg.point[0]
988 elif isinstance(funcarg, Derivative):
989 deriv = funcarg
990 # No information on this. Just assume it was x
991 old = x
992 new = funcarg.variables[0]
994 if (isinstance(deriv, Derivative) and isinstance(deriv.args[0],
995 AppliedUndef) and deriv.args[0].func == f and
996 len(deriv.args[0].args) == 1 and old == x and not
997 new.has(x) and all(i == deriv.variables[0] for i in
998 deriv.variables) and x not in ics[funcarg].free_symbols):
1000 dorder = ode_order(deriv, x)
1001 temp = 'f' + str(dorder)
1002 boundary.update({temp: new, temp + 'val': ics[funcarg]})
1003 else:
1004 raise ValueError("Invalid boundary conditions for Derivatives")
1007 # Separating functions
1008 elif isinstance(funcarg, AppliedUndef):
1009 if (funcarg.func == f and len(funcarg.args) == 1 and
1010 not funcarg.args[0].has(x) and x not in ics[funcarg].free_symbols):
1011 boundary.update({'f0': funcarg.args[0], 'f0val': ics[funcarg]})
1012 else:
1013 raise ValueError("Invalid boundary conditions for Function")
1015 else:
1016 raise ValueError("Enter boundary conditions of the form ics={f(point): value, f(x).diff(x, order).subs(x, point): value}")
1018 ode = SingleODEProblem(eq_orig, func, x, prep=prep, xi=xi, eta=eta)
1019 user_hint = kwargs.get('hint', 'default')
1020 # Used when dsolve is called without an explicit hint.
1021 # We exit early to return the first valid match
1022 early_exit = (user_hint=='default')
1023 if user_hint.endswith('_Integral'):
1024 user_hint = user_hint[:-len('_Integral')]
1025 user_map = solver_map
1026 # An explicit hint has been given to dsolve
1027 # Skip matching code for other hints
1028 if user_hint not in ['default', 'all', 'all_Integral', 'best'] and user_hint in solver_map:
1029 user_map = {user_hint: solver_map[user_hint]}
1031 for hint in user_map:
1032 solver = user_map[hint](ode)
1033 if solver.matches():
1034 matching_hints[hint] = solver
1035 if user_map[hint].has_integral:
1036 matching_hints[hint + "_Integral"] = solver
1037 if dict and early_exit:
1038 matching_hints["default"] = hint
1039 return matching_hints
1041 eq = expand(eq)
1042 # Precondition to try remove f(x) from highest order derivative
1043 reduced_eq = None
1044 if eq.is_Add:
1045 deriv_coef = eq.coeff(f(x).diff(x, order))
1046 if deriv_coef not in (1, 0):
1047 r = deriv_coef.match(a*f(x)**c1)
1048 if r and r[c1]:
1049 den = f(x)**r[c1]
1050 reduced_eq = Add(*[arg/den for arg in eq.args])
1051 if not reduced_eq:
1052 reduced_eq = eq
1054 if order == 1:
1056 # NON-REDUCED FORM OF EQUATION matches
1057 r = collect(eq, df, exact=True).match(d + e * df)
1058 if r:
1059 r['d'] = d
1060 r['e'] = e
1061 r['y'] = y
1062 r[d] = r[d].subs(f(x), y)
1063 r[e] = r[e].subs(f(x), y)
1065 # FIRST ORDER POWER SERIES WHICH NEEDS INITIAL CONDITIONS
1066 # TODO: Hint first order series should match only if d/e is analytic.
1067 # For now, only d/e and (d/e).diff(arg) is checked for existence at
1068 # at a given point.
1069 # This is currently done internally in ode_1st_power_series.
1070 point = boundary.get('f0', 0)
1071 value = boundary.get('f0val', C1)
1072 check = cancel(r[d]/r[e])
1073 check1 = check.subs({x: point, y: value})
1074 if not check1.has(oo) and not check1.has(zoo) and \
1075 not check1.has(nan) and not check1.has(-oo):
1076 check2 = (check1.diff(x)).subs({x: point, y: value})
1077 if not check2.has(oo) and not check2.has(zoo) and \
1078 not check2.has(nan) and not check2.has(-oo):
1079 rseries = r.copy()
1080 rseries.update({'terms': terms, 'f0': point, 'f0val': value})
1081 matching_hints["1st_power_series"] = rseries
1083 elif order == 2:
1084 # Homogeneous second order differential equation of the form
1085 # a3*f(x).diff(x, 2) + b3*f(x).diff(x) + c3
1086 # It has a definite power series solution at point x0 if, b3/a3 and c3/a3
1087 # are analytic at x0.
1088 deq = a3*(f(x).diff(x, 2)) + b3*df + c3*f(x)
1089 r = collect(reduced_eq,
1090 [f(x).diff(x, 2), f(x).diff(x), f(x)]).match(deq)
1091 ordinary = False
1092 if r:
1093 if not all(r[key].is_polynomial() for key in r):
1094 n, d = reduced_eq.as_numer_denom()
1095 reduced_eq = expand(n)
1096 r = collect(reduced_eq,
1097 [f(x).diff(x, 2), f(x).diff(x), f(x)]).match(deq)
1098 if r and r[a3] != 0:
1099 p = cancel(r[b3]/r[a3]) # Used below
1100 q = cancel(r[c3]/r[a3]) # Used below
1101 point = kwargs.get('x0', 0)
1102 check = p.subs(x, point)
1103 if not check.has(oo, nan, zoo, -oo):
1104 check = q.subs(x, point)
1105 if not check.has(oo, nan, zoo, -oo):
1106 ordinary = True
1107 r.update({'a3': a3, 'b3': b3, 'c3': c3, 'x0': point, 'terms': terms})
1108 matching_hints["2nd_power_series_ordinary"] = r
1110 # Checking if the differential equation has a regular singular point
1111 # at x0. It has a regular singular point at x0, if (b3/a3)*(x - x0)
1112 # and (c3/a3)*((x - x0)**2) are analytic at x0.
1113 if not ordinary:
1114 p = cancel((x - point)*p)
1115 check = p.subs(x, point)
1116 if not check.has(oo, nan, zoo, -oo):
1117 q = cancel(((x - point)**2)*q)
1118 check = q.subs(x, point)
1119 if not check.has(oo, nan, zoo, -oo):
1120 coeff_dict = {'p': p, 'q': q, 'x0': point, 'terms': terms}
1121 matching_hints["2nd_power_series_regular"] = coeff_dict
1124 # Order keys based on allhints.
1125 retlist = [i for i in allhints if i in matching_hints]
1126 if dict:
1127 # Dictionaries are ordered arbitrarily, so make note of which
1128 # hint would come first for dsolve(). Use an ordered dict in Py 3.
1129 matching_hints["default"] = retlist[0] if retlist else None
1130 matching_hints["ordered_hints"] = tuple(retlist)
1131 return matching_hints
1132 else:
1133 return tuple(retlist)
1136def classify_sysode(eq, funcs=None, **kwargs):
1137 r"""
1138 Returns a dictionary of parameter names and values that define the system
1139 of ordinary differential equations in ``eq``.
1140 The parameters are further used in
1141 :py:meth:`~sympy.solvers.ode.dsolve` for solving that system.
1143 Some parameter names and values are:
1145 'is_linear' (boolean), which tells whether the given system is linear.
1146 Note that "linear" here refers to the operator: terms such as ``x*diff(x,t)`` are
1147 nonlinear, whereas terms like ``sin(t)*diff(x,t)`` are still linear operators.
1149 'func' (list) contains the :py:class:`~sympy.core.function.Function`s that
1150 appear with a derivative in the ODE, i.e. those that we are trying to solve
1151 the ODE for.
1153 'order' (dict) with the maximum derivative for each element of the 'func'
1154 parameter.
1156 'func_coeff' (dict or Matrix) with the coefficient for each triple ``(equation number,
1157 function, order)```. The coefficients are those subexpressions that do not
1158 appear in 'func', and hence can be considered constant for purposes of ODE
1159 solving. The value of this parameter can also be a Matrix if the system of ODEs are
1160 linear first order of the form X' = AX where X is the vector of dependent variables.
1161 Here, this function returns the coefficient matrix A.
1163 'eq' (list) with the equations from ``eq``, sympified and transformed into
1164 expressions (we are solving for these expressions to be zero).
1166 'no_of_equations' (int) is the number of equations (same as ``len(eq)``).
1168 'type_of_equation' (string) is an internal classification of the type of
1169 ODE.
1171 'is_constant' (boolean), which tells if the system of ODEs is constant coefficient
1172 or not. This key is temporary addition for now and is in the match dict only when
1173 the system of ODEs is linear first order constant coefficient homogeneous. So, this
1174 key's value is True for now if it is available else it does not exist.
1176 'is_homogeneous' (boolean), which tells if the system of ODEs is homogeneous. Like the
1177 key 'is_constant', this key is a temporary addition and it is True since this key value
1178 is available only when the system is linear first order constant coefficient homogeneous.
1180 References
1181 ==========
1182 -https://eqworld.ipmnet.ru/en/solutions/sysode/sode-toc1.htm
1183 -A. D. Polyanin and A. V. Manzhirov, Handbook of Mathematics for Engineers and Scientists
1185 Examples
1186 ========
1188 >>> from sympy import Function, Eq, symbols, diff
1189 >>> from sympy.solvers.ode.ode import classify_sysode
1190 >>> from sympy.abc import t
1191 >>> f, x, y = symbols('f, x, y', cls=Function)
1192 >>> k, l, m, n = symbols('k, l, m, n', Integer=True)
1193 >>> x1 = diff(x(t), t) ; y1 = diff(y(t), t)
1194 >>> x2 = diff(x(t), t, t) ; y2 = diff(y(t), t, t)
1195 >>> eq = (Eq(x1, 12*x(t) - 6*y(t)), Eq(y1, 11*x(t) + 3*y(t)))
1196 >>> classify_sysode(eq)
1197 {'eq': [-12*x(t) + 6*y(t) + Derivative(x(t), t), -11*x(t) - 3*y(t) + Derivative(y(t), t)], 'func': [x(t), y(t)],
1198 'func_coeff': {(0, x(t), 0): -12, (0, x(t), 1): 1, (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): 1}, 'is_linear': True, 'no_of_equation': 2, 'order': {x(t): 1, y(t): 1}, 'type_of_equation': None}
1199 >>> eq = (Eq(diff(x(t),t), 5*t*x(t) + t**2*y(t) + 2), Eq(diff(y(t),t), -t**2*x(t) + 5*t*y(t)))
1200 >>> classify_sysode(eq)
1201 {'eq': [-t**2*y(t) - 5*t*x(t) + Derivative(x(t), t) - 2, t**2*x(t) - 5*t*y(t) + Derivative(y(t), t)],
1202 '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,
1203 (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,
1204 'order': {x(t): 1, y(t): 1}, 'type_of_equation': None}
1206 """
1208 # Sympify equations and convert iterables of equations into
1209 # a list of equations
1210 def _sympify(eq):
1211 return list(map(sympify, eq if iterable(eq) else [eq]))
1213 eq, funcs = (_sympify(w) for w in [eq, funcs])
1214 for i, fi in enumerate(eq):
1215 if isinstance(fi, Equality):
1216 eq[i] = fi.lhs - fi.rhs
1218 t = list(list(eq[0].atoms(Derivative))[0].atoms(Symbol))[0]
1219 matching_hints = {"no_of_equation":i+1}
1220 matching_hints['eq'] = eq
1221 if i==0:
1222 raise ValueError("classify_sysode() works for systems of ODEs. "
1223 "For scalar ODEs, classify_ode should be used")
1225 # find all the functions if not given
1226 order = {}
1227 if funcs==[None]:
1228 funcs = _extract_funcs(eq)
1230 funcs = list(set(funcs))
1231 if len(funcs) != len(eq):
1232 raise ValueError("Number of functions given is not equal to the number of equations %s" % funcs)
1234 # This logic of list of lists in funcs to
1235 # be replaced later.
1236 func_dict = {}
1237 for func in funcs:
1238 if not order.get(func, False):
1239 max_order = 0
1240 for i, eqs_ in enumerate(eq):
1241 order_ = ode_order(eqs_,func)
1242 if max_order < order_:
1243 max_order = order_
1244 eq_no = i
1245 if eq_no in func_dict:
1246 func_dict[eq_no] = [func_dict[eq_no], func]
1247 else:
1248 func_dict[eq_no] = func
1249 order[func] = max_order
1251 funcs = [func_dict[i] for i in range(len(func_dict))]
1252 matching_hints['func'] = funcs
1253 for func in funcs:
1254 if isinstance(func, list):
1255 for func_elem in func:
1256 if len(func_elem.args) != 1:
1257 raise ValueError("dsolve() and classify_sysode() work with "
1258 "functions of one variable only, not %s" % func)
1259 else:
1260 if func and len(func.args) != 1:
1261 raise ValueError("dsolve() and classify_sysode() work with "
1262 "functions of one variable only, not %s" % func)
1264 # find the order of all equation in system of odes
1265 matching_hints["order"] = order
1267 # find coefficients of terms f(t), diff(f(t),t) and higher derivatives
1268 # and similarly for other functions g(t), diff(g(t),t) in all equations.
1269 # Here j denotes the equation number, funcs[l] denotes the function about
1270 # which we are talking about and k denotes the order of function funcs[l]
1271 # whose coefficient we are calculating.
1272 def linearity_check(eqs, j, func, is_linear_):
1273 for k in range(order[func] + 1):
1274 func_coef[j, func, k] = collect(eqs.expand(), [diff(func, t, k)]).coeff(diff(func, t, k))
1275 if is_linear_ == True:
1276 if func_coef[j, func, k] == 0:
1277 if k == 0:
1278 coef = eqs.as_independent(func, as_Add=True)[1]
1279 for xr in range(1, ode_order(eqs,func) + 1):
1280 coef -= eqs.as_independent(diff(func, t, xr), as_Add=True)[1]
1281 if coef != 0:
1282 is_linear_ = False
1283 else:
1284 if eqs.as_independent(diff(func, t, k), as_Add=True)[1]:
1285 is_linear_ = False
1286 else:
1287 for func_ in funcs:
1288 if isinstance(func_, list):
1289 for elem_func_ in func_:
1290 dep = func_coef[j, func, k].as_independent(elem_func_, as_Add=True)[1]
1291 if dep != 0:
1292 is_linear_ = False
1293 else:
1294 dep = func_coef[j, func, k].as_independent(func_, as_Add=True)[1]
1295 if dep != 0:
1296 is_linear_ = False
1297 return is_linear_
1299 func_coef = {}
1300 is_linear = True
1301 for j, eqs in enumerate(eq):
1302 for func in funcs:
1303 if isinstance(func, list):
1304 for func_elem in func:
1305 is_linear = linearity_check(eqs, j, func_elem, is_linear)
1306 else:
1307 is_linear = linearity_check(eqs, j, func, is_linear)
1308 matching_hints['func_coeff'] = func_coef
1309 matching_hints['is_linear'] = is_linear
1312 if len(set(order.values())) == 1:
1313 order_eq = list(matching_hints['order'].values())[0]
1314 if matching_hints['is_linear'] == True:
1315 if matching_hints['no_of_equation'] == 2:
1316 if order_eq == 1:
1317 type_of_equation = check_linear_2eq_order1(eq, funcs, func_coef)
1318 else:
1319 type_of_equation = None
1320 # If the equation does not match up with any of the
1321 # general case solvers in systems.py and the number
1322 # of equations is greater than 2, then NotImplementedError
1323 # should be raised.
1324 else:
1325 type_of_equation = None
1327 else:
1328 if matching_hints['no_of_equation'] == 2:
1329 if order_eq == 1:
1330 type_of_equation = check_nonlinear_2eq_order1(eq, funcs, func_coef)
1331 else:
1332 type_of_equation = None
1333 elif matching_hints['no_of_equation'] == 3:
1334 if order_eq == 1:
1335 type_of_equation = check_nonlinear_3eq_order1(eq, funcs, func_coef)
1336 else:
1337 type_of_equation = None
1338 else:
1339 type_of_equation = None
1340 else:
1341 type_of_equation = None
1343 matching_hints['type_of_equation'] = type_of_equation
1345 return matching_hints
1348def check_linear_2eq_order1(eq, func, func_coef):
1349 x = func[0].func
1350 y = func[1].func
1351 fc = func_coef
1352 t = list(list(eq[0].atoms(Derivative))[0].atoms(Symbol))[0]
1353 r = {}
1354 # for equations Eq(a1*diff(x(t),t), b1*x(t) + c1*y(t) + d1)
1355 # and Eq(a2*diff(y(t),t), b2*x(t) + c2*y(t) + d2)
1356 r['a1'] = fc[0,x(t),1] ; r['a2'] = fc[1,y(t),1]
1357 r['b1'] = -fc[0,x(t),0]/fc[0,x(t),1] ; r['b2'] = -fc[1,x(t),0]/fc[1,y(t),1]
1358 r['c1'] = -fc[0,y(t),0]/fc[0,x(t),1] ; r['c2'] = -fc[1,y(t),0]/fc[1,y(t),1]
1359 forcing = [S.Zero,S.Zero]
1360 for i in range(2):
1361 for j in Add.make_args(eq[i]):
1362 if not j.has(x(t), y(t)):
1363 forcing[i] += j
1364 if not (forcing[0].has(t) or forcing[1].has(t)):
1365 # We can handle homogeneous case and simple constant forcings
1366 r['d1'] = forcing[0]
1367 r['d2'] = forcing[1]
1368 else:
1369 # Issue #9244: nonhomogeneous linear systems are not supported
1370 return None
1372 # Conditions to check for type 6 whose equations are Eq(diff(x(t),t), f(t)*x(t) + g(t)*y(t)) and
1373 # Eq(diff(y(t),t), a*[f(t) + a*h(t)]x(t) + a*[g(t) - h(t)]*y(t))
1374 p = 0
1375 q = 0
1376 p1 = cancel(r['b2']/(cancel(r['b2']/r['c2']).as_numer_denom()[0]))
1377 p2 = cancel(r['b1']/(cancel(r['b1']/r['c1']).as_numer_denom()[0]))
1378 for n, i in enumerate([p1, p2]):
1379 for j in Mul.make_args(collect_const(i)):
1380 if not j.has(t):
1381 q = j
1382 if q and n==0:
1383 if ((r['b2']/j - r['b1'])/(r['c1'] - r['c2']/j)) == j:
1384 p = 1
1385 elif q and n==1:
1386 if ((r['b1']/j - r['b2'])/(r['c2'] - r['c1']/j)) == j:
1387 p = 2
1388 # End of condition for type 6
1390 if r['d1']!=0 or r['d2']!=0:
1391 return None
1392 else:
1393 if not any(r[k].has(t) for k in 'a1 a2 b1 b2 c1 c2'.split()):
1394 return None
1395 else:
1396 r['b1'] = r['b1']/r['a1'] ; r['b2'] = r['b2']/r['a2']
1397 r['c1'] = r['c1']/r['a1'] ; r['c2'] = r['c2']/r['a2']
1398 if p:
1399 return "type6"
1400 else:
1401 # 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))
1402 return "type7"
1403def check_nonlinear_2eq_order1(eq, func, func_coef):
1404 t = list(list(eq[0].atoms(Derivative))[0].atoms(Symbol))[0]
1405 f = Wild('f')
1406 g = Wild('g')
1407 u, v = symbols('u, v', cls=Dummy)
1408 def check_type(x, y):
1409 r1 = eq[0].match(t*diff(x(t),t) - x(t) + f)
1410 r2 = eq[1].match(t*diff(y(t),t) - y(t) + g)
1411 if not (r1 and r2):
1412 r1 = eq[0].match(diff(x(t),t) - x(t)/t + f/t)
1413 r2 = eq[1].match(diff(y(t),t) - y(t)/t + g/t)
1414 if not (r1 and r2):
1415 r1 = (-eq[0]).match(t*diff(x(t),t) - x(t) + f)
1416 r2 = (-eq[1]).match(t*diff(y(t),t) - y(t) + g)
1417 if not (r1 and r2):
1418 r1 = (-eq[0]).match(diff(x(t),t) - x(t)/t + f/t)
1419 r2 = (-eq[1]).match(diff(y(t),t) - y(t)/t + g/t)
1420 if r1 and r2 and not (r1[f].subs(diff(x(t),t),u).subs(diff(y(t),t),v).has(t) \
1421 or r2[g].subs(diff(x(t),t),u).subs(diff(y(t),t),v).has(t)):
1422 return 'type5'
1423 else:
1424 return None
1425 for func_ in func:
1426 if isinstance(func_, list):
1427 x = func[0][0].func
1428 y = func[0][1].func
1429 eq_type = check_type(x, y)
1430 if not eq_type:
1431 eq_type = check_type(y, x)
1432 return eq_type
1433 x = func[0].func
1434 y = func[1].func
1435 fc = func_coef
1436 n = Wild('n', exclude=[x(t),y(t)])
1437 f1 = Wild('f1', exclude=[v,t])
1438 f2 = Wild('f2', exclude=[v,t])
1439 g1 = Wild('g1', exclude=[u,t])
1440 g2 = Wild('g2', exclude=[u,t])
1441 for i in range(2):
1442 eqs = 0
1443 for terms in Add.make_args(eq[i]):
1444 eqs += terms/fc[i,func[i],1]
1445 eq[i] = eqs
1446 r = eq[0].match(diff(x(t),t) - x(t)**n*f)
1447 if r:
1448 g = (diff(y(t),t) - eq[1])/r[f]
1449 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)):
1450 return 'type1'
1451 r = eq[0].match(diff(x(t),t) - exp(n*x(t))*f)
1452 if r:
1453 g = (diff(y(t),t) - eq[1])/r[f]
1454 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)):
1455 return 'type2'
1456 g = Wild('g')
1457 r1 = eq[0].match(diff(x(t),t) - f)
1458 r2 = eq[1].match(diff(y(t),t) - g)
1459 if r1 and r2 and not (r1[f].subs(x(t),u).subs(y(t),v).has(t) or \
1460 r2[g].subs(x(t),u).subs(y(t),v).has(t)):
1461 return 'type3'
1462 r1 = eq[0].match(diff(x(t),t) - f)
1463 r2 = eq[1].match(diff(y(t),t) - g)
1464 num, den = (
1465 (r1[f].subs(x(t),u).subs(y(t),v))/
1466 (r2[g].subs(x(t),u).subs(y(t),v))).as_numer_denom()
1467 R1 = num.match(f1*g1)
1468 R2 = den.match(f2*g2)
1469 # phi = (r1[f].subs(x(t),u).subs(y(t),v))/num
1470 if R1 and R2:
1471 return 'type4'
1472 return None
1475def check_nonlinear_2eq_order2(eq, func, func_coef):
1476 return None
1478def check_nonlinear_3eq_order1(eq, func, func_coef):
1479 x = func[0].func
1480 y = func[1].func
1481 z = func[2].func
1482 fc = func_coef
1483 t = list(list(eq[0].atoms(Derivative))[0].atoms(Symbol))[0]
1484 u, v, w = symbols('u, v, w', cls=Dummy)
1485 a = Wild('a', exclude=[x(t), y(t), z(t), t])
1486 b = Wild('b', exclude=[x(t), y(t), z(t), t])
1487 c = Wild('c', exclude=[x(t), y(t), z(t), t])
1488 f = Wild('f')
1489 F1 = Wild('F1')
1490 F2 = Wild('F2')
1491 F3 = Wild('F3')
1492 for i in range(3):
1493 eqs = 0
1494 for terms in Add.make_args(eq[i]):
1495 eqs += terms/fc[i,func[i],1]
1496 eq[i] = eqs
1497 r1 = eq[0].match(diff(x(t),t) - a*y(t)*z(t))
1498 r2 = eq[1].match(diff(y(t),t) - b*z(t)*x(t))
1499 r3 = eq[2].match(diff(z(t),t) - c*x(t)*y(t))
1500 if r1 and r2 and r3:
1501 num1, den1 = r1[a].as_numer_denom()
1502 num2, den2 = r2[b].as_numer_denom()
1503 num3, den3 = r3[c].as_numer_denom()
1504 if solve([num1*u-den1*(v-w), num2*v-den2*(w-u), num3*w-den3*(u-v)],[u, v]):
1505 return 'type1'
1506 r = eq[0].match(diff(x(t),t) - y(t)*z(t)*f)
1507 if r:
1508 r1 = collect_const(r[f]).match(a*f)
1509 r2 = ((diff(y(t),t) - eq[1])/r1[f]).match(b*z(t)*x(t))
1510 r3 = ((diff(z(t),t) - eq[2])/r1[f]).match(c*x(t)*y(t))
1511 if r1 and r2 and r3:
1512 num1, den1 = r1[a].as_numer_denom()
1513 num2, den2 = r2[b].as_numer_denom()
1514 num3, den3 = r3[c].as_numer_denom()
1515 if solve([num1*u-den1*(v-w), num2*v-den2*(w-u), num3*w-den3*(u-v)],[u, v]):
1516 return 'type2'
1517 r = eq[0].match(diff(x(t),t) - (F2-F3))
1518 if r:
1519 r1 = collect_const(r[F2]).match(c*F2)
1520 r1.update(collect_const(r[F3]).match(b*F3))
1521 if r1:
1522 if eq[1].has(r1[F2]) and not eq[1].has(r1[F3]):
1523 r1[F2], r1[F3] = r1[F3], r1[F2]
1524 r1[c], r1[b] = -r1[b], -r1[c]
1525 r2 = eq[1].match(diff(y(t),t) - a*r1[F3] + r1[c]*F1)
1526 if r2:
1527 r3 = (eq[2] == diff(z(t),t) - r1[b]*r2[F1] + r2[a]*r1[F2])
1528 if r1 and r2 and r3:
1529 return 'type3'
1530 r = eq[0].match(diff(x(t),t) - z(t)*F2 + y(t)*F3)
1531 if r:
1532 r1 = collect_const(r[F2]).match(c*F2)
1533 r1.update(collect_const(r[F3]).match(b*F3))
1534 if r1:
1535 if eq[1].has(r1[F2]) and not eq[1].has(r1[F3]):
1536 r1[F2], r1[F3] = r1[F3], r1[F2]
1537 r1[c], r1[b] = -r1[b], -r1[c]
1538 r2 = (diff(y(t),t) - eq[1]).match(a*x(t)*r1[F3] - r1[c]*z(t)*F1)
1539 if r2:
1540 r3 = (diff(z(t),t) - eq[2] == r1[b]*y(t)*r2[F1] - r2[a]*x(t)*r1[F2])
1541 if r1 and r2 and r3:
1542 return 'type4'
1543 r = (diff(x(t),t) - eq[0]).match(x(t)*(F2 - F3))
1544 if r:
1545 r1 = collect_const(r[F2]).match(c*F2)
1546 r1.update(collect_const(r[F3]).match(b*F3))
1547 if r1:
1548 if eq[1].has(r1[F2]) and not eq[1].has(r1[F3]):
1549 r1[F2], r1[F3] = r1[F3], r1[F2]
1550 r1[c], r1[b] = -r1[b], -r1[c]
1551 r2 = (diff(y(t),t) - eq[1]).match(y(t)*(a*r1[F3] - r1[c]*F1))
1552 if r2:
1553 r3 = (diff(z(t),t) - eq[2] == z(t)*(r1[b]*r2[F1] - r2[a]*r1[F2]))
1554 if r1 and r2 and r3:
1555 return 'type5'
1556 return None
1559def check_nonlinear_3eq_order2(eq, func, func_coef):
1560 return None
1563@vectorize(0)
1564def odesimp(ode, eq, func, hint):
1565 r"""
1566 Simplifies solutions of ODEs, including trying to solve for ``func`` and
1567 running :py:meth:`~sympy.solvers.ode.constantsimp`.
1569 It may use knowledge of the type of solution that the hint returns to
1570 apply additional simplifications.
1572 It also attempts to integrate any :py:class:`~sympy.integrals.integrals.Integral`\s
1573 in the expression, if the hint is not an ``_Integral`` hint.
1575 This function should have no effect on expressions returned by
1576 :py:meth:`~sympy.solvers.ode.dsolve`, as
1577 :py:meth:`~sympy.solvers.ode.dsolve` already calls
1578 :py:meth:`~sympy.solvers.ode.ode.odesimp`, but the individual hint functions
1579 do not call :py:meth:`~sympy.solvers.ode.ode.odesimp` (because the
1580 :py:meth:`~sympy.solvers.ode.dsolve` wrapper does). Therefore, this
1581 function is designed for mainly internal use.
1583 Examples
1584 ========
1586 >>> from sympy import sin, symbols, dsolve, pprint, Function
1587 >>> from sympy.solvers.ode.ode import odesimp
1588 >>> x, u2, C1= symbols('x,u2,C1')
1589 >>> f = Function('f')
1591 >>> eq = dsolve(x*f(x).diff(x) - f(x) - x*sin(f(x)/x), f(x),
1592 ... hint='1st_homogeneous_coeff_subs_indep_div_dep_Integral',
1593 ... simplify=False)
1594 >>> pprint(eq, wrap_line=False)
1595 x
1596 ----
1597 f(x)
1598 /
1599 |
1600 | / 1 \
1601 | -|u1 + -------|
1602 | | /1 \|
1603 | | sin|--||
1604 | \ \u1//
1605 log(f(x)) = log(C1) + | ---------------- d(u1)
1606 | 2
1607 | u1
1608 |
1609 /
1611 >>> pprint(odesimp(eq, f(x), 1, {C1},
1612 ... hint='1st_homogeneous_coeff_subs_indep_div_dep'
1613 ... )) #doctest: +SKIP
1614 x
1615 --------- = C1
1616 /f(x)\
1617 tan|----|
1618 \2*x /
1620 """
1621 x = func.args[0]
1622 f = func.func
1623 C1 = get_numbered_constants(eq, num=1)
1624 constants = eq.free_symbols - ode.free_symbols
1626 # First, integrate if the hint allows it.
1627 eq = _handle_Integral(eq, func, hint)
1628 if hint.startswith("nth_linear_euler_eq_nonhomogeneous"):
1629 eq = simplify(eq)
1630 if not isinstance(eq, Equality):
1631 raise TypeError("eq should be an instance of Equality")
1633 # Second, clean up the arbitrary constants.
1634 # Right now, nth linear hints can put as many as 2*order constants in an
1635 # expression. If that number grows with another hint, the third argument
1636 # here should be raised accordingly, or constantsimp() rewritten to handle
1637 # an arbitrary number of constants.
1638 eq = constantsimp(eq, constants)
1640 # Lastly, now that we have cleaned up the expression, try solving for func.
1641 # When CRootOf is implemented in solve(), we will want to return a CRootOf
1642 # every time instead of an Equality.
1644 # Get the f(x) on the left if possible.
1645 if eq.rhs == func and not eq.lhs.has(func):
1646 eq = [Eq(eq.rhs, eq.lhs)]
1648 # make sure we are working with lists of solutions in simplified form.
1649 if eq.lhs == func and not eq.rhs.has(func):
1650 # The solution is already solved
1651 eq = [eq]
1653 else:
1654 # The solution is not solved, so try to solve it
1655 try:
1656 floats = any(i.is_Float for i in eq.atoms(Number))
1657 eqsol = solve(eq, func, force=True, rational=False if floats else None)
1658 if not eqsol:
1659 raise NotImplementedError
1660 except (NotImplementedError, PolynomialError):
1661 eq = [eq]
1662 else:
1663 def _expand(expr):
1664 numer, denom = expr.as_numer_denom()
1666 if denom.is_Add:
1667 return expr
1668 else:
1669 return powsimp(expr.expand(), combine='exp', deep=True)
1671 # XXX: the rest of odesimp() expects each ``t`` to be in a
1672 # specific normal form: rational expression with numerator
1673 # expanded, but with combined exponential functions (at
1674 # least in this setup all tests pass).
1675 eq = [Eq(f(x), _expand(t)) for t in eqsol]
1677 # special simplification of the lhs.
1678 if hint.startswith("1st_homogeneous_coeff"):
1679 for j, eqi in enumerate(eq):
1680 newi = logcombine(eqi, force=True)
1681 if isinstance(newi.lhs, log) and newi.rhs == 0:
1682 newi = Eq(newi.lhs.args[0]/C1, C1)
1683 eq[j] = newi
1685 # We cleaned up the constants before solving to help the solve engine with
1686 # a simpler expression, but the solved expression could have introduced
1687 # things like -C1, so rerun constantsimp() one last time before returning.
1688 for i, eqi in enumerate(eq):
1689 eq[i] = constantsimp(eqi, constants)
1690 eq[i] = constant_renumber(eq[i], ode.free_symbols)
1692 # If there is only 1 solution, return it;
1693 # otherwise return the list of solutions.
1694 if len(eq) == 1:
1695 eq = eq[0]
1696 return eq
1699def ode_sol_simplicity(sol, func, trysolving=True):
1700 r"""
1701 Returns an extended integer representing how simple a solution to an ODE
1702 is.
1704 The following things are considered, in order from most simple to least:
1706 - ``sol`` is solved for ``func``.
1707 - ``sol`` is not solved for ``func``, but can be if passed to solve (e.g.,
1708 a solution returned by ``dsolve(ode, func, simplify=False``).
1709 - If ``sol`` is not solved for ``func``, then base the result on the
1710 length of ``sol``, as computed by ``len(str(sol))``.
1711 - If ``sol`` has any unevaluated :py:class:`~sympy.integrals.integrals.Integral`\s,
1712 this will automatically be considered less simple than any of the above.
1714 This function returns an integer such that if solution A is simpler than
1715 solution B by above metric, then ``ode_sol_simplicity(sola, func) <
1716 ode_sol_simplicity(solb, func)``.
1718 Currently, the following are the numbers returned, but if the heuristic is
1719 ever improved, this may change. Only the ordering is guaranteed.
1721 +----------------------------------------------+-------------------+
1722 | Simplicity | Return |
1723 +==============================================+===================+
1724 | ``sol`` solved for ``func`` | ``-2`` |
1725 +----------------------------------------------+-------------------+
1726 | ``sol`` not solved for ``func`` but can be | ``-1`` |
1727 +----------------------------------------------+-------------------+
1728 | ``sol`` is not solved nor solvable for | ``len(str(sol))`` |
1729 | ``func`` | |
1730 +----------------------------------------------+-------------------+
1731 | ``sol`` contains an | ``oo`` |
1732 | :obj:`~sympy.integrals.integrals.Integral` | |
1733 +----------------------------------------------+-------------------+
1735 ``oo`` here means the SymPy infinity, which should compare greater than
1736 any integer.
1738 If you already know :py:meth:`~sympy.solvers.solvers.solve` cannot solve
1739 ``sol``, you can use ``trysolving=False`` to skip that step, which is the
1740 only potentially slow step. For example,
1741 :py:meth:`~sympy.solvers.ode.dsolve` with the ``simplify=False`` flag
1742 should do this.
1744 If ``sol`` is a list of solutions, if the worst solution in the list
1745 returns ``oo`` it returns that, otherwise it returns ``len(str(sol))``,
1746 that is, the length of the string representation of the whole list.
1748 Examples
1749 ========
1751 This function is designed to be passed to ``min`` as the key argument,
1752 such as ``min(listofsolutions, key=lambda i: ode_sol_simplicity(i,
1753 f(x)))``.
1755 >>> from sympy import symbols, Function, Eq, tan, Integral
1756 >>> from sympy.solvers.ode.ode import ode_sol_simplicity
1757 >>> x, C1, C2 = symbols('x, C1, C2')
1758 >>> f = Function('f')
1760 >>> ode_sol_simplicity(Eq(f(x), C1*x**2), f(x))
1761 -2
1762 >>> ode_sol_simplicity(Eq(x**2 + f(x), C1), f(x))
1763 -1
1764 >>> ode_sol_simplicity(Eq(f(x), C1*Integral(2*x, x)), f(x))
1765 oo
1766 >>> eq1 = Eq(f(x)/tan(f(x)/(2*x)), C1)
1767 >>> eq2 = Eq(f(x)/tan(f(x)/(2*x) + f(x)), C2)
1768 >>> [ode_sol_simplicity(eq, f(x)) for eq in [eq1, eq2]]
1769 [28, 35]
1770 >>> min([eq1, eq2], key=lambda i: ode_sol_simplicity(i, f(x)))
1771 Eq(f(x)/tan(f(x)/(2*x)), C1)
1773 """
1774 # TODO: if two solutions are solved for f(x), we still want to be
1775 # able to get the simpler of the two
1777 # See the docstring for the coercion rules. We check easier (faster)
1778 # things here first, to save time.
1780 if iterable(sol):
1781 # See if there are Integrals
1782 for i in sol:
1783 if ode_sol_simplicity(i, func, trysolving=trysolving) == oo:
1784 return oo
1786 return len(str(sol))
1788 if sol.has(Integral):
1789 return oo
1791 # Next, try to solve for func. This code will change slightly when CRootOf
1792 # is implemented in solve(). Probably a CRootOf solution should fall
1793 # somewhere between a normal solution and an unsolvable expression.
1795 # First, see if they are already solved
1796 if sol.lhs == func and not sol.rhs.has(func) or \
1797 sol.rhs == func and not sol.lhs.has(func):
1798 return -2
1799 # We are not so lucky, try solving manually
1800 if trysolving:
1801 try:
1802 sols = solve(sol, func)
1803 if not sols:
1804 raise NotImplementedError
1805 except NotImplementedError:
1806 pass
1807 else:
1808 return -1
1810 # Finally, a naive computation based on the length of the string version
1811 # of the expression. This may favor combined fractions because they
1812 # will not have duplicate denominators, and may slightly favor expressions
1813 # with fewer additions and subtractions, as those are separated by spaces
1814 # by the printer.
1816 # Additional ideas for simplicity heuristics are welcome, like maybe
1817 # checking if a equation has a larger domain, or if constantsimp has
1818 # introduced arbitrary constants numbered higher than the order of a
1819 # given ODE that sol is a solution of.
1820 return len(str(sol))
1823def _extract_funcs(eqs):
1824 funcs = []
1825 for eq in eqs:
1826 derivs = [node for node in preorder_traversal(eq) if isinstance(node, Derivative)]
1827 func = []
1828 for d in derivs:
1829 func += list(d.atoms(AppliedUndef))
1830 for func_ in func:
1831 funcs.append(func_)
1832 funcs = list(uniq(funcs))
1834 return funcs
1837def _get_constant_subexpressions(expr, Cs):
1838 Cs = set(Cs)
1839 Ces = []
1840 def _recursive_walk(expr):
1841 expr_syms = expr.free_symbols
1842 if expr_syms and expr_syms.issubset(Cs):
1843 Ces.append(expr)
1844 else:
1845 if expr.func == exp:
1846 expr = expr.expand(mul=True)
1847 if expr.func in (Add, Mul):
1848 d = sift(expr.args, lambda i : i.free_symbols.issubset(Cs))
1849 if len(d[True]) > 1:
1850 x = expr.func(*d[True])
1851 if not x.is_number:
1852 Ces.append(x)
1853 elif isinstance(expr, Integral):
1854 if expr.free_symbols.issubset(Cs) and \
1855 all(len(x) == 3 for x in expr.limits):
1856 Ces.append(expr)
1857 for i in expr.args:
1858 _recursive_walk(i)
1859 return
1860 _recursive_walk(expr)
1861 return Ces
1863def __remove_linear_redundancies(expr, Cs):
1864 cnts = {i: expr.count(i) for i in Cs}
1865 Cs = [i for i in Cs if cnts[i] > 0]
1867 def _linear(expr):
1868 if isinstance(expr, Add):
1869 xs = [i for i in Cs if expr.count(i)==cnts[i] \
1870 and 0 == expr.diff(i, 2)]
1871 d = {}
1872 for x in xs:
1873 y = expr.diff(x)
1874 if y not in d:
1875 d[y]=[]
1876 d[y].append(x)
1877 for y in d:
1878 if len(d[y]) > 1:
1879 d[y].sort(key=str)
1880 for x in d[y][1:]:
1881 expr = expr.subs(x, 0)
1882 return expr
1884 def _recursive_walk(expr):
1885 if len(expr.args) != 0:
1886 expr = expr.func(*[_recursive_walk(i) for i in expr.args])
1887 expr = _linear(expr)
1888 return expr
1890 if isinstance(expr, Equality):
1891 lhs, rhs = [_recursive_walk(i) for i in expr.args]
1892 f = lambda i: isinstance(i, Number) or i in Cs
1893 if isinstance(lhs, Symbol) and lhs in Cs:
1894 rhs, lhs = lhs, rhs
1895 if lhs.func in (Add, Symbol) and rhs.func in (Add, Symbol):
1896 dlhs = sift([lhs] if isinstance(lhs, AtomicExpr) else lhs.args, f)
1897 drhs = sift([rhs] if isinstance(rhs, AtomicExpr) else rhs.args, f)
1898 for i in [True, False]:
1899 for hs in [dlhs, drhs]:
1900 if i not in hs:
1901 hs[i] = [0]
1902 # this calculation can be simplified
1903 lhs = Add(*dlhs[False]) - Add(*drhs[False])
1904 rhs = Add(*drhs[True]) - Add(*dlhs[True])
1905 elif lhs.func in (Mul, Symbol) and rhs.func in (Mul, Symbol):
1906 dlhs = sift([lhs] if isinstance(lhs, AtomicExpr) else lhs.args, f)
1907 if True in dlhs:
1908 if False not in dlhs:
1909 dlhs[False] = [1]
1910 lhs = Mul(*dlhs[False])
1911 rhs = rhs/Mul(*dlhs[True])
1912 return Eq(lhs, rhs)
1913 else:
1914 return _recursive_walk(expr)
1916@vectorize(0)
1917def constantsimp(expr, constants):
1918 r"""
1919 Simplifies an expression with arbitrary constants in it.
1921 This function is written specifically to work with
1922 :py:meth:`~sympy.solvers.ode.dsolve`, and is not intended for general use.
1924 Simplification is done by "absorbing" the arbitrary constants into other
1925 arbitrary constants, numbers, and symbols that they are not independent
1926 of.
1928 The symbols must all have the same name with numbers after it, for
1929 example, ``C1``, ``C2``, ``C3``. The ``symbolname`` here would be
1930 '``C``', the ``startnumber`` would be 1, and the ``endnumber`` would be 3.
1931 If the arbitrary constants are independent of the variable ``x``, then the
1932 independent symbol would be ``x``. There is no need to specify the
1933 dependent function, such as ``f(x)``, because it already has the
1934 independent symbol, ``x``, in it.
1936 Because terms are "absorbed" into arbitrary constants and because
1937 constants are renumbered after simplifying, the arbitrary constants in
1938 expr are not necessarily equal to the ones of the same name in the
1939 returned result.
1941 If two or more arbitrary constants are added, multiplied, or raised to the
1942 power of each other, they are first absorbed together into a single
1943 arbitrary constant. Then the new constant is combined into other terms if
1944 necessary.
1946 Absorption of constants is done with limited assistance:
1948 1. terms of :py:class:`~sympy.core.add.Add`\s are collected to try join
1949 constants so `e^x (C_1 \cos(x) + C_2 \cos(x))` will simplify to `e^x
1950 C_1 \cos(x)`;
1952 2. powers with exponents that are :py:class:`~sympy.core.add.Add`\s are
1953 expanded so `e^{C_1 + x}` will be simplified to `C_1 e^x`.
1955 Use :py:meth:`~sympy.solvers.ode.ode.constant_renumber` to renumber constants
1956 after simplification or else arbitrary numbers on constants may appear,
1957 e.g. `C_1 + C_3 x`.
1959 In rare cases, a single constant can be "simplified" into two constants.
1960 Every differential equation solution should have as many arbitrary
1961 constants as the order of the differential equation. The result here will
1962 be technically correct, but it may, for example, have `C_1` and `C_2` in
1963 an expression, when `C_1` is actually equal to `C_2`. Use your discretion
1964 in such situations, and also take advantage of the ability to use hints in
1965 :py:meth:`~sympy.solvers.ode.dsolve`.
1967 Examples
1968 ========
1970 >>> from sympy import symbols
1971 >>> from sympy.solvers.ode.ode import constantsimp
1972 >>> C1, C2, C3, x, y = symbols('C1, C2, C3, x, y')
1973 >>> constantsimp(2*C1*x, {C1, C2, C3})
1974 C1*x
1975 >>> constantsimp(C1 + 2 + x, {C1, C2, C3})
1976 C1 + x
1977 >>> constantsimp(C1*C2 + 2 + C2 + C3*x, {C1, C2, C3})
1978 C1 + C3*x
1980 """
1981 # This function works recursively. The idea is that, for Mul,
1982 # Add, Pow, and Function, if the class has a constant in it, then
1983 # we can simplify it, which we do by recursing down and
1984 # simplifying up. Otherwise, we can skip that part of the
1985 # expression.
1987 Cs = constants
1989 orig_expr = expr
1991 constant_subexprs = _get_constant_subexpressions(expr, Cs)
1992 for xe in constant_subexprs:
1993 xes = list(xe.free_symbols)
1994 if not xes:
1995 continue
1996 if all(expr.count(c) == xe.count(c) for c in xes):
1997 xes.sort(key=str)
1998 expr = expr.subs(xe, xes[0])
2000 # try to perform common sub-expression elimination of constant terms
2001 try:
2002 commons, rexpr = cse(expr)
2003 commons.reverse()
2004 rexpr = rexpr[0]
2005 for s in commons:
2006 cs = list(s[1].atoms(Symbol))
2007 if len(cs) == 1 and cs[0] in Cs and \
2008 cs[0] not in rexpr.atoms(Symbol) and \
2009 not any(cs[0] in ex for ex in commons if ex != s):
2010 rexpr = rexpr.subs(s[0], cs[0])
2011 else:
2012 rexpr = rexpr.subs(*s)
2013 expr = rexpr
2014 except IndexError:
2015 pass
2016 expr = __remove_linear_redundancies(expr, Cs)
2018 def _conditional_term_factoring(expr):
2019 new_expr = terms_gcd(expr, clear=False, deep=True, expand=False)
2021 # we do not want to factor exponentials, so handle this separately
2022 if new_expr.is_Mul:
2023 infac = False
2024 asfac = False
2025 for m in new_expr.args:
2026 if isinstance(m, exp):
2027 asfac = True
2028 elif m.is_Add:
2029 infac = any(isinstance(fi, exp) for t in m.args
2030 for fi in Mul.make_args(t))
2031 if asfac and infac:
2032 new_expr = expr
2033 break
2034 return new_expr
2036 expr = _conditional_term_factoring(expr)
2038 # call recursively if more simplification is possible
2039 if orig_expr != expr:
2040 return constantsimp(expr, Cs)
2041 return expr
2044def constant_renumber(expr, variables=None, newconstants=None):
2045 r"""
2046 Renumber arbitrary constants in ``expr`` to use the symbol names as given
2047 in ``newconstants``. In the process, this reorders expression terms in a
2048 standard way.
2050 If ``newconstants`` is not provided then the new constant names will be
2051 ``C1``, ``C2`` etc. Otherwise ``newconstants`` should be an iterable
2052 giving the new symbols to use for the constants in order.
2054 The ``variables`` argument is a list of non-constant symbols. All other
2055 free symbols found in ``expr`` are assumed to be constants and will be
2056 renumbered. If ``variables`` is not given then any numbered symbol
2057 beginning with ``C`` (e.g. ``C1``) is assumed to be a constant.
2059 Symbols are renumbered based on ``.sort_key()``, so they should be
2060 numbered roughly in the order that they appear in the final, printed
2061 expression. Note that this ordering is based in part on hashes, so it can
2062 produce different results on different machines.
2064 The structure of this function is very similar to that of
2065 :py:meth:`~sympy.solvers.ode.constantsimp`.
2067 Examples
2068 ========
2070 >>> from sympy import symbols
2071 >>> from sympy.solvers.ode.ode import constant_renumber
2072 >>> x, C1, C2, C3 = symbols('x,C1:4')
2073 >>> expr = C3 + C2*x + C1*x**2
2074 >>> expr
2075 C1*x**2 + C2*x + C3
2076 >>> constant_renumber(expr)
2077 C1 + C2*x + C3*x**2
2079 The ``variables`` argument specifies which are constants so that the
2080 other symbols will not be renumbered:
2082 >>> constant_renumber(expr, [C1, x])
2083 C1*x**2 + C2 + C3*x
2085 The ``newconstants`` argument is used to specify what symbols to use when
2086 replacing the constants:
2088 >>> constant_renumber(expr, [x], newconstants=symbols('E1:4'))
2089 E1 + E2*x + E3*x**2
2091 """
2093 # System of expressions
2094 if isinstance(expr, (set, list, tuple)):
2095 return type(expr)(constant_renumber(Tuple(*expr),
2096 variables=variables, newconstants=newconstants))
2098 # Symbols in solution but not ODE are constants
2099 if variables is not None:
2100 variables = set(variables)
2101 free_symbols = expr.free_symbols
2102 constantsymbols = list(free_symbols - variables)
2103 # Any Cn is a constant...
2104 else:
2105 variables = set()
2106 isconstant = lambda s: s.startswith('C') and s[1:].isdigit()
2107 constantsymbols = [sym for sym in expr.free_symbols if isconstant(sym.name)]
2109 # Find new constants checking that they aren't already in the ODE
2110 if newconstants is None:
2111 iter_constants = numbered_symbols(start=1, prefix='C', exclude=variables)
2112 else:
2113 iter_constants = (sym for sym in newconstants if sym not in variables)
2115 constants_found = []
2117 # make a mapping to send all constantsymbols to S.One and use
2118 # that to make sure that term ordering is not dependent on
2119 # the indexed value of C
2120 C_1 = [(ci, S.One) for ci in constantsymbols]
2121 sort_key=lambda arg: default_sort_key(arg.subs(C_1))
2123 def _constant_renumber(expr):
2124 r"""
2125 We need to have an internal recursive function
2126 """
2128 # For system of expressions
2129 if isinstance(expr, Tuple):
2130 renumbered = [_constant_renumber(e) for e in expr]
2131 return Tuple(*renumbered)
2133 if isinstance(expr, Equality):
2134 return Eq(
2135 _constant_renumber(expr.lhs),
2136 _constant_renumber(expr.rhs))
2138 if type(expr) not in (Mul, Add, Pow) and not expr.is_Function and \
2139 not expr.has(*constantsymbols):
2140 # Base case, as above. Hope there aren't constants inside
2141 # of some other class, because they won't be renumbered.
2142 return expr
2143 elif expr.is_Piecewise:
2144 return expr
2145 elif expr in constantsymbols:
2146 if expr not in constants_found:
2147 constants_found.append(expr)
2148 return expr
2149 elif expr.is_Function or expr.is_Pow:
2150 return expr.func(
2151 *[_constant_renumber(x) for x in expr.args])
2152 else:
2153 sortedargs = list(expr.args)
2154 sortedargs.sort(key=sort_key)
2155 return expr.func(*[_constant_renumber(x) for x in sortedargs])
2156 expr = _constant_renumber(expr)
2158 # Don't renumber symbols present in the ODE.
2159 constants_found = [c for c in constants_found if c not in variables]
2161 # Renumbering happens here
2162 subs_dict = {var: cons for var, cons in zip(constants_found, iter_constants)}
2163 expr = expr.subs(subs_dict, simultaneous=True)
2165 return expr
2168def _handle_Integral(expr, func, hint):
2169 r"""
2170 Converts a solution with Integrals in it into an actual solution.
2172 For most hints, this simply runs ``expr.doit()``.
2174 """
2175 if hint == "nth_linear_constant_coeff_homogeneous":
2176 sol = expr
2177 elif not hint.endswith("_Integral"):
2178 sol = expr.doit()
2179 else:
2180 sol = expr
2181 return sol
2184# XXX: Should this function maybe go somewhere else?
2187def homogeneous_order(eq, *symbols):
2188 r"""
2189 Returns the order `n` if `g` is homogeneous and ``None`` if it is not
2190 homogeneous.
2192 Determines if a function is homogeneous and if so of what order. A
2193 function `f(x, y, \cdots)` is homogeneous of order `n` if `f(t x, t y,
2194 \cdots) = t^n f(x, y, \cdots)`.
2196 If the function is of two variables, `F(x, y)`, then `f` being homogeneous
2197 of any order is equivalent to being able to rewrite `F(x, y)` as `G(x/y)`
2198 or `H(y/x)`. This fact is used to solve 1st order ordinary differential
2199 equations whose coefficients are homogeneous of the same order (see the
2200 docstrings of
2201 :obj:`~sympy.solvers.ode.single.HomogeneousCoeffSubsDepDivIndep` and
2202 :obj:`~sympy.solvers.ode.single.HomogeneousCoeffSubsIndepDivDep`).
2204 Symbols can be functions, but every argument of the function must be a
2205 symbol, and the arguments of the function that appear in the expression
2206 must match those given in the list of symbols. If a declared function
2207 appears with different arguments than given in the list of symbols,
2208 ``None`` is returned.
2210 Examples
2211 ========
2213 >>> from sympy import Function, homogeneous_order, sqrt
2214 >>> from sympy.abc import x, y
2215 >>> f = Function('f')
2216 >>> homogeneous_order(f(x), f(x)) is None
2217 True
2218 >>> homogeneous_order(f(x,y), f(y, x), x, y) is None
2219 True
2220 >>> homogeneous_order(f(x), f(x), x)
2221 1
2222 >>> homogeneous_order(x**2*f(x)/sqrt(x**2+f(x)**2), x, f(x))
2223 2
2224 >>> homogeneous_order(x**2+f(x), x, f(x)) is None
2225 True
2227 """
2229 if not symbols:
2230 raise ValueError("homogeneous_order: no symbols were given.")
2231 symset = set(symbols)
2232 eq = sympify(eq)
2234 # The following are not supported
2235 if eq.has(Order, Derivative):
2236 return None
2238 # These are all constants
2239 if (eq.is_Number or
2240 eq.is_NumberSymbol or
2241 eq.is_number
2242 ):
2243 return S.Zero
2245 # Replace all functions with dummy variables
2246 dum = numbered_symbols(prefix='d', cls=Dummy)
2247 newsyms = set()
2248 for i in [j for j in symset if getattr(j, 'is_Function')]:
2249 iargs = set(i.args)
2250 if iargs.difference(symset):
2251 return None
2252 else:
2253 dummyvar = next(dum)
2254 eq = eq.subs(i, dummyvar)
2255 symset.remove(i)
2256 newsyms.add(dummyvar)
2257 symset.update(newsyms)
2259 if not eq.free_symbols & symset:
2260 return None
2262 # assuming order of a nested function can only be equal to zero
2263 if isinstance(eq, Function):
2264 return None if homogeneous_order(
2265 eq.args[0], *tuple(symset)) != 0 else S.Zero
2267 # make the replacement of x with x*t and see if t can be factored out
2268 t = Dummy('t', positive=True) # It is sufficient that t > 0
2269 eqs = separatevars(eq.subs([(i, t*i) for i in symset]), [t], dict=True)[t]
2270 if eqs is S.One:
2271 return S.Zero # there was no term with only t
2272 i, d = eqs.as_independent(t, as_Add=False)
2273 b, e = d.as_base_exp()
2274 if b == t:
2275 return e
2278def ode_2nd_power_series_ordinary(eq, func, order, match):
2279 r"""
2280 Gives a power series solution to a second order homogeneous differential
2281 equation with polynomial coefficients at an ordinary point. A homogeneous
2282 differential equation is of the form
2284 .. math :: P(x)\frac{d^2y}{dx^2} + Q(x)\frac{dy}{dx} + R(x) y(x) = 0
2286 For simplicity it is assumed that `P(x)`, `Q(x)` and `R(x)` are polynomials,
2287 it is sufficient that `\frac{Q(x)}{P(x)}` and `\frac{R(x)}{P(x)}` exists at
2288 `x_{0}`. A recurrence relation is obtained by substituting `y` as `\sum_{n=0}^\infty a_{n}x^{n}`,
2289 in the differential equation, and equating the nth term. Using this relation
2290 various terms can be generated.
2293 Examples
2294 ========
2296 >>> from sympy import dsolve, Function, pprint
2297 >>> from sympy.abc import x
2298 >>> f = Function("f")
2299 >>> eq = f(x).diff(x, 2) + f(x)
2300 >>> pprint(dsolve(eq, hint='2nd_power_series_ordinary'))
2301 / 4 2 \ / 2\
2302 |x x | | x | / 6\
2303 f(x) = C2*|-- - -- + 1| + C1*x*|1 - --| + O\x /
2304 \24 2 / \ 6 /
2307 References
2308 ==========
2309 - https://tutorial.math.lamar.edu/Classes/DE/SeriesSolutions.aspx
2310 - George E. Simmons, "Differential Equations with Applications and
2311 Historical Notes", p.p 176 - 184
2313 """
2314 x = func.args[0]
2315 f = func.func
2316 C0, C1 = get_numbered_constants(eq, num=2)
2317 n = Dummy("n", integer=True)
2318 s = Wild("s")
2319 k = Wild("k", exclude=[x])
2320 x0 = match['x0']
2321 terms = match['terms']
2322 p = match[match['a3']]
2323 q = match[match['b3']]
2324 r = match[match['c3']]
2325 seriesdict = {}
2326 recurr = Function("r")
2328 # Generating the recurrence relation which works this way:
2329 # for the second order term the summation begins at n = 2. The coefficients
2330 # p is multiplied with an*(n - 1)*(n - 2)*x**n-2 and a substitution is made such that
2331 # the exponent of x becomes n.
2332 # For example, if p is x, then the second degree recurrence term is
2333 # an*(n - 1)*(n - 2)*x**n-1, substituting (n - 1) as n, it transforms to
2334 # an+1*n*(n - 1)*x**n.
2335 # A similar process is done with the first order and zeroth order term.
2337 coefflist = [(recurr(n), r), (n*recurr(n), q), (n*(n - 1)*recurr(n), p)]
2338 for index, coeff in enumerate(coefflist):
2339 if coeff[1]:
2340 f2 = powsimp(expand((coeff[1]*(x - x0)**(n - index)).subs(x, x + x0)))
2341 if f2.is_Add:
2342 addargs = f2.args
2343 else:
2344 addargs = [f2]
2345 for arg in addargs:
2346 powm = arg.match(s*x**k)
2347 term = coeff[0]*powm[s]
2348 if not powm[k].is_Symbol:
2349 term = term.subs(n, n - powm[k].as_independent(n)[0])
2350 startind = powm[k].subs(n, index)
2351 # Seeing if the startterm can be reduced further.
2352 # If it vanishes for n lesser than startind, it is
2353 # equal to summation from n.
2354 if startind:
2355 for i in reversed(range(startind)):
2356 if not term.subs(n, i):
2357 seriesdict[term] = i
2358 else:
2359 seriesdict[term] = i + 1
2360 break
2361 else:
2362 seriesdict[term] = S.Zero
2364 # Stripping of terms so that the sum starts with the same number.
2365 teq = S.Zero
2366 suminit = seriesdict.values()
2367 rkeys = seriesdict.keys()
2368 req = Add(*rkeys)
2369 if any(suminit):
2370 maxval = max(suminit)
2371 for term in seriesdict:
2372 val = seriesdict[term]
2373 if val != maxval:
2374 for i in range(val, maxval):
2375 teq += term.subs(n, val)
2377 finaldict = {}
2378 if teq:
2379 fargs = teq.atoms(AppliedUndef)
2380 if len(fargs) == 1:
2381 finaldict[fargs.pop()] = 0
2382 else:
2383 maxf = max(fargs, key = lambda x: x.args[0])
2384 sol = solve(teq, maxf)
2385 if isinstance(sol, list):
2386 sol = sol[0]
2387 finaldict[maxf] = sol
2389 # Finding the recurrence relation in terms of the largest term.
2390 fargs = req.atoms(AppliedUndef)
2391 maxf = max(fargs, key = lambda x: x.args[0])
2392 minf = min(fargs, key = lambda x: x.args[0])
2393 if minf.args[0].is_Symbol:
2394 startiter = 0
2395 else:
2396 startiter = -minf.args[0].as_independent(n)[0]
2397 lhs = maxf
2398 rhs = solve(req, maxf)
2399 if isinstance(rhs, list):
2400 rhs = rhs[0]
2402 # Checking how many values are already present
2403 tcounter = len([t for t in finaldict.values() if t])
2405 for _ in range(tcounter, terms - 3): # Assuming c0 and c1 to be arbitrary
2406 check = rhs.subs(n, startiter)
2407 nlhs = lhs.subs(n, startiter)
2408 nrhs = check.subs(finaldict)
2409 finaldict[nlhs] = nrhs
2410 startiter += 1
2412 # Post processing
2413 series = C0 + C1*(x - x0)
2414 for term in finaldict:
2415 if finaldict[term]:
2416 fact = term.args[0]
2417 series += (finaldict[term].subs([(recurr(0), C0), (recurr(1), C1)])*(
2418 x - x0)**fact)
2419 series = collect(expand_mul(series), [C0, C1]) + Order(x**terms)
2420 return Eq(f(x), series)
2423def ode_2nd_power_series_regular(eq, func, order, match):
2424 r"""
2425 Gives a power series solution to a second order homogeneous differential
2426 equation with polynomial coefficients at a regular point. A second order
2427 homogeneous differential equation is of the form
2429 .. math :: P(x)\frac{d^2y}{dx^2} + Q(x)\frac{dy}{dx} + R(x) y(x) = 0
2431 A point is said to regular singular at `x0` if `x - x0\frac{Q(x)}{P(x)}`
2432 and `(x - x0)^{2}\frac{R(x)}{P(x)}` are analytic at `x0`. For simplicity
2433 `P(x)`, `Q(x)` and `R(x)` are assumed to be polynomials. The algorithm for
2434 finding the power series solutions is:
2436 1. Try expressing `(x - x0)P(x)` and `((x - x0)^{2})Q(x)` as power series
2437 solutions about x0. Find `p0` and `q0` which are the constants of the
2438 power series expansions.
2439 2. Solve the indicial equation `f(m) = m(m - 1) + m*p0 + q0`, to obtain the
2440 roots `m1` and `m2` of the indicial equation.
2441 3. If `m1 - m2` is a non integer there exists two series solutions. If
2442 `m1 = m2`, there exists only one solution. If `m1 - m2` is an integer,
2443 then the existence of one solution is confirmed. The other solution may
2444 or may not exist.
2446 The power series solution is of the form `x^{m}\sum_{n=0}^\infty a_{n}x^{n}`. The
2447 coefficients are determined by the following recurrence relation.
2448 `a_{n} = -\frac{\sum_{k=0}^{n-1} q_{n-k} + (m + k)p_{n-k}}{f(m + n)}`. For the case
2449 in which `m1 - m2` is an integer, it can be seen from the recurrence relation
2450 that for the lower root `m`, when `n` equals the difference of both the
2451 roots, the denominator becomes zero. So if the numerator is not equal to zero,
2452 a second series solution exists.
2455 Examples
2456 ========
2458 >>> from sympy import dsolve, Function, pprint
2459 >>> from sympy.abc import x
2460 >>> f = Function("f")
2461 >>> eq = x*(f(x).diff(x, 2)) + 2*(f(x).diff(x)) + x*f(x)
2462 >>> pprint(dsolve(eq, hint='2nd_power_series_regular'))
2463 / 6 4 2 \
2464 | x x x |
2465 / 4 2 \ C1*|- --- + -- - -- + 1|
2466 | x x | \ 720 24 2 / / 6\
2467 f(x) = C2*|--- - -- + 1| + ------------------------ + O\x /
2468 \120 6 / x
2471 References
2472 ==========
2473 - George E. Simmons, "Differential Equations with Applications and
2474 Historical Notes", p.p 176 - 184
2476 """
2477 x = func.args[0]
2478 f = func.func
2479 C0, C1 = get_numbered_constants(eq, num=2)
2480 m = Dummy("m") # for solving the indicial equation
2481 x0 = match['x0']
2482 terms = match['terms']
2483 p = match['p']
2484 q = match['q']
2486 # Generating the indicial equation
2487 indicial = []
2488 for term in [p, q]:
2489 if not term.has(x):
2490 indicial.append(term)
2491 else:
2492 term = series(term, x=x, n=1, x0=x0)
2493 if isinstance(term, Order):
2494 indicial.append(S.Zero)
2495 else:
2496 for arg in term.args:
2497 if not arg.has(x):
2498 indicial.append(arg)
2499 break
2501 p0, q0 = indicial
2502 sollist = solve(m*(m - 1) + m*p0 + q0, m)
2503 if sollist and isinstance(sollist, list) and all(
2504 sol.is_real for sol in sollist):
2505 serdict1 = {}
2506 serdict2 = {}
2507 if len(sollist) == 1:
2508 # Only one series solution exists in this case.
2509 m1 = m2 = sollist.pop()
2510 if terms-m1-1 <= 0:
2511 return Eq(f(x), Order(terms))
2512 serdict1 = _frobenius(terms-m1-1, m1, p0, q0, p, q, x0, x, C0)
2514 else:
2515 m1 = sollist[0]
2516 m2 = sollist[1]
2517 if m1 < m2:
2518 m1, m2 = m2, m1
2519 # Irrespective of whether m1 - m2 is an integer or not, one
2520 # Frobenius series solution exists.
2521 serdict1 = _frobenius(terms-m1-1, m1, p0, q0, p, q, x0, x, C0)
2522 if not (m1 - m2).is_integer:
2523 # Second frobenius series solution exists.
2524 serdict2 = _frobenius(terms-m2-1, m2, p0, q0, p, q, x0, x, C1)
2525 else:
2526 # Check if second frobenius series solution exists.
2527 serdict2 = _frobenius(terms-m2-1, m2, p0, q0, p, q, x0, x, C1, check=m1)
2529 if serdict1:
2530 finalseries1 = C0
2531 for key in serdict1:
2532 power = int(key.name[1:])
2533 finalseries1 += serdict1[key]*(x - x0)**power
2534 finalseries1 = (x - x0)**m1*finalseries1
2535 finalseries2 = S.Zero
2536 if serdict2:
2537 for key in serdict2:
2538 power = int(key.name[1:])
2539 finalseries2 += serdict2[key]*(x - x0)**power
2540 finalseries2 += C1
2541 finalseries2 = (x - x0)**m2*finalseries2
2542 return Eq(f(x), collect(finalseries1 + finalseries2,
2543 [C0, C1]) + Order(x**terms))
2546def _frobenius(n, m, p0, q0, p, q, x0, x, c, check=None):
2547 r"""
2548 Returns a dict with keys as coefficients and values as their values in terms of C0
2549 """
2550 n = int(n)
2551 # In cases where m1 - m2 is not an integer
2552 m2 = check
2554 d = Dummy("d")
2555 numsyms = numbered_symbols("C", start=0)
2556 numsyms = [next(numsyms) for i in range(n + 1)]
2557 serlist = []
2558 for ser in [p, q]:
2559 # Order term not present
2560 if ser.is_polynomial(x) and Poly(ser, x).degree() <= n:
2561 if x0:
2562 ser = ser.subs(x, x + x0)
2563 dict_ = Poly(ser, x).as_dict()
2564 # Order term present
2565 else:
2566 tseries = series(ser, x=x0, n=n+1)
2567 # Removing order
2568 dict_ = Poly(list(ordered(tseries.args))[: -1], x).as_dict()
2569 # Fill in with zeros, if coefficients are zero.
2570 for i in range(n + 1):
2571 if (i,) not in dict_:
2572 dict_[(i,)] = S.Zero
2573 serlist.append(dict_)
2575 pseries = serlist[0]
2576 qseries = serlist[1]
2577 indicial = d*(d - 1) + d*p0 + q0
2578 frobdict = {}
2579 for i in range(1, n + 1):
2580 num = c*(m*pseries[(i,)] + qseries[(i,)])
2581 for j in range(1, i):
2582 sym = Symbol("C" + str(j))
2583 num += frobdict[sym]*((m + j)*pseries[(i - j,)] + qseries[(i - j,)])
2585 # Checking for cases when m1 - m2 is an integer. If num equals zero
2586 # then a second Frobenius series solution cannot be found. If num is not zero
2587 # then set constant as zero and proceed.
2588 if m2 is not None and i == m2 - m:
2589 if num:
2590 return False
2591 else:
2592 frobdict[numsyms[i]] = S.Zero
2593 else:
2594 frobdict[numsyms[i]] = -num/(indicial.subs(d, m+i))
2596 return frobdict
2598def _remove_redundant_solutions(eq, solns, order, var):
2599 r"""
2600 Remove redundant solutions from the set of solutions.
2602 This function is needed because otherwise dsolve can return
2603 redundant solutions. As an example consider:
2605 eq = Eq((f(x).diff(x, 2))*f(x).diff(x), 0)
2607 There are two ways to find solutions to eq. The first is to solve f(x).diff(x, 2) = 0
2608 leading to solution f(x)=C1 + C2*x. The second is to solve the equation f(x).diff(x) = 0
2609 leading to the solution f(x) = C1. In this particular case we then see
2610 that the second solution is a special case of the first and we do not
2611 want to return it.
2613 This does not always happen. If we have
2615 eq = Eq((f(x)**2-4)*(f(x).diff(x)-4), 0)
2617 then we get the algebraic solution f(x) = [-2, 2] and the integral solution
2618 f(x) = x + C1 and in this case the two solutions are not equivalent wrt
2619 initial conditions so both should be returned.
2620 """
2621 def is_special_case_of(soln1, soln2):
2622 return _is_special_case_of(soln1, soln2, eq, order, var)
2624 unique_solns = []
2625 for soln1 in solns:
2626 for soln2 in unique_solns[:]:
2627 if is_special_case_of(soln1, soln2):
2628 break
2629 elif is_special_case_of(soln2, soln1):
2630 unique_solns.remove(soln2)
2631 else:
2632 unique_solns.append(soln1)
2634 return unique_solns
2636def _is_special_case_of(soln1, soln2, eq, order, var):
2637 r"""
2638 True if soln1 is found to be a special case of soln2 wrt some value of the
2639 constants that appear in soln2. False otherwise.
2640 """
2641 # The solutions returned by dsolve may be given explicitly or implicitly.
2642 # We will equate the sol1=(soln1.rhs - soln1.lhs), sol2=(soln2.rhs - soln2.lhs)
2643 # of the two solutions.
2644 #
2645 # Since this is supposed to hold for all x it also holds for derivatives.
2646 # For an order n ode we should be able to differentiate
2647 # each solution n times to get n+1 equations.
2648 #
2649 # We then try to solve those n+1 equations for the integrations constants
2650 # in sol2. If we can find a solution that does not depend on x then it
2651 # means that some value of the constants in sol1 is a special case of
2652 # sol2 corresponding to a particular choice of the integration constants.
2654 # In case the solution is in implicit form we subtract the sides
2655 soln1 = soln1.rhs - soln1.lhs
2656 soln2 = soln2.rhs - soln2.lhs
2658 # Work for the series solution
2659 if soln1.has(Order) and soln2.has(Order):
2660 if soln1.getO() == soln2.getO():
2661 soln1 = soln1.removeO()
2662 soln2 = soln2.removeO()
2663 else:
2664 return False
2665 elif soln1.has(Order) or soln2.has(Order):
2666 return False
2668 constants1 = soln1.free_symbols.difference(eq.free_symbols)
2669 constants2 = soln2.free_symbols.difference(eq.free_symbols)
2671 constants1_new = get_numbered_constants(Tuple(soln1, soln2), len(constants1))
2672 if len(constants1) == 1:
2673 constants1_new = {constants1_new}
2674 for c_old, c_new in zip(constants1, constants1_new):
2675 soln1 = soln1.subs(c_old, c_new)
2677 # n equations for sol1 = sol2, sol1'=sol2', ...
2678 lhs = soln1
2679 rhs = soln2
2680 eqns = [Eq(lhs, rhs)]
2681 for n in range(1, order):
2682 lhs = lhs.diff(var)
2683 rhs = rhs.diff(var)
2684 eq = Eq(lhs, rhs)
2685 eqns.append(eq)
2687 # BooleanTrue/False awkwardly show up for trivial equations
2688 if any(isinstance(eq, BooleanFalse) for eq in eqns):
2689 return False
2690 eqns = [eq for eq in eqns if not isinstance(eq, BooleanTrue)]
2692 try:
2693 constant_solns = solve(eqns, constants2)
2694 except NotImplementedError:
2695 return False
2697 # Sometimes returns a dict and sometimes a list of dicts
2698 if isinstance(constant_solns, dict):
2699 constant_solns = [constant_solns]
2701 # after solving the issue 17418, maybe we don't need the following checksol code.
2702 for constant_soln in constant_solns:
2703 for eq in eqns:
2704 eq=eq.rhs-eq.lhs
2705 if checksol(eq, constant_soln) is not True:
2706 return False
2708 # If any solution gives all constants as expressions that don't depend on
2709 # x then there exists constants for soln2 that give soln1
2710 for constant_soln in constant_solns:
2711 if not any(c.has(var) for c in constant_soln.values()):
2712 return True
2714 return False
2717def ode_1st_power_series(eq, func, order, match):
2718 r"""
2719 The power series solution is a method which gives the Taylor series expansion
2720 to the solution of a differential equation.
2722 For a first order differential equation `\frac{dy}{dx} = h(x, y)`, a power
2723 series solution exists at a point `x = x_{0}` if `h(x, y)` is analytic at `x_{0}`.
2724 The solution is given by
2726 .. math:: y(x) = y(x_{0}) + \sum_{n = 1}^{\infty} \frac{F_{n}(x_{0},b)(x - x_{0})^n}{n!},
2728 where `y(x_{0}) = b` is the value of y at the initial value of `x_{0}`.
2729 To compute the values of the `F_{n}(x_{0},b)` the following algorithm is
2730 followed, until the required number of terms are generated.
2732 1. `F_1 = h(x_{0}, b)`
2733 2. `F_{n+1} = \frac{\partial F_{n}}{\partial x} + \frac{\partial F_{n}}{\partial y}F_{1}`
2735 Examples
2736 ========
2738 >>> from sympy import Function, pprint, exp, dsolve
2739 >>> from sympy.abc import x
2740 >>> f = Function('f')
2741 >>> eq = exp(x)*(f(x).diff(x)) - f(x)
2742 >>> pprint(dsolve(eq, hint='1st_power_series'))
2743 3 4 5
2744 C1*x C1*x C1*x / 6\
2745 f(x) = C1 + C1*x - ----- + ----- + ----- + O\x /
2746 6 24 60
2749 References
2750 ==========
2752 - Travis W. Walker, Analytic power series technique for solving first-order
2753 differential equations, p.p 17, 18
2755 """
2756 x = func.args[0]
2757 y = match['y']
2758 f = func.func
2759 h = -match[match['d']]/match[match['e']]
2760 point = match['f0']
2761 value = match['f0val']
2762 terms = match['terms']
2764 # First term
2765 F = h
2766 if not h:
2767 return Eq(f(x), value)
2769 # Initialization
2770 series = value
2771 if terms > 1:
2772 hc = h.subs({x: point, y: value})
2773 if hc.has(oo) or hc.has(nan) or hc.has(zoo):
2774 # Derivative does not exist, not analytic
2775 return Eq(f(x), oo)
2776 elif hc:
2777 series += hc*(x - point)
2779 for factcount in range(2, terms):
2780 Fnew = F.diff(x) + F.diff(y)*h
2781 Fnewc = Fnew.subs({x: point, y: value})
2782 # Same logic as above
2783 if Fnewc.has(oo) or Fnewc.has(nan) or Fnewc.has(-oo) or Fnewc.has(zoo):
2784 return Eq(f(x), oo)
2785 series += Fnewc*((x - point)**factcount)/factorial(factcount)
2786 F = Fnew
2787 series += Order(x**terms)
2788 return Eq(f(x), series)
2791def checkinfsol(eq, infinitesimals, func=None, order=None):
2792 r"""
2793 This function is used to check if the given infinitesimals are the
2794 actual infinitesimals of the given first order differential equation.
2795 This method is specific to the Lie Group Solver of ODEs.
2797 As of now, it simply checks, by substituting the infinitesimals in the
2798 partial differential equation.
2801 .. math:: \frac{\partial \eta}{\partial x} + \left(\frac{\partial \eta}{\partial y}
2802 - \frac{\partial \xi}{\partial x}\right)*h
2803 - \frac{\partial \xi}{\partial y}*h^{2}
2804 - \xi\frac{\partial h}{\partial x} - \eta\frac{\partial h}{\partial y} = 0
2807 where `\eta`, and `\xi` are the infinitesimals and `h(x,y) = \frac{dy}{dx}`
2809 The infinitesimals should be given in the form of a list of dicts
2810 ``[{xi(x, y): inf, eta(x, y): inf}]``, corresponding to the
2811 output of the function infinitesimals. It returns a list
2812 of values of the form ``[(True/False, sol)]`` where ``sol`` is the value
2813 obtained after substituting the infinitesimals in the PDE. If it
2814 is ``True``, then ``sol`` would be 0.
2816 """
2817 if isinstance(eq, Equality):
2818 eq = eq.lhs - eq.rhs
2819 if not func:
2820 eq, func = _preprocess(eq)
2821 variables = func.args
2822 if len(variables) != 1:
2823 raise ValueError("ODE's have only one independent variable")
2824 else:
2825 x = variables[0]
2826 if not order:
2827 order = ode_order(eq, func)
2828 if order != 1:
2829 raise NotImplementedError("Lie groups solver has been implemented "
2830 "only for first order differential equations")
2831 else:
2832 df = func.diff(x)
2833 a = Wild('a', exclude = [df])
2834 b = Wild('b', exclude = [df])
2835 match = collect(expand(eq), df).match(a*df + b)
2837 if match:
2838 h = -simplify(match[b]/match[a])
2839 else:
2840 try:
2841 sol = solve(eq, df)
2842 except NotImplementedError:
2843 raise NotImplementedError("Infinitesimals for the "
2844 "first order ODE could not be found")
2845 else:
2846 h = sol[0] # Find infinitesimals for one solution
2848 y = Dummy('y')
2849 h = h.subs(func, y)
2850 xi = Function('xi')(x, y)
2851 eta = Function('eta')(x, y)
2852 dxi = Function('xi')(x, func)
2853 deta = Function('eta')(x, func)
2854 pde = (eta.diff(x) + (eta.diff(y) - xi.diff(x))*h -
2855 (xi.diff(y))*h**2 - xi*(h.diff(x)) - eta*(h.diff(y)))
2856 soltup = []
2857 for sol in infinitesimals:
2858 tsol = {xi: S(sol[dxi]).subs(func, y),
2859 eta: S(sol[deta]).subs(func, y)}
2860 sol = simplify(pde.subs(tsol).doit())
2861 if sol:
2862 soltup.append((False, sol.subs(y, func)))
2863 else:
2864 soltup.append((True, 0))
2865 return soltup
2868def sysode_linear_2eq_order1(match_):
2869 x = match_['func'][0].func
2870 y = match_['func'][1].func
2871 func = match_['func']
2872 fc = match_['func_coeff']
2873 eq = match_['eq']
2874 r = {}
2875 t = list(list(eq[0].atoms(Derivative))[0].atoms(Symbol))[0]
2876 for i in range(2):
2877 eq[i] = Add(*[terms/fc[i,func[i],1] for terms in Add.make_args(eq[i])])
2879 # for equations Eq(a1*diff(x(t),t), a*x(t) + b*y(t) + k1)
2880 # and Eq(a2*diff(x(t),t), c*x(t) + d*y(t) + k2)
2881 r['a'] = -fc[0,x(t),0]/fc[0,x(t),1]
2882 r['c'] = -fc[1,x(t),0]/fc[1,y(t),1]
2883 r['b'] = -fc[0,y(t),0]/fc[0,x(t),1]
2884 r['d'] = -fc[1,y(t),0]/fc[1,y(t),1]
2885 forcing = [S.Zero,S.Zero]
2886 for i in range(2):
2887 for j in Add.make_args(eq[i]):
2888 if not j.has(x(t), y(t)):
2889 forcing[i] += j
2890 if not (forcing[0].has(t) or forcing[1].has(t)):
2891 r['k1'] = forcing[0]
2892 r['k2'] = forcing[1]
2893 else:
2894 raise NotImplementedError("Only homogeneous problems are supported" +
2895 " (and constant inhomogeneity)")
2897 if match_['type_of_equation'] == 'type6':
2898 sol = _linear_2eq_order1_type6(x, y, t, r, eq)
2899 if match_['type_of_equation'] == 'type7':
2900 sol = _linear_2eq_order1_type7(x, y, t, r, eq)
2901 return sol
2903def _linear_2eq_order1_type6(x, y, t, r, eq):
2904 r"""
2905 The equations of this type of ode are .
2907 .. math:: x' = f(t) x + g(t) y
2909 .. math:: y' = a [f(t) + a h(t)] x + a [g(t) - h(t)] y
2911 This is solved by first multiplying the first equation by `-a` and adding
2912 it to the second equation to obtain
2914 .. math:: y' - a x' = -a h(t) (y - a x)
2916 Setting `U = y - ax` and integrating the equation we arrive at
2918 .. math:: y - ax = C_1 e^{-a \int h(t) \,dt}
2920 and on substituting the value of y in first equation give rise to first order ODEs. After solving for
2921 `x`, we can obtain `y` by substituting the value of `x` in second equation.
2923 """
2924 C1, C2, C3, C4 = get_numbered_constants(eq, num=4)
2925 p = 0
2926 q = 0
2927 p1 = cancel(r['c']/cancel(r['c']/r['d']).as_numer_denom()[0])
2928 p2 = cancel(r['a']/cancel(r['a']/r['b']).as_numer_denom()[0])
2929 for n, i in enumerate([p1, p2]):
2930 for j in Mul.make_args(collect_const(i)):
2931 if not j.has(t):
2932 q = j
2933 if q!=0 and n==0:
2934 if ((r['c']/j - r['a'])/(r['b'] - r['d']/j)) == j:
2935 p = 1
2936 s = j
2937 break
2938 if q!=0 and n==1:
2939 if ((r['a']/j - r['c'])/(r['d'] - r['b']/j)) == j:
2940 p = 2
2941 s = j
2942 break
2944 if p == 1:
2945 equ = diff(x(t),t) - r['a']*x(t) - r['b']*(s*x(t) + C1*exp(-s*Integral(r['b'] - r['d']/s, t)))
2946 hint1 = classify_ode(equ)[1]
2947 sol1 = dsolve(equ, hint=hint1+'_Integral').rhs
2948 sol2 = s*sol1 + C1*exp(-s*Integral(r['b'] - r['d']/s, t))
2949 elif p ==2:
2950 equ = diff(y(t),t) - r['c']*y(t) - r['d']*s*y(t) + C1*exp(-s*Integral(r['d'] - r['b']/s, t))
2951 hint1 = classify_ode(equ)[1]
2952 sol2 = dsolve(equ, hint=hint1+'_Integral').rhs
2953 sol1 = s*sol2 + C1*exp(-s*Integral(r['d'] - r['b']/s, t))
2954 return [Eq(x(t), sol1), Eq(y(t), sol2)]
2956def _linear_2eq_order1_type7(x, y, t, r, eq):
2957 r"""
2958 The equations of this type of ode are .
2960 .. math:: x' = f(t) x + g(t) y
2962 .. math:: y' = h(t) x + p(t) y
2964 Differentiating the first equation and substituting the value of `y`
2965 from second equation will give a second-order linear equation
2967 .. math:: g x'' - (fg + gp + g') x' + (fgp - g^{2} h + f g' - f' g) x = 0
2969 This above equation can be easily integrated if following conditions are satisfied.
2971 1. `fgp - g^{2} h + f g' - f' g = 0`
2973 2. `fgp - g^{2} h + f g' - f' g = ag, fg + gp + g' = bg`
2975 If first condition is satisfied then it is solved by current dsolve solver and in second case it becomes
2976 a constant coefficient differential equation which is also solved by current solver.
2978 Otherwise if the above condition fails then,
2979 a particular solution is assumed as `x = x_0(t)` and `y = y_0(t)`
2980 Then the general solution is expressed as
2982 .. math:: x = C_1 x_0(t) + C_2 x_0(t) \int \frac{g(t) F(t) P(t)}{x_0^{2}(t)} \,dt
2984 .. math:: y = C_1 y_0(t) + C_2 [\frac{F(t) P(t)}{x_0(t)} + y_0(t) \int \frac{g(t) F(t) P(t)}{x_0^{2}(t)} \,dt]
2986 where C1 and C2 are arbitrary constants and
2988 .. math:: F(t) = e^{\int f(t) \,dt}, P(t) = e^{\int p(t) \,dt}
2990 """
2991 C1, C2, C3, C4 = get_numbered_constants(eq, num=4)
2992 e1 = r['a']*r['b']*r['c'] - r['b']**2*r['c'] + r['a']*diff(r['b'],t) - diff(r['a'],t)*r['b']
2993 e2 = r['a']*r['c']*r['d'] - r['b']*r['c']**2 + diff(r['c'],t)*r['d'] - r['c']*diff(r['d'],t)
2994 m1 = r['a']*r['b'] + r['b']*r['d'] + diff(r['b'],t)
2995 m2 = r['a']*r['c'] + r['c']*r['d'] + diff(r['c'],t)
2996 if e1 == 0:
2997 sol1 = dsolve(r['b']*diff(x(t),t,t) - m1*diff(x(t),t)).rhs
2998 sol2 = dsolve(diff(y(t),t) - r['c']*sol1 - r['d']*y(t)).rhs
2999 elif e2 == 0:
3000 sol2 = dsolve(r['c']*diff(y(t),t,t) - m2*diff(y(t),t)).rhs
3001 sol1 = dsolve(diff(x(t),t) - r['a']*x(t) - r['b']*sol2).rhs
3002 elif not (e1/r['b']).has(t) and not (m1/r['b']).has(t):
3003 sol1 = dsolve(diff(x(t),t,t) - (m1/r['b'])*diff(x(t),t) - (e1/r['b'])*x(t)).rhs
3004 sol2 = dsolve(diff(y(t),t) - r['c']*sol1 - r['d']*y(t)).rhs
3005 elif not (e2/r['c']).has(t) and not (m2/r['c']).has(t):
3006 sol2 = dsolve(diff(y(t),t,t) - (m2/r['c'])*diff(y(t),t) - (e2/r['c'])*y(t)).rhs
3007 sol1 = dsolve(diff(x(t),t) - r['a']*x(t) - r['b']*sol2).rhs
3008 else:
3009 x0 = Function('x0')(t) # x0 and y0 being particular solutions
3010 y0 = Function('y0')(t)
3011 F = exp(Integral(r['a'],t))
3012 P = exp(Integral(r['d'],t))
3013 sol1 = C1*x0 + C2*x0*Integral(r['b']*F*P/x0**2, t)
3014 sol2 = C1*y0 + C2*(F*P/x0 + y0*Integral(r['b']*F*P/x0**2, t))
3015 return [Eq(x(t), sol1), Eq(y(t), sol2)]
3018def sysode_nonlinear_2eq_order1(match_):
3019 func = match_['func']
3020 eq = match_['eq']
3021 fc = match_['func_coeff']
3022 t = list(list(eq[0].atoms(Derivative))[0].atoms(Symbol))[0]
3023 if match_['type_of_equation'] == 'type5':
3024 sol = _nonlinear_2eq_order1_type5(func, t, eq)
3025 return sol
3026 x = func[0].func
3027 y = func[1].func
3028 for i in range(2):
3029 eqs = 0
3030 for terms in Add.make_args(eq[i]):
3031 eqs += terms/fc[i,func[i],1]
3032 eq[i] = eqs
3033 if match_['type_of_equation'] == 'type1':
3034 sol = _nonlinear_2eq_order1_type1(x, y, t, eq)
3035 elif match_['type_of_equation'] == 'type2':
3036 sol = _nonlinear_2eq_order1_type2(x, y, t, eq)
3037 elif match_['type_of_equation'] == 'type3':
3038 sol = _nonlinear_2eq_order1_type3(x, y, t, eq)
3039 elif match_['type_of_equation'] == 'type4':
3040 sol = _nonlinear_2eq_order1_type4(x, y, t, eq)
3041 return sol
3044def _nonlinear_2eq_order1_type1(x, y, t, eq):
3045 r"""
3046 Equations:
3048 .. math:: x' = x^n F(x,y)
3050 .. math:: y' = g(y) F(x,y)
3052 Solution:
3054 .. math:: x = \varphi(y), \int \frac{1}{g(y) F(\varphi(y),y)} \,dy = t + C_2
3056 where
3058 if `n \neq 1`
3060 .. math:: \varphi = [C_1 + (1-n) \int \frac{1}{g(y)} \,dy]^{\frac{1}{1-n}}
3062 if `n = 1`
3064 .. math:: \varphi = C_1 e^{\int \frac{1}{g(y)} \,dy}
3066 where `C_1` and `C_2` are arbitrary constants.
3068 """
3069 C1, C2 = get_numbered_constants(eq, num=2)
3070 n = Wild('n', exclude=[x(t),y(t)])
3071 f = Wild('f')
3072 u, v = symbols('u, v')
3073 r = eq[0].match(diff(x(t),t) - x(t)**n*f)
3074 g = ((diff(y(t),t) - eq[1])/r[f]).subs(y(t),v)
3075 F = r[f].subs(x(t),u).subs(y(t),v)
3076 n = r[n]
3077 if n!=1:
3078 phi = (C1 + (1-n)*Integral(1/g, v))**(1/(1-n))
3079 else:
3080 phi = C1*exp(Integral(1/g, v))
3081 phi = phi.doit()
3082 sol2 = solve(Integral(1/(g*F.subs(u,phi)), v).doit() - t - C2, v)
3083 sol = []
3084 for sols in sol2:
3085 sol.append(Eq(x(t),phi.subs(v, sols)))
3086 sol.append(Eq(y(t), sols))
3087 return sol
3089def _nonlinear_2eq_order1_type2(x, y, t, eq):
3090 r"""
3091 Equations:
3093 .. math:: x' = e^{\lambda x} F(x,y)
3095 .. math:: y' = g(y) F(x,y)
3097 Solution:
3099 .. math:: x = \varphi(y), \int \frac{1}{g(y) F(\varphi(y),y)} \,dy = t + C_2
3101 where
3103 if `\lambda \neq 0`
3105 .. math:: \varphi = -\frac{1}{\lambda} log(C_1 - \lambda \int \frac{1}{g(y)} \,dy)
3107 if `\lambda = 0`
3109 .. math:: \varphi = C_1 + \int \frac{1}{g(y)} \,dy
3111 where `C_1` and `C_2` are arbitrary constants.
3113 """
3114 C1, C2 = get_numbered_constants(eq, num=2)
3115 n = Wild('n', exclude=[x(t),y(t)])
3116 f = Wild('f')
3117 u, v = symbols('u, v')
3118 r = eq[0].match(diff(x(t),t) - exp(n*x(t))*f)
3119 g = ((diff(y(t),t) - eq[1])/r[f]).subs(y(t),v)
3120 F = r[f].subs(x(t),u).subs(y(t),v)
3121 n = r[n]
3122 if n:
3123 phi = -1/n*log(C1 - n*Integral(1/g, v))
3124 else:
3125 phi = C1 + Integral(1/g, v)
3126 phi = phi.doit()
3127 sol2 = solve(Integral(1/(g*F.subs(u,phi)), v).doit() - t - C2, v)
3128 sol = []
3129 for sols in sol2:
3130 sol.append(Eq(x(t),phi.subs(v, sols)))
3131 sol.append(Eq(y(t), sols))
3132 return sol
3134def _nonlinear_2eq_order1_type3(x, y, t, eq):
3135 r"""
3136 Autonomous system of general form
3138 .. math:: x' = F(x,y)
3140 .. math:: y' = G(x,y)
3142 Assuming `y = y(x, C_1)` where `C_1` is an arbitrary constant is the general
3143 solution of the first-order equation
3145 .. math:: F(x,y) y'_x = G(x,y)
3147 Then the general solution of the original system of equations has the form
3149 .. math:: \int \frac{1}{F(x,y(x,C_1))} \,dx = t + C_1
3151 """
3152 C1, C2, C3, C4 = get_numbered_constants(eq, num=4)
3153 v = Function('v')
3154 u = Symbol('u')
3155 f = Wild('f')
3156 g = Wild('g')
3157 r1 = eq[0].match(diff(x(t),t) - f)
3158 r2 = eq[1].match(diff(y(t),t) - g)
3159 F = r1[f].subs(x(t), u).subs(y(t), v(u))
3160 G = r2[g].subs(x(t), u).subs(y(t), v(u))
3161 sol2r = dsolve(Eq(diff(v(u), u), G/F))
3162 if isinstance(sol2r, Equality):
3163 sol2r = [sol2r]
3164 for sol2s in sol2r:
3165 sol1 = solve(Integral(1/F.subs(v(u), sol2s.rhs), u).doit() - t - C2, u)
3166 sol = []
3167 for sols in sol1:
3168 sol.append(Eq(x(t), sols))
3169 sol.append(Eq(y(t), (sol2s.rhs).subs(u, sols)))
3170 return sol
3172def _nonlinear_2eq_order1_type4(x, y, t, eq):
3173 r"""
3174 Equation:
3176 .. math:: x' = f_1(x) g_1(y) \phi(x,y,t)
3178 .. math:: y' = f_2(x) g_2(y) \phi(x,y,t)
3180 First integral:
3182 .. math:: \int \frac{f_2(x)}{f_1(x)} \,dx - \int \frac{g_1(y)}{g_2(y)} \,dy = C
3184 where `C` is an arbitrary constant.
3186 On solving the first integral for `x` (resp., `y` ) and on substituting the
3187 resulting expression into either equation of the original solution, one
3188 arrives at a first-order equation for determining `y` (resp., `x` ).
3190 """
3191 C1, C2 = get_numbered_constants(eq, num=2)
3192 u, v = symbols('u, v')
3193 U, V = symbols('U, V', cls=Function)
3194 f = Wild('f')
3195 g = Wild('g')
3196 f1 = Wild('f1', exclude=[v,t])
3197 f2 = Wild('f2', exclude=[v,t])
3198 g1 = Wild('g1', exclude=[u,t])
3199 g2 = Wild('g2', exclude=[u,t])
3200 r1 = eq[0].match(diff(x(t),t) - f)
3201 r2 = eq[1].match(diff(y(t),t) - g)
3202 num, den = (
3203 (r1[f].subs(x(t),u).subs(y(t),v))/
3204 (r2[g].subs(x(t),u).subs(y(t),v))).as_numer_denom()
3205 R1 = num.match(f1*g1)
3206 R2 = den.match(f2*g2)
3207 phi = (r1[f].subs(x(t),u).subs(y(t),v))/num
3208 F1 = R1[f1]; F2 = R2[f2]
3209 G1 = R1[g1]; G2 = R2[g2]
3210 sol1r = solve(Integral(F2/F1, u).doit() - Integral(G1/G2,v).doit() - C1, u)
3211 sol2r = solve(Integral(F2/F1, u).doit() - Integral(G1/G2,v).doit() - C1, v)
3212 sol = []
3213 for sols in sol1r:
3214 sol.append(Eq(y(t), dsolve(diff(V(t),t) - F2.subs(u,sols).subs(v,V(t))*G2.subs(v,V(t))*phi.subs(u,sols).subs(v,V(t))).rhs))
3215 for sols in sol2r:
3216 sol.append(Eq(x(t), dsolve(diff(U(t),t) - F1.subs(u,U(t))*G1.subs(v,sols).subs(u,U(t))*phi.subs(v,sols).subs(u,U(t))).rhs))
3217 return set(sol)
3219def _nonlinear_2eq_order1_type5(func, t, eq):
3220 r"""
3221 Clairaut system of ODEs
3223 .. math:: x = t x' + F(x',y')
3225 .. math:: y = t y' + G(x',y')
3227 The following are solutions of the system
3229 `(i)` straight lines:
3231 .. math:: x = C_1 t + F(C_1, C_2), y = C_2 t + G(C_1, C_2)
3233 where `C_1` and `C_2` are arbitrary constants;
3235 `(ii)` envelopes of the above lines;
3237 `(iii)` continuously differentiable lines made up from segments of the lines
3238 `(i)` and `(ii)`.
3240 """
3241 C1, C2 = get_numbered_constants(eq, num=2)
3242 f = Wild('f')
3243 g = Wild('g')
3244 def check_type(x, y):
3245 r1 = eq[0].match(t*diff(x(t),t) - x(t) + f)
3246 r2 = eq[1].match(t*diff(y(t),t) - y(t) + g)
3247 if not (r1 and r2):
3248 r1 = eq[0].match(diff(x(t),t) - x(t)/t + f/t)
3249 r2 = eq[1].match(diff(y(t),t) - y(t)/t + g/t)
3250 if not (r1 and r2):
3251 r1 = (-eq[0]).match(t*diff(x(t),t) - x(t) + f)
3252 r2 = (-eq[1]).match(t*diff(y(t),t) - y(t) + g)
3253 if not (r1 and r2):
3254 r1 = (-eq[0]).match(diff(x(t),t) - x(t)/t + f/t)
3255 r2 = (-eq[1]).match(diff(y(t),t) - y(t)/t + g/t)
3256 return [r1, r2]
3257 for func_ in func:
3258 if isinstance(func_, list):
3259 x = func[0][0].func
3260 y = func[0][1].func
3261 [r1, r2] = check_type(x, y)
3262 if not (r1 and r2):
3263 [r1, r2] = check_type(y, x)
3264 x, y = y, x
3265 x1 = diff(x(t),t); y1 = diff(y(t),t)
3266 return {Eq(x(t), C1*t + r1[f].subs(x1,C1).subs(y1,C2)), Eq(y(t), C2*t + r2[g].subs(x1,C1).subs(y1,C2))}
3268def sysode_nonlinear_3eq_order1(match_):
3269 x = match_['func'][0].func
3270 y = match_['func'][1].func
3271 z = match_['func'][2].func
3272 eq = match_['eq']
3273 t = list(list(eq[0].atoms(Derivative))[0].atoms(Symbol))[0]
3274 if match_['type_of_equation'] == 'type1':
3275 sol = _nonlinear_3eq_order1_type1(x, y, z, t, eq)
3276 if match_['type_of_equation'] == 'type2':
3277 sol = _nonlinear_3eq_order1_type2(x, y, z, t, eq)
3278 if match_['type_of_equation'] == 'type3':
3279 sol = _nonlinear_3eq_order1_type3(x, y, z, t, eq)
3280 if match_['type_of_equation'] == 'type4':
3281 sol = _nonlinear_3eq_order1_type4(x, y, z, t, eq)
3282 if match_['type_of_equation'] == 'type5':
3283 sol = _nonlinear_3eq_order1_type5(x, y, z, t, eq)
3284 return sol
3286def _nonlinear_3eq_order1_type1(x, y, z, t, eq):
3287 r"""
3288 Equations:
3290 .. math:: a x' = (b - c) y z, \enspace b y' = (c - a) z x, \enspace c z' = (a - b) x y
3292 First Integrals:
3294 .. math:: a x^{2} + b y^{2} + c z^{2} = C_1
3296 .. math:: a^{2} x^{2} + b^{2} y^{2} + c^{2} z^{2} = C_2
3298 where `C_1` and `C_2` are arbitrary constants. On solving the integrals for `y` and
3299 `z` and on substituting the resulting expressions into the first equation of the
3300 system, we arrives at a separable first-order equation on `x`. Similarly doing that
3301 for other two equations, we will arrive at first order equation on `y` and `z` too.
3303 References
3304 ==========
3305 -https://eqworld.ipmnet.ru/en/solutions/sysode/sode0401.pdf
3307 """
3308 C1, C2 = get_numbered_constants(eq, num=2)
3309 u, v, w = symbols('u, v, w')
3310 p = Wild('p', exclude=[x(t), y(t), z(t), t])
3311 q = Wild('q', exclude=[x(t), y(t), z(t), t])
3312 s = Wild('s', exclude=[x(t), y(t), z(t), t])
3313 r = (diff(x(t),t) - eq[0]).match(p*y(t)*z(t))
3314 r.update((diff(y(t),t) - eq[1]).match(q*z(t)*x(t)))
3315 r.update((diff(z(t),t) - eq[2]).match(s*x(t)*y(t)))
3316 n1, d1 = r[p].as_numer_denom()
3317 n2, d2 = r[q].as_numer_denom()
3318 n3, d3 = r[s].as_numer_denom()
3319 val = solve([n1*u-d1*v+d1*w, d2*u+n2*v-d2*w, d3*u-d3*v-n3*w],[u,v])
3320 vals = [val[v], val[u]]
3321 c = lcm(vals[0].as_numer_denom()[1], vals[1].as_numer_denom()[1])
3322 b = vals[0].subs(w, c)
3323 a = vals[1].subs(w, c)
3324 y_x = sqrt(((c*C1-C2) - a*(c-a)*x(t)**2)/(b*(c-b)))
3325 z_x = sqrt(((b*C1-C2) - a*(b-a)*x(t)**2)/(c*(b-c)))
3326 z_y = sqrt(((a*C1-C2) - b*(a-b)*y(t)**2)/(c*(a-c)))
3327 x_y = sqrt(((c*C1-C2) - b*(c-b)*y(t)**2)/(a*(c-a)))
3328 x_z = sqrt(((b*C1-C2) - c*(b-c)*z(t)**2)/(a*(b-a)))
3329 y_z = sqrt(((a*C1-C2) - c*(a-c)*z(t)**2)/(b*(a-b)))
3330 sol1 = dsolve(a*diff(x(t),t) - (b-c)*y_x*z_x)
3331 sol2 = dsolve(b*diff(y(t),t) - (c-a)*z_y*x_y)
3332 sol3 = dsolve(c*diff(z(t),t) - (a-b)*x_z*y_z)
3333 return [sol1, sol2, sol3]
3336def _nonlinear_3eq_order1_type2(x, y, z, t, eq):
3337 r"""
3338 Equations:
3340 .. math:: a x' = (b - c) y z f(x, y, z, t)
3342 .. math:: b y' = (c - a) z x f(x, y, z, t)
3344 .. math:: c z' = (a - b) x y f(x, y, z, t)
3346 First Integrals:
3348 .. math:: a x^{2} + b y^{2} + c z^{2} = C_1
3350 .. math:: a^{2} x^{2} + b^{2} y^{2} + c^{2} z^{2} = C_2
3352 where `C_1` and `C_2` are arbitrary constants. On solving the integrals for `y` and
3353 `z` and on substituting the resulting expressions into the first equation of the
3354 system, we arrives at a first-order differential equations on `x`. Similarly doing
3355 that for other two equations we will arrive at first order equation on `y` and `z`.
3357 References
3358 ==========
3359 -https://eqworld.ipmnet.ru/en/solutions/sysode/sode0402.pdf
3361 """
3362 C1, C2 = get_numbered_constants(eq, num=2)
3363 u, v, w = symbols('u, v, w')
3364 p = Wild('p', exclude=[x(t), y(t), z(t), t])
3365 q = Wild('q', exclude=[x(t), y(t), z(t), t])
3366 s = Wild('s', exclude=[x(t), y(t), z(t), t])
3367 f = Wild('f')
3368 r1 = (diff(x(t),t) - eq[0]).match(y(t)*z(t)*f)
3369 r = collect_const(r1[f]).match(p*f)
3370 r.update(((diff(y(t),t) - eq[1])/r[f]).match(q*z(t)*x(t)))
3371 r.update(((diff(z(t),t) - eq[2])/r[f]).match(s*x(t)*y(t)))
3372 n1, d1 = r[p].as_numer_denom()
3373 n2, d2 = r[q].as_numer_denom()
3374 n3, d3 = r[s].as_numer_denom()
3375 val = solve([n1*u-d1*v+d1*w, d2*u+n2*v-d2*w, -d3*u+d3*v+n3*w],[u,v])
3376 vals = [val[v], val[u]]
3377 c = lcm(vals[0].as_numer_denom()[1], vals[1].as_numer_denom()[1])
3378 a = vals[0].subs(w, c)
3379 b = vals[1].subs(w, c)
3380 y_x = sqrt(((c*C1-C2) - a*(c-a)*x(t)**2)/(b*(c-b)))
3381 z_x = sqrt(((b*C1-C2) - a*(b-a)*x(t)**2)/(c*(b-c)))
3382 z_y = sqrt(((a*C1-C2) - b*(a-b)*y(t)**2)/(c*(a-c)))
3383 x_y = sqrt(((c*C1-C2) - b*(c-b)*y(t)**2)/(a*(c-a)))
3384 x_z = sqrt(((b*C1-C2) - c*(b-c)*z(t)**2)/(a*(b-a)))
3385 y_z = sqrt(((a*C1-C2) - c*(a-c)*z(t)**2)/(b*(a-b)))
3386 sol1 = dsolve(a*diff(x(t),t) - (b-c)*y_x*z_x*r[f])
3387 sol2 = dsolve(b*diff(y(t),t) - (c-a)*z_y*x_y*r[f])
3388 sol3 = dsolve(c*diff(z(t),t) - (a-b)*x_z*y_z*r[f])
3389 return [sol1, sol2, sol3]
3391def _nonlinear_3eq_order1_type3(x, y, z, t, eq):
3392 r"""
3393 Equations:
3395 .. math:: x' = c F_2 - b F_3, \enspace y' = a F_3 - c F_1, \enspace z' = b F_1 - a F_2
3397 where `F_n = F_n(x, y, z, t)`.
3399 1. First Integral:
3401 .. math:: a x + b y + c z = C_1,
3403 where C is an arbitrary constant.
3405 2. If we assume function `F_n` to be independent of `t`,i.e, `F_n` = `F_n (x, y, z)`
3406 Then, on eliminating `t` and `z` from the first two equation of the system, one
3407 arrives at the first-order equation
3409 .. math:: \frac{dy}{dx} = \frac{a F_3 (x, y, z) - c F_1 (x, y, z)}{c F_2 (x, y, z) -
3410 b F_3 (x, y, z)}
3412 where `z = \frac{1}{c} (C_1 - a x - b y)`
3414 References
3415 ==========
3416 -https://eqworld.ipmnet.ru/en/solutions/sysode/sode0404.pdf
3418 """
3419 C1 = get_numbered_constants(eq, num=1)
3420 u, v, w = symbols('u, v, w')
3421 fu, fv, fw = symbols('u, v, w', cls=Function)
3422 p = Wild('p', exclude=[x(t), y(t), z(t), t])
3423 q = Wild('q', exclude=[x(t), y(t), z(t), t])
3424 s = Wild('s', exclude=[x(t), y(t), z(t), t])
3425 F1, F2, F3 = symbols('F1, F2, F3', cls=Wild)
3426 r1 = (diff(x(t), t) - eq[0]).match(F2-F3)
3427 r = collect_const(r1[F2]).match(s*F2)
3428 r.update(collect_const(r1[F3]).match(q*F3))
3429 if eq[1].has(r[F2]) and not eq[1].has(r[F3]):
3430 r[F2], r[F3] = r[F3], r[F2]
3431 r[s], r[q] = -r[q], -r[s]
3432 r.update((diff(y(t), t) - eq[1]).match(p*r[F3] - r[s]*F1))
3433 a = r[p]; b = r[q]; c = r[s]
3434 F1 = r[F1].subs(x(t), u).subs(y(t),v).subs(z(t), w)
3435 F2 = r[F2].subs(x(t), u).subs(y(t),v).subs(z(t), w)
3436 F3 = r[F3].subs(x(t), u).subs(y(t),v).subs(z(t), w)
3437 z_xy = (C1-a*u-b*v)/c
3438 y_zx = (C1-a*u-c*w)/b
3439 x_yz = (C1-b*v-c*w)/a
3440 y_x = dsolve(diff(fv(u),u) - ((a*F3-c*F1)/(c*F2-b*F3)).subs(w,z_xy).subs(v,fv(u))).rhs
3441 z_x = dsolve(diff(fw(u),u) - ((b*F1-a*F2)/(c*F2-b*F3)).subs(v,y_zx).subs(w,fw(u))).rhs
3442 z_y = dsolve(diff(fw(v),v) - ((b*F1-a*F2)/(a*F3-c*F1)).subs(u,x_yz).subs(w,fw(v))).rhs
3443 x_y = dsolve(diff(fu(v),v) - ((c*F2-b*F3)/(a*F3-c*F1)).subs(w,z_xy).subs(u,fu(v))).rhs
3444 y_z = dsolve(diff(fv(w),w) - ((a*F3-c*F1)/(b*F1-a*F2)).subs(u,x_yz).subs(v,fv(w))).rhs
3445 x_z = dsolve(diff(fu(w),w) - ((c*F2-b*F3)/(b*F1-a*F2)).subs(v,y_zx).subs(u,fu(w))).rhs
3446 sol1 = dsolve(diff(fu(t),t) - (c*F2 - b*F3).subs(v,y_x).subs(w,z_x).subs(u,fu(t))).rhs
3447 sol2 = dsolve(diff(fv(t),t) - (a*F3 - c*F1).subs(u,x_y).subs(w,z_y).subs(v,fv(t))).rhs
3448 sol3 = dsolve(diff(fw(t),t) - (b*F1 - a*F2).subs(u,x_z).subs(v,y_z).subs(w,fw(t))).rhs
3449 return [sol1, sol2, sol3]
3451def _nonlinear_3eq_order1_type4(x, y, z, t, eq):
3452 r"""
3453 Equations:
3455 .. math:: x' = c z F_2 - b y F_3, \enspace y' = a x F_3 - c z F_1, \enspace z' = b y F_1 - a x F_2
3457 where `F_n = F_n (x, y, z, t)`
3459 1. First integral:
3461 .. math:: a x^{2} + b y^{2} + c z^{2} = C_1
3463 where `C` is an arbitrary constant.
3465 2. Assuming the function `F_n` is independent of `t`: `F_n = F_n (x, y, z)`. Then on
3466 eliminating `t` and `z` from the first two equations of the system, one arrives at
3467 the first-order equation
3469 .. math:: \frac{dy}{dx} = \frac{a x F_3 (x, y, z) - c z F_1 (x, y, z)}
3470 {c z F_2 (x, y, z) - b y F_3 (x, y, z)}
3472 where `z = \pm \sqrt{\frac{1}{c} (C_1 - a x^{2} - b y^{2})}`
3474 References
3475 ==========
3476 -https://eqworld.ipmnet.ru/en/solutions/sysode/sode0405.pdf
3478 """
3479 C1 = get_numbered_constants(eq, num=1)
3480 u, v, w = symbols('u, v, w')
3481 p = Wild('p', exclude=[x(t), y(t), z(t), t])
3482 q = Wild('q', exclude=[x(t), y(t), z(t), t])
3483 s = Wild('s', exclude=[x(t), y(t), z(t), t])
3484 F1, F2, F3 = symbols('F1, F2, F3', cls=Wild)
3485 r1 = eq[0].match(diff(x(t),t) - z(t)*F2 + y(t)*F3)
3486 r = collect_const(r1[F2]).match(s*F2)
3487 r.update(collect_const(r1[F3]).match(q*F3))
3488 if eq[1].has(r[F2]) and not eq[1].has(r[F3]):
3489 r[F2], r[F3] = r[F3], r[F2]
3490 r[s], r[q] = -r[q], -r[s]
3491 r.update((diff(y(t),t) - eq[1]).match(p*x(t)*r[F3] - r[s]*z(t)*F1))
3492 a = r[p]; b = r[q]; c = r[s]
3493 F1 = r[F1].subs(x(t),u).subs(y(t),v).subs(z(t),w)
3494 F2 = r[F2].subs(x(t),u).subs(y(t),v).subs(z(t),w)
3495 F3 = r[F3].subs(x(t),u).subs(y(t),v).subs(z(t),w)
3496 x_yz = sqrt((C1 - b*v**2 - c*w**2)/a)
3497 y_zx = sqrt((C1 - c*w**2 - a*u**2)/b)
3498 z_xy = sqrt((C1 - a*u**2 - b*v**2)/c)
3499 y_x = dsolve(diff(v(u),u) - ((a*u*F3-c*w*F1)/(c*w*F2-b*v*F3)).subs(w,z_xy).subs(v,v(u))).rhs
3500 z_x = dsolve(diff(w(u),u) - ((b*v*F1-a*u*F2)/(c*w*F2-b*v*F3)).subs(v,y_zx).subs(w,w(u))).rhs
3501 z_y = dsolve(diff(w(v),v) - ((b*v*F1-a*u*F2)/(a*u*F3-c*w*F1)).subs(u,x_yz).subs(w,w(v))).rhs
3502 x_y = dsolve(diff(u(v),v) - ((c*w*F2-b*v*F3)/(a*u*F3-c*w*F1)).subs(w,z_xy).subs(u,u(v))).rhs
3503 y_z = dsolve(diff(v(w),w) - ((a*u*F3-c*w*F1)/(b*v*F1-a*u*F2)).subs(u,x_yz).subs(v,v(w))).rhs
3504 x_z = dsolve(diff(u(w),w) - ((c*w*F2-b*v*F3)/(b*v*F1-a*u*F2)).subs(v,y_zx).subs(u,u(w))).rhs
3505 sol1 = dsolve(diff(u(t),t) - (c*w*F2 - b*v*F3).subs(v,y_x).subs(w,z_x).subs(u,u(t))).rhs
3506 sol2 = dsolve(diff(v(t),t) - (a*u*F3 - c*w*F1).subs(u,x_y).subs(w,z_y).subs(v,v(t))).rhs
3507 sol3 = dsolve(diff(w(t),t) - (b*v*F1 - a*u*F2).subs(u,x_z).subs(v,y_z).subs(w,w(t))).rhs
3508 return [sol1, sol2, sol3]
3510def _nonlinear_3eq_order1_type5(x, y, z, t, eq):
3511 r"""
3512 .. math:: x' = x (c F_2 - b F_3), \enspace y' = y (a F_3 - c F_1), \enspace z' = z (b F_1 - a F_2)
3514 where `F_n = F_n (x, y, z, t)` and are arbitrary functions.
3516 First Integral:
3518 .. math:: \left|x\right|^{a} \left|y\right|^{b} \left|z\right|^{c} = C_1
3520 where `C` is an arbitrary constant. If the function `F_n` is independent of `t`,
3521 then, by eliminating `t` and `z` from the first two equations of the system, one
3522 arrives at a first-order equation.
3524 References
3525 ==========
3526 -https://eqworld.ipmnet.ru/en/solutions/sysode/sode0406.pdf
3528 """
3529 C1 = get_numbered_constants(eq, num=1)
3530 u, v, w = symbols('u, v, w')
3531 fu, fv, fw = symbols('u, v, w', cls=Function)
3532 p = Wild('p', exclude=[x(t), y(t), z(t), t])
3533 q = Wild('q', exclude=[x(t), y(t), z(t), t])
3534 s = Wild('s', exclude=[x(t), y(t), z(t), t])
3535 F1, F2, F3 = symbols('F1, F2, F3', cls=Wild)
3536 r1 = eq[0].match(diff(x(t), t) - x(t)*F2 + x(t)*F3)
3537 r = collect_const(r1[F2]).match(s*F2)
3538 r.update(collect_const(r1[F3]).match(q*F3))
3539 if eq[1].has(r[F2]) and not eq[1].has(r[F3]):
3540 r[F2], r[F3] = r[F3], r[F2]
3541 r[s], r[q] = -r[q], -r[s]
3542 r.update((diff(y(t), t) - eq[1]).match(y(t)*(p*r[F3] - r[s]*F1)))
3543 a = r[p]; b = r[q]; c = r[s]
3544 F1 = r[F1].subs(x(t), u).subs(y(t), v).subs(z(t), w)
3545 F2 = r[F2].subs(x(t), u).subs(y(t), v).subs(z(t), w)
3546 F3 = r[F3].subs(x(t), u).subs(y(t), v).subs(z(t), w)
3547 x_yz = (C1*v**-b*w**-c)**-a
3548 y_zx = (C1*w**-c*u**-a)**-b
3549 z_xy = (C1*u**-a*v**-b)**-c
3550 y_x = dsolve(diff(fv(u), u) - ((v*(a*F3 - c*F1))/(u*(c*F2 - b*F3))).subs(w, z_xy).subs(v, fv(u))).rhs
3551 z_x = dsolve(diff(fw(u), u) - ((w*(b*F1 - a*F2))/(u*(c*F2 - b*F3))).subs(v, y_zx).subs(w, fw(u))).rhs
3552 z_y = dsolve(diff(fw(v), v) - ((w*(b*F1 - a*F2))/(v*(a*F3 - c*F1))).subs(u, x_yz).subs(w, fw(v))).rhs
3553 x_y = dsolve(diff(fu(v), v) - ((u*(c*F2 - b*F3))/(v*(a*F3 - c*F1))).subs(w, z_xy).subs(u, fu(v))).rhs
3554 y_z = dsolve(diff(fv(w), w) - ((v*(a*F3 - c*F1))/(w*(b*F1 - a*F2))).subs(u, x_yz).subs(v, fv(w))).rhs
3555 x_z = dsolve(diff(fu(w), w) - ((u*(c*F2 - b*F3))/(w*(b*F1 - a*F2))).subs(v, y_zx).subs(u, fu(w))).rhs
3556 sol1 = dsolve(diff(fu(t), t) - (u*(c*F2 - b*F3)).subs(v, y_x).subs(w, z_x).subs(u, fu(t))).rhs
3557 sol2 = dsolve(diff(fv(t), t) - (v*(a*F3 - c*F1)).subs(u, x_y).subs(w, z_y).subs(v, fv(t))).rhs
3558 sol3 = dsolve(diff(fw(t), t) - (w*(b*F1 - a*F2)).subs(u, x_z).subs(v, y_z).subs(w, fw(t))).rhs
3559 return [sol1, sol2, sol3]
3562#This import is written at the bottom to avoid circular imports.
3563from .single import SingleODEProblem, SingleODESolver, solver_map