ODE solver roadmap

[nijso.be...@gmail.com]

A couple of ideas for improving the ODE solver capabilities are listed on the GSoC project ideas page:

https://github.com/sympy/sympy/wiki/GSoC-Ideas#other-ode-ideas

I propose a roadmap for a (first order) ODE solver that leads to a structured and modular workflow where the individual functionalities can be reused for other, more advanced solvers. Currently, the first order ODE solver consists of a collection of solvers based mainly on pattern matching for well-known ODEs. More general methods for first order and second order ODEs exist that do not depend on pattern matching, but on searching for integrating factors or symmetry generators (that construct an integrating factor) of a certain type.

For first order ODEs, the key is that an integrating factor will solve the ODE. You can either search directly for integrating factors of a certain type, or you can search for a symmetry generator that will then in turn lead to an integrating factor.

Darboux methods search directly for integrating factors of a certain type. This method will find all solutions of first order ODEs with Liouvillian solutions, which is a very broad class. The idea is that you assume a polynomial shape for the integrating factor of a certain degree and then find the coefficients of using the method of undetermined coefficients. This method works well, but you need a solver for nonlinear systems of equations, so a Groebner basis computation, like the F5 algorithm would be needed in general.

Another approach is symmetry methods. A symmetry generator of the ODE is found, which directly leads to an integrating factor. Once the integrating factor is found, you can construct the first integral of the ODE, which is an implicit solution. If possible, you can then try to write the first integral as an explicit solution of the ODE. In general, we would also need Groebner basis computations for finding symmetries, but we can avoid this by implementing searches for specific symmetry generators.

So the workflow for the symmetry method is:

1. compute symmetry generator using several available options

2. compute integrating factor from symmetry generator

3. compute first integral using integrating factor

4. write first integral as explicit solution of the ODE

ODE solvers that can fit into this workflow are solvers for linear, inverse linear, exact, separable, Bernoulli, Abel (constant coefficient) and homogeneous ODEs

Some first order ODEs need special treatment, like Riccati ODEs and ODEs that have special functions as the solution (Bessel, Airy,...). That is because you need to solve the ODE first to find the symmetry generators. A solver for Rational Riccati ODEs will be important because many problems can be reduced to the solution of a Riccati ODE, like for instance second order ODEs

Performance and robustness is important. We need a regression test containing a large collection of ODEs for testing. The books of Kamke and Murphy contain such lists.

So I propose the following projects, which could be part of a GSOC project:

1. Create regression test for the Kamke and Murphy ODEs

2. Create workflow as stated above, and write auxiliary functions like computeFirstIntegralFromIntegratingFactor, testIntegratingFactor,  etc

3. Implement more structured symmetry search for well known ODEs, and make existing symmetry searches robust

4. Add more advanced symmetry search like from the paper of Kolokolnikov (https://arxiv.org/abs/math-ph/0007023)

5. Add rational Riccati solver

With this basis, methods can then be added that use this functionality:

6. Kovacic method for second order ODEs with Liouvillian solutions

7. Solver for rational first order ODEs

9. Implicit and higher degree first order ODEs

10. Computation of Janet bases for general linear second order ODEs.

11. Symmetry method for second order ODEs

12. Second order nonlinear ODEs

13. Computation of Invariant Algebraic Curves

Let me know what you think of this plan. All we need now is somebody who will do all of this :-)

Best regards,

Nijso

[Oscar Benjamin]

On Tue, 16 Mar 2021 at 09:27, nijso.be...@gmail.com
<nijso.be...@gmail.com> wrote:
>
>
> A couple of ideas for improving the ODE solver capabilities are listed on the GSoC project ideas page:
> https://github.com/sympy/sympy/wiki/GSoC-Ideas#other-ode-ideas
>
> I propose a roadmap for a (first order) ODE solver that leads to a structured and modular workflow where the individual functionalities can be reused

Hi Nijso,

I completely agree with all of your very useful suggestions. What
would be most useful is if you could write the ideas up in more detail
e.g. on the sympy wiki.

Last year I wrote something for improving the ODE systems solvers
which you can see here:
https://github.com/sympy/sympy/wiki/ODE-Systems-roadmap

That led to a successful GSOC project that massively improved the very
messy code that was there before for ODE systems. Most of the ideas
from that roadmap got implemented and thousands of lines of buggy and
repetitive code were replaced with some better more generalised code.
The key point is that we are more likely to get well-prepared GSOC
applications if the proposal for what needs doing is written up
clearly enough with examples/references. (It doesn't need to be as
explicit as the road map I wrote but we should try to be clear about
exactly which paper to read, where to see existing implementations
etc, how much would make sense for a single project, what background
needed etc)

