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

您所在的位置：网站首页 › 二维热传导方程数值解法 › 偏微分方程算法之向前欧拉法（Forward Euler）

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

2024-07-05 14:44| 来源: 网络整理| 查看: 265

一、研究对象

二、理论推导

2.1 网格剖分

2.2 建立离散方程

2.3 建立差分格式

2.4 求解差分格式

三、算例实现

一、研究对象

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

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

$\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 0x1, \space\space 0t\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 0t\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 0k\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 0k\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 0k\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 1k\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 0k\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 0x1, \space\space 0t \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 0t\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 #include #include 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

【本文地址】

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

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

今日新闻

推荐新闻