iconEuler Home

Differential Equations

We want to look at some differential equations and the numerical methods of Euler to solve these equations.

Reference for Differential Equations

To understand this tutorial, you should have an elementary knowledge about the programming language of Euler.

First, there is the Runge-Kutta method.

For an example, we take the initial value problem

11 - Differential Equations

We want to compute and plot all values of the solution between -3 and 3. So we have to split the region in two, going back and forth from 0.

We first plot the vector field using the function "vectorfield". The functions needs an expression in x and y.


11 - Differential Equations

Then we compute and plot the right side of the solution. The Runge method computes values of the solution at given points. We need to supply a function or an expression, and an initial value.

The Runge method takes a vector of x-values and an initial value at x[1]. It returns the values at all points x[i] in a vector. We can plot the solution easily with these two vectors.

>t=0:0.01:pi; s=runge("x*y*sin(y)",t,1);  ...

11 - Differential Equations

Add the solution from 0 to -3 to the plot, and the initial value.

>t=0:-0.01:-pi; s=runge("x*y*sin(y)",t,1);  ...
   plot2d(t,s,color=red,>add);  ...
   plot2d(0,1,color=red,>points,>add); ...
   textbox(latex("y'=x y \sin(y)",factor=1.5),0.4,0.2):

11 - Differential Equations

It is not necessary to return all the values in the function "runge". But then we have to tell the routine to use lots of intermediate points. In the following example, we get only two values with 500 lost intermediate values.

[1,  3.141592653589793]

Other Solvers

To demonstrate other solvers we solve

11 - Differential Equations

with the solution

11 - Differential Equations

For critical problems, there is an adaptive version of the Runge method. It does not need intermediate steps.

We test the method by counting the evaluations of f(x,y). To do this, we use a function instead of an expression. Inside the function we increase a global count.

>function f(x,y) ...
 global count;
 return -x/(1+y);

Now, let us test the adaptive Runge method.

>longestformat;  ...
   count=0; adaptiverunge("f",[0,1.99],1)[-1], count,

This is very good for 1000 evaluations of f. The correct value is the following.

>function ysol(x) := sqrt(4-x^2)-1; ...

The function has a singularity at x=2.


11 - Differential Equations

For comparison, we need 20000 function evaluations, and don't get the same accuracy with the Runge method.

>count=0; runge("f",[0,1.99],1,5000)[-1], count,

Euler has also the LSODA method, which is the default for the ode function. This is an adaptive method, and it gets very good results with very few evaluations.

>count=0; ode("f",[0,1.99],1)[-1], count,

To increase the accuracy, we can add an epsilon and intermediate points. The LSODA algorithm interpolates the intermediate points with no evaluations of f. If epsilon is too small, it will fail. But 10^-13 seems to be okay for this problem.

>count=0; ode("f",linspace(0,1.99,10000),1,eps=1e-13)[-1], count,

Guaranteed Inclusions

Euler can also get a guaranteed inclusion for the solution of this differential equations.

The function mxmidgl calls Maxima to get a Taylor expansion of the equation of very high degree. It then uses the remainder formula to get an inclusion of the solution.


The correct solution is indeed in this interval.


The function can also be used for integration by solving the differential equation

11 - Differential Equations

The function mxmiint uses this method.


Solving with Maxima

In many cases, Maxima delivers exact solutions. Let us try the initial value problem

11 - Differential Equations

The standard function of Maxima is ode2. To get the solution for an initial value, we use ic1.