Oscar

[nijso.be...@gmail.com]

Hi Oscar,

I've started with a document: https://github.com/sympy/sympy/wiki/ODE-solver-roadmap

I think the first thing to do is build a kind of regression test using e.g. the database of first order ODEs from the book of Kamke and see where we stand. Maybe somebody already did this? I started with making my own regression test a couple of months ago but other projects demanded my attention. I'm not so familiar with making regression tests, so I'll probably get back to you about how to set this up properly (I think I wanted to get much more output like time statistics I have to check).

There is a nice website listing the performance of maple and mathematica for the ODE solvers:

http://www.12000.org/my_notes/kamek/mma_12_1_maple_2020/KEse1.htm#x3-20001

An ODE solver that I wrote (for the maxima CAS) based on the symmetry papers from Cheb-Terrab and from Schwarz gave me a performance of around 86% for the first order, first degree ODEs, but very important classes are missing like solutions to Riccati and Bessel problems.

code:

https://github.com/bigfooted/maxima-odesolve/blob/master/ode1_lie.mac

with documentation here:

https://nbviewer.jupyter.org/github/bigfooted/maxima-odesolve/blob/master/Doc/ode1solve.ipynb

So by just making the existing symmetry solvers more robust we should get a similar performance without the solver getting stuck or crashing (I got 1 stuck and 1 crash, both in the Risch integration routine)

Best regards,

Nijso

[Oscar Benjamin]

On Wed, 17 Mar 2021 at 13:54, nijso.be...@gmail.com
<nijso.be...@gmail.com> wrote:
>
> I've started with a document: https://github.com/sympy/sympy/wiki/ODE-solver-roadmap

