一阶初值问题

设 $f(y,t)$ 在区域 $G:0\le t\le T,|y|<\infty$ 上连续，对于一个给定的常微分方程 $\frac{dy}{dt}=f(y,t)$ 及初值 $y(0)=y_0$ ，求解 $y(t)$ ，即
::: align-center
$\left\{\begin{matrix} \frac{dy}{dt}=f(y,t)\\ y(0)=y_0 \end{matrix}\right .$

:::
为了保证解 $y(t)$ 存在、唯一且连续依赖初值 $y_0$ ，要求 $f(y,t)$ 需要满足Lipschitz条件，即
::: align-center
存在常数 $L$ ，使得 $|f(y_1,t)-f(y_2,t)|\le L|y_1-y_2|$ 对所有 $t\in[0,T]$ 和 $y_1,y_2\in(-\infty,\infty)$ 成立

:::
假设 $f(y,t)$ 总满足以上条件，则此时微分方程有近似解法，常用的有欧拉法、改进的欧拉法与龙格—库塔法

欧拉法

将区间 $[0,T]$ 作N等分，每一小区间长度 $\tau=\frac{T}{N}$ 称为步长， $t_n=n\tau$ 称为节点，根据初值 $y(0)=y_0$ ，代入微分方程可直接解出 $t=t_0$ 的导数值 $y'=f(y_0,t_0)$ ，之后再进行迭代最终求得 $y$
::: align-center
$\begin{aligned} &t_{k+1}=t_k+\tau\\ &y_{k+1}=y_k+\tau f(y_k,t_k) \end{aligned}$

:::

推导

由泰勒公式，将 $y(t_1)=y(t_0+\tau)$ 在 $t_0$ 处展开
$y(t_1)=y(t_0)+y'(t_0)\tau+O(\tau^2)$
略去二阶无穷小量，得
$y(t_0)+y'(t_0)\tau=y_0+f(y_0,t_0)\tau\approx y_1$
类似递推即可得

此外还有数值积分推导法，将在下面介绍到

改进的欧拉法（又称预估—校正法）

顾名思义，使用此方法时

使用欧拉法计算出 $\hat{y_{k+1}}$ 进行预估
后再用 $\hat{y_{k+1}}$ 求出真正的近似值进行校正
::: align-center
$\begin{aligned} &\hat{y_{k+1}}=y_k+\tau f(y_k,t_k)\\ &y_{k+1}=y_k+\frac{\tau}{2}[f(y_k,t_k)+f(\hat{y_{k+1}},t_{k+1})] \end{aligned}$

:::

推导

此时我们使用数值积分方法进行推导
在某个区间上， $y'$ 的实际增量可理解为其下方的面积
左梯形积分.jpg
欧拉法计算的是以 $f(y_k,t_k)$ （图中为 $f(x_k,y_k)$ 为高的矩形面积），而改进的欧拉法计算的是图中所示的改进的左梯形数值积分的面积
容易知道，在步长合理时，显然左梯形的面积更接近实际函数下方面积，故改进的欧拉法实际上效果好于欧拉法，但计算较为复杂

龙格-库塔法（又称R-K法）

我们先回到欧拉法，实际上欧拉法的基本思想是使用差商代替导数，实际上，我们可以考虑
::: align-center
$\frac{y(t_{n+1})-y(t_n)}{\tau}$

:::
根据Lagrange中值定理，我们可得到
::: align-center
$\frac{y(t_{n+1})-y(t_n)}{\tau}=y'(t_n+\theta \tau),0<\theta<1$

:::
又注意到一阶初值问题中有 $y'=f(y,t)$ ，于是有
::: align-center
$y(t_{n+1})=y(t_n)+\tau f(t_n+\theta\tau,y(t_n+\theta\tau))$

:::
不妨记 $\overline{K}=f(t_n+\theta\tau,y(t_n+\theta\tau))$ ，称为区间 $[t_n,t_{n+1}]$ 上的平均斜率，可见给出一种平均斜率的计算方法，我们就能导出一种算法

在如上的欧拉公式中 $\theta=0$ ，即 $\overline{K}=f(t_n,y_n)$ ，而改进的欧拉公式取 $\overline{K}$ 为 $f(t_n.y_n)$ 和 $f(t_{n+1},\hat{y_{n+1}})$ 的平均值，而这种操作提高了精度

所以，我们不妨在区间 $[t_n,t_{n+1}]$ 内多取几个点，而将它们的加权平均作为 $\overline{K}$ ，这就有可能构造出精度更高的算法，而这就是龙格-库塔法的基本思想

二阶R-K法

首先不妨在区间内取两个点，使用
::: align-center
$\left\{\begin{aligned} &y_{n+1}=y_{n}+\tau(\lambda_1k_1+\lambda_2k_2)\\ &k_1=f(t_n,y_n)\\ &k_2=f(t_n+\alpha\tau,y_n+\beta \tau k_1) \end{aligned}\right .$

:::
其中 $\lambda_1,\lambda_2,\alpha,\beta$ 为待定系数，我们分析如何选这几个参数可以使误差最低，为此我们分析局部阶段误差(local truncation error) $y(t_{n+1})-y_{n+1}$ ，因为迭代时我们假设 $y_n=y(t_n)$ ，故上式可化为（*）
::: align-center
$\left\{\begin{aligned} &y_{n+1}=y(t_n)+\tau(\lambda_1k_1+\lambda_2k_2)\\ &k_1=f(t_n,y(t_n))=y'(t_n)\\ &k_2=f(t_n+\alpha\tau,y_n+\beta \tau k_1)\\ \space&=f(t_n,y(t_n))+\alpha\tau f_x(t_n,y(t_n))+\beta hk_1f_y(t_n,y(t_n)) +O(h^2) \end{aligned}\right .$

:::
其中 $k_2$ 的展开使用了多元函数的泰勒公式
将 $k_1,k_2$ 代入 $y_{n+1}$ 得
::: align-center
$y_{n+1}=y(t_n)+(\lambda_1+\lambda_2)\tau y'(t_n)+\lambda_2\alpha \tau^2(f_t+\frac{\beta}{\alpha}ff_y)+O(h^3)$

:::
注意到 $y(t_{n+1})$ 的泰勒展开式中
::: align-center
$y(t_{n+1})=y(t_n)+\tau y'(t_n)+\frac{\tau^2}{2}y''(t_n)+O(h^3)$

:::
的 $y'=f$ ，所以对其求全微分得 $y''=f_t+ff_y$ ，如此，两式对应，要使误差为 $y(t_{n+1})-y(t_n)=O(h^3)$ ，只需令
::: align-center
$\lambda_1+\lambda_2=1,\lambda_2\alpha=\frac{1}{2},\frac{\beta}{\alpha}=1$

:::
待定系数满足上式的（*）称为2阶龙格-库塔公式，由于有四个未知数三个方程，故四个参数不唯一，特别地，令 $\lambda_1=\lambda_2=1,\frac{1}{2},\alpha=\beta=1$ ，即得改进的欧拉公式

四阶R-K法

要进一步提高精度，必须取更多的点，如取四点构造如下形式的公式
::: align-center
$\left\{\begin{aligned} &y_{n+1}=y_n+h(\lambda_1 k_1+\lambda_2 k_2+\lambda_3 k_3+\lambda_4 k_4)\\ &k_1=f(t_n,y_n)\\ &k_2=f(t_n+\alpha_1hy_n,y_n+\beta_1 hk_1)\\ &k_3=f(t_n+\alpha_2hy_n,y_n+\beta_2 hk_1+\beta_3hk_2)\\ &k_4=f(t_n+\alpha_3hy_n,y_n+\beta_4 hk_1+\beta_5hk_2+\beta_6hk_3) \end{aligned} \right .$

:::
其中共有13个待定系数，计算更为复杂，此时局部误差小于 $O(h^5)$ ，其中满足这些方程的一组形式比较简便的系数组成的公式为
::: align-center
$\left\{\begin{aligned} &y_{n+1}=\frac{h}{6}(k_1+2k_2+2k_3+k_4)\\ &k_1=f(x_n,y_n)\\ &k_2=f(x_n+\frac{h}{2},y_n+\frac{hk_1}{2})\\ &k_3=f(x_n+\frac{h}{2},y_n+\frac{hk_2}{2})\\ &k_4=f(x_n+h,y_n+hk_3) \end{aligned} \right .$

:::
这就是常用的四阶龙格-库塔方法（简称RK方法）