# INTRODUCTORY DIFFERENTIAL EQUATIONS

    by Dr. Wolfgang Lindner, Leichlingen, Germany, 2023

## What is an ordinary Differential Equation?

## Analytical Solution Methods

## Intermezzo: iterate

## Numerical Solution Methods for IVP

fpprintprec:5$
iterate(u,x,x0,n) :=  block([numer:true],
                      cons(x0, makelist(x0:subst(x0,x,u), i,1,n)))$

iterate2(u,x,x0,n) := block([numer:true],
  cons(x0, makelist(x0:subst([x[1]=x0[1],x[2]=x0[2]],u), i,1,n)))$

iterate3(u,x,x0,n) := block([numer:true],
  cons(x0, makelist(x0:subst([x[1]=x0[1],x[2]=x0[2],x[3]=x0[3]],u),i,1,n)))$

### RK1 alias the EULER Method

ye: exp(x^2/2);      /* exact solution  ye */
yp: diff(ye, x);     /* their slope ye' */

E(x,y) := [x + 0.2, y +0.2* x*y];
       /* xn + h    Yn+  h*f(xn,Yn)  */
ratprint:false$ 
fpprintprec:5$ 
E(0.0, 1);
E(0.2, 1);
E(0.4, 1.04);
E(0.6, 1.1232);
E(0.8, 1.257984);

XY: [E(0,1),E(0.2,1),E(0.4,1.04),E(0.6,1.1232),E(0.8,1.257984)];
wxdraw2d(xaxis = true, 
    point_size=2, points_joined=true,
    xlabel="xi", ylabel="Yi",
    points(XY),
    color=red,
    explicit(exp(x^2/2),x,0,1),  /* exact solution ye */
    title="EULER polygon")$

#### The Euler Method for IVP

fpprintprec:5$
iterate2(u,x,x0,n) := block([numer:true],
  cons(x0, makelist(x0:subst([x[1]=x0[1],x[2]=x0[2]],u), i,1,n)));

f(x,y):= x*y;
xo: 0;
yo: 1;
h: 0.2;
n: 5;   /* recurrence__term    vars.   initial  steps */
iterate2( [x+h,  y+h*f(x,y)],  [x,y],  [xo,yo], n     );

U(x) := exp(x^2/2); 
U(1), numer;

/* MAXIMA : --- EULER method  --- */
EULER(x,y, xo,yo, h, n) := iterate2([x+h, y+h*f(x,y)],[x,y],[xo,yo],n);  

/* Test ODE  y'= xy with y(0)=1 */
fpprintprec: 5$ 
f(x,y):=x*y;
EULER( x,y, 0, 1, 0.2, 5);  

/* MAXIMA : --- EULER method denoted as RK1  --- */
RK1(x,y, xo,yo, h, n) := block(
    k1: f(x,y),
    iterate2([x+h, y+h*k1],[x,y],[xo,yo],n));

RK1( x,y, 0, 1, 0.2, 5); 

        ----- EXERCISE 88:

fpprintprec:5$
second(z):= z[2];
ye(x):= exp(x^2/2);
f(x,y):=x*y;
Pe: EULER( x,y, 0,1, 0.2, 5);             /* Euler method points */
L1: map(second,Pe);
L2: makelist(ye(x),x,0,1,0.2);
abs(L2-L1);
localError: lmax(abs(L2-L1));
globalError: sum(abs(L2-L1)[i],i,1,length(L1));

#### Modified Euler method alias 'RK2e of order 2'

U(x) := exp(x^2/2); 
U(1), numer;           /* exact value */

fpprintprec:5$
f(x,y) := x*y; 
k1 : f(x,y);
k2 : f(x+0.2, y + 0.2*k1);
iterate2([x+0.2, y + 0.2*(k1+k2)/2],[x,y],[0,1],5);

iterate2(u,x,x0,n) := block([numer:true],
        cons(x0, makelist(x0:subst([x[1]=x0[1],x[2]=x0[2]],u), i,1,n)))$

/* MAXIMA : --- modified EULER method --- */
modiEULER(x,y, xo,yo, h,n):= block([k1,k2], 
               k1 : f(x,y),
               k2 : f(x+h, y + h*k1),
               iterate2([x+h, y+h/2*(k1+k2)],[x,y],[xo,yo],n) )$
                                
/* Test using ODE  y'= xy with y(0)=1 */
f(x,y) := x*y;
modiEULER( x,y, 0,1, 0.2, 5);

RK2e( x,y, xo,yo, h,n) :=
       iterate2([x+h, y+h/2*( f(x,y)  + f(x+h,y+h*f(x,y)) )], 
                [x,   y], [xo,yo], n)$
 
f(x,y) := x*y;                               
RK2e( x,y, 0,1, 0.2, 5);

### RK2h alias HEUN's Method

fpprintprec: 5$ 
iterate2(u,x,x0,n) := block([numer:true],
        cons(x0, makelist(x0:subst([x[1]=x0[1],x[2]=x0[2]],u), i,1,n)))$

