ODE Systems rewrite

[Oscar Benjamin]

Hi all,

I've written a roadmap for rewriting the ODE systems code:
https://github.com/sympy/sympy/wiki/ODE-Systems-roadmap

The roadmap represents a bonfire of the existing code and proposes to
delete all 33 of the current systems solvers. Around 7 of those are
untested and entirely broken and would be removed. The remaining 26
would be replaced by 4-5 more general (and less buggy) solvers. This
would greatly simplify the code for systems of ODEs at the same time
as dramatically increasing its capabilities.

There are 23 solvers for linear systems of ODEs and the roadmap
proposes to delete 3 broken ones and replace the other 20 with 3-4
solvers for the following cases:

Constant coefficient homogeneous (vector X and matrix A):

X' = A X

Constant coefficient nonhomogeneous (vector f):

X' = A X + f(t)

Non-constant homo/nonhomo where A(t) commutes with antiderivative B(t):

X' = A(t) X + f(t)

There should also be code to reduce higher order systems to 1st order
form (reduce to the cases above) and to partition systems of ODEs into
weakly and strongly connected components.

Most of the 10 nonlinear systems solvers are broken. Those that work
mostly come under the 2-equations autonomous class

x' = f(x, y)
y' = g(x, y)

where the solution is found from dx/dy = f/g so there should be a
solver for that case.

The work described would make an excellent GSOC project for anyone
interested. It corresponds to this idea on the GSOC ideas page:
https://github.com/sympy/sympy/wiki/GSoC-2020-Ideas#systems-of-ordinary-differential-equations

Whether this gets taken up as a GSOC project or not I hope that the
roadmap will help contributors thinking about improvements to sympy's
ODE capabilities. Most contributors looking to work on ODEs in sympy
seem to focus on improvements in the single ODE case. Those are very
useful but right now the systems code needs improvements in capability
more than the single ODE code does. There are easy wins for systems
that would make sympy's ODE handling much more useful in practical
problems.

Any feedback on the roadmap is appreciated. It is part of the wiki so
anyone can edit it (do you need push access to edit the wiki?).

Oscar

[Aaron Meurer]

Thanks for the writeup. I agree matrix exponential is the right
approach. It also takes care of issues with symbolic parameters where
there are bifurcations on the geometric multiplicity. For example

>>> pprint(exp(Matrix([[x, 1], [0, 1]])))
⎡ x ⎤
⎢ x ℯ ℯ ⎥
⎢ℯ ───── - ─────⎥
⎢ x - 1 x - 1⎥
⎢ ⎥
⎣0 ℯ ⎦

If x = 1 the matrix is not diagonalizable. However, if you take the
limit as x -> 1 of the exponential expression it gives the right
answer. I believe this should always work since the exponential is
continuous and nondiagonalizability occurs at isolated points (if I'm
wrong on this please let me know).

The current systems solver, if I understand it correctly, tries to
split things into Piecewise conditions. But there are combinatorially
too many possibilities to use Piecewise for systems of n equations.
Giving a solution that is correct in the limit is sufficient.

Aaron Meurer

> --
> You received this message because you are subscribed to the Google Groups "sympy" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to sympy+un...@googlegroups.com.
> To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/CAHVvXxR%3D42pCuvK0OJnBwcoCh5ajsZ6s8e_Soj7eevwpTmyDAw%40mail.gmail.com.

[Oscar Benjamin]

I actually think that the Piecewise approach is important and is
something that we should do somehow but for it to work it needs to be
done (or doable) *everywhere* including solvers for linear systems and
polynomials:
https://github.com/sympy/sympy/issues/16861
We could have that in dsolve but it should inherit work done in
matrices/poly etc rather than implementing all the cases in the ode
module.

For ODE systems right now only the 2-equation constant coefficient
solver uses Piecewise and it has lots of bug reports. Basic things
explode because it isn't implemented very well. I won't show the
output but try this in isympy:

dsolve([f(x).diff(x)-y*f(x), g(x).diff(x)-g(x)])

Given that it doesn't work properly and is inconsistent with the rest
of solvers the Piecewise behaviour there should go and we should
return results for the generic cases.

The matrix exponential is continuous (and analytic). Whether or not
non-diagonalisability occurs at isolated points is more complicated
though as it depends which class of matrices we consider in the first
place which depends on user input in sympy's case.

Given a matrix [[a, b], [c, d]] where a, b, c, d are real we have a 4D
space of matrices. The submanifold of non-diagonalisable matrices is
3D so it is a non-generic subset of the 4D space. However if we think
of the submanifold of matrices that have repeated eigenvalues then
within *that* space non-diagonalisability is generic and
diagonalisability occurs in a 2D subspace (the symmetric case [[a, b],
[b, a]]). On the other hand if we restrict ourselves to symmetric
matrices then we get a 3D subspace [[a, b], [b, c]] within which
diagonalisability is universal.

Oscar

> To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/CAKgW%3D6K2p1P4swGCtu619TCxqys%3DfRWx4pbuFwa2o_x255obcw%40mail.gmail.com.