<s>[INST] Gain matrices of numerical schemes 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?
 [/INST] 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.
</s>