Yes, that looks good as an overall roadmap. It's written at the level
that I can understand although I'd have to take time at some point to
read some of the references. For GSOC students we might need to break
it down into more explicit tasks but that should be done on the GSOC
ideas page rather than in the roadmap though. A good description like
this is really helpful for other contributors as well. It's good if
the document can be as clear as possible with examples and can make it
clear what the current state of the solvers in sympy is (e.g. there is
a Lie group solver but I'm guessing it doesn't come close to what
you're describing in the roadmap).

> I think the first thing to do is build a kind of regression test using e.g. the database of first order ODEs from the book of Kamke and see where we stand. Maybe somebody already did this? I started with making my own regression test a couple of months ago but other projects demanded my attention. I'm not so familiar with making regression tests, so I'll probably get back to you about how to set this up properly (I think I wanted to get much more output like time statistics I have to check).

I can certainly help if you are unsure how to make regression tests!

Mohit Balwani has recently refactored the ODE test suite into more of
a database format which you can see here:
https://github.com/sympy/sympy/blob/master/sympy/solvers/ode/tests/test_single.py

The examples there are not really systematic though. Mostly they come
from either the author who added a particular method choosing examples
that happen to work with that method or otherwise bug reports from
users. There can be a reasonable rationale I think in having a
database of examples from users coming in as bug reports since it
shows what kinds of things people are trying to use the solver for but
it is not like having the Kamke examples as a database. I would be
interested just to see how dsolve currently measures.

> There is a nice website listing the performance of maple and mathematica for the ODE solvers:
> http://www.12000.org/my_notes/kamek/mma_12_1_maple_2020/KEse1.htm#x3-20001
>
> An ODE solver that I wrote (for the maxima CAS) based on the symmetry papers from Cheb-Terrab and from Schwarz gave me a performance of around 86% for the first order, first degree ODEs, but very important classes are missing like solutions to Riccati and Bessel problems.
> code:
> https://github.com/bigfooted/maxima-odesolve/blob/master/ode1_lie.mac
> with documentation here:
> https://nbviewer.jupyter.org/github/bigfooted/maxima-odesolve/blob/master/Doc/ode1solve.ipynb
>
> So by just making the existing symmetry solvers more robust we should get a similar performance without the solver getting stuck or crashing (I got 1 stuck and 1 crash, both in the Risch integration routine)

I'll need to take some time to read through these things at some point.

Your input on this would be extremely helpful as we don't have anyone
involved in the ode module who has the level of experience that you
have. At the moment most work is on improving the pattern matching
system which is itself slow and badly organised.


Oscar

[Naveen Saisreenivas Thota]

Hi Nijso,

I went through the roadmap and found that the paper  Rational and algebraic solutions of first-order algebraic ODEs has a very detailed and simple-to-understand approach for finding solutions of the Riccati Equation. However, there are 2 algorithms in the paper - Algorithm 4 (Pg 59) which computes rational general solutions and Algorithm 11 (Pg 82), which computes rational solutions of the equation. Which one is more general and is better to be implemented?

For the ODE test database, 12000.org has a listing of problems at https://www.12000.org/my_notes/kamek/maple_listing.txt for Kamke ODEs. If we require metadata about the ODEs, there is a general classification given in Table 2. For detailed info, we could even parse the HTML pages for each ODE. A similar structure exists for Murphy ODEs, so it may not be very difficult to port these examples into SymPy using parse_expr or even the latex parser.

Naveen

[nijso.be...@gmail.com]

Hi Naveen,

Algorithm 11 is more general, as it is not limited to finding the rational general solution. I'm not sure how many cases will be missed by using algorithm 4, but in any case algorithm 11 is the one that is needed also for an extension to Kovacic algorithm. Preferably, algorithm 11 should then be set up in such a way that Kovacic method can easily call it, an alg.11 cal also later be called by a solver for first order rational ODEs that are not Riccati.

the algorithm is also in the paper:

https://www3.risc.jku.at/publications/download/risc_5269/RISCReport1601.pdf

For the Kamke database, we would need a list of ODEs and a list of solutions, and we need to find out which of these ODEs are directly solvable by the method. We might need additional ODEs from the paper, thesis, or our own to cover all cases.

Best regards,

Nijso

[Naveen Saisreenivas Thota]

Hi Nijso,

I have some doubts regarding Algorithm 11. Please help me understand them -

1. In case 1 and 3 (Pg 79-80), is the value of r_i = 1?

2. How do we find the valuation v_inf(a(x))?

3. In step 7 of Algorithm 11 (Pg 82), it says we have to do the following for "all possible combinations of these vectors". Since we have 1 vector (d_-1, d_0, ..., d_N) and 2 vectors (c_i1, c_i2, ...) for each singular point, does all possible combinations mean 2^n combinations (where n is the number of singular points) as we can pick one of the two for each singular point?

4. The paper mentions that all polynomial solutions of degree m for the equation 5.16 (Pg 81) can be done using Linear Algebra. Could you please elaborate on any algorithm to do so?

Lastly, finding the coefficient vectors requires Laurent Series expansion. I'm not sure if the series module can achieve this. There seems to be a function laurent_series, but I couldn't understand the documentation. Perhaps Oscar could help us out here.

Naveen

[nijso.be...@gmail.com]

1. in case 1 r_1=1, in case 3 r_i is not relevant since it's a simple pole.
2. The definition of the valuation at infinity is mentioned in section 5.1
3. yes, I think so
4. I think he just means using the method of undetermined coefficients here

Some ideas are also explained in the first chapter of Fritz Schwarz' book: Algorithmic Lie theory for solving Ordinary Differential Equations. I think this chapter gives a good introduction, although it is going through the material at a high speed.

Best, Nijso

[Oscar Benjamin]

On Tue, 23 Mar 2021 at 12:45, Naveen Saisreenivas Thota
<naveensai...@gmail.com> wrote:
>
> I have some doubts regarding Algorithm 11. Please help me understand them -

...

>
> Lastly, finding the coefficient vectors requires Laurent Series expansion. I'm not sure if the series module can achieve this. There seems to be a function laurent_series, but I couldn't understand the documentation. Perhaps Oscar could help us out here.

Is the Laurent series actually needed?

I haven't read the paper but I looked at algorithm 11. Step 5 says
"compute the poles of a(x) and their orders". Is it not just asking
for the partial fraction expansion?

Oscar

[Aaron Meurer]

On Wed, Mar 17, 2021 at 7:54 AM nijso.be...@gmail.com
<nijso.be...@gmail.com> wrote:
>
> Hi Oscar,
>
> I've started with a document: https://github.com/sympy/sympy/wiki/ODE-solver-roadmap
>
> I think the first thing to do is build a kind of regression test using e.g. the database of first order ODEs from the book of Kamke and see where we stand. Maybe somebody already did this? I started with making my own regression test a couple of months ago but other projects demanded my attention. I'm not so familiar with making regression tests, so I'll probably get back to you about how to set this up properly (I think I wanted to get much more output like time statistics I have to check).
>
> There is a nice website listing the performance of maple and mathematica for the ODE solvers:
> http://www.12000.org/my_notes/kamek/mma_12_1_maple_2020/KEse1.htm#x3-20001
>
> An ODE solver that I wrote (for the maxima CAS) based on the symmetry papers from Cheb-Terrab and from Schwarz gave me a performance of around 86% for the first order, first degree ODEs, but very important classes are missing like solutions to Riccati and Bessel problems.
> code:
> https://github.com/bigfooted/maxima-odesolve/blob/master/ode1_lie.mac
> with documentation here:
> https://nbviewer.jupyter.org/github/bigfooted/maxima-odesolve/blob/master/Doc/ode1solve.ipynb
>
> So by just making the existing symmetry solvers more robust we should get a similar performance without the solver getting stuck or crashing (I got 1 stuck and 1 crash, both in the Risch integration routine)

At the end of the day, the power of the integrator is always going to
be a limiting factor in the ODE solver. For those solvers that call
out to integrate(), I would split the comparisons into two steps, the
computation of the solution with an unevaluated integral, and the
computation of the integral. sympy.dsolve() already does this in the
hints system with separate _Integral hints, so this will not be hard.
I know Maple has a similar feature (the dsolve hints feature is based
on Maple). I don't know about Mathematica or Maxima.

The integration itself is obviously important for the actual final
answer, but any improvements to that step will generally need to be
done in the integration code, not the ODE solving code.

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/e1bc45fc-89dd-4c79-a351-354ccf9b4694n%40googlegroups.com.

[Naveen Saisreenivas Thota]

> 1. in case 1 r_1=1, in case 3 r_i is not relevant since it's a simple pole.

So, do we add 1/(x - x_i) to y(x) when we encounter case 3?

> 2. The definition of the valuation at infinity is mentioned in section 5.1

Yes, thank you. 

> 4. I think he just means using the method of undetermined coefficients here

Oh ok. SymPy already has a solver for it, so we'd just have to call it internally.

The algorithm is very clear to me now. Are there any examples which walk through the entire algorithm? If so, could you please link them here? I'll try to write a basic function which can accomplish it and then start working on the proposal.

Naveen

[Naveen Saisreenivas Thota]

> Is the Laurent series actually needed?

> I haven't read the paper but I looked at algorithm 11. Step 5 says
> "compute the poles of a(x) and their orders". Is it not just asking
> for the partial fraction expansion?
 

I think Laurent Series is required in step 6 where we have to calculate the coefficient vectors for each singular point.

Naveen

[Naveen Saisreenivas Thota]

I checked some examples of laurent series expansion and the series function seems to give correct results. Some examples I tested where from -

1. http://courses.washington.edu/ph227814/228/W14/notes/Laurent.nb.pdf

2. https://piazza.com/class_profile/get_resource/iw9kxycftxk6su/iy5oz4r144i6rr

3. https://www.maplesoft.com/support/help/Maple/view.aspx?path=numapprox/laurent

Naveen

[Naveen Saisreenivas Thota]

Hi Nijso,

Are there any examples which walk through the entire algorithm? If so, could you please link them here? I've written a pretty basic function which should work for Case 1, but it doesn't. I could include the code for it too if required.

Naveen

[nijso.be...@gmail.com]

Hi Naveen, have a look at this report from the RISC group,where they list the algebraic solutions of the kamke ODEs. Check if there are cases that fall under case 1, or you can construct a solution yourself. Take the solution y(x) that belongs to case 1, select a couple of poles construct a solution and write it as a Riccati ode.

https://www3.risc.jku.at/publications/download/risc_5197/RISCReport15-19.pdf

Also, try splitting the solver into a number of independent functions, like a function to determine poles, checking if a pole is really a pole, etc. 

Best,

Nijso

[Naveen Saisreenivas Thota]

Hi Nijso,

Thanks for linking the report!  I'll check it out. I figured out the errors in the code. For now, it seems to be working on some examples given in Fritz's book and Kovacic's paper. The code is here - Rational Riccati Solver. Let me know what you think of it. I'll test some of the other ODEs meanwhile.

> Also, try splitting the solver into a number of independent functions, like a function to determine poles, checking if a pole is really a pole, etc. 

Sure, I'll do that.

Naveen

[Oscar Benjamin]

On Sun, 28 Mar 2021 at 16:47, Naveen Saisreenivas Thota
<naveensai...@gmail.com> wrote:
>
> Hi Nijso,
>
> Thanks for linking the report! I'll check it out. I figured out the errors in the code. For now, it seems to be working on some examples given in Fritz's book and Kovacic's paper. The code is here - Rational Riccati Solver. Let me know what you think of it. I'll test some of the other ODEs meanwhile.

That looks nice. It would be easier to comment on the details if it was a PR.

One comment I have is that if this is all working with rational
functions it might be better to work with the numerator and
denominator as Poly rather than using Expr in the calculations.

I see some checks in the code that are potentially indeterminate if
the ODE has other symbols as parameters. Those need to be handled
carefully.

For example how does this handle a case where the poles are at
"symbolic" values? I mean suppose the denominator factorises as (x -
a)*(x - b) and we don't know whether or not a == b.

> > Also, try splitting the solver into a number of independent functions, like a function to determine poles, checking if a pole is really a pole, etc.
>
> Sure, I'll do that.

Yes, I agree. With separate smaller functions you can test and
document the different parts.

Probably there should be e.g.:

1. A function for extracting b0, b1, b2 from canonical form
2. Another for converting those to normal form
3. etc

Functions working with the poles might be reusable for other methods.
The 2nd_hypergeometric solver also works with the poles and could
potentially share some code.

Oscar

[Naveen Saisreenivas Thota]

>  That looks nice. It would be easier to comment on the details if it was a PR.

I didn't make a PR yet since it is a very basic code and more things are to be changed.

> One comment I have is that if this is all working with rational

> functions it might be better to work with the numerator and
> denominator as Poly rather than using Expr in the calculations.

I didn't use Poly since there was some discussion on Poly being slow during creation.

> I see some checks in the code that are potentially indeterminate if

> the ODE has other symbols as parameters. Those need to be handled
> carefully.

I've changed some of the checks (will update it in the gist), but we will have to test more examples and figure out which checks could potentially fail.

> Functions working with the poles might be reusable for other methods.

> The 2nd_hypergeometric solver also works with the poles and could
> potentially share some code.

Yes, we can merge both of them together.

Now that a basic Riccati Solver is done, it wouldn't take more than a week or two to merge it in the codebase. I'd like to know what all we plan on doing as a GSoC project as that would help me prepare the proposal. Like, is Kovacic's algorithm more of a priority or should Computation of All Rational Solutionsof First-Order Algebraic ODEs be implemented? Ofcourse, if time permits, I would like to implement both.

Naveen

[Oscar Benjamin]

On Mon, 29 Mar 2021 at 08:13, Naveen Saisreenivas Thota
<naveensai...@gmail.com> wrote:
>
> Now that a basic Riccati Solver is done, it wouldn't take more than a week or two to merge it in the codebase. I'd like to know what all we plan on doing as a GSoC project as that would help me prepare the proposal. Like, is Kovacic's algorithm more of a priority or should Computation of All Rational Solutionsof First-Order Algebraic ODEs be implemented? Ofcourse, if time permits, I would like to implement both.

I think you underestimate how much work is involved in really making
the implementation robust and complete. Note that it's much better to
have a well-tested, complete, efficient implementation of a single
algorithm with nicely organised and documented code than it is to have
multiple half-implemented algorithms. As Nijso emphasised earlier the
most important thing first is to establish a systematic test base. We
should get the Kamke examples in and you should verify that this does
find all the rational function solutions for all of the Ricatti ODEs.

If it's possible to implement the other method as well then that's
great but quality is better than quantity here.


Oscar

[Naveen Saisreenivas Thota]

> I think you underestimate how much work is involved in really making

> the implementation robust and complete. Note that it's much better to
> have a well-tested, complete, efficient implementation of a single
> algorithm with nicely organised and documented code than it is to have
> multiple half-implemented algorithms. As Nijso emphasised earlier the
> most important thing first is to establish a systematic test base. We
> should get the Kamke examples in and you should verify that this does
> find all the rational function solutions for all of the Ricatti ODEs.

I was thinking as much, but I wanted to ask just to know your opinion as well. I did test the current code with some examples, but I am yet to test it with all of them. So, from what you say, I am planning to include Rational Riccati Solver and ODE test bank (Kamke and Murphy) as the primary items to be done and leave computation of rational solutions for a general 1st order equation as a bonus? Will this be okay?

Naveen

[Oscar Benjamin]

On Mon, 29 Mar 2021 at 10:42, Naveen Saisreenivas Thota
<naveensai...@gmail.com> wrote:
>
> > I think you underestimate how much work is involved in really making
> > the implementation robust and complete. Note that it's much better to
> > have a well-tested, complete, efficient implementation of a single
> > algorithm with nicely organised and documented code than it is to have
> > multiple half-implemented algorithms. As Nijso emphasised earlier the
> > most important thing first is to establish a systematic test base. We
> > should get the Kamke examples in and you should verify that this does
> > find all the rational function solutions for all of the Ricatti ODEs.
>
> I was thinking as much, but I wanted to ask just to know your opinion as well. I did test the current code with some examples, but I am yet to test it with all of them. So, from what you say, I am planning to include Rational Riccati Solver and ODE test bank (Kamke and Murphy) as the primary items to be done and leave computation of rational solutions for a general 1st order equation as a bonus? Will this be okay?

Yes, I think that sounds good.

Note, as I said in reply to some other queries about GSOC exactly what
you would or wouldn't achieve within the duration of the project is
less important than demonstrating that you are capable of making
significant contributions to SymPy. All tasks can turn out to be
harder or easier than expected so it's hard to estimate in advance
what is possible given a fixed timeframe.

When reviewing GSOC applications (just speaking for myself - I am not
the only reviewer) I am most interested in ensuring that we can get
the best contributors who are capable of making the most valuable
contributions to important parts of SymPy. What you are proposing here
is a significant improvement to an important part of SymPy so the main
points to focus on in your application are:
1) making it clear why this is important and how significant the improvement is
2) demonstrating that you personally understand what needs doing and
are capable of doing the necessary work