>&ode2(y+x*'diff(y,x)=log(x)+1,y,x), sol &= ic1(%,x=1,y=1)
                              x log(x) + %c
                          y = -------------

                               x log(x) + 1
                           y = ------------

If you want to make a function of this solution, use the following trick.

>function ysol(x) &= y with sol
                             x log(x) + 1

We can now plot this function.


11 - Differential Equations

The numerical solution is the same.

>res := ode("(log(x)+1-y)/x",[1,2],1); ...
 longestformat; res[2], ysol(2)

We can try to solve the equation

11 - Differential Equations

in Maxima, but it does not really work.

>&x*y*sin(y)='diff(y,x); &ode2(%,y,x)
                       /                2
                       [    1          x
                       I -------- dy = -- + %c
                       ] y sin(y)      2

Looking at the integral, we see that we need to use a numerical approximation for the integral on the left side, and a further numerical approximation to solve I(y)=x^2/2 for y. We use the fast and stable Gauss method for the integral.

>function fi (y,x) := gauss("1/(x*sin(x))",1,y,20)-x

And the super-stable bisection method for the solution.

>function fy (x) := bisect("fi",1,pi;x^2/2)

Yet we can find the solution only in [-2,2]. The result is the same as with the Runge method in this range. The Runge-Kutta method seems to work better, and the Gauss integral becomes unstable.

>longestformat; fy(2), res=runge("x*y*sin(y)",[0,2],1,500); res[2]

For a second example, we solve the problem

11 - Differential Equations

Again the solution is implicit, and we need to solve it for y.

>&ode2('diff(y,x)=-x/(1+y),y,x), sol &= rhs(solve(%,y)[2])
                            2          2
                         - y  - 2 y   x
                         ---------- = -- + %c
                             2        2

                      sqrt(- x  - 2 %c + 1) - 1

To find the constant, we can use the standard solver of Maxima instead of ic1.

>&solve((sol with x=0)=1), function ysol(x) &= sol with %
                              [%c = - -]

                           sqrt(4 - x ) - 1

Systems of Equations


11 - Differential Equations

numerically requires a rewrite into a system of equations.

11 - Differential Equations

We write a function f(x,u).

>function f(t,u) := [u[2],-u[1]^2]

The time variable t is not used in this function in this case.

Now we solve

11 - Differential Equations

with the Runge method.

>t=linspace(0,4,500); s=runge("f",t,[0,1]);

The result is a matrix with two rows y and y'. We plot both into one plot.


11 - Differential Equations

We can plot the vector field for this, since the time does not matter. Note that vector fields in the plane need two expressions in x and y, one for each component of the vector at (x,y).

>vectorfield2("y","-x^2",r=1.2); plot2d(s[1],s[2],>add):

11 - Differential Equations

We can also use a vector parameter. We solve

11 - Differential Equations

>function f(x,[y1,y2]) := [y2,-y1]; ...
 x=linspace(0,2pi,500); plot2d(x,ode("f",x,[0,1])):

11 - Differential Equations

A linear system is of the form

11 - Differential Equations

>A = [1,0.5;0.5,-2]; function f(x,u) := A.u;

We can use Maxima to derive a formula for the vector field with

11 - Differential Equations

>Axy &= @A.[x,y]
                           [  0.5 y + x  ]
                           [             ]
                           [ 0.5 x - 2 y ]

Then draw the hyperbolic vector field.


11 - Differential Equations

We now compute the track of a solution for t=0 to t=10.

>t=linspace(0,10,1000); s=runge("f",t,[0;1]); ...

11 - Differential Equations

The exact solution of such a differential equation is

11 - Differential Equations

Euler has a function to apply an analytic function to a matrix, which decomposes the matrix and applies the function to the eigenvalues.


Compare to the numerical result.


Maxima can derive a formula for this solution.

>function u(t) &= matrixexp(t*@A).[0;1];

The formula is

11 - Differential Equations

11 - Differential Equations

To print this in the comment, use

 maxima: u[1](t) = u(t)[1]
 maxima: u[2](t) = u(t)[2]

Stiff Equations

The following example

11 - Differential Equations

is very difficult to solve with the Runge method.

>c=1e-5; x=0:0.1:4; plot2d(x,ode("-y/(c+y)",x,1)):

11 - Differential Equations

But the LSODA algorithm, which is used in ode() by default, switches to a different formula if the equations becomes "stiff".

If we solve this with the ordinary Runge method for small c, we get essentially the solution of

11 - Differential Equations

>c=1e-5; x=0:0.1:4; plot2d(x,runge("-y/(c+y)",x,1)):

11 - Differential Equations

In many cases, this kind of equation simply wants to simulate a change of the behavior of f(x,y). It might be far better to use cases.

>function f(x,y) ...
   global c;
   if y>0 then return -y/(c+y);
   else return 0

The solution is not perfect, since it equals a small negative constant c>1. But it far easier to compute, and we can chop away the parts were y<0 in the plot.

>c=1e-5; x=0:0.01:4; plot2d(x,max(runge("f",x,1),0)):

11 - Differential Equations

Here is another example. We solve

11 - Differential Equations

This is a stiff equation. The simple Runge algorithm does not work. But LSODA can handle this very well.

Note that we defined f() with a vector parameter.

>function f(x,[y1,y2]) := [y2,1000*(1-y1^2)*y2-y1]; ...
 x=0:0.01:5000; y=ode("f",x,[2,0]); plot2d(x,y[1]):

11 - Differential Equations

Boundary Value Problems

Let us solve a boundary value problem numerically. The differential equation is

11 - Differential Equations

We rewrite that in two dimensions and get

11 - Differential Equations

We want the boundary values

11 - Differential Equations

We have to rewrite the equation into a system of equations.

>function f(x,y) := [y[2],x*y[1]]

We can use the LSODA algorithm to solve this for

11 - Differential Equations

>y=ode("f",[0,1],[1,0]); y[1,2]

The next step is to solve for the boundary values. We set up a function, which has a zero in the solution by taking y'(0) as a parameter and y(1)-1 as the function value.

>function g(a) ...
 return y[1,2]-1;

We know already g(0)>0, and we see g(-2)<0.

>g(0), g(-2),

Thus the secant method yields the desired value for y'(0). This is called the "Shooting Method".


A final plot of the answer.

>t=0:0.05:1; y=ode("f",t,[1,a]); plot2d(t,y[1]):

11 - Differential Equations

Here is a more complex boundary value problem, which can be solved analytically in Maxima.

11 - Differential Equations

>&ode2('diff(y,x,2)+y=sin(x),y,x),  ...
 function f(x) &= y with bc2(%,x=0,y=0,x=%pi/2,y=-1)
                                 x cos(x)
                y = %k1 sin(x) - -------- + %k2 cos(x)

                                    x cos(x)
                         - sin(x) - --------

Let us plot the solution to this boundary value problem.


11 - Differential Equations

The correct derivative in 0 is -3/2.

                                 - -

Let us try the shooting method.

>function f(x,y) := [y[2],sin(x)-y[1]]

The function g must be modified as follows.

>function g(a) ...
 return y[1,2]+1;


>g(0), g(-2),

And we do indeed find the correct derivative.

>a := bisect("g",-2,0)

Euler Home