Euler法原理
Euler法是一般数值分析课程中最基础也最重要的一种数值方法,用于求解微分方程的初值问题。一阶常微分方程的初值问题一般形式是
数值方法就是要求问题(1)的解
在若干离散点的近似值
,(n=1,2,...)的方法,建立数值解法,首先要将微分方程离散化,转化为差分方程,用数值积分公式或泰勒展开法
由此得到迭代公式
用差分方程初值问题的解来近似微分方程初值问题(1)的解, 即由(3)依次算出
的近似值
,(n=1,2,...)。这组计算式称为
向前 Euler 法。
向前 Euler 法是显式一阶方法,本质上是利用第N点的函数值和一阶导数(斜率)去近似求解第n+1点的函数值,因此它的精度不高,好处是不需要迭代求解。
如果能够利用第N+1点的一阶导数来修正误差,就能够得到更高的精度。直观上用梯形公式计算数值积分要比矩形公式好,梯形公式利用第N+1点的斜率,与第N点的斜率取平均得到等效的斜率,具有二阶精度,但是它是隐式格式,一般需用迭代法解方程
如果实际计算时精度要求不太高,用公式(3)求解时,每步可以只迭代一次,则得到改进的Euler 法,相当于将Euler公式与梯形公式结合使用:先用 Euler 公式求
的一个初步近似值
(预测值),然后用梯形公式校正求得近似值
。为了编程方便,把斜率单独提取出来计算等效斜率,迭代式如下
改进Euler法具有二阶精度,但收敛率略小于梯形法。
二维热传导方程
二维热传导方程(常物性)具有如下形式:
是热扩散系数,
为热源项,初值条件
。采用Euler法对时间项进行离散化:
Euler法只处理了时间偏导项(向前差分),进行接下来对空间二阶偏导进行离散化。用
表示第i行(y方向)第j列(x方向)的温度值,t时刻的温度场以(Ny+1,Nx+1)大小的矩阵存储,
代表x,y方向的分段数量,每一段的等分长度为
、
矩阵的每个值代表一个网格节点温度,注意y方向代表行标,x方向代表列标
每一点对
,
方向均采用中心差分格式:
把二阶偏导组装成矩阵形式,维度与
相同
写成矩阵乘积形式,记为
和
称为差分系数矩阵,均为三对角对称矩阵,但是大小不同,
为(Nx+1,Nx+1),
为(Ny+1,Ny+1)。这样迭代公式变为
这里需要注意式(9-1)和(9-2)中求得的二阶偏导矩阵的边界点是无效的,因为差分系数矩阵第一行和最后一行是不完整的[-2 1]和[1 -2],相当于边界点的温度恒为0,因此需要单独给定边界条件才能确定矩阵边界的温度值。
Python计算代码
因为给出了矩阵化的迭代式,所以编程就很简单了,直接采用numpy的向量化运算。计算10秒,前5秒加热,后5秒不加热,边界条件采用第一类边界,即给定温度值。
import
由于向前Euler法只有一阶精度,当
太大或
太小均会造成函数不收敛,得到错误数值结果。当Nt<2000,即dt>0.005时出现不收敛现象。
由于Euler法的向量化运算十分简单快速,实际操作中可以用较大的时间分段Nt,去增加总的时间精度而不会明显降低求解速度,当需要更高精度和速度时可以使用改进Euler法或Runge-Kutta法,工程上常用四阶精度的Runge-Kutta法求解。
总结
Euler法是显式一阶方法,可以快速求解ODE初值问题。
热传导方程本质上是PDE方程,在采用差分离散化空间变量后转换为时间变量的ODE方程,方可采用Euler法求解,其最终格式与热传导显式差分格式完全相同。
在处理温度U节点及其空间二阶导时采用向量化和矩阵化形式,避免编程中使用多重for-loop,而是采用向量化运算提高效率。
一般小规模问题,可用较大的时间分段Nt去提高收敛性和精度,而不必采用复杂的迭代式。
下一篇文章将讲解如何处理第二类边界条件和非矩形的边界形状。