-
-
Notifications
You must be signed in to change notification settings - Fork 4.4k
ODE Systems roadmap
This is written by Oscar Benjamin in Feb 2020 and represents my opinion of the state of SymPy at that time with a view to improvements that should be made. These would be suitable for a 2020 GSOC project but are also useful as a roadmap for the ode module in any case.
An ODE system is a system of ODEs for multiple dependent variables. Currently SymPy can solve a variety of systems of ODEs through the dsolve function e.g.:
In [12]: t = Symbol('t')
In [13]: x = Function('x')
In [14]: y = Function('y')
In [15]: eqs = [Eq(x(t).diff(t), y(t)), Eq(y(t).diff(t), x(t))]
In [16]: eqs
Out[16]:
⎡d d ⎤
⎢──(x(t)) = y(t), ──(y(t)) = x(t)⎥
⎣dt dt ⎦
In [17]: dsolve(eqs, [x(t), y(t)])
Out[17]:
⎡ t -t t -t⎤
⎣x(t) = C₁⋅ℯ + C₂⋅ℯ , y(t) = C₁⋅ℯ - C₂⋅ℯ ⎦Although some simple cases work there are significant problems with the code for handling of systems of ODEs in dsolve:
- The code is badly organised and should be refactored.
- Many basic trivial cases are not currently handled.
- There are few general purpose solvers.
- Most of the current solvers are special cases of more general solvers that could be implemented.
- The solvers that are not special cases of more general solvers are for obscure systems that are unlikely to come up in practice.
- Much of the code is untested and broken.
Combining all of these points I would argue that none of the existing code for solving systems of ODEs in dsolve needs to be kept. The only code that is worth keeping is the tests and probably those should end up being changed to some extent while doing this work. That does not mean that we should start by deleting that code: we should add new refactored general purpose solvers until the old code is not needed any more.
For an example of what I mean by refactoring the solving code see #18403 which does a similar refactor. That is for the single ODE case but we should do something similar for the systems code as well.
Currently dsolve for systems of ODEs consists mostly of solvers that handle systems of 2 or 3 equations. There is only one n-equation solver that can handle an arbitrary number of equations. It is for systems with constant coefficients and does not correctly handle defective matrices (see below). We should try to do everything in dsolve for systems using n-equation solvers. I will describe below a number of methods that can work for n-equations that should be implemented in dsolve.
I don't know if there are good textbooks for this but the following wikipedia articles would be useful for anyone looking to work on this:
https://en.wikipedia.org/wiki/Matrix_exponential
https://en.wikipedia.org/wiki/Matrix_differential_equation
https://en.wikipedia.org/wiki/Diagonalizable_matrix
https://en.wikipedia.org/wiki/Defective_matrix
https://en.wikipedia.org/wiki/Jordan_normal_form
https://en.wikipedia.org/wiki/Generalized_eigenvector
https://en.wikipedia.org/wiki/Integrating_factor
https://en.wikipedia.org/wiki/Connectivity_(graph_theory)
https://en.wikipedia.org/wiki/Strongly_connected_component
Most systems of ODEs do not have any analytic solution. SymPy's goal is to provide analytic solutions for the small number of cases where they exist. The current ODE systems code attempts to enumerate many different particular examples of systems and provide their solutions. This is a bad approach that leads to a large amount of inconsistent and buggy code. It would be better to have a small number of general solvers that can handle wide classes of problems that can be expected to be solvable.
The most basic class of systems of ODEs for which analytic solutions can be found is that of constant coefficient homogeneous systems of ODEs. For example the general 2x2 system of this class would look like
In [22]: a, b, c, d = symbols('a, b, c, d')
In [23]: eqs = [Eq(x(t).diff(t), a*x(t) + b*y(t)), Eq(y(t).diff(t), c*x(t) + d*y(t))]
In [24]: eqs
Out[24]:
⎡d d ⎤
⎢──(x(t)) = a⋅x(t) + b⋅y(t), ──(y(t)) = c⋅x(t) + d⋅y(t)⎥
⎣dt dt ⎦It is useful to think of such a system by grouping together the variables x and y in a vector X := [x, y]:
In [26]: eq = Eq(Matrix([[x(t).diff(t)], [y(t).diff(t)]]), Matrix([[a*x(t)+b*y(t)], [c*x(t)+d*y(t)]]))
In [27]: eq
Out[27]:
⎡d ⎤
⎢──(x(t))⎥
⎢dt ⎥ ⎡a⋅x(t) + b⋅y(t)⎤
⎢ ⎥ = ⎢ ⎥
⎢d ⎥ ⎣c⋅x(t) + d⋅y(t)⎦
⎢──(y(t))⎥
⎣dt ⎦ We can then pull out the derivative and factor the rhs into a matrix vector product
In [36]: eq = Eq(Derivative(Matrix([[x(t)], [y(t)]])), MatMul(Matrix([[a, b], [c, d]]), Matrix([[x(t)], [y(t)]])), evaluate=False)
In [37]: eq
Out[37]:
d ⎛⎡x(t)⎤⎞ ⎡a b⎤ ⎡x(t)⎤
──⎜⎢ ⎥⎟ = ⎢ ⎥⋅⎢ ⎥
dt⎝⎣y(t)⎦⎠ ⎣c d⎦ ⎣y(t)⎦Now finally this is in the form of the matrix ODE
X' = A * X (1)
where X is the vector [x, y], X' is its derivative wrt t and A is a square constant matrix. We can use a trial solution X = V exp(m t) where V is a constant vector and m is a constant. Substituting the trial solution into (1) gives
X' = A X --> V m exp(m t) = A V exp(m t) --> A V = m V
which is the equation for an eigenvector and eigenvalues of the matrix A. What this shows is that X = V exp(m t) is a solution of (1) if V is an eigenvector of A and m is the corresponding eigenvalue. Given the full set of eigenvectors V1, V2, ... and eigenvalues m1, m2, ... the general solution of (1) is given as
X = C1 V1 exp(m1 t) + C2 V2 exp(m2 t) + ...
where C1, C2, ... are the arbitrary constants of integration.
Taking the example from the top of this document we can find the solution using this method
In [50]: eqs = [Eq(x(t).diff(t), y(t)), Eq(y(t).diff(t), x(t))]
In [51]: eqs
Out[51]:
⎡d d ⎤
⎢──(x(t)) = y(t), ──(y(t)) = x(t)⎥
⎣dt dt ⎦
In [52]: M = Matrix([[0, 1], [1, 0]])
In [53]: M
Out[53]:
⎡0 1⎤
⎢ ⎥
⎣1 0⎦
In [54]: M.eigenvects()
Out[54]:
⎡⎛ ⎡⎡-1⎤⎤⎞ ⎛ ⎡⎡1⎤⎤⎞⎤
⎢⎜-1, 1, ⎢⎢ ⎥⎥⎟, ⎜1, 1, ⎢⎢ ⎥⎥⎟⎥
⎣⎝ ⎣⎣1 ⎦⎦⎠ ⎝ ⎣⎣1⎦⎦⎠⎦
In [55]: (m2, _, (V2,)), (m1, _, (V1,)) = M.eigenvects()
In [56]: gensol = C1*V1*exp(m1*t) + C2*V2*exp(m2*t)
In [57]: gensol
Out[57]:
⎡ t -t⎤
⎢C₁⋅ℯ - C₂⋅ℯ ⎥
⎢ ⎥
⎢ t -t⎥
⎣C₁⋅ℯ + C₂⋅ℯ ⎦
In [58]: sol = [Eq(x(t), gensol[0]), Eq(y(t), gensol[1])]
In [59]: sol
Out[59]:
⎡ t -t t -t⎤
⎣x(t) = C₁⋅ℯ - C₂⋅ℯ , y(t) = C₁⋅ℯ + C₂⋅ℯ ⎦This is the same as the solution given by dsolve apart from the sign of the arbitrary constant C2 (the sign of the constant is arbitrary so this is equally correct).
The above approach works well for diagonalisable matrices. An N x N diagonalisable matrix has N eigenvectors and corresponding eigenvalues and the general solution as formulated above will have N arbitrary constants as required for the full general solution. However when the characteristic polynomial of the matrix has repeated roots it is possible that the matrix will be defective. In this case the matrix does not have a full set of eigenvectors and it is not possible to construct the general solution as above.
Extending the approach above to handle defective matrices is possible but requires using generalised eigenvectors. I don't think that this is the right approach to take for dsolve but I will briefly illustrate how it works. Suppose that we have an ODE system with a 2 x 2 defective matrix such as
In [60]: eqs = [Eq(x(t).diff(t), x(t) + y(t)), Eq(y(t).diff(t), y(t))]
In [61]: eqs
Out[61]:
⎡d d ⎤
⎢──(x(t)) = x(t) + y(t), ──(y(t)) = y(t)⎥
⎣dt dt ⎦Here the associated defective matrix is
In [62]: M = Matrix([[1, 1], [0, 1]])
In [63]: M
Out[63]:
⎡1 1⎤
⎢ ⎥
⎣0 1⎦
In [64]: M.eigenvects()
Out[64]:
⎡⎛ ⎡⎡1⎤⎤⎞⎤
⎢⎜1, 2, ⎢⎢ ⎥⎥⎟⎥
⎣⎝ ⎣⎣0⎦⎦⎠⎦We see that the matrix only has one eigenvector V1 even though it is a 2 x 2 matrix. The matrix has one eigenvalue with algebraic multiplicity two but geometric multiplicity one and there is only one eigenvector for that eigenvalue.
Going back to (1) we can use a trial solution of the form X = (V + b t W) exp(m t) where b and m are constant scalars and V and W are constant vectors. Substituting this into (1) we find
X' = A X --> b W exp(m t) + m (V + b t W) exp(m t) = A (V + b t W) exp(m t)
--> b W + m (V + b t W) = A (V + b t W)
--> (b W + m V) + t (b m W) = A V + t (b A W)
Equating coefficients of t gives
A W = m W
so W needs to be an eigenvector of A with eigenvalue m. Equating constant coefficients gives
A V = m V + b W --> A V - m V = b W
which shows that V should be a generalised eigenvector for the eigenvalue m and constrains b for a particular choice of V and W.
The general solution for the defective matrix ODE above can be found with
In [87]: eqs = [Eq(x(t).diff(t), x(t) + y(t)), Eq(y(t).diff(t), y(t))]
In [88]: eqs
Out[88]:
⎡d d ⎤
⎢──(x(t)) = x(t) + y(t), ──(y(t)) = y(t)⎥
⎣dt dt ⎦
In [89]: m = 1
In [90]: W = Matrix([[1], [0]])
In [91]: V = Matrix([[0], [1]])
In [92]: C1*W*exp(m*t) + C2*(V+t*W)*exp(m*t)
Out[92]:
⎡ t t⎤
⎢C₁⋅ℯ + C₂⋅t⋅ℯ ⎥
⎢ ⎥
⎢ t ⎥
⎣ C₂⋅ℯ ⎦This solution is found correctly by dsolve:
In [93]: dsolve(eqs)
Out[93]:
⎡ t t⎤
⎣x(t) = (C₁ + C₂⋅t)⋅ℯ , y(t) = C₂⋅ℯ ⎦A better approach for handling constant coefficient homogeneous systems of ODEs is to use the matrix exponential. We define the exponential of a matrix using the Maclaurin series for exp as
exp(B) = I + B + B^2/2 + B^3/3! + ... + B^n/n! + ...
We can compute this in SymPy with
In [98]: M = Matrix([[1, 1], [0, 2]])
In [99]: M
Out[99]:
⎡1 1⎤
⎢ ⎥
⎣0 2⎦
In [100]: M.exp()
Out[100]:
⎡ 2⎤
⎢ℯ -ℯ + ℯ ⎥
⎢ ⎥
⎢ 2 ⎥
⎣0 ℯ ⎦Given the matrix ODE (1) (X' = A X) we can write the general solution as
X = exp(A t) C
where exp(A t) is a matrix and C is a vector of arbitrary constants. If the initial conditions are X(0) = X0 then the solution becomes
X = exp(A t) X0
Computing the matrix exponential of a diagonalisable matrix is as simple as diagonalising the matrix and then exponentiating the eigenvalues on the diagonal:
In [125]: M.diagonalize()
Out[125]:
⎛⎡1 1⎤ ⎡1 0⎤⎞
⎜⎢ ⎥, ⎢ ⎥⎟
⎝⎣0 1⎦ ⎣0 2⎦⎠
In [126]: P, D = M.diagonalize()
In [127]: P * D * P.inv()
Out[127]:
⎡1 1⎤
⎢ ⎥
⎣0 2⎦
In [128]: P * D.exp() * P.inv()
Out[128]:
⎡ 2⎤
⎢ℯ -ℯ + ℯ ⎥
⎢ ⎥
⎢ 2 ⎥
⎣0 ℯ ⎦Of course this works for the case of diagonalisable (non-defective) matrices but to handle this in general we should use the Jordan normal form. Given a defective matrix there is no diagonalisation e.g.
In [129]: M = Matrix([[1, 1], [0, 1]])
In [130]: M.diagonalize()
---------------------------------------------------------------------------
MatrixError: Matrix is not diagonalizableHowever there is a Jordan normal form
In [131]: M.jordan_form()
Out[131]:
⎛⎡1 0⎤ ⎡1 1⎤⎞
⎜⎢ ⎥, ⎢ ⎥⎟
⎝⎣0 1⎦ ⎣0 1⎦⎠
In [132]: P, J = M.jordan_form()
In [133]: P * J * P.inv()
Out[133]:
⎡1 1⎤
⎢ ⎥
⎣0 1⎦Unlike the diagonalisation the Jordan normal form exists for any square matrix over the complex numbers. The Jordan normal form decomposes into a matrix product where instead of a diagonal matrix we have a block diagonal matrix where the blocks are of a particular form known as Jordan blocks. If we know how to compute the matrix exponential of the Jordan blocks we can compute the matrix exponential of any square matrix
In [134]: P * J.exp() * P.inv()
Out[134]:
⎡ℯ ℯ⎤
⎢ ⎥
⎣0 ℯ⎦
In [135]: M.exp()
Out[135]:
⎡ℯ ℯ⎤
⎢ ⎥
⎣0 ℯ⎦This provides a way to solve any constant coefficient system of ODEs. The one possible disadvantage of this approach is that it can lead to results involving complex numbers for a system of real ODEs e.g.
In [136]: eqs = [Eq(x(t).diff(t), y(t)), Eq(y(t).diff(t), -x(t))]
In [137]: eqs
Out[137]:
⎡d d ⎤
⎢──(x(t)) = y(t), ──(y(t)) = -x(t)⎥
⎣dt dt ⎦
In [138]: dsolve(eqs)
Out[138]: [x(t) = C₁⋅cos(t) + C₂⋅sin(t), y(t) = -C₁⋅sin(t) + C₂⋅cos(t)]Computing the matrix exponential directly from the Jordan normal form gives
In [147]: (M*t).exp()
Out[147]:
⎡ ⅈ⋅t -ⅈ⋅t ⅈ⋅t -ⅈ⋅t⎤
⎢ ℯ ℯ ⅈ⋅ℯ ⅈ⋅ℯ ⎥
⎢ ──── + ───── - ────── + ───────⎥
⎢ 2 2 2 2 ⎥
⎢ ⎥
⎢ ⅈ⋅t -ⅈ⋅t ⅈ⋅t -ⅈ⋅t ⎥
⎢ⅈ⋅ℯ ⅈ⋅ℯ ℯ ℯ ⎥
⎢────── - ─────── ──── + ───── ⎥
⎣ 2 2 2 2 ⎦
In [148]: P * (J*t).exp() * P.inv()
Out[148]:
⎡ ⅈ⋅t -ⅈ⋅t ⅈ⋅t -ⅈ⋅t⎤
⎢ ℯ ℯ ⅈ⋅ℯ ⅈ⋅ℯ ⎥
⎢ ──── + ───── - ────── + ───────⎥
⎢ 2 2 2 2 ⎥
⎢ ⎥
⎢ ⅈ⋅t -ⅈ⋅t ⅈ⋅t -ⅈ⋅t ⎥
⎢ⅈ⋅ℯ ⅈ⋅ℯ ℯ ℯ ⎥
⎢────── - ─────── ──── + ───── ⎥
⎣ 2 2 2 2 ⎦There are ways to take care of this and produce real-valued results when desired e.g. https://en.wikipedia.org/wiki/Jordan_normal_form#Real_matrices which is implemented in https://github.com/sympy/sympy/pull/15449
So the idea here is that we want to find the general solution of X' = A X and we know that is given by X = exp(A t) C where C is a vector of arbitrary constants. In order to compute exp(A t) we compute it as P * (J*t).exp() * P.inv() where P, J is the Jordan decomposition of A. That means our solution is
X = P * (J*t).exp() * P.inv() * C
Remembering that C is just a vector of arbitrary constants we can define a new vector D := P.inv() * C which is also a vector of arbitrary constants which means that our general solution becomes
X = P * (J*t).exp() * D
This amounts to a redefinition of the arbitrary constants but is still equally correct as a solution of the system. The advantage though is that we no longer need to compute the inverse of P which takes time and also can introduce significant algebraic complexity into the solution.
Let's try this out for a particular system of ODEs. We take
In [156]: eqs = [Eq(x(t).diff(t), x(t) + y(t)), Eq(y(t).diff(t), y(t))]
In [157]: eqs
Out[157]:
⎡d d ⎤
⎢──(x(t)) = x(t) + y(t), ──(y(t)) = y(t)⎥
⎣dt dt ⎦
In [158]: M = Matrix([[1, 1], [0, 1]])
In [159]: M
Out[159]:
⎡1 1⎤
⎢ ⎥
⎣0 1⎦
In [160]: P, J = M.jordan_form()
In [161]: C = Matrix([[C1], [C2]])
In [162]: P * (J*t).exp() * P.inv() * C
Out[162]:
⎡ t t⎤
⎢C₁⋅ℯ + C₂⋅t⋅ℯ ⎥
⎢ ⎥
⎢ t ⎥
⎣ C₂⋅ℯ ⎦
In [163]: P * (J*t).exp() * C
Out[163]:
⎡ t t⎤
⎢C₁⋅ℯ + C₂⋅t⋅ℯ ⎥
⎢ ⎥
⎢ t ⎥
⎣ C₂⋅ℯ ⎦
In [164]: dsolve(eqs)
Out[164]:
⎡ t t⎤
⎣x(t) = (C₁ + C₂⋅t)⋅ℯ , y(t) = C₂⋅ℯ ⎦In this example using P.inv() or not makes no difference as it is the identity matrix but where there are symbolic coefficients it can make a significant difference.
The current situation is that SymPy's ODE systems code can solve a number of classes of constant coefficient ODEs using special case solvers e.g.:
https://docs.sympy.org/dev/modules/solvers/ode.html#linear-2-equations-order-1-type-1 https://docs.sympy.org/dev/modules/solvers/ode.html#linear-3-equations-order-1-type-1 https://docs.sympy.org/dev/modules/solvers/ode.html#linear-3-equations-order-1-type-2 https://docs.sympy.org/dev/modules/solvers/ode.html#linear-3-equations-order-1-type-3 https://docs.sympy.org/dev/modules/solvers/ode.html#linear-n-equations-order-1-type-1
The most commonly used of these is the 2-equation solver which has problems with equations involving symbols e.g.
eqs = [Eq(x(t).diff(t), z*x(t)), Eq(y(t).diff(t), y(t))]
(dsolve(eqs) produces a complicated Piecewise result.)
The n-eq solver in this class is the only n-equations solver that dsolve has and is broken for defective systems: https://github.com/sympy/sympy/issues/14312
There are two stalled PRs attempting to fix this: https://github.com/sympy/sympy/pull/15449 https://github.com/sympy/sympy/pull/14445
If a single n-equations solver was made to work for this then the other solvers could all be removed and sympy would have its first actually working n-equations solver for ODE systems.
I suggest the following:
- Implement a fully working
n-equations solver for constant coefficient homogeneous systems. - Implement this as a class similar to the refactored classes for single ODEs.
- Add the matching code for this right at the start of
classify_sysodeso that it takes priority over all other solvers. - At this point the special case solvers mentioned above will no longer be needed and could be deleted.
The above deals with ODE systems of the form
X' = A X
where A is a matrix of constant coefficients. Such systems are called homogeneous. The corresponding class of non-homogeneous systems is written as
X' = A X + f(t)
where the rhs is a vector that depends only on the independent variable t. For this we can use a generalisation of the method of integrating factor for single linear 1st order ODEs first rearrange
X' - A X = f(t)
Now we multiply through by exp(-A t) to get
exp(-A t) X' - A exp(-A t) X = exp(-A t) * f(t)
(Note that A commutes with exp(-A t).) We can now see that the lhs can be written as a derivative of a product
d/dt{ exp(-A t) X } = exp(-A t) X' - A exp(-A t) X
so the ODE becomes
d/dt{ exp(-A t) X } = exp(- A t) f(t)
Now the rhs depends only on t so we can integrate both sides
exp(-A t) X = int_x { exp(- A t) f(t) } + C
where C is a vector of constants. Finally we multiply through by exp(A t) to find
X = exp(A t) int_x { exp(- A t) f(t) } + exp(A t) C
We can see that this is consistent with the homogeneous case since there f(t) is zero and we would just have X = exp(A t) C as before.
Let's work through an example. Here is a system that dsolve can not currently solve:
In [286]: eqs = [Eq(x(t).diff(t), y(t)), Eq(y(t).diff(t), x(t) + t)]
In [287]: eqs
Out[287]:
⎡d d ⎤
⎢──(x(t)) = y(t), ──(y(t)) = t + x(t)⎥
⎣dt dt ⎦
In [288]: dsolve(eqs)
---------------------------------------------------------------------------
NotImplementedError We can construct the matrix and build the general solution:
In [297]: A = Matrix([[0, 1], [1, 0]])
In [298]: f = Matrix([[0], [t]])
In [299]: X = Matrix([[x(t)], [y(t)]])
In [300]: Eq(Derivative(X, t), MatAdd(MatMul(A, X, evaluate=False), f, evaluate=False), evaluate=False)
Out[300]:
d ⎛⎡x(t)⎤⎞ ⎡0 1⎤ ⎡x(t)⎤ ⎡0⎤
──⎜⎢ ⎥⎟ = ⎢ ⎥⋅⎢ ⎥ + ⎢ ⎥
dt⎝⎣y(t)⎦⎠ ⎣1 0⎦ ⎣y(t)⎦ ⎣t⎦Introducing the constants we can build the general solution and verify that it is a solution
In [22]: C1, C2 = symbols('C1, C2')
In [23]: C = Matrix([[C1], [C2]])
In [24]: gensol = (A*t).exp() * (integrate((-A*t).exp() * f) + C)
In [25]: gensol
Out[25]:
⎡⎛ t -t⎞ ⎛ -t t⎞ ⎛ t -t⎞ ⎛ t -t⎞⎤
⎢⎜ℯ ℯ ⎟ ⎜ (-2⋅t - 2)⋅ℯ (2⋅t - 2)⋅ℯ ⎟ ⎜ℯ ℯ ⎟ ⎜ (2 - 2⋅t)⋅ℯ (-2⋅t - 2)⋅ℯ ⎟⎥
⎢⎜── - ───⎟⋅⎜C₂ + ────────────── + ────────────⎟ + ⎜── + ───⎟⋅⎜C₁ + ──────────── + ──────────────⎟⎥
⎢⎝2 2 ⎠ ⎝ 4 4 ⎠ ⎝2 2 ⎠ ⎝ 4 4 ⎠⎥
⎢ ⎥
⎢⎛ t -t⎞ ⎛ t -t⎞ ⎛ t -t⎞ ⎛ -t t⎞⎥
⎢⎜ℯ ℯ ⎟ ⎜ (2 - 2⋅t)⋅ℯ (-2⋅t - 2)⋅ℯ ⎟ ⎜ℯ ℯ ⎟ ⎜ (-2⋅t - 2)⋅ℯ (2⋅t - 2)⋅ℯ ⎟⎥
⎢⎜── - ───⎟⋅⎜C₁ + ──────────── + ──────────────⎟ + ⎜── + ───⎟⋅⎜C₂ + ────────────── + ────────────⎟⎥
⎣⎝2 2 ⎠ ⎝ 4 4 ⎠ ⎝2 2 ⎠ ⎝ 4 4 ⎠⎦
In [26]: gensol = simplify(gensol.rewrite(sin))
In [27]: gensol
Out[27]:
⎡C₁⋅cosh(t) + C₂⋅sinh(t) - t⎤
⎢ ⎥
⎣C₁⋅sinh(t) + C₂⋅cosh(t) - 1⎦
In [28]: sol = [Eq(x(t), gensol[0]), Eq(y(t), gensol[1])]
In [29]: sol
Out[29]: [x(t) = C₁⋅cosh(t) + C₂⋅sinh(t) - t, y(t) = C₁⋅sinh(t) + C₂⋅cosh(t) - 1]
In [30]: from sympy.solvers.ode.subscheck import checksysodesol
In [31]: checksysodesol(eqs, sol)
Out[31]: (True, [0, 0])This means that we can express the solution of any first-order constant coefficient non-homogeneous ODE system using integrals. It might not be possible to evaluate the resulting integrals but we can at least express the solution in terms of integrals. Combined with the above results for homogeneous ODEs this means we can solve all constant coefficient first order systems both homogeneous and non-homogeneous and for any number of equations.
Above we see a solver that can express the solution for any first-order system of non-homogeneous ODEs with constant coefficients. This is a substantial generalisation of the current capability for handling non-homogeneous systems which can only handle a system of 2 first order ODEs where the function f(t) is constant:
https://docs.sympy.org/dev/modules/solvers/ode.html#linear-2-equations-order-1-type-2
That solver could be removed once a general n-equations solver using the matrix exponential is implemented.
For systems of higher order ODEs dsolve is currently limited to a few cases where there are two 2nd order ODEs. In particular what dsolve can not do is:
- Handle any system of 3 or more 2nd order ODEs.
- Handle any system involving ODEs of 3rd or higher order.
- Handle any system involving ODEs of mixed order.
However all of these cases can be reduced to 1st order and we have seen above that we can solve all 1st order ODEs with constant coefficients both homogeneous and non-homogeneous for any number of equations. Simply being able to recast higher-order systems in first order would enormously expand the class of ODE systems that can be solved.
An example:
In [40]: x, y, z = symbols('x, y, z', cls=Function)
In [41]: t = Symbol('t')
In [42]: eqs = [Eq(x(t).diff(t), y(t)), Eq(y(t).diff(t, 2), x(t))]
In [43]: eqs
Out[43]:
⎡ 2 ⎤
⎢d d ⎥
⎢──(x(t)) = y(t), ───(y(t)) = x(t)⎥
⎢dt 2 ⎥
⎣ dt ⎦
In [44]: dsolve(eqs)
---------------------------------------------------------------------------
ValueErrorIntroducing a new variable z = y' we can recast this as a system of 1st order ODEs that dsolve can solve:
In [52]: eqs2 = [Eq(x(t).diff(t), y(t)), Eq(y(t).diff(t), z(t)), Eq(z(t).diff(t), x(t))]
In [53]: dsolve(eqs2)
Out[53]:
⎡ ⎛ ⎛√3⋅t⎞ ⎛√3⋅t⎞⎞ -t ⎛ ⎛√3⋅t⎞ ⎛√3⋅t⎞⎞ -t ⎛ ⎛√3⋅t⎞ ⎛√3⋅t⎞⎞ -t
⎢ ⎜ √3⋅sin⎜────⎟ cos⎜────⎟⎟ ─── ⎜ sin⎜────⎟ √3⋅cos⎜────⎟⎟ ─── ⎜√3⋅sin⎜────⎟ cos⎜────⎟⎟ ──
⎢ t ⎜ ⎝ 2 ⎠ ⎝ 2 ⎠⎟ 2 ⎜ ⎝ 2 ⎠ ⎝ 2 ⎠⎟ 2 t ⎜ ⎝ 2 ⎠ ⎝ 2 ⎠⎟ 2
⎢x(t) = C₁⋅ℯ + C₂⋅⎜- ──────────── - ─────────⎟⋅ℯ + C₃⋅⎜- ───────── + ────────────⎟⋅ℯ , y(t) = C₁⋅ℯ + C₂⋅⎜──────────── - ─────────⎟⋅ℯ
⎣ ⎝ 2 2 ⎠ ⎝ 2 2 ⎠ ⎝ 2 2 ⎠
⎛ ⎛√3⋅t⎞ ⎛√3⋅t⎞⎞ -t -t -t ⎤
─ ⎜ sin⎜────⎟ √3⋅cos⎜────⎟⎟ ─── ─── ─── ⎥
⎜ ⎝ 2 ⎠ ⎝ 2 ⎠⎟ 2 t 2 ⎛√3⋅t⎞ 2 ⎛√3⋅t⎞⎥
+ C₃⋅⎜- ───────── - ────────────⎟⋅ℯ , z(t) = C₁⋅ℯ + C₂⋅ℯ ⋅cos⎜────⎟ + C₃⋅ℯ ⋅sin⎜────⎟⎥
⎝ 2 2 ⎠ ⎝ 2 ⎠ ⎝ 2 ⎠⎦
In [54]: xs, ys, zs = dsolve(eqs2)
In [55]: checksysodesol(eqs, [xs, ys])
Out[55]: (True, [0, 0])Implementing a solver that recasts constant coefficient systems in 1st order would generalise some of the existing special case solvers that could potentially be removed after e.g.:
https://docs.sympy.org/dev/modules/solvers/ode.html#linear-2-equations-order-2-type-1 https://docs.sympy.org/dev/modules/solvers/ode.html#linear-2-equations-order-2-type-2 https://docs.sympy.org/dev/modules/solvers/ode.html#linear-2-equations-order-2-type-3 https://docs.sympy.org/dev/modules/solvers/ode.html#linear-2-equations-order-2-type-4
Some of those solvers might still be useful as it may be more efficient to get the solution for e.g. a 2nd order matrix ODE
X'' = A X (2)
without recasting to 1st order form. Certainly on pen and paper (2) is easier to solve in 2nd order form because if A is a 2x2 matrix then converting to 1st order makes it 4x4 which increases the amount of algebra needed. For a CAS though it might still be better to handle the 4x4 matrix than have special code for this case.
If dsolve for systems could recast systems in 1st order form then it would be possible (using also the results above) to solve any system of constant coefficient ODEs, both homogeneous and non-homogeneous of any number of equations of any order including of mixed order. That would be a substantial improvement over the current situation.
The above describes a complete solution for the case of systems of ODEs with constant coefficients. Constant coefficients is an important special case for which a complete solution is possible and it is important for sympy to have a complete implementation of solvers for that case. There are however important cases involving non-constant coefficients that can be solved as well.
Here is an example of a system that does not have constant coefficients:
In [56]: eqs = [Eq(x(t).diff(t), t*y(t)), Eq(y(t).diff(t), t*x(t))]
In [57]: eqs
Out[57]:
⎡d d ⎤
⎢──(x(t)) = t⋅y(t), ──(y(t)) = t⋅x(t)⎥
⎣dt dt ⎦
In [58]: dsolve(eqs)
Out[58]:
⎡ ⌠ ⌠ ⌠ ⌠ ⎤
⎢ ⎮ t dt -⎮ t dt ⎮ t dt -⎮ t dt⎥
⎢ ⌡ ⌡ ⌡ ⌡ ⎥
⎣x(t) = C₁⋅ℯ + C₂⋅ℯ , y(t) = C₁⋅ℯ - C₂⋅ℯ ⎦Since the coefficients of x(t) and y(t) on the rhs of the ODEs includes t this system does not have constant coefficients. The coefficient matrix is
In [63]: M = Matrix([[0, t], [t, 0]])
In [64]: M
Out[64]:
⎡0 t⎤
⎢ ⎥
⎣t 0⎦
In [65]: M.is_symmetric()
Out[65]: Trueand is symmetric. This is another important class that can be solved. We have a system of ODEs represented as a matrix ODE as
X' = A(t) X (3)
where A(t) is symmetric is a function of the independent variable t but not any of the dependent variables.
An important property of the matrix exponential is that if A and B commute then A and exp(B) commute. This is easily verified from the power series definition of the matrix exponential. What this means is that if A(t) is symmetric and has symmetric antiderivative B(t) (i.e. B' = A) then A(t) commutes with exp(B(t)) and we have
d/dt{ exp(B(t)) } = d/dt{ I + B + B^2/2 + B^3/6 + ...
= 0 + A + AB/2 + BA/2 + AB^2/6 + BAB/6 + B^2A/6 ...
= A + 2AB/2 + 3AB^2/6 + ...
= A + AB + AB^2/2 + ...
= A(I + B + B^2/2 + ...)
= A exp(B)
The key result here is that if B is an antiderivative of A and A and B commute then
d/dt{ exp(B(t)) } = A(t) exp(B(t))
The condition that A and B commute is easily satisfied if A is symmetric but is in fact more general than that. If we can find a B satisfying this then the system (3) has solution
X = exp(B(t)) C
where C is a vector of constants.
Let's use this to solve the example above:
In [73]: C1, C2 = symbols('C1, C2')
In [74]: C = Matrix([[C1], [C2]])
In [75]: gensol = B.exp() * C
In [76]: gensol
Out[76]:
⎡ ⎛ 2 2 ⎞ ⎛ 2 2 ⎞⎤
⎢ ⎜ t -t ⎟ ⎜ t -t ⎟⎥
⎢ ⎜ ── ────⎟ ⎜ ── ────⎟⎥
⎢ ⎜ 2 2 ⎟ ⎜ 2 2 ⎟⎥
⎢ ⎜ℯ ℯ ⎟ ⎜ℯ ℯ ⎟⎥
⎢C₁⋅⎜─── + ─────⎟ + C₂⋅⎜─── - ─────⎟⎥
⎢ ⎝ 2 2 ⎠ ⎝ 2 2 ⎠⎥
⎢ ⎥
⎢ ⎛ 2 2 ⎞ ⎛ 2 2 ⎞⎥
⎢ ⎜ t -t ⎟ ⎜ t -t ⎟⎥
⎢ ⎜ ── ────⎟ ⎜ ── ────⎟⎥
⎢ ⎜ 2 2 ⎟ ⎜ 2 2 ⎟⎥
⎢ ⎜ℯ ℯ ⎟ ⎜ℯ ℯ ⎟⎥
⎢C₁⋅⎜─── - ─────⎟ + C₂⋅⎜─── + ─────⎟⎥
⎣ ⎝ 2 2 ⎠ ⎝ 2 2 ⎠⎦
In [77]: gensol = simplify(gensol.rewrite(sin))
In [78]: gensol
Out[78]:
⎡ ⎛ 2⎞ ⎛ 2⎞⎤
⎢ ⎜t ⎟ ⎜t ⎟⎥
⎢C₁⋅cosh⎜──⎟ + C₂⋅sinh⎜──⎟⎥
⎢ ⎝2 ⎠ ⎝2 ⎠⎥
⎢ ⎥
⎢ ⎛ 2⎞ ⎛ 2⎞⎥
⎢ ⎜t ⎟ ⎜t ⎟⎥
⎢C₁⋅sinh⎜──⎟ + C₂⋅cosh⎜──⎟⎥
⎣ ⎝2 ⎠ ⎝2 ⎠⎦
In [79]: sol = [Eq(x(t), gensol[0]), Eq(y(t), gensol[1])]
In [80]: sol
Out[80]:
⎡ ⎛ 2⎞ ⎛ 2⎞ ⎛ 2⎞ ⎛ 2⎞⎤
⎢ ⎜t ⎟ ⎜t ⎟ ⎜t ⎟ ⎜t ⎟⎥
⎢x(t) = C₁⋅cosh⎜──⎟ + C₂⋅sinh⎜──⎟, y(t) = C₁⋅sinh⎜──⎟ + C₂⋅cosh⎜──⎟⎥
⎣ ⎝2 ⎠ ⎝2 ⎠ ⎝2 ⎠ ⎝2 ⎠⎦We can verify this solution:
In [89]: xs, ys = gensol
In [90]: eqs[0].subs({x(t): xs, y(t): ys}).doit().expand()
Out[90]: True
In [91]: eqs[1].subs({x(t): xs, y(t): ys}).doit().expand()
Out[91]: True
For some reason though checksysodesol fails here:
In [81]: checksysodesol(eqs, [xs, ys])
Out[81]:
(False,
[-C1*t*exp(t) + C1*exp(t) + C2*t*exp(-t/2)*cos(sqrt(3)*t/2 + pi/3) - C2*exp(-t/2)*cos(sqrt(3)*t/2 + pi/3) + C3*t*exp(-t/2)*sin(sqrt(3)*t/2 + pi/3) - C3*exp(-t/2)*sin(sqrt(3)*t/2 + pi/3),
-C1*t*exp(t) + C1*exp(t) + C2*t*exp(-t/2)*sin(sqrt(3)*t/2 + pi/6) + C2*exp(-t/2)*cos(sqrt(3)*t/2) - C3*t*exp(-t/2)*cos(sqrt(3)*t/2 + pi/6) + C3*exp(-t/2)*sin(sqrt(3)*t/2)])
(That's probably a bug in checksysodesol).
The above provides a method for solving equations of the form of (3) where A is symmetric or more generally commutes with some antiderivative. This provides a method that can work for any number of equations. Since higher-order systems can be recast as 1st order the method can also work for those as well.
There are a number of special case solvers that can be generalised using this method (I haven't fully checked all of these but I think they can):
https://docs.sympy.org/dev/modules/solvers/ode.html#linear-2-equations-order-1-type-3 https://docs.sympy.org/dev/modules/solvers/ode.html#linear-2-equations-order-1-type-4 https://docs.sympy.org/dev/modules/solvers/ode.html#linear-2-equations-order-1-type-5 https://docs.sympy.org/dev/modules/solvers/ode.html#linear-2-equations-order-2-type-5 https://docs.sympy.org/dev/modules/solvers/ode.html#linear-2-equations-order-2-type-6 https://docs.sympy.org/dev/modules/solvers/ode.html#linear-2-equations-order-2-type-7 https://docs.sympy.org/dev/modules/solvers/ode.html#linear-2-equations-order-2-type-8 https://docs.sympy.org/dev/modules/solvers/ode.html#linear-2-equations-order-2-type-9 https://docs.sympy.org/dev/modules/solvers/ode.html#linear-2-equations-order-2-type-10 https://docs.sympy.org/dev/modules/solvers/ode.html#linear-3-equations-order-1-type-4
If an n-equation solver that handles this case was implemented then these special case solvers could be removed.
In combination with the above this would mean that there are 3 n-equation solvers for these 3 types of system:
X' = A X
X' = A X + f(t)
X' = A(t) X (where A commutes with antiderivative)
These would encompass 20 out of the 23 solvers for linear systems of ODEs in 3 n-equation solvers. This would significantly simplify the code as well as greatly expanding the capability of dsolve for systems.
The 3 remaining linear solvers that are not included are:
https://docs.sympy.org/dev/modules/solvers/ode.html#linear-2-equations-order-1-type-7 https://docs.sympy.org/dev/modules/solvers/ode.html#linear-2-equations-order-1-type-6 https://docs.sympy.org/dev/modules/solvers/ode.html#linear-2-equations-order-2-type-11
The _linear_2eq_order1_type7 solver does not make any sense to me. I don't understand what the docstring means as an explanation for the method. The only test for it returns an answer in terms of so called "particular solutions":
In [1]: x, y = symbols('x, y', cls=Function)
In [2]: eq10 = (Eq(diff(x(t),t), 5*t*x(t) + t**2*y(t)), Eq(diff(y(t),t), (1-t**2)*x(t) + (5*t+9*t**2)*y(t)))
In [3]: eq10
Out[3]:
⎛d 2 d ⎛ 2⎞ ⎛ 2 ⎞ ⎞
⎜──(x(t)) = t ⋅y(t) + 5⋅t⋅x(t), ──(y(t)) = ⎝1 - t ⎠⋅x(t) + ⎝9⋅t + 5⋅t⎠⋅y(t)⎟
⎝dt dt ⎠
In [5]: repr(dsolve(eq10))
Out[5]: '[Eq(x(t), C1*x0(t) + C2*x0(t)*Integral(t**2*exp(Integral(5*t, t))*exp(Integral(9*t**2 + 5*t, t))/x0(t)**2, t)), Eq(y(t), C1*y0(t) + C2*(y0(t)*Integral(t**2*exp(Integral(5*t, t))*exp(Integral(9*t**2 + 5*t, t))/x0(t)**2, t) + exp(Integral(5*t, t))*exp(Integral(9*t**2 + 5*t, t))/x0(t)))]'
so the method does not appear to solve any systems properly. Unless this can be found to solve anything I suggest to remove this method.
The _linear_2eq_order1_type6 method is only tested by XFAIL tests showing that it returns incorrect results. The method can probably be generalised. Essentially it amounts to this: Given a linear non-constant coefficient homogeneous system for x and y try to find a linear combination z=ax+by so that we get an equation z'=f(z) that eliminates both x and y on the rhs and can be solved as a single ODE. If that works then the solution for z gives a linear relationship between x and y that can be used to eliminate one or the other from the original ODEs to obtain either x' = g(x) or y' = h(y) which can then be solved as a single ODE. A generalised solver along these lines could be useful but the current solver is not generalised and does not work. I suggest to remove this method.
The _linear_2eq_order2_type11 method only has a single test:
In [5]: x, y = symbols('x, y', cls=Function)
In [6]: eq7 = (Eq(diff(x(t),t,t), log(t)*(t*diff(x(t),t) - x(t)) + exp(t)*(t*diff(y(t),t) - y(t))), \
...: Eq(diff(y(t),t,t), (t**2)*(t*diff(x(t),t) - x(t)) + (t)*(t*diff(y(t),t) - y(t))))
In [7]: repr(dsolve(eq7))
Out[7]: '[Eq(x(t), C3*t + t*Integral((C1*x0(t) + C2*x0(t)*Integral(t*exp(t)*exp(Integral(t**2, t))*exp(Integral(t*log(t), t))/x0(t)**2, t))/t**2, t)), Eq(y(t), C4*t + t*Integral((C1*y0(t) + C2*(y0(t)*Integral(t*exp(t)*exp(Integral(t**2, t))*exp(Integral(t*log(t), t))/x0(t)**2, t) + exp(Integral(t**2, t))*exp(Integral(t*log(t), t))/x0(t)))/t**2, t))]'
This works by reducing a particular kind of 2nd order to system to a 1st order system in the form supposedly solved by _linear_2eq_order1_type7. The solution returned comes out in terms of "particular solutions" because of the failure of the type7 solver. This method seems obscure and I am not convinced that it ever works properly. I suggest to remove this method.
If we remove these 3 questionable methods then it seems that all of the currently working linear solvers are special cases of the three general case solvers mentioned above. At that point it is clear what classes of linear ODE systems can be solved by dsolve and the documentation can be updated so that it is clear to users as well.
The natural next step is to consider
X' = A(t) X + f(t)
which is indeed possible using the same techniques as above for the constant coefficient case. Given a commuting antiderivative B (B' = A and A and B commute) we have:
X' - A(t) X = f(t)
exp(-B(t)) X' - A(t) exp(-B(t)) X = exp(-B(t)) f(t)
d/dt{ exp(-B(t)) X } = exp(-B(t)) f(t)
exp(-B(t)) X = int_x { exp(-B(t)) f(t) } + C
X = exp(B(t)) int_x { exp(-B(t)) f(t) } + exp(B(t)) C
where C is a vector of constants.
Let's follow through an example:
In [20]: x, y = symbols('x, y', cls=Function)
In [21]: t = Symbol('t')
In [22]: eqs = [Eq(x(t).diff(t), t*y(t) + t), Eq(y(t).diff(t), t*x(t))]
In [23]: eqs
Out[23]:
⎡d d ⎤
⎢──(x(t)) = t⋅y(t) + t, ──(y(t)) = t⋅x(t)⎥
⎣dt dt ⎦
In [24]: dsolve(eqs)
---------------------------------------------------------------------------
NotImplementedErrorSo dsolve can not currently solve this. However using the method above we can:
In [40]: A = Matrix([[0, t], [t, 0]])
In [41]: f = Matrix([[t], [0]])
In [42]: B = integrate(A, t)
In [43]: C1, C2 = symbols('C1, C2')
In [44]: C = Matrix([[C1], [C2]])
In [45]: gensol = B.exp() * integrate((-B).exp() * f, t) + B.exp() * C
In [46]: gensol = simplify(gensol.rewrite(sin))
In [47]: sol = [Eq(x(t), gensol[0]), Eq(y(t), gensol[1])]
In [48]: sol
Out[48]:
⎡ ⎛ 2⎞ ⎛ 2⎞ ⎛ 2⎞ ⎛ 2⎞ ⎤
⎢ ⎜t ⎟ ⎜t ⎟ ⎜t ⎟ ⎜t ⎟ ⎥
⎢x(t) = C₁⋅cosh⎜──⎟ + C₂⋅sinh⎜──⎟, y(t) = C₁⋅sinh⎜──⎟ + C₂⋅cosh⎜──⎟ - 1⎥
⎣ ⎝2 ⎠ ⎝2 ⎠ ⎝2 ⎠ ⎝2 ⎠ ⎦
In [49]: from sympy.solvers.ode.subscheck import checksysodesol
In [50]: checksysodesol(eqs, sol)
Out[50]: (True, [0, 0])Currently sympy can not solve any non-homogeneous systems with non-constant coefficients so this is a new n-equations solver for capability that sympy does not have at all right now.
The current dsolve for systems code makes no attempt to rearrange the ODEs to see if they can be put in a form that it will recognise. All solvers match ODEs in the form x' = ..., y' = ... so even this will confuse dsolve:
In [52]: eqs = [Eq(x(t).diff(t) + y(t).diff(t), x(t)), Eq(x(t).diff(t) - y(t).diff(t), y(t))]
In [53]: eqs
Out[53]:
⎡d d d d ⎤
⎢──(x(t)) + ──(y(t)) = x(t), ──(x(t)) - ──(y(t)) = y(t)⎥
⎣dt dt dt dt ⎦
In [54]: dsolve(eqs)
---------------------------------------------------------------------------
AttributeError It is trivial though to add and subtract these equations to get two equations that are in standard form. I think that dsolve should at minimum attempt linear rearrangements to get the system into canonical form. The first step in dsolve should be to solve for the derivatives algebraically:
In [55]: solve(eqs, [x(t).diff(t), y(t).diff(t)])
Out[55]:
⎧d x(t) y(t) d x(t) y(t)⎫
⎨──(x(t)): ──── + ────, ──(y(t)): ──── - ────⎬
⎩dt 2 2 dt 2 2 ⎭
In [56]: dsol = solve(eqs, [x(t).diff(t), y(t).diff(t)])
In [57]: eqs2 = [Eq(x(t).diff(t), dsol[x(t).diff(t)]), Eq(y(t).diff(t), dsol[y(t).diff(t)])]
In [58]: eqs2
Out[58]:
⎡d x(t) y(t) d x(t) y(t)⎤
⎢──(x(t)) = ──── + ────, ──(y(t)) = ──── - ────⎥
⎣dt 2 2 dt 2 2 ⎦
In [59]: dsolve(eqs2)
Out[59]:
⎡ √2⋅t -√2⋅t ⎤
⎢ ──── ────── √2⋅t -√2⋅t ⎥
⎢ 2 2 ──── ──────⎥
⎢ C₁⋅ℯ C₂⋅ℯ ⎛ 1 √2⎞ 2 ⎛ √2 1⎞ 2 ⎥
⎢x(t) = ──────── + ──────────, y(t) = C₁⋅⎜- ─ + ──⎟⋅ℯ + C₂⋅⎜- ── - ─⎟⋅ℯ ⎥
⎣ 2 2 ⎝ 2 2 ⎠ ⎝ 2 2⎠ ⎦Some systems of ODEs can be partitioned into independent systems e.g.:
In [14]: x, y = symbols('x, y', cls=Function)
In [15]: t = Symbol('t')
In [16]: eqs = [Eq(x(t).diff(t, 2), x(t)), Eq(y(t).diff(t), y(t)**2)]
In [17]: eqs
Out[17]:
⎡ 2 ⎤
⎢ d d 2 ⎥
⎢───(x(t)) = x(t), ──(y(t)) = y (t)⎥
⎢ 2 dt ⎥
⎣dt ⎦
In [18]: dsolve(eqs)
---------------------------------------------------------------------------
ValueErrorThis system is trivially solvable since we can solve it as two independent ODEs:
In [19]: sol1 = dsolve(eqs[0])
In [20]: sol2 = dsolve(eqs[1])
In [21]: sol1
Out[21]:
-t t
x(t) = C₁⋅ℯ + C₂⋅ℯ
In [22]: sol2
Out[22]:
-1
y(t) = ──────
C₁ + tPutting these together for the general solution of the system requires renumbering the constants. The solution from dsolve should be
In [23]: C1, C2, C3 = symbols('C1, C2, C3')
In [24]: sol = [sol1, sol2.subs(C1, C3)]
In [25]: sol
Out[25]:
⎡ -t t -1 ⎤
⎢x(t) = C₁⋅ℯ + C₂⋅ℯ , y(t) = ──────⎥
⎣ C₃ + t⎦
In [26]: from sympy.solvers.ode.subscheck import checksysodesol
In [27]: checksysodesol(eqs, sol)
Out[27]: (True, [0, 0])In this example it is easy to see that the equations can be separated since we have one ODE that only involves x and another that only involves y. More generally though we might have something like:
In [38]: x, y, z, w = symbols('x, y, z, w', cls=Function)
In [39]: eqs = [Eq(x(t).diff(t), y(t)), Eq(y(t).diff(t), x(t)), Eq(z(t).diff(t), w(t)), Eq(w(t).diff(t), z(t))]
In [40]: eqs
Out[40]:
⎡d d d d ⎤
⎢──(x(t)) = y(t), ──(y(t)) = x(t), ──(z(t)) = w(t), ──(w(t)) = z(t)⎥
⎣dt dt dt dt ⎦In this example we can solve the first two ODEs for x and y separately from solving the other two for z and w. The question is how to spot that the system can be partitioned in this way.
We define the idea of a directed graph between the dependent variables. We assume that the system is in canonical form x' = f(x, y, ...), y' = g(x, y, ...), .... Now if variable y appears in the rhs of the equation for x' we say that there is an edge from x to y. In this way we can represent the system above as (using dot notation):
digraph G{
x -> y
y -> x
z -> w
w -> z
}
Having done this we can partition the graph into its (weakly) connected components. There is already a function connected_components in sympy for partitioning the graph. We just need to write the code that constructs the graph in the right form:
def ode_graph(eqs):
funcs = [eq.lhs.args[0] for eq in eqs]
rhss = [eq.rhs for eq in eqs]
edges = []
for func, rhs in zip(funcs, rhss):
for other in funcs:
if other != func and rhs.has(other):
edges.append((func, other))
return funcs, edgesWith that we can easily partition the system into connected components:
In [56]: nodes, edges = ode_graph(eqs)
In [57]: connected_components((nodes, edges))
Out[57]: [[x(t), y(t)], [z(t), w(t)]]Partitioning the system in this way is always beneficial. Even if the original system (like this example) is solvable as a whole system it is still much faster to solve the smaller systems having smaller matrices than to solve the whole system with its larger matrix. This should be the first thing that dsolve does with any ODE system before attempting any kind of solving or matching against known ODE patterns.
Having split the system into connected components we know that we are working with equations where each variable appears in an ODE with each of the others. Further separation is still possible though. Consider this:
In [58]: eqs = [Eq(x(t).diff(t), x(t)**2), Eq(y(t).diff(t), x(t))]
In [59]: eqs
Out[59]:
⎡d 2 d ⎤
⎢──(x(t)) = x (t), ──(y(t)) = x(t)⎥
⎣dt dt ⎦
In [60]: dsolve(eqs)
---------------------------------------------------------------------------
TypeErrorThis system can be solved by first solving for x substituting that into the ODE for y and then solving for y. In this case x and y are in the same connected component:
In [61]: nodes, edges = ode_graph(eqs)
In [62]: connected_components((nodes, edges))
Out[62]: [[x(t), y(t)]]This is because the ODE for y involves x so they are "connected". We can identify this sort of case by partitioning the connected components further into "strongly connected" components and sympy has a function for this as well:
In [63]: from sympy.utilities.iterables import strongly_connected_components
In [64]: strongly_connected_components((nodes, edges))
Out[64]: [[x(t)], [y(t)]]The strongly_connected_components function returns results in reverse topological order which represents the order in which we should attempt to solve the subsystems:
In [65]: solx = dsolve(eqs[0], x(t))
In [66]: solx
Out[66]:
-1
x(t) = ──────
C₁ + t
In [67]: eqy = eqs[1].subs(x(t), solx.rhs)
In [68]: eqy
Out[68]:
d -1
──(y(t)) = ──────
dt C₁ + t
In [69]: soly = dsolve(eqy)
In [70]: soly
Out[70]: y(t) = C₂ - log(C₁ + t)
In [71]: sol = [solx, soly]
In [72]: checksysodesol(eqs, sol)
Out[72]: (True, [0, 0])
In [73]: sol
Out[73]:
⎡ -1 ⎤
⎢x(t) = ──────, y(t) = C₂ - log(C₁ + t)⎥
⎣ C₁ + t ⎦We can use this in more complex examples:
In [74]: eqs = [Eq(x(t).diff(t), x(t)**2), Eq(y(t).diff(t), x(t)+z(t)), Eq(z(t).diff(t), y(t)), Eq(w(t).diff(t), z(t))]
In [75]: eqs
Out[75]:
⎡d 2 d d d ⎤
⎢──(x(t)) = x (t), ──(y(t)) = x(t) + z(t), ──(z(t)) = y(t), ──(w(t)) = z(t)⎥
⎣dt dt dt dt ⎦
In [76]: nodes, edges = ode_graph(eqs)
In [77]: strongly_connected_components((nodes, edges))
Out[77]: [[x(t)], [y(t), z(t)], [w(t)]]This shows that we can solve this system by:
- solving the ODE for
xand substituting that into the remaining equations - solving the ODEs for
yandzas a system and then substituting into the remaining equation - solving the ODE for
w.
Using this approach there can be many cases of a linear system of ODEs coupled to a nonlinear ODE that we can solve separately as a single ODE. Also a large linear system with a sparse matrix can often be solved efficiently in stages.
There is one situation where we might not want to use this approach which is if the original system is constant coefficient homogeneous then it can be solved without integrals. Breaking into strongly connected components makes subsequent systems become non-homogeneous requiring the use of integrals to find a solution. Such a system (being constant coefficient) would necessarily have manageable integrals but it might still be preferable to keep it as a whole system rather than partition in this way. Leaving aside this caveat breaking into strongly connected components should always be the next step after breaking into connected components so this should probably be done before any attempt to match or solve the system of ODEs.
There are 10 nonlinear solvers in dsolve for systems. The first 4 of these are:
Each of these boils down to the same thing so I'll explain how to generalise that here.
In general a system of two ODEs in two dependent variables looks like:
dx/dt = f(x, y, t)
dy/dt = g(x, y, t)
If the RHSs do not explicitly depend on t then we have an autonomous system which looks like
dx/dt = f(x, y)
dy/dt = g(x, y)
In any such autonomous system of two ODEs we can divide one ODE by the other to obtain a single ODE in which one or other of the original dependent variables becomes the independent variable. This gives two possibilities:
dx/dy = f(x, y)/g(x, y) := h(x, y)
or
dy/dx = g(x, y)/f(x, y) := p(x, y)
Either way that gives a single ODE that can be solved on its own. The solution gives x in terms of y or vice-versa and we can use that to eliminate x or y from the original ODEs to obtain a single ODE for x over t or for y over t. That means that the whole autonomous system of 2 coupled ODEs can be solved by solving 2 separate single ODEs one after the other.
More generally we don't need the ODEs to be autonomous provided the dependence on t cancels out during division e.g.:
dx/dt = f(x, y)*h(t)
dy/dt = g(x, y)*h(t)
The examples from the tests are:
eq1 = (Eq(diff(x(t),t),x(t)*y(t)**3), Eq(diff(y(t),t),y(t)**5))
eq2 = (Eq(diff(x(t),t), exp(3*x(t))*y(t)**3),Eq(diff(y(t),t), y(t)**5))
eq3 = (Eq(diff(x(t),t), y(t)*x(t)), Eq(diff(y(t),t), x(t)**3))
eq4 = (Eq(diff(x(t),t),x(t)*y(t)*sin(t)**2), Eq(diff(y(t),t),y(t)**2*sin(t)**2))
eq6 = (Eq(diff(x(t),t),x(t)**2*y(t)**3), Eq(diff(y(t),t),y(t)**5))
Note firstly that eq1, eq2, and eq6 can all be solved in stages (see strongly connected components above). Secondly the solutions currently given by dsolve for eq3 and eq4 fail the checksysodesol test and do not appear to be correct.
Let's take a simple example that works to demonstrate the method:
In [132]: x, y = symbols('x, y', cls=Function)
In [133]: z, t = symbols('z, t')
In [134]: eqs = (Eq(diff(x(t),t), y(t)*x(t)), Eq(diff(y(t),t), x(t)*y(t)))
In [135]: eqs
Out[135]:
⎛d d ⎞
⎜──(x(t)) = x(t)⋅y(t), ──(y(t)) = x(t)⋅y(t)⎟
⎝dt dt ⎠
In [136]: neweq_yz = Eq(y(z).diff(z), (eqs[1].rhs/eqs[0].rhs).subs(y(t), y(z)).subs(x(t), z))
In [138]: neweq_yz
Out[138]:
d
──(y(z)) = 1
dz
In [139]: yz = dsolve(neweq_yz, y(z))
In [140]: yz
Out[140]: y(z) = C₁ + z
In [141]: yt = yz.subs(y(z), y(t)).subs(z, x(t))
In [142]: yt
Out[142]: y(t) = C₁ + x(t)
In [143]: neweq_x = eqs[0].subs(yt.lhs, yt.rhs)
In [144]: neweq_x
Out[144]:
d
──(x(t)) = (C₁ + x(t))⋅x(t)
dt
In [145]: solx = dsolve(neweq_x, x(t))
In [146]: solx
Out[146]:
C₁⋅(C₂ + t)
C₁⋅ℯ
x(t) = ────────────────
C₁⋅(C₂ + t)
1 - ℯ
In [147]: soly = yt.subs(solx.lhs, solx.rhs)
In [148]: soly
Out[148]:
C₁⋅(C₂ + t)
C₁⋅ℯ
y(t) = C₁ + ────────────────
C₁⋅(C₂ + t)
1 - ℯ
In [24]: checksysodesol(eqs, [solx, soly])
Out[24]: (True, [0, 0])Currently dsolve hits a bug for this example (in the type3 solver):
In [3]: dsolve(eqs)
---------------------------------------------------------------------------
TypeError: 'Equality' object is not iterableThese 4 nonlinear 2 equation solvers could all be replaced by a single more general solver along the lines above. This method can also be generalised to 3 or more ODEs. The way it works is if we have
dx/dt = f(x, y, z) p(t)
dy/dt = g(x, y, z) p(t)
dz/dt = h(x, y, z) p(t)
Then we can eliminate t from any pair e.g.
dx/dy = f/g
This is only useful though if we end up with a rhs that depends only on x and y so that we can solve a single ODE for x and y. Otherwise we can in general at least reduce the system of 3 ODEs to 2 ODEs e.g.:
dx/dy = f(x, y, z)/g(x, y, z)
dz/dy = h(x, y, z)/g(x, y, z)
This system would then potentially be solvable. There are 3 ways to reduce this system since we can choose any of x, y, or z to become the independent variable.
We could use this method to also handle some of the 3-eq nonlinear solvers. Altogether this amounts to replacing the following nonlinear solvers with a single method:
- https://docs.sympy.org/dev/modules/solvers/ode.html#nonlinear-2-equations-order-1-type-1
- https://docs.sympy.org/dev/modules/solvers/ode.html#nonlinear-2-equations-order-1-type-2
- https://docs.sympy.org/dev/modules/solvers/ode.html#nonlinear-2-equations-order-1-type-3
- https://docs.sympy.org/dev/modules/solvers/ode.html#nonlinear-2-equations-order-1-type-4
- https://docs.sympy.org/dev/modules/solvers/ode.html#nonlinear-3-equations-order-1-type-1
- https://docs.sympy.org/dev/modules/solvers/ode.html#nonlinear-3-equations-order-1-type-2
These 3 methods work along similar lines but require first identifying an invariant (all 3 of these solvers are untested and broken):
- https://docs.sympy.org/dev/modules/solvers/ode.html#nonlinear-3-equations-order-1-type-3
- https://docs.sympy.org/dev/modules/solvers/ode.html#nonlinear-3-equations-order-1-type-4
- https://docs.sympy.org/dev/modules/solvers/ode.html#nonlinear-3-equations-order-1-type-5
The only remaining nonlinear solver is this one:
That has a single test:
eq5 = (Eq(x(t),t*diff(x(t),t)+diff(x(t),t)*diff(y(t),t)), Eq(y(t),t*diff(y(t),t)+diff(y(t),t)**2))
which could in principle be solved in stages (strongly connected components).
It is worth noting that many of these nonlinear solvers are untested and don't work. Most also represent relatively obscure classes of ODE system that are unlikely to come up in practical problems.