PROGRAM DiffDemo;
{ Program to demonstrate numerical solutions to ordinary, first order
differential equations (initial value problems). Starting from
(x, y) = (0, y0), the solution to y' = dy/dx is advanced through
N steps of size h to a final point (Xf, Yf).
Two "self-starting" methods are used here. The primitive Euler method
introduces the topic, while the 4'th order Runge-Kutta formula is
useful in practice.
Note that many other more sophisticated and accurate methods of numerical
solution of differential equations exist.
This program was compiled with Turbo Pascal V5.5, but should be
compilable with Turbo Pascal versions 3.0 and above. Deviations from
standard Pascal are also minimal, particularly in the coding of the
algorithms themselves. }
VAR
I, N: Integer;
X, Xf, h,
Y0, Y_True,
Y_Euler,
Y_RK4: Real;
FUNCTION F(x, y: Real): Real;
{ The function to be integrated. I.e. y' = dy/dx = F(x, y) }
BEGIN
F:= y;
END; { FUNCTION F }
FUNCTION Y_solution(x, y0: Real): Real;
{ The analytical solution y = y(x) to the differential equation
y' = F(x, y) with initial condition y(0) = y0. }
BEGIN
Y_solution:= exp(x) + Y0;
END; { FUNCTION Y_solution }
(* FUNCTION F(x, y: Real): Real;
{ The function to be integrated. I.e. y' = dy/dx = F(x, y) }
BEGIN
F:= -y*y;
END; { FUNCTION F }
FUNCTION Y_solution(x, y0: Real): Real;
{ The analytical solution y = y(x) to the differential equation
y' = F(x, y) with initial condition y(0) = y0. }
BEGIN
Y_solution:= 1/(x + 1/y0);
END; { FUNCTION Y_solution } *)
FUNCTION Euler(x, y, h: Real): Real;
{ Uses F(x, y) to advance the solution of the first order
ordinary differential equation from x to x+h. Note that
the "truncation error" in the first order Euler method is
O(h), so the estimated solution is not accurate. }
BEGIN
Euler:= y + h*F(x, y);
END; { FUNCTION Euler }
FUNCTION Runge_Kutta4(x, y, h: Real): Real;
{ Uses F(x, y) to advance the solution of the first order
ordinary differential equation from x to x+h. The truncation
error is O(h^4), so the method is reasonable accurate despite
the simplicity of coding. }
VAR
h_half,
x_half,
y_half_1,
y_half_2,
y_1: Real;
BEGIN
h_half := h/2;
x_half := x + h_half;
{ Form various estimates of the solution at x+h/2 and x+h, then
combine them using the formula below. }
y_half_1 := Euler(x, y, h_half);
y_half_2 := Euler(x_half, y_half_1, h_half);
y_1 := Euler(x_half, y_half_2, h);
Runge_Kutta4:= y + h/3 *
((f(x, y) + f(x+h, y_1))/2 +
f(x_half, y_half_1) + f(x_half, y_half_2));
END; { FUNCTION Runge_Kutta4 }
BEGIN { Main }
{ Get input from user---target x, Xf, and number of steps, N }
REPEAT
Write('Enter final (high) x value: ');
ReadLn(Xf);
UNTIL (Xf > 0);
REPEAT
Write('Enter number of steps to take: ');
ReadLn(N);
UNTIL (N > 0);
h:= Xf/N; { Step size }
X:= 0.0; { Initial X }
Y0:= 1; { Initial Y (given as part of initial value problem) }
WriteLn; WriteLn;
WriteLn('The differential equation will be solved ');
WriteLn('from x = 0 to x = ', Xf:7:2);
WriteLn('using ', N:1, ' steps of size ', h:7:4);
WriteLn; WriteLn;
WriteLn('Step':5, 'x value':10, 'True y':10,
'Euler':10, 'Error':10, 'RK-4':10, 'Error':10);
Write(0:5, X:10:4, Y0:10:4);
FOR I:= 1 TO 4 DO Write('--':10); WriteLn;
Y_Euler:= Y0;
Y_RK4:= Y0;
{ Compare true solution with Euler and Runge-Kutta 4'th order
at each step }
FOR I:= 1 TO N DO
BEGIN
Y_Euler:= Euler(X, Y_Euler, h);
Y_RK4:= Runge_Kutta4(X, Y_RK4, h);
X:= X + h;
Y_True:= Y_solution(X, Y0);
WriteLn(I:5, X:10:4, Y_True:10:4,
Y_Euler:10:4, (Y_Euler - Y_True):10:4,
Y_RK4:10:4, (Y_RK4 - Y_True):10:4);
END; { FOR I }
END.