[Discussion] GSoC 2020: Improving and extending ODE

[Milan Jolly]

Hello everyone,

My name is Milan Jolly and I am an undergraduate student at Indian Institute of Technology, Patna. For the past 2 month, I have been learning and exploring sympy through either contributions, reading documentation or trying examples out. This last month I have learned a lot of new things thanks to the well designed code-base, the structured way this community works and most importantly the maintainers who make it work. It has been a pleasure to be a part of the community. 

I am interested in participating for GSoC this year and I would like to work for this org during the summers if I am lucky. I particularly want to work on improving the current ODE module as it is given in the idea list. There is a lot of work that needs to be taken care of like:

1. Implementing solvers for solving constant coefficient non-homogeneous systems

2. Solving mixed order ODEs

3. Adding rearrangements to solve the system

These are not my ideas but I have taken inspiration from the ideas page but I am up for working on these. If someone can guide me regarding this then it would be really helpful.

[Aaron Meurer]

Are you familiar with Oscar's efforts to refactor the ODE module?
https://github.com/sympy/sympy/pull/18403

You should at least be aware of that, if not including it in your
proposal. The ODE module is starting to collapse under the weight of
itself.

I'm not clear what you mean by rearrangements, if this is the same
thing, but another idea to be aware of is the idea of splitting a
system using its strongly connected components
https://github.com/sympy/sympy/issues/16174. However, this is the sort
of thing that could get very messy the way the code is currently
rewritten, which is why a refactor is so important.

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/39adc638-2da1-4d4b-ada8-3cff7cf188bc%40googlegroups.com.

[Oscar Benjamin]

I have started writing a roadmap and intro to the maths for dsolve for
systems here:
https://github.com/sympy/sympy/wiki/ODE-Systems-roadmap

> To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/CAKgW%3D6KUGaBU-k0z9iEEK%2BauZqf5zpjCWS-E_rtBNjMAyAS82w%40mail.gmail.com.

[Milan Jolly]

Thanks Aaron Meurer and Oscar Benjamin. I will look into the links you shared and start working towards them right away.

To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/CAHVvXxSKa9mmNSrFeV8ZCLrUGGTddXL48EmdZo2KBcWPgEz6RQ%40mail.gmail.com.

[Milan Jolly]

Ok so I have gone through the links suggested and I have realised that as far as ODE module is concerned, refactoring is the most important task. But, as far as that is concerned, I think Mohit Balwani is working on this for a while and I want to limit any collisions with my co-contributors. So, I have couple of ideas to work on:

1. Helping to extend the solvers, i.e.implementing a fully working n-equations solver for constant coefficient homogeneous systems. This is from the ODE systems map. I am interested in working on this but I understand that it might be hard to work upon it while refactoring takes place. Still, if its possible to work on this and if no one else has started to work in this direction yet then I am willing to work for this.

2. Using connected components function implemented by Oscar Benjamin in https://github.com/sympy/sympy/pull/16225 to enhance ODE solvers and computing eigen values faster as mentioned here https://github.com/sympy/sympy/issues/16207 .

3. This idea is not mentioned in the ideas page and is something of my own. If there is anything possible, then I can also work on extending functions like maximum, minimum, argmax, argmin, etc in calculus module. I have been working on the issue https://github.com/sympy/sympy/pull/18550 and I think there is some scope to extend these functionalities.

[Oscar Benjamin]

The current refactoring effort applies only to the case of solving
*single* ODEs. The ODE systems code also needs to be refactored but
(in my opinion) needs a complete rewrite. That is what the roadmap is
about (it describes how to rewrite everything). The code for systems
of ODEs should also get refactored in the process but there is no need
to "refactor" it in its current form if it is in fact being
*completely* rewritten: we can just make sure that the new code is
written the way we want it to be.

> --
> 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/1033f581-abbb-4be5-a5b2-1988f4261535%40googlegroups.com.

[Milan Jolly]

So, I am interested in rewriting parts of the current ODE as discussed in the roadmap. Is there any work started in that direction and if not then can I create a PR for the same?

To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/CAHVvXxTeWturK6WmtHKxakLqbV1yhp5_KoTPs8vPtmbu8%3D2VxQ%40mail.gmail.com.

[Oscar Benjamin]

