iconEuler Home

Differential Equations

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

  ../reference/dgl

To understand this notebeook, 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 - Introduction to 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.

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

We first plot the vector field using the function "vectorfield".
>vectorfield("x*y*sin(y)",-3,3,0,4);
Then we compute and plot the right side of the solution.
>t=0:0.01:3; s=runge("x*y*sin(y)",t,1); plot2d(t,s,a=-3,b=3,c=0,d=4,add=1);
And add the left side.
>t=0:-0.01:-3; s=runge("x*y*sin(y)",t,1); plot2d(t,s,add=1):

11 - Introduction to 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.
>runge("x*y*sin(y)",[0,10],1,500)
[ 1  3.14159265359 ]
We can try to solve this equation 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 Gauß 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 instable.
>longestformat; fy(2), res=runge("x*y*sin(y)",[0,2],1,500); res[2]
      3.103165883964501 
      3.103165883778106 

Higher Order Solvers

The LSODA solver is a special solver for differential equations,
which switches to a specials state in stiff cases. So it can be used
for all differential equations. It is implemented in the ode()
function of Euler.

This method is very effective, and should be the method of choice.
Since it is an adaptive procedure, intermediate steps are automatic,
and need not be provided.
>longestformat; fy(2), res=ode("x*y*sin(y)",[0,2],1); res[-1]
      3.103165883964501 
      3.103165883948223 
There is also an adaptive Runge format, which might be more effective
sometimes, but is not most of the time.
>fy(2), res=adaptiverunge("x*y*sin(y)",[0,2],1); res[-1]
      3.103165883964501 
      3.103165883957769 
Using Maxima, we can get high order numerical solvers for
differential equations, which deliver guaranteed interval inclusions
of the solution.

The function "mxmidgl" computes the derivatives of an expression, so
that it can use a Taylor expansion of given degree. In the example,
we use degree 5.
>v:=mxmidgl("x*y*sin(y)",0:0.01:2,1,deg=5); v[-1]
              ~3.1031658633,3.1031659046~ 

Solving with Maxima

In other cases, Maxima is more successful with exact solutions.
Let us try the initial value problem y+x*y'=log(x)+1, y(1)=1.
>&ode2(y+x*'diff(y,x)=log(x)+1,y,x), &ic1(%,x=1,y=1)
                              x log(x) + %c
                          y = -------------
                                    x


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

The numerical solution is the same.
>res := ode("(log(x)+1-y)/x",[1,2],1); ...
 longestformat; res[2], (2*log(2)+1)/2
      1.193147180564709 
      1.193147180559945 

Stiff Equations

The following example

11 - Introduction to 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 - Introduction to Differential Equations

A Boundary Value Problem

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

11 - Introduction to Differential Equations

We rewrite that in two dimensions and get

11 - Introduction to Differential Equations

We want the boundary values

11 - Introduction to 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 - Introduction to Differential Equations


>y=ode("f",[0,1],[1,0]); y[1,2]
      1.172299970082665 
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) ...
t=0:0.001:1;
y=ode("f",[0,1],[1,a]);
return y[1,2]-1;
endfunction
We know already g(0)>0, and we see g(-2)<0.
>g(0), g(-2),
     0.1722999700826648 
      -1.99837932614611 
Thus the secant method yields the desired value for y'(0).
>a=secant("g",-2,2)
    -0.1587521200144439 
A final plot of the answer.
>t=0:0.05:1; y=ode("f",t,[1,a]); plot2d(t,y[1]):

11 - Introduction to Differential Equations

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

11 - Introduction to Differential Equations

>&ode2('diff(y,x,2)+y=sin(x),y,x),  ...
                                 x cos(x)
                y = %k1 sin(x) - -------- + %k2 cos(x)
                                    2

>  function f(x) &= y with bc2(%,x=0,y=0,x=%pi/2,y=-1)
                                    x cos(x)
                         - sin(x) - --------
                                       2

Let us plot the solution to this boundary value problem.
>plot2d(f,0,pi/2):

11 - Introduction to Differential Equations

The correct derivative in 0 is -3/2.
>&diffat(f(x),x=0)
                                   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) ...
y=ode("f",[0,pi/2],[0,a]);
return y[1,2]+1;
endfunction
Now
>g(0), g(-2),
      1.500000000004344 
    -0.5000000000005538 
And we do indeed find the correct derivative.
>a := bisect("g",-2,0)
      -1.50000000000091 

Euler Home