Introduction
Say I have a linear multiple degree of freedom system, as below: \begin{equation*} \frac{dx}{dt}=ax+by \end{equation*} \begin{equation*} \frac{dy}{dt}=dx+ey \end{equation*} This can also be written as $\frac{dr}{dt}=Ar$ with $A$ given by \begin{equation*} A= \begin{bmatrix} a & b \\ c & d \\ \end{bmatrix} \end{equation*} Matrix $A$ can be used to write the numerical schemes (Euler, Heun and Runge-Kutta) using gain matrices.
For Euler is this \begin{equation*} r_{n+1}=r_{n}+\frac{dr}{dt}\Delta t=r_{n}+Ar_{n}\Delta t \end{equation*} and rewriting gives \begin{equation*} r_{n+1}=\left(I+A\Delta t\right)r_{n}. \end{equation*} Now we can multiply old vector $r_{n}$ with $G=(I+A\Delta t)$ to get the new vector $r_{n+1}$.
In a similar way, the gain matrices $(G)$ can be determined for the Heun and Runge-Kutta schemes.
For the Heun scheme follows \begin{equation*} \begin{split} r_{n+1}&=r_{n}+\frac{1}{2}\left(\frac{dr}{dt}\rvert_{(r_{n},t_{n})}+\frac{dr}{dt}\rvert_{(r_{p},t_{n+1})}\right)\Delta t \\ &=r_{n}+\frac{1}{2}\left(Ar_{n}+Ar_{p}\right) \\ &=r_{n}+\frac{1}{2}Ar_{n}\Delta t+\frac{1}{2}Ar_{p}\Delta t \\ &=\left(I+\frac{1}{2}A\Delta t\right)r_{n}+\left(\frac{1}{2}A\Delta t\right)r_{p} \\ \end{split} \end{equation*}
Question
Runge Kutta is defined as \begin{equation*} r_{n+1}=r_{n}+\frac{1}{6}(k_{1}+2k_{2}+2k_{3}+k_{4})\Delta t \end{equation*} with \begin{align*} k_{1}&=f\left(r_{n},t_{n}\right) \\ k_{2}&=f\left(r_{n}+\frac{1}{2}k_{1}\Delta t,t_{n}+\frac{1}{2}\Delta t\right) \\ k_{3}&=f\left(r_{n}+\frac{1}{2}k_{2}\Delta t,t_{n}+\frac{1}{2}\Delta t\right) \\ k_{4}&=f\left(r_{n}+k_{3}\Delta t,t_{n}+\Delta t\right) \end{align*}
For the Runge-Kutta scheme I can't figure out how to get gain matrices. Somebody who knows how to do this?
You just insert the $f$ of the homogeneous linear equation, \begin{alignat}{3} k_1&=f(y)&&&&=Ay \\ k_2&=f(y+0.5Δt·k_1)&&=A(y+0.5Δt·k_1))&&=Ay+0.5Δt·A^2y \\ k_3&=f(y+0.5Δt·k_3)&&=A(y+0.5Δt·k_3))&&=Ay+0.5Δt·A^2y+0.25Δt^2·A^3y \\ k_4&=f(y+Δt·k_4)&&=A(y+Δt·k_4)&&=Ay+Δt·A^2y+0.5Δt^2·A^3y+0.25Δt^3·A^4y \\\hline &k_1+2k_2+2k_3+k_4&&&&=6Ay+3Δt·A^2y+Δt^2·A^3y+0.25Δt^3·A^4y \end{alignat} which ends up giving the exact 4th order Taylor polynomial for $e^{Δt⋅A}y$ in the RK4 step.