向前欧拉+中心差商离散边界法，向前欧拉+左边界方程用向前差分，右边界用向后差分法和Crank-Nicolson格式+边界方程用中心差分法解决抛物方程边值问题

抛物方程的一般形式为：

$$\frac{\partial u}{\partial t} = \frac{\partial^2 u}{\partial x^2} + f(x,t),\quad a \leq x \leq b,\quad t_0 \leq t \leq T$$

其中 $f(x,t)$ 是已知函数，$u(x,t)$ 是未知函数，$a$ 和 $b$ 是空间区间的两个端点，$t_0$ 和 $T$ 分别是时间区间的起始和结束时间。

我们需要解决此方程的边值问题，即在空间区间的两个端点处给出边界条件：

$$u(a,t) = g_1(t),\quad u(b,t) = g_2(t),\quad t_0 \leq t \leq T$$

现在我们来看三种求解抛物方程边值问题的方法。

1. 向前欧拉+中心差商离散边界法

在求解抛物方程时，向前欧拉+中心差商离散边界法通常被用来处理空间区间的边界条件。为了方便表示，我们将空间区间 $[a,b]$ 离散化为 $N+1$ 个节点，即

$$x_i = a + ih,\quad i=0,1,\cdots,N,$$

其中 $h=(b-a)/N$ 是节点间隔。将时间区间 $[t_0,T]$ 离散化为 $M+1$ 个时间步长，即

$$t_n = t_0 + nk,\quad n=0,1,\cdots,M,$$

其中 $k=(T-t_0)/M$ 是时间步长。

现在我们用向前欧拉差分公式对时间导数进行离散化：

$$\frac{u_{i,n+1}-u_{i,n}}{k} = \frac{u_{i+1,n}-2u_{i,n}+u_{i-1,n}}{h^2} + f_{i,n},\quad i=1,2,\cdots,N-1,\quad n=0,1,\cdots,M-1,$$

其中 $u_{i,n} \approx u(x_i,t_n)$，$f_{i,n} \approx f(x_i,t_n)$。注意到在上式中，$u_{0,n}$ 和 $u_{N,n}$ 没有直接表示出来，这意味着我们需要在边界处采用中心差商离散化方法。

对于左边界 $x=a$，我们将 $u_{0,n}$ 用中心差商离散化：

$$\frac{u_{1,n} - u_{-1,n}}{2h} = \frac{u_{0,n} - g_1(t_n)}{h}$$

将 $u_{-1,n}$ 替换为 $2u_{0,n}-u_{1,n}$，得到

$$u_{0,n} = \frac{1}{2}(u_{1,n} + g_1(t_n) - h\Delta u_{0,n}),$$

其中 $\Delta u_{0,n} = \frac{u_{1,n} - g_1(t_n)}{h}$。

对于右边界 $x=b$，我们将 $u_{N,n}$ 用中心差商离散化：

$$\frac{u_{N+1,n} - u_{N-1,n}}{2h} = \frac{g_2(t_n) - u_{N,n}}{h}$$

将 $u_{N+1,n}$ 替换为 $2u_{N,n}-u_{N-1,n}$，得到

$$u_{N,n} = \frac{1}{2}(u_{N-1,n} + g_2(t_n) + h\Delta u_{N,n}),$$

其中 $\Delta u_{N,n} = \frac{g_2(t_n) - u_{N-1,n}}{h}$。

在上述公式中，$\Delta u_{0,n}$ 和 $\Delta u_{N,n}$ 分别表示左右边界处的一阶导数。

2. 向前欧拉+左边界方程用向前差分，右边界用向后差分法

与第一种方法类似，我们将空间区间 $[a,b]$ 离散化为 $N+1$ 个节点，时间区间 $[t_0,T]$ 离散化为 $M+1$ 个时间步长。我们仍然采用向前欧拉差分公式对时间导数进行离散化：

