# Using Runge-Kutta-Fehlberg 4-5 for higher dimension systems

When applying RKF45 algorithm to a first order ODE with higher dimensions, e.g. $f(t,y_1,y_2)$ and $f(t,y_1,y_2,y_3)$, is the procedure simply a matter of applying RKF45 to each dimension in turn? I.e. $f(t,y_1;y_2,y_3)$, $f(t,y_2;y_1,y_3)$, and $f(t,y_3;y_1,y_2)$, where constants appear after the semi-colon, would give me the corresponding component values. Or is there more to it than this, perhaps a specific algorithm for such cases?

RK4 in 3 functions/variables

After doing a bit more searching and thinking I found this. Consider just the the RK-4 algorithm in 3 functions $y_1,y_2,y_3$ for now. Let

$$\left.\begin{array}{l}y_1'(t)=f_1(t,y_1,y_2,y_3)\\y_2'(t)=f_2(t,y_1,y_2,y_3)\\y_3'(t)=f_3(t,y_1,y_2,y_3)\end{array}\right\}.$$

Given we know $\Phi_n=\left(y_1^{(n)},y_2^{(n)},y_3^{(n)}\right)$, e.g. the initial value $\Phi_0$, we wish to find $\Phi_{n+1}=\left(y_1^{(n+1)},y_2^{(n+1)},y_3^{(n+1)}\right)$. We have

$$\begin{array}{l}a_j^{(n)} = f_j\left(y_1^{(n)},y_2^{(n)},y_3^{(n)}\right)\\ b_j^{(n)}=f_j\left(y_1^{(n)}+\frac{h}{2}a_1^{(n)},y_2^{(n)}+\frac{h}{2}a_2^{(n)},y_3^{(n)}+\frac{h}{2}a_3^{(n)}\right)\\ c_j^{(n)}=f_j\left(y_1^{(n)}+\frac{h}{2}b_1^{(n)},y_2^{(n)}+\frac{h}{2}b_2^{(n)},y_3^{(n)}+\frac{h}{2}b_3^{(n)}\right)\\ d_j^{(n)}=f_j\left(y_1^{(n)}+hc_1^{(n)},y_2^{(n)}+hc_2^{(n)},y_3^{(n)}+hc_3^{(n)}\right)\\ \end{array},$$

and then

$$y_j^{(n+1)}=y_j^{(n)}+\frac{h}{6}\left(a_j^{(n)}+2b_j^{(n)}+2c_j^{(n)}+d_j^{(n)}\right),$$ where $j\in\{1,2,3\}$.

Is this the right way ?

If so I could adapt to RK45 quite easily I think, e.g. adapt to RK5, then use the following in the code below instead of R=abs(y1-y2) / h:

$$R=\frac{1}{h}\left\|\mathbf{y}_{RK4}^{(n+1)}-\mathbf{y}_{RK5}^{(n+1)}\right\|_2.$$

RKF45 code for 1 function

Just to show what I did originally, this is my implementation for the 1 function case:

    y = initial_y0;

print("step " + step + ", t = " + t + ", w = " + y);

while(t < t1)
{
h = min(h, t1 - t);

k1 = h * f(t, y);
k2 = h * f(t + h/4, y + k1/4);
k3 = h * f(t + 3*h/8, y + 3*k1/32 + 9*k2/32);
k4 = h * f(t + 12*h/13, y + 1932*k1/2197 - 7200*k2/2197 + 7296*k3/2197);
k5 = h * f(t + h, y + 439*k1/216 - 8*k2 + 3680*k3/513 - 845*k4/4104);
k6 = h * f(t + h/2, y - 8*k1/27 + 2*k2 - 3544*k3/2565 + 1859*k4/4104 - 11*k5/40);

y1 = y + 25*k1/216 + 1408*k3/2565 + 2197*k4/4104 - k5/5;
y2 = y + 16*k1/135 + 6656*k3/12825 + 28561*k4/56430 - 9*k5/50 + 2*k6/55;

R = abs(y1 - y2) / h;
delta = 0.84  * power(epsilon/R, 0.25);

if(R <= epsilon)        // if error < required
{
t += h;             // next t value
y = y1;             // use the RK4 approx

print("step " + step + ", t = " + t + ", w = " + y);
step = step + 1
}

h = h * delta;          // adapt the step size
}


