龙格库塔算法原理详解

龙格库塔算法

前言: 本文从简单欧拉法讲起，详细探究了龙格库塔算法的原理。

0. 问题提出

求解下面的一阶微分方程
$$\{\begin{array}{rl}& y^{\prime }(x)=\frac{dy}{dx}=f(x,y)\\ & y(x_0)=y_0\end{array}$$
即已知上述条件，则 y(x) 如何求出？

1. 欧拉法

1.1 基本原理

上述方程组给出了某点的坐标，以及该点的斜率值。可以使用欧拉法对 y(x) 进行迭代的近似求解。求解分析如下：

已知 $x=x_0$ 时 $y=y_0$ ，且斜率为 $y^{\prime }(x_0)=f(x_0,y_0)$ ，那么就可以用过该点且斜率为$y^{\prime }(x_0)$的直线近似求解$x_0$附近一点$x_1$的函数值。

构建直线：
$$y(x)=y_0+f(x_0,y_0)(x-x_0)$$
那么，在$x=x_1$时，求得
$$y(x_1)=y_0+f(x_0,y_0)(x_1-x_0)$$
即在 $x=x_1$ 时，$y=y(x_1)$，且斜率为 $y^{\prime }(x_1)=f(x_1,y_1)$。重复以上步骤，可以求得 $y(x_2),y(x_3),\cdots ,y(x_n)$

如果等间隔取 x 值，即 $x_1-x_0=x_2-x_1=x_3-x_2=\cdots =x_n-x_{n-1}=h$，也即$x_n=x_0+n\cdot h$可以得到
$$\begin{gathered} y(x_1)=y_0+f(x_0,y_0)\cdot h \\ y(x_2)=y_1+f(x_1,y_1)\cdot h \\ ⋮ \\ y(x_n)=y_{n-1}+f(x_{n-1},y_{n-1})\cdot h \end{gathered}$$
总结为：
$$\{\begin{array}{rl}& y(x_n)=y_{n-1}+f(x_{n-1},y_{n-1})\cdot h\\ & x_n=x_0+n\cdot h\end{array}$$
该公式就是欧拉公式。

1.2 误差分析

值得注意的是，如上求解 $y(x_1)$ 是有误差的，而在求解 $y(x_2)$ 的时候，虽然斜率是准确的，但函数值是基于 $y(x_1)$ 的，即在第一次误差的基础上计算的，因此会累积误差。

可以使用泰勒展开分析每次迭代产生的误差。对 $y(x)$ 进行泰勒展开
y ( x ) = y ( x k ) + y ′ ( x k ) ( x − x k ) + y ′ ′ ( x k ) 2 ! ( x − x k ) 2 + ⋯ + y n + 1 ( ϵ ) ( n + 1 ) ! ( x − x k ) n + 1 y(x) = y(x_k) +y'(x_k)(x-x_k)+\frac{y''(x_k)}{2!}(x-x_k)^2+\cdots+\frac{y^{n+1}(\epsilon)}{(n+1)!}(x-x_k)^{n+1} \\ y(x)=y(xk​)+y′(x