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
}
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:
I've tested it against Mathematica, e.g.
Works very well !