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

On Mon, 24 Feb 2020 at 23:13, Aaron Meurer <[email protected]> wrote:
>
> 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.

-- 
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/CAHVvXxRVXfo4Ov4KcbBsjM9TY5BaXWVoYO_qZV0DNidmDk%3DGxQ%40mail.gmail.com.

Reply via email to