Server : Apache/2.4.41 (Ubuntu) System : Linux absol.cf 5.4.0-198-generic #218-Ubuntu SMP Fri Sep 27 20:18:53 UTC 2024 x86_64 User : www-data ( 33) PHP Version : 7.4.33 Disable Function : pcntl_alarm,pcntl_fork,pcntl_waitpid,pcntl_wait,pcntl_wifexited,pcntl_wifstopped,pcntl_wifsignaled,pcntl_wifcontinued,pcntl_wexitstatus,pcntl_wtermsig,pcntl_wstopsig,pcntl_signal,pcntl_signal_get_handler,pcntl_signal_dispatch,pcntl_get_last_error,pcntl_strerror,pcntl_sigprocmask,pcntl_sigwaitinfo,pcntl_sigtimedwait,pcntl_exec,pcntl_getpriority,pcntl_setpriority,pcntl_async_signals,pcntl_unshare, Directory : /usr/local/lib/python3.6/dist-packages/sympy/solvers/ode/
Upload File :
Current File : //usr/local/lib/python3.6/dist-packages/sympy/solvers/ode/single.py
#
# This is the module for ODE solver classes for single ODEs.
#
import typing
if typing.TYPE_CHECKING:
from typing import ClassVar
from typing import Dict, Type
from typing import Iterator, List, Optional
from sympy.core import S
from sympy.core.exprtools import factor_terms
from sympy.core.expr import Expr
from sympy.core.function import AppliedUndef, Derivative, Function, expand
from sympy.core.numbers import Float
from sympy.core.relational import Equality, Eq
from sympy.core.symbol import Symbol, Dummy, Wild
from sympy.core.mul import Mul
from sympy.functions import exp, sqrt, tan, log
from sympy.integrals import Integral
from sympy.polys.polytools import cancel, factor
from sympy.simplify.simplify import simplify
from sympy.simplify.radsimp import fraction
from sympy.utilities import numbered_symbols
from sympy.solvers.solvers import solve
from sympy.solvers.deutils import ode_order, _preprocess
class ODEMatchError(NotImplementedError):
"""Raised if a SingleODESolver is asked to solve an ODE it does not match"""
pass
def cached_property(func):
'''Decorator to cache property method'''
attrname = '_' + func.__name__
def propfunc(self):
val = getattr(self, attrname, None)
if val is None:
val = func(self)
setattr(self, attrname, val)
return val
return property(propfunc)
class SingleODEProblem:
"""Represents an ordinary differential equation (ODE)
This class is used internally in the by dsolve and related
functions/classes so that properties of an ODE can be computed
efficiently.
Examples
========
This class is used internally by dsolve. To instantiate an instance
directly first define an ODE problem:
>>> from sympy import Function, Symbol
>>> x = Symbol('x')
>>> f = Function('f')
>>> eq = f(x).diff(x, 2)
Now you can create a SingleODEProblem instance and query its properties:
>>> from sympy.solvers.ode.single import SingleODEProblem
>>> problem = SingleODEProblem(f(x).diff(x), f(x), x)
>>> problem.eq
Derivative(f(x), x)
>>> problem.func
f(x)
>>> problem.sym
x
"""
# Instance attributes:
eq = None # type: Expr
func = None # type: AppliedUndef
sym = None # type: Symbol
_order = None # type: int
_eq_expanded = None # type: Expr
_eq_preprocessed = None # type: Expr
def __init__(self, eq, func, sym, prep=True):
assert isinstance(eq, Expr)
assert isinstance(func, AppliedUndef)
assert isinstance(sym, Symbol)
assert isinstance(prep, bool)
self.eq = eq
self.func = func
self.sym = sym
self.prep = prep
@cached_property
def order(self) -> int:
return ode_order(self.eq, self.func)
@cached_property
def eq_preprocessed(self) -> Expr:
return self._get_eq_preprocessed()
@cached_property
def eq_expanded(self) -> Expr:
return expand(self.eq_preprocessed)
def _get_eq_preprocessed(self) -> Expr:
if self.prep:
process_eq, process_func = _preprocess(self.eq, self.func)
if process_func != self.func:
raise ValueError
else:
process_eq = self.eq
return process_eq
def get_numbered_constants(self, num=1, start=1, prefix='C') -> List[Symbol]:
"""
Returns a list of constants that do not occur
in eq already.
"""
ncs = self.iter_numbered_constants(start, prefix)
Cs = [next(ncs) for i in range(num)]
return Cs
def iter_numbered_constants(self, start=1, prefix='C') -> Iterator[Symbol]:
"""
Returns an iterator of constants that do not occur
in eq already.
"""
atom_set = self.eq.free_symbols
func_set = self.eq.atoms(Function)
if func_set:
atom_set |= {Symbol(str(f.func)) for f in func_set}
return numbered_symbols(start=start, prefix=prefix, exclude=atom_set)
# TODO: Add methods that can be used by many ODE solvers:
# order
# is_linear()
# get_linear_coefficients()
# eq_prepared (the ODE in prepared form)
class SingleODESolver:
"""
Base class for Single ODE solvers.
Subclasses should implement the _matches and _get_general_solution
methods. This class is not intended to be instantiated directly but its
subclasses are as part of dsolve.
Examples
========
You can use a subclass of SingleODEProblem to solve a particular type of
ODE. We first define a particular ODE problem:
>>> from sympy import Function, Symbol
>>> x = Symbol('x')
>>> f = Function('f')
>>> eq = f(x).diff(x, 2)
Now we solve this problem using the NthAlgebraic solver which is a
subclass of SingleODESolver:
>>> from sympy.solvers.ode.single import NthAlgebraic, SingleODEProblem
>>> problem = SingleODEProblem(eq, f(x), x)
>>> solver = NthAlgebraic(problem)
>>> solver.get_general_solution()
[Eq(f(x), _C*x + _C)]
The normal way to solve an ODE is to use dsolve (which would use
NthAlgebraic and other solvers internally). When using dsolve a number of
other things are done such as evaluating integrals, simplifying the
solution and renumbering the constants:
>>> from sympy import dsolve
>>> dsolve(eq, hint='nth_algebraic')
Eq(f(x), C1 + C2*x)
"""
# Subclasses should store the hint name (the argument to dsolve) in this
# attribute
hint = None # type: ClassVar[str]
# Subclasses should define this to indicate if they support an _Integral
# hint.
has_integral = None # type: ClassVar[bool]
# The ODE to be solved
ode_problem = None # type: SingleODEProblem
# Cache whether or not the equation has matched the method
_matched = None # type: Optional[bool]
# Subclasses should store in this attribute the list of order(s) of ODE
# that subclass can solve or leave it to None if not specific to any order
order = None # type: Optional[list]
def __init__(self, ode_problem):
self.ode_problem = ode_problem
def matches(self) -> bool:
if self.order is not None and self.ode_problem.order not in self.order:
self._matched = False
return self._matched
if self._matched is None:
self._matched = self._matches()
return self._matched
def get_general_solution(self, *, simplify: bool = True) -> List[Equality]:
if not self.matches():
msg = "%s solver can not solve:\n%s"
raise ODEMatchError(msg % (self.hint, self.ode_problem.eq))
return self._get_general_solution()
def _matches(self) -> bool:
msg = "Subclasses of SingleODESolver should implement matches."
raise NotImplementedError(msg)
def _get_general_solution(self, *, simplify: bool = True) -> List[Equality]:
msg = "Subclasses of SingleODESolver should implement get_general_solution."
raise NotImplementedError(msg)
class SinglePatternODESolver(SingleODESolver):
'''Superclass for ODE solvers based on pattern matching'''
def wilds(self):
prob = self.ode_problem
f = prob.func.func
x = prob.sym
order = prob.order
return self._wilds(f, x, order)
def wilds_match(self):
match = self._wilds_match
return [match.get(w, S.Zero) for w in self.wilds()]
def _matches(self):
eq = self.ode_problem.eq_expanded
f = self.ode_problem.func.func
x = self.ode_problem.sym
order = self.ode_problem.order
df = f(x).diff(x)
if order != 1:
return False
pattern = self._equation(f(x), x, 1)
if not pattern.coeff(df).has(Wild):
eq = expand(eq / eq.coeff(df))
eq = eq.collect(f(x), func = cancel)
self._wilds_match = match = eq.match(pattern)
if match is not None:
return self._verify(f(x))
return False
def _verify(self, fx) -> bool:
return True
def _wilds(self, f, x, order):
msg = "Subclasses of SingleODESolver should implement _wilds"
raise NotImplementedError(msg)
def _equation(self, fx, x, order):
msg = "Subclasses of SingleODESolver should implement _equation"
raise NotImplementedError(msg)
class NthAlgebraic(SingleODESolver):
r"""
Solves an `n`\th order ordinary differential equation using algebra and
integrals.
There is no general form for the kind of equation that this can solve. The
the equation is solved algebraically treating differentiation as an
invertible algebraic function.
Examples
========
>>> from sympy import Function, dsolve, Eq
>>> from sympy.abc import x
>>> f = Function('f')
>>> eq = Eq(f(x) * (f(x).diff(x)**2 - 1), 0)
>>> dsolve(eq, f(x), hint='nth_algebraic')
[Eq(f(x), 0), Eq(f(x), C1 - x), Eq(f(x), C1 + x)]
Note that this solver can return algebraic solutions that do not have any
integration constants (f(x) = 0 in the above example).
"""
hint = 'nth_algebraic'
has_integral = True # nth_algebraic_Integral hint
def _matches(self):
r"""
Matches any differential equation that nth_algebraic can solve. Uses
`sympy.solve` but teaches it how to integrate derivatives.
This involves calling `sympy.solve` and does most of the work of finding a
solution (apart from evaluating the integrals).
"""
eq = self.ode_problem.eq
func = self.ode_problem.func
var = self.ode_problem.sym
# Derivative that solve can handle:
diffx = self._get_diffx(var)
# Replace derivatives wrt the independent variable with diffx
def replace(eq, var):
def expand_diffx(*args):
differand, diffs = args[0], args[1:]
toreplace = differand
for v, n in diffs:
for _ in range(n):
if v == var:
toreplace = diffx(toreplace)
else:
toreplace = Derivative(toreplace, v)
return toreplace
return eq.replace(Derivative, expand_diffx)
# Restore derivatives in solution afterwards
def unreplace(eq, var):
return eq.replace(diffx, lambda e: Derivative(e, var))
subs_eqn = replace(eq, var)
try:
# turn off simplification to protect Integrals that have
# _t instead of fx in them and would otherwise factor
# as t_*Integral(1, x)
solns = solve(subs_eqn, func, simplify=False)
except NotImplementedError:
solns = []
solns = [simplify(unreplace(soln, var)) for soln in solns]
solns = [Equality(func, soln) for soln in solns]
self.solutions = solns
return len(solns) != 0
def _get_general_solution(self, *, simplify: bool = True):
return self.solutions
# This needs to produce an invertible function but the inverse depends
# which variable we are integrating with respect to. Since the class can
# be stored in cached results we need to ensure that we always get the
# same class back for each particular integration variable so we store these
# classes in a global dict:
_diffx_stored = {} # type: Dict[Symbol, Type[Function]]
@staticmethod
def _get_diffx(var):
diffcls = NthAlgebraic._diffx_stored.get(var, None)
if diffcls is None:
# A class that behaves like Derivative wrt var but is "invertible".
class diffx(Function):
def inverse(self):
# don't use integrate here because fx has been replaced by _t
# in the equation; integrals will not be correct while solve
# is at work.
return lambda expr: Integral(expr, var) + Dummy('C')
diffcls = NthAlgebraic._diffx_stored.setdefault(var, diffx)
return diffcls
class FirstLinear(SinglePatternODESolver):
r"""
Solves 1st order linear differential equations.
These are differential equations of the form
.. math:: dy/dx + P(x) y = Q(x)\text{.}
These kinds of differential equations can be solved in a general way. The
integrating factor `e^{\int P(x) \,dx}` will turn the equation into a
separable equation. The general solution is::
>>> from sympy import Function, dsolve, Eq, pprint, diff, sin
>>> from sympy.abc import x
>>> f, P, Q = map(Function, ['f', 'P', 'Q'])
>>> genform = Eq(f(x).diff(x) + P(x)*f(x), Q(x))
>>> pprint(genform)
d
P(x)*f(x) + --(f(x)) = Q(x)
dx
>>> pprint(dsolve(genform, f(x), hint='1st_linear_Integral'))
/ / \
| | |
| | / | /
| | | | |
| | | P(x) dx | - | P(x) dx
| | | | |
| | / | /
f(x) = |C1 + | Q(x)*e dx|*e
| | |
\ / /
Examples
========
>>> f = Function('f')
>>> pprint(dsolve(Eq(x*diff(f(x), x) - f(x), x**2*sin(x)),
... f(x), '1st_linear'))
f(x) = x*(C1 - cos(x))
References
==========
- https://en.wikipedia.org/wiki/Linear_differential_equation#First_order_equation
- M. Tenenbaum & H. Pollard, "Ordinary Differential Equations",
Dover 1963, pp. 92
# indirect doctest
"""
hint = '1st_linear'
has_integral = True
order = [1]
def _wilds(self, f, x, order):
P = Wild('P', exclude=[f(x)])
Q = Wild('Q', exclude=[f(x), f(x).diff(x)])
return P, Q
def _equation(self, fx, x, order):
P, Q = self.wilds()
return fx.diff(x) + P*fx - Q
def _get_general_solution(self, *, simplify: bool = True):
P, Q = self.wilds_match()
fx = self.ode_problem.func
x = self.ode_problem.sym
(C1,) = self.ode_problem.get_numbered_constants(num=1)
gensol = Eq(fx, (((C1 + Integral(Q*exp(Integral(P, x)),x))
* exp(-Integral(P, x)))))
return [gensol]
class AlmostLinear(SinglePatternODESolver):
r"""
Solves an almost-linear differential equation.
The general form of an almost linear differential equation is
.. math:: a(x) g'(f(x)) f'(x) + b(x) g(f(x)) + c(x)
Here `f(x)` is the function to be solved for (the dependent variable).
The substitution `g(f(x)) = u(x)` leads to a linear differential equation
for `u(x)` of the form `a(x) u' + b(x) u + c(x) = 0`. This can be solved
for `u(x)` by the `first_linear` hint and then `f(x)` is found by solving
`g(f(x)) = u(x)`.
See Also
========
:meth:`sympy.solvers.ode.single.FirstLinear`
Examples
========
>>> from sympy import Function, pprint, sin, cos
>>> from sympy.solvers.ode import dsolve
>>> from sympy.abc import x
>>> f = Function('f')
>>> d = f(x).diff(x)
>>> eq = x*d + x*f(x) + 1
>>> dsolve(eq, f(x), hint='almost_linear')
Eq(f(x), (C1 - Ei(x))*exp(-x))
>>> pprint(dsolve(eq, f(x), hint='almost_linear'))
-x
f(x) = (C1 - Ei(x))*e
>>> example = cos(f(x))*f(x).diff(x) + sin(f(x)) + 1
>>> pprint(example)
d
sin(f(x)) + cos(f(x))*--(f(x)) + 1
dx
>>> pprint(dsolve(example, f(x), hint='almost_linear'))
/ -x \ / -x \
[f(x) = pi - asin\C1*e - 1/, f(x) = asin\C1*e - 1/]
References
==========
- Joel Moses, "Symbolic Integration - The Stormy Decade", Communications
of the ACM, Volume 14, Number 8, August 1971, pp. 558
"""
hint = "almost_linear"
has_integral = True
order = [1]
def _wilds(self, f, x, order):
P = Wild('P', exclude=[f(x).diff(x)])
Q = Wild('Q', exclude=[f(x).diff(x)])
return P, Q
def _equation(self, fx, x, order):
P, Q = self.wilds()
return P*fx.diff(x) + Q
def _verify(self, fx):
a, b = self.wilds_match()
c, b = b.as_independent(fx) if b.is_Add else (S.Zero, b)
# a, b and c are the function a(x), b(x) and c(x) respectively.
# c(x) is obtained by separating out b as terms with and without fx i.e, l(y)
# The following conditions checks if the given equation is an almost-linear differential equation using the fact that
# a(x)*(l(y))' / l(y)' is independent of l(y)
if b.diff(fx) != 0 and not simplify(b.diff(fx)/a).has(fx):
self.ly = factor_terms(b).as_independent(fx, as_Add=False)[1] # Gives the term containing fx i.e., l(y)
self.ax = a / self.ly.diff(fx)
self.cx = -c # cx is taken as -c(x) to simplify expression in the solution integral
self.bx = factor_terms(b) / self.ly
return True
return False
def _get_general_solution(self, *, simplify: bool = True):
x = self.ode_problem.sym
(C1,) = self.ode_problem.get_numbered_constants(num=1)
gensol = Eq(self.ly, (((C1 + Integral((self.cx/self.ax)*exp(Integral(self.bx/self.ax, x)),x))
* exp(-Integral(self.bx/self.ax, x)))))
return [gensol]
class Bernoulli(SinglePatternODESolver):
r"""
Solves Bernoulli differential equations.
These are equations of the form
.. math:: dy/dx + P(x) y = Q(x) y^n\text{, }n \ne 1`\text{.}
The substitution `w = 1/y^{1-n}` will transform an equation of this form
into one that is linear (see the docstring of
:py:meth:`~sympy.solvers.ode.single.FirstLinear`). The general solution is::
>>> from sympy import Function, dsolve, Eq, pprint
>>> from sympy.abc import x, n
>>> f, P, Q = map(Function, ['f', 'P', 'Q'])
>>> genform = Eq(f(x).diff(x) + P(x)*f(x), Q(x)*f(x)**n)
>>> pprint(genform)
d n
P(x)*f(x) + --(f(x)) = Q(x)*f (x)
dx
>>> pprint(dsolve(genform, f(x), hint='Bernoulli_Integral'), num_columns=110)
-1
-----
n - 1
// / / \ \
|| | | | |
|| | / | / | / |
|| | | | | | | |
|| | (1 - n)* | P(x) dx | (1 - n)* | P(x) dx | (n - 1)* | P(x) dx|
|| | | | | | | |
|| | / | / | / |
f(x) = ||C1 - n* | Q(x)*e dx + | Q(x)*e dx|*e |
|| | | | |
\\ / / / /
Note that the equation is separable when `n = 1` (see the docstring of
:py:meth:`~sympy.solvers.ode.ode.ode_separable`).
>>> pprint(dsolve(Eq(f(x).diff(x) + P(x)*f(x), Q(x)*f(x)), f(x),
... hint='separable_Integral'))
f(x)
/
| /
| 1 |
| - dy = C1 + | (-P(x) + Q(x)) dx
| y |
| /
/
Examples
========
>>> from sympy import Function, dsolve, Eq, pprint, log
>>> from sympy.abc import x
>>> f = Function('f')
>>> pprint(dsolve(Eq(x*f(x).diff(x) + f(x), log(x)*f(x)**2),
... f(x), hint='Bernoulli'))
1
f(x) = -----------------
C1*x + log(x) + 1
References
==========
- https://en.wikipedia.org/wiki/Bernoulli_differential_equation
- M. Tenenbaum & H. Pollard, "Ordinary Differential Equations",
Dover 1963, pp. 95
# indirect doctest
"""
hint = "Bernoulli"
has_integral = True
order = [1]
def _wilds(self, f, x, order):
P = Wild('P', exclude=[f(x)])
Q = Wild('Q', exclude=[f(x)])
n = Wild('n', exclude=[x, f(x), f(x).diff(x)])
return P, Q, n
def _equation(self, fx, x, order):
P, Q, n = self.wilds()
return fx.diff(x) + P*fx - Q*fx**n
def _get_general_solution(self, *, simplify: bool = True):
P, Q, n = self.wilds_match()
fx = self.ode_problem.func
x = self.ode_problem.sym
(C1,) = self.ode_problem.get_numbered_constants(num=1)
if n==1:
gensol = Eq(log(fx), (
(C1 + Integral((-P + Q),x)
)))
else:
gensol = Eq(fx**(1-n), (
(C1 - (n - 1) * Integral(Q*exp(-n*Integral(P, x))
* exp(Integral(P, x)), x)
) * exp(-(1 - n)*Integral(P, x)))
)
return [gensol]
class Factorable(SingleODESolver):
r"""
Solves equations having a solvable factor.
This function is used to solve the equation having factors. Factors may be of type algebraic or ode. It
will try to solve each factor independently. Factors will be solved by calling dsolve. We will return the
list of solutions.
Examples
========
>>> from sympy import Function, dsolve, pprint
>>> from sympy.abc import x
>>> f = Function('f')
>>> eq = (f(x)**2-4)*(f(x).diff(x)+f(x))
>>> pprint(dsolve(eq, f(x)))
-x
[f(x) = 2, f(x) = -2, f(x) = C1*e ]
"""
hint = "factorable"
has_integral = False
def _matches(self):
eq = self.ode_problem.eq
f = self.ode_problem.func.func
x = self.ode_problem.sym
order =self.ode_problem.order
df = f(x).diff(x)
self.eqs = []
eq = eq.collect(f(x), func = cancel)
eq = fraction(factor(eq))[0]
factors = Mul.make_args(factor(eq))
roots = [fac.as_base_exp() for fac in factors if len(fac.args)!=0]
if len(roots)>1 or roots[0][1]>1:
for base,expo in roots:
if base.has(f(x)):
self.eqs.append(base)
if len(self.eqs)>0:
return True
roots = solve(eq, df)
if len(roots)>0:
self.eqs = [(df - root) for root in roots]
if len(self.eqs)==1:
if order>1:
return False
if self.eqs[0].has(Float):
return False
return fraction(factor(self.eqs[0]))[0]-eq!=0
return True
return False
def _get_general_solution(self, *, simplify: bool = True):
func = self.ode_problem.func.func
x = self.ode_problem.sym
eqns = self.eqs
sols = []
for eq in eqns:
try:
sol = dsolve(eq, func(x))
except NotImplementedError:
continue
else:
if isinstance(sol, list):
sols.extend(sol)
else:
sols.append(sol)
if sols == []:
raise NotImplementedError("The given ODE " + str(eq) + " cannot be solved by"
+ " the factorable group method")
return sols
class RiccatiSpecial(SinglePatternODESolver):
r"""
The general Riccati equation has the form
.. math:: dy/dx = f(x) y^2 + g(x) y + h(x)\text{.}
While it does not have a general solution [1], the "special" form, `dy/dx
= a y^2 - b x^c`, does have solutions in many cases [2]. This routine
returns a solution for `a(dy/dx) = b y^2 + c y/x + d/x^2` that is obtained
by using a suitable change of variables to reduce it to the special form
and is valid when neither `a` nor `b` are zero and either `c` or `d` is
zero.
>>> from sympy.abc import x, a, b, c, d
>>> from sympy.solvers.ode import dsolve, checkodesol
>>> from sympy import pprint, Function
>>> f = Function('f')
>>> y = f(x)
>>> genform = a*y.diff(x) - (b*y**2 + c*y/x + d/x**2)
>>> sol = dsolve(genform, y)
>>> pprint(sol, wrap_line=False)
/ / __________________ \\
| __________________ | / 2 ||
| / 2 | \/ 4*b*d - (a + c) *log(x)||
-|a + c - \/ 4*b*d - (a + c) *tan|C1 + ----------------------------||
\ \ 2*a //
f(x) = ------------------------------------------------------------------------
2*b*x
>>> checkodesol(genform, sol, order=1)[0]
True
References
==========
1. http://www.maplesoft.com/support/help/Maple/view.aspx?path=odeadvisor/Riccati
2. http://eqworld.ipmnet.ru/en/solutions/ode/ode0106.pdf -
http://eqworld.ipmnet.ru/en/solutions/ode/ode0123.pdf
"""
hint = "Riccati_special_minus2"
has_integral = False
order = [1]
def _wilds(self, f, x, order):
a = Wild('a', exclude=[x, f(x), f(x).diff(x), 0])
b = Wild('b', exclude=[x, f(x), f(x).diff(x), 0])
c = Wild('c', exclude=[x, f(x), f(x).diff(x)])
d = Wild('d', exclude=[x, f(x), f(x).diff(x)])
return a, b, c, d
def _equation(self, fx, x, order):
a, b, c, d = self.wilds()
return a*fx.diff(x) + b*fx**2 + c*fx/x + d/x**2
def _get_general_solution(self, *, simplify: bool = True):
a, b, c, d = self.wilds_match()
fx = self.ode_problem.func
x = self.ode_problem.sym
(C1,) = self.ode_problem.get_numbered_constants(num=1)
mu = sqrt(4*d*b - (a - c)**2)
gensol = Eq(fx, (a - c - mu*tan(mu/(2*a)*log(x) + C1))/(2*b*x))
return [gensol]
# Avoid circular import:
from .ode import dsolve
Server : Apache/2.4.41 (Ubuntu)