请详细推导四阶的龙格库塔

四阶龙格-库塔法是一种数值解常微分方程的方法，它的推导过程如下：

考虑求解一阶常微分方程

$$
\frac{dy}{dt} = f(t,y),\quad y(t_0) = y_0
$$

将时间 $[t_n, t_{n+1}]$ 均分成 $N$ 个子区间，每个子区间的长度为 $h = (t_{n+1} - t_n)/N$。则 $y_{n+1}$ 可以通过 $y_n$ 和 $f(t_n, y_n)$ 进行逼近，我们可以使用以下的龙格-库塔公式：

$$
\begin{aligned}
k_1 &= f(t_n, y_n) \\
k_2 &= f(t_n + h/2, y_n + h/2 k_1) \\
k_3 &= f(t_n + h/2, y_n + h/2 k_2) \\
k_4 &= f(t_n + h, y_n + h k_3) \\
y_{n+1} &= y_n + \frac{h}{6}(k_1 + 2k_2 + 2k_3 + k_4)
\end{aligned}
$$

其中，$k_1$ 是在 $t_n$ 处的斜率，$k_2$ 是在 $t_n + h/2$ 处的斜率，$k_3$ 是在 $t_n + h/2$ 处的斜率，$k_4$ 是在 $t_n + h$ 处的斜率。这个公式就是四阶龙格-库塔法的表达式。

该公式的推导依赖于泰勒展开式，将 $y(t_{n+1})$ 在 $y(t_n)$ 处展开：

$$
\begin{aligned}
y(t_{n+1}) &= y(t_n) + h\frac{dy}{dt}\bigg|_{t_n} + \frac{1}{2}h^2\frac{d^2y}{dt^2}\bigg|_{t_n} + \frac{1}{6}h^3\frac{d^3y}{dt^3}\bigg|_{t_n} + \frac{1}{24}h^4\frac{d^4y}{dt^4}\bigg|_{t_n} + O(h^5) \\
&= y(t_n) + hf(t_n, y_n) + \frac{1}{2}h^2\frac{d}{dt}f(t_n, y_n) + \frac{1}{6}h^3\frac{d^2}{dt^2}f(t_n, y_n) + \frac{1}{24}h^4\frac{d^3}{dt^3}f(t_n, y_n) + O(h^5)
\end{aligned}
$$

然后，我们使用中点法近似第二项和第三项，四阶龙格-库塔法的表达式就得到了。