Then if your application is successful and it turns out that (based on
the work you have already done) it is not hard to complete some of the
tasks listed then there is no shortage of other things to be done for
ODEs in SymPy. On the other hand if one of the tasks turns out to be
more involved than expected then it is better to limit the scope of
the project later and make sure that the parts that are implemented
are done well.

A general point that I often make to students is that (usually) it is
better to do half a job well than to do the whole job badly. If half a
job is done well then it makes a good starting point for someone in
future to finish that work. If the whole job is done badly it
potentially makes it more difficult for someone else to improve that
work than it would be for them if starting from scratch.

Oscar

[Naveen Saisreenivas Thota]

> When reviewing GSOC applications (just speaking for myself - I am not

> the only reviewer) I am most interested in ensuring that we can get

> the best contributors who are capable of making the most valuable

> contributions to important parts of SymPy. What you are proposing here

> is a significant improvement to an important part of SymPy so the main

> points to focus on in your application are:

> 1) making it clear why this is important and how significant the improvement is

> 2) demonstrating that you personally understand what needs doing and

> are capable of doing the necessary work

Okay, thank you for the advice, Oscar! I'll make the proposal and post it here so that you and others can review it.

Naveen

[Naveen Saisreenivas Thota]

Hi Oscar,

Should we add all the ODEs of Kamke and Murphy or only Riccati ODEs? In either case, how do we plan on parsing the solution from Maple/Mathematica? I could see that there is a Mathematica Parser, but even that seems to be very basic and is not parsing some complex expressions.

