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 On Mon, Feb 24, 2020 at 3:31 PM Oscar Benjamin <[email protected]> wrote: > > 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 > > -- > 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 [email protected]. > To view this discussion on the web visit > https://groups.google.com/d/msgid/sympy/CAHVvXxR%3D42pCuvK0OJnBwcoCh5ajsZ6s8e_Soj7eevwpTmyDAw%40mail.gmail.com. -- 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 [email protected]. To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/CAKgW%3D6K2p1P4swGCtu619TCxqys%3DfRWx4pbuFwa2o_x255obcw%40mail.gmail.com.
