Coverage for /usr/lib/python3/dist-packages/sympy/solvers/ode/single.py: 23%
1114 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-14 15:55 +0200
1#
2# This is the module for ODE solver classes for single ODEs.
3#
5from __future__ import annotations
6from typing import ClassVar, Iterator
8from .riccati import match_riccati, solve_riccati
9from sympy.core import Add, S, Pow, Rational
10from sympy.core.cache import cached_property
11from sympy.core.exprtools import factor_terms
12from sympy.core.expr import Expr
13from sympy.core.function import AppliedUndef, Derivative, diff, Function, expand, Subs, _mexpand
14from sympy.core.numbers import zoo
15from sympy.core.relational import Equality, Eq
16from sympy.core.symbol import Symbol, Dummy, Wild
17from sympy.core.mul import Mul
18from sympy.functions import exp, tan, log, sqrt, besselj, bessely, cbrt, airyai, airybi
19from sympy.integrals import Integral
20from sympy.polys import Poly
21from sympy.polys.polytools import cancel, factor, degree
22from sympy.simplify import collect, simplify, separatevars, logcombine, posify # type: ignore
23from sympy.simplify.radsimp import fraction
24from sympy.utilities import numbered_symbols
25from sympy.solvers.solvers import solve
26from sympy.solvers.deutils import ode_order, _preprocess
27from sympy.polys.matrices.linsolve import _lin_eq2dict
28from sympy.polys.solvers import PolyNonlinearError
29from .hypergeometric import equivalence_hypergeometric, match_2nd_2F1_hypergeometric, \
30 get_sol_2F1_hypergeometric, match_2nd_hypergeometric
31from .nonhomogeneous import _get_euler_characteristic_eq_sols, _get_const_characteristic_eq_sols, \
32 _solve_undetermined_coefficients, _solve_variation_of_parameters, _test_term, _undetermined_coefficients_match, \
33 _get_simplified_sol
34from .lie_group import _ode_lie_group
37class ODEMatchError(NotImplementedError):
38 """Raised if a SingleODESolver is asked to solve an ODE it does not match"""
39 pass
42class SingleODEProblem:
43 """Represents an ordinary differential equation (ODE)
45 This class is used internally in the by dsolve and related
46 functions/classes so that properties of an ODE can be computed
47 efficiently.
49 Examples
50 ========
52 This class is used internally by dsolve. To instantiate an instance
53 directly first define an ODE problem:
55 >>> from sympy import Function, Symbol
56 >>> x = Symbol('x')
57 >>> f = Function('f')
58 >>> eq = f(x).diff(x, 2)
60 Now you can create a SingleODEProblem instance and query its properties:
62 >>> from sympy.solvers.ode.single import SingleODEProblem
63 >>> problem = SingleODEProblem(f(x).diff(x), f(x), x)
64 >>> problem.eq
65 Derivative(f(x), x)
66 >>> problem.func
67 f(x)
68 >>> problem.sym
69 x
70 """
72 # Instance attributes:
73 eq = None # type: Expr
74 func = None # type: AppliedUndef
75 sym = None # type: Symbol
76 _order = None # type: int
77 _eq_expanded = None # type: Expr
78 _eq_preprocessed = None # type: Expr
79 _eq_high_order_free = None
81 def __init__(self, eq, func, sym, prep=True, **kwargs):
82 assert isinstance(eq, Expr)
83 assert isinstance(func, AppliedUndef)
84 assert isinstance(sym, Symbol)
85 assert isinstance(prep, bool)
86 self.eq = eq
87 self.func = func
88 self.sym = sym
89 self.prep = prep
90 self.params = kwargs
92 @cached_property
93 def order(self) -> int:
94 return ode_order(self.eq, self.func)
96 @cached_property
97 def eq_preprocessed(self) -> Expr:
98 return self._get_eq_preprocessed()
100 @cached_property
101 def eq_high_order_free(self) -> Expr:
102 a = Wild('a', exclude=[self.func])
103 c1 = Wild('c1', exclude=[self.sym])
104 # Precondition to try remove f(x) from highest order derivative
105 reduced_eq = None
106 if self.eq.is_Add:
107 deriv_coef = self.eq.coeff(self.func.diff(self.sym, self.order))
108 if deriv_coef not in (1, 0):
109 r = deriv_coef.match(a*self.func**c1)
110 if r and r[c1]:
111 den = self.func**r[c1]
112 reduced_eq = Add(*[arg/den for arg in self.eq.args])
113 if not reduced_eq:
114 reduced_eq = expand(self.eq)
115 return reduced_eq
117 @cached_property
118 def eq_expanded(self) -> Expr:
119 return expand(self.eq_preprocessed)
121 def _get_eq_preprocessed(self) -> Expr:
122 if self.prep:
123 process_eq, process_func = _preprocess(self.eq, self.func)
124 if process_func != self.func:
125 raise ValueError
126 else:
127 process_eq = self.eq
128 return process_eq
130 def get_numbered_constants(self, num=1, start=1, prefix='C') -> list[Symbol]:
131 """
132 Returns a list of constants that do not occur
133 in eq already.
134 """
135 ncs = self.iter_numbered_constants(start, prefix)
136 Cs = [next(ncs) for i in range(num)]
137 return Cs
139 def iter_numbered_constants(self, start=1, prefix='C') -> Iterator[Symbol]:
140 """
141 Returns an iterator of constants that do not occur
142 in eq already.
143 """
144 atom_set = self.eq.free_symbols
145 func_set = self.eq.atoms(Function)
146 if func_set:
147 atom_set |= {Symbol(str(f.func)) for f in func_set}
148 return numbered_symbols(start=start, prefix=prefix, exclude=atom_set)
150 @cached_property
151 def is_autonomous(self):
152 u = Dummy('u')
153 x = self.sym
154 syms = self.eq.subs(self.func, u).free_symbols
155 return x not in syms
157 def get_linear_coefficients(self, eq, func, order):
158 r"""
159 Matches a differential equation to the linear form:
161 .. math:: a_n(x) y^{(n)} + \cdots + a_1(x)y' + a_0(x) y + B(x) = 0
163 Returns a dict of order:coeff terms, where order is the order of the
164 derivative on each term, and coeff is the coefficient of that derivative.
165 The key ``-1`` holds the function `B(x)`. Returns ``None`` if the ODE is
166 not linear. This function assumes that ``func`` has already been checked
167 to be good.
169 Examples
170 ========
172 >>> from sympy import Function, cos, sin
173 >>> from sympy.abc import x
174 >>> from sympy.solvers.ode.single import SingleODEProblem
175 >>> f = Function('f')
176 >>> eq = f(x).diff(x, 3) + 2*f(x).diff(x) + \
177 ... x*f(x).diff(x, 2) + cos(x)*f(x).diff(x) + x - f(x) - \
178 ... sin(x)
179 >>> obj = SingleODEProblem(eq, f(x), x)
180 >>> obj.get_linear_coefficients(eq, f(x), 3)
181 {-1: x - sin(x), 0: -1, 1: cos(x) + 2, 2: x, 3: 1}
182 >>> eq = f(x).diff(x, 3) + 2*f(x).diff(x) + \
183 ... x*f(x).diff(x, 2) + cos(x)*f(x).diff(x) + x - f(x) - \
184 ... sin(f(x))
185 >>> obj = SingleODEProblem(eq, f(x), x)
186 >>> obj.get_linear_coefficients(eq, f(x), 3) == None
187 True
189 """
190 f = func.func
191 x = func.args[0]
192 symset = {Derivative(f(x), x, i) for i in range(order+1)}
193 try:
194 rhs, lhs_terms = _lin_eq2dict(eq, symset)
195 except PolyNonlinearError:
196 return None
198 if rhs.has(func) or any(c.has(func) for c in lhs_terms.values()):
199 return None
200 terms = {i: lhs_terms.get(f(x).diff(x, i), S.Zero) for i in range(order+1)}
201 terms[-1] = rhs
202 return terms
204 # TODO: Add methods that can be used by many ODE solvers:
205 # order
206 # is_linear()
207 # get_linear_coefficients()
208 # eq_prepared (the ODE in prepared form)
211class SingleODESolver:
212 """
213 Base class for Single ODE solvers.
215 Subclasses should implement the _matches and _get_general_solution
216 methods. This class is not intended to be instantiated directly but its
217 subclasses are as part of dsolve.
219 Examples
220 ========
222 You can use a subclass of SingleODEProblem to solve a particular type of
223 ODE. We first define a particular ODE problem:
225 >>> from sympy import Function, Symbol
226 >>> x = Symbol('x')
227 >>> f = Function('f')
228 >>> eq = f(x).diff(x, 2)
230 Now we solve this problem using the NthAlgebraic solver which is a
231 subclass of SingleODESolver:
233 >>> from sympy.solvers.ode.single import NthAlgebraic, SingleODEProblem
234 >>> problem = SingleODEProblem(eq, f(x), x)
235 >>> solver = NthAlgebraic(problem)
236 >>> solver.get_general_solution()
237 [Eq(f(x), _C*x + _C)]
239 The normal way to solve an ODE is to use dsolve (which would use
240 NthAlgebraic and other solvers internally). When using dsolve a number of
241 other things are done such as evaluating integrals, simplifying the
242 solution and renumbering the constants:
244 >>> from sympy import dsolve
245 >>> dsolve(eq, hint='nth_algebraic')
246 Eq(f(x), C1 + C2*x)
247 """
249 # Subclasses should store the hint name (the argument to dsolve) in this
250 # attribute
251 hint: ClassVar[str]
253 # Subclasses should define this to indicate if they support an _Integral
254 # hint.
255 has_integral: ClassVar[bool]
257 # The ODE to be solved
258 ode_problem = None # type: SingleODEProblem
260 # Cache whether or not the equation has matched the method
261 _matched: bool | None = None
263 # Subclasses should store in this attribute the list of order(s) of ODE
264 # that subclass can solve or leave it to None if not specific to any order
265 order: list | None = None
267 def __init__(self, ode_problem):
268 self.ode_problem = ode_problem
270 def matches(self) -> bool:
271 if self.order is not None and self.ode_problem.order not in self.order:
272 self._matched = False
273 return self._matched
275 if self._matched is None:
276 self._matched = self._matches()
277 return self._matched
279 def get_general_solution(self, *, simplify: bool = True) -> list[Equality]:
280 if not self.matches():
281 msg = "%s solver cannot solve:\n%s"
282 raise ODEMatchError(msg % (self.hint, self.ode_problem.eq))
283 return self._get_general_solution(simplify_flag=simplify)
285 def _matches(self) -> bool:
286 msg = "Subclasses of SingleODESolver should implement matches."
287 raise NotImplementedError(msg)
289 def _get_general_solution(self, *, simplify_flag: bool = True) -> list[Equality]:
290 msg = "Subclasses of SingleODESolver should implement get_general_solution."
291 raise NotImplementedError(msg)
294class SinglePatternODESolver(SingleODESolver):
295 '''Superclass for ODE solvers based on pattern matching'''
297 def wilds(self):
298 prob = self.ode_problem
299 f = prob.func.func
300 x = prob.sym
301 order = prob.order
302 return self._wilds(f, x, order)
304 def wilds_match(self):
305 match = self._wilds_match
306 return [match.get(w, S.Zero) for w in self.wilds()]
308 def _matches(self):
309 eq = self.ode_problem.eq_expanded
310 f = self.ode_problem.func.func
311 x = self.ode_problem.sym
312 order = self.ode_problem.order
313 df = f(x).diff(x, order)
315 if order not in [1, 2]:
316 return False
318 pattern = self._equation(f(x), x, order)
320 if not pattern.coeff(df).has(Wild):
321 eq = expand(eq / eq.coeff(df))
322 eq = eq.collect([f(x).diff(x), f(x)], func = cancel)
324 self._wilds_match = match = eq.match(pattern)
325 if match is not None:
326 return self._verify(f(x))
327 return False
329 def _verify(self, fx) -> bool:
330 return True
332 def _wilds(self, f, x, order):
333 msg = "Subclasses of SingleODESolver should implement _wilds"
334 raise NotImplementedError(msg)
336 def _equation(self, fx, x, order):
337 msg = "Subclasses of SingleODESolver should implement _equation"
338 raise NotImplementedError(msg)
341class NthAlgebraic(SingleODESolver):
342 r"""
343 Solves an `n`\th order ordinary differential equation using algebra and
344 integrals.
346 There is no general form for the kind of equation that this can solve. The
347 the equation is solved algebraically treating differentiation as an
348 invertible algebraic function.
350 Examples
351 ========
353 >>> from sympy import Function, dsolve, Eq
354 >>> from sympy.abc import x
355 >>> f = Function('f')
356 >>> eq = Eq(f(x) * (f(x).diff(x)**2 - 1), 0)
357 >>> dsolve(eq, f(x), hint='nth_algebraic')
358 [Eq(f(x), 0), Eq(f(x), C1 - x), Eq(f(x), C1 + x)]
360 Note that this solver can return algebraic solutions that do not have any
361 integration constants (f(x) = 0 in the above example).
362 """
364 hint = 'nth_algebraic'
365 has_integral = True # nth_algebraic_Integral hint
367 def _matches(self):
368 r"""
369 Matches any differential equation that nth_algebraic can solve. Uses
370 `sympy.solve` but teaches it how to integrate derivatives.
372 This involves calling `sympy.solve` and does most of the work of finding a
373 solution (apart from evaluating the integrals).
374 """
375 eq = self.ode_problem.eq
376 func = self.ode_problem.func
377 var = self.ode_problem.sym
379 # Derivative that solve can handle:
380 diffx = self._get_diffx(var)
382 # Replace derivatives wrt the independent variable with diffx
383 def replace(eq, var):
384 def expand_diffx(*args):
385 differand, diffs = args[0], args[1:]
386 toreplace = differand
387 for v, n in diffs:
388 for _ in range(n):
389 if v == var:
390 toreplace = diffx(toreplace)
391 else:
392 toreplace = Derivative(toreplace, v)
393 return toreplace
394 return eq.replace(Derivative, expand_diffx)
396 # Restore derivatives in solution afterwards
397 def unreplace(eq, var):
398 return eq.replace(diffx, lambda e: Derivative(e, var))
400 subs_eqn = replace(eq, var)
401 try:
402 # turn off simplification to protect Integrals that have
403 # _t instead of fx in them and would otherwise factor
404 # as t_*Integral(1, x)
405 solns = solve(subs_eqn, func, simplify=False)
406 except NotImplementedError:
407 solns = []
409 solns = [simplify(unreplace(soln, var)) for soln in solns]
410 solns = [Equality(func, soln) for soln in solns]
412 self.solutions = solns
413 return len(solns) != 0
415 def _get_general_solution(self, *, simplify_flag: bool = True):
416 return self.solutions
418 # This needs to produce an invertible function but the inverse depends
419 # which variable we are integrating with respect to. Since the class can
420 # be stored in cached results we need to ensure that we always get the
421 # same class back for each particular integration variable so we store these
422 # classes in a global dict:
423 _diffx_stored: dict[Symbol, type[Function]] = {}
425 @staticmethod
426 def _get_diffx(var):
427 diffcls = NthAlgebraic._diffx_stored.get(var, None)
429 if diffcls is None:
430 # A class that behaves like Derivative wrt var but is "invertible".
431 class diffx(Function):
432 def inverse(self):
433 # don't use integrate here because fx has been replaced by _t
434 # in the equation; integrals will not be correct while solve
435 # is at work.
436 return lambda expr: Integral(expr, var) + Dummy('C')
438 diffcls = NthAlgebraic._diffx_stored.setdefault(var, diffx)
440 return diffcls
443class FirstExact(SinglePatternODESolver):
444 r"""
445 Solves 1st order exact ordinary differential equations.
447 A 1st order differential equation is called exact if it is the total
448 differential of a function. That is, the differential equation
450 .. math:: P(x, y) \,\partial{}x + Q(x, y) \,\partial{}y = 0
452 is exact if there is some function `F(x, y)` such that `P(x, y) =
453 \partial{}F/\partial{}x` and `Q(x, y) = \partial{}F/\partial{}y`. It can
454 be shown that a necessary and sufficient condition for a first order ODE
455 to be exact is that `\partial{}P/\partial{}y = \partial{}Q/\partial{}x`.
456 Then, the solution will be as given below::
458 >>> from sympy import Function, Eq, Integral, symbols, pprint
459 >>> x, y, t, x0, y0, C1= symbols('x,y,t,x0,y0,C1')
460 >>> P, Q, F= map(Function, ['P', 'Q', 'F'])
461 >>> pprint(Eq(Eq(F(x, y), Integral(P(t, y), (t, x0, x)) +
462 ... Integral(Q(x0, t), (t, y0, y))), C1))
463 x y
464 / /
465 | |
466 F(x, y) = | P(t, y) dt + | Q(x0, t) dt = C1
467 | |
468 / /
469 x0 y0
471 Where the first partials of `P` and `Q` exist and are continuous in a
472 simply connected region.
474 A note: SymPy currently has no way to represent inert substitution on an
475 expression, so the hint ``1st_exact_Integral`` will return an integral
476 with `dy`. This is supposed to represent the function that you are
477 solving for.
479 Examples
480 ========
482 >>> from sympy import Function, dsolve, cos, sin
483 >>> from sympy.abc import x
484 >>> f = Function('f')
485 >>> dsolve(cos(f(x)) - (x*sin(f(x)) - f(x)**2)*f(x).diff(x),
486 ... f(x), hint='1st_exact')
487 Eq(x*cos(f(x)) + f(x)**3/3, C1)
489 References
490 ==========
492 - https://en.wikipedia.org/wiki/Exact_differential_equation
493 - M. Tenenbaum & H. Pollard, "Ordinary Differential Equations",
494 Dover 1963, pp. 73
496 # indirect doctest
498 """
499 hint = "1st_exact"
500 has_integral = True
501 order = [1]
503 def _wilds(self, f, x, order):
504 P = Wild('P', exclude=[f(x).diff(x)])
505 Q = Wild('Q', exclude=[f(x).diff(x)])
506 return P, Q
508 def _equation(self, fx, x, order):
509 P, Q = self.wilds()
510 return P + Q*fx.diff(x)
512 def _verify(self, fx) -> bool:
513 P, Q = self.wilds()
514 x = self.ode_problem.sym
515 y = Dummy('y')
517 m, n = self.wilds_match()
519 m = m.subs(fx, y)
520 n = n.subs(fx, y)
521 numerator = cancel(m.diff(y) - n.diff(x))
523 if numerator.is_zero:
524 # Is exact
525 return True
526 else:
527 # The following few conditions try to convert a non-exact
528 # differential equation into an exact one.
529 # References:
530 # 1. Differential equations with applications
531 # and historical notes - George E. Simmons
532 # 2. https://math.okstate.edu/people/binegar/2233-S99/2233-l12.pdf
534 factor_n = cancel(numerator/n)
535 factor_m = cancel(-numerator/m)
536 if y not in factor_n.free_symbols:
537 # If (dP/dy - dQ/dx) / Q = f(x)
538 # then exp(integral(f(x))*equation becomes exact
539 factor = factor_n
540 integration_variable = x
541 elif x not in factor_m.free_symbols:
542 # If (dP/dy - dQ/dx) / -P = f(y)
543 # then exp(integral(f(y))*equation becomes exact
544 factor = factor_m
545 integration_variable = y
546 else:
547 # Couldn't convert to exact
548 return False
550 factor = exp(Integral(factor, integration_variable))
551 m *= factor
552 n *= factor
553 self._wilds_match[P] = m.subs(y, fx)
554 self._wilds_match[Q] = n.subs(y, fx)
555 return True
557 def _get_general_solution(self, *, simplify_flag: bool = True):
558 m, n = self.wilds_match()
559 fx = self.ode_problem.func
560 x = self.ode_problem.sym
561 (C1,) = self.ode_problem.get_numbered_constants(num=1)
562 y = Dummy('y')
564 m = m.subs(fx, y)
565 n = n.subs(fx, y)
567 gen_sol = Eq(Subs(Integral(m, x)
568 + Integral(n - Integral(m, x).diff(y), y), y, fx), C1)
569 return [gen_sol]
572class FirstLinear(SinglePatternODESolver):
573 r"""
574 Solves 1st order linear differential equations.
576 These are differential equations of the form
578 .. math:: dy/dx + P(x) y = Q(x)\text{.}
580 These kinds of differential equations can be solved in a general way. The
581 integrating factor `e^{\int P(x) \,dx}` will turn the equation into a
582 separable equation. The general solution is::
584 >>> from sympy import Function, dsolve, Eq, pprint, diff, sin
585 >>> from sympy.abc import x
586 >>> f, P, Q = map(Function, ['f', 'P', 'Q'])
587 >>> genform = Eq(f(x).diff(x) + P(x)*f(x), Q(x))
588 >>> pprint(genform)
589 d
590 P(x)*f(x) + --(f(x)) = Q(x)
591 dx
592 >>> pprint(dsolve(genform, f(x), hint='1st_linear_Integral'))
593 / / \
594 | | |
595 | | / | /
596 | | | | |
597 | | | P(x) dx | - | P(x) dx
598 | | | | |
599 | | / | /
600 f(x) = |C1 + | Q(x)*e dx|*e
601 | | |
602 \ / /
605 Examples
606 ========
608 >>> f = Function('f')
609 >>> pprint(dsolve(Eq(x*diff(f(x), x) - f(x), x**2*sin(x)),
610 ... f(x), '1st_linear'))
611 f(x) = x*(C1 - cos(x))
613 References
614 ==========
616 - https://en.wikipedia.org/wiki/Linear_differential_equation#First-order_equation_with_variable_coefficients
617 - M. Tenenbaum & H. Pollard, "Ordinary Differential Equations",
618 Dover 1963, pp. 92
620 # indirect doctest
622 """
623 hint = '1st_linear'
624 has_integral = True
625 order = [1]
627 def _wilds(self, f, x, order):
628 P = Wild('P', exclude=[f(x)])
629 Q = Wild('Q', exclude=[f(x), f(x).diff(x)])
630 return P, Q
632 def _equation(self, fx, x, order):
633 P, Q = self.wilds()
634 return fx.diff(x) + P*fx - Q
636 def _get_general_solution(self, *, simplify_flag: bool = True):
637 P, Q = self.wilds_match()
638 fx = self.ode_problem.func
639 x = self.ode_problem.sym
640 (C1,) = self.ode_problem.get_numbered_constants(num=1)
641 gensol = Eq(fx, ((C1 + Integral(Q*exp(Integral(P, x)), x))
642 * exp(-Integral(P, x))))
643 return [gensol]
646class AlmostLinear(SinglePatternODESolver):
647 r"""
648 Solves an almost-linear differential equation.
650 The general form of an almost linear differential equation is
652 .. math:: a(x) g'(f(x)) f'(x) + b(x) g(f(x)) + c(x)
654 Here `f(x)` is the function to be solved for (the dependent variable).
655 The substitution `g(f(x)) = u(x)` leads to a linear differential equation
656 for `u(x)` of the form `a(x) u' + b(x) u + c(x) = 0`. This can be solved
657 for `u(x)` by the `first_linear` hint and then `f(x)` is found by solving
658 `g(f(x)) = u(x)`.
660 See Also
661 ========
662 :obj:`sympy.solvers.ode.single.FirstLinear`
664 Examples
665 ========
667 >>> from sympy import dsolve, Function, pprint, sin, cos
668 >>> from sympy.abc import x
669 >>> f = Function('f')
670 >>> d = f(x).diff(x)
671 >>> eq = x*d + x*f(x) + 1
672 >>> dsolve(eq, f(x), hint='almost_linear')
673 Eq(f(x), (C1 - Ei(x))*exp(-x))
674 >>> pprint(dsolve(eq, f(x), hint='almost_linear'))
675 -x
676 f(x) = (C1 - Ei(x))*e
677 >>> example = cos(f(x))*f(x).diff(x) + sin(f(x)) + 1
678 >>> pprint(example)
679 d
680 sin(f(x)) + cos(f(x))*--(f(x)) + 1
681 dx
682 >>> pprint(dsolve(example, f(x), hint='almost_linear'))
683 / -x \ / -x \
684 [f(x) = pi - asin\C1*e - 1/, f(x) = asin\C1*e - 1/]
687 References
688 ==========
690 - Joel Moses, "Symbolic Integration - The Stormy Decade", Communications
691 of the ACM, Volume 14, Number 8, August 1971, pp. 558
692 """
693 hint = "almost_linear"
694 has_integral = True
695 order = [1]
697 def _wilds(self, f, x, order):
698 P = Wild('P', exclude=[f(x).diff(x)])
699 Q = Wild('Q', exclude=[f(x).diff(x)])
700 return P, Q
702 def _equation(self, fx, x, order):
703 P, Q = self.wilds()
704 return P*fx.diff(x) + Q
706 def _verify(self, fx):
707 a, b = self.wilds_match()
708 c, b = b.as_independent(fx) if b.is_Add else (S.Zero, b)
709 # a, b and c are the function a(x), b(x) and c(x) respectively.
710 # c(x) is obtained by separating out b as terms with and without fx i.e, l(y)
711 # The following conditions checks if the given equation is an almost-linear differential equation using the fact that
712 # a(x)*(l(y))' / l(y)' is independent of l(y)
714 if b.diff(fx) != 0 and not simplify(b.diff(fx)/a).has(fx):
715 self.ly = factor_terms(b).as_independent(fx, as_Add=False)[1] # Gives the term containing fx i.e., l(y)
716 self.ax = a / self.ly.diff(fx)
717 self.cx = -c # cx is taken as -c(x) to simplify expression in the solution integral
718 self.bx = factor_terms(b) / self.ly
719 return True
721 return False
723 def _get_general_solution(self, *, simplify_flag: bool = True):
724 x = self.ode_problem.sym
725 (C1,) = self.ode_problem.get_numbered_constants(num=1)
726 gensol = Eq(self.ly, ((C1 + Integral((self.cx/self.ax)*exp(Integral(self.bx/self.ax, x)), x))
727 * exp(-Integral(self.bx/self.ax, x))))
729 return [gensol]
732class Bernoulli(SinglePatternODESolver):
733 r"""
734 Solves Bernoulli differential equations.
736 These are equations of the form
738 .. math:: dy/dx + P(x) y = Q(x) y^n\text{, }n \ne 1`\text{.}
740 The substitution `w = 1/y^{1-n}` will transform an equation of this form
741 into one that is linear (see the docstring of
742 :obj:`~sympy.solvers.ode.single.FirstLinear`). The general solution is::
744 >>> from sympy import Function, dsolve, Eq, pprint
745 >>> from sympy.abc import x, n
746 >>> f, P, Q = map(Function, ['f', 'P', 'Q'])
747 >>> genform = Eq(f(x).diff(x) + P(x)*f(x), Q(x)*f(x)**n)
748 >>> pprint(genform)
749 d n
750 P(x)*f(x) + --(f(x)) = Q(x)*f (x)
751 dx
752 >>> pprint(dsolve(genform, f(x), hint='Bernoulli_Integral'), num_columns=110)
753 -1
754 -----
755 n - 1
756 // / / \ \
757 || | | | |
758 || | / | / | / |
759 || | | | | | | |
760 || | -(n - 1)* | P(x) dx | -(n - 1)* | P(x) dx | (n - 1)* | P(x) dx|
761 || | | | | | | |
762 || | / | / | / |
763 f(x) = ||C1 - n* | Q(x)*e dx + | Q(x)*e dx|*e |
764 || | | | |
765 \\ / / / /
768 Note that the equation is separable when `n = 1` (see the docstring of
769 :obj:`~sympy.solvers.ode.single.Separable`).
771 >>> pprint(dsolve(Eq(f(x).diff(x) + P(x)*f(x), Q(x)*f(x)), f(x),
772 ... hint='separable_Integral'))
773 f(x)
774 /
775 | /
776 | 1 |
777 | - dy = C1 + | (-P(x) + Q(x)) dx
778 | y |
779 | /
780 /
783 Examples
784 ========
786 >>> from sympy import Function, dsolve, Eq, pprint, log
787 >>> from sympy.abc import x
788 >>> f = Function('f')
790 >>> pprint(dsolve(Eq(x*f(x).diff(x) + f(x), log(x)*f(x)**2),
791 ... f(x), hint='Bernoulli'))
792 1
793 f(x) = -----------------
794 C1*x + log(x) + 1
796 References
797 ==========
799 - https://en.wikipedia.org/wiki/Bernoulli_differential_equation
801 - M. Tenenbaum & H. Pollard, "Ordinary Differential Equations",
802 Dover 1963, pp. 95
804 # indirect doctest
806 """
807 hint = "Bernoulli"
808 has_integral = True
809 order = [1]
811 def _wilds(self, f, x, order):
812 P = Wild('P', exclude=[f(x)])
813 Q = Wild('Q', exclude=[f(x)])
814 n = Wild('n', exclude=[x, f(x), f(x).diff(x)])
815 return P, Q, n
817 def _equation(self, fx, x, order):
818 P, Q, n = self.wilds()
819 return fx.diff(x) + P*fx - Q*fx**n
821 def _get_general_solution(self, *, simplify_flag: bool = True):
822 P, Q, n = self.wilds_match()
823 fx = self.ode_problem.func
824 x = self.ode_problem.sym
825 (C1,) = self.ode_problem.get_numbered_constants(num=1)
826 if n==1:
827 gensol = Eq(log(fx), (
828 C1 + Integral((-P + Q), x)
829 ))
830 else:
831 gensol = Eq(fx**(1-n), (
832 (C1 - (n - 1) * Integral(Q*exp(-n*Integral(P, x))
833 * exp(Integral(P, x)), x)
834 ) * exp(-(1 - n)*Integral(P, x)))
835 )
836 return [gensol]
839class Factorable(SingleODESolver):
840 r"""
841 Solves equations having a solvable factor.
843 This function is used to solve the equation having factors. Factors may be of type algebraic or ode. It
844 will try to solve each factor independently. Factors will be solved by calling dsolve. We will return the
845 list of solutions.
847 Examples
848 ========
850 >>> from sympy import Function, dsolve, pprint
851 >>> from sympy.abc import x
852 >>> f = Function('f')
853 >>> eq = (f(x)**2-4)*(f(x).diff(x)+f(x))
854 >>> pprint(dsolve(eq, f(x)))
855 -x
856 [f(x) = 2, f(x) = -2, f(x) = C1*e ]
859 """
860 hint = "factorable"
861 has_integral = False
863 def _matches(self):
864 eq_orig = self.ode_problem.eq
865 f = self.ode_problem.func.func
866 x = self.ode_problem.sym
867 df = f(x).diff(x)
868 self.eqs = []
869 eq = eq_orig.collect(f(x), func = cancel)
870 eq = fraction(factor(eq))[0]
871 factors = Mul.make_args(factor(eq))
872 roots = [fac.as_base_exp() for fac in factors if len(fac.args)!=0]
873 if len(roots)>1 or roots[0][1]>1:
874 for base, expo in roots:
875 if base.has(f(x)):
876 self.eqs.append(base)
877 if len(self.eqs)>0:
878 return True
879 roots = solve(eq, df)
880 if len(roots)>0:
881 self.eqs = [(df - root) for root in roots]
882 # Avoid infinite recursion
883 matches = self.eqs != [eq_orig]
884 return matches
885 for i in factors:
886 if i.has(f(x)):
887 self.eqs.append(i)
888 return len(self.eqs)>0 and len(factors)>1
890 def _get_general_solution(self, *, simplify_flag: bool = True):
891 func = self.ode_problem.func.func
892 x = self.ode_problem.sym
893 eqns = self.eqs
894 sols = []
895 for eq in eqns:
896 try:
897 sol = dsolve(eq, func(x))
898 except NotImplementedError:
899 continue
900 else:
901 if isinstance(sol, list):
902 sols.extend(sol)
903 else:
904 sols.append(sol)
906 if sols == []:
907 raise NotImplementedError("The given ODE " + str(eq) + " cannot be solved by"
908 + " the factorable group method")
909 return sols
912class RiccatiSpecial(SinglePatternODESolver):
913 r"""
914 The general Riccati equation has the form
916 .. math:: dy/dx = f(x) y^2 + g(x) y + h(x)\text{.}
918 While it does not have a general solution [1], the "special" form, `dy/dx
919 = a y^2 - b x^c`, does have solutions in many cases [2]. This routine
920 returns a solution for `a(dy/dx) = b y^2 + c y/x + d/x^2` that is obtained
921 by using a suitable change of variables to reduce it to the special form
922 and is valid when neither `a` nor `b` are zero and either `c` or `d` is
923 zero.
925 >>> from sympy.abc import x, a, b, c, d
926 >>> from sympy import dsolve, checkodesol, pprint, Function
927 >>> f = Function('f')
928 >>> y = f(x)
929 >>> genform = a*y.diff(x) - (b*y**2 + c*y/x + d/x**2)
930 >>> sol = dsolve(genform, y, hint="Riccati_special_minus2")
931 >>> pprint(sol, wrap_line=False)
932 / / __________________ \\
933 | __________________ | / 2 ||
934 | / 2 | \/ 4*b*d - (a + c) *log(x)||
935 -|a + c - \/ 4*b*d - (a + c) *tan|C1 + ----------------------------||
936 \ \ 2*a //
937 f(x) = ------------------------------------------------------------------------
938 2*b*x
940 >>> checkodesol(genform, sol, order=1)[0]
941 True
943 References
944 ==========
946 - https://www.maplesoft.com/support/help/Maple/view.aspx?path=odeadvisor/Riccati
947 - https://eqworld.ipmnet.ru/en/solutions/ode/ode0106.pdf -
948 https://eqworld.ipmnet.ru/en/solutions/ode/ode0123.pdf
949 """
950 hint = "Riccati_special_minus2"
951 has_integral = False
952 order = [1]
954 def _wilds(self, f, x, order):
955 a = Wild('a', exclude=[x, f(x), f(x).diff(x), 0])
956 b = Wild('b', exclude=[x, f(x), f(x).diff(x), 0])
957 c = Wild('c', exclude=[x, f(x), f(x).diff(x)])
958 d = Wild('d', exclude=[x, f(x), f(x).diff(x)])
959 return a, b, c, d
961 def _equation(self, fx, x, order):
962 a, b, c, d = self.wilds()
963 return a*fx.diff(x) + b*fx**2 + c*fx/x + d/x**2
965 def _get_general_solution(self, *, simplify_flag: bool = True):
966 a, b, c, d = self.wilds_match()
967 fx = self.ode_problem.func
968 x = self.ode_problem.sym
969 (C1,) = self.ode_problem.get_numbered_constants(num=1)
970 mu = sqrt(4*d*b - (a - c)**2)
972 gensol = Eq(fx, (a - c - mu*tan(mu/(2*a)*log(x) + C1))/(2*b*x))
973 return [gensol]
976class RationalRiccati(SinglePatternODESolver):
977 r"""
978 Gives general solutions to the first order Riccati differential
979 equations that have atleast one rational particular solution.
981 .. math :: y' = b_0(x) + b_1(x) y + b_2(x) y^2
983 where `b_0`, `b_1` and `b_2` are rational functions of `x`
984 with `b_2 \ne 0` (`b_2 = 0` would make it a Bernoulli equation).
986 Examples
987 ========
989 >>> from sympy import Symbol, Function, dsolve, checkodesol
990 >>> f = Function('f')
991 >>> x = Symbol('x')
993 >>> eq = -x**4*f(x)**2 + x**3*f(x).diff(x) + x**2*f(x) + 20
994 >>> sol = dsolve(eq, hint="1st_rational_riccati")
995 >>> sol
996 Eq(f(x), (4*C1 - 5*x**9 - 4)/(x**2*(C1 + x**9 - 1)))
997 >>> checkodesol(eq, sol)
998 (True, 0)
1000 References
1001 ==========
1003 - Riccati ODE: https://en.wikipedia.org/wiki/Riccati_equation
1004 - N. Thieu Vo - Rational and Algebraic Solutions of First-Order Algebraic ODEs:
1005 Algorithm 11, pp. 78 - https://www3.risc.jku.at/publications/download/risc_5387/PhDThesisThieu.pdf
1006 """
1007 has_integral = False
1008 hint = "1st_rational_riccati"
1009 order = [1]
1011 def _wilds(self, f, x, order):
1012 b0 = Wild('b0', exclude=[f(x), f(x).diff(x)])
1013 b1 = Wild('b1', exclude=[f(x), f(x).diff(x)])
1014 b2 = Wild('b2', exclude=[f(x), f(x).diff(x)])
1015 return (b0, b1, b2)
1017 def _equation(self, fx, x, order):
1018 b0, b1, b2 = self.wilds()
1019 return fx.diff(x) - b0 - b1*fx - b2*fx**2
1021 def _matches(self):
1022 eq = self.ode_problem.eq_expanded
1023 f = self.ode_problem.func.func
1024 x = self.ode_problem.sym
1025 order = self.ode_problem.order
1027 if order != 1:
1028 return False
1030 match, funcs = match_riccati(eq, f, x)
1031 if not match:
1032 return False
1033 _b0, _b1, _b2 = funcs
1034 b0, b1, b2 = self.wilds()
1035 self._wilds_match = match = {b0: _b0, b1: _b1, b2: _b2}
1036 return True
1038 def _get_general_solution(self, *, simplify_flag: bool = True):
1039 # Match the equation
1040 b0, b1, b2 = self.wilds_match()
1041 fx = self.ode_problem.func
1042 x = self.ode_problem.sym
1043 return solve_riccati(fx, x, b0, b1, b2, gensol=True)
1046class SecondNonlinearAutonomousConserved(SinglePatternODESolver):
1047 r"""
1048 Gives solution for the autonomous second order nonlinear
1049 differential equation of the form
1051 .. math :: f''(x) = g(f(x))
1053 The solution for this differential equation can be computed
1054 by multiplying by `f'(x)` and integrating on both sides,
1055 converting it into a first order differential equation.
1057 Examples
1058 ========
1060 >>> from sympy import Function, symbols, dsolve
1061 >>> f, g = symbols('f g', cls=Function)
1062 >>> x = symbols('x')
1064 >>> eq = f(x).diff(x, 2) - g(f(x))
1065 >>> dsolve(eq, simplify=False)
1066 [Eq(Integral(1/sqrt(C1 + 2*Integral(g(_u), _u)), (_u, f(x))), C2 + x),
1067 Eq(Integral(1/sqrt(C1 + 2*Integral(g(_u), _u)), (_u, f(x))), C2 - x)]
1069 >>> from sympy import exp, log
1070 >>> eq = f(x).diff(x, 2) - exp(f(x)) + log(f(x))
1071 >>> dsolve(eq, simplify=False)
1072 [Eq(Integral(1/sqrt(-2*_u*log(_u) + 2*_u + C1 + 2*exp(_u)), (_u, f(x))), C2 + x),
1073 Eq(Integral(1/sqrt(-2*_u*log(_u) + 2*_u + C1 + 2*exp(_u)), (_u, f(x))), C2 - x)]
1075 References
1076 ==========
1078 - https://eqworld.ipmnet.ru/en/solutions/ode/ode0301.pdf
1079 """
1080 hint = "2nd_nonlinear_autonomous_conserved"
1081 has_integral = True
1082 order = [2]
1084 def _wilds(self, f, x, order):
1085 fy = Wild('fy', exclude=[0, f(x).diff(x), f(x).diff(x, 2)])
1086 return (fy, )
1088 def _equation(self, fx, x, order):
1089 fy = self.wilds()[0]
1090 return fx.diff(x, 2) + fy
1092 def _verify(self, fx):
1093 return self.ode_problem.is_autonomous
1095 def _get_general_solution(self, *, simplify_flag: bool = True):
1096 g = self.wilds_match()[0]
1097 fx = self.ode_problem.func
1098 x = self.ode_problem.sym
1099 u = Dummy('u')
1100 g = g.subs(fx, u)
1101 C1, C2 = self.ode_problem.get_numbered_constants(num=2)
1102 inside = -2*Integral(g, u) + C1
1103 lhs = Integral(1/sqrt(inside), (u, fx))
1104 return [Eq(lhs, C2 + x), Eq(lhs, C2 - x)]
1107class Liouville(SinglePatternODESolver):
1108 r"""
1109 Solves 2nd order Liouville differential equations.
1111 The general form of a Liouville ODE is
1113 .. math:: \frac{d^2 y}{dx^2} + g(y) \left(\!
1114 \frac{dy}{dx}\!\right)^2 + h(x)
1115 \frac{dy}{dx}\text{.}
1117 The general solution is:
1119 >>> from sympy import Function, dsolve, Eq, pprint, diff
1120 >>> from sympy.abc import x
1121 >>> f, g, h = map(Function, ['f', 'g', 'h'])
1122 >>> genform = Eq(diff(f(x),x,x) + g(f(x))*diff(f(x),x)**2 +
1123 ... h(x)*diff(f(x),x), 0)
1124 >>> pprint(genform)
1125 2 2
1126 /d \ d d
1127 g(f(x))*|--(f(x))| + h(x)*--(f(x)) + ---(f(x)) = 0
1128 \dx / dx 2
1129 dx
1130 >>> pprint(dsolve(genform, f(x), hint='Liouville_Integral'))
1131 f(x)
1132 / /
1133 | |
1134 | / | /
1135 | | | |
1136 | - | h(x) dx | | g(y) dy
1137 | | | |
1138 | / | /
1139 C1 + C2* | e dx + | e dy = 0
1140 | |
1141 / /
1143 Examples
1144 ========
1146 >>> from sympy import Function, dsolve, Eq, pprint
1147 >>> from sympy.abc import x
1148 >>> f = Function('f')
1149 >>> pprint(dsolve(diff(f(x), x, x) + diff(f(x), x)**2/f(x) +
1150 ... diff(f(x), x)/x, f(x), hint='Liouville'))
1151 ________________ ________________
1152 [f(x) = -\/ C1 + C2*log(x) , f(x) = \/ C1 + C2*log(x) ]
1154 References
1155 ==========
1157 - Goldstein and Braun, "Advanced Methods for the Solution of Differential
1158 Equations", pp. 98
1159 - https://www.maplesoft.com/support/help/Maple/view.aspx?path=odeadvisor/Liouville
1161 # indirect doctest
1163 """
1164 hint = "Liouville"
1165 has_integral = True
1166 order = [2]
1168 def _wilds(self, f, x, order):
1169 d = Wild('d', exclude=[f(x).diff(x), f(x).diff(x, 2)])
1170 e = Wild('e', exclude=[f(x).diff(x)])
1171 k = Wild('k', exclude=[f(x).diff(x)])
1172 return d, e, k
1174 def _equation(self, fx, x, order):
1175 # Liouville ODE in the form
1176 # f(x).diff(x, 2) + g(f(x))*(f(x).diff(x))**2 + h(x)*f(x).diff(x)
1177 # See Goldstein and Braun, "Advanced Methods for the Solution of
1178 # Differential Equations", pg. 98
1179 d, e, k = self.wilds()
1180 return d*fx.diff(x, 2) + e*fx.diff(x)**2 + k*fx.diff(x)
1182 def _verify(self, fx):
1183 d, e, k = self.wilds_match()
1184 self.y = Dummy('y')
1185 x = self.ode_problem.sym
1186 self.g = simplify(e/d).subs(fx, self.y)
1187 self.h = simplify(k/d).subs(fx, self.y)
1188 if self.y in self.h.free_symbols or x in self.g.free_symbols:
1189 return False
1190 return True
1192 def _get_general_solution(self, *, simplify_flag: bool = True):
1193 d, e, k = self.wilds_match()
1194 fx = self.ode_problem.func
1195 x = self.ode_problem.sym
1196 C1, C2 = self.ode_problem.get_numbered_constants(num=2)
1197 int = Integral(exp(Integral(self.g, self.y)), (self.y, None, fx))
1198 gen_sol = Eq(int + C1*Integral(exp(-Integral(self.h, x)), x) + C2, 0)
1200 return [gen_sol]
1203class Separable(SinglePatternODESolver):
1204 r"""
1205 Solves separable 1st order differential equations.
1207 This is any differential equation that can be written as `P(y)
1208 \tfrac{dy}{dx} = Q(x)`. The solution can then just be found by
1209 rearranging terms and integrating: `\int P(y) \,dy = \int Q(x) \,dx`.
1210 This hint uses :py:meth:`sympy.simplify.simplify.separatevars` as its back
1211 end, so if a separable equation is not caught by this solver, it is most
1212 likely the fault of that function.
1213 :py:meth:`~sympy.simplify.simplify.separatevars` is
1214 smart enough to do most expansion and factoring necessary to convert a
1215 separable equation `F(x, y)` into the proper form `P(x)\cdot{}Q(y)`. The
1216 general solution is::
1218 >>> from sympy import Function, dsolve, Eq, pprint
1219 >>> from sympy.abc import x
1220 >>> a, b, c, d, f = map(Function, ['a', 'b', 'c', 'd', 'f'])
1221 >>> genform = Eq(a(x)*b(f(x))*f(x).diff(x), c(x)*d(f(x)))
1222 >>> pprint(genform)
1223 d
1224 a(x)*b(f(x))*--(f(x)) = c(x)*d(f(x))
1225 dx
1226 >>> pprint(dsolve(genform, f(x), hint='separable_Integral'))
1227 f(x)
1228 / /
1229 | |
1230 | b(y) | c(x)
1231 | ---- dy = C1 + | ---- dx
1232 | d(y) | a(x)
1233 | |
1234 / /
1236 Examples
1237 ========
1239 >>> from sympy import Function, dsolve, Eq
1240 >>> from sympy.abc import x
1241 >>> f = Function('f')
1242 >>> pprint(dsolve(Eq(f(x)*f(x).diff(x) + x, 3*x*f(x)**2), f(x),
1243 ... hint='separable', simplify=False))
1244 / 2 \ 2
1245 log\3*f (x) - 1/ x
1246 ---------------- = C1 + --
1247 6 2
1249 References
1250 ==========
1252 - M. Tenenbaum & H. Pollard, "Ordinary Differential Equations",
1253 Dover 1963, pp. 52
1255 # indirect doctest
1257 """
1258 hint = "separable"
1259 has_integral = True
1260 order = [1]
1262 def _wilds(self, f, x, order):
1263 d = Wild('d', exclude=[f(x).diff(x), f(x).diff(x, 2)])
1264 e = Wild('e', exclude=[f(x).diff(x)])
1265 return d, e
1267 def _equation(self, fx, x, order):
1268 d, e = self.wilds()
1269 return d + e*fx.diff(x)
1271 def _verify(self, fx):
1272 d, e = self.wilds_match()
1273 self.y = Dummy('y')
1274 x = self.ode_problem.sym
1275 d = separatevars(d.subs(fx, self.y))
1276 e = separatevars(e.subs(fx, self.y))
1277 # m1[coeff]*m1[x]*m1[y] + m2[coeff]*m2[x]*m2[y]*y'
1278 self.m1 = separatevars(d, dict=True, symbols=(x, self.y))
1279 self.m2 = separatevars(e, dict=True, symbols=(x, self.y))
1280 if self.m1 and self.m2:
1281 return True
1282 return False
1284 def _get_match_object(self):
1285 fx = self.ode_problem.func
1286 x = self.ode_problem.sym
1287 return self.m1, self.m2, x, fx
1289 def _get_general_solution(self, *, simplify_flag: bool = True):
1290 m1, m2, x, fx = self._get_match_object()
1291 (C1,) = self.ode_problem.get_numbered_constants(num=1)
1292 int = Integral(m2['coeff']*m2[self.y]/m1[self.y],
1293 (self.y, None, fx))
1294 gen_sol = Eq(int, Integral(-m1['coeff']*m1[x]/
1295 m2[x], x) + C1)
1296 return [gen_sol]
1299class SeparableReduced(Separable):
1300 r"""
1301 Solves a differential equation that can be reduced to the separable form.
1303 The general form of this equation is
1305 .. math:: y' + (y/x) H(x^n y) = 0\text{}.
1307 This can be solved by substituting `u(y) = x^n y`. The equation then
1308 reduces to the separable form `\frac{u'}{u (\mathrm{power} - H(u))} -
1309 \frac{1}{x} = 0`.
1311 The general solution is:
1313 >>> from sympy import Function, dsolve, pprint
1314 >>> from sympy.abc import x, n
1315 >>> f, g = map(Function, ['f', 'g'])
1316 >>> genform = f(x).diff(x) + (f(x)/x)*g(x**n*f(x))
1317 >>> pprint(genform)
1318 / n \
1319 d f(x)*g\x *f(x)/
1320 --(f(x)) + ---------------
1321 dx x
1322 >>> pprint(dsolve(genform, hint='separable_reduced'))
1323 n
1324 x *f(x)
1325 /
1326 |
1327 | 1
1328 | ------------ dy = C1 + log(x)
1329 | y*(n - g(y))
1330 |
1331 /
1333 See Also
1334 ========
1335 :obj:`sympy.solvers.ode.single.Separable`
1337 Examples
1338 ========
1340 >>> from sympy import dsolve, Function, pprint
1341 >>> from sympy.abc import x
1342 >>> f = Function('f')
1343 >>> d = f(x).diff(x)
1344 >>> eq = (x - x**2*f(x))*d - f(x)
1345 >>> dsolve(eq, hint='separable_reduced')
1346 [Eq(f(x), (1 - sqrt(C1*x**2 + 1))/x), Eq(f(x), (sqrt(C1*x**2 + 1) + 1)/x)]
1347 >>> pprint(dsolve(eq, hint='separable_reduced'))
1348 ___________ ___________
1349 / 2 / 2
1350 1 - \/ C1*x + 1 \/ C1*x + 1 + 1
1351 [f(x) = ------------------, f(x) = ------------------]
1352 x x
1354 References
1355 ==========
1357 - Joel Moses, "Symbolic Integration - The Stormy Decade", Communications
1358 of the ACM, Volume 14, Number 8, August 1971, pp. 558
1359 """
1360 hint = "separable_reduced"
1361 has_integral = True
1362 order = [1]
1364 def _degree(self, expr, x):
1365 # Made this function to calculate the degree of
1366 # x in an expression. If expr will be of form
1367 # x**p*y, (wheare p can be variables/rationals) then it
1368 # will return p.
1369 for val in expr:
1370 if val.has(x):
1371 if isinstance(val, Pow) and val.as_base_exp()[0] == x:
1372 return (val.as_base_exp()[1])
1373 elif val == x:
1374 return (val.as_base_exp()[1])
1375 else:
1376 return self._degree(val.args, x)
1377 return 0
1379 def _powers(self, expr):
1380 # this function will return all the different relative power of x w.r.t f(x).
1381 # expr = x**p * f(x)**q then it will return {p/q}.
1382 pows = set()
1383 fx = self.ode_problem.func
1384 x = self.ode_problem.sym
1385 self.y = Dummy('y')
1386 if isinstance(expr, Add):
1387 exprs = expr.atoms(Add)
1388 elif isinstance(expr, Mul):
1389 exprs = expr.atoms(Mul)
1390 elif isinstance(expr, Pow):
1391 exprs = expr.atoms(Pow)
1392 else:
1393 exprs = {expr}
1395 for arg in exprs:
1396 if arg.has(x):
1397 _, u = arg.as_independent(x, fx)
1398 pow = self._degree((u.subs(fx, self.y), ), x)/self._degree((u.subs(fx, self.y), ), self.y)
1399 pows.add(pow)
1400 return pows
1402 def _verify(self, fx):
1403 num, den = self.wilds_match()
1404 x = self.ode_problem.sym
1405 factor = simplify(x/fx*num/den)
1406 # Try representing factor in terms of x^n*y
1407 # where n is lowest power of x in factor;
1408 # first remove terms like sqrt(2)*3 from factor.atoms(Mul)
1409 num, dem = factor.as_numer_denom()
1410 num = expand(num)
1411 dem = expand(dem)
1412 pows = self._powers(num)
1413 pows.update(self._powers(dem))
1414 pows = list(pows)
1415 if(len(pows)==1) and pows[0]!=zoo:
1416 self.t = Dummy('t')
1417 self.r2 = {'t': self.t}
1418 num = num.subs(x**pows[0]*fx, self.t)
1419 dem = dem.subs(x**pows[0]*fx, self.t)
1420 test = num/dem
1421 free = test.free_symbols
1422 if len(free) == 1 and free.pop() == self.t:
1423 self.r2.update({'power' : pows[0], 'u' : test})
1424 return True
1425 return False
1426 return False
1428 def _get_match_object(self):
1429 fx = self.ode_problem.func
1430 x = self.ode_problem.sym
1431 u = self.r2['u'].subs(self.r2['t'], self.y)
1432 ycoeff = 1/(self.y*(self.r2['power'] - u))
1433 m1 = {self.y: 1, x: -1/x, 'coeff': 1}
1434 m2 = {self.y: ycoeff, x: 1, 'coeff': 1}
1435 return m1, m2, x, x**self.r2['power']*fx
1438class HomogeneousCoeffSubsDepDivIndep(SinglePatternODESolver):
1439 r"""
1440 Solves a 1st order differential equation with homogeneous coefficients
1441 using the substitution `u_1 = \frac{\text{<dependent
1442 variable>}}{\text{<independent variable>}}`.
1444 This is a differential equation
1446 .. math:: P(x, y) + Q(x, y) dy/dx = 0
1448 such that `P` and `Q` are homogeneous and of the same order. A function
1449 `F(x, y)` is homogeneous of order `n` if `F(x t, y t) = t^n F(x, y)`.
1450 Equivalently, `F(x, y)` can be rewritten as `G(y/x)` or `H(x/y)`. See
1451 also the docstring of :py:meth:`~sympy.solvers.ode.homogeneous_order`.
1453 If the coefficients `P` and `Q` in the differential equation above are
1454 homogeneous functions of the same order, then it can be shown that the
1455 substitution `y = u_1 x` (i.e. `u_1 = y/x`) will turn the differential
1456 equation into an equation separable in the variables `x` and `u`. If
1457 `h(u_1)` is the function that results from making the substitution `u_1 =
1458 f(x)/x` on `P(x, f(x))` and `g(u_2)` is the function that results from the
1459 substitution on `Q(x, f(x))` in the differential equation `P(x, f(x)) +
1460 Q(x, f(x)) f'(x) = 0`, then the general solution is::
1462 >>> from sympy import Function, dsolve, pprint
1463 >>> from sympy.abc import x
1464 >>> f, g, h = map(Function, ['f', 'g', 'h'])
1465 >>> genform = g(f(x)/x) + h(f(x)/x)*f(x).diff(x)
1466 >>> pprint(genform)
1467 /f(x)\ /f(x)\ d
1468 g|----| + h|----|*--(f(x))
1469 \ x / \ x / dx
1470 >>> pprint(dsolve(genform, f(x),
1471 ... hint='1st_homogeneous_coeff_subs_dep_div_indep_Integral'))
1472 f(x)
1473 ----
1474 x
1475 /
1476 |
1477 | -h(u1)
1478 log(x) = C1 + | ---------------- d(u1)
1479 | u1*h(u1) + g(u1)
1480 |
1481 /
1483 Where `u_1 h(u_1) + g(u_1) \ne 0` and `x \ne 0`.
1485 See also the docstrings of
1486 :obj:`~sympy.solvers.ode.single.HomogeneousCoeffBest` and
1487 :obj:`~sympy.solvers.ode.single.HomogeneousCoeffSubsIndepDivDep`.
1489 Examples
1490 ========
1492 >>> from sympy import Function, dsolve
1493 >>> from sympy.abc import x
1494 >>> f = Function('f')
1495 >>> pprint(dsolve(2*x*f(x) + (x**2 + f(x)**2)*f(x).diff(x), f(x),
1496 ... hint='1st_homogeneous_coeff_subs_dep_div_indep', simplify=False))
1497 / 3 \
1498 |3*f(x) f (x)|
1499 log|------ + -----|
1500 | x 3 |
1501 \ x /
1502 log(x) = log(C1) - -------------------
1503 3
1505 References
1506 ==========
1508 - https://en.wikipedia.org/wiki/Homogeneous_differential_equation
1509 - M. Tenenbaum & H. Pollard, "Ordinary Differential Equations",
1510 Dover 1963, pp. 59
1512 # indirect doctest
1514 """
1515 hint = "1st_homogeneous_coeff_subs_dep_div_indep"
1516 has_integral = True
1517 order = [1]
1519 def _wilds(self, f, x, order):
1520 d = Wild('d', exclude=[f(x).diff(x), f(x).diff(x, 2)])
1521 e = Wild('e', exclude=[f(x).diff(x)])
1522 return d, e
1524 def _equation(self, fx, x, order):
1525 d, e = self.wilds()
1526 return d + e*fx.diff(x)
1528 def _verify(self, fx):
1529 self.d, self.e = self.wilds_match()
1530 self.y = Dummy('y')
1531 x = self.ode_problem.sym
1532 self.d = separatevars(self.d.subs(fx, self.y))
1533 self.e = separatevars(self.e.subs(fx, self.y))
1534 ordera = homogeneous_order(self.d, x, self.y)
1535 orderb = homogeneous_order(self.e, x, self.y)
1536 if ordera == orderb and ordera is not None:
1537 self.u = Dummy('u')
1538 if simplify((self.d + self.u*self.e).subs({x: 1, self.y: self.u})) != 0:
1539 return True
1540 return False
1541 return False
1543 def _get_match_object(self):
1544 fx = self.ode_problem.func
1545 x = self.ode_problem.sym
1546 self.u1 = Dummy('u1')
1547 xarg = 0
1548 yarg = 0
1549 return [self.d, self.e, fx, x, self.u, self.u1, self.y, xarg, yarg]
1551 def _get_general_solution(self, *, simplify_flag: bool = True):
1552 d, e, fx, x, u, u1, y, xarg, yarg = self._get_match_object()
1553 (C1,) = self.ode_problem.get_numbered_constants(num=1)
1554 int = Integral(
1555 (-e/(d + u1*e)).subs({x: 1, y: u1}),
1556 (u1, None, fx/x))
1557 sol = logcombine(Eq(log(x), int + log(C1)), force=True)
1558 gen_sol = sol.subs(fx, u).subs(((u, u - yarg), (x, x - xarg), (u, fx)))
1559 return [gen_sol]
1562class HomogeneousCoeffSubsIndepDivDep(SinglePatternODESolver):
1563 r"""
1564 Solves a 1st order differential equation with homogeneous coefficients
1565 using the substitution `u_2 = \frac{\text{<independent
1566 variable>}}{\text{<dependent variable>}}`.
1568 This is a differential equation
1570 .. math:: P(x, y) + Q(x, y) dy/dx = 0
1572 such that `P` and `Q` are homogeneous and of the same order. A function
1573 `F(x, y)` is homogeneous of order `n` if `F(x t, y t) = t^n F(x, y)`.
1574 Equivalently, `F(x, y)` can be rewritten as `G(y/x)` or `H(x/y)`. See
1575 also the docstring of :py:meth:`~sympy.solvers.ode.homogeneous_order`.
1577 If the coefficients `P` and `Q` in the differential equation above are
1578 homogeneous functions of the same order, then it can be shown that the
1579 substitution `x = u_2 y` (i.e. `u_2 = x/y`) will turn the differential
1580 equation into an equation separable in the variables `y` and `u_2`. If
1581 `h(u_2)` is the function that results from making the substitution `u_2 =
1582 x/f(x)` on `P(x, f(x))` and `g(u_2)` is the function that results from the
1583 substitution on `Q(x, f(x))` in the differential equation `P(x, f(x)) +
1584 Q(x, f(x)) f'(x) = 0`, then the general solution is:
1586 >>> from sympy import Function, dsolve, pprint
1587 >>> from sympy.abc import x
1588 >>> f, g, h = map(Function, ['f', 'g', 'h'])
1589 >>> genform = g(x/f(x)) + h(x/f(x))*f(x).diff(x)
1590 >>> pprint(genform)
1591 / x \ / x \ d
1592 g|----| + h|----|*--(f(x))
1593 \f(x)/ \f(x)/ dx
1594 >>> pprint(dsolve(genform, f(x),
1595 ... hint='1st_homogeneous_coeff_subs_indep_div_dep_Integral'))
1596 x
1597 ----
1598 f(x)
1599 /
1600 |
1601 | -g(u1)
1602 | ---------------- d(u1)
1603 | u1*g(u1) + h(u1)
1604 |
1605 /
1606 <BLANKLINE>
1607 f(x) = C1*e
1609 Where `u_1 g(u_1) + h(u_1) \ne 0` and `f(x) \ne 0`.
1611 See also the docstrings of
1612 :obj:`~sympy.solvers.ode.single.HomogeneousCoeffBest` and
1613 :obj:`~sympy.solvers.ode.single.HomogeneousCoeffSubsDepDivIndep`.
1615 Examples
1616 ========
1618 >>> from sympy import Function, pprint, dsolve
1619 >>> from sympy.abc import x
1620 >>> f = Function('f')
1621 >>> pprint(dsolve(2*x*f(x) + (x**2 + f(x)**2)*f(x).diff(x), f(x),
1622 ... hint='1st_homogeneous_coeff_subs_indep_div_dep',
1623 ... simplify=False))
1624 / 2 \
1625 | 3*x |
1626 log|----- + 1|
1627 | 2 |
1628 \f (x) /
1629 log(f(x)) = log(C1) - --------------
1630 3
1632 References
1633 ==========
1635 - https://en.wikipedia.org/wiki/Homogeneous_differential_equation
1636 - M. Tenenbaum & H. Pollard, "Ordinary Differential Equations",
1637 Dover 1963, pp. 59
1639 # indirect doctest
1641 """
1642 hint = "1st_homogeneous_coeff_subs_indep_div_dep"
1643 has_integral = True
1644 order = [1]
1646 def _wilds(self, f, x, order):
1647 d = Wild('d', exclude=[f(x).diff(x), f(x).diff(x, 2)])
1648 e = Wild('e', exclude=[f(x).diff(x)])
1649 return d, e
1651 def _equation(self, fx, x, order):
1652 d, e = self.wilds()
1653 return d + e*fx.diff(x)
1655 def _verify(self, fx):
1656 self.d, self.e = self.wilds_match()
1657 self.y = Dummy('y')
1658 x = self.ode_problem.sym
1659 self.d = separatevars(self.d.subs(fx, self.y))
1660 self.e = separatevars(self.e.subs(fx, self.y))
1661 ordera = homogeneous_order(self.d, x, self.y)
1662 orderb = homogeneous_order(self.e, x, self.y)
1663 if ordera == orderb and ordera is not None:
1664 self.u = Dummy('u')
1665 if simplify((self.e + self.u*self.d).subs({x: self.u, self.y: 1})) != 0:
1666 return True
1667 return False
1668 return False
1670 def _get_match_object(self):
1671 fx = self.ode_problem.func
1672 x = self.ode_problem.sym
1673 self.u1 = Dummy('u1')
1674 xarg = 0
1675 yarg = 0
1676 return [self.d, self.e, fx, x, self.u, self.u1, self.y, xarg, yarg]
1678 def _get_general_solution(self, *, simplify_flag: bool = True):
1679 d, e, fx, x, u, u1, y, xarg, yarg = self._get_match_object()
1680 (C1,) = self.ode_problem.get_numbered_constants(num=1)
1681 int = Integral(simplify((-d/(e + u1*d)).subs({x: u1, y: 1})), (u1, None, x/fx)) # type: ignore
1682 sol = logcombine(Eq(log(fx), int + log(C1)), force=True)
1683 gen_sol = sol.subs(fx, u).subs(((u, u - yarg), (x, x - xarg), (u, fx)))
1684 return [gen_sol]
1687class HomogeneousCoeffBest(HomogeneousCoeffSubsIndepDivDep, HomogeneousCoeffSubsDepDivIndep):
1688 r"""
1689 Returns the best solution to an ODE from the two hints
1690 ``1st_homogeneous_coeff_subs_dep_div_indep`` and
1691 ``1st_homogeneous_coeff_subs_indep_div_dep``.
1693 This is as determined by :py:meth:`~sympy.solvers.ode.ode.ode_sol_simplicity`.
1695 See the
1696 :obj:`~sympy.solvers.ode.single.HomogeneousCoeffSubsIndepDivDep`
1697 and
1698 :obj:`~sympy.solvers.ode.single.HomogeneousCoeffSubsDepDivIndep`
1699 docstrings for more information on these hints. Note that there is no
1700 ``ode_1st_homogeneous_coeff_best_Integral`` hint.
1702 Examples
1703 ========
1705 >>> from sympy import Function, dsolve, pprint
1706 >>> from sympy.abc import x
1707 >>> f = Function('f')
1708 >>> pprint(dsolve(2*x*f(x) + (x**2 + f(x)**2)*f(x).diff(x), f(x),
1709 ... hint='1st_homogeneous_coeff_best', simplify=False))
1710 / 2 \
1711 | 3*x |
1712 log|----- + 1|
1713 | 2 |
1714 \f (x) /
1715 log(f(x)) = log(C1) - --------------
1716 3
1718 References
1719 ==========
1721 - https://en.wikipedia.org/wiki/Homogeneous_differential_equation
1722 - M. Tenenbaum & H. Pollard, "Ordinary Differential Equations",
1723 Dover 1963, pp. 59
1725 # indirect doctest
1727 """
1728 hint = "1st_homogeneous_coeff_best"
1729 has_integral = False
1730 order = [1]
1732 def _verify(self, fx):
1733 if HomogeneousCoeffSubsIndepDivDep._verify(self, fx) and HomogeneousCoeffSubsDepDivIndep._verify(self, fx):
1734 return True
1735 return False
1737 def _get_general_solution(self, *, simplify_flag: bool = True):
1738 # There are two substitutions that solve the equation, u1=y/x and u2=x/y
1739 # # They produce different integrals, so try them both and see which
1740 # # one is easier
1741 sol1 = HomogeneousCoeffSubsIndepDivDep._get_general_solution(self)
1742 sol2 = HomogeneousCoeffSubsDepDivIndep._get_general_solution(self)
1743 fx = self.ode_problem.func
1744 if simplify_flag:
1745 sol1 = odesimp(self.ode_problem.eq, *sol1, fx, "1st_homogeneous_coeff_subs_indep_div_dep")
1746 sol2 = odesimp(self.ode_problem.eq, *sol2, fx, "1st_homogeneous_coeff_subs_dep_div_indep")
1747 return min([sol1, sol2], key=lambda x: ode_sol_simplicity(x, fx, trysolving=not simplify))
1750class LinearCoefficients(HomogeneousCoeffBest):
1751 r"""
1752 Solves a differential equation with linear coefficients.
1754 The general form of a differential equation with linear coefficients is
1756 .. math:: y' + F\left(\!\frac{a_1 x + b_1 y + c_1}{a_2 x + b_2 y +
1757 c_2}\!\right) = 0\text{,}
1759 where `a_1`, `b_1`, `c_1`, `a_2`, `b_2`, `c_2` are constants and `a_1 b_2
1760 - a_2 b_1 \ne 0`.
1762 This can be solved by substituting:
1764 .. math:: x = x' + \frac{b_2 c_1 - b_1 c_2}{a_2 b_1 - a_1 b_2}
1766 y = y' + \frac{a_1 c_2 - a_2 c_1}{a_2 b_1 - a_1
1767 b_2}\text{.}
1769 This substitution reduces the equation to a homogeneous differential
1770 equation.
1772 See Also
1773 ========
1774 :obj:`sympy.solvers.ode.single.HomogeneousCoeffBest`
1775 :obj:`sympy.solvers.ode.single.HomogeneousCoeffSubsIndepDivDep`
1776 :obj:`sympy.solvers.ode.single.HomogeneousCoeffSubsDepDivIndep`
1778 Examples
1779 ========
1781 >>> from sympy import dsolve, Function, pprint
1782 >>> from sympy.abc import x
1783 >>> f = Function('f')
1784 >>> df = f(x).diff(x)
1785 >>> eq = (x + f(x) + 1)*df + (f(x) - 6*x + 1)
1786 >>> dsolve(eq, hint='linear_coefficients')
1787 [Eq(f(x), -x - sqrt(C1 + 7*x**2) - 1), Eq(f(x), -x + sqrt(C1 + 7*x**2) - 1)]
1788 >>> pprint(dsolve(eq, hint='linear_coefficients'))
1789 ___________ ___________
1790 / 2 / 2
1791 [f(x) = -x - \/ C1 + 7*x - 1, f(x) = -x + \/ C1 + 7*x - 1]
1794 References
1795 ==========
1797 - Joel Moses, "Symbolic Integration - The Stormy Decade", Communications
1798 of the ACM, Volume 14, Number 8, August 1971, pp. 558
1799 """
1800 hint = "linear_coefficients"
1801 has_integral = True
1802 order = [1]
1804 def _wilds(self, f, x, order):
1805 d = Wild('d', exclude=[f(x).diff(x), f(x).diff(x, 2)])
1806 e = Wild('e', exclude=[f(x).diff(x)])
1807 return d, e
1809 def _equation(self, fx, x, order):
1810 d, e = self.wilds()
1811 return d + e*fx.diff(x)
1813 def _verify(self, fx):
1814 self.d, self.e = self.wilds_match()
1815 a, b = self.wilds()
1816 F = self.d/self.e
1817 x = self.ode_problem.sym
1818 params = self._linear_coeff_match(F, fx)
1819 if params:
1820 self.xarg, self.yarg = params
1821 u = Dummy('u')
1822 t = Dummy('t')
1823 self.y = Dummy('y')
1824 # Dummy substitution for df and f(x).
1825 dummy_eq = self.ode_problem.eq.subs(((fx.diff(x), t), (fx, u)))
1826 reps = ((x, x + self.xarg), (u, u + self.yarg), (t, fx.diff(x)), (u, fx))
1827 dummy_eq = simplify(dummy_eq.subs(reps))
1828 # get the re-cast values for e and d
1829 r2 = collect(expand(dummy_eq), [fx.diff(x), fx]).match(a*fx.diff(x) + b)
1830 if r2:
1831 self.d, self.e = r2[b], r2[a]
1832 orderd = homogeneous_order(self.d, x, fx)
1833 ordere = homogeneous_order(self.e, x, fx)
1834 if orderd == ordere and orderd is not None:
1835 self.d = self.d.subs(fx, self.y)
1836 self.e = self.e.subs(fx, self.y)
1837 return True
1838 return False
1839 return False
1841 def _linear_coeff_match(self, expr, func):
1842 r"""
1843 Helper function to match hint ``linear_coefficients``.
1845 Matches the expression to the form `(a_1 x + b_1 f(x) + c_1)/(a_2 x + b_2
1846 f(x) + c_2)` where the following conditions hold:
1848 1. `a_1`, `b_1`, `c_1`, `a_2`, `b_2`, `c_2` are Rationals;
1849 2. `c_1` or `c_2` are not equal to zero;
1850 3. `a_2 b_1 - a_1 b_2` is not equal to zero.
1852 Return ``xarg``, ``yarg`` where
1854 1. ``xarg`` = `(b_2 c_1 - b_1 c_2)/(a_2 b_1 - a_1 b_2)`
1855 2. ``yarg`` = `(a_1 c_2 - a_2 c_1)/(a_2 b_1 - a_1 b_2)`
1858 Examples
1859 ========
1861 >>> from sympy import Function, sin
1862 >>> from sympy.abc import x
1863 >>> from sympy.solvers.ode.single import LinearCoefficients
1864 >>> f = Function('f')
1865 >>> eq = (-25*f(x) - 8*x + 62)/(4*f(x) + 11*x - 11)
1866 >>> obj = LinearCoefficients(eq)
1867 >>> obj._linear_coeff_match(eq, f(x))
1868 (1/9, 22/9)
1869 >>> eq = sin((-5*f(x) - 8*x + 6)/(4*f(x) + x - 1))
1870 >>> obj = LinearCoefficients(eq)
1871 >>> obj._linear_coeff_match(eq, f(x))
1872 (19/27, 2/27)
1873 >>> eq = sin(f(x)/x)
1874 >>> obj = LinearCoefficients(eq)
1875 >>> obj._linear_coeff_match(eq, f(x))
1877 """
1878 f = func.func
1879 x = func.args[0]
1880 def abc(eq):
1881 r'''
1882 Internal function of _linear_coeff_match
1883 that returns Rationals a, b, c
1884 if eq is a*x + b*f(x) + c, else None.
1885 '''
1886 eq = _mexpand(eq)
1887 c = eq.as_independent(x, f(x), as_Add=True)[0]
1888 if not c.is_Rational:
1889 return
1890 a = eq.coeff(x)
1891 if not a.is_Rational:
1892 return
1893 b = eq.coeff(f(x))
1894 if not b.is_Rational:
1895 return
1896 if eq == a*x + b*f(x) + c:
1897 return a, b, c
1899 def match(arg):
1900 r'''
1901 Internal function of _linear_coeff_match that returns Rationals a1,
1902 b1, c1, a2, b2, c2 and a2*b1 - a1*b2 of the expression (a1*x + b1*f(x)
1903 + c1)/(a2*x + b2*f(x) + c2) if one of c1 or c2 and a2*b1 - a1*b2 is
1904 non-zero, else None.
1905 '''
1906 n, d = arg.together().as_numer_denom()
1907 m = abc(n)
1908 if m is not None:
1909 a1, b1, c1 = m
1910 m = abc(d)
1911 if m is not None:
1912 a2, b2, c2 = m
1913 d = a2*b1 - a1*b2
1914 if (c1 or c2) and d:
1915 return a1, b1, c1, a2, b2, c2, d
1917 m = [fi.args[0] for fi in expr.atoms(Function) if fi.func != f and
1918 len(fi.args) == 1 and not fi.args[0].is_Function] or {expr}
1919 m1 = match(m.pop())
1920 if m1 and all(match(mi) == m1 for mi in m):
1921 a1, b1, c1, a2, b2, c2, denom = m1
1922 return (b2*c1 - b1*c2)/denom, (a1*c2 - a2*c1)/denom
1924 def _get_match_object(self):
1925 fx = self.ode_problem.func
1926 x = self.ode_problem.sym
1927 self.u1 = Dummy('u1')
1928 u = Dummy('u')
1929 return [self.d, self.e, fx, x, u, self.u1, self.y, self.xarg, self.yarg]
1932class NthOrderReducible(SingleODESolver):
1933 r"""
1934 Solves ODEs that only involve derivatives of the dependent variable using
1935 a substitution of the form `f^n(x) = g(x)`.
1937 For example any second order ODE of the form `f''(x) = h(f'(x), x)` can be
1938 transformed into a pair of 1st order ODEs `g'(x) = h(g(x), x)` and
1939 `f'(x) = g(x)`. Usually the 1st order ODE for `g` is easier to solve. If
1940 that gives an explicit solution for `g` then `f` is found simply by
1941 integration.
1944 Examples
1945 ========
1947 >>> from sympy import Function, dsolve, Eq
1948 >>> from sympy.abc import x
1949 >>> f = Function('f')
1950 >>> eq = Eq(x*f(x).diff(x)**2 + f(x).diff(x, 2), 0)
1951 >>> dsolve(eq, f(x), hint='nth_order_reducible')
1952 ... # doctest: +NORMALIZE_WHITESPACE
1953 Eq(f(x), C1 - sqrt(-1/C2)*log(-C2*sqrt(-1/C2) + x) + sqrt(-1/C2)*log(C2*sqrt(-1/C2) + x))
1955 """
1956 hint = "nth_order_reducible"
1957 has_integral = False
1959 def _matches(self):
1960 # Any ODE that can be solved with a substitution and
1961 # repeated integration e.g.:
1962 # `d^2/dx^2(y) + x*d/dx(y) = constant
1963 #f'(x) must be finite for this to work
1964 eq = self.ode_problem.eq_preprocessed
1965 func = self.ode_problem.func
1966 x = self.ode_problem.sym
1967 r"""
1968 Matches any differential equation that can be rewritten with a smaller
1969 order. Only derivatives of ``func`` alone, wrt a single variable,
1970 are considered, and only in them should ``func`` appear.
1971 """
1972 # ODE only handles functions of 1 variable so this affirms that state
1973 assert len(func.args) == 1
1974 vc = [d.variable_count[0] for d in eq.atoms(Derivative)
1975 if d.expr == func and len(d.variable_count) == 1]
1976 ords = [c for v, c in vc if v == x]
1977 if len(ords) < 2:
1978 return False
1979 self.smallest = min(ords)
1980 # make sure func does not appear outside of derivatives
1981 D = Dummy()
1982 if eq.subs(func.diff(x, self.smallest), D).has(func):
1983 return False
1984 return True
1986 def _get_general_solution(self, *, simplify_flag: bool = True):
1987 eq = self.ode_problem.eq
1988 f = self.ode_problem.func.func
1989 x = self.ode_problem.sym
1990 n = self.smallest
1991 # get a unique function name for g
1992 names = [a.name for a in eq.atoms(AppliedUndef)]
1993 while True:
1994 name = Dummy().name
1995 if name not in names:
1996 g = Function(name)
1997 break
1998 w = f(x).diff(x, n)
1999 geq = eq.subs(w, g(x))
2000 gsol = dsolve(geq, g(x))
2002 if not isinstance(gsol, list):
2003 gsol = [gsol]
2005 # Might be multiple solutions to the reduced ODE:
2006 fsol = []
2007 for gsoli in gsol:
2008 fsoli = dsolve(gsoli.subs(g(x), w), f(x)) # or do integration n times
2009 fsol.append(fsoli)
2011 return fsol
2014class SecondHypergeometric(SingleODESolver):
2015 r"""
2016 Solves 2nd order linear differential equations.
2018 It computes special function solutions which can be expressed using the
2019 2F1, 1F1 or 0F1 hypergeometric functions.
2021 .. math:: y'' + A(x) y' + B(x) y = 0\text{,}
2023 where `A` and `B` are rational functions.
2025 These kinds of differential equations have solution of non-Liouvillian form.
2027 Given linear ODE can be obtained from 2F1 given by
2029 .. math:: (x^2 - x) y'' + ((a + b + 1) x - c) y' + b a y = 0\text{,}
2031 where {a, b, c} are arbitrary constants.
2033 Notes
2034 =====
2036 The algorithm should find any solution of the form
2038 .. math:: y = P(x) _pF_q(..; ..;\frac{\alpha x^k + \beta}{\gamma x^k + \delta})\text{,}
2040 where pFq is any of 2F1, 1F1 or 0F1 and `P` is an "arbitrary function".
2041 Currently only the 2F1 case is implemented in SymPy but the other cases are
2042 described in the paper and could be implemented in future (contributions
2043 welcome!).
2046 Examples
2047 ========
2049 >>> from sympy import Function, dsolve, pprint
2050 >>> from sympy.abc import x
2051 >>> f = Function('f')
2052 >>> eq = (x*x - x)*f(x).diff(x,2) + (5*x - 1)*f(x).diff(x) + 4*f(x)
2053 >>> pprint(dsolve(eq, f(x), '2nd_hypergeometric'))
2054 _
2055 / / 4 \\ |_ /-1, -1 | \
2056 |C1 + C2*|log(x) + -----||* | | | x|
2057 \ \ x + 1// 2 1 \ 1 | /
2058 f(x) = --------------------------------------------
2059 3
2060 (x - 1)
2063 References
2064 ==========
2066 - "Non-Liouvillian solutions for second order linear ODEs" by L. Chan, E.S. Cheb-Terrab
2068 """
2069 hint = "2nd_hypergeometric"
2070 has_integral = True
2072 def _matches(self):
2073 eq = self.ode_problem.eq_preprocessed
2074 func = self.ode_problem.func
2075 r = match_2nd_hypergeometric(eq, func)
2076 self.match_object = None
2077 if r:
2078 A, B = r
2079 d = equivalence_hypergeometric(A, B, func)
2080 if d:
2081 if d['type'] == "2F1":
2082 self.match_object = match_2nd_2F1_hypergeometric(d['I0'], d['k'], d['sing_point'], func)
2083 if self.match_object is not None:
2084 self.match_object.update({'A':A, 'B':B})
2085 # We can extend it for 1F1 and 0F1 type also.
2086 return self.match_object is not None
2088 def _get_general_solution(self, *, simplify_flag: bool = True):
2089 eq = self.ode_problem.eq
2090 func = self.ode_problem.func
2091 if self.match_object['type'] == "2F1":
2092 sol = get_sol_2F1_hypergeometric(eq, func, self.match_object)
2093 if sol is None:
2094 raise NotImplementedError("The given ODE " + str(eq) + " cannot be solved by"
2095 + " the hypergeometric method")
2097 return [sol]
2100class NthLinearConstantCoeffHomogeneous(SingleODESolver):
2101 r"""
2102 Solves an `n`\th order linear homogeneous differential equation with
2103 constant coefficients.
2105 This is an equation of the form
2107 .. math:: a_n f^{(n)}(x) + a_{n-1} f^{(n-1)}(x) + \cdots + a_1 f'(x)
2108 + a_0 f(x) = 0\text{.}
2110 These equations can be solved in a general manner, by taking the roots of
2111 the characteristic equation `a_n m^n + a_{n-1} m^{n-1} + \cdots + a_1 m +
2112 a_0 = 0`. The solution will then be the sum of `C_n x^i e^{r x}` terms,
2113 for each where `C_n` is an arbitrary constant, `r` is a root of the
2114 characteristic equation and `i` is one of each from 0 to the multiplicity
2115 of the root - 1 (for example, a root 3 of multiplicity 2 would create the
2116 terms `C_1 e^{3 x} + C_2 x e^{3 x}`). The exponential is usually expanded
2117 for complex roots using Euler's equation `e^{I x} = \cos(x) + I \sin(x)`.
2118 Complex roots always come in conjugate pairs in polynomials with real
2119 coefficients, so the two roots will be represented (after simplifying the
2120 constants) as `e^{a x} \left(C_1 \cos(b x) + C_2 \sin(b x)\right)`.
2122 If SymPy cannot find exact roots to the characteristic equation, a
2123 :py:class:`~sympy.polys.rootoftools.ComplexRootOf` instance will be return
2124 instead.
2126 >>> from sympy import Function, dsolve
2127 >>> from sympy.abc import x
2128 >>> f = Function('f')
2129 >>> dsolve(f(x).diff(x, 5) + 10*f(x).diff(x) - 2*f(x), f(x),
2130 ... hint='nth_linear_constant_coeff_homogeneous')
2131 ... # doctest: +NORMALIZE_WHITESPACE
2132 Eq(f(x), C5*exp(x*CRootOf(_x**5 + 10*_x - 2, 0))
2133 + (C1*sin(x*im(CRootOf(_x**5 + 10*_x - 2, 1)))
2134 + C2*cos(x*im(CRootOf(_x**5 + 10*_x - 2, 1))))*exp(x*re(CRootOf(_x**5 + 10*_x - 2, 1)))
2135 + (C3*sin(x*im(CRootOf(_x**5 + 10*_x - 2, 3)))
2136 + C4*cos(x*im(CRootOf(_x**5 + 10*_x - 2, 3))))*exp(x*re(CRootOf(_x**5 + 10*_x - 2, 3))))
2138 Note that because this method does not involve integration, there is no
2139 ``nth_linear_constant_coeff_homogeneous_Integral`` hint.
2141 Examples
2142 ========
2144 >>> from sympy import Function, dsolve, pprint
2145 >>> from sympy.abc import x
2146 >>> f = Function('f')
2147 >>> pprint(dsolve(f(x).diff(x, 4) + 2*f(x).diff(x, 3) -
2148 ... 2*f(x).diff(x, 2) - 6*f(x).diff(x) + 5*f(x), f(x),
2149 ... hint='nth_linear_constant_coeff_homogeneous'))
2150 x -2*x
2151 f(x) = (C1 + C2*x)*e + (C3*sin(x) + C4*cos(x))*e
2153 References
2154 ==========
2156 - https://en.wikipedia.org/wiki/Linear_differential_equation section:
2157 Nonhomogeneous_equation_with_constant_coefficients
2158 - M. Tenenbaum & H. Pollard, "Ordinary Differential Equations",
2159 Dover 1963, pp. 211
2161 # indirect doctest
2163 """
2164 hint = "nth_linear_constant_coeff_homogeneous"
2165 has_integral = False
2167 def _matches(self):
2168 eq = self.ode_problem.eq_high_order_free
2169 func = self.ode_problem.func
2170 order = self.ode_problem.order
2171 x = self.ode_problem.sym
2172 self.r = self.ode_problem.get_linear_coefficients(eq, func, order)
2173 if order and self.r and not any(self.r[i].has(x) for i in self.r if i >= 0):
2174 if not self.r[-1]:
2175 return True
2176 else:
2177 return False
2178 return False
2180 def _get_general_solution(self, *, simplify_flag: bool = True):
2181 fx = self.ode_problem.func
2182 order = self.ode_problem.order
2183 roots, collectterms = _get_const_characteristic_eq_sols(self.r, fx, order)
2184 # A generator of constants
2185 constants = self.ode_problem.get_numbered_constants(num=len(roots))
2186 gsol = Add(*[i*j for (i, j) in zip(constants, roots)])
2187 gsol = Eq(fx, gsol)
2188 if simplify_flag:
2189 gsol = _get_simplified_sol([gsol], fx, collectterms)
2191 return [gsol]
2194class NthLinearConstantCoeffVariationOfParameters(SingleODESolver):
2195 r"""
2196 Solves an `n`\th order linear differential equation with constant
2197 coefficients using the method of variation of parameters.
2199 This method works on any differential equations of the form
2201 .. math:: f^{(n)}(x) + a_{n-1} f^{(n-1)}(x) + \cdots + a_1 f'(x) + a_0
2202 f(x) = P(x)\text{.}
2204 This method works by assuming that the particular solution takes the form
2206 .. math:: \sum_{x=1}^{n} c_i(x) y_i(x)\text{,}
2208 where `y_i` is the `i`\th solution to the homogeneous equation. The
2209 solution is then solved using Wronskian's and Cramer's Rule. The
2210 particular solution is given by
2212 .. math:: \sum_{x=1}^n \left( \int \frac{W_i(x)}{W(x)} \,dx
2213 \right) y_i(x) \text{,}
2215 where `W(x)` is the Wronskian of the fundamental system (the system of `n`
2216 linearly independent solutions to the homogeneous equation), and `W_i(x)`
2217 is the Wronskian of the fundamental system with the `i`\th column replaced
2218 with `[0, 0, \cdots, 0, P(x)]`.
2220 This method is general enough to solve any `n`\th order inhomogeneous
2221 linear differential equation with constant coefficients, but sometimes
2222 SymPy cannot simplify the Wronskian well enough to integrate it. If this
2223 method hangs, try using the
2224 ``nth_linear_constant_coeff_variation_of_parameters_Integral`` hint and
2225 simplifying the integrals manually. Also, prefer using
2226 ``nth_linear_constant_coeff_undetermined_coefficients`` when it
2227 applies, because it does not use integration, making it faster and more
2228 reliable.
2230 Warning, using simplify=False with
2231 'nth_linear_constant_coeff_variation_of_parameters' in
2232 :py:meth:`~sympy.solvers.ode.dsolve` may cause it to hang, because it will
2233 not attempt to simplify the Wronskian before integrating. It is
2234 recommended that you only use simplify=False with
2235 'nth_linear_constant_coeff_variation_of_parameters_Integral' for this
2236 method, especially if the solution to the homogeneous equation has
2237 trigonometric functions in it.
2239 Examples
2240 ========
2242 >>> from sympy import Function, dsolve, pprint, exp, log
2243 >>> from sympy.abc import x
2244 >>> f = Function('f')
2245 >>> pprint(dsolve(f(x).diff(x, 3) - 3*f(x).diff(x, 2) +
2246 ... 3*f(x).diff(x) - f(x) - exp(x)*log(x), f(x),
2247 ... hint='nth_linear_constant_coeff_variation_of_parameters'))
2248 / / / x*log(x) 11*x\\\ x
2249 f(x) = |C1 + x*|C2 + x*|C3 + -------- - ----|||*e
2250 \ \ \ 6 36 ///
2252 References
2253 ==========
2255 - https://en.wikipedia.org/wiki/Variation_of_parameters
2256 - https://planetmath.org/VariationOfParameters
2257 - M. Tenenbaum & H. Pollard, "Ordinary Differential Equations",
2258 Dover 1963, pp. 233
2260 # indirect doctest
2262 """
2263 hint = "nth_linear_constant_coeff_variation_of_parameters"
2264 has_integral = True
2266 def _matches(self):
2267 eq = self.ode_problem.eq_high_order_free
2268 func = self.ode_problem.func
2269 order = self.ode_problem.order
2270 x = self.ode_problem.sym
2271 self.r = self.ode_problem.get_linear_coefficients(eq, func, order)
2273 if order and self.r and not any(self.r[i].has(x) for i in self.r if i >= 0):
2274 if self.r[-1]:
2275 return True
2276 else:
2277 return False
2278 return False
2280 def _get_general_solution(self, *, simplify_flag: bool = True):
2281 eq = self.ode_problem.eq_high_order_free
2282 f = self.ode_problem.func.func
2283 x = self.ode_problem.sym
2284 order = self.ode_problem.order
2285 roots, collectterms = _get_const_characteristic_eq_sols(self.r, f(x), order)
2286 # A generator of constants
2287 constants = self.ode_problem.get_numbered_constants(num=len(roots))
2288 homogen_sol = Add(*[i*j for (i, j) in zip(constants, roots)])
2289 homogen_sol = Eq(f(x), homogen_sol)
2290 homogen_sol = _solve_variation_of_parameters(eq, f(x), roots, homogen_sol, order, self.r, simplify_flag)
2291 if simplify_flag:
2292 homogen_sol = _get_simplified_sol([homogen_sol], f(x), collectterms)
2293 return [homogen_sol]
2296class NthLinearConstantCoeffUndeterminedCoefficients(SingleODESolver):
2297 r"""
2298 Solves an `n`\th order linear differential equation with constant
2299 coefficients using the method of undetermined coefficients.
2301 This method works on differential equations of the form
2303 .. math:: a_n f^{(n)}(x) + a_{n-1} f^{(n-1)}(x) + \cdots + a_1 f'(x)
2304 + a_0 f(x) = P(x)\text{,}
2306 where `P(x)` is a function that has a finite number of linearly
2307 independent derivatives.
2309 Functions that fit this requirement are finite sums functions of the form
2310 `a x^i e^{b x} \sin(c x + d)` or `a x^i e^{b x} \cos(c x + d)`, where `i`
2311 is a non-negative integer and `a`, `b`, `c`, and `d` are constants. For
2312 example any polynomial in `x`, functions like `x^2 e^{2 x}`, `x \sin(x)`,
2313 and `e^x \cos(x)` can all be used. Products of `\sin`'s and `\cos`'s have
2314 a finite number of derivatives, because they can be expanded into `\sin(a
2315 x)` and `\cos(b x)` terms. However, SymPy currently cannot do that
2316 expansion, so you will need to manually rewrite the expression in terms of
2317 the above to use this method. So, for example, you will need to manually
2318 convert `\sin^2(x)` into `(1 + \cos(2 x))/2` to properly apply the method
2319 of undetermined coefficients on it.
2321 This method works by creating a trial function from the expression and all
2322 of its linear independent derivatives and substituting them into the
2323 original ODE. The coefficients for each term will be a system of linear
2324 equations, which are be solved for and substituted, giving the solution.
2325 If any of the trial functions are linearly dependent on the solution to
2326 the homogeneous equation, they are multiplied by sufficient `x` to make
2327 them linearly independent.
2329 Examples
2330 ========
2332 >>> from sympy import Function, dsolve, pprint, exp, cos
2333 >>> from sympy.abc import x
2334 >>> f = Function('f')
2335 >>> pprint(dsolve(f(x).diff(x, 2) + 2*f(x).diff(x) + f(x) -
2336 ... 4*exp(-x)*x**2 + cos(2*x), f(x),
2337 ... hint='nth_linear_constant_coeff_undetermined_coefficients'))
2338 / / 3\\
2339 | | x || -x 4*sin(2*x) 3*cos(2*x)
2340 f(x) = |C1 + x*|C2 + --||*e - ---------- + ----------
2341 \ \ 3 // 25 25
2343 References
2344 ==========
2346 - https://en.wikipedia.org/wiki/Method_of_undetermined_coefficients
2347 - M. Tenenbaum & H. Pollard, "Ordinary Differential Equations",
2348 Dover 1963, pp. 221
2350 # indirect doctest
2352 """
2353 hint = "nth_linear_constant_coeff_undetermined_coefficients"
2354 has_integral = False
2356 def _matches(self):
2357 eq = self.ode_problem.eq_high_order_free
2358 func = self.ode_problem.func
2359 order = self.ode_problem.order
2360 x = self.ode_problem.sym
2361 self.r = self.ode_problem.get_linear_coefficients(eq, func, order)
2362 does_match = False
2363 if order and self.r and not any(self.r[i].has(x) for i in self.r if i >= 0):
2364 if self.r[-1]:
2365 eq_homogeneous = Add(eq, -self.r[-1])
2366 undetcoeff = _undetermined_coefficients_match(self.r[-1], x, func, eq_homogeneous)
2367 if undetcoeff['test']:
2368 self.trialset = undetcoeff['trialset']
2369 does_match = True
2370 return does_match
2372 def _get_general_solution(self, *, simplify_flag: bool = True):
2373 eq = self.ode_problem.eq
2374 f = self.ode_problem.func.func
2375 x = self.ode_problem.sym
2376 order = self.ode_problem.order
2377 roots, collectterms = _get_const_characteristic_eq_sols(self.r, f(x), order)
2378 # A generator of constants
2379 constants = self.ode_problem.get_numbered_constants(num=len(roots))
2380 homogen_sol = Add(*[i*j for (i, j) in zip(constants, roots)])
2381 homogen_sol = Eq(f(x), homogen_sol)
2382 self.r.update({'list': roots, 'sol': homogen_sol, 'simpliy_flag': simplify_flag})
2383 gsol = _solve_undetermined_coefficients(eq, f(x), order, self.r, self.trialset)
2384 if simplify_flag:
2385 gsol = _get_simplified_sol([gsol], f(x), collectterms)
2386 return [gsol]
2389class NthLinearEulerEqHomogeneous(SingleODESolver):
2390 r"""
2391 Solves an `n`\th order linear homogeneous variable-coefficient
2392 Cauchy-Euler equidimensional ordinary differential equation.
2394 This is an equation with form `0 = a_0 f(x) + a_1 x f'(x) + a_2 x^2 f''(x)
2395 \cdots`.
2397 These equations can be solved in a general manner, by substituting
2398 solutions of the form `f(x) = x^r`, and deriving a characteristic equation
2399 for `r`. When there are repeated roots, we include extra terms of the
2400 form `C_{r k} \ln^k(x) x^r`, where `C_{r k}` is an arbitrary integration
2401 constant, `r` is a root of the characteristic equation, and `k` ranges
2402 over the multiplicity of `r`. In the cases where the roots are complex,
2403 solutions of the form `C_1 x^a \sin(b \log(x)) + C_2 x^a \cos(b \log(x))`
2404 are returned, based on expansions with Euler's formula. The general
2405 solution is the sum of the terms found. If SymPy cannot find exact roots
2406 to the characteristic equation, a
2407 :py:obj:`~.ComplexRootOf` instance will be returned
2408 instead.
2410 >>> from sympy import Function, dsolve
2411 >>> from sympy.abc import x
2412 >>> f = Function('f')
2413 >>> dsolve(4*x**2*f(x).diff(x, 2) + f(x), f(x),
2414 ... hint='nth_linear_euler_eq_homogeneous')
2415 ... # doctest: +NORMALIZE_WHITESPACE
2416 Eq(f(x), sqrt(x)*(C1 + C2*log(x)))
2418 Note that because this method does not involve integration, there is no
2419 ``nth_linear_euler_eq_homogeneous_Integral`` hint.
2421 The following is for internal use:
2423 - ``returns = 'sol'`` returns the solution to the ODE.
2424 - ``returns = 'list'`` returns a list of linearly independent solutions,
2425 corresponding to the fundamental solution set, for use with non
2426 homogeneous solution methods like variation of parameters and
2427 undetermined coefficients. Note that, though the solutions should be
2428 linearly independent, this function does not explicitly check that. You
2429 can do ``assert simplify(wronskian(sollist)) != 0`` to check for linear
2430 independence. Also, ``assert len(sollist) == order`` will need to pass.
2431 - ``returns = 'both'``, return a dictionary ``{'sol': <solution to ODE>,
2432 'list': <list of linearly independent solutions>}``.
2434 Examples
2435 ========
2437 >>> from sympy import Function, dsolve, pprint
2438 >>> from sympy.abc import x
2439 >>> f = Function('f')
2440 >>> eq = f(x).diff(x, 2)*x**2 - 4*f(x).diff(x)*x + 6*f(x)
2441 >>> pprint(dsolve(eq, f(x),
2442 ... hint='nth_linear_euler_eq_homogeneous'))
2443 2
2444 f(x) = x *(C1 + C2*x)
2446 References
2447 ==========
2449 - https://en.wikipedia.org/wiki/Cauchy%E2%80%93Euler_equation
2450 - C. Bender & S. Orszag, "Advanced Mathematical Methods for Scientists and
2451 Engineers", Springer 1999, pp. 12
2453 # indirect doctest
2455 """
2456 hint = "nth_linear_euler_eq_homogeneous"
2457 has_integral = False
2459 def _matches(self):
2460 eq = self.ode_problem.eq_preprocessed
2461 f = self.ode_problem.func.func
2462 order = self.ode_problem.order
2463 x = self.ode_problem.sym
2464 match = self.ode_problem.get_linear_coefficients(eq, f(x), order)
2465 self.r = None
2466 does_match = False
2468 if order and match:
2469 coeff = match[order]
2470 factor = x**order / coeff
2471 self.r = {i: factor*match[i] for i in match}
2472 if self.r and all(_test_term(self.r[i], f(x), i) for i in
2473 self.r if i >= 0):
2474 if not self.r[-1]:
2475 does_match = True
2476 return does_match
2478 def _get_general_solution(self, *, simplify_flag: bool = True):
2479 fx = self.ode_problem.func
2480 eq = self.ode_problem.eq
2481 homogen_sol = _get_euler_characteristic_eq_sols(eq, fx, self.r)[0]
2482 return [homogen_sol]
2485class NthLinearEulerEqNonhomogeneousVariationOfParameters(SingleODESolver):
2486 r"""
2487 Solves an `n`\th order linear non homogeneous Cauchy-Euler equidimensional
2488 ordinary differential equation using variation of parameters.
2490 This is an equation with form `g(x) = a_0 f(x) + a_1 x f'(x) + a_2 x^2 f''(x)
2491 \cdots`.
2493 This method works by assuming that the particular solution takes the form
2495 .. math:: \sum_{x=1}^{n} c_i(x) y_i(x) {a_n} {x^n} \text{, }
2497 where `y_i` is the `i`\th solution to the homogeneous equation. The
2498 solution is then solved using Wronskian's and Cramer's Rule. The
2499 particular solution is given by multiplying eq given below with `a_n x^{n}`
2501 .. math:: \sum_{x=1}^n \left( \int \frac{W_i(x)}{W(x)} \, dx
2502 \right) y_i(x) \text{, }
2504 where `W(x)` is the Wronskian of the fundamental system (the system of `n`
2505 linearly independent solutions to the homogeneous equation), and `W_i(x)`
2506 is the Wronskian of the fundamental system with the `i`\th column replaced
2507 with `[0, 0, \cdots, 0, \frac{x^{- n}}{a_n} g{\left(x \right)}]`.
2509 This method is general enough to solve any `n`\th order inhomogeneous
2510 linear differential equation, but sometimes SymPy cannot simplify the
2511 Wronskian well enough to integrate it. If this method hangs, try using the
2512 ``nth_linear_constant_coeff_variation_of_parameters_Integral`` hint and
2513 simplifying the integrals manually. Also, prefer using
2514 ``nth_linear_constant_coeff_undetermined_coefficients`` when it
2515 applies, because it does not use integration, making it faster and more
2516 reliable.
2518 Warning, using simplify=False with
2519 'nth_linear_constant_coeff_variation_of_parameters' in
2520 :py:meth:`~sympy.solvers.ode.dsolve` may cause it to hang, because it will
2521 not attempt to simplify the Wronskian before integrating. It is
2522 recommended that you only use simplify=False with
2523 'nth_linear_constant_coeff_variation_of_parameters_Integral' for this
2524 method, especially if the solution to the homogeneous equation has
2525 trigonometric functions in it.
2527 Examples
2528 ========
2530 >>> from sympy import Function, dsolve, Derivative
2531 >>> from sympy.abc import x
2532 >>> f = Function('f')
2533 >>> eq = x**2*Derivative(f(x), x, x) - 2*x*Derivative(f(x), x) + 2*f(x) - x**4
2534 >>> dsolve(eq, f(x),
2535 ... hint='nth_linear_euler_eq_nonhomogeneous_variation_of_parameters').expand()
2536 Eq(f(x), C1*x + C2*x**2 + x**4/6)
2538 """
2539 hint = "nth_linear_euler_eq_nonhomogeneous_variation_of_parameters"
2540 has_integral = True
2542 def _matches(self):
2543 eq = self.ode_problem.eq_preprocessed
2544 f = self.ode_problem.func.func
2545 order = self.ode_problem.order
2546 x = self.ode_problem.sym
2547 match = self.ode_problem.get_linear_coefficients(eq, f(x), order)
2548 self.r = None
2549 does_match = False
2551 if order and match:
2552 coeff = match[order]
2553 factor = x**order / coeff
2554 self.r = {i: factor*match[i] for i in match}
2555 if self.r and all(_test_term(self.r[i], f(x), i) for i in
2556 self.r if i >= 0):
2557 if self.r[-1]:
2558 does_match = True
2560 return does_match
2562 def _get_general_solution(self, *, simplify_flag: bool = True):
2563 eq = self.ode_problem.eq
2564 f = self.ode_problem.func.func
2565 x = self.ode_problem.sym
2566 order = self.ode_problem.order
2567 homogen_sol, roots = _get_euler_characteristic_eq_sols(eq, f(x), self.r)
2568 self.r[-1] = self.r[-1]/self.r[order]
2569 sol = _solve_variation_of_parameters(eq, f(x), roots, homogen_sol, order, self.r, simplify_flag)
2571 return [Eq(f(x), homogen_sol.rhs + (sol.rhs - homogen_sol.rhs)*self.r[order])]
2574class NthLinearEulerEqNonhomogeneousUndeterminedCoefficients(SingleODESolver):
2575 r"""
2576 Solves an `n`\th order linear non homogeneous Cauchy-Euler equidimensional
2577 ordinary differential equation using undetermined coefficients.
2579 This is an equation with form `g(x) = a_0 f(x) + a_1 x f'(x) + a_2 x^2 f''(x)
2580 \cdots`.
2582 These equations can be solved in a general manner, by substituting
2583 solutions of the form `x = exp(t)`, and deriving a characteristic equation
2584 of form `g(exp(t)) = b_0 f(t) + b_1 f'(t) + b_2 f''(t) \cdots` which can
2585 be then solved by nth_linear_constant_coeff_undetermined_coefficients if
2586 g(exp(t)) has finite number of linearly independent derivatives.
2588 Functions that fit this requirement are finite sums functions of the form
2589 `a x^i e^{b x} \sin(c x + d)` or `a x^i e^{b x} \cos(c x + d)`, where `i`
2590 is a non-negative integer and `a`, `b`, `c`, and `d` are constants. For
2591 example any polynomial in `x`, functions like `x^2 e^{2 x}`, `x \sin(x)`,
2592 and `e^x \cos(x)` can all be used. Products of `\sin`'s and `\cos`'s have
2593 a finite number of derivatives, because they can be expanded into `\sin(a
2594 x)` and `\cos(b x)` terms. However, SymPy currently cannot do that
2595 expansion, so you will need to manually rewrite the expression in terms of
2596 the above to use this method. So, for example, you will need to manually
2597 convert `\sin^2(x)` into `(1 + \cos(2 x))/2` to properly apply the method
2598 of undetermined coefficients on it.
2600 After replacement of x by exp(t), this method works by creating a trial function
2601 from the expression and all of its linear independent derivatives and
2602 substituting them into the original ODE. The coefficients for each term
2603 will be a system of linear equations, which are be solved for and
2604 substituted, giving the solution. If any of the trial functions are linearly
2605 dependent on the solution to the homogeneous equation, they are multiplied
2606 by sufficient `x` to make them linearly independent.
2608 Examples
2609 ========
2611 >>> from sympy import dsolve, Function, Derivative, log
2612 >>> from sympy.abc import x
2613 >>> f = Function('f')
2614 >>> eq = x**2*Derivative(f(x), x, x) - 2*x*Derivative(f(x), x) + 2*f(x) - log(x)
2615 >>> dsolve(eq, f(x),
2616 ... hint='nth_linear_euler_eq_nonhomogeneous_undetermined_coefficients').expand()
2617 Eq(f(x), C1*x + C2*x**2 + log(x)/2 + 3/4)
2619 """
2620 hint = "nth_linear_euler_eq_nonhomogeneous_undetermined_coefficients"
2621 has_integral = False
2623 def _matches(self):
2624 eq = self.ode_problem.eq_high_order_free
2625 f = self.ode_problem.func.func
2626 order = self.ode_problem.order
2627 x = self.ode_problem.sym
2628 match = self.ode_problem.get_linear_coefficients(eq, f(x), order)
2629 self.r = None
2630 does_match = False
2632 if order and match:
2633 coeff = match[order]
2634 factor = x**order / coeff
2635 self.r = {i: factor*match[i] for i in match}
2636 if self.r and all(_test_term(self.r[i], f(x), i) for i in
2637 self.r if i >= 0):
2638 if self.r[-1]:
2639 e, re = posify(self.r[-1].subs(x, exp(x)))
2640 undetcoeff = _undetermined_coefficients_match(e.subs(re), x)
2641 if undetcoeff['test']:
2642 does_match = True
2643 return does_match
2645 def _get_general_solution(self, *, simplify_flag: bool = True):
2646 f = self.ode_problem.func.func
2647 x = self.ode_problem.sym
2648 chareq, eq, symbol = S.Zero, S.Zero, Dummy('x')
2649 for i in self.r.keys():
2650 if i >= 0:
2651 chareq += (self.r[i]*diff(x**symbol, x, i)*x**-symbol).expand()
2653 for i in range(1, degree(Poly(chareq, symbol))+1):
2654 eq += chareq.coeff(symbol**i)*diff(f(x), x, i)
2656 if chareq.as_coeff_add(symbol)[0]:
2657 eq += chareq.as_coeff_add(symbol)[0]*f(x)
2658 e, re = posify(self.r[-1].subs(x, exp(x)))
2659 eq += e.subs(re)
2661 self.const_undet_instance = NthLinearConstantCoeffUndeterminedCoefficients(SingleODEProblem(eq, f(x), x))
2662 sol = self.const_undet_instance.get_general_solution(simplify = simplify_flag)[0]
2663 sol = sol.subs(x, log(x))
2664 sol = sol.subs(f(log(x)), f(x)).expand()
2666 return [sol]
2669class SecondLinearBessel(SingleODESolver):
2670 r"""
2671 Gives solution of the Bessel differential equation
2673 .. math :: x^2 \frac{d^2y}{dx^2} + x \frac{dy}{dx} y(x) + (x^2-n^2) y(x)
2675 if `n` is integer then the solution is of the form ``Eq(f(x), C0 besselj(n,x)
2676 + C1 bessely(n,x))`` as both the solutions are linearly independent else if
2677 `n` is a fraction then the solution is of the form ``Eq(f(x), C0 besselj(n,x)
2678 + C1 besselj(-n,x))`` which can also transform into ``Eq(f(x), C0 besselj(n,x)
2679 + C1 bessely(n,x))``.
2681 Examples
2682 ========
2684 >>> from sympy.abc import x
2685 >>> from sympy import Symbol
2686 >>> v = Symbol('v', positive=True)
2687 >>> from sympy import dsolve, Function
2688 >>> f = Function('f')
2689 >>> y = f(x)
2690 >>> genform = x**2*y.diff(x, 2) + x*y.diff(x) + (x**2 - v**2)*y
2691 >>> dsolve(genform)
2692 Eq(f(x), C1*besselj(v, x) + C2*bessely(v, x))
2694 References
2695 ==========
2697 https://math24.net/bessel-differential-equation.html
2699 """
2700 hint = "2nd_linear_bessel"
2701 has_integral = False
2703 def _matches(self):
2704 eq = self.ode_problem.eq_high_order_free
2705 f = self.ode_problem.func
2706 order = self.ode_problem.order
2707 x = self.ode_problem.sym
2708 df = f.diff(x)
2709 a = Wild('a', exclude=[f,df])
2710 b = Wild('b', exclude=[x, f,df])
2711 a4 = Wild('a4', exclude=[x,f,df])
2712 b4 = Wild('b4', exclude=[x,f,df])
2713 c4 = Wild('c4', exclude=[x,f,df])
2714 d4 = Wild('d4', exclude=[x,f,df])
2715 a3 = Wild('a3', exclude=[f, df, f.diff(x, 2)])
2716 b3 = Wild('b3', exclude=[f, df, f.diff(x, 2)])
2717 c3 = Wild('c3', exclude=[f, df, f.diff(x, 2)])
2718 deq = a3*(f.diff(x, 2)) + b3*df + c3*f
2719 r = collect(eq,
2720 [f.diff(x, 2), df, f]).match(deq)
2721 if order == 2 and r:
2722 if not all(r[key].is_polynomial() for key in r):
2723 n, d = eq.as_numer_denom()
2724 eq = expand(n)
2725 r = collect(eq,
2726 [f.diff(x, 2), df, f]).match(deq)
2728 if r and r[a3] != 0:
2729 # leading coeff of f(x).diff(x, 2)
2730 coeff = factor(r[a3]).match(a4*(x-b)**b4)
2732 if coeff:
2733 # if coeff[b4] = 0 means constant coefficient
2734 if coeff[b4] == 0:
2735 return False
2736 point = coeff[b]
2737 else:
2738 return False
2740 if point:
2741 r[a3] = simplify(r[a3].subs(x, x+point))
2742 r[b3] = simplify(r[b3].subs(x, x+point))
2743 r[c3] = simplify(r[c3].subs(x, x+point))
2745 # making a3 in the form of x**2
2746 r[a3] = cancel(r[a3]/(coeff[a4]*(x)**(-2+coeff[b4])))
2747 r[b3] = cancel(r[b3]/(coeff[a4]*(x)**(-2+coeff[b4])))
2748 r[c3] = cancel(r[c3]/(coeff[a4]*(x)**(-2+coeff[b4])))
2749 # checking if b3 is of form c*(x-b)
2750 coeff1 = factor(r[b3]).match(a4*(x))
2751 if coeff1 is None:
2752 return False
2753 # c3 maybe of very complex form so I am simply checking (a - b) form
2754 # if yes later I will match with the standerd form of bessel in a and b
2755 # a, b are wild variable defined above.
2756 _coeff2 = r[c3].match(a - b)
2757 if _coeff2 is None:
2758 return False
2759 # matching with standerd form for c3
2760 coeff2 = factor(_coeff2[a]).match(c4**2*(x)**(2*a4))
2761 if coeff2 is None:
2762 return False
2764 if _coeff2[b] == 0:
2765 coeff2[d4] = 0
2766 else:
2767 coeff2[d4] = factor(_coeff2[b]).match(d4**2)[d4]
2769 self.rn = {'n':coeff2[d4], 'a4':coeff2[c4], 'd4':coeff2[a4]}
2770 self.rn['c4'] = coeff1[a4]
2771 self.rn['b4'] = point
2772 return True
2773 return False
2775 def _get_general_solution(self, *, simplify_flag: bool = True):
2776 f = self.ode_problem.func.func
2777 x = self.ode_problem.sym
2778 n = self.rn['n']
2779 a4 = self.rn['a4']
2780 c4 = self.rn['c4']
2781 d4 = self.rn['d4']
2782 b4 = self.rn['b4']
2783 n = sqrt(n**2 + Rational(1, 4)*(c4 - 1)**2)
2784 (C1, C2) = self.ode_problem.get_numbered_constants(num=2)
2785 return [Eq(f(x), ((x**(Rational(1-c4,2)))*(C1*besselj(n/d4,a4*x**d4/d4)
2786 + C2*bessely(n/d4,a4*x**d4/d4))).subs(x, x-b4))]
2789class SecondLinearAiry(SingleODESolver):
2790 r"""
2791 Gives solution of the Airy differential equation
2793 .. math :: \frac{d^2y}{dx^2} + (a + b x) y(x) = 0
2795 in terms of Airy special functions airyai and airybi.
2797 Examples
2798 ========
2800 >>> from sympy import dsolve, Function
2801 >>> from sympy.abc import x
2802 >>> f = Function("f")
2803 >>> eq = f(x).diff(x, 2) - x*f(x)
2804 >>> dsolve(eq)
2805 Eq(f(x), C1*airyai(x) + C2*airybi(x))
2806 """
2807 hint = "2nd_linear_airy"
2808 has_integral = False
2810 def _matches(self):
2811 eq = self.ode_problem.eq_high_order_free
2812 f = self.ode_problem.func
2813 order = self.ode_problem.order
2814 x = self.ode_problem.sym
2815 df = f.diff(x)
2816 a4 = Wild('a4', exclude=[x,f,df])
2817 b4 = Wild('b4', exclude=[x,f,df])
2818 match = self.ode_problem.get_linear_coefficients(eq, f, order)
2819 does_match = False
2820 if order == 2 and match and match[2] != 0:
2821 if match[1].is_zero:
2822 self.rn = cancel(match[0]/match[2]).match(a4+b4*x)
2823 if self.rn and self.rn[b4] != 0:
2824 self.rn = {'b':self.rn[a4],'m':self.rn[b4]}
2825 does_match = True
2826 return does_match
2828 def _get_general_solution(self, *, simplify_flag: bool = True):
2829 f = self.ode_problem.func.func
2830 x = self.ode_problem.sym
2831 (C1, C2) = self.ode_problem.get_numbered_constants(num=2)
2832 b = self.rn['b']
2833 m = self.rn['m']
2834 if m.is_positive:
2835 arg = - b/cbrt(m)**2 - cbrt(m)*x
2836 elif m.is_negative:
2837 arg = - b/cbrt(-m)**2 + cbrt(-m)*x
2838 else:
2839 arg = - b/cbrt(-m)**2 + cbrt(-m)*x
2841 return [Eq(f(x), C1*airyai(arg) + C2*airybi(arg))]
2844class LieGroup(SingleODESolver):
2845 r"""
2846 This hint implements the Lie group method of solving first order differential
2847 equations. The aim is to convert the given differential equation from the
2848 given coordinate system into another coordinate system where it becomes
2849 invariant under the one-parameter Lie group of translations. The converted
2850 ODE can be easily solved by quadrature. It makes use of the
2851 :py:meth:`sympy.solvers.ode.infinitesimals` function which returns the
2852 infinitesimals of the transformation.
2854 The coordinates `r` and `s` can be found by solving the following Partial
2855 Differential Equations.
2857 .. math :: \xi\frac{\partial r}{\partial x} + \eta\frac{\partial r}{\partial y}
2858 = 0
2860 .. math :: \xi\frac{\partial s}{\partial x} + \eta\frac{\partial s}{\partial y}
2861 = 1
2863 The differential equation becomes separable in the new coordinate system
2865 .. math :: \frac{ds}{dr} = \frac{\frac{\partial s}{\partial x} +
2866 h(x, y)\frac{\partial s}{\partial y}}{
2867 \frac{\partial r}{\partial x} + h(x, y)\frac{\partial r}{\partial y}}
2869 After finding the solution by integration, it is then converted back to the original
2870 coordinate system by substituting `r` and `s` in terms of `x` and `y` again.
2872 Examples
2873 ========
2875 >>> from sympy import Function, dsolve, exp, pprint
2876 >>> from sympy.abc import x
2877 >>> f = Function('f')
2878 >>> pprint(dsolve(f(x).diff(x) + 2*x*f(x) - x*exp(-x**2), f(x),
2879 ... hint='lie_group'))
2880 / 2\ 2
2881 | x | -x
2882 f(x) = |C1 + --|*e
2883 \ 2 /
2886 References
2887 ==========
2889 - Solving differential equations by Symmetry Groups,
2890 John Starrett, pp. 1 - pp. 14
2892 """
2893 hint = "lie_group"
2894 has_integral = False
2896 def _has_additional_params(self):
2897 return 'xi' in self.ode_problem.params and 'eta' in self.ode_problem.params
2899 def _matches(self):
2900 eq = self.ode_problem.eq
2901 f = self.ode_problem.func.func
2902 order = self.ode_problem.order
2903 x = self.ode_problem.sym
2904 df = f(x).diff(x)
2905 y = Dummy('y')
2906 d = Wild('d', exclude=[df, f(x).diff(x, 2)])
2907 e = Wild('e', exclude=[df])
2908 does_match = False
2909 if self._has_additional_params() and order == 1:
2910 xi = self.ode_problem.params['xi']
2911 eta = self.ode_problem.params['eta']
2912 self.r3 = {'xi': xi, 'eta': eta}
2913 r = collect(eq, df, exact=True).match(d + e * df)
2914 if r:
2915 r['d'] = d
2916 r['e'] = e
2917 r['y'] = y
2918 r[d] = r[d].subs(f(x), y)
2919 r[e] = r[e].subs(f(x), y)
2920 self.r3.update(r)
2921 does_match = True
2922 return does_match
2924 def _get_general_solution(self, *, simplify_flag: bool = True):
2925 eq = self.ode_problem.eq
2926 x = self.ode_problem.sym
2927 func = self.ode_problem.func
2928 order = self.ode_problem.order
2929 df = func.diff(x)
2931 try:
2932 eqsol = solve(eq, df)
2933 except NotImplementedError:
2934 eqsol = []
2936 desols = []
2937 for s in eqsol:
2938 sol = _ode_lie_group(s, func, order, match=self.r3)
2939 if sol:
2940 desols.extend(sol)
2942 if desols == []:
2943 raise NotImplementedError("The given ODE " + str(eq) + " cannot be solved by"
2944 + " the lie group method")
2945 return desols
2948solver_map = {
2949 'factorable': Factorable,
2950 'nth_linear_constant_coeff_homogeneous': NthLinearConstantCoeffHomogeneous,
2951 'nth_linear_euler_eq_homogeneous': NthLinearEulerEqHomogeneous,
2952 'nth_linear_constant_coeff_undetermined_coefficients': NthLinearConstantCoeffUndeterminedCoefficients,
2953 'nth_linear_euler_eq_nonhomogeneous_undetermined_coefficients': NthLinearEulerEqNonhomogeneousUndeterminedCoefficients,
2954 'separable': Separable,
2955 '1st_exact': FirstExact,
2956 '1st_linear': FirstLinear,
2957 'Bernoulli': Bernoulli,
2958 'Riccati_special_minus2': RiccatiSpecial,
2959 '1st_rational_riccati': RationalRiccati,
2960 '1st_homogeneous_coeff_best': HomogeneousCoeffBest,
2961 '1st_homogeneous_coeff_subs_indep_div_dep': HomogeneousCoeffSubsIndepDivDep,
2962 '1st_homogeneous_coeff_subs_dep_div_indep': HomogeneousCoeffSubsDepDivIndep,
2963 'almost_linear': AlmostLinear,
2964 'linear_coefficients': LinearCoefficients,
2965 'separable_reduced': SeparableReduced,
2966 'nth_linear_constant_coeff_variation_of_parameters': NthLinearConstantCoeffVariationOfParameters,
2967 'nth_linear_euler_eq_nonhomogeneous_variation_of_parameters': NthLinearEulerEqNonhomogeneousVariationOfParameters,
2968 'Liouville': Liouville,
2969 '2nd_linear_airy': SecondLinearAiry,
2970 '2nd_linear_bessel': SecondLinearBessel,
2971 '2nd_hypergeometric': SecondHypergeometric,
2972 'nth_order_reducible': NthOrderReducible,
2973 '2nd_nonlinear_autonomous_conserved': SecondNonlinearAutonomousConserved,
2974 'nth_algebraic': NthAlgebraic,
2975 'lie_group': LieGroup,
2976 }
2978# Avoid circular import:
2979from .ode import dsolve, ode_sol_simplicity, odesimp, homogeneous_order