#20 (04/07/2025)

Statement of Differential Equations

= 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(t₀) = y₀,

(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

y_n+1−y_n

≈ f(t, y).

(3)

Equation (3) can be written as

y_n+1 = y_n + h f(t, y),

(4)

which can be used to predict y_{n + 1} from y_n and the slope at that point.

Example:

= y, y(0) = 1.

(5)

(Analytical solution) Separation of variables:

y = e^t.

(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=%f %f %f\n", t, y, exp(t));
   y = y + h * f(t,y);
   t = t + h;
}
 return 0;
}

Runge-Kutta method

k₁

f(t, y)

k₂

f(t +

, y +

k₁)

k₃

f(t +

, y +

k₂)

k₄

f(t + h, y + h k₃),

y_n+1 = y_n + h

k₁ + 2 k₂ + 2 k₃ + k₄

(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

d² y

dt²

= − y, y(0)=0, y′(0)=1.

(8)

By setting

y₁ ≡ y, y₂ ≡

dy₁

the equation above is converted to

dy₁

y₂

(9)

dy₂

−y₁,

(10)

with

y₁(0)=0, y₂(0)=1.

The exact solutions are y₁=sin(t) and y₂=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)

p (v − u)

(11)

−u w + R u − v

(12)

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');

Matalb cheat for differential equations

The following code numerically solves the differential equation:

= −t² y, y(0) = 1.

f=@(t,y) -t*t*y;
tspan=linspace(0,5,100);
y0=1;
[t,y]=ode45(f, tspan, y0);
plot(t,y);

Octave cheat for differential equations

The following code numerically solves the differential equation:

= −t² y, y(0) = 1.

f=@(y,t) -t*t*y;
tspan=linspace(0, 3, 100);
y0=1;
sol=lsode(f, y0, tspan);
plot(tspan, sol);

File translated from T_EX by T_TH, version 4.03.
On 08 Apr 2025, 21:09.