Naveen

[nijso.be...@gmail.com]

Hi Naveen, Oscar,

I agree with Oscar that there is still much work to be done, even having 'only' the rational Riccati solver implemented and tested is still a lot of work, after all you have to implement and test all subcases (movable poles, fixed poles, poles at infinity). And the solver should not only work on the solvable ODEs, but should also not get stuck on ODEs that are not solvable using this method. Also sometimes we will find only special solutions and we would need to construct the general solution from them. Since this is a core solver for possibly many other solvers, it is important that it functions very well.

I have put the first 367 kamke odes (the ODEs of first order and first degree) in sympy format in a list, it is here: https://github.com/bigfooted/sympy_ode

I think a modular approach is best and the solver should consist of a number of independent functions that can be re-used elsewhere. For instance:

isRiccati(ode)  : return True if the input ODE is Riccati, False otherwise,

Riccati2Normal(ode) : returns the Riccati ode in normal form,

Poles(ode) : returns the poles, and their order and multiplicities,

etc..

But I think you have made a good start.

Regarding the GSoC application: I would focus a bit more on that now, it's always good to be able to get some feedback before submitting (not sure how that works regarding independent reviewing, though).

Best regards,

Nijso

...

[nijso.be...@gmail.com]

By the way, Naveen, could you create a repository for your work? It will be easier to work with.

