Adaptive RK Methods

We implemented the Runge-Kutta-Fehlberg adaptive method

   0   |
  1/4  |    1/4
  3/8  |    3/32       9/32
 12/13 | 1932/219  -7200/2197  7296/2197
   1   |  439/216      -8      3680/513     -845/4104
  1/2  |   -8/27        2     -3544/2565    1859/4104  -11/40
 ------------------------------------------------------------------
       |   16/135       0      6656/12825  28561/56430  -9/50  2/55
       |   25/216       0      1408/2565    2197/4104   -1/5    0

to solve the differential equation

y' = f(t,y) where y − y sin(yt)

with initial condition y(0)=1 on the interval [0,3].

Checking the order of the Method

The order of the first method is 5 and the order of the second method is 4. The following Reduce script verifies the order.

load_package "dfpart";
load_package "taylor";
operator y;
generic_function f(t,y);
for all x let df(y(x),x)=f(x,y(x));
k1:=f(t,y(t));
k2:=f(t+h/4,y(t)+h/4*k1);
k3:=f(t+3*h/8,y(t)+3*h/32*k1+9*h/32*k2);
k4:=f(t+12*h/13,y(t)+1932*h/2197*k1-7200*h/2197*k2+7296*h/2197*k3);
k5:=f(t+h,y(t)+439*h/216*k1-8*h*k2+3680*h/513*k3-845*h/4104*k4);
k6:=f(t+h/2,y(t)-8*h/27*k1+2*h*k2-3544*h/2565*k3+1859*h/4104*k4-11*h/40*k5);
y1:=y(t)+h*(16*k1/135+6656*k3/12825+28561*k4/56430-9*k5/50+2*k6/55);
y2:=y(t)+h*(25*k1/216+1408*k3/2565+2197*k4/4104-k5/5);
taylor(y(t+h)-y1,h,0,5);
taylor(y(t+h)-y2,h,0,4);
quit;

The last three lines of output should have been

14: taylor(y(t+h)-y1,h,0,5);

       6
0 + O(h )

15: taylor(y(t+h)-y2,h,0,4);

       5
0 + O(h )

16: quit;

however, the script got stuck because the O(h⁶) Taylor's expansion for y(t+h)−y1 took too long to evaluate. Rerunning the script after class resulted in an execution time of 10 minutes 54 seconds. In class we, at least, checked that y1 was correct to O(h⁵) using

taylor(y(t+h)-y1,h,0,4);

which provided a partial verification of our typing.

The Adaptive Integrator

We wrote the program fehlberg.c to implement adaptive time steps. The first verson of the code only decreased the size of the step if it was too big. The output looks like

0.3 1.28015
0.6 1.45258
0.9 1.50166
1.2 1.50553
1.5 1.56443
1.8 1.91984
# step too big--halve and try again!
1.95 2.49396
2.1 3.10106
2.25 3.27486
2.4 3.29115
2.55 3.3212
2.7 3.64521
# step too big--halve and try again!
# step too big--halve and try again!
2.7375 3.88399
2.775 4.17937
2.8125 4.41972
2.85 4.56037
2.8875 4.63164
2.925 4.66469
2.9625 4.67768
3 4.68094

The improved version rkadapt.c also increases the size of the step if is too small. Note that if we increase the time step, we also need to check that the next time step doesn't go beyond the end time T. The output now looks like

0.3 1.28015
# step too small--doubling next one!
0.9 1.50202
1.5 1.56518
# step too small--doubling next one!
# step too big--halve and try again!
2.1 3.14774
# step too big--halve and try again!
# step too big--halve and try again!
2.25 3.29843
2.4 3.31041
2.55 3.35044
# step too small--doubling next one!
# step too big--halve and try again!
2.7 3.75055
# step too big--halve and try again!
2.775 4.31224
2.85 4.61962
2.925 4.69585
3 4.70526

A more accurate answer can be obtained by setting

double error=0.000001;

at the top of the program to obtain y(3)≈4.68228. Note that both of our original approximations were within 0.01 of this more accurate answer.

Last Updated: Fri Nov 21 13:26:45 PST 2014