Geology Reference
In-Depth Information
The subroutine REL is also used in monitoring stepsize in the Runge-Kutta integ-
ration to ensure that relative error is kept below the specified value of EPS. As out-
lined in Press
et al.
(1992, pp. 708-712) the Runge-Kutta integration is performed
forasinglestepofsize2
h
and then for two steps of size
h
. Since the truncation
error is
O
h
5
, the latter gives a truncation error of order 2
h
5
, the former a trunca-
tion error of order (2
h
)
5
. The di
ff
erence of the second solution minus the first is
then of order 30
h
5
. The second solution can then be improved by adding 1/15 of
this di
erence, and the subroutine REL is used to determine the maximum relative
error. If it is larger than the specified value, EPS, the stepsize needs to be reduced; if
smaller, the stepsize can be increased. In the latter case, the stepsize is increased in
proportion to the fifth root of the ratio EPS
ff
ERRMAX. If the stepsize is decreased,
there is a potential that additional error could accumulate in proportion to the num-
ber of extra steps. As a precaution, the stepsize is only decreased in proportion to
the fourth root of the ratio EPS
/
ERRMAX. The constant of proportionality for both
increases and decreases of stepsize is taken as the conservative value of 0.9. The
subroutine REL, listed below, finds the maximum ratio ERRMAX of elements of
the array DELTA(6,3) to those of the array YIC(6,3), for solutions of degree N.
/
SUBROUTINE REL(ERRMAX,DELTA,YIC,N)
C
C This subroutine is used in the variable stepsize Runge-Kutta
C integration of the fundamental solutions in the inner core.
C It finds the maximum relative error in the elements of
C the difference array, DELTA(6,3), compared to the elements of
C the fundamental solution array, YIC(6,3), for degree N.
C
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
DIMENSION DELTA(6,3),YIC(6,3)
C Zero relative error.
ERRMAX=0.D0
DO 10 I=1,6
C Branch for degree zero.
IF(N.EQ.0)GO TO 11
GO TO 12
C If degree is zero, skip y3 and y4.
11 IF(I.EQ.3.OR.I.EQ.4)GO TO 10
12 DO 13 J=1,3
C Third solution is only fundamental solution for degree zero.
IF(N.EQ.0.AND.J.NE.3)GO TO 13
IF(YIC(I,J).EQ.0.D0)GO TO 13
RATIO=DABS(DELTA(I,J)/YIC(I,J))
ERRMAX=DMAX1(ERRMAX,RATIO)
13
CONTINUE
10
CONTINUE
RETURN
END
The subroutine RK4, listed below, propagates the solution vector, over step-
size H, by the fourth-order Runge-Kutta scheme (3.297), with calls to the sub-
routine YPRIME, described in Section 3.6, for the calculation of derivatives.
Search WWH ::
Custom Search