$$\frac{u_{i,n+1}-u_{i,n}}{k} = \frac{u_{i+1,n}-2u_{i,n}+u_{i-1,n}}{h^2} + f_{i,n},\quad i=1,2,\cdots,N-1,\quad n=0,1,\cdots,M-1,$$

其中 $u_{i,n} \approx u(x_i,t_n)$，$f_{i,n} \approx f(x_i,t_n)$。

对于左边界 $x=a$，我们采用向前差分法：

$$\frac{u_{1,n}-u_{0,n}}{h} = \Delta u_{0,n},$$

其中 $\Delta u_{0,n} = \frac{u_{1,n} - g_1(t_n)}{h}$。移项得到

$$u_{0,n} = u_{1,n} - h\Delta u_{0,n}.$$

对于右边界 $x=b$，我们采用向后差分法：

$$\frac{u_{N,n}-u_{N-1,n}}{h} = \Delta u_{N,n},$$

其中 $\Delta u_{N,n} = \frac{g_2(t_n) - u_{N-1,n}}{h}$。移项得到

$$u_{N,n} = u_{N-1,n} + h\Delta u_{N,n}.$$

3. Crank-Nicolson格式+边界方程用中心差分法

Crank-Nicolson格式是一种半离散格式，即先将空间变量离散化，再将时间变量离散化。它克服了显式和隐式格式的一些缺点，既能够保证数值稳定，又能够提高计算精度。

在求解抛物方程时，我们采用以下离散化公式：

$$\frac{u_{i,n+1}-u_{i,n}}{k} = \frac{1}{2}(\frac{\partial^2 u}{\partial x^2}\big|_{i,n+1} + \frac{\partial^2 u}{\partial x^2}\big|_{i,n}) + f_{i,n+1/2},\quad i=1,2,\cdots,N-1,\quad n=0,1,\cdots,M-1,$$

其中 $\frac{\partial^2 u}{\partial x^2}\big|_{i,n+1}$ 和 $\frac{\partial^2 u}{\partial x^2}\big|_{i,n}$ 分别表示 $u(x_i,t_{n+1})$ 和 $u(x_i,t_n)$ 在 $x_i$ 处的二阶导数，$f_{i,n+1/2} \approx f(x_i,t_{n+1/2})$。

我们可以使用中心差分公式来离散化边界条件：

对于左边界 $x=a$，我们有

$$\frac{u_{1,n+1}-u_{-1,n+1}}{2h} = \frac{u_{1,n}-u_{-1,n}}{2h} + f_{0,n+1/2}$$

将 $u_{-1,n+1}$ 替换为 $2u_{0,n+1}-u_{1,n+1}$，得到

$$u_{0,n+1} = \frac{1}{2}(u_{1,n+1} + g_1(t_{n+1}) - \frac{k}{h}\Delta u_{0,n+1/2}),$$

其中 $\Delta u_{0,n+1/2} = \frac{u_{1,n}-g_1(t_{n+1/2})}{h}$。

对于右边界 $x=b$，我们有

$$\frac{u_{N+1,n+1}-u_{N-1,n+1}}{2h} = \frac{u_{N+1,n}-u_{N-1,n}}{2h} + f_{N,n+1/2}$$

将 $u_{N+1,n+1}$ 替换为 $2u_{N,n+1}-u_{N-1,n+1}$，得到

$$u_{N,n+1} = \frac{1}{2}(u_{N-1,n+1} + g_2(t_{n+1}) + \frac{k}{h}\Delta u_{N,n+1/2}),$$

其中 $\Delta u_{N,n+1/2} = \frac{g_2(t_{n+1/2})-u_{N-1,n}}{h}$。

在上述公式中，$\Delta u_{0,n+1/2}$ 和 $\Delta u_{N,n+1/2}$ 分别表示左右边界处的一阶导数。

以上三种方法都可以用来求解抛物方程边值问题，它们的实现方法略有不同，但本质上都是基于离散化方法来求解微分方程。