{ 26/01/84 011141605 HLL RKF45 PASCAL PASCAL } const NEQN = 20; type workarray = array[1..NEQN] of real; procedure rkf45 ( n : integer; { number of equations } function f { compute y_prime = dy / dx; return true if ok; else false } ( n : integer; { number of equations } x : real; { independent variable } y : workarray; { dependent variable } var y_prime : workarray { dy / dx } ) : boolean; var y_prime : workarray; { derivatives of y computed by f } x0 : real; { initial value of independent variable } var y0 : workarray; { initial value of dependent variable } x1 : real; { final value of independent variable } epsilon : workarray; { integration tolerance, expressed as an absolute value, not a fraction } firsth, { initial step size } minh { smallest step size tolerated } : real; var y1 : workarray; { final value of dependent variable; only returned if success is true } var lasth : real; { last value of step size; recommended for restart if desired } var failure : integer { was integration successful? } ); {----------------------------------------------------------------- | Runge-Kutta-Fehlberg integration of the system of n ordinary | | differential equations | | | | dy(x) [i] | | ----- = y_prime(x) [i]; | | dx | | | | y(x0) [i] = y0 [i]; --> y(x1) [i] = y1 [i]; | | | | (i = 1...n) | | | | using the 5-stage fourth order Runge-Kutta method of Fehlberg. | | | | The function f(n, x, y, y_prime) computes y_prime from the | | values of x and y(x) and returns true if the calculation is | | successful, else false. | | | | In application rkf45 should be used in a series of steps | | to cover the desired total interval; the method will take | | finer steps within this interval as required. Internally | | the step size is halved if the error estimate per unit step | | exceeds epsilon. The step size is doubled if the error | | estimate is less than epsilon / 32, unless this exceeds the | | range available. The process terminates if the step size | | becomes less than minh. The variable "failure" is 0 if the | | integration is successful, 1 if the number of equations | | exceeds NEQN, 2 if the step size became less than minh without | | reaching the requested error tolerance, 3 if the derivative | | calculation failed. | | | | Taken from An Introduction to Numerical Methods with Pascal | | p. 283 by L.V. Atkinson and P.J. Harley, (Addison Wesley). | | | | Modified | | (1) to protect the routine from bad step size and collapsed | | x interval, | | (2) to allow stepping in the negative direction, | | (3) to return several error status indications, | | (4) to work in n-dimensions. | | (5) better value of nearzero (roundoff handling) | | (6) pass y0 and y_prime as var parameters. | | | | 26 January 1984 | | H. L. Lynch | | | ------------------------------------------------------------------} const mach_acc = 1.0e-15; { relative machine precision for real } var yt, { temporary storage for calc. kn } k1, k2, { interpolated devivative values } k3, k4, { interpolated derivative values } k5, k6 { interpolated derivative values } : workarray; nearzero, { within machine accuracy of end point } error, { error estimate } x, { curr. value of independent variable} h, { current step size } ah, { absolute value of current step size} tolerance { temp. var. in error tolerance calc.} : real; state { current state of program } : (stepping, x1reached, htoosmall, bad_derivative); good_deriv, { good derivative returned from f() } err_too_large, { need to reduce step size } err_too_small, { need to increase step size } shortrange { flag to prevent step beyond x1 and return good estimated h } : boolean; i { dummy index for arrays } : integer; begin { Runge-Kutta-Fehlberg rkf45 } if (n > NEQN) then { Attempt to over run internal arrays } begin failure := 1; return; end; x := x0; for i := 1 to n do y1[i] := y0[i]; shortrange := false; state := stepping; h := firsth; { check for bad input } if ((h * (x1 - x0)) < 0.0) then { step is in wrong direction } h := -h; if ((((x1 - x0) / h) < 1) or (h = 0)) then { bad step size } begin h := x1 - x0; if (abs(h) < minh) then { interval has collapsed } return; end; ah := abs(h); nearzero := mach_acc * abs(x1 - x0); repeat { while state = stepping, checked at the end } good_deriv := f(n, x, y1, y_prime); for i := 1 to n do k1[i] := y_prime[i]; for i := 1 to n do yt[i] := y1[i] + h * k1[i] / 4.0; good_deriv := f(n, x + h / 4.0, yt, y_prime) and good_deriv; for i := 1 to n do k2[i] := y_prime[i]; for i := 1 to n do yt[i] := y1[i] + h * (3.0 * k1[i] + 9.0 * k2[i]) / 32.0; good_deriv := f(n, x + 3.0 * h / 8.0, yt, y_prime) and good_deriv; for i := 1 to n do k3[i] := y_prime[i]; for i := 1 to n do yt[i] := y1[i] + h * (1932.0 * k1[i] - 7200.0 * k2[i] + 7296.0 * k3[i]) / 2197.0; good_deriv := f(n, x + 12.0 * h / 13.0, yt, y_prime) and good_deriv; for i := 1 to n do k4[i] := y_prime[i]; for i := 1 to n do yt[i] := y1[i] + h * (439.0 * k1[i] / 216.0 - 8.0 * k2[i] + 3680.0 * k3[i] / 513.0 - 845.0 * k4[i] / 4104.0); good_deriv := f(n, x + h, yt, y_prime) and good_deriv; for i := 1 to n do k5[i] := y_prime[i]; for i := 1 to n do yt[i] := y1[i] + h * (-8.0 * k1[i] / 27.0 + 2.0 * k2[i] - 3544.0 * k3[i] / 2565.0 + 1859.0 * k4[i] / 4104.0 - 11.0 * k5[i] / 40.0); good_deriv := f(n, x + h / 2.0, yt, y_prime) and good_deriv; for i := 1 to n do k6[i] := y_prime[i]; if (not good_deriv) then state := bad_derivative; { Set up err_too_large := at least one error is too large; err_too_small := all errors are too small. } err_too_large := false; err_too_small := true; for i := 1 to n do begin error := ah * abs(k1[i] / 360.0 - 128.0 * k3[i] / 4275.0 - 2197.0 * k4[i] / 75240.0 + k5[i] / 50.0 + 2 * k6[i] / 55.0); tolerance := ah * epsilon[i]; err_too_large := (error > tolerance) or err_too_large; err_too_small := (error <= tolerance / 32.0) and err_too_small; end; { end i loop } if (err_too_large) then begin h := h / 2.0; ah := abs(h); if (ah < minh) then state := htoosmall; end else begin { Good result, take the step. } x := x + h; for i := 1 to n do begin y1[i] := y1[i] + h * (25.0 * k1[i] / 216.0 + 1408.0 * k3[i] / 2565.0 + 2197.0 * k4[i] / 4104.0 - k5[i] / 5.0); end; if (err_too_small) then begin h := h * 2.0; ah := abs(h); end; if (abs(x1 - x) <= nearzero) then state := x1reached else begin shortrange := (ah > abs(x1 - x)); if (shortrange) then begin lasth := h; h := x1 - x; ah := abs(h); end; end; { end x1 not yet reached } end; { end err not too large } until (state <> stepping); if (not shortrange) then lasth := h; if (state = x1reached) then failure := 0 else if (state = htoosmall) then failure := 2 else if (state = bad_derivative) then failure := 3; end { end Runge-Kutta-Fehlberg rkf45 } ;