This section in the roadmap refers to existing stalled PRs trying to
fix the n-equations solver for constant coefficient homogeneous ODEs
which is the first step:
https://github.com/sympy/sympy/wiki/ODE-Systems-roadmap#constant-coefficients---current-status

A first step would be to attempt to revive one or both of those PRs
and finish them off.

> To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/CAMrWc1BZKcwjFvVQwZZ80P6s-CP62Tn_VN30zBEdjmqhrN44DA%40mail.gmail.com.

[Milan Jolly]

I will go through the roadmap. Also, I will work on reviving and finishing the stalled PRs namely the matrix exponential one for now as I am interested in working towards this. Thanks.

To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/CAHVvXxRGFJC%2BGeowaNJcVFQsQsGq%2Bh0SW-jmYezsSzU5%3DvipcA%40mail.gmail.com.

[Milan Jolly]

I have  mostly read and understood matrix exponentials and Jordan forms along with the ODE systems roadmap. But I am unclear as to what has already been done when it comes to implementing the general solvers. For example: The matrix exponentials part has already been implemented and now I have a PR that has revived the matrix exponential code.

I want to make a proposal and contribute to make these general solvers during this summer if my proposal gets accepted. But I am unclear what should be the parts I need to work during community bonding period, phase 1, phase 2 and phase 3 as I am unaware how much time each part of the general solvers would take. 

If someone can help me in this regard(helping me with these 2 questions) then it would be great.

[Oscar Benjamin]

It's hard to say how much time each of these would take. The roadmap
aims to completely replace all of the existing code for systems of
ODEs. How much of that you think you would be able to do is up to you
if making a proposal.

None of the other things described in the roadmap is implemented
anywhere as far as I know. Following the roadmap it should be possible
to close all of these issues I think:
https://github.com/sympy/sympy/issues?q=is%3Aopen+is%3Aissue+label%3Asolvers.dsolve.system

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

[Nicolas Guarin]

I have been working in a function that turn a higher-order system of ODEs into a system of first order equations. So you think that it might help?

[Oscar Benjamin]

On Sun, 15 Mar 2020 at 14:45, Nicolas Guarin <nicog...@gmail.com> wrote:
>
> I have been working in a function that turn a higher-order system of ODEs into a system of first order equations. So you think that it might help?

Absolutely. This is mentioned in the roadmap:
https://github.com/sympy/sympy/wiki/ODE-Systems-roadmap#higher-and-mixed-order-systems-of-odes

I think that being able to rewrite single higher-order ODEs in the
form of a first order system is also useful since that's the form
expected by most numerical integration routines.

--
Oscar

[Milan Jolly]

Thanks for your reply. I have planned a rough layout for the phases. I took a lot of time this past month to understand all the mathematics that will be involved and have grasped some part of it.

If I am lucky and get selected for GSOC'20 for this organisation, then the below is the rough plan. Please comment on suggestions if necessary.

Community Bonding phase:

1. Using matrix exponential to solve first order linear constant coefficient homogeneous systems(n equations).

2. Adding new test cases and/or updating old ones.

3. Removing and closing related issues if they are solved by the addition of this general solver. Identifying and removing the special cases solvers which are covered by this general solver.

Phase I:

1. Adding technique to solve first order constant coefficient non-homogeneous systems(n equations).

2. Adding the functionality that reduces higher order linear ODEs to first order linear ODEs(if not done already, and if done, then incorporating it to solve higher order ODEs).

3. Adding a special case solver when non-constant linear first order ODE has symmetric coefficient matrix.

Phase II:

1. Adding technique to solve non-constant non-homogeneous linear ODE based off the solver added by the end of Phase I.

2. Evaluating and eliminating unnecessary solvers.

3. Closing related issues solved by the general solvers and identifying and removing unwanted solvers.

4. Adding basic rearrangements to simplify the system of ODEs.

Phase III:

1. Dividing the ODEs by evaluating which sub-systems are weakly and strongly connected and handling both of these cases accordingly.

2. Adding a special case solver where the independent variable can be eliminated and thus solving the system becomes easier.

3. Wrapping things up: adding test cases, eliminating unwanted solvers and updating documentation.

This is the rough layout and my plan for summer if I get selected. If this plan seems ok then I would include this plan in my proposal.

>>> >> > To unsubscribe from this group and stop receiving emails from it, send an email to sy...@googlegroups.com.

