[WIP] Rational Riccati solver

[naveensaigit]

This PR introduces a Rational Riccati ODE solver. This code can also be used to solve a subclass of 2nd order ODEs solvable by Kovacic's Algorithm

References to other Issues or PRs

Partial Fix #19183, #21503

Brief description of what is fixed or changed

The algorithm present here is implemented. The examples used to test are here.

Other comments

I will soon add the Kovacic solver.

Release Notes

• solvers
  • ode: Added a new solver for Riccati equation

[sympy-bot]

✅

Hi, I am the SymPy bot (v161). I'm here to help you write a release notes entry. Please read the guide on how to write release notes.

Your release notes are in good order.

Here is what the release notes will look like:

• solvers
  • ode: Added a new solver for Riccati equation (#21459 by @naveensaigit)

This will be added to https://github.com/sympy/sympy/wiki/Release-Notes-for-1.9.

Click here to see the pull request description that was parsed.

This PR introduces a Rational Riccati ODE solver. This code can also be used to solve a subclass of 2nd order ODEs solvable by [Kovacic's Algorithm](https://cs.uwaterloo.ca/research/tr/1984/CS-84-35.pdf)

#### References to other Issues or PRs
Partial Fix #19183, #21503

#### Brief description of what is fixed or changed
The algorithm present [here](https://www3.risc.jku.at/publications/download/risc_5387/PhDThesisThieu.pdf) is implemented. The examples used to test are [here](https://www3.risc.jku.at/publications/download/risc_5197/RISCReport15-19.pdf).

#### Other comments
I will soon add the Kovacic solver.

#### Release Notes
<!-- BEGIN RELEASE NOTES -->
* solvers
  * ode:  Added a new solver for Riccati equation
<!-- END RELEASE NOTES -->

Update

The release notes on the wiki have been updated.

[oscarbenjamin]

I haven't looked through this in detail but it's good to see this coming.

I would put the algorithm in a separate file like ricatti.py. There can be a solver in single.py but it should just be a wrapper around the basic algorithm routines from ricatti.py. The basic routines should each be tested separately from dsolve.

The code needs to give much more explanation about the algorithm. Imagine someone in future who hasn't read the papers on which the algorithm is based trying to fix a simple bug. It's too much for them to go read a PhD thesis to understand what your code is doing. The top of the ricatti.py module should outline the algorithm and should also make it very clear what the reference for the algorithm is. Any assumptions you make and anything that you are not 100% sure about should be clearly documented within the code so that someone in future can understand what you were thinking when you wrote that code (e.g. if the assumption turns out to be incorrect in some cases).

[sympy-bot]

🟠

Hi, I am the SymPy bot (v161). I've noticed that some of your commits add or delete files. Since this is sometimes done unintentionally, I wanted to alert you about it.

This is an experimental feature of SymPy Bot. If you have any feedback on it, please comment at sympy/sympy-bot#75.

The following commits add new files:

• 382ea94:
  • sympy/solvers/ode/riccati.py

If these files were added/deleted on purpose, you can ignore this message.

[oscarbenjamin]

Is the current structure okay?

The basic structure looks okay.

I'm going to hold off on detailed comments for now because there is one big question that needs to be resolved which is whether this algorithm should be formulated in terms of Poly rather than Expr.

There are lots of places which are using cancel, together, as_numer_denom etc which suggests that you should really be using a pair of Poly rather than an Expr. Every time you call cancel the whole expression is converted into a pair of Poly and then the Poly.cancel method is used and then the results are converted back to Expr. This is very inefficient compared to just working with Poly and using the Poly.cancel method.

Poly expressions are always in a canonical form which simplifies a lot of things and makes algebra more robust, predictable and efficient for an algorithm such as this.

Currently this uses series to compute Laurent series but the series function is slow because it is intended for more general cases. I'm pretty sure that the relevant Laurent coefficients for a rational function can be computed more efficiently using only polynomial operations.

If this can be written in terms of Poly and if it is better to do so then it is better to make that decision now as it will only be harder later. This is especially important if this is just the first of many related algorithms that should all share code and use the same kind of internal representation.

There are certain limitations of Poly such as the fact that symbolic powers of the independent variable would not be possible e.g. x/(x + 1)**n where n is a symbol. If this algorithm also has the same kind of limitations then I think it would be better to use Poly.

We can discuss this when we meet.

[naveensaigit]

I think using Poly would be best since finding the degree is not possible for polynomials with symbolic powers (which means the algorithm won't work for equations with symbolic powers in x).

As for Laurent series, I'll have to see how to get the coefficients in an efficient way without using the inbuilt function. I think that there is something mentioned in Smith's implementation of Kovavic's algorithm about finding the coefficients.

[asmeurer]

If we use Poly, these helper functions in risch will be useful.

sympy/sympy/integrals/risch.py

Lines 811 to 886 in 31d469a

	def frac_in(f, t, *, cancel=False, **kwargs):
	"""
	Returns the tuple (fa, fd), where fa and fd are Polys in t.
	Explanation
	===========
	This is a common idiom in the Risch Algorithm functions, so we abstract
	it out here. ``f`` should be a basic expression, a Poly, or a tuple (fa, fd),
	where fa and fd are either basic expressions or Polys, and f == fa/fd.
	**kwargs are applied to Poly.
	"""
	if type(f) is tuple:
	fa, fd = f
	f = fa.as_expr()/fd.as_expr()
	fa, fd = f.as_expr().as_numer_denom()
	fa, fd = fa.as_poly(t, **kwargs), fd.as_poly(t, **kwargs)
	if cancel:
	fa, fd = fa.cancel(fd, include=True)
	if fa is None or fd is None:
	raise ValueError("Could not turn %s into a fraction in %s." % (f, t))
	return (fa, fd)
	def as_poly_1t(p, t, z):
	"""
	(Hackish) way to convert an element ``p`` of K[t, 1/t] to K[t, z].
	In other words, ``z == 1/t`` will be a dummy variable that Poly can handle
	better.
	See issue 5131.
	Examples
	========
	>>> from sympy import random_poly
	>>> from sympy.integrals.risch import as_poly_1t
	>>> from sympy.abc import x, z
	>>> p1 = random_poly(x, 10, -10, 10)
	>>> p2 = random_poly(x, 10, -10, 10)
	>>> p = p1 + p2.subs(x, 1/x)
	>>> as_poly_1t(p, x, z).as_expr().subs(z, 1/x) == p
	True
	"""
	# TODO: Use this on the final result. That way, we can avoid answers like
	# (...)*exp(-x).
	pa, pd = frac_in(p, t, cancel=True)
	if not pd.is_monomial:
	# XXX: Is there a better Poly exception that we could raise here?
	# Either way, if you see this (from the Risch Algorithm) it indicates
	# a bug.
	raise PolynomialError("%s is not an element of K[%s, 1/%s]." % (p, t, t))
	d = pd.degree(t)
	one_t_part = pa.slice(0, d + 1)
	r = pd.degree() - pa.degree()
	t_part = pa - one_t_part
	try:
	t_part = t_part.to_field().exquo(pd)
	except DomainError as e:
	# issue 4950
	raise NotImplementedError(e)
	# Compute the negative degree parts.
	one_t_part = Poly.from_list(reversed(one_t_part.rep.rep), *one_t_part.gens,
	domain=one_t_part.domain)
	if 0 < r < oo:
	one_t_part *= Poly(t**r, t)
	one_t_part = one_t_part.replace(t, z) # z will be 1/t
	if pd.nth(d):
	one_t_part *= Poly(1/pd.nth(d), z, expand=False)
	ans = t_part.as_poly(t, z, expand=False) + one_t_part.as_poly(t, z,
	expand=False)
	return ans

[oscarbenjamin]

To make the substitution x -> 1/x in a rational function represented as two Polys you can use num.transform:

In [11]: f = (x**2 - 1)/(x**2 + 1)

In [12]: f
Out[12]: 
 2    
x  - 1
──────
 2    
x  + 1

In [13]: num, den = [Poly(e, x) for e in f.as_numer_denom()]

In [14]: num
Out[14]: Poly(x**2 - 1, x, domain='ZZ')

In [15]: den
Out[15]: Poly(x**2 + 1, x, domain='ZZ')

In [18]: num.transform(Poly(1, x), Poly(x, x))
Out[18]: Poly(-x**2 + 1, x, domain='ZZ')

In [19]: den.transform(Poly(1, x), Poly(x, x))
Out[19]: Poly(x**2 + 1, x, domain='ZZ')

In [20]: num.transform(Poly(1, x), Poly(x, x)) / den.transform(Poly(1, x), Poly(x, x))
Out[20]: 
     2
1 - x 
──────
 2    
x  + 1

If the degree of the numerator and denominator are different then some factors of x need to be added somewhere.

[naveensaigit]

If the degree of the numerator and denominator are different then some factors of x need to be added somewhere.

I think this should work.

In [3]: f = (x**2 - 1)/x**3

In [4]: f.subs(x, 1/x)
Out[4]: 
 3 ⎛     1 ⎞
x ⋅⎜-1 + ──⎟
   ⎜      2⎟
   ⎝     x ⎠

In [5]: num, den = [Poly(e, x) for e in f.as_numer_denom()]

In [6]: num.transform(Poly(1, x), Poly(x, x)) / den.transform(Poly(1, x), Poly(x, x))
Out[6]: 
     2
1 - x

In [7]: (num.transform(Poly(1, x), Poly(x, x)) / den.transform(Poly(1, x), Poly(x, x)))*x**(den.degree(x) - num.degree(x))
Out[7]: 
  ⎛     2⎞
x⋅⎝1 - x ⎠

[naveensaigit]

We were discussing yesterday that it's not possible to use series with Poly. Should we use Expr for this, or should this be implemented?

[oscarbenjamin]

We were discussing yesterday that it's not possible to use series with Poly. Should we use Expr for this, or should this be implemented?

For now you could just just use Expr for that part.

I think it should be straight-forward to make a routine to compute the series for rational functions though. For example the coefficients in the Laurent part of the series can be calculated like:

In [12]: f = (1 + x)/(4 - x**2)                                                                      

In [13]: r1, r2 = roots(f.as_numer_denom()[1])                                                       

In [14]: f.series(x, r1)                                                                             
Out[14]: 
                            2            3            4            5                            
      1       9    3⋅(x + 2)    3⋅(x + 2)    3⋅(x + 2)    3⋅(x + 2)    3⋅x    ⎛       6        ⎞
- ───────── + ── + ────────── + ────────── + ────────── + ────────── + ─── + O⎝(x + 2) ; x → -2⎠
  4⋅(x + 2)   32      256          1024         4096        16384       64                      

In [15]: cancel(f*(x - r1)).subs(x, r1)                                                              
Out[15]: -1/4

In general for the series you have a recursion relation:
https://en.wikipedia.org/wiki/Rational_function#Taylor_series

For the f here we can rearrange by substituting x -> x + r1 multiplying out the x from the denominator:

In [35]: cancel(f.subs(x, x+r1)*x)                                                                   
Out[35]: 
1 - x
─────
x - 4

Now we want the Taylor series of this around x=0. We write this as

In [39]: a = IndexedBase('a')                                                                        

In [40]: Eq(1 - x, (x - 4)*Sum(a[m]*x**m, (m, 0, oo)))                                               
Out[40]: 
                  ∞          
                 ___         
                 ╲           
                  ╲    m     
1 - x = (x - 4)⋅  ╱   x ⋅a[m]
                 ╱           
                 ‾‾‾         
                m = 0 

Multiply out and rearrange the indices of the sums:

1 - x = -4*a0 + (a0 - 4*a1)*x + Sum((a[i-1] - 4*a[i])*x**i, (x, 2, oo)

Then equate coefficients of powers of x and you get the recurrence problem:

a0 = -1/4
a1 = 3/16
a[m] = a[m-1]/4

The part of this that would be potentially tricky though is needing to be able to compute in an extension field that includes the roots of the polynomials when the roots are not rational e.g.

In [43]: f = (1 + x)/(2 - x**2)                                                                      

In [44]: r1, r2 = roots(f.as_numer_denom()[1])                                                       

In [45]: f.series(x, r1)                                                                             
Out[45]: 
  1   √2                                                                                           
- ─ + ──                                                                                           
  2   4    √2   1   ⎛√2   1 ⎞            ⎛1    √2⎞         2   ⎛ √2    1 ⎞         3   ⎛ 1     √2⎞ 
──────── + ── + ─ + ⎜── + ──⎟⋅(x + √2) + ⎜── + ──⎟⋅(x + √2)  + ⎜─── + ───⎟⋅(x + √2)  + ⎜─── + ───⎟⋅
 x + √2    8    8   ⎝32   16⎠            ⎝64   64⎠             ⎝256   128⎠             ⎝512   512⎠ 

                                                           
                                                           
        4   ⎛ √2     1  ⎞         5    ⎛        6         ⎞
(x + √2)  + ⎜──── + ────⎟⋅(x + √2)  + O⎝(x + √2) ; x → -√2⎠
            ⎝2048   1024⎠                                  

To handle something like that you want the domain to be QQ<sqrt(2)> rather than EX ideally:

In [49]: f = (1 + x)/(2 - x**2)                                                                      

In [50]: num, den = [Poly(e, x) for e in f.as_numer_denom()]                                         

In [51]: r1, r2 = roots(den)                                                                         

In [52]: r1                                                                                          
Out[52]: -√2

In [55]: Poly(r1, x)                                                                                 
Out[55]: Poly(-sqrt(2), x, domain='EX')

In [56]: Poly(r1, x, extension=True)                                                                 
Out[56]: Poly(-sqrt(2), x, domain='QQ<sqrt(2)>')

In [57]: r1p = Poly(r1, x, extension=True)                                                           

In [58]: num.transform(x + r1p, 1)                                                                   
Out[58]: Poly(x + 1 - sqrt(2), x, domain='QQ<sqrt(2)>')

In [59]: den.transform(x + r1p, 1)                                                                   
Out[59]: Poly(-x**2 + 2*sqrt(2)*x, x, domain='QQ<sqrt(2)>')

In [60]: cancel(f.subs(x, x + r1)*x)                                                                 
Out[60]: 
-x - 1 + √2
───────────
  x - 2⋅√2 

In general we might want to use FiniteExtension as a domain for these calculations if the roots have symbols in.

[oscarbenjamin]

One of the advantages of having separate functions for the different steps of the algorithm is that we can make separate tests for them. There should be direct tests for the functions in ricatti.py separately from the dsolve tests.

[naveensaigit]

One of the advantages of having separate functions for the different steps of the algorithm is that we can make separate tests for them. There should be direct tests for the functions in ricatti.py separately from the dsolve tests.

Should this be in a separate file, like test_riccati.py?

[oscarbenjamin]

Should this be in a separate file, like test_riccati.py?

Yes

[naveensaigit]

I've made most of the changes that you suggested. I will now work on creating the RationalFunction class.

[naveensaigit]

@oscarbenjamin Anything else that is left to do here? Right now, only the sphinx build is failing. I think its because I should add a newline before the docstring ends. I will change that along with any others if required.

[github-actions]

Benchmark results from GitHub Actions

Lower numbers are good, higher numbers are bad. A ratio less than 1
means a speed up and greater than 1 means a slowdown. Green lines
beginning with + are slowdowns (the PR is slower then master or
master is slower than the previous release). Red lines beginning
with - are speedups.

Significantly changed benchmark results (PR vs master)

Significantly changed benchmark results (master vs previous release)

       before           after         ratio
     [ed9a550f]       [fb2d2999]
     <sympy-1.8^0>                 
+         954±3ms          1.62±0s     1.70  dsolve.TimeDsolve01.time_dsolve
-      7.83±0.01s          3.73±0s     0.48  integrate.TimeIntegrationRisch02.time_doit_risch(100)
-      13.3±0.8ms      8.38±0.02ms     0.63  logic.LogicSuite.time_load_file
+        69.5±4μs      2.27±0.01ms    32.59  matrices.TimeDiagonalEigenvals.time_eigenvals
-     4.56±0.01ms         2.46±0ms     0.54  solve.TimeRationalSystem.time_linsolve(10)
-         922±2μs          574±2μs     0.62  solve.TimeRationalSystem.time_linsolve(5)
-        1.07±0ms          704±5μs     0.66  solve.TimeSparseSystem.time_linsolve_eqs(10)
-        2.02±0ms         1.28±0ms     0.63  solve.TimeSparseSystem.time_linsolve_eqs(20)
-     2.99±0.01ms         1.87±0ms     0.63  solve.TimeSparseSystem.time_linsolve_eqs(30)

Full benchmark results can be found as artifacts in GitHub Actions
(click on checks at the top of the PR).

[oscarbenjamin]

Looks good. Let's get this in!

[oscarbenjamin]

Actually it would be good to expand the release note a bit. Can you expand the release note a little to explain what the significance of this solver is:
https://github.com/sympy/sympy/wiki/Release-Notes-for-1.9

[oscarbenjamin]

It would be good to add some benchmark cases for the algorithm to the benchmarks repo:
https://github.com/sympy/sympy_benchmarks

[naveensaigit]

Actually it would be good to expand the release note a bit. Can you expand the release note a little to explain what the significance of this solver is:
https://github.com/sympy/sympy/wiki/Release-Notes-for-1.9

I've added some notes for it. Let me know if something else must also be mentioned.

It would be good to add some benchmark cases for the algorithm to the benchmarks repo:
https://github.com/sympy/sympy_benchmarks

Sure. Should I put it in benchmarks/dsolve.py ?

[oscarbenjamin]

Sure. Should I put it in benchmarks/dsolve.py ?

Yes, I think so. It would be good to test both how long dsolve takes but also how long the solver takes without dsolve.

[elieraphael]

Hi,
On 07/21/2021, Riccati's equation (#19183) was closed via #21459. However, solving Riccati's equation y’ = -y**2 + x*y + 1 using sympy still does not work (sympy does not find the closed-form solution y(x) = E^(-x^2/2)/ (Sqrt(Pi/2) Erf(x/Sqrt(2)) + c)) + x (with Erf the error function, and c an arbitrary constant). Thanks for your help!
Elie.

[bigfooted]

The Riccati solver can only find rational solutions. To solve this ODE, the kovacic method or something similar needs to be implemented as well.

[elieraphael]

Is such an implementation planned in the near future? That would be really useful.

[elieraphael]

Any progress in that direction? Thanks a lot. Sympy is great!

[oscarbenjamin]

There hasn't been any progress. I think we need to spell out somewhere what should be done so that it can be maybe another GSOC project.

[elieraphael]

Thanks a lot!

…

Le 30 juin 2022 à 14:51, Oscar Benjamin ***@***.***> a écrit : There hasn't been any progress. I think we need to spell out somewhere what should be done so that it can be maybe another GSOC project. — Reply to this email directly, view it on GitHub <#21459 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/APKF3JFGQAQ3ROK3OQO5L23VRWJ4FANCNFSM44WHY3ZQ>. You are receiving this because you commented.