偏微分方程算法之二阶双曲型方程交替方向隐格式

本文主要是介绍偏微分方程算法之二阶双曲型方程交替方向隐格式，希望对大家解决编程问题提供一定的参考价值，需要的开发者们随着小编来一起学习吧！

一、研究目标

二、理论推导

2.1 显差分格式

2.2 交替方向隐格式

三、算例实现

四、结论

一、研究目标

本节的研究对象为如下的双曲型偏微分方程初边值问题：

$\left\{\begin{matrix} \frac{\partial^{2}u(x,y,t)}{\partial t^{2}}-(\frac{\partial^{2}u(x,y,t)}{\partial x^{2}}+\frac{\partial^{2}u(x,y,t)}{\partial y^{2}})=f(x,y,t),(x,y)\in \Omega=[0,a]\times [0,b],0<t\leqslant T,\\ u(x,y,0)=\varphi(x,y),\frac{\partial u}{\partial t}(x,y,0)=\Psi(x,y),(x,y)\in \Omega \space\space\space\space(1),\\ u(0,y,t)=g_{1}(y,t),u(a,y,t)=g_{2}(y,t),0\leqslant y \leqslant b,0<t\leqslant T,\\ u(x,0,t)=g_{3}(x,t),u(x,b,t)=g_{4}(x,t),0\leqslant x \leqslant a,0<t \leqslant T \end{matrix}\right.$

为了保证连续性，公式（1）中相关函数满足：

$g_{1}(0,t)=g_{3}(0,t),g_{1}(b,t)=g_{4}(b,t)\\g_{2}(0,t)=g_{3}(a,t),g_{2}(b,t)=g_{4}(a,t)$

$\varphi(0,y)=g_{1}(y,0),\varphi(a,y)=g_{2}(y,0)\\\varphi(x,0)=g_{3}(x,0),\varphi(x,b)=g_{4}(x,0)$

二、理论推导

2.1 显差分格式

第一步：网格剖分。在三维长方体空间（二维平面加一维时间轴）进行网格剖分，将区域[0,a]等分m份，区域[0,b]等分n份，区域[0,T]等分l份，即

$x_{i}=i\cdot\Delta x=\frac{ia}{m},0\leqslant i\leqslant m,y_{j}=j\cdot\Delta y=\frac{jb}{n},0\leqslant j\leqslant n,t_{k}=k\cdot\Delta t=\frac{kT}{l},0\leqslant k\leqslant l$

得到网格节点坐标为 $(x_{i},y_{j},t_{k})$ 。我们利用数值方法获得精确解 $u(x,y,t)$ 在网格节点 $(x_{i},y_{j},t_{k})$ 处的近似值，即数值解 $u^{k}_{i,j}$ 。

第二步：弱化原方程，仅在离散点处成立。即

$\left\{\begin{matrix} \frac{\partial^{2}u}{\partial t^{2}}|_{(x_{i},y_{j},t_{k})}-(\frac{\partial^{2}u}{\partial x^{2}}+\frac{\partial^{2}}{\partial y^{2}})|_{(x_{i},y_{j},t_{k})}=f(x_{i},y_{j},t_{k}),0\leqslant i\leqslant m,0\leqslant j\leqslant n,0<k\leqslant l,\\ u(x_{i},y_{j},t_{0})=\varphi(x_{i},y_{j}),\frac{\partial u}{\partial t}(x_{i},y_{j},t_{0})=\Psi(x_{i},y_{j}),0\leqslant i\leqslant m,0\leqslant j\leqslant n,\space\space(2)\\ u(x_{0},y_{j},t_{k})=g_{1}(y_{j},t_{k}),u(x_{m},y_{j},t_{k})=g_{2}(y_{j},t_{k}),0\leqslant j\leqslant n,0<k\leqslant l,\\ u(x_{i},y_{0},t_{k})=g_{3}(x_{i},t_{k}),u(x_{i},y_{n},t_{k})=g_{4}(x_{i},t_{k}),0\leqslant i\leqslant m,0<k\leqslant l \end{matrix}\right.$

第三步：差商代替微商处理偏导数。方程组二阶偏导数使用中心差商来近似，即

$\frac{\partial^{2}u}{\partial t^{2}}|_{(x_{i},y_{j},t_{k})}\approx\frac{u(x_{i},y_{j},t_{k-1})-2u(x_{i},y_{j},t_{k})+u(x_{i},y_{j},t_{k+1})}{\Delta t^{2}}$

$\frac{\partial^{2}u}{\partial x^{2}}|_{(x_{i},y_{j},t_{k})}\approx\frac{u(x_{i-1},y_{j},t_{k})-2u(x_{i},y_{j},t_{k})+u(x_{i+1},y_{j},t_{k})}{\Delta x^{2}}$

$\frac{\partial^{2}u}{\partial y^{2}}|_{(x_{i},y_{j},t_{k})}\approx\frac{u(x_{i},y_{j-1},t_{k})-2u(x_{i},y_{j},t_{k})+u(x_{i},y_{j+1},t_{k})}{\Delta y^{2}}$

对于初始条件中的一阶偏导数，也采用二阶中心差商，即

$\frac{\partial u}{\partial t}|_{(x_{i},y_{j},t_{k})}\approx\frac{u(x_{i},y_{j},t_{1})-u(x_{i},y_{j},t_{-1})}{2\Delta t}$

将差商公式带入公式（2），用数值解代替精确解并忽略高阶项，可得对应的数值格式为

$\left\{\begin{matrix} \frac{u^{k-1}_{i,j}-2u^{k}_{i,j}+u^{k+1}_{i,j}}{\Delta t^{2}}-(\frac{u^{k}_{i-1,j}-2u^{k}_{i,j}+u^{k}_{i+1,j}}{\Delta x^{2}}+\frac{u^{k}_{i,j-1}-2u^{k}_{i,j}+u^{k}_{i,j+1}}{\Delta y^{2}})=f(x_{i},y_{j},t_{k}),1\leqslant i\leqslant m-1,1\leqslant j\leqslant n-1,1\leqslant k\leqslant l-1,\\ u^{0}_{i,j}=\varphi(x_{i},y_{j}),u^{1}_{i,j}=u^{-1}_{i,j}+2\Delta t\Psi(x_{i},y_{j}),0\leqslant i\leqslant m,0\leqslant j\leqslant n,\space\space(3)\\ u^{k}_{0,j}=g_{1}(y_{j},t_{k}),u^{k}_{m,j}=g_{2}(y_{j},t_{k}),0\leqslant j\leqslant n,0<k\leqslant l,\\ u^{k}_{i,0}=g_{3}(x_{i},t_{k}),u^{k}_{i,n}=g_{4}(x_{i},t_{k}),0\leqslant i\leqslant m,0<k\leqslant l \end{matrix}\right.$

为消去虚拟边界项 $u^{-1}_{i,j}$ ，让公式（3）中第1式对k=0成立，可得

$\frac{u^{-1}_{i,j}-2u^{0}_{i,j}+u^{1}_{i,j}}{\Delta t^{2}}-(\frac{u^{0}_{i-1,j}-2u^{0}_{i,j}+u^{0}_{i+1,j}}{\Delta x^{2}}+\frac{u^{0}_{i,j-1}-2u^{0}_{i,j}+u^{0}_{i,j+1}}{\Delta y^{2}})=f(x_{i},y_{j},t_{0})$

再结合公式（3）中第3式中越界的 $u^{-1}_{i,j}=u^{1}_{i,j}-2\Delta t\Psi(x_{i},y_{j})$ ，消去虚拟边界项 $u^{-1}_{i,j}$ ，可得

$2u^{1}_{i,j}-2\Delta t\Psi(x_{i},y_{j})-2u^{0}_{i,j}=\Delta t^{2}(\frac{u^{0}_{i-1,j}-2u^{0}_{i,j}+u^{0}_{i+1,j}}{\Delta x^{2}}+\frac{u^{0}_{i,j-1}-2u^{0}_{i,j}+u^{0}_{i,j+1}}{\Delta y^{2}}+f(x_{i},y_{j},t_{0}))$

记 $r_{1}=\Delta t^{2}/\Delta x^{2}, r_{2}=\Delta t^{2}/\Delta y^{2}$ ，整理上式可得

$u^{1}_{i,j}=\Delta t\Psi(x_{i},y_{j})+(r_{1}(u^{0}_{i-1,j}+u^{0}_{i+1,j})+r_{2}(u^{0}_{i,j-1}+u^{0}_{i,j+1})+f(x_{i},y_{j},t_{0})\Delta t^{2})/2+(1-r_{1}-r_{2})u^{0}_{i,j} \space\space\space\space(4)$

由公式（4）知，公式（3）即为三层显格式：

$\left\{\begin{matrix} u^{k+1}_{i,j}=r_{1}(u^{k}_{i-1,j}+u^{k}_{i+1,j})+r_{2}(u^{0}_{i,j-1}+u^{k}_{i,j+1})+2(1-r_{1}-r_{2})u^{k}_{i,j}-u^{k-1}_{i,j}=f(x_{i},y_{j},t_{k})\Delta t^{2},1\leqslant i\leqslant m-1,1\leqslant j\leqslant n-1,1\leqslant k\leqslant l-1,\\ u^{0}_{i,j}=\varphi(x_{i},y_{j}),0\leqslant i\leqslant m,0\leqslant j\leqslant n,\\ u^{1}_{i,j}=\Delta t\Psi(x_{i},y_{j})+(r_{1}(u^{0}_{i-1,j}+u^{0}_{i+1,j})+r_{2}(u^{0}_{i,j-1}+u^{0}_{i,j+1})+f(x_{i},y_{j},t_{0})\Delta t^{2})/2+(1-r_{1}-r_{2})u^{0}_{i,j},1\leqslant i\leqslant m-1,1\leqslant j\leqslant n-1,\\ u^{k}_{0,j}=g_{1}(y_{j},t_{k}),u^{k}_{m,j}=g_{2}(y_{j},t_{k}),0\leqslant j\leqslant n,0<k\leqslant l,\\ u^{k}_{i,0}=g_{3}(x_{i},t_{k}),u^{k}_{i,n}=g_{4}(x_{i},t_{k}),0\leqslant i\leqslant m,0<k\leqslant l \end{matrix}\right.$

该格式的局部截断误差为 $O(\Delta t^{2}+\Delta x^{2}+\Delta y^{2})$ ，考虑到其稳定性限制，将其改为隐格式。

2.2 交替方向隐格式

由于

$\frac{\partial ^{2}u}{\partial x^{2}}|_{(x_{i,}y_{j},t_{k})}\approx\frac{1}{2}(\frac{\partial ^{2}u}{\partial x^{2}}|_{(x_{i},y_{j},t_{k-1})}+\frac{\partial^{2}u}{\partial x^{2}}|_{(x_{i},y_{j},t_{k+1})})$

$\frac{\partial ^{2}u}{\partial y^{2}}|_{(x_{i,}y_{j},t_{k})}\approx\frac{1}{2}(\frac{\partial ^{2}u}{\partial y^{2}}|_{(x_{i},y_{j},t_{k-1})}+\frac{\partial^{2}u}{\partial y^{2}}|_{(x_{i},y_{j},t_{k+1})})$

由上面两式可将原来在点 $(x_{i},y_{j},t_{k})$ 处的连续方程

$\frac{\partial^{2}u}{\partial t^{2}}|_{(x_{i},y_{j},t_{k})}-(\frac{\partial^{2}u}{\partial x^{2}}+\frac{\partial^{2}u}{\partial y^{2}})|_{(x_{i},y_{j},t_{k})}=f(x_{i},y_{j},t_{k})$

改写为

$\frac{\partial^{2}u}{\partial t^{2}}|_{(x_{i},y_{j},t_{k})}-\frac{1}{2}((\frac{\partial^{2}u}{\partial x^{2}}+\frac{\partial^{2}u}{\partial y^{2}})|_{(x_{i},y_{j},t_{k-1})}+(\frac{\partial^{2}u}{\partial x^{2}}+\frac{\partial^{2}u}{\partial y^{2}})|_{(x_{i},y_{j},t_{k+1})})=f(x_{i},y_{j},t_{k})+O(\Delta t^{2})$

再用中心差商处理二阶偏导数，将数值解代替精确解并忽略高阶项，可得数值隐格式为：

$\frac{u^{k-1}_{i,j}-2u^{k}_{i,j}+u^{k+1}_{i,j}}{\Delta t^{2}}-\frac{1}{2}(\frac{u^{k-1}_{i-1,j}-2u^{k-1}_{i,j}+u^{k-1}_{i+1,j}}{\Delta x^{2}}+\frac{u^{k+1}_{i-1,j}-2u^{k+1}_{i,j}+u^{k+1}_{i+1,j}}{\Delta x^{2}})-\frac{1}{2}(\frac{u^{k-1}_{i,j-1}-2u^{k-1}_{i,j}+u^{k-1}_{i,j+1}}{\Delta y^{2}}+\frac{u^{k+1}_{i,j-1}-2u^{k+1}_{i,j}+u^{k+1}_{i,j+1}}{\Delta y^{2}})=f(x_{i},y_{j},t_{k}) \space\space(5)$

公式（5）在求解上存在困难，故参照二维抛物型偏微分方程初边值问题的交替方向隐格式设计思想，先将公式（5）写成算子的形式，可改写为：

$\frac{u^{k-1}_{i,j}-2u^{k}_{i,j}+u^{k+1}_{i,j}}{\Delta t^{2}}-\frac{1}{2}(\frac{\delta^{2}_{x}u^{k-1}_{i,j}}{\Delta x^{2}}+\frac{\delta^{2}_{x}u^{k+1}_{i,j}}{\Delta x^{2}})-\frac{1}{2}(\frac{\delta^{2}_{y}u^{k-1}_{i,j}}{\Delta y^{2}}+\frac{\delta^{2}_{y}u^{k+1}_{i,j}}{\Delta y^{2}})=f(x_{i},y_{j},t_{k})$

也就是 $u^{k+1}_{i,j}-2^{k}_{i,j}+u^{k+1}_{i,j}=\frac{r_{1}}{2}(\delta^{2}_{x}u^{k-1}_{i,j}+\delta^{2}_{x}u^{k+1}_{i,j})+\frac{r_{2}}{2}(\delta^{2}_{y}u^{k-1}_{i,j}+\delta^{2}_{y}u^{k+1}_{i,j})+f(x_{i},y_{j},t_{k})\Delta t^{2}$

即

$(1-\frac{r_{1}}{2}\delta^{2}_{x}-\frac{r_{2}}{2}\delta^{2}_{y})u^{k+1}_{i,j}=(\frac{r_{1}}{2}\delta^{2}_{x}+\frac{r_{2}}{2}\delta^{2}_{y}-1)u^{k-1}_{i,j}+2u^{k}_{i,j}+f(x_{i},y_{j},t_{k})\Delta t^{2}$

对上式左端添加一个辅助项，即

$(1-\frac{r_{1}}{2}\delta^{2}_{x}-\frac{r_{2}}{2}\delta^{2}_{y}+\frac{r_{1}r_{2}}{4}\delta^{2}_{x}\delta^{2}_{y})u^{k+1}_{i,j}=(\frac{r_{1}}{2}\delta^{2}_{x}+\frac{r_{2}}{2}\delta^{2}_{y}-1)u^{k-1}_{i,j}+2u^{k}_{i,j}+f(x_{i},y_{j},t_{k})\Delta t^{2} \space\space(6)$

对公式（6）进行分解，可以写为

$(1-\frac{r_{1}}{2}\delta^{2}_{x})(1-\frac{r_{2}}{2}\delta^{2}_{y})u^{k+1}_{i,j}=(\frac{r_{1}}{2}\delta^{2}_{x}+\frac{r_{2}}{2}\delta^{2}_{y}-1)u^{k-1}_{i,j}+2u^{k}_{i,j}+f(x_{i},y_{j},t_{k})\Delta t^{2} \space\space(6)$

通过引入中间变量 $V_{i,j}$ ，得

$\left\{\begin{matrix} (1-\frac{r_{1}}{2}\delta^{2}_{x})V_{i,j}=(\frac{r_{1}}{2}\delta^{2}_{x}+\frac{r_{2}}{2}\delta^{2}_{y}-1)u^{k-1}_{i,j}+2u^{k}_{i,j}+f(x_{i},y_{j},t_{k})\\ (1-\frac{r_{2}}{2}\delta^{2}_{y})u^{k+1}_{i,j}=V_{i,j} \space\space(7)\end{matrix}\right.$

上式是一个x方向的三对角线性方程组，计算第1式中的 $V_{i,j}$ 时需要用到 $V_{0,j}$ 和 $V_{m,j}$ ，这可以直接从第2式中解得，即

$V_{0,j}=(1-\frac{r_{2}}{2}\delta^{2}_{y})u^{k+1}_{0,j},V_{m,j}=(1-\frac{r_{2}}{2}\delta^{2}_{y})u^{k+1}_{m,j} \space\space(8)$

我们可以直接将公式（7）、（8）写成差分形式，即

$\left\{\begin{matrix} -\frac{r_{1}}{2}V_{i-1,j}+(1+r_{1})V_{i,j}-\frac{r_{1}}{2}V_{i+1,j}=2u^{k}_{i,j}+f(x_{i},y_{j},t_{k})\Delta t^{2}+\\\frac{r_{1}}{2}(u^{k-1}_{i-1,j}+u^{k-1}_{i+1,j})+\frac{r_{2}}{2}(u^{k-1}_{i,j-1}+u^{k-1}_{i,j+1})-(1+r_{1}+r_{2})u^{k-1}_{i,j}\\ -\frac{r_{2}}{2}u^{k+1}_{i,j-1}+(1+r_{2})u^{k+1}_{i,j}-\frac{r_{2}}{2}u^{k+1}_{i,j+1}=V_{i,j} \space\space(9)\end{matrix}\right.$

其中， $V_{0,j}=-\frac{r_{2}}{2}(u^{k+1}_{0,j-1}+u^{k+1}_{0,j+1})+(1+r_{2})u^{k+1}_{0,j}$

$V_{m,j}=-\frac{r_{2}}{2}(u^{k+1}_{m,j-1}+u^{k+1}_{m,j+1})+(1+r_{2})u^{k+1}_{m,j}$

假设第k和k-1个时间层的信息已知，这样在公式（9）的第1式中，先固定某个j，这时公式（9）的第1式就是一个m-1阶线性方程组，其系数矩阵是三对角矩阵，可以用追赶法求解。每解一个方程组就可以得到一个在[0,a]x[0,b]矩形区域内行网格点（即行节点 $(x_{i},y_{j})$ ）上的中间量V的信息，所以要得到矩形区域内所有的网格点信息，需要求解n-1个这样的线性系统。同样，公式（9）的第2式是一个y方向的三对角线性方程组时，求解时要先固定某个i，这时公式（9）的第2式是一个n-1阶线性方程组，系数矩阵是三对角矩阵，可以用追赶法求解，每解一个方程组得到第k+1个时间层上[0,a]x[0,b]矩形区域内列网格点（即 $(x_{i},y_{j},t_{k+1})$ ）上的数值解的信息，所以要得到第k+1时间层上所有网格点的信息，需要求解m-1个这样的线性系统。这样就实现了x方向和y方向交替求解的隐格式。

三、算例实现

用交替方向隐格式求解下列二维双曲型方程初边值问题：

$\left\{\begin{matrix} \frac{\partial^{2}u}{\partial t^{2}}-(\frac{\partial^{2}u}{\partial x^{2}}+\frac{\partial^{2}u}{\partial y^{2}})=\frac{4t}{(1+x^{2}+y^{2})^{2}}(1-\frac{2(x^{2}+y^{2})}{1+x^{2}+y^{2}}),0<x,y<1,0<t\leqslant 1,\\ u(x,y,0)=0,\frac{\partial u}{\partial t}(x,y,0)=\frac{1}{1+x^{2}+y^{2}},0\leqslant x,y\leqslant 1,\\ u(0,y,t)=\frac{t}{1+y^{2}},u(1,y,t)=\frac{t}{2+y^{2}},0\leqslant y\leqslant 1,0<t\leqslant 1,\\ u(x,0,t)=\frac{t}{1+x^{2}},u(x,1,t)=\frac{t}{2+x^{2}},0<x<1,0<t\leqslant 1. \end{matrix}\right.$

已知精确解为 $u(x,y,t)=\frac{t}{1+x^{2}+y^{2}}$ 。分别取步长为 $\Delta x=1/10,\Delta y=1/20,\Delta t=1/40$ 和 $\Delta x=1/20,\Delta y=1/40,\Delta t=1/80$ ，给出在节点 $(0.2i,0.25j,1.00),i=1,2,3,4,j=1,2,3$ 处的数值解和误差。

代码如下：

#include<cmath>
#include<stdlib.h>
#include<stdio.h>int main(int argc, char* argv[])
{int i, j, k, m, n, L, gap_i, gap_j;double a, b, T, r1, r2, dx, dy, dt, *x, *y, *t, ***u, **v;double *a1, *b1, *c1, *d1, *a2, *b2, *c2, *d2, *ans, temp;double f(double x, double y, double t);double phi(double x, double y);double psi(double x, double y);double g1(double y, double t);double g2(double y, double t);double g3(double x, double t);double g4(double x, double t);double *chase_algorithm(double *a, double *b, double *c, double *d, int n);double exact(double x, double y, double t);a = 1.0;b = 1.0;T = 1.0;m = 10;n = 20;L = 40;dx = a / m;dy = b / n;dt = T / L;temp = dt / dx;r1 = temp*temp;temp = dt / dy;r2 = temp*temp;printf("m=%d, n=%d, L=%d.\n", m, n, L);printf("r1=%.4f,   r2=%.4f.\n", r1, r2);x = (double *)malloc(sizeof(double)*(m + 1));for (i = 0; i <= m; i++)x[i] = i*dx;y = (double *)malloc(sizeof(double)*(n + 1));for (j = 0; j <= n; j++)y[j] = j*dy;t = (double *)malloc(sizeof(double)*(L + 1));for (k = 0; k <= L; k++)t[k] = k*dt;u = (double ***)malloc(sizeof(double **)*((m + 1)*(n + 1)*(L + 1)));v = (double **)malloc(sizeof(double *)*((m + 1)*(n + 1)));for (i = 0; i <= m; i++){u[i] = (double **)malloc(sizeof(double *)*((n + 1)*(L + 1)));v[i] = (double *)malloc(sizeof(double)*(n + 1));}for (i = 0; i <= m; i++){for (j = 0; j <= n; j++)u[i][j] = (double *)malloc(sizeof(double)*(L + 1));}for (i = 0; i <= m; i++){for (j = 0; j <= n; j++){u[i][j][0] = phi(x[i], y[j]);}}for (i = 1; i< m; i++){for (j = 1; j < n; j++){u[i][j][1] = u[i][j][0] + dt*psi(x[i], y[j]) + (r1*(u[i - 1][j][0] - 2 * u[i][j][0] + u[i + 1][j][0]) + r2*(u[i][j - 1][0] - 2 * u[i][j][0] + u[i][j + 1][0]) + dt*dt*f(x[i], y[j], t[0])) / 2.0;}}for (k = 1; k <= L; k++){for (j = 0; j <= n; j++){u[0][j][k] = g1(y[j], t[k]);u[m][j][k] = g2(y[j], t[k]);}for (i = 1; i <= m - 1; i++){u[i][0][k] = g3(x[i], t[k]);u[i][n][k] = g4(x[i], t[k]);}}a1 = (double *)malloc(sizeof(double)*(m - 1));b1 = (double *)malloc(sizeof(double)*(m - 1));c1 = (double *)malloc(sizeof(double)*(m - 1));d1 = (double *)malloc(sizeof(double)*(m - 1));for (i = 0; i< m - 1; i++){a1[i] = -r1 / 2.0;b1[i] = 1.0 + r1;c1[i] = a1[i];}a2 = (double *)malloc(sizeof(double)*(n - 1));b2 = (double *)malloc(sizeof(double)*(n - 1));c2 = (double *)malloc(sizeof(double)*(n - 1));d2 = (double *)malloc(sizeof(double)*(n - 1));for (j = 0; j < n - 1; j++){a2[j] = -r2 / 2.0;b2[j] = 1.0 + r2;c2[j] = a2[j];}for (k = 1; k < L; k++){for (j = 1; j <= n - 1; j++){for (i = 1; i <= m - 1; i++){d1[i-1]=2 * u[i][j][k] + dt*dt*f(x[i], y[j], t[k])+r1*(u[i-1][j][k-1]+u[i+1][j][k-1])/2.0+r2*(u[i][j-1][k-1]+u[i][j+1][k-1])/2.0-(1+r1+r2)*u[i][j][k-1];}v[0][j] = -r2*(u[0][j - 1][k + 1] + u[0][j + 1][k + 1]) / 2.0 + (1 + r2)*u[0][j][k + 1];d1[0] = d1[0] + r1*v[0][j] / 2.0;v[m][j] = -r2*(u[m][j - 1][k + 1] + u[m][j + 1][k + 1]) / 2.0 + (1 + r2)*u[m][j][k + 1];d1[m - 2] = d1[m - 2] + r1*v[m][j] / 2.0;ans = chase_algorithm(a1, b1, c1, d1, m - 1);for (i = 1; i <= m - 1; i++)v[i][j] = ans[i - 1];free(ans);}for (i = 1; i <= m - 1; i++){for (j = 1; j <= n - 1; j++)d2[j - 1] = v[i][j];d2[0] = d2[0] + r2*u[i][0][k + 1] / 2.0;d2[n - 2] = d2[n - 2] + r2*u[i][n][k + 1] / 2.0;ans = chase_algorithm(a2, b2, c2, d2, n - 1);for (j = 1; j <= n - 1; j++)u[i][j][k + 1] = ans[j - 1];free(ans);}}k = L ;gap_i = m / 5;gap_j = n / 4;for (i = gap_i; i <= m - 1; i = i + gap_i){for (j = gap_j; j <= n - 1; j = j + gap_j){temp = fabs(exact(x[i], y[j], t[k]) - u[i][j][k]);printf("(%.2f, %.2f, %.2f)   y=%f,  err=%.4e\n", x[i], y[j], t[k], u[i][j][k], temp);}}free(x); free(y); free(t);free(a1); free(b1); free(c1); free(d1);free(a2); free(b2); free(c2); free(d2);for (i = 0; i <= m; i++){for (j = 0; j <= n; j++)free(u[i][j]);}for (i = 0; i <= m; i++){free(u[i]);}free(u);return 0;
}double f(double x, double y, double t)
{double z,temp;z = x*x + y*y;temp = 1 - 2 * z / (1 + z);return 4*t*temp/((1+z)*(1+z));
}
double phi(double x, double y)
{return 0.0;
}
double psi(double x, double y)
{return 1.0/(1+x*x+y*y);
}
double g1(double y, double t)
{return t/(1+y*y);
}
double g2(double y, double t)
{return t/(2+y*y);
}
double g3(double x, double t)
{return t/(1+x*x);
}
double g4(double x, double t)
{return t/(2+x*x);
}
double *chase_algorithm(double *a, double *b, double *c, double *d, int n)
{int i;double * ans, *g, *w, p;ans = (double *)malloc(sizeof(double)*n);g = (double *)malloc(sizeof(double)*n);w = (double *)malloc(sizeof(double)*n);g[0] = d[0] / b[0];w[0] = c[0] / b[0];for (i = 1; i<n; i++){p = b[i] - a[i] * w[i - 1];g[i] = (d[i] - a[i] * g[i - 1]) / p;w[i] = c[i] / p;}ans[n - 1] = g[n - 1];i = n - 2;do{ans[i] = g[i] - w[i] * ans[i + 1];i = i - 1;} while (i >= 0);free(g);free(w);return ans;
}
double exact(double x, double y, double t)
{return t/(1+x*x+y*y);
}

当 $\Delta x=1/10,\Delta y=1/20,\Delta t=1/40$ 时，计算结果如下：

m=10, n=20, L=40.
r1=0.0625,   r2=0.2500.
(0.20, 0.25, 1.00)   y=0.907133,  err=1.0342e-04
(0.20, 0.50, 1.00)   y=0.775225,  err=3.0722e-05
(0.20, 0.75, 1.00)   y=0.624030,  err=4.5462e-06
(0.40, 0.25, 1.00)   y=0.817921,  err=7.4625e-05
(0.40, 0.50, 1.00)   y=0.709103,  err=1.1636e-04
(0.40, 0.75, 1.00)   y=0.580475,  err=7.6977e-05
(0.60, 0.25, 1.00)   y=0.702808,  err=1.7953e-04
(0.60, 0.50, 1.00)   y=0.620917,  err=2.0139e-04
(0.60, 0.75, 1.00)   y=0.520023,  err=1.3270e-04
(0.80, 0.25, 1.00)   y=0.587236,  err=1.3524e-04
(0.80, 0.50, 1.00)   y=0.528946,  err=1.5470e-04
(0.80, 0.75, 1.00)   y=0.453919,  err=1.1003e-04

当 $\Delta x=1/20,\Delta y=1/40,\Delta t=1/80$ 时，计算结果如下：

m=20, n=40, L=80.
r1=0.0625,   r2=0.2500.
(0.20, 0.25, 1.00)   y=0.907056,  err=2.6327e-05
(0.20, 0.50, 1.00)   y=0.775202,  err=7.7268e-06
(0.20, 0.75, 1.00)   y=0.624026,  err=9.2774e-07
(0.40, 0.25, 1.00)   y=0.817978,  err=1.8204e-05
(0.40, 0.50, 1.00)   y=0.709191,  err=2.8662e-05
(0.40, 0.75, 1.00)   y=0.580532,  err=1.9294e-05
(0.60, 0.25, 1.00)   y=0.702943,  err=4.4834e-05
(0.60, 0.50, 1.00)   y=0.621068,  err=5.0469e-05
(0.60, 0.75, 1.00)   y=0.520123,  err=3.2859e-05
(0.80, 0.25, 1.00)   y=0.587338,  err=3.3376e-05
(0.80, 0.50, 1.00)   y=0.529063,  err=3.7709e-05
(0.80, 0.75, 1.00)   y=0.454003,  err=2.6860e-05