There are many ways to try to solve differential equations. Some of them belong to computer algebra. Some do not. This text aims at comparing these methods and presenting the contributions of the computer algebra team of the LIFL. Here is an example of an ordinary differential equation with an initial condition. The unknown function is x(t).

x' = x(3 - x),
x(0) = 1.

Numerical integration

The numerical integration of the equation with its initial condition consists in approximating the graph of its solution by a finite number of points. Numerical integration is not considered as a branch of computer algebra. Let us however see how to do it with the MAPLE computer algebra software. One enters the differential equation

> equadiff := diff(x(t),t) - x(t)*(3-x(t));
                                /d      \ 
                    equadiff := |-- x(t)| - x(t) (3 - x(t)) 
                                \dt     /

One solves it, providing its initial condition and declaring that a numerical integration is desired.

> solution := dsolve ({equadiff, x(0)=1}, x(t), numeric);
                   solution := proc(x_rkf45)  ...  end proc 

The solution, in the MAPLE sense, is a function which permits to evaluate (at least approximatively) the solution of the differential equation for many different values of the time.

> solution (0);
                              [t = 0., x(t) = 1.]

> solution (0.5);
                     [t = 0.5, x(t) = 2.07431460567341386]

The simplest of the numerical integration methods is Euler's method. In general, the more sophisticated Runge-Kutta schemes are used. Designing these schemes implies to solve some systems of equations. In practice, this can only be achieved by means of non commutative algebra methods to which |Michel Petitot contributed.

Closed form integration

Integrating a differential equation under closed form consists in computing its solution as a formula. This is the type of integration usually taught in calculus lectures. Closed form integration is a branch of computer algebra. The computer algebra team has only indirectly contributed to this topic by developing methods able to determine whether two differential equations are equivalent up to some change of coordinates (see below). In general, closed form integration is unfortunately not possible: solutions of differential equations do not need to be writable as finite formulas involving the well known functions: sine, exponentials, ... Over the example, close form integration is however possible, the following commands show.

One solves the differential equation under closed form. Since the initial condition is not precised, one gets a formula depending on an arbitrary constant, denoted _C1 by MAPLE.

> dsolve (equadiff, x(t));
                                          3
                          x(t) = -------------------
                                 1 + 3 exp(-3 t) _C1

Qualitative analysis

The qualitative analysis of a differential equation consists in determining the shape of its solutions, without solving the equation. Here is the sort of result one could get by qualitative analysis, over the example.

There are two constant solutions: the functions x(t) = 0 and x(t) = 3. If the initial condition is negative then the curve is decreasing and tends towards minus infinity (as the time tends to plus infinity). If the initial condition is positive and smaller than 3, the curve is increasing and tends towards the limit value 3. If the initial condition is greater than 3, the curve is decreasing and tends towards the limit value 3. The constant solution x(t) = 0 is unstable. The constant solution x(t) = 3 is stable.

The qualitative analysis of an equation, and more generally of a system of differential equations can be performed by means of numerical or symbolic methods. The MAPLE package RegularChains, designed by some members of the team (Marc Moreno Maza, François Lemaire), is a tool which may help performing the qualitative analysis of a polynomial differential system. This package is dedicated to polynomial systems solving.

Differential elimination

Differential elimination is a computer algebra method our team contributed to. To understand its usefulness, it is necessary to consider a system of differential equations. Here is an example.

x' = 0.7 y + sin (2.5 z),
y' = 1.4 x + cos (2.5 z),
1 = x^2 + y^2.

Even a reader not familiar with this type of system, which mixes differential and non differential equations, may understand that numerically integrating this system is problematic. Let us try to integrate it using Euler's method. There are three unknown functions x(t), y(t) and z(t). Assume that three initial conditions x(0), y(0) and z(0) as well as some step size h are fixed. Replacing x, y and z by initial conditions in the two first equations, one can easily compute the values of the derivatives x'(0) and y'(0) and thereby some approximations of x(h) and y(h). However, one does see any simple way to compute z'(0) and the approximation of z(h). It seems that some differential equation of the following form is missing.

z' = something.

The point is that this equation is not missing: it is hidden in the system. Differential equation permits to reveal such hidden equations. Let us see how to do it using the MAPLE DifferentialAlgebra package, written by François Boulier, with the help of Edgardo Cheb-Terrab. First we load the package.

