偏微分方程算法之向前欧拉法（Forward Euler）

猿核试Bug愁

已于 2024-04-09 16:22:24 修改

阅读量2.3k

点赞数 29

分类专栏：偏微分方程（抛物型）文章标签：算法

于 2024-04-09 16:12:10 首次发布

本文链接：https://blog.csdn.net/L_peanut/article/details/137541192

版权

偏微分方程（抛物型）专栏收录该内容

15 篇文章 12 订阅 ¥9.90 ¥99.00

订阅专栏

超级会员免费看

本文介绍了使用向前欧拉法求解一维非齐次热传导方程的初边值问题。通过网格剖分、离散方程建立、差分格式以及算例实现，详细阐述了该方法的理论推导和数值求解过程，揭示了稳定性条件对计算结果的影响。

摘要由CSDN通过智能技术生成

一、研究对象

以一维（空间维度）非齐次热传导方程的定解问题研究抛物型方程的差分方法为研究对象（热传导方程描述的是存在热源的一个区域内，与周围介质存在热交换的物体内部温度的分布及温度如何随时间变化）。

研究对象为抛物型方程初边值问题：

$\left\{\begin{matrix} \frac{\partial u(x,t)}{\partial t}-a\frac{\partial^{2}u(x,t)}{\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(1)\\ u(0,t)=\alpha(t), \space\space\space\space u(1,t)=\beta(t),\space\space 0<t\leqslant T \end{matrix}\right.$

其中，a>0为常数。

二、理论推导

类比常微分方程差分法求解，公式（1）的差分法可以通过多步实现：

2.1 网格剖分

对二维求解区域 $\Omega ={(x,t)|0\leqslant x\leqslant 1,0\leqslant t}$ 作矩形网格剖分，采取空间域、时间域的等距剖分。将空间[0,1]等分m份，节点为 $x_{i}=0+ih$ ，步长h=1/m。将空间域[a,b]等分m份，节点为 $x_{i}=a+ih,0 \leqslant i\leqslant m$ 且 $h=(b-a)/m$ 。对时间域[0,T]等分n份，节点为 $t_{k}=0+k\tau ,o\leqslant k\leqslant n$ 且 $\tau =T/n$ 。于是得到 $(m+1)*(n+1)$ 个网格节点 $(x_{i},t_{k}),0\leqslant i\leqslant m,0\leqslant k\leqslant n$ 。如图所示：

2.2 建立离散方程

在剖分好的网格节点上建立离散方程。本质上是将求解域内的微分方程弱化为节点上的离散方程，即：

$\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}), \space\space 1\leqslant i\leqslant m-1,\space\space 0<k\leqslant n,\\ u(x_{i},t_{0})=\varphi (x_{i}),\space\space\space\space 0\leqslant i\leqslant m,\space\space\space\space (2)\\ u(x_{0},t_{k})=\alpha(t_{k}),\space\space u(x_{m},t_{k})=\beta(t_{k}),\space\space 0<k\leqslant n \end{matrix}\right.$

2.3 建立差分格式

由一阶向前差商公式：

$\left\{\begin{matrix} \frac{\partial u}{\partial x}(x_{i},y_{j})\approx \frac{u(x_{i+1},y_{j})-u(x_{i},y_{j})}{h_{x}}\\ \frac{\partial u}{\partial y}(x_{i},y_{j})\approx \frac{u(x_{i},y_{j+1})-u(x_{i},y_{j})}{h_{y}} \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})}{\Delta t}\approx \frac{u(x_{i},t_{k+1})-u(x_{i},t_{k})}{\tau } \space\space\space\space (3)$

由二阶偏导数的二阶中心差商公式：

