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
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):
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
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~
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
The following example
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)):
Let us solve a boundary value problem numerically. The differential equation is
We rewrite that in two dimensions and get
We want the boundary values
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
>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]):
Here is a more complex boundary value problem, which can be solved analytically in Maxima.
>&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):
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