> with (DifferentialAlgebra):

One assigns the differential system to the variable syst. Since differential elimination methods only apply to polynomial systems and since the input system is not polynomial, one encodes the sine and the cosine by introducing two extra variables and three extra equations: s(t) encodes the sine, c(t) the cosine. The so obtained system is polynomial. It is equivalent to the input system.

> syst := [diff(x(t),t) - 7/10*y(t) - s(t),
>              diff(y(t),t) - 14/10*x(t) - c(t),
>              x(t)^2 + y(t)^2 - 1,
>              diff(s(t),t) - 25/10*diff(z(t),t)*c(t),
>              diff(c(t),t) + 25/10*diff(z(t),t)*s(t),
>              s(t)^2 + c(t)^2 - 1]:

One assigns to the variable R the context in which the elimination must be performed: there is only one derivation (w.r.t. to the time) and the ranking is precised: without entering details, the fact that the unknown function z stands on the rightmost place in the list of the unknown function, indicates to the software that we are looking for a differential equation of the form

z' = something

> R := DifferentialRing (derivations = [t], blocks = [[s,c,x,y,z]]):

The differential elimination algorithm is applied to syst and R. The result is stored in the variable ideal.

> ideal := RosenfeldGroebner (syst, R):

The computation being complete, one may ask for the list of the computed equations. The desired equation is a bit large. It stands on the second place of the list.

> Equations (ideal [1], solved);
 d                          d                                 4
[-- y(t) = c(t) + 7/5 x(t), -- z(t) = - (-13230 c(t) x(t) y(t)
 dt                         dt

                           2                              6             4
     + 14700 c(t) x(t) y(t)  - 3940 c(t) x(t) + 12348 y(t)  - 25809 y(t)

                 2           /
     + 16961 y(t)  - 3500)  /  (
                           /

              6             4             2
    11025 y(t)  - 22050 y(t)  + 13525 y(t)  - 2500),

                                         3
             -10 c(t) x(t) y(t) + 21 y(t)  - 21 y(t)
    s(t) = - ---------------------------------------,
                                 2
                          10 y(t)  - 10

        2                       2   441     4   541     2
    c(t)  = -21/5 c(t) x(t) y(t)  + --- y(t)  - --- y(t)  + 1,
                                    100         100

        2        2
    x(t)  = -y(t)  + 1]

Equivalence problems

Two ordinary differential equations (ODE) are said to be equivalent if they are equal up to some change of coordinates. An equivalence problem solver has been developed in our team by Michel Petitot and Raouf Dridi. This solver relies on preliminary works of Sylvain Neut (MAPLE Cartan package) and, more generally, on Élie Cartan's equivalence method. Here is its principle:

  • Precalculate a table of known ordinary differential equations (Van Der Pol equation, Painlevé equations, Airy equation ... and even all the ones in Kamke's book).
  • Given a differential equation to be solved, the solver searches an ODE in the list and a change of coordinates which performs the equivalence.

The difficult problem consists in computing the change of coordinates. In the general setting, it is itself solution of a partial differential equations system. Our solver focuses on changes of coordinates which do not imply any integration. Roughly speaking, the idea consists in restricting the authorized group of transformations in such a way that the symmetry group of the equation gets discrete.

Example of Van Der Pol equation

Let us consider the following second order ordinary differential equation:

 2
d
--- y(x) =
  2
dx

     /d      \2   /d      \        2          4 /d      \          2
    -|-- y(x)|  + |-- y(x)| epsilon  - epsilon  |-- y(x)| exp(y(x))  - 1
     \dx     /    \dx     /                     \dx     /

For this equation, the classical dsolve solver of MAPLE returns a complicated first order ODE and an unsolved quadrature. Our solver finds out that this equation is equivalent to Van Der Pol's

            / 2      \
            |d       |                 2      /d      \
            |--- y(x)| - epsilon (-y(x)  + 1) |-- y(x)| + y(x) = 0
            |  2     |                        \dx     /
            \dx      /

with respect to the change of coordinates (which preserves periodic solutions):

                                                         2
           (x, y, epsilon) -> (x, exp(y) epsilon, epsilon )