Introduction

Below is a roadmap for developing a robust and general first order ODE solver for SymPy. I will focus for now on explicit first order ODEs of the form . A first order ODE solver is important as it can function as a basis for other ODE solvers and even PDE solvers. Many ODEs can be reduced to the problem of solving a first order ODE, and even PDEs can be transformed into first order second order ODEs (method of characteristics, separation of variables). It is therefore important that components of the solver can be reused by other solvers and have flexible calling behaviour. Furthermore, one should keep in mind that it can be notoriously difficult to find compact (in terms of expression length), explicit analytic solutions of any ODE. Sometimes one is not even interested in such an expression but in, say, the integrating factor, a first integral, the invariant algebraic curves or a symmetry generator of an ODE. A solver that returns any or all of these 'solutions' is more valuable than a solver that just returns a (sometimes very ugly) expression, or nothing at all. So what is a viable approach?

Fundamental solvers for first order ODEs (1-ODEs)

Some ODE solver methods are more fundamental than others, in the sense that they can be considered specializations of a more general method. Some ODEs do not need any solving, as they can be integrated immediately, like separable ODEs of the form . However, it can be difficult sometimes to see if an ODE is separable, for instance the ODE

is separable. One could try to manipulate the expression on the right-hand side until it has the form but structured methods of finding out when an ODE (or a function F(x,y)) is separable or not exist, see e.g. the paper of J.A. Cid A simple method to find out if an ordinary differential equation is separable

As for the separable ODE, also exact 1-ODEs and linear 1-ODEs can be solved almost directly by integration. Because in the end, the solution of the ODE is obtained by integration, the performance of an ODE solver depends on the performance of the integration routines.

Many well-known 1-ODEs have their own special solution methods, and a natural way to extend capabilities of the ODE solver is by simply adding more and more specialized solution methods for specific ODEs. But woudldn't it be better if there was one general method that can solve all first order ODEs, or at least very large classes of ODEs, so we do not need all these dedicated solvers? So which solution methods exist that can solve the largest class of 1-ODEs?

Symmetry methods

A very general method to tackle differential equations are symmetry methods. A symmetry is basically a coordinate transformation that leaves the equation unchanged (an invariant transformation). For a quick introduction, see Peter E. Hydon, Symmetry Methods for Differential Equations, or H. Stephani, Differential Equations. The high level idea is to find a Lie symmetry (after Sophus Lie, pronounced Lee), and with this Lie symmetry construct an integrating factor. The integrating factor transforms the ODE into an exact ODE, which we know how to solve. So when you find a symmetry of a first order ODE, you solve the ODE. How do we find symmetries? For first order ODEs, we have to find a solution of the linearized symmetry condition:

This is a linear PDE in the two variables . Since most 1-ODEs have very simple symmetries, we could try a couple of symmetries to see if we can find one that works. The simplest symmetry is . If we fill that into our symmetry generator, we get:

The most general expression for a symmetry that we can still manage to solve is and a solution method is given by Kolokolnikov and Cheb-Terrab, First order ODEs, Symmetries and Linear Transformations.

To demonstrate the power of this method, a large collection of ODEs from the book of E. Kamke, Differentialgleichungen was tested and this solver alone could solve 75% of the first order ODEs, including linear, separable, and exact ODEs. A robust symmetry method is the most general and the single most important algorithm for solving 1-ODEs. The keyword here is of course robust, as intermediate expressions during the computation can grow and a successful solution depends much on the quality of SymPy to be able to simplify such expressions.

What is available in SymPy?

Sympy has a number of symmetry methods implemented already, but they are almost all specializations of the method of Kolokolnikov. Of course, Specializations exist because they are simpler to implement and usually faster to run. However, these methods could all benefit from a shared solver hierarchy to reduce code duplication. Furthermore, it does not make much sense to expose the user to the different choices available here: The solver method can simply be: "Lie Solver" and internally either a list of simple Lie solvers or the one size fits all method of Kolokolnikov can be called.

