Coverage for /usr/lib/python3/dist-packages/sympy/solvers/ode/subscheck.py: 9%
179 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-14 15:55 +0200
1from sympy.core import S, Pow
2from sympy.core.function import (Derivative, AppliedUndef, diff)
3from sympy.core.relational import Equality, Eq
4from sympy.core.symbol import Dummy
5from sympy.core.sympify import sympify
7from sympy.logic.boolalg import BooleanAtom
8from sympy.functions import exp
9from sympy.series import Order
10from sympy.simplify.simplify import simplify, posify, besselsimp
11from sympy.simplify.trigsimp import trigsimp
12from sympy.simplify.sqrtdenest import sqrtdenest
13from sympy.solvers import solve
14from sympy.solvers.deutils import _preprocess, ode_order
15from sympy.utilities.iterables import iterable, is_sequence
18def sub_func_doit(eq, func, new):
19 r"""
20 When replacing the func with something else, we usually want the
21 derivative evaluated, so this function helps in making that happen.
23 Examples
24 ========
26 >>> from sympy import Derivative, symbols, Function
27 >>> from sympy.solvers.ode.subscheck import sub_func_doit
28 >>> x, z = symbols('x, z')
29 >>> y = Function('y')
31 >>> sub_func_doit(3*Derivative(y(x), x) - 1, y(x), x)
32 2
34 >>> sub_func_doit(x*Derivative(y(x), x) - y(x)**2 + y(x), y(x),
35 ... 1/(x*(z + 1/x)))
36 x*(-1/(x**2*(z + 1/x)) + 1/(x**3*(z + 1/x)**2)) + 1/(x*(z + 1/x))
37 ...- 1/(x**2*(z + 1/x)**2)
38 """
39 reps= {func: new}
40 for d in eq.atoms(Derivative):
41 if d.expr == func:
42 reps[d] = new.diff(*d.variable_count)
43 else:
44 reps[d] = d.xreplace({func: new}).doit(deep=False)
45 return eq.xreplace(reps)
48def checkodesol(ode, sol, func=None, order='auto', solve_for_func=True):
49 r"""
50 Substitutes ``sol`` into ``ode`` and checks that the result is ``0``.
52 This works when ``func`` is one function, like `f(x)` or a list of
53 functions like `[f(x), g(x)]` when `ode` is a system of ODEs. ``sol`` can
54 be a single solution or a list of solutions. Each solution may be an
55 :py:class:`~sympy.core.relational.Equality` that the solution satisfies,
56 e.g. ``Eq(f(x), C1), Eq(f(x) + C1, 0)``; or simply an
57 :py:class:`~sympy.core.expr.Expr`, e.g. ``f(x) - C1``. In most cases it
58 will not be necessary to explicitly identify the function, but if the
59 function cannot be inferred from the original equation it can be supplied
60 through the ``func`` argument.
62 If a sequence of solutions is passed, the same sort of container will be
63 used to return the result for each solution.
65 It tries the following methods, in order, until it finds zero equivalence:
67 1. Substitute the solution for `f` in the original equation. This only
68 works if ``ode`` is solved for `f`. It will attempt to solve it first
69 unless ``solve_for_func == False``.
70 2. Take `n` derivatives of the solution, where `n` is the order of
71 ``ode``, and check to see if that is equal to the solution. This only
72 works on exact ODEs.
73 3. Take the 1st, 2nd, ..., `n`\th derivatives of the solution, each time
74 solving for the derivative of `f` of that order (this will always be
75 possible because `f` is a linear operator). Then back substitute each
76 derivative into ``ode`` in reverse order.
78 This function returns a tuple. The first item in the tuple is ``True`` if
79 the substitution results in ``0``, and ``False`` otherwise. The second
80 item in the tuple is what the substitution results in. It should always
81 be ``0`` if the first item is ``True``. Sometimes this function will
82 return ``False`` even when an expression is identically equal to ``0``.
83 This happens when :py:meth:`~sympy.simplify.simplify.simplify` does not
84 reduce the expression to ``0``. If an expression returned by this
85 function vanishes identically, then ``sol`` really is a solution to
86 the ``ode``.
88 If this function seems to hang, it is probably because of a hard
89 simplification.
91 To use this function to test, test the first item of the tuple.
93 Examples
94 ========
96 >>> from sympy import (Eq, Function, checkodesol, symbols,
97 ... Derivative, exp)
98 >>> x, C1, C2 = symbols('x,C1,C2')
99 >>> f, g = symbols('f g', cls=Function)
100 >>> checkodesol(f(x).diff(x), Eq(f(x), C1))
101 (True, 0)
102 >>> assert checkodesol(f(x).diff(x), C1)[0]
103 >>> assert not checkodesol(f(x).diff(x), x)[0]
104 >>> checkodesol(f(x).diff(x, 2), x**2)
105 (False, 2)
107 >>> eqs = [Eq(Derivative(f(x), x), f(x)), Eq(Derivative(g(x), x), g(x))]
108 >>> sol = [Eq(f(x), C1*exp(x)), Eq(g(x), C2*exp(x))]
109 >>> checkodesol(eqs, sol)
110 (True, [0, 0])
112 """
113 if iterable(ode):
114 return checksysodesol(ode, sol, func=func)
116 if not isinstance(ode, Equality):
117 ode = Eq(ode, 0)
118 if func is None:
119 try:
120 _, func = _preprocess(ode.lhs)
121 except ValueError:
122 funcs = [s.atoms(AppliedUndef) for s in (
123 sol if is_sequence(sol, set) else [sol])]
124 funcs = set().union(*funcs)
125 if len(funcs) != 1:
126 raise ValueError(
127 'must pass func arg to checkodesol for this case.')
128 func = funcs.pop()
129 if not isinstance(func, AppliedUndef) or len(func.args) != 1:
130 raise ValueError(
131 "func must be a function of one variable, not %s" % func)
132 if is_sequence(sol, set):
133 return type(sol)([checkodesol(ode, i, order=order, solve_for_func=solve_for_func) for i in sol])
135 if not isinstance(sol, Equality):
136 sol = Eq(func, sol)
137 elif sol.rhs == func:
138 sol = sol.reversed
140 if order == 'auto':
141 order = ode_order(ode, func)
142 solved = sol.lhs == func and not sol.rhs.has(func)
143 if solve_for_func and not solved:
144 rhs = solve(sol, func)
145 if rhs:
146 eqs = [Eq(func, t) for t in rhs]
147 if len(rhs) == 1:
148 eqs = eqs[0]
149 return checkodesol(ode, eqs, order=order,
150 solve_for_func=False)
152 x = func.args[0]
154 # Handle series solutions here
155 if sol.has(Order):
156 assert sol.lhs == func
157 Oterm = sol.rhs.getO()
158 solrhs = sol.rhs.removeO()
160 Oexpr = Oterm.expr
161 assert isinstance(Oexpr, Pow)
162 sorder = Oexpr.exp
163 assert Oterm == Order(x**sorder)
165 odesubs = (ode.lhs-ode.rhs).subs(func, solrhs).doit().expand()
167 neworder = Order(x**(sorder - order))
168 odesubs = odesubs + neworder
169 assert odesubs.getO() == neworder
170 residual = odesubs.removeO()
172 return (residual == 0, residual)
174 s = True
175 testnum = 0
176 while s:
177 if testnum == 0:
178 # First pass, try substituting a solved solution directly into the
179 # ODE. This has the highest chance of succeeding.
180 ode_diff = ode.lhs - ode.rhs
182 if sol.lhs == func:
183 s = sub_func_doit(ode_diff, func, sol.rhs)
184 s = besselsimp(s)
185 else:
186 testnum += 1
187 continue
188 ss = simplify(s.rewrite(exp))
189 if ss:
190 # with the new numer_denom in power.py, if we do a simple
191 # expansion then testnum == 0 verifies all solutions.
192 s = ss.expand(force=True)
193 else:
194 s = 0
195 testnum += 1
196 elif testnum == 1:
197 # Second pass. If we cannot substitute f, try seeing if the nth
198 # derivative is equal, this will only work for odes that are exact,
199 # by definition.
200 s = simplify(
201 trigsimp(diff(sol.lhs, x, order) - diff(sol.rhs, x, order)) -
202 trigsimp(ode.lhs) + trigsimp(ode.rhs))
203 # s2 = simplify(
204 # diff(sol.lhs, x, order) - diff(sol.rhs, x, order) - \
205 # ode.lhs + ode.rhs)
206 testnum += 1
207 elif testnum == 2:
208 # Third pass. Try solving for df/dx and substituting that into the
209 # ODE. Thanks to Chris Smith for suggesting this method. Many of
210 # the comments below are his, too.
211 # The method:
212 # - Take each of 1..n derivatives of the solution.
213 # - Solve each nth derivative for d^(n)f/dx^(n)
214 # (the differential of that order)
215 # - Back substitute into the ODE in decreasing order
216 # (i.e., n, n-1, ...)
217 # - Check the result for zero equivalence
218 if sol.lhs == func and not sol.rhs.has(func):
219 diffsols = {0: sol.rhs}
220 elif sol.rhs == func and not sol.lhs.has(func):
221 diffsols = {0: sol.lhs}
222 else:
223 diffsols = {}
224 sol = sol.lhs - sol.rhs
225 for i in range(1, order + 1):
226 # Differentiation is a linear operator, so there should always
227 # be 1 solution. Nonetheless, we test just to make sure.
228 # We only need to solve once. After that, we automatically
229 # have the solution to the differential in the order we want.
230 if i == 1:
231 ds = sol.diff(x)
232 try:
233 sdf = solve(ds, func.diff(x, i))
234 if not sdf:
235 raise NotImplementedError
236 except NotImplementedError:
237 testnum += 1
238 break
239 else:
240 diffsols[i] = sdf[0]
241 else:
242 # This is what the solution says df/dx should be.
243 diffsols[i] = diffsols[i - 1].diff(x)
245 # Make sure the above didn't fail.
246 if testnum > 2:
247 continue
248 else:
249 # Substitute it into ODE to check for self consistency.
250 lhs, rhs = ode.lhs, ode.rhs
251 for i in range(order, -1, -1):
252 if i == 0 and 0 not in diffsols:
253 # We can only substitute f(x) if the solution was
254 # solved for f(x).
255 break
256 lhs = sub_func_doit(lhs, func.diff(x, i), diffsols[i])
257 rhs = sub_func_doit(rhs, func.diff(x, i), diffsols[i])
258 ode_or_bool = Eq(lhs, rhs)
259 ode_or_bool = simplify(ode_or_bool)
261 if isinstance(ode_or_bool, (bool, BooleanAtom)):
262 if ode_or_bool:
263 lhs = rhs = S.Zero
264 else:
265 lhs = ode_or_bool.lhs
266 rhs = ode_or_bool.rhs
267 # No sense in overworking simplify -- just prove that the
268 # numerator goes to zero
269 num = trigsimp((lhs - rhs).as_numer_denom()[0])
270 # since solutions are obtained using force=True we test
271 # using the same level of assumptions
272 ## replace function with dummy so assumptions will work
273 _func = Dummy('func')
274 num = num.subs(func, _func)
275 ## posify the expression
276 num, reps = posify(num)
277 s = simplify(num).xreplace(reps).xreplace({_func: func})
278 testnum += 1
279 else:
280 break
282 if not s:
283 return (True, s)
284 elif s is True: # The code above never was able to change s
285 raise NotImplementedError("Unable to test if " + str(sol) +
286 " is a solution to " + str(ode) + ".")
287 else:
288 return (False, s)
291def checksysodesol(eqs, sols, func=None):
292 r"""
293 Substitutes corresponding ``sols`` for each functions into each ``eqs`` and
294 checks that the result of substitutions for each equation is ``0``. The
295 equations and solutions passed can be any iterable.
297 This only works when each ``sols`` have one function only, like `x(t)` or `y(t)`.
298 For each function, ``sols`` can have a single solution or a list of solutions.
299 In most cases it will not be necessary to explicitly identify the function,
300 but if the function cannot be inferred from the original equation it
301 can be supplied through the ``func`` argument.
303 When a sequence of equations is passed, the same sequence is used to return
304 the result for each equation with each function substituted with corresponding
305 solutions.
307 It tries the following method to find zero equivalence for each equation:
309 Substitute the solutions for functions, like `x(t)` and `y(t)` into the
310 original equations containing those functions.
311 This function returns a tuple. The first item in the tuple is ``True`` if
312 the substitution results for each equation is ``0``, and ``False`` otherwise.
313 The second item in the tuple is what the substitution results in. Each element
314 of the ``list`` should always be ``0`` corresponding to each equation if the
315 first item is ``True``. Note that sometimes this function may return ``False``,
316 but with an expression that is identically equal to ``0``, instead of returning
317 ``True``. This is because :py:meth:`~sympy.simplify.simplify.simplify` cannot
318 reduce the expression to ``0``. If an expression returned by each function
319 vanishes identically, then ``sols`` really is a solution to ``eqs``.
321 If this function seems to hang, it is probably because of a difficult simplification.
323 Examples
324 ========
326 >>> from sympy import Eq, diff, symbols, sin, cos, exp, sqrt, S, Function
327 >>> from sympy.solvers.ode.subscheck import checksysodesol
328 >>> C1, C2 = symbols('C1:3')
329 >>> t = symbols('t')
330 >>> x, y = symbols('x, y', cls=Function)
331 >>> eq = (Eq(diff(x(t),t), x(t) + y(t) + 17), Eq(diff(y(t),t), -2*x(t) + y(t) + 12))
332 >>> sol = [Eq(x(t), (C1*sin(sqrt(2)*t) + C2*cos(sqrt(2)*t))*exp(t) - S(5)/3),
333 ... Eq(y(t), (sqrt(2)*C1*cos(sqrt(2)*t) - sqrt(2)*C2*sin(sqrt(2)*t))*exp(t) - S(46)/3)]
334 >>> checksysodesol(eq, sol)
335 (True, [0, 0])
336 >>> eq = (Eq(diff(x(t),t),x(t)*y(t)**4), Eq(diff(y(t),t),y(t)**3))
337 >>> sol = [Eq(x(t), C1*exp(-1/(4*(C2 + t)))), Eq(y(t), -sqrt(2)*sqrt(-1/(C2 + t))/2),
338 ... Eq(x(t), C1*exp(-1/(4*(C2 + t)))), Eq(y(t), sqrt(2)*sqrt(-1/(C2 + t))/2)]
339 >>> checksysodesol(eq, sol)
340 (True, [0, 0])
342 """
343 def _sympify(eq):
344 return list(map(sympify, eq if iterable(eq) else [eq]))
345 eqs = _sympify(eqs)
346 for i in range(len(eqs)):
347 if isinstance(eqs[i], Equality):
348 eqs[i] = eqs[i].lhs - eqs[i].rhs
349 if func is None:
350 funcs = []
351 for eq in eqs:
352 derivs = eq.atoms(Derivative)
353 func = set().union(*[d.atoms(AppliedUndef) for d in derivs])
354 funcs.extend(func)
355 funcs = list(set(funcs))
356 if not all(isinstance(func, AppliedUndef) and len(func.args) == 1 for func in funcs)\
357 and len({func.args for func in funcs})!=1:
358 raise ValueError("func must be a function of one variable, not %s" % func)
359 for sol in sols:
360 if len(sol.atoms(AppliedUndef)) != 1:
361 raise ValueError("solutions should have one function only")
362 if len(funcs) != len({sol.lhs for sol in sols}):
363 raise ValueError("number of solutions provided does not match the number of equations")
364 dictsol = {}
365 for sol in sols:
366 func = list(sol.atoms(AppliedUndef))[0]
367 if sol.rhs == func:
368 sol = sol.reversed
369 solved = sol.lhs == func and not sol.rhs.has(func)
370 if not solved:
371 rhs = solve(sol, func)
372 if not rhs:
373 raise NotImplementedError
374 else:
375 rhs = sol.rhs
376 dictsol[func] = rhs
377 checkeq = []
378 for eq in eqs:
379 for func in funcs:
380 eq = sub_func_doit(eq, func, dictsol[func])
381 ss = simplify(eq)
382 if ss != 0:
383 eq = ss.expand(force=True)
384 if eq != 0:
385 eq = sqrtdenest(eq).simplify()
386 else:
387 eq = 0
388 checkeq.append(eq)
389 if len(set(checkeq)) == 1 and list(set(checkeq))[0] == 0:
390 return (True, checkeq)
391 else:
392 return (False, checkeq)