4

I'm currently looking at the implementation of the Dormand-Prince 5(4) Runge-Kutta algorithm (also known as Dopr5, or ode45) in Numerical Recipes Chapter 17. I understand the integration part, but I have trouble with the dense output function, as it is poorly documented.

Equation 17.2.16 in the book, says that I can interpolate between integration steps ($0\le\theta\le1$), using the following formula:

$$ y(x_n+\theta h)=y_n+b_1(\theta)k_1+b_2(\theta)k_2+b_3(\theta)k_3+b_4(\theta)k_4+b_5(\theta)k_5+b_6(\theta)k_6 $$

However, looking at the implementation, I find that Numerical Recipes actually uses the following equation:

$$ y(x_n+\theta h)=r_1+\theta(r_2+(1-\theta)(r_3+\theta(r_4+(1-\theta)r_5))) $$ Or expanded... $$ y(x_n+\theta h)=r_1+\theta(r_2+r_3)-\theta^2(r_3-r_4-r_5)-\theta^3(r_4+2r_5)+\theta^4r_5 $$

With: $$ r_1=y_n $$ $$ r_2=y_{n+1}-y_n $$ $$ r_3=y_n+h f_n-y_{n+1} $$ $$ r_4=2(y_{n+1}-y_n)-h(f_n+f_{n+1}) $$ $$ r_5=d_1hf_n+d_3k_3+d_4k_4+d_5k_5+d_6k_6+d_7hf_{n+1} $$ And: $$ d_1=-\frac{12715105075}{11282082432} $$ $$ d_3=\frac{87487479700}{32700410799} $$ $$ d_4=-\frac{10690763975}{1880347072} $$ $$ d_5=\frac{701980252875}{199316789632} $$ $$ d_6=-\frac{1453857185}{822651844} $$ $$ d_7=\frac{69997945}{29380423} $$ My question is: How do I get the values for the coefficients $d_1$,...,$d_7$?

Numerical Recipes cites Shampine (1986) and Dormand+Prince (1986), however none of these two derives the coefficients found in the code.

Dormand+Prince (1986) at least derive the polynomials used for equation 17.2.16 (see their equation 25), and their polynomials share the same denominator with Numerical Recipes. However, I have no idea how to derive the coefficients from them.

I tried expanding the formulas using the definition of $y_{n+1}$, and the fact, that $hf_n=k_1$. However, the polynomials I get are different to the ones given by Dormand+Prince (1986).

I would appreciate any help solving this conundrum. Thanks!

jan.sende
  • 193

1 Answers1

2

Okay, I played around with the formulas, and by chance found the translation formula for the coefficients...

One needs to calculate $y(x_n+\frac{1}{2}h)$ of the expanded implementation, using the fact that: $hf_n=k_1$, $hf_{n+1}=k_7$, plus the definition $y_{n+1}=y_n+b_1k_1+b_3k_3+b_4k_4+b_5k_5+b_6k_6$ (for ode45, $b_2=0$).

If one reorders the result in terms of $k$s, one gets:

$$ y(x_n+\frac{1}{2}h)=y_n+\left(\frac{1}{8}+\frac{b_1}{2}+\frac{d_1}{16}\right)k_1+\left(\frac{b_3}{2}+\frac{d_3}{16}\right)k_3+\left(\frac{b_4}{2}+\frac{d_4}{16}\right)k_4+\left(\frac{b_5}{2}+\frac{d_5}{16}\right)k_5+\left(\frac{b_6}{2}+\frac{d_6}{16}\right)k_6+\left(-\frac{1}{8}+\frac{d_7}{16}\right)k_7 $$

Each of these factors is the same as the $c^*$ given in Shampine (1986) (note, that he factored out a $\frac{1}{2}$). Thus, we can derive the following translation:

$$ d_1=-2(1+4b_1-4c^*_1) $$ $$ d_3=-8(b_3-c^*_3) $$ $$ d_4=-8(b_4-c^*_4) $$ $$ d_5=-8(b_5-c^*_5) $$ $$ d_6=-8(b_6-c^*_6) $$ $$ d_7=2(1+4c^*_7) $$ Inserting the numerical values, one can show, that the $d$s used in Numerical Recipes are right.

However, I could not find any justification on why to use the equations Numerical Recipes uses, or how to derive them from scratch. (I tried reversing the steps given above, but the procedure seems arbitrary.) I guess, they try to reduce the number of evaluations, but in modern hardware that should not matter anymore, as caching is way more important.

Sidenote 1:

The coefficients given in Shampine (1986) can be derived from Dormand+Prince (1986), by inserting $\frac{1}{2}$ into the given polynomials for $b^*$ (the same values are denoted with different symbols in both articles).

Sidenote 2:

Both the procedures in Shampine (1986) and Numerical Recipes interpolate by constructing a quartic function over the interval $0\le\theta\le 1$ from the following:

$$ y(x_n+0h)=y_{n} $$ $$ y(x_n+\frac{1/2}h)=y_{n}+c^*_1k_1+c^*_3k_3+c^*_4k_4+c^*_5k_5+c^*_6k_6+c^*_7k_7 $$ $$ y(x_n+1h)=y_{n+1} $$ $$ y'(x_n+0h)=f_{n} $$ $$ y'(x_n+1h)=f_{n+1} $$

jan.sende
  • 193