Of course, it is unnecessary to call such a powerful method when most of the input equations from users are linear or separable ODEs. But instead of calling the ODE solver and going through a list of solution methods for ODEs, continuing to the next one when the current solver fails, one could go through a list of solution methods for symmetries of the ODE. This has three advantages. Firstly, the symmetry method is independent on pattern matching, making the method very robust in the sense that no equations that could be solved with the method are missed (assuming that all other functionalities like expression simplifiers work satisfyingly). A second advantage is that much larger classes of ODEs can be solved using symmetries, making its range of applicability much larger. And finally, this method can do more than just try to find a solution to the ODE: the symmetry generator can also return the solvable ODE class for a given symmetry, or we can compute an integrating factor from the symmetry, giving us the transformation to an exact ODE. All this additional information can give valuable insights into the behavior of the ODE, more than just returning (or not) the solution.

Riccati ODE

The Riccati ODE is given by

This ODE actually has symmetries of the form given in the paper of Kolokolnikov, but to find them we must solve a Riccati ODE. So Riccati ODEs need a special treatment. Kolokolnikov mentions a treatment for a couple of subclasses of Riccati ODEs, but a very important class is the rational Riccati ODE that has rational solutions. Many other rational ODEs can be solved by solving a rational Riccati ODE, see Grasegger, Vo, Winkler, Computation of all rational solutions of first order algebraic ODEs Also, second order Linear ODEs with Liouvillian solutions can be solved by solving a Riccati ODE using Kovacic' algorithm. A discussion on how to solve rational Riccati ODEs can be found in Yuan, A Note on the Riccati Differential Equation and in Rational and algebraic solutions of first-order algebraic ODEs. Since Riccati equations are important in many applications, like the analysis of algebraic curves (parametrization!), this functionality is an important addition to SymPy's ODE solvers.

Nonliouvillian solutions

A number of 1-ODEs have Nonliouvillian solutions, and their solution can be written in terms of e.g. Bessel functions or hypergeometric functions. We have seen previously that a symmetry method is basically the search for a transformation that makes the ODE exact. Sometimes it is easier to search for a transformation into another equation that has a known solution. We could for instance transform the 1-ODE into a linear second order ODE and see if it is actually a Bessel ODE in disguise. For this we would then need a solver for second order ODEs that can deal with Nonliouvillian solutions, like Chan, Cheb-Terrab, Non-Liouvillian solution for second order linear ODEs. The success of this implementation relies heavily on the capabilities of the SymPy to deal with such special functions.

higher degree ODEs

A naive strategy for solving first order ODEs of higher degree of the form f(y',y,x) is to solve the ODE for the independent variable y' and then solve each of the ODEs individually. However, there are many cases where this strategy either leads to unattractive solutions or to no solution at all when the solver fails to find explicit solutions in terms of y'.

How to continue

Here are some first steps:

A number of methods to find Lie symmetries and solve 1-ODEs already exist in SymPy. However, they can all be brought into the same symmetry solver class, thus greatly reducing code duplication. A clean-up and validation for solver performance (does it find and solve the solvable ODEs?) and robustness (will it not crash or get stuck during computations?) of the symmetry methods should be performed.

Next to that, a solver for Rational Riccati ODEs should be implemented. The implementation should be modular, with separate routines available to determine the genus of the ODE, the movable and fixed poles, and the poles at infinity. These functions can then also be used separately and called by other solvers that have a need for them.

Finally, we need the ODEs from Kamke and Murphy as the standard ODEs for solver performance testing.

With these first steps, one can then proceed to implement solvers for general rational first order ODEs and for second order linear ODEs. With Riccati's solver available, one could also start with a module to construct the Janet base of second order linear ODEs, see e.g. F. Schwarz, Algorithmic Lie theory for solving ordinary differential equations. With such techniques, second order ODEs can be classified and then solved (or not) based on their Janet base.