>>> >> > To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/1033f581-abbb-4be5-a5b2-1988f4261535%40googlegroups.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 sy...@googlegroups.com.

>>> >> To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/CAHVvXxTeWturK6WmtHKxakLqbV1yhp5_KoTPs8vPtmbu8%3D2VxQ%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 sy...@googlegroups.com.

>>> > To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/CAMrWc1BZKcwjFvVQwZZ80P6s-CP62Tn_VN30zBEdjmqhrN44DA%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 sy...@googlegroups.com.

>>> To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/CAHVvXxRGFJC%2BGeowaNJcVFQsQsGq%2Bh0SW-jmYezsSzU5%3DvipcA%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 sy...@googlegroups.com.

[Oscar Benjamin]

That seems reasonable to me. Since the plan is a total rewrite I think
that it would be good to put some time in at the beginning for
designing how all of these pieces would fit together. For example even
if the connected components part comes at the end it would be good to
think about how that code would fit in from the beginning and to
clearly document it both in issues and in the code.

Getting a good design is actually more important than implementing all
of the pieces. If the groundwork is done then other contributors in
future can easily implement the remaining features one by one. Right
now it is not easy to improve the code for systems because of the way
that it is structured.

> 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/b2d61327-da27-4baf-a753-cf4069abfc78%40googlegroups.com.

[Oscar Benjamin]

There are possibilities to go from nonlinear to linear e.g.:

In [6]: x, y = symbols('x, y', cls=Function)

In [7]: eqs = [x(t).diff(t)**2 - y(t)**2, y(t).diff(t)**2 - x(t)**2]

In [8]: eqs
Out[8]:
⎡ 2 2⎤
⎢ 2 ⎛d ⎞ 2 ⎛d ⎞ ⎥
⎢- y (t) + ⎜──(x(t))⎟ , - x (t) + ⎜──(y(t))⎟ ⎥
⎣ ⎝dt ⎠ ⎝dt ⎠ ⎦

In [9]: solve(eqs, [x(t).diff(t), y(t).diff(t)], dict=True)
Out[9]:
⎡⎧d d ⎫ ⎧d d ⎫
⎧d d ⎫ ⎧d d ⎫⎤
⎢⎨──(x(t)): -y(t), ──(y(t)): -x(t)⎬, ⎨──(x(t)): -y(t), ──(y(t)):
x(t)⎬, ⎨──(x(t)): y(t), ──(y(t)): -x(t)⎬, ⎨──(x(t)): y(t), ──(y(t)):
x(t)⎬⎥
⎣⎩dt dt ⎭ ⎩dt dt ⎭
⎩dt dt ⎭ ⎩dt dt ⎭⎦

On Mon, 16 Mar 2020 at 15:48, Milan Jolly <milan....@gmail.com> wrote:
>
> Thanks for the suggestion, I have started with the design for these solvers. But I have one doubt, namely since now we are using linear_eq_to_matrix function to check if the system of ODEs is linear or not, would we require the canonical rearrangements part? Or rather are there other cases when we can reduce non-linear ODEs into linear ODEs.

>> > To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/b2d61327-da27-4baf-a753-cf4069abfc78%40googlegroups.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 sympy+un...@googlegroups.com.

> To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/38ad9558-b8e0-465a-aeb4-5ea782144b55%40googlegroups.com.

[Milan Jolly]

So, I have made a rough layout of the main function that will be used to solve ODEs with the methods like neq_nth_order_linear_constant_coeff_homogeneous/nonhomogeneous, neq_nth_linear_symmetric_coeff_homogeneous/nonhomogeneous, special case non-linear solvers, etc.

Some notations used:

   eqs: Equations, funcs: dependent variables, t: independent variable, wcc: weakly connected component, scc: strongly connected component

Introduction to helper functions that will be used(these are temporary names, parameters and return elements and may be changed if required):

1. match_ode:- 

        Parameters: eqs, funcs, t

        Returns: dictionary which has important keys like: order(a dict that has func as a key and maximum order found as value), is_linear, is_constant, is_homogeneous, eqs, funcs.

2. component_division:-

        Paramters: eqs, funcs

        Returns: A 3D list where the eqs are first divided into its wccs and then into its sccs.

        This function is suggested to be implemented later. So, until all the other solvers are not ready(tested and working), this function will just take eqs and return [[eqs]].

