有限元之抛物型方程初边值问题解法

一、原方程的变分形式

二、有限元法进行空间半离散

三、差分法进行时间全离散

四、相关量的数值计算

五、编程时的说明

六、算例实现

6.1 C++代码

6.2 计算结果

本节我们将采用有限元法联合差分法来数值求解抛物型方程的初边值问题：

$\left\{\begin{matrix} \frac{\partial u}{\partial t}-a\frac{\partial^{2}u}{\partial x^{2}}=f(x,t),\;0<x<1,0<t\leqslant T,\\ u(x,0)=\Psi(x),\;0\leqslant x \leqslant1,\\ u(0,t)=u(1,t)=0,\;0<t\leqslant T. \end{matrix}\right. \;\;\;(1)$

其中常数 $a>0$ 。

一、原方程的变分形式

将公式（1）中的第1个方程两边同时乘以函数 $v(x)$ 并在区间[0,1]上积分，可得

$\int_{0}^{1}\frac{\partial u(x,t)}{\partial t}v(x)dx-a\int_{0}^{1}\frac{\partial^{2}u(x,t)}{\partial x^{2}}v(x)dx=\int_{0}^{1}f(x,t)v(x)dx$

若取 $v(x)\in H_{0}^{1}([0,1])$ ，从而 $v(0)=v(1)=0$ ，再利用分部积分，上式变成为

$\int_{0}^{1}\frac{\partial u(x,t)}{\partial t}v(x)dx+a\int_{0}^{1}\frac{\partial u(x,t)}{\partial x}v^{'}(x)dx=\int_{0}^{1}f(x,t)v(x)dx$

同样，由公式（1）第2个方程可得

$\int_{0}^{1}u(x,0)v(x)dx=\int_{0}^{1}\Psi(x)v(x)dx$

这样，就得到原问题的变分形式：

在任意时刻 $t\in (0,T]$ ，求函数 $u(x,t)\in H_{0}^{1}([0,1])$ （指的是在固定的时刻t，关于x的函数 $u(x,t)\in H_{0}^{1}([0,1])$ ），使得

$\int_{0}^{1}\frac{\partial u(x,t)}{\partial t}v(x)dx+a\int_{0}^{1}\frac{\partial u(x,t)}{\partial x}v^{'}(x)dx=\int_{0}^{1}f(x,t)v(x)dx,\; \forall v \in H_{0}^{1}([0,1]) \;\;(2)$

$\int_{0}^{1}u(x,0)v(x)dx=\int_{0}^{1}\Psi(x)v(x)dx,\;\forall v \in H_{0}^{1}([0,1])\;\;(3)$

二、有限元法进行空间半离散

首先对空间区域进行剖分离散，为简单起见，对[0,1]进行等距剖分，步长为h，节点坐标为 $x_{i}=ih,0\leqslant i\leqslant m\;\;and\;\;h=\frac{1}{m}$ 。设 $V_{h}$ 为 $H_{0}^{1}([0,1])$ 中的分片连续线性函数空间，其基函数为 $\phi_{0}(x),\phi_{1}(x),\cdots,\phi_{m}(x)\in V_{h}$ ，其定义式为

$\varphi_{i}(x)=\left\{\begin{matrix} \frac{x-x_{i-1}}{h},x_{i-1}\leqslant x \leqslant x_{i},\\ \frac{x_{i+1}-x}{h},x_{i}<x\leqslant x_{i+1}=\left\{\begin{matrix} 1-\frac{|x-x_{i}|}{h},x_{i-1}\leqslant x \leqslant x_{i+1},\\ 0,\;else \end{matrix}\right. 1 \leqslant i \leqslant m-1\\ 0,\;else \end{matrix}\right.$

且 $\varphi_{0}(x)=\left\{\begin{matrix} \frac{x_{1}-x}{h},x\in[x_{0},x_{1}]\\ 0,\;else \end{matrix}\right.,\varphi_{m}(x)=\left\{\begin{matrix} \frac{x-x_{m-1}}{h},x\in[x_{m-1},x_{m}]\\ 0,\;else \end{matrix}\right.$

可知有

$\varphi_{i}(x_{j})=\left\{\begin{matrix} 1,i=j,\\ 0,i\neq j. \end{matrix}\right.$

现在假设数值解具有可分离变量的性质，即

$u_{h}(x,t)=\sum_{i=0}^{m}\alpha_{i}(t)\phi_{i}(x)$

则原来的连续型问题（公式2、公式3）可以相应的离散为：

$\int_{0}^{1}\frac{\partial u_{h}(x,t)}{\partial t}v_{h}(x)dx+a\int_{0}^{1}\frac{\partial u_{h}(x,t)}{\partial x}v^{'}_{h}(x)dx=\int_{0}^{1}f(x,t)v_{h}(x)dx,$

$\int_{0}^{1}u_{h}(x,0)v_{h}(x)dx=\int_{0}^{1}\Psi(x)v_{h}(x)dx,\;\forall v \in V_{h}$

即

$\sum_{i=0}^{m}\alpha^{'}_{i}(t)\int_{0}^{1}\varphi_{i}(x)\varphi_{j}(x)dx+a\sum_{i=0}^{m}\alpha_{i}(t)\int_{0}^{1}\varphi^{'}_{i}(x)\varphi^{'}_{j}(x)dx=\int_{0}^{1}f(x,t)\varphi_{j}(x)dx,\;\;(4)$

$\sum_{i=0}^{m}\alpha_{i}(0)\int_{0}^{1}\varphi_{i}(x)\varphi_{j}(x)dx=\int_{0}^{1}\Psi(x)\varphi_{j}(x)dx,j=0,1,2,\cdots,m.\;\;(5)$

如果引入记号

$a_{i,j}=\int_{0}^{1}\varphi_{i}(x)\varphi_{j}(x)dx$ $b_{i,j}=\int_{0}^{1}\varphi^{'}_{i}(x)\varphi^{'}_{j}(x)dx$ $c_{j}=\int_{0}^{1}\Psi(x)\varphi_{j}(x)dx$

$F_{j}(t)=\int_{0}^{1}f(x,t)\varphi_{j}(x)dx$ , $j=0,1,2,\cdots,m\;\;\;\;(6)$

则有 $a_{i,j}=a_{j,i},b_{i,j}=b_{j,i}$ ，且公式（4）和公式（5）可以简写为

$\left\{\begin{matrix} A\mathbf{\alpha}^{'}(t)+aB\mathbf{\alpha}(t)=\mathbf{F}(t),\\ A\mathbf{\alpha}(0)=\mathbf{C}. \end{matrix}\right.\;\;\;(7)$

其中m+1方阵A，B定义为 $A=(a_{i,j}),B=(b_{i,j})$ ，且向量 $\mathbf{C},\mathbf{F}(t),\mathbf{\alpha}(t)$ 分别定义为

$\mathbf{C}=\begin{pmatrix} c_{0}\\ c_{1}\\ \vdots\\ c_{m-1}\\ c_{m} \end{pmatrix}$ $\mathbf{F}(t)=\begin{pmatrix} F_{0}(t)\\ F_{1}(t)\\ \vdots\\ F_{m-1}(t)\\ F_{m} (t)\end{pmatrix}$ $\mathbf{\alpha}(t)=\begin{pmatrix} \alpha_{0}(t)\\ \alpha_{1}(t)\\ \vdots\\ \alpha_{m-1}(t)\\ \alpha_{m}(t) \end{pmatrix}$

在结构力学中，A称为总质量矩阵，B称为总刚度矩阵，F（t）称为总荷载。

三、差分法进行时间全离散

经过有限元法离散后得到的离散结构式公式（7）只是在空间上进行了离散，时间变量依然是连续的，这种离散称为半离散，还需要对时间变量进行离散，这样最后得到的离散形式称为全离散。此处采用差分离散，即关于时间的导数用差商近似代替，于是需要对时间区域进行剖分，仍采用等距剖分，对时间区间[0,T]进行剖分。设时间步长为τ，节点坐标为 $t_{k}=k\tau,0\leqslant k\leqslant n\;\;\;and\;\;\; \tau=\frac{T}{n}$ 。对于半离散方程系统公式（7），在离散节点上应有 $A\mathbf{\alpha}^{'}(t_{\beta})+aB\alpha(t_{\beta})=\mathbf{F}(t_{\beta})$ 。

特别的：

（1）若取 $\beta=k,\alpha^{'}(t_{k})=\frac{\mathbf{\alpha}(t_{k+1})-\mathbf{\alpha}(t_{k})}{\tau}+O(\tau)$ ，用数值解 $\alpha^{k}=(\alpha_{0}^{k},\alpha_{1}^{k},\cdots,\alpha_{m}^{k})^{T}$ 代替精确解 $\alpha(t_{k})$ 并忽略高阶小项，可得向前欧拉法：

$A\frac{\mathbf{\alpha}^{k+1}-\mathbf{\alpha}^{k}}{\tau}+aB\mathbf{\alpha}^{k}=\mathbf{F}(t_{k})$ ，局部截断误差为 $O(\tau)$

（2）若取 $\beta=k+\frac{1}{2}$ ，其中

$t_{k+\frac{1}{2}}=\frac{t_{k}+t_{k+1}}{2}\;\;and\;\;\mathbf{\alpha}^{'}(t_{k+\frac{1}{2}})=\frac{\mathbf{\alpha}(t_{k+1})-\mathbf{\alpha}(t_{k})}{\tau}+O(\tau^{2})$

再利用数值解 $\alpha^{k}=(\alpha_{0}^{k},\alpha_{1}^{k},\cdots,\alpha_{m}^{k})^{T}$ 代替精确解 $\alpha(t_{k})$ 并忽略高阶小项，可得Crank-Nicolson方法：

$A\frac{\mathbf{\alpha}^{k+1}-\mathbf{\alpha}^{k}}{\tau}+aB\frac{\mathbf{\alpha}^{k+1}+\mathbf{\alpha}^{k}}{2}=\mathbf{F}(\frac{t_{k}+t_{k+1}}{2})$ ，局部截断误差为 $O(\tau^{2})$

（3）若取 $\beta=k+1,\alpha^{'}(t_{k+1})=\frac{\mathbf{\alpha}(t_{k+1})-\mathbf{\alpha}(t_{k})}{\tau}+O(\tau)$ ，再利用数值解 $\alpha^{k}=(\alpha_{0}^{k},\alpha_{1}^{k},\cdots,\alpha_{m}^{k})^{T}$ 代替精确解 $\alpha(t_{k})$ 并忽略高阶小项，可得向后欧拉法：

$A\frac{\mathbf{\alpha}^{k+1}-\mathbf{\alpha}^{k}}{\tau}+aB\mathbf{\alpha}^{k+1}=\mathbf{F}(t_{k+1})$ ，局部截断误差为 $O(\tau)$

整理以后可以得到全离散格式：

对 $k=0,1,\cdots,n-1,$

$(A+\theta a \tau B)\mathbf{\alpha}^{k+1}=(A-(1-\theta)a\tau B)\mathbf{\alpha}^{k}+\tau \mathbf{F}((1-\theta)t_{k}+\theta t_{k+1}) \;\;(8)$

取 $\theta =0,\frac{1}{2},1$ 分别对应向前欧拉有限元格式、Crank-Nicolson有限元格式、向后欧拉有限元格式。

由公式（7）第2式，上述格式都可以取初始值：

$\mathbf{\alpha}^{0}=\mathbf{\alpha}(0)=A^{-1}\mathbf{C}\;\;\;(9)$

采用追赶法求解公式（8）即可。

四、相关量的数值计算

1、 $a_{i,j}$ 的计算。

直接计算可得总质量矩阵

$A=(a_{ij})=\frac{h}{6}\begin{pmatrix} 2 & 1 & & 0 & \\ 1 & 4 & 1 & & \\ & \ddots & \ddots & \ddots & \\ & & 1 & 4 & 1\\ & 0 & & 1 & 2 \end{pmatrix}_{(m+1)\times(m+1)}$

2、 $b_{i,j}$ 的计算。

直接计算可得总刚度矩阵

$B=(b_{ij})=\frac{1}{h}\begin{pmatrix} 1 & -1 & & 0 & \\ -1 & 2 & -1 & & \\ & \ddots & \ddots & \ddots & \\ & & -1 & 2 & -1\\ & 0 & & -1 & 1 \end{pmatrix}_{(m+1)\times(m+1)}$

3、 $c_{j}$ 和 $F_{j}(t)$ 的计算。

可以利用数值节分公式 $\int_{a}^{b}g(x)dx \approx \frac{b-a}{2}(g(\frac{a+b}{2}-\frac{b-a}{2\sqrt{3}})+g(\frac{a+b}{2}+\frac{b-a}{2\sqrt{3}}))$ 近似计算得到。

4、边界条件的处理。

由于对连续问题有边界条件 $u(0,t)=u(1,t)=0$ ，离散情况下则为 $u_{h}(0,t)=u_{h}(1,t)=0$ ，从而利用

$u_{h}(x_{i},t)=\sum_{j=0}^{m}\alpha_{j}(t)\varphi_{j}(x)|_{x=x_{i}}=\alpha_{i}(t)$

得到 $\alpha_{0}(t)=\alpha_{m}(t)=0$ ，即有 $\alpha_{0}(t_{k})=\alpha_{m}(t_{k})=0,k=1,2,\cdots,n$ ，因此可以取边界条件为 $\alpha_{0}^{k}=\alpha_{m}^{k}=0,k=1,2,\cdots,n$ 。

五、编程时的说明

1、初始值。

上文中初始值取成公式（9），因此在实际计算中需要求解一个三对角线性方程组，计算比较复杂。事实上还可以这样考虑，由于初始条件为 $u(x,0)=\Psi(x)$ ，所以

$u_{h}(x_{i},0)=\sum_{j=0}^{m}\alpha_{j}(0)\varphi_{j}(x_{i})=\alpha_{i}(0)\approx\Psi(x)$

可以将 $\alpha_{i}^{0}$ 取成 $\Psi(x_{i}),i=0,1,\cdots,m$ 从而简化初始值的计算。

2、线性方程组的系数矩阵。

原则上线性方程组（公式 8）中的矩阵A,B和向量F的计算都可以利用小的剖分单元上的单元质量矩阵、单元刚度矩阵和单元荷载来合成。具体操作可以参考专栏第一篇文章内容。

有限元之初探门径-CSDN博客https://blog.csdn.net/L_peanut/article/details/139091636 3、数值格式（公式 8）是一个时间渐进格式，在第0个时间层，初始信息为 $\mathbf{\alpha}^{0}$ ，可以用公式（9）来计算，也可以直接取 $\alpha_{i}^{0}$ 为 $\Psi(x_{i}),i=0,1,\cdots,m$ ，然后利用公式（8）计算第1个时间层上的信息 $\mathbf{\alpha}^{1}$ ，但是在计算 $\mathbf{\alpha}^{1}$ 时，由于边界信息已知，即 $\mathbf{\alpha}_{0}^{1}=\alpha_{m}^{1}=0$ ，所以需要修改线性方程组（公式 8）左边的系数矩阵 $A+\theta a \tau B$ ，使它的第0行元素和第0列元素除了第0行第0列这一个元素取成1以外，其它元素都取成0。同样，第m行元素和第m列元素除了第m行第m列这一个元素取成1以外，其它元素都取成0。最后修改公式（8）右边的列向量 $(A-(1-\theta)a\tau B)\mathbf{\alpha}^{k}+\tau \mathbf{F}((1-\theta)t_{k}+\theta t_{k+1})$ ，使其第0个元素及第m个元素均为0。这样，线性方程组（公式 8）求解完以后可以保证边界条件 $\mathbf{\alpha}_{0}^{1}=\alpha_{m}^{1}=0$ 成立，这样得到完整的第1个时间层上的信息 $\mathbf{\alpha}^{1}$ 。然后按照同样的方法继续计算后面第2,3，…，n个时间层上的信息。最终，可以得到数值解在第k，k=0,1，…，n个时间层上的近似式

$u_{u}(x,t_{k})=\sum_{j=0}^{m}\alpha_{j}(t_{k})\varphi_{j}(x)\approx\sum_{j=0}^{m}\alpha^{k}_{j}\varphi_{j}(x)$

显然，有限元离散后x方向每个位置都有数值解，但差分离散后t方向只在节点处有信息。

六、算例实现

用有限元求解抛物型方程的初边值问题：

$\left\{\begin{matrix} \frac{\partial u}{\partial t}-2\frac{\partial ^{2}u}{\partial x^{2}}=0,0<x<\pi,0<t\leqslant1,\\ u(x,0) =5sin(2x)-30sin(3x),0\leqslant x\leqslant \pi,\\ u(0,t)=u(1,t)=0,0<t\leqslant 1. \end{matrix}\right.$

已知其精确解为 $u(x,t)=5e^{-8t}sin(2x)-30e^{-18t}sin(3x)$ 。分别取第一种剖分数m=n=16和第二种剖分数m=n=32，输出在节点 $(\frac{i\pi}{8},\frac{1}{2}),i=1,2,\cdots,7$ 处的数值解并给出误差。

6.1 C++代码

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

int m,n;
double h,tau,T,*xco,**alpha;
double gauss=0.5773502692;   //数值积分中的高斯点

int main(int argc, char *argv[])
{
        int i,k;
        double a1,*t,*a,*b,*c,*d,*ans,tmid,temp1,temp2;
        double f(double x, double t);
        double phi(int i, double x);
        double psi(double x);
        double fun1(int i, double x, double t);
        double integral(double a, double b, int i, double t, double(*fun)(int, double, double));
        double exact(double x, double t);
        double *chase_algorithm(int n, double *a, double *b, double *c, double *d);

        m=16;
        n=16;
        printf("m=%d, n=%d.\n",m,n);

        h=pi/m;    //空间步长
        T=1.0;    //时间终点
        tau=T/n;   //时间步长
        a1=2.0;   //热的扩散率,原方程中的常系数a

        xco=(double *)malloc(sizeof(double)*(m+1));   //位移x方向上的剖分节点坐标
        for(i=0;i<=m;i++)
                xco[i]=i*h;

        t=(double *)malloc(sizeof(double)*(n+1));   //时间t方向上的剖分节点坐标
        for(k=0;k<=n;k++)
                t[k]=k*tau;

        alpha=(double **)malloc(sizeof(double *)*(m+1));   //设置二维数组alpha(i,k), i表示空间分量,k表示时间分量
        for(i=0;i<=m;i++)
                alpha[i]=(double *)malloc(sizeof(double)*(n+1));
        for(i=0;i<=m;i++)
                alpha[i][0]=psi(xco[i]);   //用简化的方法计算初始值

        a=(double *)malloc(sizeof(double)*(m+1));   //计算线性方程组的系数矩阵（三对角）
        b=(double *)malloc(sizeof(double)*(m+1));
        c=(double *)malloc(sizeof(double)*(m+1));
        temp1=h/6.0-a1*tau/(2*h);
        temp2=2*h/3.0+a1*tau/h;
        for(i=0;i<=m;i++)
        {
                a[i]=temp1;   //下次对角线上的元素
                b[i]=temp2;   //主对角线上的元素
                c[i]=temp1;   //上次对角线上的元素
        }
        b[0]=b[0]/2.0;
        b[m]=b[m]/2.0;

        d=(double *)malloc(sizeof(double)*(m+1));   //计算线性方程组的右端项
        for(k=0;k<n;k++)
        {
                tmid=(t[k]+t[k+1])/2.0;   //计算向量tau*F(tmid)
                for(i=1;i<m;i++)
                {
                        d[i]=integral(xco[i-1],xco[i],i,tmid,fun1)+integral(xco[i],xco[i+1],i,tmid,fun1);
                }
                d[0]=integral(xco[0],xco[1],0,tmid,fun1);
                d[m]=integral(xco[m-1],xco[m],m,tmid,fun1);

                temp1=2.0*h/3.0-a1*tau/h;
                temp2=h/6.0+a1*tau/(2*h);
                //计算(A-a*tau*B/2.0)*alpha[k]
                for(i=1;i<=m-1;i++)
                {
                        d[i]=d[i]+temp2*(alpha[i-1][k]+alpha[i+1][k])+temp1*alpha[i][k];
                }
                d[0]=d[0]+temp1*alpha[0][k]/2.0+temp2*alpha[1][k];
                d[m]=d[m]+temp2*alpha[m-1][k]+temp1*alpha[m][k]/2.0;

                //处理零边界条件
                d[0]=0.0; d[m]=0.0;
                b[0]=1.0; c[0]=0.0;
                a[1]=0.0; a[m]=0.0;
                c[m-1]=0.0; b[m]=1.0;

                ans=chase_algorithm(m+1,a,b,c,d);
                for(i=0;i<=m;i++)
                        alpha[i][k+1]=ans[i];
                free(ans);
        }

        free(a);free(b);free(c);free(d);
        k=m/8;
        for(i=k;i<m;i=i+k)
        {
                printf("alpha[%d][%d]=%f, err=%.4e.\n", i, n/2, alpha[i][n/2], fabs(exact(xco[i], T/2.0)-alpha[i][n/2]));
        }
        for(i=0;i<=m;i++)
                free(alpha[i]);
        free(xco);free(alpha);free(t);

        return 0;
}



//右端项函数f(x,t)
double f(double x, double t)
{
        return 0.0;
}
//基函数
double phi(int i, double x)
{
        double temp, z;

        temp=fabs(x-xco[i]);
        if(temp<=h)
                z=1.0-temp/h;
        else
                z=0.0;

        return z;
}
double psi(double x)
{
        return 5*sin(2*x)-30*sin(3*x);
}
//计算单元荷载时的被积函数
double fun1(int i, double x, double t)
{
        return f(x,t)*phi(i,x);
}
//在区间[a,b]上对被积函数fun(i,x,t)进行数值积分（两点高斯公式）
double integral(double a, double b, int i, double t, double(*fun)(int, double, double))
{
        double mid, w, ans;
        mid=(b+a)/2.0;
        w=(b-a)/2.0;
        ans=w*((*fun)(i,mid+w*gauss,t)+(*fun)(i,mid-w*gauss,t));

        return ans;
}
double exact(double x, double t)
{
        return 5*exp(-8.0*t)*sin(2*x)-30*exp(-18.0*t)*sin(3*x);
}
//追赶法
double * chase_algorithm(int n, double *a, double *b, double *c, double *d)
{
        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;
}

6.2 计算结果

（1）当m=n=16时，计算结果如下：

m=16, n=16.
alpha[2][8]=0.055488, err=5.8471e-03.
alpha[4][8]=0.078932, err=1.0029e-02.
alpha[6][8]=0.056490, err=9.6821e-03.
alpha[8][8]=0.000767, err=2.9351e-03.
alpha[10][8]=-0.055903, err=7.4356e-03.
alpha[12][8]=-0.080017, err=1.4179e-02.
alpha[14][8]=-0.056905, err=1.1271e-02.

（2）当m=n=32时，计算结果如下：

m=32, n=32.
alpha[4][16]=0.060101, err=1.2341e-03.
alpha[8][16]=0.086607, err=2.3533e-03.
alpha[12][16]=0.063613, err=2.5593e-03.
alpha[16][16]=0.002688, err=1.0143e-03.
alpha[20][16]=-0.061556, err=1.7830e-03.
alpha[24][16]=-0.090408, err=3.7876e-03.
alpha[28][16]=-0.065068, err=3.1082e-03.