偏微分方程算法之混合边界差分

一、研究对象

二、差分格式

2.1 向前欧拉格式

1. 中心差商

1.1.1 理论推导

1.1.2 算例实现

2. x=0处向前差商，x=1处向后差商

1.2.1 理论推导

1.2.2 算例实现

2.2 Crank-Nicolson格式

2.2.1 理论推导

2.2.2 算例实现

一、研究对象

这里我们以混合边界（导数边界）条件下的抛物型方程初边值问题：

$\left\{\begin{matrix} \frac{\partial u}{\partial t}-a\frac{\partial^{2}}{\partial x^{2}}=f(x,t), \space\space 0<x<1,\space\space 0<t \leqslant T,\\ u(x,0)=\varphi(x),\space\space\space\space 0 \leqslant x \leqslant 1, \space\space\space\space(1)\\ \frac{\partial u}{\partial x}(0,t)-\lambda u(0,t)=\alpha(t), \frac{\partial u}{\partial x}(1,t)+\mu u(1,t)=\beta(t),0<t\leqslant T \end{matrix}\right.$

其中， $\lambda \geqslant 0,\mu \geqslant 0$ 且当 $\lambda,\mu$ 同时为0时公式（1）中的边界条件是诺依曼条件。

二、差分格式

这里我们用向前欧拉法显格式和Crank-Nicolson格式进行差分格式建立。

2.1 向前欧拉格式

1. 中心差商

1.1.1 理论推导

网格剖分参照偏微分方程算法之向前欧拉法（Forward Euler）-CSDN博客。在节点 $(x_{i},t_{k})$ 处得到节点离散方程：

$\left\{\begin{matrix} \frac{\partial u}{\partial t}|_{({x_{i},t_{k}})}-a\frac{\partial^{2}u}{\partial x^{2}}|_{({x_{i},t_{k}})}=f(x_{i},t_{k}),\\ u(x_{i},t_{0})=\varphi(x_{i}),\space\space\space\space (2)\\ \frac{\partial u}{\partial x}(x_{0},t_{k})-\lambda u(x_{0},t_{k})=\alpha(t_{k}),\frac{\partial u}{\partial x}(x_{m},t_{k})+\mu u(x_{m},t_{k})=\beta(t_{k}) \end{matrix}\right.$

利用一阶向前差商代替微商，可得:

$\frac{\partial u}{\partial t}|_{(x_{i},t_{k})}\approx \frac{u(x_{i},t_{k+1})-u(x_{i},t_{k})}{\tau}$

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

边界条件采用中心差商

$\frac{\partial u}{\partial x}|_{(x_{0},t_{k})}\approx \frac{u(x_{1},t_{k})-u(x_{-1},t_{k})}{2h},\frac{\partial u}{\partial x}|_{(x_{m},t_{k})}\approx \frac{u(x_{m+1},t_{k})-u(x_{m-1},t_{k})}{2h}$

其中 $u(x_{-1},t_{k}),u(x_{m+1},t_{k})$ 中x变量都已经越界，属于虚拟数值，将在下文单独处理。将上面各式带入公式（2）中，将数值解代替精确解并忽略高阶项，可得到离散差分格式：

$\left\{\begin{matrix} \frac{u^{k+1}_{i}-u^{k}_{i}}{\tau}-a\frac{u^{k}_{i+1}-2u^{k}_{i}+u^{k}_{i-1}}{h^{2}}=f(x_{i},t_{k}),\space\space 1\leqslant i \leqslant m-1,\space\space 0 \leqslant k \leqslant n-1,\\ u^{0}_{i}=\varphi(x_{i}),\space\space\space\space(3) \\ \frac{u^{k}_{1}-u^{k}_{-1}}{2h}-\lambda u^{k}_{0}=\alpha(t_k),\frac{u^{k}_{m+1}-u^{k}_{m-1}}{2h}+\mu u^{k}_{m}=\beta(t_k), 1 \leqslant k \leqslant n \end{matrix}\right.$

公式（3）中第1式可以写成：

$u^{k+1}_{i}-u^{k}_{i}-r(u^{k}_{i+1}-2u^{k}_{i}+u^{k}_{i-1})=\tau f(x_{i},t_{k}),\space\space 1 \leqslant i \leqslant m-1,0 \leqslant k \leqslant n-1 \space\space\space\space (4)$

其中 $r=\frac{a\tau}{h^{2}}$ 。为处理越界问题，设公式（4）对i=0和i=1都成立，即：

$u^{k+1}_{0}-u^{k}_{0}-r(u^{k}_{1}-2u^{k}_{0}+u^{k}_{-1})=\tau f(x_{0},t_{k}),u^{k+1}_{m}-u^{k}_{m}-r(u^{k}_{m+1}-2u^{k}_{m}+u^{k}_{m-1})=\tau f(x_{m},t_{k}) \space\space\space\space(5)$

将上式与公式（3）中的第3式以及 $u^{k}_{-1}=u^{k}_{1}-2\lambda h u^{k}_{0}-2h\alpha(t_{k})$ 、 $u^{k}_{m+1}=u^{k}_{m-1}-2\mu h u^{k}_{m}+2h\beta(t_{k})$ 联立，可得：

$\left\{\begin{matrix} u^{k+1}_{0}=(1-2r-2r\lambda h)u^{k}_{0}+2ru^{k}_{1}-2rh\alpha(t_{k})+\tau f(x_{0},t_{k}), \space\space (6)\\ u^{k+1}_{m}=2ru^{k}_{m-1}+(1-2r-2r\mu h)u^{k}_{m}+2rh\beta(t_{k})+\tau f(x_{m},t_{k}) \end{matrix}\right.$

联合公式（5）、（6）可得：

$\left\{\begin{matrix} u^{k+1}_{0}=(1-2r-2r\lambda h)u^{k}_{0}+2ru^{k}_{1}-2rh\alpha(t_{k})+\tau f(x_{0},t_{k}) , [1]\\ u^{k+1}_{i}=ru^{k}_{i-1}+(1-2r)u^{k}_{i}+ru^{k}_{i+1}+\tau f(x_{i},t_{k}),1\leqslant i\leqslant m-1,0\leqslant k\leqslant n-1,[2]\\ u^{k+1}_{m}=2ru^{k}_{m-1}+(1-2r-2r\mu h)u^{k}_{m}+2rh\beta(t_{k})+\tau f(x_{m},t_{k}),[3]\\ u^{0}_{i}=\varphi(x_{i}),0\leqslant i\leqslant, m[4] \end{matrix}\right.$

1.1.2 算例实现

抛物型初边值问题：

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

已知精确解为 $u(x,t)=4\sum_{n=1}^{\infty}[\frac{sec\alpha_{n}}{(3+4\alpha^{2}_{n})}e^{-4\alpha^{2}_{n}t}cos2\alpha_{n}(x-\frac{1}{2})]$ ，其中 $\alpha_{n}$ 是方程 $\alpha tan\alpha=\frac{1}{2}$ 的根。取 $h=0.1,\tau=0.0025$ 。

代码如下：

#include<cmath>
#include<stdio.h>
#include<stdlib.h>


int main(int argc, char* argv[])
{
        int m, n, i, k;
        double h, tau, a,lambda,mu,r;
        double *x, *t,**u;
        double f(double x, double t);
        double phi(double x);
        double alpha(double t);
        double beta(double t);

        m=10;
        n=400;
        h=1.0/m;
        tau=1.0/n;
        a=1.0;
        lambda=1.0;
        mu=1.0;
        r=a*tau/(h*h);
        printf("r=%.4f.\n", r);

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

        t=(double *)malloc(sizeof(double)*(n+1));
        for(k=0;k<=n;k++)
                t[k]=k*tau;

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

        for(i=0;i<=m;i++)
                u[i][0]=phi(x[i]);

        for(k=0;k<n;k++)
        {
                u[0][k+1]=(1.0-2*r-2*r*lambda*h)*u[0][k]+2*r*u[1][k]-2*r*h*alpha(t[k])+tau*f(x[0], t[k]);
                for(i=1;i<m;i++)
                        u[i][k+1]=r*u[i-1][k]+(1-2*r)*u[i][k]+r*u[i+1][k]+tau*f(x[i],t[k]);
                u[m][k+1]=2*r*u[m-1][k]+(1.0-2*r-2*r*mu*h)*u[m][k]+2*r*h*beta(t[k])+tau*f(x[m],t[k]);
        }

        printf("t/x     0          0.1       0.2       0.3       0.4       0.5\n");

        for(k=1;k<=8;k++)
        {
                printf("%.4f  ", t[k]);
                for(i=0;i<=m/2;i++)
                        printf("%.4f    ", u[i][k]);
                printf("\n");
        }

        printf("\n");
        printf("……\n");
        printf("\n");
        printf("0.1000  ");

        for(i=0;i<=m/2;i++)
                printf("%.4f    ", u[i][40]);
        printf("\n");

        for(k=1; k<=4; k=2*k)
        {
                printf("%.4f  ", t[k*100]);
                for(i=0;i<=m/2;i++)
                         printf("%.4f    ", u[i][k*100]);
                printf("\n");
        }

        return 0;
}


double f(double x, double t)
{
        return 0;
}
double phi(double x)
{
        return 1.0;
}
double alpha(double t)
{
        return 0.0;
}
double beta(double t)
{
        return 0.0;
}

结果如下：

r=0.2500.
t/x     0          0.1       0.2       0.3       0.4       0.5
0.0025  0.9500    1.0000    1.0000    1.0000    1.0000    1.0000
0.0050  0.9275    0.9875    1.0000    1.0000    1.0000    1.0000
0.0075  0.9111    0.9756    0.9969    1.0000    1.0000    1.0000
0.0100  0.8978    0.9648    0.9923    0.9992    1.0000    1.0000
0.0125  0.8864    0.9549    0.9872    0.9977    0.9998    1.0000
0.0150  0.8764    0.9459    0.9818    0.9956    0.9993    0.9999
0.0175  0.8673    0.9375    0.9762    0.9931    0.9985    0.9996
0.0200  0.8590    0.9296    0.9708    0.9902    0.9974    0.9991

……

0.1000  0.7175    0.7829    0.8345    0.8718    0.8942    0.9017
0.2500  0.5541    0.6048    0.6452    0.6745    0.6923    0.6983
0.5000  0.3612    0.3942    0.4205    0.4396    0.4512    0.4551
1.0000  0.1534    0.1674    0.1786    0.1867    0.1917    0.1933

2. x=0处向前差商，x=1处向后差商

1.2.1 理论推导

利用一阶向前差商代替微商，可得:

$\frac{\partial u}{\partial t}|_{(x_{i},t_{k})}\approx \frac{u(x_{i},t_{k+1})-u(x_{i},t_{k})}{\tau}$

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

边界条件处理如下：

$\frac{\partial u}{\partial x}|_{(x_{0},t_{k})}\approx \frac{u(x_{1},t_{k})-u(x_{0},t_{k})}{h},\frac{\partial u}{\partial x}|_{(x_{m},t_{k})}\approx \frac{u(x_{m},t_{k})-u(x_{m-1},t_{k})}{h}$

将上式带入公式（2），将数值解代替精确解并忽略高阶项，可得离散格式：

$\left\{\begin{matrix} \frac{u^{k+1}_{i}-u^{k}_{i}}{\tau}-a\frac{u^{k}_{i+1}-2u^{k}_{i}+u^{k}_{i-1}}{h^{2}}=f(x_{i},t_{k}),\space\space 1\leqslant i \leqslant m-1,\space\space 0 \leqslant k \leqslant n-1,\\ u^{0}_{i}=\varphi(x_{i}),0 \leqslant i \leqslant m,\space\space\space\space(7) \\ \frac{u^{k}_{1}-u^{k}_{0}}{h}-\lambda u^{k}_{0}=\alpha(t_k),\frac{u^{k}_{m}-u^{k}_{m-1}}{h}+\mu u^{k}_{m}=\beta(t_k), 1 \leqslant k \leqslant n \end{matrix}\right.$

整理可得：

$\left\{\begin{matrix} u^{k+1}_{i}=ru^{k}_{i-1}+(1-2r)u^{k}_{i}+ru^{k}_{i+1}+\tau f(x_{i},t_{k}),1 \leqslant i\leqslant m-1,0\leqslant k\leqslant n-1,[1] \\ u^{0}_{i}=\varphi(x_{i}),0\leqslant i\leqslant m,[2] \\u^{k}_{0}=\frac{u^{k}_{1}-h\alpha(t_{k})}{1+\lambda h},[3] \\ u^{k}_{m}=\frac{u^{k}_{m-1}+h\beta(t_{k})}{1+\mu h},1\leqslant k\leqslant n,[4] \end{matrix}\right.$

1.2.2 算例实现

抛物型初边值问题：

代码如下：

#include<cmath>
#include<stdio.h>
#include<stdlib.h>


int main(int argc, char* argv[])
{
        int m, n, i, k;
        double h, tau, a,lambda,mu,r;
        double *x, *t,**u;
        double f(double x, double t);
        double phi(double x);
        double alpha(double t);
        double beta(double t);

        m=10;
        n=400;
        h=1.0/m;
        tau=1.0/n;
        a=1.0;
        lambda=1.0;
        mu=1.0;
        r=a*tau/(h*h);
        printf("r=%.4f.\n", r);

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

        t=(double *)malloc(sizeof(double)*(n+1));
        for(k=0;k<=n;k++)
                t[k]=k*tau;

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

        for(i=0;i<=m;i++)
                u[i][0]=phi(x[i]);

        for(k=0;k<n;k++)
        {
                 for(i=1;i<m;i++)
                     u[i][k+1]=r*u[i-1][k]+(1-2*r)*u[i][k]+r*u[i+1][k]+tau*f(x[i],t[k]);
                 u[0][k+1]=(u[1][k+1]-h*alpha(t[k]))/(1.0+lambda*h);
                 u[m][k+1]=(u[m-1][k+1]+h*beta(t[k]))/(1.0+mu*h);
        }

        printf("t/x     0          0.1       0.2       0.3       0.4       0.5\n");

        for(k=1;k<=8;k++)
        {
                printf("%.4f  ", t[k]);
                for(i=0;i<=m/2;i++)
                        printf("%.4f    ", u[i][k]);
                printf("\n");
        }

        printf("\n");
        printf("……\n");
        printf("\n");
        printf("0.1000  ");

        for(i=0;i<=m/2;i++)
                printf("%.4f    ", u[i][40]);
        printf("\n");

        for(k=1; k<=4; k=2*k)
        {
                printf("%.4f  ", t[k*100]);
                for(i=0;i<=m/2;i++)
                         printf("%.4f    ", u[i][k*100]);
                printf("\n");
        }

        return 0;
}


double f(double x, double t)
{
        return 0;
}
double phi(double x)
{
        return 1.0;
}
double alpha(double t)
{
        return 0.0;
}
double beta(double t)
{
        return 0.0;
}

结果如下：

r=0.2500.
t/x     0          0.1       0.2       0.3       0.4       0.5
0.0025  0.9091    1.0000    1.0000    1.0000    1.0000    1.0000
0.0050  0.8884    0.9773    1.0000    1.0000    1.0000    1.0000
0.0075  0.8734    0.9607    0.9943    1.0000    1.0000    1.0000
0.0100  0.8612    0.9473    0.9873    0.9986    1.0000    1.0000
0.0125  0.8507    0.9358    0.9801    0.9961    0.9996    1.0000
0.0150  0.8415    0.9256    0.9730    0.9930    0.9989    0.9998
0.0175  0.8331    0.9164    0.9662    0.9895    0.9976    0.9993
0.0200  0.8255    0.9080    0.9596    0.9857    0.9960    0.9985

……

0.1000  0.6901    0.7591    0.8140    0.8537    0.8778    0.8859
0.2500  0.5230    0.5753    0.6170    0.6474    0.6658    0.6720
0.5000  0.3298    0.3627    0.3890    0.4082    0.4198    0.4237
1.0000  0.1311    0.1442    0.1547    0.1623    0.1669    0.1685

2.2 Crank-Nicolson格式

边界条件采用中心差商。

2.2.1 理论推导

在虚拟节点 $(x_{i},t_{k+\frac{1}2{}})$ 处得离散方程：

$\left\{\begin{matrix} \frac{\partial u}{\partial t}|_{({x_{i},t_{k+\frac{1}{2}}})}-a\frac{\partial^{2}u}{\partial x^{2}}|_{({x_{i},t_{k+\frac{1}{2}}})}=f(x_{i},t_{k+\frac{1}{2}}),\\ u(x_{i},t_{0})=\varphi(x_{i}),\space\space\space\space (8)\\ \frac{\partial u}{\partial x}(x_{0},t_{k})-\lambda u(x_{0},t_{k})=\alpha(t_{k}),\frac{\partial u}{\partial x}(x_{m},t_{k})+\mu u(x_{m},t_{k})=\beta(t_{k}) \end{matrix}\right.$

利用差商代替微商：

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

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

$\frac{\partial u}{\partial x}|_{(x_{0},t_{k})}\approx \frac{u(x_{1},t_{k})-u(x_{-1},t_{k})}{2h},\frac{\partial u}{\partial x}|_{(x_{m},t_{k})}\approx \frac{u(x_{m+1},t_{k})-u(x_{m-1},t_{k})}{2h}$

其中 $u(x_{-1},t_{k}),u(x_{m+1},t_{k})$ 同样越界，将上式代入公式（8），用数值解代替精确解并忽略高阶项，可得离散格式：

$\left\{\begin{matrix} \frac{u^{k+1}_{i}-u^{k}_{i}}{\tau}-a\frac{u^{k}_{i+1}-2u^{k}_{i}+u^{k}_{i-1}+u^{k+1}_{i+1}-2u^{k+1}_{i}+u^{k+1}_{i-1}}{2h^{2}}=f(x_{i},t_{k+\frac{1}{2}}),\space\space 1\leqslant i \leqslant m-1,\space\space 0 \leqslant k \leqslant n-1,\\ u^{0}_{i}=\varphi(x_{i}),0 \leqslant i \leqslant m,\space\space\space\space(9) \\ \frac{u^{k}_{1}-u^{k}_{-1}}{2h}-\lambda u^{k}_{0}=\alpha(t_k),\frac{u^{k}_{m+1}-u^{k}_{m-1}}{2h}+\mu u^{k}_{m}=\beta(t_k), 1 \leqslant k \leqslant n \end{matrix}\right.$

公式（9）中第1式可写为

$-\frac{r}{2}u^{k+1}_{i-1}+(1+r)u^{k+1}_{i}-\frac{r}{2}u^{k+1}_{i+1}=\frac{r}{2}u^{k}_{i-1}+(1-r)u^{k}_{i}+\frac{r}{2}u^{k}_{i+1}+\tau f(x_{i},t_{k+\frac{1}{2}}),\space\space 1 \leqslant i \leqslant m-1,0 \leqslant k \leqslant n-1 \space\space\space\space (10)$

为处理越界问题，设公式（10）对i=0和i=m都成立，即：

$-\frac{r}{2}u^{k+1}_{-1}+(1+r)u^{k+1}_{0}-\frac{r}{2}u^{k+1}_{1}=\frac{r}{2}u^{k}_{-1}+(1-r)u^{k}_{0}+\frac{r}{2}u^{k}_{1}+\tau f(x_{0},t_{k+\frac{1}{2}})$

$-\frac{r}{2}u^{k+1}_{m-1}+(1+r)u^{k+1}_{m}-\frac{r}{2}u^{k+1}_{m+1}=\frac{r}{2}u^{k}_{m-1}+(1-r)u^{k}_{m}+\frac{r}{2}u^{k}_{m+1}+\tau f(x_{m},t_{k+\frac{1}{2}})$

将上式与公式（9）中的第3式以及 $u^{k}_{-1}=u^{k}_{1}-2\lambda h u^{k}_{0}-2h\alpha(t_{k})$ 、 $u^{k}_{m+1}=u^{k}_{m-1}-2\mu h u^{k}_{m}+2h\beta(t_{k})$ 联立，可得：

$\left\{\begin{matrix} (1+r+r\lambda h)u^{k+1}_{0}-ru^{k+1}_{1} =(1-r-r\lambda h)u^{k}_{0}+ru^{k}_{1}-rh\alpha(t_{k})-rh\alpha(t_{k+1})+\tau f(x_{0},t_{k+\frac{1}{2}})\\ (1+r+r\mu h)u^{k+1}_{m}-ru^{k+1}_{m-1} =(1-r-r\mu h)u^{k}_{m}+ru^{k}_{m-1}+rh\beta(t_{k})+rh\beta(t_{k+1})+\tau f(x_{m},t_{k+\frac{1}{2}}) \end{matrix}\right.$

联合上面两式与公式（10）可得：

$\left\{\begin{matrix} (1+r+r\lambda h)u^{k+1}_{0}-ru^{k+1}_{1}=(1-r-r\lambda h)u^{k}_{0}+ru^{k}_{1}-rh\alpha(t_{k})-rh\alpha(t_{k+1})+\tau f(x_{0},t_{k+\frac{1}{2}}),0\leqslant k\leqslant n-1\\ -\frac{r}{2}u^{k+1}_{i-1}+(1+r)u^{k+1}_{i}-\frac{r}{2}u^{k+1}_{i+1}=\frac{r}{2}u^{k}_{i-1}+(1-r)u^{k}_{i}+\frac{r}{2}u^{k}_{i+1}+\tau f(x_{i},t_{k+\frac{1}{2}}),1\leqslant i\leqslant m-1,0\leqslant k\leqslant n-1,\\ -ru^{k+1}_{m-1}+(1+r+r\mu h)u^{k+1}_{m}=ru^{k}_{m-1}+(1-r-r\mu h)u^{k}_{m}+rh\beta(t_{k})+rh\beta(t_{k+1})+\tau f(x_{m},t_{k+\frac{1}{2}}),0\leqslant k\leqslant n-1,\\ u^{0}_{i}=\varphi(x_{i}),0\leqslant i\leqslant m\end{matrix}\right.$

上式可写出矩阵形式：

$\begin{pmatrix} 1+r+r\lambda h & -r & & & & \\ -\frac{r}{2} & 1+r & -\frac{r}{2} & & 0 & \\ & & \ddots & & & \\ & & & \ddots & & \\ & 0 & & -\frac{r}{2} & 1+r &-\frac{r}{2} \\ & & & & -r & 1+r+r\mu h \end{pmatrix}$ $\begin{pmatrix} u^{k+1}_{0}\\ u^{k+1}_{1}\\ \vdots\\ \vdots\\ u^{k+1}_{m-1}\\ u^{k+1}_{m} \end{pmatrix}=$

$\begin{pmatrix} 1-r-r\lambda h & r & & & & \\ \frac{r}{2} & 1-r & \frac{r}{2} & & 0 & \\ & & \ddots & & & \\ & & & \ddots & & \\ & 0 & & \frac{r}{2} & 1-r &\frac{r}{2} \\ & & & & r & 1-r-r\mu h \end{pmatrix}$ $\begin{pmatrix} u^{k}_{0}\\ u^{k}_{1}\\ \vdots\\ \vdots\\ u^{k}_{m-1}\\ u^{k}_{m} \end{pmatrix}+$ $\begin{pmatrix} -rh\alpha(t_{k})-rh\alpha(t_{k+1})+\tau f(x_{0},t_{k+\frac{1}{2}})\\ \tau f(x_{1},t_{k+\frac{1}{2}})\\ \vdots\\ \vdots\\ \tau f(x_{m-1},t_{k+\frac{1}{2}})\\ rh\beta(t_{k})+rh\beta(t_{k+1})+\tau f(x_{m},t_{k+\frac{1}{2}}) \end{pmatrix}$

上式可用追赶法求解。

2.2.2 算例实现

抛物型初边值问题：

代码如下：

#include<cmath>
#include<stdio.h>
#include<stdlib.h>


int main(int argc, char* argv[])
{
        int m, n, i, k;
        double h, tau, a, lambda,mu,r;
        double *x, *t, *a1, *b, *c, *d, *ans, **u, tkmid;
        double f(double x, double t);
        double phi(double x);
        double alpha(double t);
        double beta(double t);
        double * chase_algorithm(double *a, double *b, double *c, double *d, int n);

        m=10;
        n=400;
        h=1.0/m;
        tau=1.0/n;
        a=1.0;
        lambda=1.0;
        mu=1.0;
        r=a*tau/(h*h);
        printf("r=%.4f\n", r);

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

        t=(double *)malloc(sizeof(double)*(n+1));
        for(k=0;k<=n;k++)
                t[k] = k*tau;

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

        for(i=0;i<=m;i++)
                u[i][0]=phi(x[i]);

        a1=(double *)malloc(sizeof(double)*(m+1));
        b=(double *)malloc(sizeof(double)*(m+1));
        c=(double *)malloc(sizeof(double)*(m+1));
        d=(double *)malloc(sizeof(double)*(m+1));
        ans=(double *)malloc(sizeof(double)*(m+1));

        for(k=0;k<n;k++)
        {
                tkmid=(t[k]+t[k+1])/2.0;
                for(i=1;i<m;i++)
                {
                         d[i]=r*u[i-1][k]/2.0+(1.0-r)*u[i][k]+r*u[i+1][k]/2.0+tau*f(x[i],tkmid);
                         a1[i]=-r/2.0;
                         b[i]=1.0+r;
                         c[i]=a1[i];
                }
                b[0]=1.0+r+r*lambda*h;
                b[m]=1.0+r+r*mu*h;
                c[0]=-r;
                a1[m]=-r;

                d[0]=(1.0-r-r*lambda*h)*u[0][k]+r*u[1][k]-r*h*alpha(t[k])-r*h*alpha(t[k+1])+tau*f(x[0],tkmid);
                d[m]=r*u[m-1][k]+(1.0-r-r*mu*h)*u[m][k]+r*h*beta(t[k])+r*h*beta(t[k+1])+tau*f(x[m],tkmid);
                ans=chase_algorithm(a1,b,c,d,m+1);
                for(i=0;i<=m;i++)
                         u[i][k+1]=ans[i];
        }

        free(a1);free(b);free(c);free(d);

        printf("t/x     0          0.1       0.2       0.3       0.4       0.5\n");
        for(k=1;k<=8;k++)
        {
                printf("%.4f  ", t[k]);
                for(i=0;i<=m/2;i++)
                         printf("%.4f    ", u[i][k]);
                printf("\n");
        }
        printf("\n");
        printf("……\n");
        printf("\n");
        printf("0.1000  ");
        for(i=0;i<=m/2;i++)
                printf("%.4f    ", u[i][40]);
        printf("\n");

        for(k=1;k<=4;k=2*k)
        {
                printf("%.4f  ", t[k*100]);
                for(i=0;i<=m/2;i++)
                         printf("%.4f    ", u[i][k*100]);
                printf("\n");
        }

        return 0;
}


double f(double x, double t)
{
        return 0;
}
double phi(double x)
{
        return 1.0;
}
double alpha(double t)
{
        return 0.0;
}
double beta(double t)
{
        return 0.0;
}
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;
}

结果如下：

r=0.2500
t/x     0          0.1       0.2       0.3       0.4       0.5
0.0025  0.9600    0.9960    0.9996    1.0000    1.0000    1.0000
0.0050  0.9347    0.9868    0.9980    0.9997    1.0000    1.0000
0.0075  0.9164    0.9765    0.9950    0.9991    0.9999    1.0000
0.0100  0.9021    0.9663    0.9910    0.9980    0.9996    0.9999
0.0125  0.8900    0.9567    0.9864    0.9964    0.9992    0.9997
0.0150  0.8795    0.9478    0.9813    0.9944    0.9985    0.9993
0.0175  0.8701    0.9394    0.9762    0.9920    0.9975    0.9988
0.0200  0.8616    0.9315    0.9709    0.9893    0.9963    0.9981

……

0.1000  0.7180    0.7834    0.8350    0.8720    0.8943    0.9017
0.2500  0.5547    0.6054    0.6458    0.6751    0.6929    0.6989
0.5000  0.3618    0.3949    0.4213    0.4404    0.4520    0.4559
1.0000  0.1540    0.1681    0.1793    0.1874    0.1924    0.1940