偏微分方程算法之混合边界条件下的差分法

一、研究目标

二、理论推导

三、算例实现

四、结论

一、研究目标

我们在前几节中介绍了Poisson方程的边值问题，接下来对椭圆型偏微分方程的混合边值问题进行探讨，研究对象为：

$\left\{\begin{matrix} -(\frac{\partial^{2}u(x,y)}{\partial x^{2}}+\frac{\partial^{2}u(x,y)}{\partial y^{2}})=f(x,y),(x,y)\in\Omega,\space\space(1)\\ \frac{\partial u(x,y)}{\partial \mathbf{n}}+\lambda(x,y)u=\mu(x,y),(x,y)\in\partial \Omega=\Gamma \end{matrix}\right.$

其中， $\Omega$ 为矩形区域 $\Omega=[(x,y)|a\leqslant x\leqslant b,c\leqslant y\leqslant d]$ ， $f(x,y)$ 为 $\Omega$ 上的连续函数， $\mathit{\mathbf{n}}$ 是 $\Gamma$ 上的单位法向量，从而 $\frac{\partial u(x,y)}{\partial\mathbf{n}}$ 表示方向导数， $\lambda(x,y),\mu(x,y)$ 为 $\Gamma$ 的连续函数且 $\lambda(x,y)$ 非负。对于矩形区域 $\Omega$ 而言，其边界上的法向量没有统一的表达式，需要对四条边界线段分别讨论。可见：

在左边界 $\mathbf{n}_{1}=(-1,0)$ ，边界条件为 $(-\frac{\partial u}{\partial x}+\lambda(x,y)u)|_{(a,y)}=\mu_{1}(a,y)$ ；

在右边界 $\mathbf{n}_{2}=(1,0)$ ，边界条件为 $(\frac{\partial u}{\partial x}+\lambda(x,y)u)|_{(b,y)}=\mu_{2}(b,y)$ ；

在下边界 $\mathbf{n}_{3}=(0,-1)$ ，边界条件为 $(-\frac{\partial u}{\partial y}+\lambda(x,y)u)|_{(x,c)}=\mu_{3}(x,c)$ ；

在上边界 $\mathbf{n}_{4}=(0,1)$ ，边界条件为 $(\frac{\partial u}{\partial y}+\lambda(x,y)u)|_{(x,d)}=\mu_{4}(x,d)$ 。

于是，公式（1）可以具体写为：

$\left\{\begin{matrix} -(\frac{\partial^{2}u(x,y)}{\partial x^{2}}+\frac{\partial^{2}u(x,y)}{\partial y^{2}})=f(x,y),(x,y)\in(a,b)\times (c,d),\\ (\frac{\partial u(x,y)}{\partial x}-\lambda(x,y)u)|_{(a,y)}=\varphi_{1}(y),c\leqslant y\leqslant d,\\ (\frac{\partial u(x,y)}{\partial x}+\lambda(x,y)u)|_{(b,y)}=\varphi_{2}(y),c\leqslant y\leqslant d,\\ (\frac{\partial u(x,y)}{\partial y}-\lambda(x,y)u)|_{(x,c)}=\Psi_{1}(x),a\leqslant x\leqslant b,\\ (\frac{\partial u(x,y)}{\partial y})+\lambda(x,y)u)|_{(x,d)},\Psi_{2}(x),a\leqslant x\leqslant b \end{matrix}\right.$

二、理论推导

首先进行矩形区域 $\Omega$ 的等距剖分，得到各网格节点 $(x_{i},y_{j})$ ，且 $x_{i}=a+i\cdot\Delta x,i=0,1,\cdot\cdot\cdot,m$ ， $\Delta x=(b-a)/m$ ； $y_{j}=c+j\cdot\Delta y,j=0,1,\cdot\cdot\cdot,n$ ， $\Delta y=(c-d)/n$ 。然后弱化原方程，使其仅在离散节点上成立，从而有

$\left\{\begin{matrix} -(\frac{\partial^{2}u(x,y)}{\partial x^{2}}+\frac{\partial^{2}u(x,y)}{\partial y^{2}})|_{(x_{i},y_{j})}=f(x_{i},y_{j}),1\leqslant i\leqslant m-1,1\leqslant j\leqslant n-1,\\ \frac{\partial u(x,y)}{\partial x}|_{(x_{0},y_{j})}-\lambda(x_{0},y_{j})u(x_{0},y_{j})=\varphi_{1}(y_{j}),0\leqslant j\leqslant n, \\ \frac{\partial u(x,y)}{\partial x}|_{(x_{m},y_{j})}+\lambda(x_{m},y_{j})u(x_{m},y_{j})=\varphi_{2}(y_{j}),0\leqslant j\leqslant n, \\ \frac{\partial u(x,y)}{\partial y}|_{(x_{i},y_{0})}-\lambda(x_{i},y_{0})u(x_{i},y_{0})=\Psi_{1}(x_{i}),0\leqslant i\leqslant m, \\ \frac{\partial u(x,y)}{\partial y}|_{(x_{i},y_{n})}+\lambda(x_{i},y_{n})u(x_{i},y_{n})=\Psi_{2}(x_{i}),0\leqslant i\leqslant m. \end{matrix}\right.$

将上式中的一阶、二阶偏导数分别用关于一阶导数的中心差商和关于二阶导数的中心差商来近似，得

$\left\{\begin{matrix} -(\frac{u(x_{i-1},y_{j})-2u(x_{i},y_{j})+u(x_{i+1},y_{j})}{\Delta x^{2}}+\frac{u(x_{i},y_{j-1})-2u(x_{i},y_{j})+u(x_{i},y_{j+1})}{\Delta y^{2}})\\ =f(x_{i},y_{j})+O(\Delta x^{2}+\Delta y^{2}),1\leqslant i\leqslant m-1,1\leqslant j\leqslant n-1,\\ \frac{u(x_{1},y_{j})-u(x_{-1},y_{j})}{2\Delta x}-\lambda(x_{0},y_{j})u(x_{0},y_{j})=\varphi_{1}(y_{j})+O(\Delta x^{2}),0\leqslant j\leqslant n,\\ \frac{u(x_{m+1},y_{j})-u(x_{m-1},y_{j})}{2\Delta x}+\lambda(x_{m},y_{j})u(x_{m},y_{j})=\varphi(y_{j})+O(\Delta x^{2}),0\leqslant j\leqslant n,\\ \frac{u(x_{i},y_{1})-u(x_{i},y_{-1})}{2\Delta y}-\lambda(x_{i},y_{0})u(x_{i},y_{0})=\Psi_{1}(x_{i})+O(\Delta y^{2}),0\leqslant i\leqslant m,\\ \frac{u(x_{i},y_{n+1})-u(x_{i},y_{n-1})}{2\Delta y}+\lambda(x_{i},y_{n})u(x_{i},y_{n})=\Psi_{2}(x_{i})+O(\Delta y^{2}),0\leqslant i\leqslant m. \end{matrix}\right.$

用数值解代替精确解并忽略高阶项，得到以下数值格式：

$\left\{\begin{matrix} -(\frac{u_{i-1,j}-2u_{i,j}+u_{i+1,j}}{\Delta x^{2}}+\frac{u_{i,j-1}-2u_{i,j}+u_{i,j+1}}{\Delta y^{2}})=f_{i,j},1\leqslant i\leqslant m-1,1\leqslant j\geqslant n-1,\\ u_{1,j}-u_{-1,j}=2\Delta x(\varphi_{1}(y_{j})+\lambda_{0,j}u_{0,j}),0\leqslant j\leqslant n,\\ u_{m+1,j}-u_{m-1,j}=2\Delta x(\varphi_{2}(y_{j})-\lambda_{m,j}u_{m,j}),0\leqslant j\leqslant n,\\ u_{i,1}-u_{i,-1}=2\Delta y(\Psi_{1}(x_{i})+\lambda_{i,0}u_{i,0}),0\leqslant i\leqslant m,\\ u_{i,n+1}-u_{i,n-1}=2\Delta y(\Psi_{2}(x_{i})-\lambda_{i,n}u_{i,n}),0\leqslant i\leqslant m. \end{matrix}\right.$ $(2)$

其中， $\lambda_{i,j}=\lambda(x_{i},y_{j}),f_{i,j}=f(x_{i},y_{j})$ 。该格式的局部截断误差为 $O(\Delta x^{2}+\Delta y^{2})$ ，同时需要处理其中的下标越界问题。由于函数的连续性，我们认为公式（2）中的第1式对i=0也成立，即有

$-(\frac{u_{-1,j}-2u_{0,j}+u_{1,j}}{\Delta x^{2}}+\frac{u_{0,j-1}-2u_{0,j}+u_{0,j+1}}{\Delta y^{2}})=f_{i,j},1\leqslant j\leqslant n-1$

再成公式（2）的第2式中解出 $u_{-1,j}=u_{1,j}-2\Delta x(\varphi_{1}(y_{j})+\lambda_{0,j}u_{0,j})$ 代入上式，得

$-\frac{2u_{1,j}-2(1+\Delta x\lambda_{0,j})u_{0,j}}{\Delta x^{2}}-\frac{u_{0,j-1}-2u_{0,j}+u_{0,j+1}}{\Delta y^{2}}=f_{i,j}-\frac{2}{\Delta x}\varphi_{1}(y_{j}),1\leqslant j\leqslant n-1$

同样，设公式（2）中的第1式分别对 $i=m,j=0,j=n$ 成立，再分别从公式（2）的第3、4、5式中解出 $u_{m+1,j},u_{i,-1},u_{i,n+1}$ 代入前面刚得到的方程，就可以处理掉越界下标，得到以下格式：

$\left\{\begin{matrix} -\frac{u_{i-1,j}-2u_{i,j}+u_{i+1,j}}{\Delta x^{2}}-\frac{u_{i,j-1}-2u_{i,j}+u_{i,j+1}}{\Delta y^{2}}=f_{i,j},1\leqslant i\leqslant m-1,1\leqslant j\leqslant n-1,\\ -\frac{2u_{1,j}-2(1+\Delta x\lambda_{0,j})u_{0,j}}{\Delta x^{2}}-\frac{u_{0,j-1}-2u_{0,j}+u_{0,j+1}}{\Delta y^{2}}=f_{0,j}-\frac{2}{\Delta x}\varphi_{1}(y_{j}),1\leqslant j\leqslant n-1,\\ -\frac{2u_{m-1,j}-2(1+\Delta x\lambda_{m,j})u_{m,j}}{\Delta x^{2}}-\frac{u_{m,j-1}-2u_{m,j}+u_{m,j+1}}{\Delta y^{2}}=f_{m,j}+\frac{2}{\Delta x}\varphi_{2}(y_{j}),1\leqslant j\leqslant n-1,\\ -\frac{u_{i-1,0}-2u_{i,0}+u_{i+1,0}}{\Delta x^{2}}-\frac{2u_{i,1}-2(1+\Delta y\lambda_{i,0})u_{i,0}}{\Delta y^{2}}=f_{i,0}-\frac{2}{\Delta y}\Psi_{1}(x_{i}),1\leqslant i\leqslant m-1,\\ -\frac{u_{i-1,n}-2u_{i,n}+u_{i+1,n}}{\Delta x^{2}}-\frac{2u_{i,n-1}-2(1+\Delta y\lambda_{i,n})u_{i,n}}{\Delta y^{2}}=f_{i,n}+\frac{2}{\Delta y}\Psi_{2}(x_{i}),1\leqslant i\leqslant m-1 \end{matrix}\right.$

至此我们共有（m+1）x（n+1）个待求量 $u_{i,j},0\leqslant i\leqslant m,0\leqslant j\leqslant n$ ,而现有（m-1）x（n-1）个关于内节点的方程，2（n-1）个关于左、右边界上的节点（不含端点）的方程及2（m-1）个关于上、下边界上的节点（也不含端点）的方程，还需要补充

$(m+1)(n+1)-(m-1)(n-1)-2(n-1)-2(m-1)=4$

个方程，也就是关于矩形区域的4个顶点的方程。为此，设公式（2）中第1式对 $i=0,j=0$ 成立，即

$-\frac{u_{-1,0}-2u_{0,0}+u_{1,0}}{\Delta x^{2}}-\frac{u_{0,-1}-2u_{0,0}+u_{0,1}}{\Delta y^{2}}=f_{0,0} \; (3)$

然后再从公式（2）的第2式和第4式中分别解出 $u_{-1,0}=u_{1,0}-2\Delta x(\varphi_{1}(y_{0})+\lambda_{0,0}u_{0,0})$ 和 $u_{0,-1}=u_{0,1}-2\Delta y(\Psi_{1}(x_{0})+\lambda_{0,0}u_{0,0})$ ，代入公式（3）中，得到 $\Omega$ 左下顶点处的方程：

$-\frac{2u_{1,0}-2(1+\Delta x\lambda_{0,0})u_{0,0}}{\Delta x^{2}}-\frac{2u_{0,1}-2(1+\Delta y\lambda_{0,0})u_{0,0}}{\Delta y^{2}}=f_{0,0}-\frac{2}{\Delta x}\varphi_{1}(y_{0})-\frac{2}{\Delta y}\Psi_{1}(x_{0})$

上式可以简化为

$-\frac{2u_{1,0}}{\Delta x^{2}}+[\frac{2(1+\Delta x\lambda_{0,0})}{\Delta x^{2}}+\frac{2(1+\Delta y\lambda_{0,0})}{\Delta y^{2}}]u_{0,0}-\frac{2u_{0,1}}{\Delta y^{2}}=f_{0,0}-\frac{2}{\Delta x}\varphi_{1}(y_{0})-\frac{2}{\Delta y}\Psi_{1}(x_{0})$

同样，设公式（2）中第1式分别对 $i=m,j=0$ 、 $i=0,j=n$ 和 $i=m,j=n$ 成立，然后再从公式（2）中第3式和第4式解出 $u_{m+1,0},u_{m,-1},u_{-1,n},u_{0,n+1},u_{m+1,n},u_{m,n+1}$ 分别代入刚才得到的3个方程，就得到 $\Omega$ 的右下顶点、左上顶点、和右上顶点处的方程，这样，我们就有了完整的处理带导数边界条件的椭圆型方程的数值格式：

$\left\{\begin{matrix} -\frac{u_{i,j-1}}{\Delta y^{2}}-\frac{u_{i-1,j}}{\Delta x^{2}}+2[\frac{1}{\Delta x^{2}}+\frac{1}{\Delta y^{2}}]u_{i,j}-\frac{u_{i+1,j}}{\Delta x^{2}}-\frac{u_{i,j+1}}{\Delta y^{2}}=f_{i,j},1\leqslant i\leqslant m-1,1\leqslant j\leqslant n-1,\\ -\frac{u_{0,j-1}}{\Delta y^{2}}+[\frac{2(1+\Delta x\lambda_{0,j})}{\Delta x^{2}}+\frac{2}{\Delta y^{2}}]u_{0,j}-\frac{2u_{1,j}}{\Delta x^{2}}-\frac{u_{0,j+1}}{\Delta y^{2}}=f_{0,j}-\frac{2}{\Delta x}\varphi_{1}(y_{j}),1\leqslant j\leqslant n-1,\\ -\frac{u_{m,j-1}}{\Delta y^{2}}-\frac{2u_{m-1,j}}{\Delta x^{2}}+[\frac{2(1+\Delta x\lambda_{m,j})}{\Delta x^{2}}+\frac{2}{\Delta y^{2}}]u_{m,j}-\frac{u_{m,j+1}}{\Delta y^{2}}=f_{m,j}+\frac{2}{\Delta x}\varphi_{2}(y_{j}),1\leqslant j\leqslant n-1,\\ -\frac{u_{i-1,0}}{\Delta x^{2}}+[\frac{2}{\Delta x^{2}}+\frac{2(1+\Delta y\lambda_{i,0})}{\Delta y^{2}}]u_{i,0}-\frac{u_{i+1,0}}{\Delta x^{2}}-\frac{2u_{i,1}}{\Delta y^{2}}=f_{i,0}-\frac{2}{\Delta y}\Psi_{1}(x_{i}),1\leqslant i\leqslant m-1,\\ -\frac{u_{i-1,n}}{\Delta x^{2}}-\frac{2u_{i,n-1}}{\Delta y^{2}}+[\frac{2}{\Delta x^{2}}+\frac{2(1+\Delta y\lambda_{i,n})}{\Delta y^{2}}]u_{i,n-\frac{u_{i+1,n}}{\Delta x^{2}}}=f_{i,n}+\frac{2}{\Delta y}\Psi_{2}(x_{i}),1\leqslant i\leqslant m-1,\\ [\frac{2(1+\Delta x\lambda_{0,0})}{\Delta x^{2}}+\frac{2(1+\Delta y \lambda_{0,0})}{\Delta y^{2}}]u_{0,0}-\frac{2u_{1,0}}{\Delta x^{2}}-\frac{2u_{0,1}}{\Delta y^{2}}=f_{0,0}-\frac{2}{\Delta x}\varphi_{1}(y_{0})-\frac{2}{\Delta y}\Psi_{1}(x_{0}),\\ -\frac{2u_{m-1,0}}{\Delta x^{2}}+[\frac{2(1+\Delta x\lambda_{m,0})}{\Delta x^{2}}+\frac{2(1+\Delta y \lambda_{m,0})}{\Delta y^{2}}]u_{m,0}-\frac{2u_{m,1}}{\Delta y^{2}}=f_{m,0}+\frac{2}{\Delta x}\varphi_{2}(y_{0})-\frac{2}{\Delta y}\Psi_{1}(x_{m}),\\ -\frac{2u_{0,n-1}}{\Delta y^{2}}+[\frac{2(1+\Delta x\lambda_{0,n})}{\Delta x^{2}}+\frac{2(1+\Delta y\lambda_{0,n})}{\Delta y^{2}}]u_{0,n}-\frac{2u_{1,n}}{\Delta x^{2}}=f_{0,n}-\frac{2}{\Delta x}\varphi_{1}(y_{n})+\frac{2}{\Delta y}\Psi_{2}(x_{0}),\\ -\frac{2u_{m,n-1}}{\Delta y^{2}}-\frac{2u_{m-1,n}}{\Delta x^{2}}+[\frac{2(1+\Delta x\lambda_{m,n})}{\Delta x^{2}}+\frac{2(1+\Delta y\lambda_{m,n})}{\Delta y^{2}}]u_{m,n}=f_{m,n}+\frac{2}{\Delta x}\varphi_{2}(y_{n})+\frac{2}{\Delta y}\Psi_{2}(x_{m}). \end{matrix}\right.$

记 $\beta =\frac{1}{\Delta x^{2}},\gamma=\frac{1}{\Delta y^{2}},\alpha=2(\frac{1}{\Delta x^{2}}+\frac{1}{\Delta y^{2}}),\xi=\frac{2}{\Delta x},\eta=\frac{2}{\Delta y}$ ，有

$\begin{Bmatrix} -\gamma u_{i,j-1}-\beta u_{i-1,j}+\alpha u_{i,j}-\beta u_{i+1,j}-\gamma u_{i,j+1}=f_{i,j},1\leqslant i\leqslant m-1,1\leqslant j\leqslant n-1,\\ -\gamma u_{0,j-1}+[\alpha+\xi \lambda_{0,j}]u_{0,j}-2\beta u_{1,j}-\gamma u_{0,j+1}=f_{0,j}-\xi\varphi_{1}(y_{j}),1\leqslant j\leqslant n-1,\\ -\gamma u_{m,j-1}-2\beta u_{m-1,j}+[\alpha+\xi\lambda_{m,j}]u_{m,j}-\gamma u_{m,j+1}=f_{m,j}+\xi \varphi_{2}(y_{j}),1\leqslant j\leqslant n-1,\\ -\beta u_{i-1,0}+[\alpha+\eta\lambda_{i,0}]u_{i,0}-\beta u_{i+1,0}-2\gamma u_{i,1}=f_{i,0}-\eta\Psi_{1}(x_{i}),1\leqslant i\leqslant m-1,\\ -2\gamma u_{i,n-1}-\beta u_{i-1,n}+[\alpha+\eta\lambda_{i,n}]u_{i,n}-\beta u_{i+1,n}=f_{i,n}+\eta\Psi_{2}(x_{i}),1\leqslant i\leqslant m-1,\\ [\alpha+(\xi+\eta)\lambda_{0,0}]u_{0,0}-2\beta u_{1,0}-2\gamma u_{0,1}=f_{0,0}-\xi\varphi_{1}(y_{0})-\eta\Psi_{1}(x_{0}),\\ -2\beta u_{m-1,0}+[\alpha+(\xi+\eta)\lambda_{m,0}]u_{m,0}-2\gamma u_{m,1}=f_{m,0}+\xi\varphi_{2}(y_{0})-\eta\Psi_{1}(x_{m}),\\ -2\gamma u_{0,n-1}+[\alpha+(\xi+\eta)\lambda_{0,n}]u_{0,n}-2\beta u_{1,n}=f_{0,n}-\xi \varphi_{1}(y_{n})+\eta\Psi_{2}(x_{0}),\\ -2\gamma u_{m,n-1}-2\beta u_{m-1,n}+[\alpha+(\xi+\eta)\lambda_{m,n}]u_{m,n}=f_{m,n}+\xi\varphi_{2}(y_{n})+\eta\Psi_{2}(x_{m}). \end{Bmatrix}\;(4)$

为计算分别，可将上述方程组写成矩阵形式。

首先，将公式（4）中第4、6、7式写成：

$\begin{pmatrix} \alpha+(\xi+\eta)\lambda_{0,0} & -2\beta & & & & & \\ -\beta & \alpha+\eta\lambda_{1,0} & -\beta & & & & \\ & -\beta & \alpha+\eta\lambda_{2,0} & -\beta & & & \\ & & \ddots & \ddots & \ddots & & \\ & & & -\beta & \alpha+\eta\lambda_{m-2,0} & -\beta & \\ & & & & -\beta & \alpha+\eta\lambda_{m-1,0} & -\beta\\ & & & & & -2\beta & \alpha+(\xi+\eta)\lambda_{m,0} \end{pmatrix}$ $+\begin{pmatrix} u_{0,0}\\ u_{1,0}\\ u_{2,0}\\ \vdots\\ u_{m-2,0}\\ u_{m-1,0}\\ u_{m,0} \end{pmatrix}$ $+\begin{pmatrix} -2\gamma & & & & & & \\ & -2\gamma & & & & & \\ & & -2\gamma & & & & \\ & & & -2\gamma & & & \\ & & & & -2\gamma & & \\ & & & & & -2\gamma & \\ & & & & & & -2\gamma \end{pmatrix}$ $\begin{pmatrix} u_{0,1}\\ u_{1,1}\\ u_{2,1}\\ \vdots\\ u_{m-2,1}\\ u_{m-1,1}\\ u_{m,1} \end{pmatrix}$ $=\begin{pmatrix} f_{0,0}-\eta\Psi_{1}(x_{0})-\xi\varphi_{1}(y_{0})\\ f_{1,0}-\eta\Psi_{1}(x_{1})\\ f_{2,0}-\eta\Psi_{1}(x_{2})\\ \vdots\\ f_{m-2,0}-\eta\Psi_{1}(x_{m-2})\\ f_{m-1,0}-\eta\Psi_{1}(x_{m-1})\\ f_{m,0}-\eta\Psi_{1}(x_{m})-\xi\varphi_{2}(y_{0}) \end{pmatrix} \; (5)$

上面的式子可以简记为

$C\mathbf{u_{0}}+2A\mathbf{u}_{1}=\mathbf{f}_{0}$

其中， $A=-\gamma I$ ， $I$ 为m+1阶单位矩阵，且C为公式（5）最左端的三对角矩阵， $\mathbf{f}_{0}$ 为公式（5）右端的向量， $\mathbf{u}_{j}=(u_{0,j},u_{1,j},\cdot\cdot\cdot ,u_{m-1,j},u_{m,j})^{T},0\leqslant j\leqslant n$ 。接着，公式（4）中的第1,2,3式可以写成

$\begin{pmatrix} -\gamma & & & & & & \\ & -\gamma & & & & & \\ & & -\gamma & & & & \\ & &\ddots &\ddots &\ddots & & \\ & & & & -\gamma & & \\ & & & & & -\gamma & \\ & & & & & & -\gamma \end{pmatrix}$ $\begin{pmatrix} u_{0,j-1}\\ u_{1,j-1}\\ u_{2,j-1}\\ \vdots\\ u_{m-2,j-1}\\ u_{m-1,j-1}\\ u_{m,j-1} \end{pmatrix}+$ $\begin{pmatrix} \alpha+\xi\lambda_{0,j} & -2\beta & & & & & \\ -\beta & \alpha & -\beta & & & & \\ & -\beta & \alpha & -\beta & & & \\ & & \ddots & \ddots & \ddots & & \\ & & & -\beta & \alpha & -\beta & \\ & & & & -\beta & \alpha & -\beta\\ & & & & & -2\beta & \alpha+\xi\lambda_{m,j} \end{pmatrix}$ $\begin{pmatrix} u_{0,j}\\ u_{1,j}\\ u_{2,j}\\ \vdots\\ u_{m-2,j}\\ u_{m-1,j}\\ u_{m,j} \end{pmatrix}+$ $\begin{pmatrix} -\gamma & & & & & & \\ & -\gamma & & & & & \\ & & -\gamma & & & & \\ & &\ddots &\ddots &\ddots & & \\ & & & & -\gamma & & \\ & & & & & -\gamma & \\ & & & & & & -\gamma \end{pmatrix}$ $\begin{pmatrix} u_{0,j+1}\\ u_{1,j+1}\\ u_{2,j+1}\\ \vdots\\ u_{m-2,j+1}\\ u_{m-1,j+1}\\ u_{m,j+1} \end{pmatrix}=$ $\begin{pmatrix} f_{0,j}-\xi \varphi_{1}(y_{j})\\ f_{1,j}\\ f_{2,j}\\ \vdots\\ f_{m-2,j}\\ f_{m-1,j}\\ f_{m,j}+\xi \varphi_{2}(y_{j}) \end{pmatrix}\;(6)$ , $1\leqslant j\leqslant n-1$

该式可以简记为

$A\mathbf{u}_{j-1}+B_{j}\mathbf{u}_{j}+A\mathbf{u}_{j+1}=\mathbf{f}_{j},1\leqslant j\leqslant n-1$

其中， $B_{j}$ 为公式（6）中的三对角矩阵， $\mathbf{f}_{j}$ 为公式（6）右端的向量。最后，公式（4）的第5、8、9式可以写成

$\begin{pmatrix} -2\gamma & & & & & & \\ & -2\gamma & & & & & \\ & & -2\gamma & & & & \\ & & \ddots& \ddots& \ddots& & \\ & & & & -2\gamma & & \\ & & & & & -2\gamma & \\ & & & & & & -2\gamma \end{pmatrix}$ $\begin{pmatrix} u_{0,n-1}\\ u_{1,n-1}\\ u_{2,n-1}\\ \vdots\\ u_{m-2,n-1}\\ u_{m-1,n-1}\\ u_{m,n-1} \end{pmatrix}+$ $\begin{pmatrix} \alpha+(\xi+\eta)\lambda_{0,n} & -2\beta & & & & & \\ -\beta & \alpha+\eta\lambda_{1,n} & -\beta & & & & \\ & -\beta & \alpha+\eta\lambda_{2,n} & -\beta & & & \\ & & \ddots & \ddots & \ddots & & \\ & & & -\beta & \alpha+\eta\lambda_{m-2,n} & -\beta & \\ & & & & -\beta & \alpha+\eta\lambda_{m-1,n} & -\beta\\ & & & & & -2\beta & \alpha+(\xi+\eta)\lambda_{m,n} \end{pmatrix}$ $\begin{pmatrix} u_{0,n}\\ u_{1,n}\\ u_{2,n}\\ \vdots\\ u_{m-2,n}\\ u_{m-1,n}\\ u_{m,n} \end{pmatrix}=$ $=\begin{pmatrix} f_{0,n}+\eta\Psi_{1}(x_{0})-\xi\varphi_{1}(y_{n})\\ f_{1,n}+\eta\Psi_{2}(x_{1})\\ f_{2,n}-\eta\Psi_{2}(x_{2})\\ \vdots\\ f_{m-2,n}+\eta\Psi_{2}(x_{m-2})\\ f_{m-1,n}+\eta\Psi_{2}(x_{m-1})\\ f_{m,n}+\eta\Psi_{2}(x_{m})+\xi\varphi_{2}(y_{n}) \end{pmatrix} \; (7)$

该式可以简记为

$2A\mathbf{u}_{n-1}+D\mathbf{u}_{n}=\mathbf{f}_{n}$

其中，D为公式（7）中的三对角矩阵， $\mathbf{f}_{n}$ 为公式（7）右端的向量。于是，由公式（5）、（6）、（7）可得到公式（4）写成块三对角矩阵为

$\begin{pmatrix} C & 2A & & & & & \\ A & B_{1} & A & & & & \\ & A & B_{2} & A & & & \\ & & \ddots & \ddots & \ddots & & \\ & & & A & B_{m-2} & A & \\ & & & & A & B_{m-1} & A\\ & & & & & 2A & D \end{pmatrix}$ $\begin{pmatrix} \mathbf{u}_{0}\\ \mathbf{u}_{1}\\ \mathbf{u}_{2}\\ \vdots\\ \mathbf{u}_{m-2}\\ \mathbf{u}_{m-1}\\ \mathbf{u}_{m} \end{pmatrix}$ $=\begin{pmatrix} \mathbf{f}_{0}\\ \mathbf{f}_{1}\\ \mathbf{f}_{2}\\ \vdots\\ \mathbf{f}_{m-2}\\ \mathbf{f}_{m-1}\\ \mathbf{f}_{m} \end{pmatrix}\;(8)$

采用Gauss-Seidel迭代法求解公式（8）。

三、算例实现

用差分格式-公式（8）求解椭圆型方程混合边值问题：

$\left\{\begin{matrix} -(\frac{\partial^{2}u}{\partial x^{2}}+\frac{\partial^{2}u}{\partial y^{2}})=(\pi^{2}-1)e^{x}sin(\pi y),0<x<2,0<y<1,\\ (\frac{\partial u(x,y)}{\partial x}-(x^{2}+y^{2})u)|_{(0,y)}=(1-y^{2})sin(\pi y),0\leqslant y\leqslant 1,\\ (\frac{\partial u(x,y)}{\partial x}+(x^{2}+y^{2})u)|_{(2,y)}=(5+y^{2})e^{2}sin(\pi y),0\leqslant y\leqslant 1,\\ (\frac{\partial u(x,y)}{\partial y}-(x^{2}+y^{2})u)|_{(x,c)}=\pi e^{x},0\leqslant x\leqslant 1,\\ (\frac{\partial u(x,y)}{\partial y}+(x^{2}+y^{2})u)|_{(x,d)}=-\pi e^{x},0\leqslant x\leqslant 1 \end{matrix}\right.$

已知精确解为 $u(x,y)=e^{x}sin(\pi y)$ 。分别取步长为 $\Delta x=\Delta y=\frac{1}{16}$ 和 $\Delta x=\Delta y=\frac{1}{32}$ ，输出6个节点 $(0.5i,0.25)$ 和 $(0.5i,0.50),i=1,2,3$ 处的数值解，并给出误差，要求在各节点处最大误差的迭代误差限为 $0.5\times10^{-10}$ 。

代码如下：

#include <cmath>
#include <stdlib.h>
#include <stdio.h>
#define pi 3.14159265359


int main(int argc, char*argv[])
{
        int m,n,i,j,k;
        double xa,xb,ya,yb,dx,dy,alpha,beta,gamma,maxerr;
        double *x,*y,**u,**v,**lambda,*d,kexi,eta,temp;
        double f(double x, double y);
        double lambda_function(double x, double y);
        double phi1(double y);
        double phi2(double y);
        double psi1(double x);
        double psi2(double x);
        double exact(double x, double y);

        xa=0.0;
        xb=2.0;
        ya=0.0;
        yb=1.0;
        m=64;
        n=32;
        printf("m=%d,n=%d.\n",m,n);

        dx=(xb-xa)/m;
        dy=(yb-ya)/n;
        beta=1.0/(dx*dx);
        gamma=1.0/(dy*dy);
        alpha=2*(beta+gamma);
        kexi=2.0/dx;
        eta=2.0/dy;

        x=(double*)malloc(sizeof(double)*(m+1));
        for(i=0;i<=m;i++)
                x[i]=xa+i*dx;

        y=(double*)malloc(sizeof(double)*(n+1));
        for(j=0;j<=n;j++)
                y[j]=ya+j*dy;

        u=(double**)malloc(sizeof(double*)*(m+1));
        v=(double**)malloc(sizeof(double*)*(m+1));
        lambda=(double**)malloc(sizeof(double*)*(m+1));
        for(i=0;i<=m;i++)
        {
                u[i]=(double*)malloc(sizeof(double)*(n+1));
                v[i]=(double*)malloc(sizeof(double)*(n+1));
                lambda[i]=(double*)malloc(sizeof(double)*(n+1));
        }

        for(i=0;i<=m;i++)
        {
                for(j=0;j<=n;j++)
                {
                        u[i][j]=0.0;
                        v[i][j]=0.0;
                        lambda[i][j]=lambda_function(x[i],y[j]);
                }
        }

        d=(double*)malloc(sizeof(double)*(m+1));
        k=0;
        do
        {
                maxerr=0.0;
                for(i=0;i<=m;i++)
                        d[i]=f(x[i],y[0])-eta*psi1(x[i]);
                d[0]=d[0]-kexi*phi1(y[0]);
                d[m]=d[m]+kexi*phi2(y[0]);
                v[0][0]=(d[0]+2*gamma*u[0][1]+2*beta*u[1][0])/(alpha+(kexi+eta)*lambda[0][0]);
                for(i=1;i<m;i++)
                        v[i][0]=(d[i]+2*gamma*u[i][1]+beta*(v[i-1][0]+u[i+1][0]))/(alpha+eta*lambda[i][0]);
                v[m][0]=(d[m]+2*gamma*u[m][1]+2*beta*v[m-1][0])/(alpha+(kexi+eta)*lambda[m][0]);

                for(j=1;j<n;j++)
                {
                        for(i=0;i<=m;i++)
                        {
                                d[i]=f(x[i],y[j]);
                        }
                        d[0]=d[0]-kexi*phi1(y[j]);
                        d[m]=d[m]+kexi*phi2(y[j]);
                        v[0][j]=(d[0]+gamma*(u[0][j+1]+v[0][j-1])+2*beta*u[1][j])/(alpha+kexi*lambda[0][j]);
                        for(i=1;i<m;i++)
                        {
                                v[i][j]=(d[i]+gamma*(v[i][j-1]+u[i][j+1])+beta*(v[i-1][j]+u[i+1][j]))/alpha;
                        }
                        v[m][j]=(d[m]+gamma*(v[m][j-1]+u[m][j+1])+2*beta*v[m-1][j])/(alpha+kexi*lambda[m][j]);
                }

                for(i=0;i<=m;i++)
                        d[i]=f(x[i],y[n])+eta*psi2(x[i]);
                d[0]=d[0]-kexi*phi1(y[n]);
                d[m]=d[m]+kexi*phi2(y[n]);
                v[0][n]=(d[0]+2*gamma*v[0][n-1]+2*beta*u[1][n])/(alpha+(kexi+eta)*lambda[0][n]);
                for(i=1;i<m;i++)
                        v[i][n]=(d[i]+2*gamma*v[i][n-1]+beta*(v[i-1][n]+u[i+1][n]))/(alpha+eta*lambda[i][n]);
                v[m][n]=(d[m]+2*gamma*v[m][n-1]+2*beta*v[m-1][n])/(alpha+(kexi+eta)*lambda[m][n]);

                for(i=0;i<=m;i++)
                {
                        for(j=0;j<=n;j++)
                        {
                                temp=fabs(u[i][j]-v[i][j]);
                                if(temp>maxerr)
                                        maxerr=temp;
                                u[i][j]=v[i][j];
                        }
                }
                k=k+1;
        }while((maxerr>0.5*1e-10)&&(k<=1e+8));
        printf("k=%d\n",k);

        k=m/4;
        for(i=k;i<m;i=i+k)
        {
                printf("(%.2f,0.25), y=%f, err=%.4e.\n",x[i],u[i][n/4],fabs(exact(x[i],y[n/4])-u[i][n/4]));
        }

        k=m/4;
        for(i=k;i<m;i=i+k)
        {
                printf("(%.2f,0.50), y=%f, err=%.4e.\n",x[i],u[i][n/2],fabs(exact(x[i],y[n/2])-u[i][n/2]));
        }

        for(i=0;i<=m;i++)
        {
                free(u[i]);
                free(v[i]);
                free(lambda[i]);
        }
        free(u);free(v);free(lambda);
        free(x);free(y);free(d);

        return 0;
}



double f(double x, double y)
{
        return (pi*pi-1)*exp(x)*sin(pi*y);
}
double lambda_function(double x, double y)
{
        return x*x+y*y;
}
double phi1(double y)
{
        return (1-y*y)*sin(pi*y);
}
double phi2(double y)
{
        return (5.0+y*y)*exp(2.0)*sin(pi*y);
}
double psi1(double x)
{
        return pi*exp(x);
}
double psi2(double x)
{
        return -pi*exp(x);
}
double exact(double x, double y)
{
        return exp(x)*sin(pi*y);
}

当 $\Delta x=\Delta y=\frac{1}{16}$ 时，计算结果如下：

m=32,n=16.
k=4323
(0.50,0.25), y=1.152179, err=1.3643e-02.
(1.00,0.25), y=1.911016, err=1.1100e-02.
(1.50,0.25), y=3.162159, err=6.8738e-03.
(0.50,0.50), y=1.638607, err=1.0115e-02.
(1.00,0.50), y=2.711255, err=7.0265e-03.
(1.50,0.50), y=4.479936, err=1.7526e-03.

当 $\Delta x=\Delta y=\frac{1}{32}$ 时，计算结果如下：

m=64,n=32.
k=16007
(0.50,0.25), y=1.162412, err=3.4097e-03.
(1.00,0.25), y=1.919341, err=2.7745e-03.
(1.50,0.25), y=3.167313, err=1.7193e-03.
(0.50,0.50), y=1.646193, err=2.5286e-03.
(1.00,0.50), y=2.716524, err=1.7575e-03.
(1.50,0.50), y=4.481249, err=4.3972e-04.