#18 (04/15/2024)

Statement of Differential Equations


dy
dt
 = f(t, y),
(1)
where t is the independent variable and f(t, y) is a function of t and y. Equation (1) along with the initial boundary condition of
y(t0) = y0,
(2)
is called the initial value problem.

Euler's method

As most of the useful and important differential equations cannot be solved analytically, one has to resort to numerical techniques.
The Euler method is the most primitive numerical technique and should be used as a first approximation.
In the Euler method, d y / d t is approximated by its forward difference as

yn+1−yn
h
 ≈ f(t, y).
(3)
euler.gif
Equation (3) can be written as

yn+1 = yn + h f(t, y),
(4)
which can be used to predict yn + 1 from yn and the slope at that point.
Example:
dy
dt
 = y,     y(0) = 1.
(5)
(Analytical solution) Separation of variables:
y = et.
(6)
#include <stdio.h>
#include <math.h>

double f(double t, double y)
{
return y;
}

int main()
{
double h=0.1, y, t;
int i;

t=0.0; y = 1.0;

for (i=0; i<= 10; i++)
{
  printf("t= %lf %lf %lf\n", t, y, exp(t));
  y= y+h*f(t,y);
  t=t+h;
}
return 0;
}

Runge-Kutta method


k1
=
f(t, y)
k2
=
f(t + h
2
 , y + h
2
 k1)
k3
=
f(t + h
2
 , y + h
2
 k2)
k4
=
f(t + h, y + h k3),

yn+1 = yn + h k1 + 2 k2 + 2 k3 + k4
6
 .
(7)
#include <stdio.h>
#include <math.h>

double f(double t, double y)
{return y;}

main()
{
double h=0.1, t, y, k1,k2,k3,k4;
int i;

/* initial value */
t=0.0; y=1.0;

for (i=0; i<=10; i++)
{

 printf("t= %lf rk= %lf exact=%lf\n", t, y, exp(t));
 k1=h*f(t,y);
 k2=h*f(t+h/2, y+k1/2.0);
 k3=h*f(t+h/2, y+k2/2.0);
 k4=h*f(t+h, y+k3);
 y= y+(k1+2.0*k2+2.0*k3+k4)/6.0;
 t=t+h;
}
return 0;
}

Higher order ODE

Example
d2 y
dt2
 = − y,     y(0)=0,     y′(0)=1.
(8)
By setting
y1 ≡ y,     y2 ≡ dy1
dt
 ,
the equation above is converted to

dy1
dt
=
y2
(9)
dy2
dt
=
−y1,
(10)
with
y1(0)=0,     y2(0)=1.
The exact solutions are y1=sin(t) and y2=cos(t). 
#include <stdio.h>
#include <math.h>

double f1(double t, double y1, double y2)
{return y2;}

double f2(double t, double y1, double y2)
{return -y1;}

int main()
{
double h=0.02, y1, y2, t;
int i;

y1=0.0; y2=1.0;
t=0.0;

for (i=0; i<=100; i++)
{
 /*printf("t= %lf %lf %lf\n", t, y1, sin(t));*/
 printf("%lf %lf %lf\n", t, y1, sin(t));
 y1=y1+h*f1(t,y1,y2);
 y2=y2+h*f2(t,y1,y2);
 t=t+h;
}
return 0;
}


Example (Lorenz equations/Strange Attracor/Chaos)

du
dt
=
p (v − u)
(11)
dv
dt
=
−u w + R u − v
(12)
dw
dt
=
u v − b w
(13)
#include <stdio.h>
#include <math.h>
#define P 16.0
#define b 4.0
#define R 35.0

double f1(double u, double v, double w)
{return P*(v-u);}

double f2(double u, double v, double w)
{return -u*w+R*u-v;}

double f3(double u, double v, double w)
{return u*v-b*w;}

int main()
{
double h, t, u, v, w;
int i;

/* initial value */
t=0.0; h=0.01;
u=5.0; v=5.0; w=5.0;

for (i=0; i< 3000; i++)
{
 u=u+h*f1(u,v,w);
 v=v+h*f2(u,v,w);
 w=w+h*f3(u,v,w);
 printf("%lf %lf\n", u, w);
 t=t+h;
}
 return 0;
}

c:\tmp\gcc -o lorenz.exe lorenz.c
c:\tmp\lorenz.exe > lorenz.dat
(Launch gnuplot)
Change directory to c:\tmp
plot "lorenz.dat" with line
Here is a code in Matlab/Octave. 
P=16;b=4;R=35;
f1=@(u,v,w) P*(v-u);
f2=@(u,v,w) -u*w+R*u-v;
f3=@(u,v,w) u*v-b*w;

t=0;h=0.01;
u=5;v=5;w=5;

for i=1:1:3000
 u=u+h*f1(u,v,w);
 v=v+h*f2(u,v,w);
 w=w+h*f3(u,v,w);
 uu(i)=u;vv(i)=v; 
 t=t+h;
end;

plot(uu',vv');

Octave cheat for differential equations

The following code numerically solves the differential equation:

dy
dt
 = −t2 y,     y(0) = 1.
f=@(y,t) -t*t*y;
t=linspace(0, 3, 100);
y0=1;
sol=lsode(f, y0, t);
plot(t, sol);


File translated from TEX by TTH, version 4.03.
On 16 Apr 2024, 15:12.