$\left\{\begin{matrix} \frac{\partial^{2}u}{\partial x^{2}}(x_{i},y_{j})\approx \frac{u(x_{i+1},y_{j})-2u(x_{i},y_{j})+u(x_{i-1},y_{j})}{h^{2}_{x}}\\ \frac{\partial^{2}u}{\partial y^{2}}(x_{i},y_{j})\approx \frac{u(x_{i},y_{j+1})-2u(x_{i},y_{j})+u(x_{i},y_{j-1})}{h^{2}_{y}} \end{matrix}\right.$

得到关于空间的二阶偏导数的中心差商公式：

$\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})}{\Delta x^{2}}=\frac{u(x_{i+1},t_{k})-2u(x_{i},t_{k})+u(x_{i-1},t_{k})}{h^{2}}$

将关于时间的一阶偏导数的差商形式公式、关于空间的二阶偏导数的中心差商公式代入公式（2），可得：

$\frac{u(x_{i},t_{k+1})-u(x_{i},t_{k})}{\tau}-a\frac{u(x_{i+1},t_{k})-2u(x_{i},t_{k})+u(x_{i-1},t_{k})}{h^{2}}=f(x_{i},t_{k})+C(\tau+h^{2})$

其中，C为常数，满足 $C=\underset{(x,t)\in \Omega}{max} \begin{Bmatrix} \frac{\partial^{2}u(x,t)}{\partial t^{2}},\frac{\partial^{4}u(x,t)}{\partial x^{4}} \end{Bmatrix}$ 。利用数值解 $u^{k}_{i}$ 代替精确解 $u(x_{i},t_{k})$ 并忽略高阶项，即可得到向前欧拉格式：

$\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<k\leqslant n-1, \\ u^{0}_{i}=\varphi (x_{i}),\space\space\space\space 0\leqslant i\leqslant m,\space\space\space\space (4) \\ u^{k}_{0}=\alpha(t_{k}),\space\space u^{k}_{m}=\beta(t_{k}),\space\space 1<k\leqslant n \end{matrix}\right.$

向前欧拉格式的截断误差为 $O(\tau+h^{2})$ ，整理公式（4）后得：

$\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}) ,\space\space 1\leqslant i \leqslant m-1,\space\space 0\leqslant k \leqslant n-1,\\ u^{0}_{i}=\varphi (x_{i}),\space\space 0\leqslant i \leqslant m, \space\space\space\space (5)\\ u^{k}_{0}=\alpha(t_{k}),\space\space u^{k}_{m}=\beta(t_{k}),\space\space 0<k\leqslant n \end{matrix}\right.$

其中 $r=\frac{a\tau}{h^{2}}$ 称为网比，是与时间、空间相关的一个值。

2.4 求解差分格式

公式（5）显示第k+1个时间层 $t=t_{k+1}$ 上的数值解 $u^{k+1}_{i}(1 \leqslant i\leqslant m-1)$ 可以由第k层上的已知信息 $u^k_{i-1},u^{k}_{i},u^{k}_{i+1}$ 表示，这是一个时间层进（Time marching），如图所示：

意味着可以通过第0层上的初始信息 $u^{0}_{i},0 \leqslant i \leqslant m$ ，利用公式（5）中第1式计算出第1层上的信息 $u^{1}_{i},1 \leqslant i \leqslant m-1$ ，结合边界条件 $u^{1}_{0},u^{1}_{m}$ 已知（即公式5中第2、3式），即第1层上所有信息 $u^{1}_{i},0 \leqslant i \leqslant m$ 便全部已知。再通过公式（5）中的第1、3式计算第2层的信息，如此反复，即可得到所有网格点的信息，即获得所有网格点的数值解。

可将公式（5）改写为矩阵形式，即：

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

三、算例实现

使用向前欧拉法计算抛物型方程初边值问题：

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

已知精确解为 $u(x,t)=x(x^{2}+e^{t})$ 。分别取 $h_{1}=\frac{1}{5},\tau_{1}=\frac{1}{100},h_{2}=\frac{1}{5},\tau_{2}=\frac{1}{200},h_{3}=\frac{1}{10},\tau_{3}=\frac{1}{200},h_{4}=\frac{1}{10},\tau_{4}=\frac{1}{100}$ 。列出节点 $(0.4,0.2i),i=1,2,3,4,5$ 处的结果及误差。

