fortran用adi方法解二维热方程_Python数值优化：使用Euler法求解二维热传导方程1...

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去提高收敛性和精度，而不必采用复杂的迭代式。

下一篇文章将讲解如何处理第二类边界条件和非矩形的边界形状。