# 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

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.

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 – John Alexiou 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 !