/* MAXIMA --- HEUN's method --- */
HEUN( x,y, xo,yo, h,n) := block(  
           k1 : f(x,y),
               k2 : f( x+2*h/3,  y+2*h/3*k1 ),
               iterate2([x+h, y+h/4*(k1+3*k2)],[x,y],[xo,yo],n) );
                                 
/* Test with ODE  y'= xy with y(0)=1 */
f(x,y):= x*y;                
HEUN( x,y, 0,1, 0.2, 5);

fpprintprec:5$
k(z):= z[2];
ye(x):=exp(x^2/2);
f(x,y):= x*y;   
P : HEUN(x,y, 0,1, 0.2, 5);
L1: map(k,P);
L2: makelist(ye(x),x,0,1,0.2);
abs(L2-L1);
globalErrorHEUN: sum(abs(L2-L1)[i],i,1,length(L1));

### RK2m alias the Midpoint Method

fpprintprec: 5$ 
iterate2(u,x,x0,n) := block([numer:true],
        cons(x0, makelist(x0:subst([x[1]=x0[1],x[2]=x0[2]],u), i,1,n)))$

/* MAXIMA : --- MIDPOINT METHOD --- */
midpoint(x,y, xo,yo, h,n) := 
         iterate2([x+h, y+h*f( x+h/2, y+h/2*f(x,y) ) ],[x,y],[xo,yo],n);
	
/* Test with ODE  y'= xy with y(0)=1 */
f(x,y):=x*y;
midpoint( x,y, 0,1, 0.2, 5);

fpprintprec:5$
k(z):= z[2];
ye(x):=exp(x^2/2);
f(x,y):=x*y;
P : midpoint(x,y, 0,1, 0.2, 5);
L1: map(k,P);
L2: makelist(ye(x),x,0,1,0.2);
abs(L2-L1);
globalErrorMidpoint: sum(abs(L2-L1)[i],i,1,length(L1));

        ----- EXERCISE 105:

/* MAXIMA : --- MIDPOINT METHOD as RK2 method --- */
midpoint1( x,y, xo,yo, h,n) :=  block(  
               k1 : f(x ,y),  
               k2 : f(x, y + h/2*k1), 
	            iterate2([x+h, y+h*k2], [x,y], [xo,yo], n));

/* Test with ODE  y'= xy with y(0)=1 */
f(x,y):=x*y;
midpoint1(x,y, 0,1, 0.2, 5);

/* MAXIMA --- Midpoint method variant--- */
midpoint2(x,y, xo,yo, h,n) := 
         iterate2([x+h, 
                   y+1/2*h*( f(x,y)+f(x+h,y+h)) ], [x,y], [xo,yo], n);
	
f(x,y):=x*y;
midpoint2( x,y, 0,1, 0.2, 5);

### RK4 alias the classic RUNGE-KUTTA method

kill(all)$
fpprintprec: 5$ 

iterate2(u,x,x0,n) := block([numer:true],
        cons(x0, makelist(x0:subst([x[1]=x0[1],x[2]=x0[2]],u), i,1,n)))$

#### the Runge–Kutta RK4 recurrence

/* MAXIMA : --- classic RUNGE-KUTTA  method of order 4 --- */
RK4(x,y,xo,yo,h,n):= block(
                     k1 : f(x,y),
  	                k2 : f(x+h/2, y + h/2*k1),
	                  k3 : f(x+h/2, y + h/2*k2),
	                  k4 : f(x+h  , y + h*k3),
	                     iterate2([x+h, y+h/6*(1*k1+2*k2+2*k3+1*k4)], 
	                              [x  , y], [xo,yo], n) )$
	                
/* Test ODE: y'= xy with y(0)=1 */
f(x,y):= x*y;
RK4( x,y, 0,1, 0.2, 5);

#### the build-in rk(.) method

results: rk(x-y^2, y,1, [x, 0,8, 0.1])$
wxdraw2d(xaxis = true, 
    points_joined=true,
    xlabel="xi", ylabel="Yi",
    points(results),
    title="rk(.) polygon")$

       ----- EXERCISE 111:

iterate2(u,x,x0,n) := block([numer:true],
        cons(x0, makelist(x0:subst([x[1]=x0[1],x[2]=x0[2]],u), i,1,n)))$
/* MAXIMA --- RUNGE-KUTTA order 3 --- */
RK3(x,y,xo,yo,h,n) := block(
                   k1 :  f(x,y),
	                k2 : f(x+h/2, y + h/2*k1),
	                k3 : f(x+h, y + 2*h*k2-h*k1),
	                iterate2([x+h, y+h/6*(1*k1+4*k2+1*k3)],[x,y],[xo,yo],n))$
	                 
/* Test ODE: y'= xy with y(0)=1 */
f(x,y):= x*y;
RK3(x,y, 0,1, 0.2, 5);