No. Adaptng the notation from the Wikipedia Runge-Kutta article, you have $y^n_{ijk}$ and you want $y^{n+1}_{ijk}$. For each $ijk$ you have to construct a $k_{ijk}$. Let's say we're doing the diffusion equation with central differencing, diffusion coeff $\alpha$. Then

$\displaystyle k_{ijk} = \alpha\times (y^n_{i+1jk}+y^n_{i-1jk}+y^n_{ij+1k}+y^n_{ij-1k}+ y^n_{ijk+1}+y^n_{ijk-1}-6y^n_{ijk})/\Delta x^2$

with maybe an overall minus sign, depending on how you define things.

The way I think of it is, I take my multidimensional grid and arrange all my grid points in a single vector. Then, for each grid point I ask, "What is the $k$ value for this point?"

• Thanks. Do you know of a reference which contains RK45 in 2 and 3 dimension spaces? – Pixel Nov 17 '15 at 20:55
• No, I had to "pick it up off the streets". :) But, you can write down the algorithm and take a look at each "k" and just ask yourself, "OK, what should go where this 'k' is?" – bob.sacamento Nov 17 '15 at 21:00
• I've updated my question - but I'm not clear if it relates to your answer. Does my update agree with your answer? Perhaps your $k_{ijk}$ are my $a,b,c,d$, e.g. $k_{1jn}=a_j^{(n)}$ – Pixel Nov 18 '15 at 9:02
• Yes, I think you're getting it. My $k$ does correspond to $a$, $b$, etc. If any of the $f$'s involve derivatives, you have to be a little careful, but it is conceptually the same thing. Regarding error, it becomes as much an art as a science. You have to make some decisions about what kind of error you are interested in. Do you want the error to be the maximum error of all the $y$'s? Do you want it to be the sum of all the $y$ errors, or the Euclidean "distance"? There's no blanket prescription for this. It depends on the problem you are working. – bob.sacamento Nov 18 '15 at 14:19

In term of programming you need to apply each step to all the variables before you move to the next step.

Consider an array function double[] f(double t, double[] y) for the slopes. You algorithm should follow these general steps I outline below in C#

public Program()
{
double t=0, h=0.2;
double[] y=new double[] { 10, 0 };
do
{
y=Rk4Step(t, y, h);
t+=h;
// y contains the solution for the step
} while(t<10);
}

public double[] f(double t, double[] y)
{
return new double[] { y[1], -10*y[0]-1*y[1] };
}

public double[] Rk4Step(double t, double[] y, double h)
{
int n=y.Length;
double[] yn=new double[n];

double[] K0=f(t, y);
for(int i=0; i<n; i++)
{
yn[i]=y[i]+(h/2)*K0[i];
}
double[] K1=f(t+h/2, yn);
for(int i=0; i<n; i++)
{
yn[i]=y[i]+(h/2)*K1[i];
}
double[] K2=f(t+h/2, yn);
for(int i=0; i<n; i++)
{
yn[i]=y[i]+h*K2[i];
}
double[] K3=f(t+h, yn);

for(int i=0; i<n; i++)
{
yn[i]=y[i]+(h/6)*(K0[i]+2*K1[i]+2*K2[i]+K3[i]);
}

return yn;
}


Of course you need the proper coefficients for RK45.

• Thanks - I missed that - the $k_i$ should be computed for all components before moving onto $k_{i+1}$ because to compute the $k_{i+1}$ for all components we need $k_i$ for all components. – Pixel Nov 19 '15 at 10:02
• I thought so because that is what I missed the first time I did multivariate RK – ja72 Nov 19 '15 at 12:56

Ok, I think I got it now. Based on articles on the internet, comments from other answers, and the little grey cells, I think the following is correct. But please - comments still welcomed !

