边界离散

考虑如下的一维初边值问题 \[ \left\{ \begin{aligned} & u_t = u_{xx} + f(x,t), \,\, x \in (0,1),t>0\\ & u(x,0) = u_0(x),\,\, x \in (0,1) \end{aligned} \right. \] 我们取左侧为整网格，即 \(x_0 = 0\)，\(x_1 = \Delta x\)，边界条件为 \[ B_1(u) = a u + b u_x = g(t),\,(x = 0) \] 取右侧网格为半网格，即 \(x_{N-1} = 1-\frac{\Delta x}2\)，\(x_N = 1 + \frac{\Delta x}2\)，边界条件为 \[ B_2(u) = c u + d u_x = h(t),\,(x=1) \] 这里 \(a,b,c,d\) 均为常数。

对于内部直接采用二阶离散 \[ \frac{v_{j}^{n+1}-v_j^n}{\Delta t} = \frac{v_{j+1}^n-2v_j^n + v_{j-1}^n}{\Delta x^2} + f_j^n \]

对于左侧整网格的边界离散有很多方案，一阶方案例如（A1） \[ a v_0^{n+1} + b \frac{v_1^{n+1}-v_0^{n+1}}{\Delta x} = g(t^{n+1}) \] 二阶方案通常需要引入虚拟点 \(x_{-1} = -\Delta x\)，例如第一种方案为（A2） \[ a v_0^{n+1} + b\frac{v_1^{n+1}-v_{-1}^{n+1}}{2\Delta x} = g(t^{n+1}) \] 第二种方案为（A3） \[ a \frac{v_{-1}^{n+1}+v_1^{n+1}}2 + b\frac{v_1^{n+1}-v_{-1}^{n+1}}{2\Delta x} = g(t^{n+1}) \]

对于右侧半网格的边界离散例如（B） \[ c \frac{v_{N-1}^{n+1}+v_N^{n+1}}{2} + d \frac{v_N^{n+1}-v_{N-1}^{n+1}}{\Delta x} = h(t^{n+1}) \]

下面我们依次探究这几种边界处理方案的相容性。

逐点相容性

这几种离散都满足逐点相容性，内部的逐点相容性是显然的 \[ \begin{aligned} \tau_j^n &= \frac{u_{j}^{n+1}-u_j^n}{\Delta t} - \frac{u_{j+1}^n-2u_j^n + u_{j-1}^n}{\Delta x^2} - f_j^n\\ &= (u_t)_j^n + \frac{\Delta t}2 (u_{tt})_j^n + O(\Delta t^2) - (u_{xx})_j^n - \frac{\Delta x^2}{12}(u_{xxxx})_j^n + O(\Delta x^4) - f_j^n\\ &= \frac{\Delta t}2 (u_{tt})_j^n - \frac{\Delta x^2}{12}(u_{xxxx})_j^n + O(\Delta t^2+\Delta x^4) \\ \end{aligned} \]

为了便于区分记号，我们使用 \(\tau\) 代表内部离散所对应的局部截断误差，\(\sigma\) 代表边界离散所对应的局部截断误差，\(\Phi\) 代表模相容性验证中的余项。

对（A1）至少为一阶 \[ \begin{aligned} \sigma_0^{n+1} &= a u_0^{n+1} + b \frac{u_1^{n+1}-u_0^{n+1}}{\Delta x} - g(t^{n+1}) \\ &= a u_0^{n+1} + b (u_x)_0^{n+1} + \frac{b}2 \Delta x (u_{xx})_0^{n+1} + O(\Delta x^2) - g(t^{n+1}) \\ &= \frac{b}2 (u_{xx})_0^{n+1} \Delta x + O(\Delta x^2) \end{aligned} \] 对（A2）至少为二阶 \[ \begin{aligned} \sigma_0^{n+1} &= a u_0^{n+1} + b\frac{u_1^{n+1}-u_{-1}^{n+1}}{2\Delta x} - g(t^{n+1}) \\ &= \frac{b}6 (u_{xxx})_0^{n+1} \Delta x^2 + O(\Delta x^4) \end{aligned} \] 对（A3）至少为二阶 \[ \begin{aligned} \sigma_0^{n+1} &= a\frac{u_{-1}^{n+1}+u_1^{n+1}}2 + b\frac{u_1^{n+1}-u_{-1}^{n+1}}{2\Delta x} - g(t^{n+1}) \\ &= (\frac{a}2 (u_{xx})_0^{n+1} + \frac{b}6 (u_{xxx})_0^{n+1}) \Delta x^2 + O(\Delta x^4) \end{aligned} \] 对（B）至少为二阶 \[ \begin{aligned} \sigma_{N-\frac12}^{n+1} ={}& c \frac{u_{N-1}^{n+1}+u_N^{n+1}}{2} + d \frac{u_N^{n+1}-u_{N-1}^{n+1}}{\Delta x} - h(t^{n+1}) \\ ={}& c (u_{N-\frac12}^{n+1} + \frac{1}8 (u_{xx})_{N-\frac12}^{n+1} \Delta x^2+ O(\Delta x^4)) \\ &+ d((u_x)_{N-\frac12}^{n+1} + \frac{1}{24}(u_{xxx})_{N-\frac12}^{n+1}\Delta x^2 + O(\Delta x^4)) - h(t^{n+1})\\ ={}&(\frac{c}8 (u_{xx})_{N-\frac12}^{n+1} + \frac{d}{24}(u_{xxx})_{N-\frac12}^{n+1} )\Delta x^2 + O(\Delta x^4) \end{aligned} \]

这里的讨论都是基于最一般的情形，可能在 \(a,b,c,d\) 取某些特殊值时达到更高的阶数。

最大模相容性

这几种离散方式的最大模相容性则是未必的，因为边界和内部的离散模板会相互影响，影响交界处的近似效果。

(A1) 整网格一阶离散

例如在左边界处，我们关注 \(x_1\) 处的实际计算格式 \[ \begin{aligned} v_1^{n+1} = v_1^n + \frac{\Delta t}{\Delta x^2}(v_{2}^n-2v_1^n + v_{0}^n) + \Delta t f_1^n \end{aligned} \] 这里必须将边界离散耦合在一起分析，例如（A1） \[ a v_0^{n} + b \frac{v_1^{n}-v_0^{n}}{\Delta x} = g(t^{n}) \] 可以得到 \[ \begin{aligned} \Phi^n_1 ={}& \frac{u_1^{n+1} - u_1^n}{\Delta t} - \frac{u_{2}^n-2 u_1^n + \widetilde{u}_{0}^n}{\Delta x^2} + f_1^n \\ & \text{where}\,\, a \widetilde{u}_0^{n} + b \frac{u_1^{n}-\widetilde{u}_0^{n}}{\Delta x} = g(t^{n}) \end{aligned} \] 改写为上述形式可以最大程度地减少计算量，利用前面的逐点相容性结果易得 \[ \begin{aligned} \Phi^n_1 = \tau_1^n - \frac{\widetilde{u}_{0}^n-u_{0}^n}{\Delta x^2},\,\,\, \sigma_0^n + \left(a - \frac{b}{\Delta x}\right)(\widetilde{u}_{0}^n-u_{0}^n) = 0 \end{aligned} \] 代入即可得到 \[ \begin{aligned} \Phi^n_1 &= \tau_1^n + \frac{\sigma_0^n}{(a \Delta x - b)\Delta x} \end{aligned} \] 只需要注意到 \[ \frac{\sigma_0^n}{(a \Delta x - b)\Delta x} = \frac{1}{a\Delta x - b}\left(\frac{b}2 (u_{xx})_0^n + O(\Delta x)\right) = -\frac12 (u_{xx})_0^n + O(\Delta x) \] 就可以得到：使用边界处理（A1）的格式是最大模不相容的！

(A2),(A3) 整网格二阶离散

对于在左边界引入虚拟点的（A2）和（A3），我们需要关注的是 \(x_0\) 处的离散 \[ \begin{aligned} v_0^{n+1} = v_0^n + \frac{\Delta t}{\Delta x^2}(v_{1}^n-2v_0^n + v_{-1}^n) + \Delta t f_0^n \end{aligned} \] 将边界离散（A2）耦合在一起 \[ a v_0^{n} + b\frac{v_1^{n}-v_{-1}^{n}}{2\Delta x} = g(t^{n}) \] 可以得到 \[ \begin{aligned} \Phi^n_0 ={}& \frac{u_0^{n+1} - u_0^n}{\Delta t} - \frac{u_{1}^n-2 u_0^n + \widetilde{u}_{-1}^n}{\Delta x^2} + f_0^n \\ & \text{where}\,\, a u_0^{n} + b\frac{u_1^{n}-\widetilde{u}_{-1}^{n}}{2\Delta x} = g(t^{n}) \end{aligned} \] 利用前面的逐点相容性结果易得 \[ \begin{aligned} \Phi^n_0 = \tau_0^n - \frac{\widetilde{u}_{-1}^n-u_{-1}^n}{\Delta x^2}, \,\,\, \sigma_0^n - \frac{b}{2\Delta x} (\widetilde{u}_{-1}^n-u_{-1}^n) = 0 \end{aligned} \] 代入即可得到 \[ \begin{aligned} \Phi^n_1 &= \tau_1^n - \frac{2\sigma_0^n}{b \Delta x} \end{aligned} \] 注意到 \[ \begin{aligned} \sigma_0^{n} = \frac{b}6 (u_{xxx})_0^{n} \Delta x^2 + O(\Delta x^4) \end{aligned} \] 因此使用边界处理（A2）的格式是最大模相容的，但是空间阶数相比于逐点会掉一阶。（如果 \(a u + b u_x = g(t)\) 取某些特殊情况的话，可能让对应系数为零，阶数提高一点）

对于边界离散（A3）的分析也是类似的 \[ a \frac{v_{-1}^{n}+v_1^{n}}2 + b\frac{v_1^{n}-v_{-1}^{n}}{2\Delta x} = g(t^{n}) \] 耦合在一起可以得到 \[ \begin{aligned} \Phi^n_0 ={}& \frac{u_0^{n+1} - u_0^n}{\Delta t} - \frac{u_{1}^n-2 u_0^n + \widetilde{u}_{-1}^n}{\Delta x^2} + f_0^n \\ & \text{where}\,\, a \frac{\widetilde{u}_{-1}^{n}+u_1^{n}}2 + b\frac{u_1^{n}-\widetilde{u}_{-1}^{n}}{2\Delta x} = g(t^{n}) \end{aligned} \] 利用前面的逐点相容性结果易得 \[ \begin{aligned} \Phi^n_0 = \tau_0^n - \frac{\widetilde{u}_{-1}^n-u_{-1}^n}{\Delta x^2}, \,\,\, \sigma_0^n + \frac12 \left(a - \frac{b}{\Delta x}\right) (\widetilde{u}_{-1}^n-u_{-1}^n) = 0 \end{aligned} \] 代入即可得到 \[ \begin{aligned} \Phi^n_1 &= \tau_1^n + \frac{2\sigma_0^n}{(a \Delta x - b)\Delta x} \end{aligned} \] 注意到 \[ \begin{aligned} \frac{2\sigma_0^n}{(a \Delta x - b)\Delta x} ={}& \frac{2}{a\Delta x - b}\left(\left(\frac{a}2 (u_{xx})_0^{n} + \frac{b}6 (u_{xxx})_0^{n}\right) \Delta x + O(\Delta x^3)\right) \\ ={}& -\left(\frac{a}b (u_{xx})_0^{n} + \frac13 (u_{xxx})_0^{n}\right) \Delta x + O(\Delta x^2) \end{aligned} \] 因此使用边界处理（A3）的格式和（A2）一样，虽然是最大模相容的，但是空间阶数相比于逐点会掉一阶。

(B) 半网格二阶离散

对于在右边界的边界离散（B），我们需要关注的是 \(x_{N-1}\) 处的实际计算格式 \[ \begin{aligned} v_{N-1}^{n+1} = v_{N-1}^n + \frac{\Delta t}{\Delta x^2}(v_{N}^n-2v_{N-1}^n + v_{N-2}^n) + \Delta t f_{N-1}^n \end{aligned} \] 将边界离散（B）耦合在一起 \[ c \frac{v_{N-1}^{n}+v_N^{n}}{2} + d \frac{v_N^{n}-v_{N-1}^{n}}{\Delta x} = h(t^{n}) \] 可以得到 \[ \begin{aligned} \Phi^n_{N-1} ={}& \frac{u_{N-1}^{n+1} - u_{N-1}^n}{\Delta t} - \frac{\widetilde{u}_{N}^n-2 u_{N-1}^n + u_{N-2}^n}{\Delta x^2} + f_{N-1}^n \\ & \text{where}\,\, c \frac{u_{N-1}^{n}+\widetilde{u}_N^{n}}{2} + d \frac{\widetilde{u}_N^{n}-u_{N-1}^{n}}{\Delta x} = h(t^{n}) \end{aligned} \] 利用前面的逐点相容性结果易得 \[ \begin{aligned} \Phi^n_{N-1} = \tau_{N-1}^n - \frac{\widetilde{u}_{N}^n-u_{N}^n}{\Delta x^2}, \,\,\, \sigma_{N-1}^n + \left(\frac{c}2 + \frac{d}{\Delta x}\right) (\widetilde{u}_{N}^n-u_{N}^n) = 0 \end{aligned} \] 代入即可得到 \[ \begin{aligned} \Phi^n_{N-1} &= \tau_{N-1}^n + \frac{2\sigma_{N-1}^n}{(c \Delta x + 2 d)\Delta x} \end{aligned} \] 注意到 \[ \begin{aligned} \frac{2\sigma_{N-1}^n}{(c \Delta x + 2 d)\Delta x} ={}& \frac{2}{c \Delta x + 2 d}\left( \left(\frac{c}8 (u_{xx})_{N-\frac12}^{n} + \frac{d}{24}(u_{xxx})_{N-\frac12}^{n} \right)\Delta x + O(\Delta x^3) \right) \\ ={}& \left(\frac{c}{8d} (u_{xx})_{N-\frac12}^{n} + \frac{1}{24}(u_{xxx})_{N-\frac12}^{n} \right)\Delta x + O(\Delta x^2) \end{aligned} \] 因此使用边界处理（B）是最大模相容的，但是空间阶数相比于逐点会掉一阶。