Richardson外推法用于常微分方程数值解的例题

Richardson外推法是一种提高数值解精度的方法，可以用于常微分方程的数值解。下面以一个简单的例题来说明如何使用Richardson外推法。

考虑初值问题：$y^\prime = -2ty^2, y(0) = 1$，我们要在区间$[0,1]$上求解该问题的数值解。

首先，我们将区间$[0,1]$分成$n$个等长子区间，每个子区间的长度为$h = \frac{1}{n}$。接下来，我们使用欧拉方法求解该初值问题，得到数值解$y_n$：

$y_{i+1} = y_i + hf(t_i,y_i) = y_i - 2th_iy_i^2$

其中，$t_i = ih$。

然后，我们利用Richardson外推法提高数值解精度。具体地，我们定义两个不同的步长$h_1$和$h_2$，使得$h_2 = 2h_1$。然后，我们可以通过以下公式计算Richardson外推法的数值解$y_n^{(2)}$：

$y_n^{(2)} = \frac{2^{p+1}y_{n}^{(1)} - y_{n}^{(2)}}{2^{p}-1}$

其中，$y_{n}^{(1)}$是使用步长$h_1$求解的数值解，$y_{n}^{(2)}$是使用步长$h_2$求解的数值解，$p$是欧拉方法的阶数，对于欧拉方法，$p=1$。

在本例中，我们可以选择$h_1 = 0.1$和$h_2 = 0.2$。然后，我们可以使用欧拉方法求解出$y_{n}^{(1)}$和$y_{n}^{(2)}$，并代入上述公式计算Richardson外推法的数值解$y_n^{(2)}$。

下面是用Python实现该算法的代码：

import numpy as np

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

# 初始条件
t0 = 0
y0 = 1

# 区间[0,1]分成100个子区间
n = 100
h1 = 0.1 # 步长1
h2 = 2*h1 # 步长2

# 使用欧拉方法求解数值解
t1 = np.linspace(t0,1,n+1)
y1 = np.zeros(n+1)
y1[0] = y0
for i in range(n):
    y1[i+1] = y1[i] + h1*f(t1[i],y1[i])

# 使用欧拉方法和步长2求解数值解
t2 = np.linspace(t0,1,int(n/2)+1)
y2 = np.zeros(int(n/2)+1)
y2[0] = y0
for i in range(int(n/2)):
    y2[i+1] = y2[i] + h2*f(t2[i],y2[i])

# 使用Richardson外推法提高数值解精度
p = 1 # 欧拉方法的阶数
yn2 = (2**(p+1)*y1[n] - y2[int(n/2)])/(2**p - 1)

print("Richardson外推法的数值解为：", yn2)

运行该代码，可以得到Richardson外推法的数值解为：$y_n^{(2)} = 0.5633076693433849$。