$$a_j^{(n)}=hf_j(t,y_1,y_2,y_3).$$ $$b_j^{(n)}=hf_j\left(t+\frac{h}{4},y_1^{(n)}+\frac{1}{4}a_1^{(n)} ,y_"^{(n)}+\frac{1}{4}a_2^{(n)},y_3^{(n)}+\frac{1}{4}a_3^{(n)}\right).$$ $$c_j^{(n)}=hf_j\left(t+\frac{3}{8}h,y_1^{(n)}+\frac{3}{32}a_1^{(n)}+\frac{9}{32}b_1^{(n)},y_2^{(n)}+\frac{3}{32}a_2^{(n)}+\frac{9}{32}b_2^{(n)},y_3^{(n)}+\frac{3}{32}a_3^{(n)}+\frac{9}{32}b_3^{(n)},\right).$$ $$d_j^{(n)}=hf_j\left(t+\frac{12}{13}h,y_1^{(n)}+\frac{1932}{2197}a_1^{(n)}-\frac{7200}{2197}b_1^{(n)}+\frac{7296}{2197}c_1^{(n)},y_2^{(n)}+\frac{1932}{2197}a_2^{(n)}-\frac{7200}{2197}b_2^{(n)}+\frac{7296}{2197}c_2^{(n)},y_3^{(n)}+\frac{1932}{2197}a_3^{(n)}-\frac{7200}{2197}b_3^{(n)}+\frac{7296}{2197}c_3^{(n)}\right).$$ $$e_j^{(n)}=hf_j\left(t+h,y_1^{(n)}+\frac{439}{216}a_1^{(n)}-8b_1^{(n)}+\frac{3680}{513}c_1^{(n)}-\frac{845}{4104}d_1^{(n)},y_2^{(n)}+\frac{439}{216}a_2^{(n)}-8b_2^{(n)}+\frac{3680}{513}c_2^{(n)}-\frac{845}{4104}d_2^{(n)},y_3^{(n)}+\frac{439}{216}a_3^{(n)}-8b_3^{(n)}+\frac{3680}{513}c_3^{(n)}-\frac{845}{4104}d_3^{(n)}\right).$$ $$f_j^{(n)}=f_j\left(t+\frac{h}{2},y_1^{(n)}-\frac{8}{27}a_1^{(n)}+2b_1^{(n)}-\frac{3544}{2565}c_1^{(n)}+\frac{1859}{4104}d_1^{(n)}-\frac{11}{40}e_1^{(n)},y_2^{(n)}-\frac{8}{27}a_2^{(n)}+2b_2^{(n)}-\frac{3544}{2565}c_2^{(n)}+\frac{1859}{4104}d_2^{(n)}-\frac{11}{40}e_2^{(n)},y_3^{(n)}-\frac{8}{27}a_3^{(n)}+2b_3^{(n)}-\frac{3544}{2565}c_3^{(n)}+\frac{1859}{4104}d_3^{(n)}-\frac{11}{40}e_3^{(n)}\right).$$

Then we compute the RK4 and RK5 values $\bar{y}$ and $\hat{y}$ respectively:

$$\bar{y}_j^{(n+1)}=y_j^{(n)}+\frac{25}{216}a_j^{(n)}+\frac{1408}{2565}c_j^{(n)}+\frac{2197}{4101}d_j^{(n)}-\frac{1}{5}e_j^{(n)}.$$

$$\hat{y}_j^{(n+1)}=y_j^{(n)}+\frac{16}{135}a_j^{(n)}+\frac{6656}{12825}c_j^{(n)}+\frac{28561}{56430}d_j^{(n)}-\frac{9}{50}e_j^{(n)}+\frac{2}{55}f_j^{(n)}.$$

Next we compute the error,

$$R=\frac{1}{h}\left\|\left(\bar{y}_j^{(n+1)}\right)_{j=1}^3-\left(\hat{y}_j^{(n+1)}\right)_{j=1}^3\right\|,$$ where $\|\cdot\|$ is some suitable norm.

Using $R$ we may decide that the tolerance is good in which case we let $y_j^{(n+1)}=\bar{y}_j^{(n+1)}$, and increment $t$, otherwise we reduce the step size $h$. Repeat until $t\geq t_b$, at which point we stop.

Here's the code for RK45 in 3 dimensions:

    while(t < t1)
{
h = min(h, t1 - t);
diffmax = 0;

for(int j = 0; j < 3; j++)
k1[j] = h * fj(j, t, y[0], y[1], y[2]);

for(int j = 0; j < 3; j++)
k2[j] = h * fj(j, t + h/4, y[0] + k1[0]/4,
y[1] + k1[1]/4,
y[2] + k1[2]/4);

for(int j = 0; j < 3; j++)
k3[j] = h * fj(j, t + 3*h/8, y[0] + 3*k1[0]/32 + 9*k2[0]/32,
y[1] + 3*k1[1]/32 + 9*k2[1]/32,
y[2] + 3*k1[2]/32 + 9*k2[2]/32);

for(int j = 0; j < 3; j++)
k4[j] = h * fj(j, t + 12*h/13, y[0] + 1932*k1[0]/2197 - 7200*k2[0]/2197 + 7296*k3[0]/2197,
y[1] + 1932*k1[1]/2197 - 7200*k2[1]/2197 + 7296*k3[1]/2197,
y[2] + 1932*k1[2]/2197 - 7200*k2[2]/2197 + 7296*k3[2]/2197);

for(int j = 0; j < 3; j++)
k5[j] = h * fj(j, t + h, y[0] + 439*k1[0]/216 - 8*k2[0] + 3680*k3[0]/513 - 845*k4[0]/4104,
y[1] + 439*k1[1]/216 - 8*k2[1] + 3680*k3[1]/513 - 845*k4[1]/4104,
y[2] + 439*k1[2]/216 - 8*k2[2] + 3680*k3[2]/513 - 845*k4[2]/4104);

for(int j = 0; j < 3; j++)
k6[j] = h * fj(j, t + h/2, y[0] - 8*k1[0]/27 + 2*k2[0] - 3544*k3[0]/2565 + 1859*k4[0]/4104 - 11*k5[0]/40,
y[1] - 8*k1[1]/27 + 2*k2[1] - 3544*k3[1]/2565 + 1859*k4[1]/4104 - 11*k5[1]/40,
y[2] - 8*k1[2]/27 + 2*k2[2] - 3544*k3[2]/2565 + 1859*k4[2]/4104 - 11*k5[2]/40);

for(int j = 0; j < 3; j++)
{
y_rk4[j] = y[j] + 25*k1[j]/216 + 1408*k3[j]/2565 + 2197*k4[j]/4104 - k5[j]/5;                       // compute RK4
y_rk5 = y[j] + 16*k1[j]/135 + 6656*k3[j]/12825 + 28561*k4[j]/56430 - 9*k5[j]/50 + 2*k6[j]/55;       // compute RK5

// choose the largest error - maybe this bit could be improved...
if(abs(y_rk5-y_rk4[j]) > diffmax)
diffmax = Math.abs(y_rk5-y_rk4[j]);
}

// compute error
R = diffmax / h;
delta = 0.84  * pow(epsilon / R, 0.25);

if(R <= epsilon)
{
t += h;
y[0] = y_rk4[0];
y[1] = y_rk4[1];
y[2] = y_rk4[2];

print("step " + step + ", t = " + t + ", y = (" + y[0] + ", " + y[1] + ", " + y[2] + ")");
step = step + 1;
}

h = h * delta;
}


I've tested it against Mathematica, e.g.

sol = NDSolve[{y1'[t] == y2[t]*y3[t] + t*y1[t] - y3[t] + 1, y2'[t] == t*y2[t] + y2[t] + 1, y3'[t] == y1[t] t + t*y3[t] - 3, y1[0] == 0.5, y2[0] == 0.8, y3[0] == 0.2}, {y1, y2, y3}, {t, 2}]

{y1[2] /. sol, y2[2] /. sol, y3[2] /. sol}


Works very well !

What usually is done that you assemble the y123 into a vector y=[y1,y2,y3] and have the ODE function F return the vector of the scalar results, F(t,y)=[ f1(t,y), f2(t,y), f3(t,y) ]. After all operations are translated into vector arithmetic, the vector case does not look (much) different from the scalar case.