Coverage for /usr/lib/python3/dist-packages/sympy/solvers/ode/nonhomogeneous.py: 8%
261 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-14 15:55 +0200
1r"""
2This File contains helper functions for nth_linear_constant_coeff_undetermined_coefficients,
3nth_linear_euler_eq_nonhomogeneous_undetermined_coefficients,
4nth_linear_constant_coeff_variation_of_parameters,
5and nth_linear_euler_eq_nonhomogeneous_variation_of_parameters.
7All the functions in this file are used by more than one solvers so, instead of creating
8instances in other classes for using them it is better to keep it here as separate helpers.
10"""
11from collections import defaultdict
12from sympy.core import Add, S
13from sympy.core.function import diff, expand, _mexpand, expand_mul
14from sympy.core.relational import Eq
15from sympy.core.sorting import default_sort_key
16from sympy.core.symbol import Dummy, Wild
17from sympy.functions import exp, cos, cosh, im, log, re, sin, sinh, \
18 atan2, conjugate
19from sympy.integrals import Integral
20from sympy.polys import (Poly, RootOf, rootof, roots)
21from sympy.simplify import collect, simplify, separatevars, powsimp, trigsimp # type: ignore
22from sympy.utilities import numbered_symbols
23from sympy.solvers.solvers import solve
24from sympy.matrices import wronskian
25from .subscheck import sub_func_doit
26from sympy.solvers.ode.ode import get_numbered_constants
29def _test_term(coeff, func, order):
30 r"""
31 Linear Euler ODEs have the form K*x**order*diff(y(x), x, order) = F(x),
32 where K is independent of x and y(x), order>= 0.
33 So we need to check that for each term, coeff == K*x**order from
34 some K. We have a few cases, since coeff may have several
35 different types.
36 """
37 x = func.args[0]
38 f = func.func
39 if order < 0:
40 raise ValueError("order should be greater than 0")
41 if coeff == 0:
42 return True
43 if order == 0:
44 if x in coeff.free_symbols:
45 return False
46 return True
47 if coeff.is_Mul:
48 if coeff.has(f(x)):
49 return False
50 return x**order in coeff.args
51 elif coeff.is_Pow:
52 return coeff.as_base_exp() == (x, order)
53 elif order == 1:
54 return x == coeff
55 return False
58def _get_euler_characteristic_eq_sols(eq, func, match_obj):
59 r"""
60 Returns the solution of homogeneous part of the linear euler ODE and
61 the list of roots of characteristic equation.
63 The parameter ``match_obj`` is a dict of order:coeff terms, where order is the order
64 of the derivative on each term, and coeff is the coefficient of that derivative.
66 """
67 x = func.args[0]
68 f = func.func
70 # First, set up characteristic equation.
71 chareq, symbol = S.Zero, Dummy('x')
73 for i in match_obj:
74 if i >= 0:
75 chareq += (match_obj[i]*diff(x**symbol, x, i)*x**-symbol).expand()
77 chareq = Poly(chareq, symbol)
78 chareqroots = [rootof(chareq, k) for k in range(chareq.degree())]
79 collectterms = []
81 # A generator of constants
82 constants = list(get_numbered_constants(eq, num=chareq.degree()*2))
83 constants.reverse()
85 # Create a dict root: multiplicity or charroots
86 charroots = defaultdict(int)
87 for root in chareqroots:
88 charroots[root] += 1
89 gsol = S.Zero
90 ln = log
91 for root, multiplicity in charroots.items():
92 for i in range(multiplicity):
93 if isinstance(root, RootOf):
94 gsol += (x**root) * constants.pop()
95 if multiplicity != 1:
96 raise ValueError("Value should be 1")
97 collectterms = [(0, root, 0)] + collectterms
98 elif root.is_real:
99 gsol += ln(x)**i*(x**root) * constants.pop()
100 collectterms = [(i, root, 0)] + collectterms
101 else:
102 reroot = re(root)
103 imroot = im(root)
104 gsol += ln(x)**i * (x**reroot) * (
105 constants.pop() * sin(abs(imroot)*ln(x))
106 + constants.pop() * cos(imroot*ln(x)))
107 collectterms = [(i, reroot, imroot)] + collectterms
109 gsol = Eq(f(x), gsol)
111 gensols = []
112 # Keep track of when to use sin or cos for nonzero imroot
113 for i, reroot, imroot in collectterms:
114 if imroot == 0:
115 gensols.append(ln(x)**i*x**reroot)
116 else:
117 sin_form = ln(x)**i*x**reroot*sin(abs(imroot)*ln(x))
118 if sin_form in gensols:
119 cos_form = ln(x)**i*x**reroot*cos(imroot*ln(x))
120 gensols.append(cos_form)
121 else:
122 gensols.append(sin_form)
123 return gsol, gensols
126def _solve_variation_of_parameters(eq, func, roots, homogen_sol, order, match_obj, simplify_flag=True):
127 r"""
128 Helper function for the method of variation of parameters and nonhomogeneous euler eq.
130 See the
131 :py:meth:`~sympy.solvers.ode.single.NthLinearConstantCoeffVariationOfParameters`
132 docstring for more information on this method.
134 The parameter are ``match_obj`` should be a dictionary that has the following
135 keys:
137 ``list``
138 A list of solutions to the homogeneous equation.
140 ``sol``
141 The general solution.
143 """
144 f = func.func
145 x = func.args[0]
146 r = match_obj
147 psol = 0
148 wr = wronskian(roots, x)
150 if simplify_flag:
151 wr = simplify(wr) # We need much better simplification for
152 # some ODEs. See issue 4662, for example.
153 # To reduce commonly occurring sin(x)**2 + cos(x)**2 to 1
154 wr = trigsimp(wr, deep=True, recursive=True)
155 if not wr:
156 # The wronskian will be 0 iff the solutions are not linearly
157 # independent.
158 raise NotImplementedError("Cannot find " + str(order) +
159 " solutions to the homogeneous equation necessary to apply " +
160 "variation of parameters to " + str(eq) + " (Wronskian == 0)")
161 if len(roots) != order:
162 raise NotImplementedError("Cannot find " + str(order) +
163 " solutions to the homogeneous equation necessary to apply " +
164 "variation of parameters to " +
165 str(eq) + " (number of terms != order)")
166 negoneterm = S.NegativeOne**(order)
167 for i in roots:
168 psol += negoneterm*Integral(wronskian([sol for sol in roots if sol != i], x)*r[-1]/wr, x)*i/r[order]
169 negoneterm *= -1
171 if simplify_flag:
172 psol = simplify(psol)
173 psol = trigsimp(psol, deep=True)
174 return Eq(f(x), homogen_sol.rhs + psol)
177def _get_const_characteristic_eq_sols(r, func, order):
178 r"""
179 Returns the roots of characteristic equation of constant coefficient
180 linear ODE and list of collectterms which is later on used by simplification
181 to use collect on solution.
183 The parameter `r` is a dict of order:coeff terms, where order is the order of the
184 derivative on each term, and coeff is the coefficient of that derivative.
186 """
187 x = func.args[0]
188 # First, set up characteristic equation.
189 chareq, symbol = S.Zero, Dummy('x')
191 for i in r.keys():
192 if isinstance(i, str) or i < 0:
193 pass
194 else:
195 chareq += r[i]*symbol**i
197 chareq = Poly(chareq, symbol)
198 # Can't just call roots because it doesn't return rootof for unsolveable
199 # polynomials.
200 chareqroots = roots(chareq, multiple=True)
201 if len(chareqroots) != order:
202 chareqroots = [rootof(chareq, k) for k in range(chareq.degree())]
204 chareq_is_complex = not all(i.is_real for i in chareq.all_coeffs())
206 # Create a dict root: multiplicity or charroots
207 charroots = defaultdict(int)
208 for root in chareqroots:
209 charroots[root] += 1
210 # We need to keep track of terms so we can run collect() at the end.
211 # This is necessary for constantsimp to work properly.
212 collectterms = []
213 gensols = []
214 conjugate_roots = [] # used to prevent double-use of conjugate roots
215 # Loop over roots in theorder provided by roots/rootof...
216 for root in chareqroots:
217 # but don't repoeat multiple roots.
218 if root not in charroots:
219 continue
220 multiplicity = charroots.pop(root)
221 for i in range(multiplicity):
222 if chareq_is_complex:
223 gensols.append(x**i*exp(root*x))
224 collectterms = [(i, root, 0)] + collectterms
225 continue
226 reroot = re(root)
227 imroot = im(root)
228 if imroot.has(atan2) and reroot.has(atan2):
229 # Remove this condition when re and im stop returning
230 # circular atan2 usages.
231 gensols.append(x**i*exp(root*x))
232 collectterms = [(i, root, 0)] + collectterms
233 else:
234 if root in conjugate_roots:
235 collectterms = [(i, reroot, imroot)] + collectterms
236 continue
237 if imroot == 0:
238 gensols.append(x**i*exp(reroot*x))
239 collectterms = [(i, reroot, 0)] + collectterms
240 continue
241 conjugate_roots.append(conjugate(root))
242 gensols.append(x**i*exp(reroot*x) * sin(abs(imroot) * x))
243 gensols.append(x**i*exp(reroot*x) * cos( imroot * x))
245 # This ordering is important
246 collectterms = [(i, reroot, imroot)] + collectterms
247 return gensols, collectterms
250# Ideally these kind of simplification functions shouldn't be part of solvers.
251# odesimp should be improved to handle these kind of specific simplifications.
252def _get_simplified_sol(sol, func, collectterms):
253 r"""
254 Helper function which collects the solution on
255 collectterms. Ideally this should be handled by odesimp.It is used
256 only when the simplify is set to True in dsolve.
258 The parameter ``collectterms`` is a list of tuple (i, reroot, imroot) where `i` is
259 the multiplicity of the root, reroot is real part and imroot being the imaginary part.
261 """
262 f = func.func
263 x = func.args[0]
264 collectterms.sort(key=default_sort_key)
265 collectterms.reverse()
266 assert len(sol) == 1 and sol[0].lhs == f(x)
267 sol = sol[0].rhs
268 sol = expand_mul(sol)
269 for i, reroot, imroot in collectterms:
270 sol = collect(sol, x**i*exp(reroot*x)*sin(abs(imroot)*x))
271 sol = collect(sol, x**i*exp(reroot*x)*cos(imroot*x))
272 for i, reroot, imroot in collectterms:
273 sol = collect(sol, x**i*exp(reroot*x))
274 sol = powsimp(sol)
275 return Eq(f(x), sol)
278def _undetermined_coefficients_match(expr, x, func=None, eq_homogeneous=S.Zero):
279 r"""
280 Returns a trial function match if undetermined coefficients can be applied
281 to ``expr``, and ``None`` otherwise.
283 A trial expression can be found for an expression for use with the method
284 of undetermined coefficients if the expression is an
285 additive/multiplicative combination of constants, polynomials in `x` (the
286 independent variable of expr), `\sin(a x + b)`, `\cos(a x + b)`, and
287 `e^{a x}` terms (in other words, it has a finite number of linearly
288 independent derivatives).
290 Note that you may still need to multiply each term returned here by
291 sufficient `x` to make it linearly independent with the solutions to the
292 homogeneous equation.
294 This is intended for internal use by ``undetermined_coefficients`` hints.
296 SymPy currently has no way to convert `\sin^n(x) \cos^m(y)` into a sum of
297 only `\sin(a x)` and `\cos(b x)` terms, so these are not implemented. So,
298 for example, you will need to manually convert `\sin^2(x)` into `[1 +
299 \cos(2 x)]/2` to properly apply the method of undetermined coefficients on
300 it.
302 Examples
303 ========
305 >>> from sympy import log, exp
306 >>> from sympy.solvers.ode.nonhomogeneous import _undetermined_coefficients_match
307 >>> from sympy.abc import x
308 >>> _undetermined_coefficients_match(9*x*exp(x) + exp(-x), x)
309 {'test': True, 'trialset': {x*exp(x), exp(-x), exp(x)}}
310 >>> _undetermined_coefficients_match(log(x), x)
311 {'test': False}
313 """
314 a = Wild('a', exclude=[x])
315 b = Wild('b', exclude=[x])
316 expr = powsimp(expr, combine='exp') # exp(x)*exp(2*x + 1) => exp(3*x + 1)
317 retdict = {}
319 def _test_term(expr, x):
320 r"""
321 Test if ``expr`` fits the proper form for undetermined coefficients.
322 """
323 if not expr.has(x):
324 return True
325 elif expr.is_Add:
326 return all(_test_term(i, x) for i in expr.args)
327 elif expr.is_Mul:
328 if expr.has(sin, cos):
329 foundtrig = False
330 # Make sure that there is only one trig function in the args.
331 # See the docstring.
332 for i in expr.args:
333 if i.has(sin, cos):
334 if foundtrig:
335 return False
336 else:
337 foundtrig = True
338 return all(_test_term(i, x) for i in expr.args)
339 elif expr.is_Function:
340 if expr.func in (sin, cos, exp, sinh, cosh):
341 if expr.args[0].match(a*x + b):
342 return True
343 else:
344 return False
345 else:
346 return False
347 elif expr.is_Pow and expr.base.is_Symbol and expr.exp.is_Integer and \
348 expr.exp >= 0:
349 return True
350 elif expr.is_Pow and expr.base.is_number:
351 if expr.exp.match(a*x + b):
352 return True
353 else:
354 return False
355 elif expr.is_Symbol or expr.is_number:
356 return True
357 else:
358 return False
360 def _get_trial_set(expr, x, exprs=set()):
361 r"""
362 Returns a set of trial terms for undetermined coefficients.
364 The idea behind undetermined coefficients is that the terms expression
365 repeat themselves after a finite number of derivatives, except for the
366 coefficients (they are linearly dependent). So if we collect these,
367 we should have the terms of our trial function.
368 """
369 def _remove_coefficient(expr, x):
370 r"""
371 Returns the expression without a coefficient.
373 Similar to expr.as_independent(x)[1], except it only works
374 multiplicatively.
375 """
376 term = S.One
377 if expr.is_Mul:
378 for i in expr.args:
379 if i.has(x):
380 term *= i
381 elif expr.has(x):
382 term = expr
383 return term
385 expr = expand_mul(expr)
386 if expr.is_Add:
387 for term in expr.args:
388 if _remove_coefficient(term, x) in exprs:
389 pass
390 else:
391 exprs.add(_remove_coefficient(term, x))
392 exprs = exprs.union(_get_trial_set(term, x, exprs))
393 else:
394 term = _remove_coefficient(expr, x)
395 tmpset = exprs.union({term})
396 oldset = set()
397 while tmpset != oldset:
398 # If you get stuck in this loop, then _test_term is probably
399 # broken
400 oldset = tmpset.copy()
401 expr = expr.diff(x)
402 term = _remove_coefficient(expr, x)
403 if term.is_Add:
404 tmpset = tmpset.union(_get_trial_set(term, x, tmpset))
405 else:
406 tmpset.add(term)
407 exprs = tmpset
408 return exprs
410 def is_homogeneous_solution(term):
411 r""" This function checks whether the given trialset contains any root
412 of homogeneous equation"""
413 return expand(sub_func_doit(eq_homogeneous, func, term)).is_zero
415 retdict['test'] = _test_term(expr, x)
416 if retdict['test']:
417 # Try to generate a list of trial solutions that will have the
418 # undetermined coefficients. Note that if any of these are not linearly
419 # independent with any of the solutions to the homogeneous equation,
420 # then they will need to be multiplied by sufficient x to make them so.
421 # This function DOES NOT do that (it doesn't even look at the
422 # homogeneous equation).
423 temp_set = set()
424 for i in Add.make_args(expr):
425 act = _get_trial_set(i, x)
426 if eq_homogeneous is not S.Zero:
427 while any(is_homogeneous_solution(ts) for ts in act):
428 act = {x*ts for ts in act}
429 temp_set = temp_set.union(act)
431 retdict['trialset'] = temp_set
432 return retdict
435def _solve_undetermined_coefficients(eq, func, order, match, trialset):
436 r"""
437 Helper function for the method of undetermined coefficients.
439 See the
440 :py:meth:`~sympy.solvers.ode.single.NthLinearConstantCoeffUndeterminedCoefficients`
441 docstring for more information on this method.
443 The parameter ``trialset`` is the set of trial functions as returned by
444 ``_undetermined_coefficients_match()['trialset']``.
446 The parameter ``match`` should be a dictionary that has the following
447 keys:
449 ``list``
450 A list of solutions to the homogeneous equation.
452 ``sol``
453 The general solution.
455 """
456 r = match
457 coeffs = numbered_symbols('a', cls=Dummy)
458 coefflist = []
459 gensols = r['list']
460 gsol = r['sol']
461 f = func.func
462 x = func.args[0]
464 if len(gensols) != order:
465 raise NotImplementedError("Cannot find " + str(order) +
466 " solutions to the homogeneous equation necessary to apply" +
467 " undetermined coefficients to " + str(eq) +
468 " (number of terms != order)")
470 trialfunc = 0
471 for i in trialset:
472 c = next(coeffs)
473 coefflist.append(c)
474 trialfunc += c*i
476 eqs = sub_func_doit(eq, f(x), trialfunc)
478 coeffsdict = dict(list(zip(trialset, [0]*(len(trialset) + 1))))
480 eqs = _mexpand(eqs)
482 for i in Add.make_args(eqs):
483 s = separatevars(i, dict=True, symbols=[x])
484 if coeffsdict.get(s[x]):
485 coeffsdict[s[x]] += s['coeff']
486 else:
487 coeffsdict[s[x]] = s['coeff']
489 coeffvals = solve(list(coeffsdict.values()), coefflist)
491 if not coeffvals:
492 raise NotImplementedError(
493 "Could not solve `%s` using the "
494 "method of undetermined coefficients "
495 "(unable to solve for coefficients)." % eq)
497 psol = trialfunc.subs(coeffvals)
499 return Eq(f(x), gsol.rhs + psol)