使用外推法提高一阶双曲方程的迎风格式的精度，并与原格式的误差进行对比

一阶双曲方程的迎风格式如下：

$\frac{U_i^{n+1}-U_i^n}{\Delta t}+a\frac{U_i^n-U_{i-1}^n}{\Delta x}=0$

其中，$U_i^n$ 表示在时刻 $n$ 时位置 $x_i$ 上的解，$a$ 表示波速，$\Delta t$ 和 $\Delta x$ 分别为时间和空间的离散步长。

使用外推法可以提高该格式的精度。具体来说，我们可以通过引入一个虚拟网格点 $U_{i+1/2}^n$，并在真实网格点 $U_i^n$ 和 $U_{i+1}^n$ 上分别使用迎风格式，得到：

$\frac{U_i^{n+1}-U_i^n}{\Delta t}+a\frac{U_i^n-U_{i-1}^n}{\Delta x}=0$

$\frac{U_{i+1/2}^{n+1}-U_{i+1/2}^n}{\Delta t}+a\frac{U_{i+1/2}^n-U_i^n}{\Delta x}=0$

然后，我们可以将两个方程相加，并对 $U_{i+1/2}^n$ 进行外推：

$U_{i+1/2}^n=\frac{1}{2}(U_i^n+U_{i+1}^n)-\frac{\Delta t}{2\Delta x}a(U_i^n-U_{i+1}^n)$

将 $U_{i+1/2}^n$ 代入原式，可以得到：

$\frac{U_i^{n+1}-U_i^n}{\Delta t}+a\frac{3}{2} \frac{U_i^n-U_{i-1}^n}{\Delta x}-a\frac{1}{2}\frac{U_{i-2}^n-2U_{i-1}^n+U_i^n}{\Delta x}=0$

可以看到，使用外推法后，迎风格式的精度提高了一阶。为了与原格式的误差进行对比，我们可以考虑一个简单的例子：$a=1$，$\Delta x=\Delta t=0.1$，初始条件为 $U_i^0=\sin(\pi x_i)$。在 $t=1$ 处，真实解为 $U_i^1=\sin(\pi (x_i-t))$。我们可以比较外推后的格式与原格式的数值解和误差：

import numpy as np
import matplotlib.pyplot as plt

# 迎风格式
def upwind(U, a, dt, dx):
    U_n = np.zeros_like(U)
    U_n[1:-1] = U[1:-1] - a*dt/dx*(U[1:-1]-U[:-2])
    return U_n

# 外推迎风格式
def upwind_extrap(U, a, dt, dx):
    U_n = np.zeros_like(U)
    U_half = 0.5*(U[1:]+U[:-1]) - 0.5*a*dt/dx*(U[1:]-U[:-1])
    U_n[1:-1] = U[1:-1] - a*dt/dx*(3/2*U[1:-1]-1/2*U[:-2]) - \
                a*dt/dx/2*(U_half[1:]-U_half[:-1])
    return U_n

# 初始条件
def initial(x):
    return np.sin(np.pi*x)

a = 1
dx = 0.1
dt = 0.1
x = np.arange(0, 1+dx, dx)
U = initial(x)
U_true = initial(x-dt)

# 迭代
for i in range(10):
    U = upwind(U, a, dt, dx)
    U_true = initial(x-(i+2)*dt)

U_extra = upwind_extrap(initial(x), a, dt, dx)

# 绘图
plt.plot(x, U, label='Upwind')
plt.plot(x, U_extra, label='Upwind with Extrapolation')
plt.plot(x, U_true, label='True')
plt.legend()
plt.show()

# 计算误差
error = np.abs(U-U_true)
error_extra = np.abs(U_extra-U_true)
print('Error of Upwind:', np.max(error))
print('Error of Upwind with Extrapolation:', np.max(error_extra))

从结果可以看出，外推后的迎风格式的误差比原格式更小，精度提高了一阶。