3. get_coeff_matrix:-

        Parameters: eqs, funcs

        Returns: coefficient matrix A(t) and f(t)

        This function takes in a first order linear ODE and returns matrix A(t) and f(t) from X' = A(t) * X + f(t).

4. nth_order_to_first_order:-

        Parameters: eqs, order

        Returns: first order ODE with new introduced dependent variables.

And all the first order linear solvers mentioned above.

Now, besides the main function, there are two separate functions depending on whether the system of ODEs is linear or not, namely _linear_ode_sol and _non_linear_ode_sol.

1. _first_order_linear_ode_sol:-

        Parameters: match dict(obtained earlier and maybe modified in ode_sol)

        Returns: Dict with keys as func and value as its solution that solves the ODE.

        Working: First, extracts A(t) and f(t) using get_coeff_matrix, then using match dict, identify which solver is required and solve the ODE if it is possible to do so. For example: Case where A(t) is not symmetric isn't solved.

2. _non_linear_ode_sol has similar Parameters and Returns but the function operates differently that's why it is essential to use a different function. But I don't have a clear understanding

   of how to design _non_linear_ode_sol yet but here is what I have came up with: First match the condition where it is possible seperate out the independent variable to get a relationship

   between the dependent variables and then finally, just use the special solver to solve the ODE.

