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
}
Best Answer
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?"