S EC . 9.5

RUNGE -K UTTA M ETHODS

497

Runge-Kutta-Fehlberg Method (RKF45)
One way to guarantee accuracy in the solution of an I.V.P. is to solve the problem twice
using step sizes h and h/2 and compare answers at the mesh points corresponding to
the larger step size. But this requires a significant amount of computation for thesmaller step size and must be repeated if it is determined that the agreement is notgood enough.
The Runge-Kutta-Fehlberg method (denoted RKF45) is one way to try to resolve
this problem. It has a procedure to determine if the proper step size h is being used. At
each step, two different approximations for the solution are made and compared. If thetwo answers are in close agreement, the approximation is accepted. If the two answers
do not agree to a specified accuracy, the step size is reduced. If the answers agree tomore significant digits than required, the step size is increased.
Each step requires the use of the following six values:(28)

k1 = h f (tk , yk ),11
k2 = h f tk + h, yk + k1 ,44339
k3 = h f tk + h, yk + k1 + k2 ,
832
32
12
1932
7200
7296
k1 −
k2 +k3 ,k4 = h f tk + h, yk +132197
2197
2197
439
3680
845
k1 − 8k2 +
k3 −
k4 ,
k5 = h f tk + h, yk +
216
513
4104
1835441859
11
k6 = h f tk + h, yk − k1 + 2k2 −k3 +k4 − k5 .
22725654104
40

Then an approximation to the solution of the I.V.P. is made using a Runge-Kutta
method of order 4:(29)yk+1 = yk +251408
2197
1k1 +
k3 +
k4 − k5 ,21625654101
5where the four function values f 1 , f 3 , f 4 , and f 5 are used. Notice that f 2 is not usedin formula (29). A better value for the solution is determined using a Runge-Kutta
method of order 5:(30)z k+1 = yk +16665628,561
9
2
k1 +
k3 +
k4 − k5 + k6 .
135
12,82556,4305055

The optimal step size sh can be determined by multiplying the scalar s times thecurrent step size h. The scalar s is(31)s=tol h2|z k+1 − yk+1 |1/4≈ 0.84tol h|z k+1 − yk+1 |

1/4

498

C HAP. 9S OLUTION OF D IFFERENTIAL E QUATIONS

Table 9.10

RKF45 Solution to y  = 1 + y 2 , y(0) = 0ktk

RK45 approximationyk

True solution,y(tk ) = tan(tk )Errory(tk ) − yk

01
2
3
4
56
78
910

0.0
0.2
0.4
0.6
0.8
1.0
1.1
1.2
1.3
1.35
1.4

0.0000000
0.2027100
0.4227933
0.6841376
1.0296434
1.5574398
1.9648085
2.5722408
3.6023295
4.4555714
5.7985045

0.0000000
0.2027100
0.4227931
0.6841368
1.0296386
1.5774077
1.9647597
2.5721516
3.6021024
4.4552218
5.7978837

0.0000000
0.0000000−0.0000002−0.0000008−0.0000048−0.0000321−0.0000488−0.0000892−0.0002271−0.0003496−0.0006208

where tol is the specified error control tolerance.
The derivation of formula (31) can be found in advanced books on numerical analysis. It is important to learn that a fixed step size is not the best strategy even though
it would give a nicer-appearing table of values. If values are needed that are not in thetable, polynomial interpolation should be used.
Example 9.12.

Compare RKF45 and RK4 solutions to the I.V.P.
y = 1 + y2

withy(0) = 0on[0, 1.4].

An RKF45 program was used with the value tol = 2 × 10−5 for the error control
tolerance. It changed the step size automatically and generated the 10 approximations tothe solution in Table 9.10. An RK4 program was used with the a priori step size of h = 0.1,which required the computer to generate 14 approximations at the equally spaced points inTable 9.11. The approximations at the right endpoint arey(1.4) ≈ y10 = 5.7985045

andy(1.4) ≈ y14 = 5.7919748

andE 14 = 0.0059089

and the errors are
E 10 = −0.0006208

for the RKF45 and RK4 methods, respectively. The RKF45 method has the smallererror.

S EC . 9.5

RUNGE -K UTTA M ETHODS

Table 9.11

499

RK4 Solution to y  = 1 + y 2 , y(0) = 0ktk

RK4 approximationyk

True solution,y(tk ) = tan(tk )Errory(tk ) − yk

0
12
34
5
6
7
8
9
10
11
12
13
14

0.0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1.0
1.1
1.2
1.3
1.4

0.0000000
0.1003346
0.2027099
0.3093360
0.4227930
0.5463023
0.6841368
0.8422886
1.0296391
1.2601588
1.5574064
1.9647466
2.5720718
3.6015634
5.7919748

0.0000000
0.1003347
0.2027100
0.3093362
0.4227932
0.5463025
0.6841368
0.8422884
1.0296386
1.2601582
1.5574077
1.9647597
2.5721516
3.6021024
5.7978837

0.0000000
0.0000001
0.0000001
0.0000002
0.0000002
0.0000002
0.0000000−0.0000002−0.0000005−0.0000006
0.0000013
0.0000131
0.0000798
0.0005390
0.0059089

Numerical Methods Using Matlab, 4th Edition, 2004
John H. Mathews and Kurtis K. Fink
ISBN: 0-13-065248-2
Prentice-Hall Inc.
Upper Saddle River, New Jersey, USA
http://vig.prenhall.com/