Re: [LF] Runge-Kutta 4 (or 4-5 adaptive) Fortran 90 code

[Lahey Support]

I looked at DERKF, which someone suggested, and it was a real
dog's breakfast of old Fortran 66 code. The following version
is adapted from the NSWC Math Library (which has been placed
in the public domain).


SUBROUTINE rk(n, t, h, a, f)

! Code converted using TO_F90 by Alan Miller
! Date: 2000-02-19 Time: 15:41:29

! ************************************************** ****************
! FOURTH ORDER RUNGE-KUTTA PROCEDURE FOR SOLVING DY = F(T,Y)
! ************************************************** ****************

! Let y'(t) = f(t,y(t)) denote a system of n ordinary first order nonstiff
! differential equations where f(t,y) = {f1(t,y), ..., fn(t,y)} =
! {y1(y), ..., yn(t)}. Assume that y(t0) is known. Then for a small
! real number h, the routine RK is available for computing y(t0+h).

IMPLICIT NONE
INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(12, 60)

INTEGER, INTENT(IN) :: n
REAL (dp), INTENT(IN OUT) :: t
REAL (dp), INTENT(IN) :: h
REAL (dp), INTENT(OUT) :: a(:)

! EXTERNAL f

INTERFACE
SUBROUTINE f(t, z)
IMPLICIT NONE
INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(12, 60)
REAL (dp), INTENT(IN) :: t
REAL (dp), INTENT(OUT) :: z(:)
END SUBROUTINE f
END INTERFACE

REAL (dp) :: ha, ta
INTEGER :: k, m, m1, mk, nk, np1

np1 = n + 1
IF (h /= 0.0_dp) THEN

ha = .5_dp * h
ta = t + ha
m = n + n
m1 = m + 1

DO k = 1, n
nk = n + k
mk = m + k
a(nk) = ha * a(nk)
a(mk) = a(k) + a(nk)
END DO
CALL f(ta, a(m1:))

DO k = 1, n
nk = n + k
mk = m + k
a(nk) = a(nk) + h * a(mk)
a(mk) = a(k) + ha * a(mk)
END DO
CALL f(ta, a(m1:))

t = t + h
DO k = 1, n
nk = n + k
mk = m + k
a(mk) = h * a(mk)
a(nk) = a(nk) + a(mk)
a(mk) = a(k) + a(mk)
END DO
CALL f(t, a(m1:))

DO k = 1, n
nk = n + k
mk = m + k
a(k) = a(k) + (a(nk) + ha*a(mk)) / 3.0_dp
a(nk) = a(k)
END DO
CALL f(t, a(np1:))
RETURN
END IF

DO k = 1, n
nk = n + k
a(nk) = a(k)
END DO
CALL f(t, a(np1:))

RETURN
END SUBROUTINE rk


Alan.Miller @ vic.cmis.CSIRO.AU = Alan Miller, Honorary Research Fellow (Retired)
CSIRO Mathematical & Information Sciences, Melbourne, Australia
Phone: +61 3 9545-8036 Fax: +61 3 9545-8080 URL: http://www.ozemail.com.au/~milleraj
Mail: CSIRO CMIS, Private Bag 10, Clayton South MDC, Vic. 3169, Australia
-----------------------------------------------------
To unsubscribe from the Fortran Forum, send a message
to [address removed] with the following command as
the first and only line of the message body:
unsubscribe fortran
-----------------------------------------------------