代码如下：

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

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

        n=100;  //时间域n等分
        m=5;    //空间域m等分
        a=1.0;
        h=1.0/m;  //空间步长
        tau=1.0/n;  //时间步长
        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(i=0;i<=n;i++)
                t[i]=i*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(i=1;i<=n;i++)
        {
                u[0][i]=alpha(t[i]);
                u[m][i]=beta(t[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]);
                }
        }

        j=int(0.2/tau);
        number=int(0.4/h);
        for(k=j;k<=n;k=k+j)
        {
                printf("(x,t)=(%.1f,%.1f), y=%f, exact=%f, err=%.4e.\n",x[number],t[k],u[number][k],exact(x[number],t[k]),fabs(u[number][k]-exact(x[number],t[k])));
        }

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

        return 0;
}


double phi(double x)
{
        return x*x*x+x;
}

double alpha(double t)
{
        return 0.0;
}
double beta(double t)
{
        return 1.0+exp(t);
}
double f(double x, double t)
{
        return x*exp(t)-6*x;
}
double exact(double x, double t)
{
        return x*(x*x+exp(t));
}

当m=5,n=100,r=1/4时，运行结果如下：

r=0.2500.
(x,t)=(0.4,0.2), y=0.552290, exact=0.552561, err=2.7087e-04.
(x,t)=(0.4,0.4), y=0.660358, exact=0.660730, err=3.7148e-04.
(x,t)=(0.4,0.6), y=0.792388, exact=0.792848, err=4.5918e-04.
(x,t)=(0.4,0.8), y=0.953655, exact=0.954216, err=5.6158e-04.
(x,t)=(0.4,1.0), y=1.150627, exact=1.151313, err=6.8601e-04.

当m=5,n=200,r=1/8时，运行结果如下：

r=0.1250.
(x,t)=(0.4,0.2), y=0.552427, exact=0.552561, err=1.3428e-04.
(x,t)=(0.4,0.4), y=0.660545, exact=0.660730, err=1.8521e-04.
(x,t)=(0.4,0.6), y=0.792618, exact=0.792848, err=2.2921e-04.
(x,t)=(0.4,0.8), y=0.953936, exact=0.954216, err=2.8038e-04.
(x,t)=(0.4,1.0), y=1.150970, exact=1.151313, err=3.4252e-04.

当m=10,n=200,r=1/2时，运行结果如下：

r=0.5000.
(x,t)=(0.4,0.2), y=0.552426, exact=0.552561, err=1.3555e-04.
(x,t)=(0.4,0.4), y=0.660544, exact=0.660730, err=1.8589e-04.
(x,t)=(0.4,0.6), y=0.792618, exact=0.792848, err=2.2978e-04.
(x,t)=(0.4,0.8), y=0.953935, exact=0.954216, err=2.8102e-04.
(x,t)=(0.4,1.0), y=1.150969, exact=1.151313, err=3.4329e-04.

当m=10,n=100,r=1时，运行结果如下：

r=1.0000.
(x,t)=(0.4,0.2), y=-283.280312, exact=0.552561, err=2.8383e+02.
(x,t)=(0.4,0.4), y=-609165439303.370361, exact=0.660730, err=6.0917e+11.
(x,t)=(0.4,0.6), y=-1115707124224570556416.000000, exact=0.792848, err=1.1157e+21.
(x,t)=(0.4,0.8), y=-2008676972462646683560975532032.000000, exact=0.954216, err=2.0087e+30.
(x,t)=(0.4,1.0), y=-3608553481656707469335928453820848275456.000000, exact=1.151313, err=3.6086e+39.

从四种不同时空步长的计算结果可知，前三种情况的计算都可以收敛，最后一种步长的网比不满足稳定性条件，导致误差会指数级增长，无法收敛。同时，在空间步长不变的情况下，时间步长缩小一半，最后结果的误差也会缩小一半；在时间步长不变的情况下，空间步长缩小一半，误差结果无明显变化；如果时间步长、空间步长都缩小一半，最后的结果误差也会缩小一半。