Now, coming to the main function ode_sol(for now, I haven't considered initial values):-

  Parameters: eqs, funcs, t

  Returns: Solution in a dict form where func is the key and value is the solution for that corresponding func.

  Working: 

  The steps of its working-

      1. Preprocess the equations.

      2. Get the match dict using match_ode function.

      3. Convert nth order equations to first order equations using nth_order_to_first_order while storing the funcs seperately so that we can later filter out the dependent variables that were introduced in this step.

      4. Get the 3D list of equations using component_division function.

      5. Iterate through the wccs and solve and store solutions seperately but for sccs, first solve the first set of equations in a scc, then substitute the solutions found in the first set of the current scc to the second 

         set of current scc. Keep doing this until the all the sets for a particular scc is solved.

      6. For solving a component, choose either _linear_ode_sol or _non_linear_ode_sol depending upon the set of equations to be solved.

      7. Return a dict by taking out values from the solution obtained using all the dependent variables in funcs as there may be more variables introduced when we made the system into first order.

For now, this is what I have came up with. Obviously the order in which we will proceed is, build the basic layout of the main function and component_division will just increase the number of dimensions to 3 rudimentarily as we

will have to first ensure that the general solvers work well since working on both of them simultaneously will make it tough to pinpoint the errors. Along with that, non-linear solvers can be implemented later, we can just raise a 

NotImplementedError for now till we have completed both the general linear solvers and the component_division and then add the special case solvers.

> To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/38ad9558-b8e0-465a-aeb4-5ea782144b55%40googlegroups.com.

[Oscar Benjamin]

That sounds reasonable.

Note that we can't start raising NotImplementedError yet. You will
need to think about how to introduce the new code gradually while
still ensuring that dsolve falls back on the old code for cases not
yet handled by the new code.

> 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/debf5cf4-2817-43ab-b4ad-64c60bf94c1e%40googlegroups.com.

[Milan Jolly]

In ODE systems roadmap, you have mentioned that for system of ODEs where the coefficient matrix is non-constant, if the coefficient matrix A(t) is symmetric, then A(t) and its anti derivative B(t) commute and thus we get the solution based on this fact. But it is also mentioned that if A and B commuting is more general than when A is symmetric, that is, it is possible that A is not symmetric but A and B commute. So, for that solver, should we first compute its anti derivative and test it that commutes with A or just check if A is symmetric and use the solution?

> To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/debf5cf4-2817-43ab-b4ad-64c60bf94c1e%40googlegroups.com.

[Oscar Benjamin]

It's not always the case that symmetric matrices commute so actually
checking if it is symmetric is not sufficient e.g.:

In [83]: M = Matrix([[2*x**2, x], [x, x**2]])

In [84]: M.is_symmetric()
Out[84]: True

In [85]: M*M.diff(x) == M.diff(x)*M
Out[85]: False

Maybe there is something that can be said more generally about
`exp(M(t)).diff(t)` when `M` is symmetric but does not necessarily
commute with `M.diff(t)`...

> 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/07fab7ec-3854-4a29-8e33-6e0804177d49%40googlegroups.com.

[Milan Jolly]

Thanks for clearing my doubt.

Now, I have started preparing my GSOC proposal and it will be ready soon. But, I wanted to know that will it be ok that I don't give details about the implementations of the helper functions and solvers and simply state what they do, which parameters they take,  what they return and how they fit in the solving process while I give more details about how they fit together more generally. I would like to elucidate more on how the main function ode_sol handles the system of equations using the helper functions and various solvers as it is the only thing that is not clearly mentioned in the roadmap.

> To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/07fab7ec-3854-4a29-8e33-6e0804177d49%40googlegroups.com.

[Oscar Benjamin]

Stating clearly what the different parts do in high-level terms should
be sufficient.

> 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/c69dc11d-d2bc-4e57-83e2-e06c483886aa%40googlegroups.com.

[Milan Jolly]

I am almost done with my proposal but I wanted to ask something regarding the same. Now, I have made a theory section and an implementation section separately. In the theory section, I have given sufficient examples but the implementation section became very big and I haven't added examples of their working as I have just given the rough plan for all the functions that will be implemented. Is it mandatory to give some demonstrations for each of the functions(solvers, helpers, and the main ones)? 

To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/CAHVvXxSxwZevbwy%3DiOMe_awP%3Djoj6wuLQRoFuhmjjLfeeFv%3DWw%40mail.gmail.com.

[Milan Jolly]

Here is the link to my proposal: https://docs.google.com/document/d/12QN19LSjwEvYoSukyq-BWd76ZrI24FQuU0CGIOIx6Ww/edit?usp=sharing

> To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/c69dc11d-d2bc-4e57-83e2-e06c483886aa%40googlegroups.com.

[Oscar Benjamin]

I took a quick look. It's long so I didn't read it fully but it looks
good. There is a lot of detail about what you would do but perhaps the
motivation section can be strengthened. What does all of this mean for
end users etc? If you completed the work described then sympy's
capabilities for systems of ODEs would be expanded enormously.

> 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/b802dab5-7280-4744-b099-e6c5562ffe8b%40googlegroups.com.

[Milan Jolly]

Thank you for your feedback. I have added another paragraph in the Motivation section where the I have added how these new solvers are advantageous to the end users.

> To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/b802dab5-7280-4744-b099-e6c5562ffe8b%40googlegroups.com.

[Nicolas Guarin]

You can check what I have done about it in the following notebook:

https://nbviewer.jupyter.org/github/nicoguaro/notebooks_examples/blob/master/ode2sys.ipynb

Let me know how can I help with the development.

Best,
Nicolás

PS: Sorry for taking that much to answer.

[Oscar Benjamin]

Hi Nicolás,

There are some good things in that notebook. Qualitative analysis of
nonlinear ODEs is orthogonal to the work here but absolutely something
that would be great to have in sympy.

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 sympy+un...@googlegroups.com.

> To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/76c616d8-e55d-4770-a0c8-66769d9e0fcf%40googlegroups.com.

[Milan Jolly]

I took a quick look. I found the logic of ode2sys and ode_rewrite very resourceful since its required for component division and for reducing higher order ODEs. 

[Milan Jolly]

If there are no issues with the proposal or the timeline mentioned(which I will update soon due to the GSoC program timeline changes), then I am planning on submitting the proposal in 2 days. Feedback would be appreciated if possible.

[Oscar Benjamin]

I had a quick look and it seems reasonable.

> 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/40081db9-7195-4dd8-9925-102280f78ed2%40googlegroups.com.

[Milan Jolly]

Thanks for your feedback. I will submit the proposal soon.

To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/CAHVvXxT5EGs0wtY-St4j0XpAd29Oz86ujsJMpyPwhm_NAY_4aQ%40mail.gmail.com.

[Milan Jolly]

Hello mentors, 

I have submitted my final proposal. This past month I was preoccupied with proposal making, learning theory and designing details of the project but I will resume my contributions soon. I hope this proposal helps in further designing of the module regardless of my selection for the program.

doc link: https://docs.google.com/document/d/1L7WgLGU_4oKfOivEEFk93MKbueNEZj2_IWFqUsEJmgQ/edit?usp=sharing