26 lines
1.1 KiB
Text
26 lines
1.1 KiB
Text
BEGIN
|
|
PROC rk4 = (PROC (REAL, REAL) REAL f, REAL y, x, dx) REAL :
|
|
BEGIN CO Fourth-order Runge-Kutta method CO
|
|
REAL dy1 = dx * f(x, y);
|
|
REAL dy2 = dx * f(x + dx / 2.0, y + dy1 / 2.0);
|
|
REAL dy3 = dx * f(x + dx / 2.0, y + dy2 / 2.0);
|
|
REAL dy4 = dx * f(x + dx, y + dy3);
|
|
y + (dy1 + 2.0 * dy2 + 2.0 * dy3 + dy4) / 6.0
|
|
END;
|
|
REAL x0 = 0, x1 = 10, y0 = 1.0; CO Boundary conditions. CO
|
|
REAL dx = 0.1; CO Step size. CO
|
|
INT num points = ENTIER ((x1 - x0) / dx + 0.5); CO Add 0.5 for rounding errors. CO
|
|
[0:num points]REAL y; y[0] := y0; CO Grid and starting point.CO
|
|
PROC dy by dx = (REAL x, y) REAL : x * sqrt(y); CO Differential equation. CO
|
|
FOR i TO num points
|
|
DO
|
|
y[i] := rk4 (dy by dx, y[i-1], x0 + dx * (i - 1), dx)
|
|
OD;
|
|
print ((" x true y calc y relative error", newline));
|
|
FOR i FROM 0 BY 10 TO num points
|
|
DO
|
|
REAL x = x0 + dx * i;
|
|
REAL true y = (x * x + 4.0) ^ 2 / 16.0;
|
|
printf (($3(-zzd.7dxxx), -d.4de-ddl$, x, true y, y[i], y[i] / true y - 1.0))
|
|
OD
|
|
END
|