Editing Runge–Kutta methods (section)

==Derivation of the Runge–Kutta fourth-order method==
In general a Runge–Kutta method of order <math>s</math> can be written as:
:<math>y_{t + h} = y_t + h \cdot \sum_{i=1}^s a_i k_i +\mathcal{O}(h^{s+1}),</math>
where:
:<math>k_i =\sum_{j = 1}^s \beta_{ij} f(k_j,\ t_n + \alpha_i h)</math>
are increments obtained evaluating the derivatives of <math>y_t</math> at the <math>i</math>-th order.

We develop the derivation<ref>{{cite web |last1=Lyu |first1=Ling-Hsiao |title=Appendix C. Derivation of the Numerical Integration Formulae |url=http://www.ss.ncu.edu.tw/~lyu/lecture_files_en/lyu_NSSP_Notes/Lyu_NSSP_AppendixC.pdf |website=Numerical Simulation of Space Plasmas (I) Lecture Notes |publisher=Institute of Space Science, National Central University |access-date=17 April 2022 |date=August 2016}}</ref> for the Runge–Kutta fourth-order method using the general formula with <math>s=4</math> evaluated, as explained above, at the starting point, the midpoint and the end point of any interval <math>(t,\ t +h)</math>; thus, we choose:

: <math>
\begin{align}
&\alpha_i & &\beta_{ij} \\
\alpha_1 &= 0 & \beta_{21} &= \frac{1}{2} \\
\alpha_2 &= \frac{1}{2} &  \beta_{32} &= \frac{1}{2} \\
\alpha_3 &= \frac{1}{2} &  \beta_{43} &= 1 \\
\alpha_4 &= 1 & &\\
\end{align}
</math>

and <math>\beta_{ij} = 0</math> otherwise. We begin by defining the following quantities:
:<math>\begin{align}
y^1_{t+h} &= y_t + hf\left(y_t,\ t\right) \\
y^2_{t+h} &= y_t + hf\left(y^1_{t+h/2},\ t+\frac{h}{2}\right) \\
y^3_{t+h} &= y_t + hf\left(y^2_{t+h/2},\ t+\frac{h}{2}\right)
\end{align}</math>
where <math>y^1_{t+h/2} = \dfrac{y_t + y^1_{t+h}}{2}</math> and <math>y^2_{t+h/2} = \dfrac{y_t + y^2_{t+h}}{2}. </math>
If we define:
:<math>\begin{align}
k_1 &= f(y_t,\ t) \\
k_2 &= f\left(y^1_{t+h/2},\ t + \frac{h}{2}\right) = f\left(y_t + \frac{h}{2} k_1,\ t + \frac{h}{2}\right) \\
k_3 &= f\left(y^2_{t+h/2},\ t + \frac{h}{2}\right) = f\left(y_t + \frac{h}{2} k_2,\ t + \frac{h}{2}\right) \\
k_4 &= f\left(y^3_{t+h},\ t + h\right) = f\left(y_t + h k_3,\ t + h\right)
\end{align}</math>
and for the previous relations we can show that the following equalities hold up to <math>\mathcal{O}(h^2)</math>:<math display=block>\begin{align}
k_2 &= f\left(y^1_{t+h/2},\ t + \frac{h}{2}\right) = f\left(y_t + \frac{h}{2} k_1,\ t + \frac{h}{2}\right) \\
&= f\left(y_t,\ t\right) + \frac{h}{2} \frac{d}{dt}f\left(y_t,\ t\right) \\
k_3 &= f\left(y^2_{t+h/2},\ t + \frac{h}{2}\right) = f\left(y_t + \frac{h}{2} f\left(y_t + \frac{h}{2} k_1,\ t + \frac{h}{2}\right),\ t + \frac{h}{2}\right) \\
&= f\left(y_t,\ t\right) + \frac{h}{2} \frac{d}{dt} \left[ f\left(y_t,\ t\right) + \frac{h}{2} \frac{d}{dt}f\left(y_t,\ t\right) \right] \\
k_4 &= f\left(y^3_{t+h},\ t + h\right) = f\left(y_t + h f\left(y_t + \frac{h}{2} k_2,\ t + \frac{h}{2}\right),\ t + h\right) \\
&= f\left(y_t + h f\left(y_t + \frac{h}{2} f\left(y_t + \frac{h}{2} f\left(y_t,\ t\right),\ t + \frac{h}{2}\right),\ t + \frac{h}{2}\right),\ t + h\right)  \\
&= f\left(y_t,\ t\right) + h \frac{d}{dt} \left[ f\left(y_t,\ t\right) + \frac{h}{2} \frac{d}{dt}\left[ f\left(y_t,\ t\right) + \frac{h}{2} \frac{d}{dt}f\left(y_t,\ t\right) \right]\right]
\end{align}</math>
where:<math display=block>\frac{d}{dt} f(y_t,\ t) = \frac{\partial}{\partial y} f(y_t,\ t) \dot y_t + \frac{\partial}{\partial t} f(y_t,\ t) = f_y(y_t,\ t) \dot y_t + f_t(y_t,\ t) := \ddot y_t</math>
is the total derivative of <math>f</math> with respect to time.

If we now express the general formula using what we just derived we obtain:<math display=block>\begin{align}
y_{t+h} = {} & y_t + h  \left\lbrace a \cdot f(y_t,\ t) + b \cdot \left[ f(y_t,\ t) + \frac{h}{2} \frac{d}{dt}f(y_t,\ t) \right] \right.+ \\
& {}+ c \cdot \left[ f(y_t,\ t) + \frac{h}{2} \frac{d}{dt} \left[ f\left(y_t,\ t\right) + \frac{h}{2} \frac{d}{dt}f(y_t,\ t) \right] \right] + \\
&{}+ d \cdot \left[f(y_t,\ t) + h \frac{d}{dt} \left[ f(y_t,\ t) + \frac{h}{2} \frac{d}{dt}\left[ f(y_t,\ t)
+ \left. \frac{h}{2} \frac{d}{dt}f(y_t,\ t) \right]\right]\right]\right\rbrace + \mathcal{O}(h^5) \\
= {} & y_t + a \cdot h f_t + b \cdot h f_t + b \cdot \frac{h^2}{2} \frac{df_t}{dt}   + c \cdot h f_t+ c \cdot \frac{h^2}{2} \frac{df_t}{dt} + \\
&{}+ c \cdot \frac{h^3}{4} \frac{d^2f_t}{dt^2} + d \cdot h f_t + d \cdot h^2 \frac{df_t}{dt} + d \cdot \frac{h^3}{2} \frac{d^2f_t}{dt^2} + d \cdot \frac{h^4}{4} \frac{d^3f_t}{dt^3} + \mathcal{O}(h^5)
\end{align}</math>

and comparing this with the [[Taylor series]] of <math>y_{t+h}</math> around <math>t</math>:<math display=block>\begin{align}
 y_{t+h} &= y_t + h \dot y_t + \frac{h^2}{2} \ddot y_t + \frac{h^3}{6} y^{(3)}_t + \frac{h^4}{24} y^{(4)}_t + \mathcal{O}(h^5) = \\
&= y_t + h f(y_t,\ t) + \frac{h^2}{2} \frac{d}{dt}f(y_t,\ t) + \frac{h^3}{6} \frac{d^2}{dt^2}f(y_t,\ t) + \frac{h^4}{24} \frac{d^3}{dt^3}f(y_t,\ t)
\end{align}</math>

we obtain a system of constraints on the coefficients:

:<math>
 \begin{cases}
 & a + b + c + d = 1 \\[6pt]
 & \frac{1}{2} b + \frac{1}{2} c  + d = \frac{1}{2} \\[6pt]
 & \frac{1}{4} c + \frac{1}{2} d = \frac{1}{6} \\[6pt]
 & \frac{1}{4} d = \frac{1}{24}
 \end{cases}</math>
which when solved gives <math>a = \frac{1}{6}, b = \frac{1}{3}, c = \frac{1}{3}, d = \frac{1}{6}</math> as stated above.