Coverage for /usr/lib/python3/dist-packages/sympy/solvers/ode/systems.py: 9%
703 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-14 15:55 +0200
1from sympy.core import Add, Mul, S
2from sympy.core.containers import Tuple
3from sympy.core.exprtools import factor_terms
4from sympy.core.numbers import I
5from sympy.core.relational import Eq, Equality
6from sympy.core.sorting import default_sort_key, ordered
7from sympy.core.symbol import Dummy, Symbol
8from sympy.core.function import (expand_mul, expand, Derivative,
9 AppliedUndef, Function, Subs)
10from sympy.functions import (exp, im, cos, sin, re, Piecewise,
11 piecewise_fold, sqrt, log)
12from sympy.functions.combinatorial.factorials import factorial
13from sympy.matrices import zeros, Matrix, NonSquareMatrixError, MatrixBase, eye
14from sympy.polys import Poly, together
15from sympy.simplify import collect, radsimp, signsimp # type: ignore
16from sympy.simplify.powsimp import powdenest, powsimp
17from sympy.simplify.ratsimp import ratsimp
18from sympy.simplify.simplify import simplify
19from sympy.sets.sets import FiniteSet
20from sympy.solvers.deutils import ode_order
21from sympy.solvers.solveset import NonlinearError, solveset
22from sympy.utilities.iterables import (connected_components, iterable,
23 strongly_connected_components)
24from sympy.utilities.misc import filldedent
25from sympy.integrals.integrals import Integral, integrate
28def _get_func_order(eqs, funcs):
29 return {func: max(ode_order(eq, func) for eq in eqs) for func in funcs}
32class ODEOrderError(ValueError):
33 """Raised by linear_ode_to_matrix if the system has the wrong order"""
34 pass
37class ODENonlinearError(NonlinearError):
38 """Raised by linear_ode_to_matrix if the system is nonlinear"""
39 pass
42def _simpsol(soleq):
43 lhs = soleq.lhs
44 sol = soleq.rhs
45 sol = powsimp(sol)
46 gens = list(sol.atoms(exp))
47 p = Poly(sol, *gens, expand=False)
48 gens = [factor_terms(g) for g in gens]
49 if not gens:
50 gens = p.gens
51 syms = [Symbol('C1'), Symbol('C2')]
52 terms = []
53 for coeff, monom in zip(p.coeffs(), p.monoms()):
54 coeff = piecewise_fold(coeff)
55 if isinstance(coeff, Piecewise):
56 coeff = Piecewise(*((ratsimp(coef).collect(syms), cond) for coef, cond in coeff.args))
57 else:
58 coeff = ratsimp(coeff).collect(syms)
59 monom = Mul(*(g ** i for g, i in zip(gens, monom)))
60 terms.append(coeff * monom)
61 return Eq(lhs, Add(*terms))
64def _solsimp(e, t):
65 no_t, has_t = powsimp(expand_mul(e)).as_independent(t)
67 no_t = ratsimp(no_t)
68 has_t = has_t.replace(exp, lambda a: exp(factor_terms(a)))
70 return no_t + has_t
73def simpsol(sol, wrt1, wrt2, doit=True):
74 """Simplify solutions from dsolve_system."""
76 # The parameter sol is the solution as returned by dsolve (list of Eq).
77 #
78 # The parameters wrt1 and wrt2 are lists of symbols to be collected for
79 # with those in wrt1 being collected for first. This allows for collecting
80 # on any factors involving the independent variable before collecting on
81 # the integration constants or vice versa using e.g.:
82 #
83 # sol = simpsol(sol, [t], [C1, C2]) # t first, constants after
84 # sol = simpsol(sol, [C1, C2], [t]) # constants first, t after
85 #
86 # If doit=True (default) then simpsol will begin by evaluating any
87 # unevaluated integrals. Since many integrals will appear multiple times
88 # in the solutions this is done intelligently by computing each integral
89 # only once.
90 #
91 # The strategy is to first perform simple cancellation with factor_terms
92 # and then multiply out all brackets with expand_mul. This gives an Add
93 # with many terms.
94 #
95 # We split each term into two multiplicative factors dep and coeff where
96 # all factors that involve wrt1 are in dep and any constant factors are in
97 # coeff e.g.
98 # sqrt(2)*C1*exp(t) -> ( exp(t), sqrt(2)*C1 )
99 #
100 # The dep factors are simplified using powsimp to combine expanded
101 # exponential factors e.g.
102 # exp(a*t)*exp(b*t) -> exp(t*(a+b))
103 #
104 # We then collect coefficients for all terms having the same (simplified)
105 # dep. The coefficients are then simplified using together and ratsimp and
106 # lastly by recursively applying the same transformation to the
107 # coefficients to collect on wrt2.
108 #
109 # Finally the result is recombined into an Add and signsimp is used to
110 # normalise any minus signs.
112 def simprhs(rhs, rep, wrt1, wrt2):
113 """Simplify the rhs of an ODE solution"""
114 if rep:
115 rhs = rhs.subs(rep)
116 rhs = factor_terms(rhs)
117 rhs = simp_coeff_dep(rhs, wrt1, wrt2)
118 rhs = signsimp(rhs)
119 return rhs
121 def simp_coeff_dep(expr, wrt1, wrt2=None):
122 """Split rhs into terms, split terms into dep and coeff and collect on dep"""
123 add_dep_terms = lambda e: e.is_Add and e.has(*wrt1)
124 expandable = lambda e: e.is_Mul and any(map(add_dep_terms, e.args))
125 expand_func = lambda e: expand_mul(e, deep=False)
126 expand_mul_mod = lambda e: e.replace(expandable, expand_func)
127 terms = Add.make_args(expand_mul_mod(expr))
128 dc = {}
129 for term in terms:
130 coeff, dep = term.as_independent(*wrt1, as_Add=False)
131 # Collect together the coefficients for terms that have the same
132 # dependence on wrt1 (after dep is normalised using simpdep).
133 dep = simpdep(dep, wrt1)
135 # See if the dependence on t cancels out...
136 if dep is not S.One:
137 dep2 = factor_terms(dep)
138 if not dep2.has(*wrt1):
139 coeff *= dep2
140 dep = S.One
142 if dep not in dc:
143 dc[dep] = coeff
144 else:
145 dc[dep] += coeff
146 # Apply the method recursively to the coefficients but this time
147 # collecting on wrt2 rather than wrt2.
148 termpairs = ((simpcoeff(c, wrt2), d) for d, c in dc.items())
149 if wrt2 is not None:
150 termpairs = ((simp_coeff_dep(c, wrt2), d) for c, d in termpairs)
151 return Add(*(c * d for c, d in termpairs))
153 def simpdep(term, wrt1):
154 """Normalise factors involving t with powsimp and recombine exp"""
155 def canonicalise(a):
156 # Using factor_terms here isn't quite right because it leads to things
157 # like exp(t*(1+t)) that we don't want. We do want to cancel factors
158 # and pull out a common denominator but ideally the numerator would be
159 # expressed as a standard form polynomial in t so we expand_mul
160 # and collect afterwards.
161 a = factor_terms(a)
162 num, den = a.as_numer_denom()
163 num = expand_mul(num)
164 num = collect(num, wrt1)
165 return num / den
167 term = powsimp(term)
168 rep = {e: exp(canonicalise(e.args[0])) for e in term.atoms(exp)}
169 term = term.subs(rep)
170 return term
172 def simpcoeff(coeff, wrt2):
173 """Bring to a common fraction and cancel with ratsimp"""
174 coeff = together(coeff)
175 if coeff.is_polynomial():
176 # Calling ratsimp can be expensive. The main reason is to simplify
177 # sums of terms with irrational denominators so we limit ourselves
178 # to the case where the expression is polynomial in any symbols.
179 # Maybe there's a better approach...
180 coeff = ratsimp(radsimp(coeff))
181 # collect on secondary variables first and any remaining symbols after
182 if wrt2 is not None:
183 syms = list(wrt2) + list(ordered(coeff.free_symbols - set(wrt2)))
184 else:
185 syms = list(ordered(coeff.free_symbols))
186 coeff = collect(coeff, syms)
187 coeff = together(coeff)
188 return coeff
190 # There are often repeated integrals. Collect unique integrals and
191 # evaluate each once and then substitute into the final result to replace
192 # all occurrences in each of the solution equations.
193 if doit:
194 integrals = set().union(*(s.atoms(Integral) for s in sol))
195 rep = {i: factor_terms(i).doit() for i in integrals}
196 else:
197 rep = {}
199 sol = [Eq(s.lhs, simprhs(s.rhs, rep, wrt1, wrt2)) for s in sol]
200 return sol
203def linodesolve_type(A, t, b=None):
204 r"""
205 Helper function that determines the type of the system of ODEs for solving with :obj:`sympy.solvers.ode.systems.linodesolve()`
207 Explanation
208 ===========
210 This function takes in the coefficient matrix and/or the non-homogeneous term
211 and returns the type of the equation that can be solved by :obj:`sympy.solvers.ode.systems.linodesolve()`.
213 If the system is constant coefficient homogeneous, then "type1" is returned
215 If the system is constant coefficient non-homogeneous, then "type2" is returned
217 If the system is non-constant coefficient homogeneous, then "type3" is returned
219 If the system is non-constant coefficient non-homogeneous, then "type4" is returned
221 If the system has a non-constant coefficient matrix which can be factorized into constant
222 coefficient matrix, then "type5" or "type6" is returned for when the system is homogeneous or
223 non-homogeneous respectively.
225 Note that, if the system of ODEs is of "type3" or "type4", then along with the type,
226 the commutative antiderivative of the coefficient matrix is also returned.
228 If the system cannot be solved by :obj:`sympy.solvers.ode.systems.linodesolve()`, then
229 NotImplementedError is raised.
231 Parameters
232 ==========
234 A : Matrix
235 Coefficient matrix of the system of ODEs
236 b : Matrix or None
237 Non-homogeneous term of the system. The default value is None.
238 If this argument is None, then the system is assumed to be homogeneous.
240 Examples
241 ========
243 >>> from sympy import symbols, Matrix
244 >>> from sympy.solvers.ode.systems import linodesolve_type
245 >>> t = symbols("t")
246 >>> A = Matrix([[1, 1], [2, 3]])
247 >>> b = Matrix([t, 1])
249 >>> linodesolve_type(A, t)
250 {'antiderivative': None, 'type_of_equation': 'type1'}
252 >>> linodesolve_type(A, t, b=b)
253 {'antiderivative': None, 'type_of_equation': 'type2'}
255 >>> A_t = Matrix([[1, t], [-t, 1]])
257 >>> linodesolve_type(A_t, t)
258 {'antiderivative': Matrix([
259 [ t, t**2/2],
260 [-t**2/2, t]]), 'type_of_equation': 'type3'}
262 >>> linodesolve_type(A_t, t, b=b)
263 {'antiderivative': Matrix([
264 [ t, t**2/2],
265 [-t**2/2, t]]), 'type_of_equation': 'type4'}
267 >>> A_non_commutative = Matrix([[1, t], [t, -1]])
268 >>> linodesolve_type(A_non_commutative, t)
269 Traceback (most recent call last):
270 ...
271 NotImplementedError:
272 The system does not have a commutative antiderivative, it cannot be
273 solved by linodesolve.
275 Returns
276 =======
278 Dict
280 Raises
281 ======
283 NotImplementedError
284 When the coefficient matrix does not have a commutative antiderivative
286 See Also
287 ========
289 linodesolve: Function for which linodesolve_type gets the information
291 """
293 match = {}
294 is_non_constant = not _matrix_is_constant(A, t)
295 is_non_homogeneous = not (b is None or b.is_zero_matrix)
296 type = "type{}".format(int("{}{}".format(int(is_non_constant), int(is_non_homogeneous)), 2) + 1)
298 B = None
299 match.update({"type_of_equation": type, "antiderivative": B})
301 if is_non_constant:
302 B, is_commuting = _is_commutative_anti_derivative(A, t)
303 if not is_commuting:
304 raise NotImplementedError(filldedent('''
305 The system does not have a commutative antiderivative, it cannot be solved
306 by linodesolve.
307 '''))
309 match['antiderivative'] = B
310 match.update(_first_order_type5_6_subs(A, t, b=b))
312 return match
315def _first_order_type5_6_subs(A, t, b=None):
316 match = {}
318 factor_terms = _factor_matrix(A, t)
319 is_homogeneous = b is None or b.is_zero_matrix
321 if factor_terms is not None:
322 t_ = Symbol("{}_".format(t))
323 F_t = integrate(factor_terms[0], t)
324 inverse = solveset(Eq(t_, F_t), t)
326 # Note: A simple way to check if a function is invertible
327 # or not.
328 if isinstance(inverse, FiniteSet) and not inverse.has(Piecewise)\
329 and len(inverse) == 1:
331 A = factor_terms[1]
332 if not is_homogeneous:
333 b = b / factor_terms[0]
334 b = b.subs(t, list(inverse)[0])
335 type = "type{}".format(5 + (not is_homogeneous))
336 match.update({'func_coeff': A, 'tau': F_t,
337 't_': t_, 'type_of_equation': type, 'rhs': b})
339 return match
342def linear_ode_to_matrix(eqs, funcs, t, order):
343 r"""
344 Convert a linear system of ODEs to matrix form
346 Explanation
347 ===========
349 Express a system of linear ordinary differential equations as a single
350 matrix differential equation [1]. For example the system $x' = x + y + 1$
351 and $y' = x - y$ can be represented as
353 .. math:: A_1 X' = A_0 X + b
355 where $A_1$ and $A_0$ are $2 \times 2$ matrices and $b$, $X$ and $X'$ are
356 $2 \times 1$ matrices with $X = [x, y]^T$.
358 Higher-order systems are represented with additional matrices e.g. a
359 second-order system would look like
361 .. math:: A_2 X'' = A_1 X' + A_0 X + b
363 Examples
364 ========
366 >>> from sympy import Function, Symbol, Matrix, Eq
367 >>> from sympy.solvers.ode.systems import linear_ode_to_matrix
368 >>> t = Symbol('t')
369 >>> x = Function('x')
370 >>> y = Function('y')
372 We can create a system of linear ODEs like
374 >>> eqs = [
375 ... Eq(x(t).diff(t), x(t) + y(t) + 1),
376 ... Eq(y(t).diff(t), x(t) - y(t)),
377 ... ]
378 >>> funcs = [x(t), y(t)]
379 >>> order = 1 # 1st order system
381 Now ``linear_ode_to_matrix`` can represent this as a matrix
382 differential equation.
384 >>> (A1, A0), b = linear_ode_to_matrix(eqs, funcs, t, order)
385 >>> A1
386 Matrix([
387 [1, 0],
388 [0, 1]])
389 >>> A0
390 Matrix([
391 [1, 1],
392 [1, -1]])
393 >>> b
394 Matrix([
395 [1],
396 [0]])
398 The original equations can be recovered from these matrices:
400 >>> eqs_mat = Matrix([eq.lhs - eq.rhs for eq in eqs])
401 >>> X = Matrix(funcs)
402 >>> A1 * X.diff(t) - A0 * X - b == eqs_mat
403 True
405 If the system of equations has a maximum order greater than the
406 order of the system specified, a ODEOrderError exception is raised.
408 >>> eqs = [Eq(x(t).diff(t, 2), x(t).diff(t) + x(t)), Eq(y(t).diff(t), y(t) + x(t))]
409 >>> linear_ode_to_matrix(eqs, funcs, t, 1)
410 Traceback (most recent call last):
411 ...
412 ODEOrderError: Cannot represent system in 1-order form
414 If the system of equations is nonlinear, then ODENonlinearError is
415 raised.
417 >>> eqs = [Eq(x(t).diff(t), x(t) + y(t)), Eq(y(t).diff(t), y(t)**2 + x(t))]
418 >>> linear_ode_to_matrix(eqs, funcs, t, 1)
419 Traceback (most recent call last):
420 ...
421 ODENonlinearError: The system of ODEs is nonlinear.
423 Parameters
424 ==========
426 eqs : list of SymPy expressions or equalities
427 The equations as expressions (assumed equal to zero).
428 funcs : list of applied functions
429 The dependent variables of the system of ODEs.
430 t : symbol
431 The independent variable.
432 order : int
433 The order of the system of ODEs.
435 Returns
436 =======
438 The tuple ``(As, b)`` where ``As`` is a tuple of matrices and ``b`` is the
439 the matrix representing the rhs of the matrix equation.
441 Raises
442 ======
444 ODEOrderError
445 When the system of ODEs have an order greater than what was specified
446 ODENonlinearError
447 When the system of ODEs is nonlinear
449 See Also
450 ========
452 linear_eq_to_matrix: for systems of linear algebraic equations.
454 References
455 ==========
457 .. [1] https://en.wikipedia.org/wiki/Matrix_differential_equation
459 """
460 from sympy.solvers.solveset import linear_eq_to_matrix
462 if any(ode_order(eq, func) > order for eq in eqs for func in funcs):
463 msg = "Cannot represent system in {}-order form"
464 raise ODEOrderError(msg.format(order))
466 As = []
468 for o in range(order, -1, -1):
469 # Work from the highest derivative down
470 syms = [func.diff(t, o) for func in funcs]
472 # Ai is the matrix for X(t).diff(t, o)
473 # eqs is minus the remainder of the equations.
474 try:
475 Ai, b = linear_eq_to_matrix(eqs, syms)
476 except NonlinearError:
477 raise ODENonlinearError("The system of ODEs is nonlinear.")
479 Ai = Ai.applyfunc(expand_mul)
481 As.append(Ai if o == order else -Ai)
483 if o:
484 eqs = [-eq for eq in b]
485 else:
486 rhs = b
488 return As, rhs
491def matrix_exp(A, t):
492 r"""
493 Matrix exponential $\exp(A*t)$ for the matrix ``A`` and scalar ``t``.
495 Explanation
496 ===========
498 This functions returns the $\exp(A*t)$ by doing a simple
499 matrix multiplication:
501 .. math:: \exp(A*t) = P * expJ * P^{-1}
503 where $expJ$ is $\exp(J*t)$. $J$ is the Jordan normal
504 form of $A$ and $P$ is matrix such that:
506 .. math:: A = P * J * P^{-1}
508 The matrix exponential $\exp(A*t)$ appears in the solution of linear
509 differential equations. For example if $x$ is a vector and $A$ is a matrix
510 then the initial value problem
512 .. math:: \frac{dx(t)}{dt} = A \times x(t), x(0) = x0
514 has the unique solution
516 .. math:: x(t) = \exp(A t) x0
518 Examples
519 ========
521 >>> from sympy import Symbol, Matrix, pprint
522 >>> from sympy.solvers.ode.systems import matrix_exp
523 >>> t = Symbol('t')
525 We will consider a 2x2 matrix for comupting the exponential
527 >>> A = Matrix([[2, -5], [2, -4]])
528 >>> pprint(A)
529 [2 -5]
530 [ ]
531 [2 -4]
533 Now, exp(A*t) is given as follows:
535 >>> pprint(matrix_exp(A, t))
536 [ -t -t -t ]
537 [3*e *sin(t) + e *cos(t) -5*e *sin(t) ]
538 [ ]
539 [ -t -t -t ]
540 [ 2*e *sin(t) - 3*e *sin(t) + e *cos(t)]
542 Parameters
543 ==========
545 A : Matrix
546 The matrix $A$ in the expression $\exp(A*t)$
547 t : Symbol
548 The independent variable
550 See Also
551 ========
553 matrix_exp_jordan_form: For exponential of Jordan normal form
555 References
556 ==========
558 .. [1] https://en.wikipedia.org/wiki/Jordan_normal_form
559 .. [2] https://en.wikipedia.org/wiki/Matrix_exponential
561 """
562 P, expJ = matrix_exp_jordan_form(A, t)
563 return P * expJ * P.inv()
566def matrix_exp_jordan_form(A, t):
567 r"""
568 Matrix exponential $\exp(A*t)$ for the matrix *A* and scalar *t*.
570 Explanation
571 ===========
573 Returns the Jordan form of the $\exp(A*t)$ along with the matrix $P$ such that:
575 .. math::
576 \exp(A*t) = P * expJ * P^{-1}
578 Examples
579 ========
581 >>> from sympy import Matrix, Symbol
582 >>> from sympy.solvers.ode.systems import matrix_exp, matrix_exp_jordan_form
583 >>> t = Symbol('t')
585 We will consider a 2x2 defective matrix. This shows that our method
586 works even for defective matrices.
588 >>> A = Matrix([[1, 1], [0, 1]])
590 It can be observed that this function gives us the Jordan normal form
591 and the required invertible matrix P.
593 >>> P, expJ = matrix_exp_jordan_form(A, t)
595 Here, it is shown that P and expJ returned by this function is correct
596 as they satisfy the formula: P * expJ * P_inverse = exp(A*t).
598 >>> P * expJ * P.inv() == matrix_exp(A, t)
599 True
601 Parameters
602 ==========
604 A : Matrix
605 The matrix $A$ in the expression $\exp(A*t)$
606 t : Symbol
607 The independent variable
609 References
610 ==========
612 .. [1] https://en.wikipedia.org/wiki/Defective_matrix
613 .. [2] https://en.wikipedia.org/wiki/Jordan_matrix
614 .. [3] https://en.wikipedia.org/wiki/Jordan_normal_form
616 """
618 N, M = A.shape
619 if N != M:
620 raise ValueError('Needed square matrix but got shape (%s, %s)' % (N, M))
621 elif A.has(t):
622 raise ValueError('Matrix A should not depend on t')
624 def jordan_chains(A):
625 '''Chains from Jordan normal form analogous to M.eigenvects().
626 Returns a dict with eignevalues as keys like:
627 {e1: [[v111,v112,...], [v121, v122,...]], e2:...}
628 where vijk is the kth vector in the jth chain for eigenvalue i.
629 '''
630 P, blocks = A.jordan_cells()
631 basis = [P[:,i] for i in range(P.shape[1])]
632 n = 0
633 chains = {}
634 for b in blocks:
635 eigval = b[0, 0]
636 size = b.shape[0]
637 if eigval not in chains:
638 chains[eigval] = []
639 chains[eigval].append(basis[n:n+size])
640 n += size
641 return chains
643 eigenchains = jordan_chains(A)
645 # Needed for consistency across Python versions
646 eigenchains_iter = sorted(eigenchains.items(), key=default_sort_key)
647 isreal = not A.has(I)
649 blocks = []
650 vectors = []
651 seen_conjugate = set()
652 for e, chains in eigenchains_iter:
653 for chain in chains:
654 n = len(chain)
655 if isreal and e != e.conjugate() and e.conjugate() in eigenchains:
656 if e in seen_conjugate:
657 continue
658 seen_conjugate.add(e.conjugate())
659 exprt = exp(re(e) * t)
660 imrt = im(e) * t
661 imblock = Matrix([[cos(imrt), sin(imrt)],
662 [-sin(imrt), cos(imrt)]])
663 expJblock2 = Matrix(n, n, lambda i,j:
664 imblock * t**(j-i) / factorial(j-i) if j >= i
665 else zeros(2, 2))
666 expJblock = Matrix(2*n, 2*n, lambda i,j: expJblock2[i//2,j//2][i%2,j%2])
668 blocks.append(exprt * expJblock)
669 for i in range(n):
670 vectors.append(re(chain[i]))
671 vectors.append(im(chain[i]))
672 else:
673 vectors.extend(chain)
674 fun = lambda i,j: t**(j-i)/factorial(j-i) if j >= i else 0
675 expJblock = Matrix(n, n, fun)
676 blocks.append(exp(e * t) * expJblock)
678 expJ = Matrix.diag(*blocks)
679 P = Matrix(N, N, lambda i,j: vectors[j][i])
681 return P, expJ
684# Note: To add a docstring example with tau
685def linodesolve(A, t, b=None, B=None, type="auto", doit=False,
686 tau=None):
687 r"""
688 System of n equations linear first-order differential equations
690 Explanation
691 ===========
693 This solver solves the system of ODEs of the following form:
695 .. math::
696 X'(t) = A(t) X(t) + b(t)
698 Here, $A(t)$ is the coefficient matrix, $X(t)$ is the vector of n independent variables,
699 $b(t)$ is the non-homogeneous term and $X'(t)$ is the derivative of $X(t)$
701 Depending on the properties of $A(t)$ and $b(t)$, this solver evaluates the solution
702 differently.
704 When $A(t)$ is constant coefficient matrix and $b(t)$ is zero vector i.e. system is homogeneous,
705 the system is "type1". The solution is:
707 .. math::
708 X(t) = \exp(A t) C
710 Here, $C$ is a vector of constants and $A$ is the constant coefficient matrix.
712 When $A(t)$ is constant coefficient matrix and $b(t)$ is non-zero i.e. system is non-homogeneous,
713 the system is "type2". The solution is:
715 .. math::
716 X(t) = e^{A t} ( \int e^{- A t} b \,dt + C)
718 When $A(t)$ is coefficient matrix such that its commutative with its antiderivative $B(t)$ and
719 $b(t)$ is a zero vector i.e. system is homogeneous, the system is "type3". The solution is:
721 .. math::
722 X(t) = \exp(B(t)) C
724 When $A(t)$ is commutative with its antiderivative $B(t)$ and $b(t)$ is non-zero i.e. system is
725 non-homogeneous, the system is "type4". The solution is:
727 .. math::
728 X(t) = e^{B(t)} ( \int e^{-B(t)} b(t) \,dt + C)
730 When $A(t)$ is a coefficient matrix such that it can be factorized into a scalar and a constant
731 coefficient matrix:
733 .. math::
734 A(t) = f(t) * A
736 Where $f(t)$ is a scalar expression in the independent variable $t$ and $A$ is a constant matrix,
737 then we can do the following substitutions:
739 .. math::
740 tau = \int f(t) dt, X(t) = Y(tau), b(t) = b(f^{-1}(tau))
742 Here, the substitution for the non-homogeneous term is done only when its non-zero.
743 Using these substitutions, our original system becomes:
745 .. math::
746 Y'(tau) = A * Y(tau) + b(tau)/f(tau)
748 The above system can be easily solved using the solution for "type1" or "type2" depending
749 on the homogeneity of the system. After we get the solution for $Y(tau)$, we substitute the
750 solution for $tau$ as $t$ to get back $X(t)$
752 .. math::
753 X(t) = Y(tau)
755 Systems of "type5" and "type6" have a commutative antiderivative but we use this solution
756 because its faster to compute.
758 The final solution is the general solution for all the four equations since a constant coefficient
759 matrix is always commutative with its antidervative.
761 An additional feature of this function is, if someone wants to substitute for value of the independent
762 variable, they can pass the substitution `tau` and the solution will have the independent variable
763 substituted with the passed expression(`tau`).
765 Parameters
766 ==========
768 A : Matrix
769 Coefficient matrix of the system of linear first order ODEs.
770 t : Symbol
771 Independent variable in the system of ODEs.
772 b : Matrix or None
773 Non-homogeneous term in the system of ODEs. If None is passed,
774 a homogeneous system of ODEs is assumed.
775 B : Matrix or None
776 Antiderivative of the coefficient matrix. If the antiderivative
777 is not passed and the solution requires the term, then the solver
778 would compute it internally.
779 type : String
780 Type of the system of ODEs passed. Depending on the type, the
781 solution is evaluated. The type values allowed and the corresponding
782 system it solves are: "type1" for constant coefficient homogeneous
783 "type2" for constant coefficient non-homogeneous, "type3" for non-constant
784 coefficient homogeneous, "type4" for non-constant coefficient non-homogeneous,
785 "type5" and "type6" for non-constant coefficient homogeneous and non-homogeneous
786 systems respectively where the coefficient matrix can be factorized to a constant
787 coefficient matrix.
788 The default value is "auto" which will let the solver decide the correct type of
789 the system passed.
790 doit : Boolean
791 Evaluate the solution if True, default value is False
792 tau: Expression
793 Used to substitute for the value of `t` after we get the solution of the system.
795 Examples
796 ========
798 To solve the system of ODEs using this function directly, several things must be
799 done in the right order. Wrong inputs to the function will lead to incorrect results.
801 >>> from sympy import symbols, Function, Eq
802 >>> from sympy.solvers.ode.systems import canonical_odes, linear_ode_to_matrix, linodesolve, linodesolve_type
803 >>> from sympy.solvers.ode.subscheck import checkodesol
804 >>> f, g = symbols("f, g", cls=Function)
805 >>> x, a = symbols("x, a")
806 >>> funcs = [f(x), g(x)]
807 >>> eqs = [Eq(f(x).diff(x) - f(x), a*g(x) + 1), Eq(g(x).diff(x) + g(x), a*f(x))]
809 Here, it is important to note that before we derive the coefficient matrix, it is
810 important to get the system of ODEs into the desired form. For that we will use
811 :obj:`sympy.solvers.ode.systems.canonical_odes()`.
813 >>> eqs = canonical_odes(eqs, funcs, x)
814 >>> eqs
815 [[Eq(Derivative(f(x), x), a*g(x) + f(x) + 1), Eq(Derivative(g(x), x), a*f(x) - g(x))]]
817 Now, we will use :obj:`sympy.solvers.ode.systems.linear_ode_to_matrix()` to get the coefficient matrix and the
818 non-homogeneous term if it is there.
820 >>> eqs = eqs[0]
821 >>> (A1, A0), b = linear_ode_to_matrix(eqs, funcs, x, 1)
822 >>> A = A0
824 We have the coefficient matrices and the non-homogeneous term ready. Now, we can use
825 :obj:`sympy.solvers.ode.systems.linodesolve_type()` to get the information for the system of ODEs
826 to finally pass it to the solver.
828 >>> system_info = linodesolve_type(A, x, b=b)
829 >>> sol_vector = linodesolve(A, x, b=b, B=system_info['antiderivative'], type=system_info['type_of_equation'])
831 Now, we can prove if the solution is correct or not by using :obj:`sympy.solvers.ode.checkodesol()`
833 >>> sol = [Eq(f, s) for f, s in zip(funcs, sol_vector)]
834 >>> checkodesol(eqs, sol)
835 (True, [0, 0])
837 We can also use the doit method to evaluate the solutions passed by the function.
839 >>> sol_vector_evaluated = linodesolve(A, x, b=b, type="type2", doit=True)
841 Now, we will look at a system of ODEs which is non-constant.
843 >>> eqs = [Eq(f(x).diff(x), f(x) + x*g(x)), Eq(g(x).diff(x), -x*f(x) + g(x))]
845 The system defined above is already in the desired form, so we do not have to convert it.
847 >>> (A1, A0), b = linear_ode_to_matrix(eqs, funcs, x, 1)
848 >>> A = A0
850 A user can also pass the commutative antiderivative required for type3 and type4 system of ODEs.
851 Passing an incorrect one will lead to incorrect results. If the coefficient matrix is not commutative
852 with its antiderivative, then :obj:`sympy.solvers.ode.systems.linodesolve_type()` raises a NotImplementedError.
853 If it does have a commutative antiderivative, then the function just returns the information about the system.
855 >>> system_info = linodesolve_type(A, x, b=b)
857 Now, we can pass the antiderivative as an argument to get the solution. If the system information is not
858 passed, then the solver will compute the required arguments internally.
860 >>> sol_vector = linodesolve(A, x, b=b)
862 Once again, we can verify the solution obtained.
864 >>> sol = [Eq(f, s) for f, s in zip(funcs, sol_vector)]
865 >>> checkodesol(eqs, sol)
866 (True, [0, 0])
868 Returns
869 =======
871 List
873 Raises
874 ======
876 ValueError
877 This error is raised when the coefficient matrix, non-homogeneous term
878 or the antiderivative, if passed, are not a matrix or
879 do not have correct dimensions
880 NonSquareMatrixError
881 When the coefficient matrix or its antiderivative, if passed is not a
882 square matrix
883 NotImplementedError
884 If the coefficient matrix does not have a commutative antiderivative
886 See Also
887 ========
889 linear_ode_to_matrix: Coefficient matrix computation function
890 canonical_odes: System of ODEs representation change
891 linodesolve_type: Getting information about systems of ODEs to pass in this solver
893 """
895 if not isinstance(A, MatrixBase):
896 raise ValueError(filldedent('''\
897 The coefficients of the system of ODEs should be of type Matrix
898 '''))
900 if not A.is_square:
901 raise NonSquareMatrixError(filldedent('''\
902 The coefficient matrix must be a square
903 '''))
905 if b is not None:
906 if not isinstance(b, MatrixBase):
907 raise ValueError(filldedent('''\
908 The non-homogeneous terms of the system of ODEs should be of type Matrix
909 '''))
911 if A.rows != b.rows:
912 raise ValueError(filldedent('''\
913 The system of ODEs should have the same number of non-homogeneous terms and the number of
914 equations
915 '''))
917 if B is not None:
918 if not isinstance(B, MatrixBase):
919 raise ValueError(filldedent('''\
920 The antiderivative of coefficients of the system of ODEs should be of type Matrix
921 '''))
923 if not B.is_square:
924 raise NonSquareMatrixError(filldedent('''\
925 The antiderivative of the coefficient matrix must be a square
926 '''))
928 if A.rows != B.rows:
929 raise ValueError(filldedent('''\
930 The coefficient matrix and its antiderivative should have same dimensions
931 '''))
933 if not any(type == "type{}".format(i) for i in range(1, 7)) and not type == "auto":
934 raise ValueError(filldedent('''\
935 The input type should be a valid one
936 '''))
938 n = A.rows
940 # constants = numbered_symbols(prefix='C', cls=Dummy, start=const_idx+1)
941 Cvect = Matrix([Dummy() for _ in range(n)])
943 if b is None and any(type == typ for typ in ["type2", "type4", "type6"]):
944 b = zeros(n, 1)
946 is_transformed = tau is not None
947 passed_type = type
949 if type == "auto":
950 system_info = linodesolve_type(A, t, b=b)
951 type = system_info["type_of_equation"]
952 B = system_info["antiderivative"]
954 if type in ("type5", "type6"):
955 is_transformed = True
956 if passed_type != "auto":
957 if tau is None:
958 system_info = _first_order_type5_6_subs(A, t, b=b)
959 if not system_info:
960 raise ValueError(filldedent('''
961 The system passed isn't {}.
962 '''.format(type)))
964 tau = system_info['tau']
965 t = system_info['t_']
966 A = system_info['A']
967 b = system_info['b']
969 intx_wrtt = lambda x: Integral(x, t) if x else 0
970 if type in ("type1", "type2", "type5", "type6"):
971 P, J = matrix_exp_jordan_form(A, t)
972 P = simplify(P)
974 if type in ("type1", "type5"):
975 sol_vector = P * (J * Cvect)
976 else:
977 Jinv = J.subs(t, -t)
978 sol_vector = P * J * ((Jinv * P.inv() * b).applyfunc(intx_wrtt) + Cvect)
979 else:
980 if B is None:
981 B, _ = _is_commutative_anti_derivative(A, t)
983 if type == "type3":
984 sol_vector = B.exp() * Cvect
985 else:
986 sol_vector = B.exp() * (((-B).exp() * b).applyfunc(intx_wrtt) + Cvect)
988 if is_transformed:
989 sol_vector = sol_vector.subs(t, tau)
991 gens = sol_vector.atoms(exp)
993 if type != "type1":
994 sol_vector = [expand_mul(s) for s in sol_vector]
996 sol_vector = [collect(s, ordered(gens), exact=True) for s in sol_vector]
998 if doit:
999 sol_vector = [s.doit() for s in sol_vector]
1001 return sol_vector
1004def _matrix_is_constant(M, t):
1005 """Checks if the matrix M is independent of t or not."""
1006 return all(coef.as_independent(t, as_Add=True)[1] == 0 for coef in M)
1009def canonical_odes(eqs, funcs, t):
1010 r"""
1011 Function that solves for highest order derivatives in a system
1013 Explanation
1014 ===========
1016 This function inputs a system of ODEs and based on the system,
1017 the dependent variables and their highest order, returns the system
1018 in the following form:
1020 .. math::
1021 X'(t) = A(t) X(t) + b(t)
1023 Here, $X(t)$ is the vector of dependent variables of lower order, $A(t)$ is
1024 the coefficient matrix, $b(t)$ is the non-homogeneous term and $X'(t)$ is the
1025 vector of dependent variables in their respective highest order. We use the term
1026 canonical form to imply the system of ODEs which is of the above form.
1028 If the system passed has a non-linear term with multiple solutions, then a list of
1029 systems is returned in its canonical form.
1031 Parameters
1032 ==========
1034 eqs : List
1035 List of the ODEs
1036 funcs : List
1037 List of dependent variables
1038 t : Symbol
1039 Independent variable
1041 Examples
1042 ========
1044 >>> from sympy import symbols, Function, Eq, Derivative
1045 >>> from sympy.solvers.ode.systems import canonical_odes
1046 >>> f, g = symbols("f g", cls=Function)
1047 >>> x, y = symbols("x y")
1048 >>> funcs = [f(x), g(x)]
1049 >>> eqs = [Eq(f(x).diff(x) - 7*f(x), 12*g(x)), Eq(g(x).diff(x) + g(x), 20*f(x))]
1051 >>> canonical_eqs = canonical_odes(eqs, funcs, x)
1052 >>> canonical_eqs
1053 [[Eq(Derivative(f(x), x), 7*f(x) + 12*g(x)), Eq(Derivative(g(x), x), 20*f(x) - g(x))]]
1055 >>> system = [Eq(Derivative(f(x), x)**2 - 2*Derivative(f(x), x) + 1, 4), Eq(-y*f(x) + Derivative(g(x), x), 0)]
1057 >>> canonical_system = canonical_odes(system, funcs, x)
1058 >>> canonical_system
1059 [[Eq(Derivative(f(x), x), -1), Eq(Derivative(g(x), x), y*f(x))], [Eq(Derivative(f(x), x), 3), Eq(Derivative(g(x), x), y*f(x))]]
1061 Returns
1062 =======
1064 List
1066 """
1067 from sympy.solvers.solvers import solve
1069 order = _get_func_order(eqs, funcs)
1071 canon_eqs = solve(eqs, *[func.diff(t, order[func]) for func in funcs], dict=True)
1073 systems = []
1074 for eq in canon_eqs:
1075 system = [Eq(func.diff(t, order[func]), eq[func.diff(t, order[func])]) for func in funcs]
1076 systems.append(system)
1078 return systems
1081def _is_commutative_anti_derivative(A, t):
1082 r"""
1083 Helper function for determining if the Matrix passed is commutative with its antiderivative
1085 Explanation
1086 ===========
1088 This function checks if the Matrix $A$ passed is commutative with its antiderivative with respect
1089 to the independent variable $t$.
1091 .. math::
1092 B(t) = \int A(t) dt
1094 The function outputs two values, first one being the antiderivative $B(t)$, second one being a
1095 boolean value, if True, then the matrix $A(t)$ passed is commutative with $B(t)$, else the matrix
1096 passed isn't commutative with $B(t)$.
1098 Parameters
1099 ==========
1101 A : Matrix
1102 The matrix which has to be checked
1103 t : Symbol
1104 Independent variable
1106 Examples
1107 ========
1109 >>> from sympy import symbols, Matrix
1110 >>> from sympy.solvers.ode.systems import _is_commutative_anti_derivative
1111 >>> t = symbols("t")
1112 >>> A = Matrix([[1, t], [-t, 1]])
1114 >>> B, is_commuting = _is_commutative_anti_derivative(A, t)
1115 >>> is_commuting
1116 True
1118 Returns
1119 =======
1121 Matrix, Boolean
1123 """
1124 B = integrate(A, t)
1125 is_commuting = (B*A - A*B).applyfunc(expand).applyfunc(factor_terms).is_zero_matrix
1127 is_commuting = False if is_commuting is None else is_commuting
1129 return B, is_commuting
1132def _factor_matrix(A, t):
1133 term = None
1134 for element in A:
1135 temp_term = element.as_independent(t)[1]
1136 if temp_term.has(t):
1137 term = temp_term
1138 break
1140 if term is not None:
1141 A_factored = (A/term).applyfunc(ratsimp)
1142 can_factor = _matrix_is_constant(A_factored, t)
1143 term = (term, A_factored) if can_factor else None
1145 return term
1148def _is_second_order_type2(A, t):
1149 term = _factor_matrix(A, t)
1150 is_type2 = False
1152 if term is not None:
1153 term = 1/term[0]
1154 is_type2 = term.is_polynomial()
1156 if is_type2:
1157 poly = Poly(term.expand(), t)
1158 monoms = poly.monoms()
1160 if monoms[0][0] in (2, 4):
1161 cs = _get_poly_coeffs(poly, 4)
1162 a, b, c, d, e = cs
1164 a1 = powdenest(sqrt(a), force=True)
1165 c1 = powdenest(sqrt(e), force=True)
1166 b1 = powdenest(sqrt(c - 2*a1*c1), force=True)
1168 is_type2 = (b == 2*a1*b1) and (d == 2*b1*c1)
1169 term = a1*t**2 + b1*t + c1
1171 else:
1172 is_type2 = False
1174 return is_type2, term
1177def _get_poly_coeffs(poly, order):
1178 cs = [0 for _ in range(order+1)]
1179 for c, m in zip(poly.coeffs(), poly.monoms()):
1180 cs[-1-m[0]] = c
1181 return cs
1184def _match_second_order_type(A1, A0, t, b=None):
1185 r"""
1186 Works only for second order system in its canonical form.
1188 Type 0: Constant coefficient matrix, can be simply solved by
1189 introducing dummy variables.
1190 Type 1: When the substitution: $U = t*X' - X$ works for reducing
1191 the second order system to first order system.
1192 Type 2: When the system is of the form: $poly * X'' = A*X$ where
1193 $poly$ is square of a quadratic polynomial with respect to
1194 *t* and $A$ is a constant coefficient matrix.
1196 """
1197 match = {"type_of_equation": "type0"}
1198 n = A1.shape[0]
1200 if _matrix_is_constant(A1, t) and _matrix_is_constant(A0, t):
1201 return match
1203 if (A1 + A0*t).applyfunc(expand_mul).is_zero_matrix:
1204 match.update({"type_of_equation": "type1", "A1": A1})
1206 elif A1.is_zero_matrix and (b is None or b.is_zero_matrix):
1207 is_type2, term = _is_second_order_type2(A0, t)
1208 if is_type2:
1209 a, b, c = _get_poly_coeffs(Poly(term, t), 2)
1210 A = (A0*(term**2).expand()).applyfunc(ratsimp) + (b**2/4 - a*c)*eye(n, n)
1211 tau = integrate(1/term, t)
1212 t_ = Symbol("{}_".format(t))
1213 match.update({"type_of_equation": "type2", "A0": A,
1214 "g(t)": sqrt(term), "tau": tau, "is_transformed": True,
1215 "t_": t_})
1217 return match
1220def _second_order_subs_type1(A, b, funcs, t):
1221 r"""
1222 For a linear, second order system of ODEs, a particular substitution.
1224 A system of the below form can be reduced to a linear first order system of
1225 ODEs:
1226 .. math::
1227 X'' = A(t) * (t*X' - X) + b(t)
1229 By substituting:
1230 .. math:: U = t*X' - X
1232 To get the system:
1233 .. math:: U' = t*(A(t)*U + b(t))
1235 Where $U$ is the vector of dependent variables, $X$ is the vector of dependent
1236 variables in `funcs` and $X'$ is the first order derivative of $X$ with respect to
1237 $t$. It may or may not reduce the system into linear first order system of ODEs.
1239 Then a check is made to determine if the system passed can be reduced or not, if
1240 this substitution works, then the system is reduced and its solved for the new
1241 substitution. After we get the solution for $U$:
1243 .. math:: U = a(t)
1245 We substitute and return the reduced system:
1247 .. math::
1248 a(t) = t*X' - X
1250 Parameters
1251 ==========
1253 A: Matrix
1254 Coefficient matrix($A(t)*t$) of the second order system of this form.
1255 b: Matrix
1256 Non-homogeneous term($b(t)$) of the system of ODEs.
1257 funcs: List
1258 List of dependent variables
1259 t: Symbol
1260 Independent variable of the system of ODEs.
1262 Returns
1263 =======
1265 List
1267 """
1269 U = Matrix([t*func.diff(t) - func for func in funcs])
1271 sol = linodesolve(A, t, t*b)
1272 reduced_eqs = [Eq(u, s) for s, u in zip(sol, U)]
1273 reduced_eqs = canonical_odes(reduced_eqs, funcs, t)[0]
1275 return reduced_eqs
1278def _second_order_subs_type2(A, funcs, t_):
1279 r"""
1280 Returns a second order system based on the coefficient matrix passed.
1282 Explanation
1283 ===========
1285 This function returns a system of second order ODE of the following form:
1287 .. math::
1288 X'' = A * X
1290 Here, $X$ is the vector of dependent variables, but a bit modified, $A$ is the
1291 coefficient matrix passed.
1293 Along with returning the second order system, this function also returns the new
1294 dependent variables with the new independent variable `t_` passed.
1296 Parameters
1297 ==========
1299 A: Matrix
1300 Coefficient matrix of the system
1301 funcs: List
1302 List of old dependent variables
1303 t_: Symbol
1304 New independent variable
1306 Returns
1307 =======
1309 List, List
1311 """
1312 func_names = [func.func.__name__ for func in funcs]
1313 new_funcs = [Function(Dummy("{}_".format(name)))(t_) for name in func_names]
1314 rhss = A * Matrix(new_funcs)
1315 new_eqs = [Eq(func.diff(t_, 2), rhs) for func, rhs in zip(new_funcs, rhss)]
1317 return new_eqs, new_funcs
1320def _is_euler_system(As, t):
1321 return all(_matrix_is_constant((A*t**i).applyfunc(ratsimp), t) for i, A in enumerate(As))
1324def _classify_linear_system(eqs, funcs, t, is_canon=False):
1325 r"""
1326 Returns a dictionary with details of the eqs if the system passed is linear
1327 and can be classified by this function else returns None
1329 Explanation
1330 ===========
1332 This function takes the eqs, converts it into a form Ax = b where x is a vector of terms
1333 containing dependent variables and their derivatives till their maximum order. If it is
1334 possible to convert eqs into Ax = b, then all the equations in eqs are linear otherwise
1335 they are non-linear.
1337 To check if the equations are constant coefficient, we need to check if all the terms in
1338 A obtained above are constant or not.
1340 To check if the equations are homogeneous or not, we need to check if b is a zero matrix
1341 or not.
1343 Parameters
1344 ==========
1346 eqs: List
1347 List of ODEs
1348 funcs: List
1349 List of dependent variables
1350 t: Symbol
1351 Independent variable of the equations in eqs
1352 is_canon: Boolean
1353 If True, then this function will not try to get the
1354 system in canonical form. Default value is False
1356 Returns
1357 =======
1359 match = {
1360 'no_of_equation': len(eqs),
1361 'eq': eqs,
1362 'func': funcs,
1363 'order': order,
1364 'is_linear': is_linear,
1365 'is_constant': is_constant,
1366 'is_homogeneous': is_homogeneous,
1367 }
1369 Dict or list of Dicts or None
1370 Dict with values for keys:
1371 1. no_of_equation: Number of equations
1372 2. eq: The set of equations
1373 3. func: List of dependent variables
1374 4. order: A dictionary that gives the order of the
1375 dependent variable in eqs
1376 5. is_linear: Boolean value indicating if the set of
1377 equations are linear or not.
1378 6. is_constant: Boolean value indicating if the set of
1379 equations have constant coefficients or not.
1380 7. is_homogeneous: Boolean value indicating if the set of
1381 equations are homogeneous or not.
1382 8. commutative_antiderivative: Antiderivative of the coefficient
1383 matrix if the coefficient matrix is non-constant
1384 and commutative with its antiderivative. This key
1385 may or may not exist.
1386 9. is_general: Boolean value indicating if the system of ODEs is
1387 solvable using one of the general case solvers or not.
1388 10. rhs: rhs of the non-homogeneous system of ODEs in Matrix form. This
1389 key may or may not exist.
1390 11. is_higher_order: True if the system passed has an order greater than 1.
1391 This key may or may not exist.
1392 12. is_second_order: True if the system passed is a second order ODE. This
1393 key may or may not exist.
1394 This Dict is the answer returned if the eqs are linear and constant
1395 coefficient. Otherwise, None is returned.
1397 """
1399 # Error for i == 0 can be added but isn't for now
1401 # Check for len(funcs) == len(eqs)
1402 if len(funcs) != len(eqs):
1403 raise ValueError("Number of functions given is not equal to the number of equations %s" % funcs)
1405 # ValueError when functions have more than one arguments
1406 for func in funcs:
1407 if len(func.args) != 1:
1408 raise ValueError("dsolve() and classify_sysode() work with "
1409 "functions of one variable only, not %s" % func)
1411 # Getting the func_dict and order using the helper
1412 # function
1413 order = _get_func_order(eqs, funcs)
1414 system_order = max(order[func] for func in funcs)
1415 is_higher_order = system_order > 1
1416 is_second_order = system_order == 2 and all(order[func] == 2 for func in funcs)
1418 # Not adding the check if the len(func.args) for
1419 # every func in funcs is 1
1421 # Linearity check
1422 try:
1424 canon_eqs = canonical_odes(eqs, funcs, t) if not is_canon else [eqs]
1425 if len(canon_eqs) == 1:
1426 As, b = linear_ode_to_matrix(canon_eqs[0], funcs, t, system_order)
1427 else:
1429 match = {
1430 'is_implicit': True,
1431 'canon_eqs': canon_eqs
1432 }
1434 return match
1436 # When the system of ODEs is non-linear, an ODENonlinearError is raised.
1437 # This function catches the error and None is returned.
1438 except ODENonlinearError:
1439 return None
1441 is_linear = True
1443 # Homogeneous check
1444 is_homogeneous = True if b.is_zero_matrix else False
1446 # Is general key is used to identify if the system of ODEs can be solved by
1447 # one of the general case solvers or not.
1448 match = {
1449 'no_of_equation': len(eqs),
1450 'eq': eqs,
1451 'func': funcs,
1452 'order': order,
1453 'is_linear': is_linear,
1454 'is_homogeneous': is_homogeneous,
1455 'is_general': True
1456 }
1458 if not is_homogeneous:
1459 match['rhs'] = b
1461 is_constant = all(_matrix_is_constant(A_, t) for A_ in As)
1463 # The match['is_linear'] check will be added in the future when this
1464 # function becomes ready to deal with non-linear systems of ODEs
1466 if not is_higher_order:
1467 A = As[1]
1468 match['func_coeff'] = A
1470 # Constant coefficient check
1471 is_constant = _matrix_is_constant(A, t)
1472 match['is_constant'] = is_constant
1474 try:
1475 system_info = linodesolve_type(A, t, b=b)
1476 except NotImplementedError:
1477 return None
1479 match.update(system_info)
1480 antiderivative = match.pop("antiderivative")
1482 if not is_constant:
1483 match['commutative_antiderivative'] = antiderivative
1485 return match
1486 else:
1487 match['type_of_equation'] = "type0"
1489 if is_second_order:
1490 A1, A0 = As[1:]
1492 match_second_order = _match_second_order_type(A1, A0, t)
1493 match.update(match_second_order)
1495 match['is_second_order'] = True
1497 # If system is constant, then no need to check if its in euler
1498 # form or not. It will be easier and faster to directly proceed
1499 # to solve it.
1500 if match['type_of_equation'] == "type0" and not is_constant:
1501 is_euler = _is_euler_system(As, t)
1502 if is_euler:
1503 t_ = Symbol('{}_'.format(t))
1504 match.update({'is_transformed': True, 'type_of_equation': 'type1',
1505 't_': t_})
1506 else:
1507 is_jordan = lambda M: M == Matrix.jordan_block(M.shape[0], M[0, 0])
1508 terms = _factor_matrix(As[-1], t)
1509 if all(A.is_zero_matrix for A in As[1:-1]) and terms is not None and not is_jordan(terms[1]):
1510 P, J = terms[1].jordan_form()
1511 match.update({'type_of_equation': 'type2', 'J': J,
1512 'f(t)': terms[0], 'P': P, 'is_transformed': True})
1514 if match['type_of_equation'] != 'type0' and is_second_order:
1515 match.pop('is_second_order', None)
1517 match['is_higher_order'] = is_higher_order
1519 return match
1521def _preprocess_eqs(eqs):
1522 processed_eqs = []
1523 for eq in eqs:
1524 processed_eqs.append(eq if isinstance(eq, Equality) else Eq(eq, 0))
1526 return processed_eqs
1529def _eqs2dict(eqs, funcs):
1530 eqsorig = {}
1531 eqsmap = {}
1532 funcset = set(funcs)
1533 for eq in eqs:
1534 f1, = eq.lhs.atoms(AppliedUndef)
1535 f2s = (eq.rhs.atoms(AppliedUndef) - {f1}) & funcset
1536 eqsmap[f1] = f2s
1537 eqsorig[f1] = eq
1538 return eqsmap, eqsorig
1541def _dict2graph(d):
1542 nodes = list(d)
1543 edges = [(f1, f2) for f1, f2s in d.items() for f2 in f2s]
1544 G = (nodes, edges)
1545 return G
1548def _is_type1(scc, t):
1549 eqs, funcs = scc
1551 try:
1552 (A1, A0), b = linear_ode_to_matrix(eqs, funcs, t, 1)
1553 except (ODENonlinearError, ODEOrderError):
1554 return False
1556 if _matrix_is_constant(A0, t) and b.is_zero_matrix:
1557 return True
1559 return False
1562def _combine_type1_subsystems(subsystem, funcs, t):
1563 indices = [i for i, sys in enumerate(zip(subsystem, funcs)) if _is_type1(sys, t)]
1564 remove = set()
1565 for ip, i in enumerate(indices):
1566 for j in indices[ip+1:]:
1567 if any(eq2.has(funcs[i]) for eq2 in subsystem[j]):
1568 subsystem[j] = subsystem[i] + subsystem[j]
1569 remove.add(i)
1570 subsystem = [sys for i, sys in enumerate(subsystem) if i not in remove]
1571 return subsystem
1574def _component_division(eqs, funcs, t):
1576 # Assuming that each eq in eqs is in canonical form,
1577 # that is, [f(x).diff(x) = .., g(x).diff(x) = .., etc]
1578 # and that the system passed is in its first order
1579 eqsmap, eqsorig = _eqs2dict(eqs, funcs)
1581 subsystems = []
1582 for cc in connected_components(_dict2graph(eqsmap)):
1583 eqsmap_c = {f: eqsmap[f] for f in cc}
1584 sccs = strongly_connected_components(_dict2graph(eqsmap_c))
1585 subsystem = [[eqsorig[f] for f in scc] for scc in sccs]
1586 subsystem = _combine_type1_subsystems(subsystem, sccs, t)
1587 subsystems.append(subsystem)
1589 return subsystems
1592# Returns: List of equations
1593def _linear_ode_solver(match):
1594 t = match['t']
1595 funcs = match['func']
1597 rhs = match.get('rhs', None)
1598 tau = match.get('tau', None)
1599 t = match['t_'] if 't_' in match else t
1600 A = match['func_coeff']
1602 # Note: To make B None when the matrix has constant
1603 # coefficient
1604 B = match.get('commutative_antiderivative', None)
1605 type = match['type_of_equation']
1607 sol_vector = linodesolve(A, t, b=rhs, B=B,
1608 type=type, tau=tau)
1610 sol = [Eq(f, s) for f, s in zip(funcs, sol_vector)]
1612 return sol
1615def _select_equations(eqs, funcs, key=lambda x: x):
1616 eq_dict = {e.lhs: e.rhs for e in eqs}
1617 return [Eq(f, eq_dict[key(f)]) for f in funcs]
1620def _higher_order_ode_solver(match):
1621 eqs = match["eq"]
1622 funcs = match["func"]
1623 t = match["t"]
1624 sysorder = match['order']
1625 type = match.get('type_of_equation', "type0")
1627 is_second_order = match.get('is_second_order', False)
1628 is_transformed = match.get('is_transformed', False)
1629 is_euler = is_transformed and type == "type1"
1630 is_higher_order_type2 = is_transformed and type == "type2" and 'P' in match
1632 if is_second_order:
1633 new_eqs, new_funcs = _second_order_to_first_order(eqs, funcs, t,
1634 A1=match.get("A1", None), A0=match.get("A0", None),
1635 b=match.get("rhs", None), type=type,
1636 t_=match.get("t_", None))
1637 else:
1638 new_eqs, new_funcs = _higher_order_to_first_order(eqs, sysorder, t, funcs=funcs,
1639 type=type, J=match.get('J', None),
1640 f_t=match.get('f(t)', None),
1641 P=match.get('P', None), b=match.get('rhs', None))
1643 if is_transformed:
1644 t = match.get('t_', t)
1646 if not is_higher_order_type2:
1647 new_eqs = _select_equations(new_eqs, [f.diff(t) for f in new_funcs])
1649 sol = None
1651 # NotImplementedError may be raised when the system may be actually
1652 # solvable if it can be just divided into sub-systems
1653 try:
1654 if not is_higher_order_type2:
1655 sol = _strong_component_solver(new_eqs, new_funcs, t)
1656 except NotImplementedError:
1657 sol = None
1659 # Dividing the system only when it becomes essential
1660 if sol is None:
1661 try:
1662 sol = _component_solver(new_eqs, new_funcs, t)
1663 except NotImplementedError:
1664 sol = None
1666 if sol is None:
1667 return sol
1669 is_second_order_type2 = is_second_order and type == "type2"
1671 underscores = '__' if is_transformed else '_'
1673 sol = _select_equations(sol, funcs,
1674 key=lambda x: Function(Dummy('{}{}0'.format(x.func.__name__, underscores)))(t))
1676 if match.get("is_transformed", False):
1677 if is_second_order_type2:
1678 g_t = match["g(t)"]
1679 tau = match["tau"]
1680 sol = [Eq(s.lhs, s.rhs.subs(t, tau) * g_t) for s in sol]
1681 elif is_euler:
1682 t = match['t']
1683 tau = match['t_']
1684 sol = [s.subs(tau, log(t)) for s in sol]
1685 elif is_higher_order_type2:
1686 P = match['P']
1687 sol_vector = P * Matrix([s.rhs for s in sol])
1688 sol = [Eq(f, s) for f, s in zip(funcs, sol_vector)]
1690 return sol
1693# Returns: List of equations or None
1694# If None is returned by this solver, then the system
1695# of ODEs cannot be solved directly by dsolve_system.
1696def _strong_component_solver(eqs, funcs, t):
1697 from sympy.solvers.ode.ode import dsolve, constant_renumber
1699 match = _classify_linear_system(eqs, funcs, t, is_canon=True)
1700 sol = None
1702 # Assuming that we can't get an implicit system
1703 # since we are already canonical equations from
1704 # dsolve_system
1705 if match:
1706 match['t'] = t
1708 if match.get('is_higher_order', False):
1709 sol = _higher_order_ode_solver(match)
1711 elif match.get('is_linear', False):
1712 sol = _linear_ode_solver(match)
1714 # Note: For now, only linear systems are handled by this function
1715 # hence, the match condition is added. This can be removed later.
1716 if sol is None and len(eqs) == 1:
1717 sol = dsolve(eqs[0], func=funcs[0])
1718 variables = Tuple(eqs[0]).free_symbols
1719 new_constants = [Dummy() for _ in range(ode_order(eqs[0], funcs[0]))]
1720 sol = constant_renumber(sol, variables=variables, newconstants=new_constants)
1721 sol = [sol]
1723 # To add non-linear case here in future
1725 return sol
1728def _get_funcs_from_canon(eqs):
1729 return [eq.lhs.args[0] for eq in eqs]
1732# Returns: List of Equations(a solution)
1733def _weak_component_solver(wcc, t):
1735 # We will divide the systems into sccs
1736 # only when the wcc cannot be solved as
1737 # a whole
1738 eqs = []
1739 for scc in wcc:
1740 eqs += scc
1741 funcs = _get_funcs_from_canon(eqs)
1743 sol = _strong_component_solver(eqs, funcs, t)
1744 if sol:
1745 return sol
1747 sol = []
1749 for j, scc in enumerate(wcc):
1750 eqs = scc
1751 funcs = _get_funcs_from_canon(eqs)
1753 # Substituting solutions for the dependent
1754 # variables solved in previous SCC, if any solved.
1755 comp_eqs = [eq.subs({s.lhs: s.rhs for s in sol}) for eq in eqs]
1756 scc_sol = _strong_component_solver(comp_eqs, funcs, t)
1758 if scc_sol is None:
1759 raise NotImplementedError(filldedent('''
1760 The system of ODEs passed cannot be solved by dsolve_system.
1761 '''))
1763 # scc_sol: List of equations
1764 # scc_sol is a solution
1765 sol += scc_sol
1767 return sol
1770# Returns: List of Equations(a solution)
1771def _component_solver(eqs, funcs, t):
1772 components = _component_division(eqs, funcs, t)
1773 sol = []
1775 for wcc in components:
1777 # wcc_sol: List of Equations
1778 sol += _weak_component_solver(wcc, t)
1780 # sol: List of Equations
1781 return sol
1784def _second_order_to_first_order(eqs, funcs, t, type="auto", A1=None,
1785 A0=None, b=None, t_=None):
1786 r"""
1787 Expects the system to be in second order and in canonical form
1789 Explanation
1790 ===========
1792 Reduces a second order system into a first order one depending on the type of second
1793 order system.
1794 1. "type0": If this is passed, then the system will be reduced to first order by
1795 introducing dummy variables.
1796 2. "type1": If this is passed, then a particular substitution will be used to reduce the
1797 the system into first order.
1798 3. "type2": If this is passed, then the system will be transformed with new dependent
1799 variables and independent variables. This transformation is a part of solving
1800 the corresponding system of ODEs.
1802 `A1` and `A0` are the coefficient matrices from the system and it is assumed that the
1803 second order system has the form given below:
1805 .. math::
1806 A2 * X'' = A1 * X' + A0 * X + b
1808 Here, $A2$ is the coefficient matrix for the vector $X''$ and $b$ is the non-homogeneous
1809 term.
1811 Default value for `b` is None but if `A1` and `A0` are passed and `b` is not passed, then the
1812 system will be assumed homogeneous.
1814 """
1815 is_a1 = A1 is None
1816 is_a0 = A0 is None
1818 if (type == "type1" and is_a1) or (type == "type2" and is_a0)\
1819 or (type == "auto" and (is_a1 or is_a0)):
1820 (A2, A1, A0), b = linear_ode_to_matrix(eqs, funcs, t, 2)
1822 if not A2.is_Identity:
1823 raise ValueError(filldedent('''
1824 The system must be in its canonical form.
1825 '''))
1827 if type == "auto":
1828 match = _match_second_order_type(A1, A0, t)
1829 type = match["type_of_equation"]
1830 A1 = match.get("A1", None)
1831 A0 = match.get("A0", None)
1833 sys_order = {func: 2 for func in funcs}
1835 if type == "type1":
1836 if b is None:
1837 b = zeros(len(eqs))
1838 eqs = _second_order_subs_type1(A1, b, funcs, t)
1839 sys_order = {func: 1 for func in funcs}
1841 if type == "type2":
1842 if t_ is None:
1843 t_ = Symbol("{}_".format(t))
1844 t = t_
1845 eqs, funcs = _second_order_subs_type2(A0, funcs, t_)
1846 sys_order = {func: 2 for func in funcs}
1848 return _higher_order_to_first_order(eqs, sys_order, t, funcs=funcs)
1851def _higher_order_type2_to_sub_systems(J, f_t, funcs, t, max_order, b=None, P=None):
1853 # Note: To add a test for this ValueError
1854 if J is None or f_t is None or not _matrix_is_constant(J, t):
1855 raise ValueError(filldedent('''
1856 Correctly input for args 'A' and 'f_t' for Linear, Higher Order,
1857 Type 2
1858 '''))
1860 if P is None and b is not None and not b.is_zero_matrix:
1861 raise ValueError(filldedent('''
1862 Provide the keyword 'P' for matrix P in A = P * J * P-1.
1863 '''))
1865 new_funcs = Matrix([Function(Dummy('{}__0'.format(f.func.__name__)))(t) for f in funcs])
1866 new_eqs = new_funcs.diff(t, max_order) - f_t * J * new_funcs
1868 if b is not None and not b.is_zero_matrix:
1869 new_eqs -= P.inv() * b
1871 new_eqs = canonical_odes(new_eqs, new_funcs, t)[0]
1873 return new_eqs, new_funcs
1876def _higher_order_to_first_order(eqs, sys_order, t, funcs=None, type="type0", **kwargs):
1877 if funcs is None:
1878 funcs = sys_order.keys()
1880 # Standard Cauchy Euler system
1881 if type == "type1":
1882 t_ = Symbol('{}_'.format(t))
1883 new_funcs = [Function(Dummy('{}_'.format(f.func.__name__)))(t_) for f in funcs]
1884 max_order = max(sys_order[func] for func in funcs)
1885 subs_dict = {func: new_func for func, new_func in zip(funcs, new_funcs)}
1886 subs_dict[t] = exp(t_)
1888 free_function = Function(Dummy())
1890 def _get_coeffs_from_subs_expression(expr):
1891 if isinstance(expr, Subs):
1892 free_symbol = expr.args[1][0]
1893 term = expr.args[0]
1894 return {ode_order(term, free_symbol): 1}
1896 if isinstance(expr, Mul):
1897 coeff = expr.args[0]
1898 order = list(_get_coeffs_from_subs_expression(expr.args[1]).keys())[0]
1899 return {order: coeff}
1901 if isinstance(expr, Add):
1902 coeffs = {}
1903 for arg in expr.args:
1905 if isinstance(arg, Mul):
1906 coeffs.update(_get_coeffs_from_subs_expression(arg))
1908 else:
1909 order = list(_get_coeffs_from_subs_expression(arg).keys())[0]
1910 coeffs[order] = 1
1912 return coeffs
1914 for o in range(1, max_order + 1):
1915 expr = free_function(log(t_)).diff(t_, o)*t_**o
1916 coeff_dict = _get_coeffs_from_subs_expression(expr)
1917 coeffs = [coeff_dict[order] if order in coeff_dict else 0 for order in range(o + 1)]
1918 expr_to_subs = sum(free_function(t_).diff(t_, i) * c for i, c in
1919 enumerate(coeffs)) / t**o
1920 subs_dict.update({f.diff(t, o): expr_to_subs.subs(free_function(t_), nf)
1921 for f, nf in zip(funcs, new_funcs)})
1923 new_eqs = [eq.subs(subs_dict) for eq in eqs]
1924 new_sys_order = {nf: sys_order[f] for f, nf in zip(funcs, new_funcs)}
1926 new_eqs = canonical_odes(new_eqs, new_funcs, t_)[0]
1928 return _higher_order_to_first_order(new_eqs, new_sys_order, t_, funcs=new_funcs)
1930 # Systems of the form: X(n)(t) = f(t)*A*X + b
1931 # where X(n)(t) is the nth derivative of the vector of dependent variables
1932 # with respect to the independent variable and A is a constant matrix.
1933 if type == "type2":
1934 J = kwargs.get('J', None)
1935 f_t = kwargs.get('f_t', None)
1936 b = kwargs.get('b', None)
1937 P = kwargs.get('P', None)
1938 max_order = max(sys_order[func] for func in funcs)
1940 return _higher_order_type2_to_sub_systems(J, f_t, funcs, t, max_order, P=P, b=b)
1942 # Note: To be changed to this after doit option is disabled for default cases
1943 # new_sysorder = _get_func_order(new_eqs, new_funcs)
1944 #
1945 # return _higher_order_to_first_order(new_eqs, new_sysorder, t, funcs=new_funcs)
1947 new_funcs = []
1949 for prev_func in funcs:
1950 func_name = prev_func.func.__name__
1951 func = Function(Dummy('{}_0'.format(func_name)))(t)
1952 new_funcs.append(func)
1953 subs_dict = {prev_func: func}
1954 new_eqs = []
1956 for i in range(1, sys_order[prev_func]):
1957 new_func = Function(Dummy('{}_{}'.format(func_name, i)))(t)
1958 subs_dict[prev_func.diff(t, i)] = new_func
1959 new_funcs.append(new_func)
1961 prev_f = subs_dict[prev_func.diff(t, i-1)]
1962 new_eq = Eq(prev_f.diff(t), new_func)
1963 new_eqs.append(new_eq)
1965 eqs = [eq.subs(subs_dict) for eq in eqs] + new_eqs
1967 return eqs, new_funcs
1970def dsolve_system(eqs, funcs=None, t=None, ics=None, doit=False, simplify=True):
1971 r"""
1972 Solves any(supported) system of Ordinary Differential Equations
1974 Explanation
1975 ===========
1977 This function takes a system of ODEs as an input, determines if the
1978 it is solvable by this function, and returns the solution if found any.
1980 This function can handle:
1981 1. Linear, First Order, Constant coefficient homogeneous system of ODEs
1982 2. Linear, First Order, Constant coefficient non-homogeneous system of ODEs
1983 3. Linear, First Order, non-constant coefficient homogeneous system of ODEs
1984 4. Linear, First Order, non-constant coefficient non-homogeneous system of ODEs
1985 5. Any implicit system which can be divided into system of ODEs which is of the above 4 forms
1986 6. Any higher order linear system of ODEs that can be reduced to one of the 5 forms of systems described above.
1988 The types of systems described above are not limited by the number of equations, i.e. this
1989 function can solve the above types irrespective of the number of equations in the system passed.
1990 But, the bigger the system, the more time it will take to solve the system.
1992 This function returns a list of solutions. Each solution is a list of equations where LHS is
1993 the dependent variable and RHS is an expression in terms of the independent variable.
1995 Among the non constant coefficient types, not all the systems are solvable by this function. Only
1996 those which have either a coefficient matrix with a commutative antiderivative or those systems which
1997 may be divided further so that the divided systems may have coefficient matrix with commutative antiderivative.
1999 Parameters
2000 ==========
2002 eqs : List
2003 system of ODEs to be solved
2004 funcs : List or None
2005 List of dependent variables that make up the system of ODEs
2006 t : Symbol or None
2007 Independent variable in the system of ODEs
2008 ics : Dict or None
2009 Set of initial boundary/conditions for the system of ODEs
2010 doit : Boolean
2011 Evaluate the solutions if True. Default value is True. Can be
2012 set to false if the integral evaluation takes too much time and/or
2013 is not required.
2014 simplify: Boolean
2015 Simplify the solutions for the systems. Default value is True.
2016 Can be set to false if simplification takes too much time and/or
2017 is not required.
2019 Examples
2020 ========
2022 >>> from sympy import symbols, Eq, Function
2023 >>> from sympy.solvers.ode.systems import dsolve_system
2024 >>> f, g = symbols("f g", cls=Function)
2025 >>> x = symbols("x")
2027 >>> eqs = [Eq(f(x).diff(x), g(x)), Eq(g(x).diff(x), f(x))]
2028 >>> dsolve_system(eqs)
2029 [[Eq(f(x), -C1*exp(-x) + C2*exp(x)), Eq(g(x), C1*exp(-x) + C2*exp(x))]]
2031 You can also pass the initial conditions for the system of ODEs:
2033 >>> dsolve_system(eqs, ics={f(0): 1, g(0): 0})
2034 [[Eq(f(x), exp(x)/2 + exp(-x)/2), Eq(g(x), exp(x)/2 - exp(-x)/2)]]
2036 Optionally, you can pass the dependent variables and the independent
2037 variable for which the system is to be solved:
2039 >>> funcs = [f(x), g(x)]
2040 >>> dsolve_system(eqs, funcs=funcs, t=x)
2041 [[Eq(f(x), -C1*exp(-x) + C2*exp(x)), Eq(g(x), C1*exp(-x) + C2*exp(x))]]
2043 Lets look at an implicit system of ODEs:
2045 >>> eqs = [Eq(f(x).diff(x)**2, g(x)**2), Eq(g(x).diff(x), g(x))]
2046 >>> dsolve_system(eqs)
2047 [[Eq(f(x), C1 - C2*exp(x)), Eq(g(x), C2*exp(x))], [Eq(f(x), C1 + C2*exp(x)), Eq(g(x), C2*exp(x))]]
2049 Returns
2050 =======
2052 List of List of Equations
2054 Raises
2055 ======
2057 NotImplementedError
2058 When the system of ODEs is not solvable by this function.
2059 ValueError
2060 When the parameters passed are not in the required form.
2062 """
2063 from sympy.solvers.ode.ode import solve_ics, _extract_funcs, constant_renumber
2065 if not iterable(eqs):
2066 raise ValueError(filldedent('''
2067 List of equations should be passed. The input is not valid.
2068 '''))
2070 eqs = _preprocess_eqs(eqs)
2072 if funcs is not None and not isinstance(funcs, list):
2073 raise ValueError(filldedent('''
2074 Input to the funcs should be a list of functions.
2075 '''))
2077 if funcs is None:
2078 funcs = _extract_funcs(eqs)
2080 if any(len(func.args) != 1 for func in funcs):
2081 raise ValueError(filldedent('''
2082 dsolve_system can solve a system of ODEs with only one independent
2083 variable.
2084 '''))
2086 if len(eqs) != len(funcs):
2087 raise ValueError(filldedent('''
2088 Number of equations and number of functions do not match
2089 '''))
2091 if t is not None and not isinstance(t, Symbol):
2092 raise ValueError(filldedent('''
2093 The independent variable must be of type Symbol
2094 '''))
2096 if t is None:
2097 t = list(list(eqs[0].atoms(Derivative))[0].atoms(Symbol))[0]
2099 sols = []
2100 canon_eqs = canonical_odes(eqs, funcs, t)
2102 for canon_eq in canon_eqs:
2103 try:
2104 sol = _strong_component_solver(canon_eq, funcs, t)
2105 except NotImplementedError:
2106 sol = None
2108 if sol is None:
2109 sol = _component_solver(canon_eq, funcs, t)
2111 sols.append(sol)
2113 if sols:
2114 final_sols = []
2115 variables = Tuple(*eqs).free_symbols
2117 for sol in sols:
2119 sol = _select_equations(sol, funcs)
2120 sol = constant_renumber(sol, variables=variables)
2122 if ics:
2123 constants = Tuple(*sol).free_symbols - variables
2124 solved_constants = solve_ics(sol, funcs, constants, ics)
2125 sol = [s.subs(solved_constants) for s in sol]
2127 if simplify:
2128 constants = Tuple(*sol).free_symbols - variables
2129 sol = simpsol(sol, [t], constants, doit=doit)
2131 final_sols.append(sol)
2133 sols = final_sols
2135 return sols