25 lines
781 B
Text
25 lines
781 B
Text
Runge_Kutta: procedure options (main); /* 10 March 2014 */
|
|
declare (y, dy1, dy2, dy3, dy4) float (18);
|
|
declare t fixed decimal (10,1);
|
|
declare dt float (18) static initial (0.1);
|
|
|
|
y = 1;
|
|
do t = 0 to 10 by 0.1;
|
|
dy1 = dt * ydash(t, y);
|
|
dy2 = dt * ydash(t + dt/2, y + dy1/2);
|
|
dy3 = dt * ydash(t + dt/2, y + dy2/2);
|
|
dy4 = dt * ydash(t + dt, y + dy3);
|
|
|
|
if mod(t, 1.0) = 0 then
|
|
put skip edit('y(', trim(t), ')=', y, ', error = ', abs(y - (t**2 + 4)**2 / 16 ))
|
|
(3 a, column(9), f(16,10), a, f(13,10));
|
|
y = y + (dy1 + 2*dy2 + 2*dy3 + dy4)/6;
|
|
end;
|
|
|
|
|
|
ydash: procedure (t, y) returns (float(18));
|
|
declare (t, y) float (18) nonassignable;
|
|
return ( t*sqrt(y) );
|
|
end ydash;
|
|
|
|
end Runge_kutta;
|