Best regards,

Nijso

[nijso.be...@gmail.com]

Hi Naveen,

This ODE comes from Banks, Gundersen, Laine: Meromorphic Solutions of the Riccati Differential Equation (example 4.2 http://www.acadsci.fi/mathematica/Vol06/vol06pp369-398.pdf)

ode = y(x).diff(x) - y(x)**2 - 21/2 + 9/4*x**2

Your current implementation fails to find the solution (-1/(x-1)   -1/x - 1/(x+1) + 3x/2)

This ODE is example 4.6 from the thesis of Thieu:

ode = y(x).diff(x) - ( (-3*x**2+2*x-2)/(x*(x-1)**2) - (6*x**2-x+3)/(x*(x-1))*y(x) - (3*x**2+1)/(x)*y(x)*y(x))

Your implementation finds 3 complex solutions while it should return 1 rational solution.

Best regards,

Nijso

[Naveen Saisreenivas Thota]

Hi Nijso, 

I've made some changes to the code - split it up into functions and added all the Riccati ODEs with Rational Solution from the report. The present code can solve all the ODEs, but not completely in the sense that I've substituted for constants like a, b, c with values in the equation.

I have a doubt in Step 11 of Algorithm 11.

I'm unable to understand how to compute solutions symbolically when the value of "m" is unknown. Is it possible? I can use sympy's assumptions to make sure that m.is_integer gives True, but how do I find polynomial solutions of degree "m" for the auxiliary equation when "m" is unknown?

Naveen

[nijso.be...@gmail.com]

Hi Naveen, 

Great to see that you are making so much progress in such a short time! Regarding your question:

>I'm unable to understand how to compute solutions symbolically when the value of "m" is unknown. Is it possible? I can use sympy's assumptions to make sure that m.is_integer gives True, but how do I find polynomial solutions of degree "m" >for the auxiliary equation when "m" is unknown?

You must have an upper bound for m, or else you have to check all polynomials of degree 1..infinity. The upper bound is computed in step 8, so that means that all d_i and c_i must be known values.  There might be special solutions that can be found by searching for a low order polynomial solution (P=1,2,3,4), and from that you could then compute the general solution. But it is not guaranteed that this will give you a solution.

Can you make a repository for your implementation? Then it will also be easier to collaborate. You can then also put all ODEs in a separate regression test. Since it is at this moment an independent software package, it needs a license, it would be good to choose the same as sympy (new bsd).

The solution should be in Q, so rational, and some solutions are now in C, so complex. This should be looked into. 

You previously mentioned that you would also like to work on either Kovacic or the generalization of the rational Riccati ODE to rational first order ODEs if time permits. This Riccati solver is a base solver for these methods and should therefore be thoroughly tested. That is why I think that before adding more functionality, the focus should be on a good regression test. You could make a time planning, with time reserved for making the regression test, and finding and fixing problems for the Riccati solver. In my experience, successfully going through the entire ODE database takes some time, especially if you care about time efficiency and returning 'nice' answers.

Then, looking at Kovacic' method, you will see that it also depends on the characterization of the poles of f(x) in the ODE y"=f(x)*y, so making a function for pole computation in such a way that it can be used by Kovacic' method will allow code sharing. Kovacic' method can easily give you very large expressions that can take a very long time to compute, so it takes much more time to debug and test this algorithm. That is something that you should take into account.

Best regards,

Nijso

[Naveen Saisreenivas Thota]

Hi Nijso,

> Great to see that you are making so much progress in such a short time!

Thank you! I have to say that you are really helping a lot creating a roadmap that can help the ODE module a lot. Infact, the Riccati solver wasn't even being considered until you talked about it!

> You must have an upper bound for m, or else you have to check all polynomials of degree 1..infinity. The upper bound is computed in step 8, so that means that all d_i and c_i must be >known values.  There might be special solutions that can be found by searching for a low order polynomial solution (P=1,2,3,4), and from that you could then compute the general solution. > But it is not guaranteed that this will give you a solution.

So there is no way to generate a general solution if there are symbolic constants in the equation? For example, if we consider Example 1.24, can we not get the solution as given in the table, i.e. in terms of "a" and "b" without substituting numbers for them like I'm presently doing?

> Can you make a repository for your implementation? Then it will also be easier to collaborate. You can then also put all ODEs in a separate regression test. Since it is at this moment an  >independent software package, it needs a license, it would be good to choose the same as sympy (new bsd).

Sure, I will do it by today or tomorrow.

> The solution should be in Q, so rational, and some solutions are now in C, so complex. This should be looked into. 

To prevent this, should I discard imaginary poles?

> You previously mentioned that you would also like to work on either Kovacic or the generalization of the rational Riccati ODE to rational first order ODEs if time permits. This Riccati solver > is a base solver for these methods and should therefore be thoroughly tested. That is why I think that before adding more functionality, the focus should be on a good regression test. You > could make a time planning, with time reserved for making the regression test, and finding and fixing problems for the Riccati solver. In my experience, successfully going through the  >entire ODE database takes some time, especially if you care about time efficiency and returning 'nice' answers.

Yes, in fact testing is taking more time than writing the code for the solver, so I understand what you mean! So, I would definitely focus on testing the Riccati solver on as many examples as possible from Kamke, Murphy, etc. If time still remains, I will go ahead and implement either Kovacic's algorithm or the algorithm to get all rational solutions of a general first order ODE depending on what feels more necessary for SymPy.

> Then, looking at Kovacic' method, you will see that it also depends on the characterization of the poles of f(x) in the ODE y"=f(x)*y, so making a function for pole computation in such a  >way that it can be used by Kovacic' method will allow code sharing. Kovacic' method can easily give you very large expressions that can take a very long time to compute, so it takes much > more time to debug and test this algorithm. That is something that you should take into account.

I've seen some papers which use Differential Galois theory and come up with an algorithm that essentially solves the same problem, but the solvers are much faster. We can look into all of the available options and then decide.

Naveen

[Naveen Saisreenivas Thota]

Hi Nijso and Oscar,

I just created a GitHub repo for the solver - RationalRiccatiSolver.

Naveen

[Naveen Saisreenivas Thota]

Hi Nijso and Oscar,

I just finished making my proposal. Please go through it and let me know what you feel.

Regards

Naveen

[nijso.be...@gmail.com]

Hi Naveen, 

Great, very detailed. But maybe first do the testing of the solver with the ODE database before merging with dsolve, unless your idea is to use the database for all ODE solvers. That would also considerably shift the focus towards testing of all subsolvers and if this is your plan (you mention something along those lines) then you should place more focus on this part, maybe even change the title of the proposal as this will then be most of your work.

You could also sketch how the Riccati ODE solver would be extended to e.g. second order linear ODEs. Actually, if you simply add a function SecondOrder2Riccati to transform second order linear ODEs and test if it has a rational solution, then you basically have the first class (of 3) of Liouvillian solutions for the second order linear ODE. The second class is also not difficult (as you can see in Bronsteins implementation in axiom, or in my implementation of Smith's maple code in maxima), you just have to construct different Riccati ODEs from the original one.

After implementing Kovacic' method, Riccati ODEs that do not have rational solutions can be transformed into a second order linear ODE and with Kovacic' solver we check if it has a Liouvillian solution.

By the way, the Kovacic method will solve more than 50% of all Kamke ODEs that are second order linear, so pretty good.

Best regards,

Nijso

[Naveen Saisreenivas Thota]

> But maybe first do the testing of the solver with the ODE database before merging with dsolve, unless your idea is to use the database for all ODE solvers. That would also considerably shift the focus towards testing of all subsolvers and if this is your plan (you mention something along those lines) then you should place more focus on this part, maybe even change the title of the proposal as this will then be most of your work.

I was kind of confused whether we need to parse all ODEs or not. It may be helpful in the future, but there are only about 1000 ODEs which are either first or second order. All others are third or higher order. We can maybe skip these. I was thinking of changing the schedule a bit. Two possibilities that crossed my mind are - 

1. I open the PRs for the solver and parser in Phase I. The parser can be merged after appropriate testing. Then, when we parse all the required ODEs, we simultaneously test the solver and make changes as and when required. We can finally merge the solver in the end.

2. I start with only parser and port all the required ODEs in Phase I (and a week of Phase II if required). In Phase II, I start work on the solver and test it thoroughly before merging it.

Which one do you think is better?

> You could also sketch how the Riccati ODE solver would be extended to e.g. second order linear ODEs. Actually, if you simply add a function SecondOrder2Riccati to transform second order linear ODEs and test if it has a rational solution, then you basically have the first class (of 3) of Liouvillian solutions for the second order linear ODE. The second class is also not difficult (as you can see in Bronsteins implementation in axiom, or in my implementation of Smith's maple code in maxima), you just have to construct different Riccati ODEs from the original one.

Sure. Here again, I am confused about whether we plan to implement Smith's paper, or Kovacic's paper. It is mentioned in Smith that the common part in Kovacic's algorithm for all the 3 steps is unified, whereas what we are implementing now would be a separate Step 1 procedure. Could you please let me know about this?

Naveen