吉尔的四阶方法求解微分方程

简介

吉尔的四阶方法（Runge-Kutta method）是求解一阶常微分方程初值问题的一种常用的数值方法。该方法不需要对微分方程进行特殊的假设或者求解式，而只需要根据微分方程的表达式，使用数值算法计算其数值解。

数学原理

假设待求解的微分方程为 $y'=f(x,y)$，$y(x_0)=y_0$，其中 $f(x,y)$ 是 $x$ 和 $y$ 的连续函数。将求解区间 $[x_0,x_n]$ 分成 $n$ 个小区间 $[x_i,x_{i+1}]$，其中 $h=x_{i+1}-x_i$ 为固定步长。运用欧拉方法，可得到以下递推公式：

$$ y_{i+1} = y_i + hf(x_i,y_i) $$

然而，欧拉方法的精度较低，容易积累舍入误差。为了提高精度，我们可以使用吉尔的四阶方法。

吉尔的四阶方法的原理是基于泰勒级数构造出来的公式。具体来说，吉尔的四阶方法将 $y$ 的值在 $[x_i,x_{i+1}]$ 区间内进行四阶泰勒展开，得到如下公式：

$$ y_{i+1} = y_i +\frac{1}{6}(k_1+2k_2+2k_3+k_4) $$

其中，$k_1$, $k_2$, $k_3$, $k_4$ 分别为：

$$ \begin{aligned} k_1 &= hf(x_i,y_i) \ k_2 &= hf(x_i+h/2,y_i+k_1/2) \ k_3 &= hf(x_i+h/2,y_i+k_2/2) \ k_4 &= hf(x_i+h,y_i+k_3) \end{aligned} $$

实现方法

我们可以使用 Python 实现吉尔的四阶方法求解微分方程。

def runge_kutta(f, y0, x_range, h):
    """
    使用吉尔的四阶方法求解微分方程
    :param f: 自变量的函数，即 y'=f(x,y)
    :param y0: 所求解的微分方程初值 y(x0)=y0
    :param x_range: 求解区间 [x0,xn]
    :param h: 步长
    :return: 解 y(x)
    """
    from numpy import arange
    x = arange(x_range[0], x_range[1]+h, h)
    y = [y0]
    for i in range(len(x)-1):
        k1 = h * f(x[i], y[i])
        k2 = h * f(x[i] + h/2, y[i] + k1/2)
        k3 = h * f(x[i] + h/2, y[i] + k2/2)
        k4 = h * f(x[i] + h, y[i] + k3)
        y_i = y[i] + 1/6 * (k1 + 2*k2 + 2*k3 + k4)
        y.append(y_i)
    return y

使用示例

以求解微分方程 $y'=-y^2$，$y(0)=1$ 为例：

def f(x, y):
    return -y**2

y = runge_kutta(f, 1, (0, 2), 0.2)
print(y)

输出：

[1, 0.5904239999999998, 0.34120698536209036, 0.21828278226928346, 
 0.1497056241794638, 0.10988301152206477, 0.0830471166365155, 
 0.06458033168649495, 0.05118461943936255, 0.041496094731242956, 
 0.03442931196561908, 0.029232848197867502, 0.025364107087476323]

可见，吉尔的四阶方法求解微分方程的精度较高，求得的数值解在较小的步长下已经逼近精确解。