有限差分法.docx

有限差分法
例题1：求常微分方程
(1,1)
定解条件为，
将求解区域离散为个点。子区间长度，分点，由

忽略余项，将差商代替微商，用近似表示，记，对内点建立差分方程，
(1,2)
由边界值确定。整理后得到关于的方程组，
(1,3)
为三对角方程组，可使用追赶法求解。
为验证计算结果，求出的精确解
(1,4)
并引入相对误差，
(1,5)
标准误差，
(1,6)
取，计算结果如下，

图(1,1) 方程(1,1)差分解及精确解结果

图(1,2) 差分法求解结果相对误差分布

由图(1,1)可见差分法结果与精确解基本吻合。图(1,2)显示，在边界附近误差较小，在求解区域内部误差较大。最大误差不足，总体标准误差，满足一般工程精度需求。

例题2；如图(2,1)厚度为的无限大平板，其热扩散系数为，初始温度为，求在边界温度作用下，平板内温度场变化。

图(2,1) 例题2问题描述

问题控制方程为，
(2,1)
边界条件
(2,2)
初始条件
(2,3)
同样先求出精确解和差分数值解。其精确解，
(2,4)
其中
差分法求解。对时间及空间离散为有限个差分点。对时间取向后差，对空间取中心差。差分格式，
(2,5)
其中为时间步长，为空间步长。
取点如下。

图(2,2) 例2差分格式取点示意

对内点使用格式，边界点值有边界条件得到。除层外，每一时间层由上一层隐式求得。探讨误差传播及总体误差情况，计算相对误差和标准误差
相对误差
(2,6)
标准误差
(2,7)
其中，为空间离散点数，为时间离散点数。
求解结果如下，

图(2,3) 温度场随时间变化示意 图(2,4) 误差随时间步、空间步变化

为探究精确解、数值解和误差在空间分布上的关系，将三者画在同一幅图中。将每一时刻的结果连起来可构成一幅动态图（见附件）。现分别取结果如下

图(2.5) t=0.09s计算结果及误差分布

图(2.6) t=1.07s计算结果及误差分布

图(2.7) t=1.19s计算结果及误差分布

图(2.3)为差分结果，实际上与精确解结果图别无二致。场内无热源（井），初始温度在边界温度的作用下，最终使场内各点温度区域边界温度。
误差分布相关讨论：
1. 观察图(2,4)与图(2,5)。将空间划分三等份，在初始阶段，两边的三分之一误差明显大于中间三分之一。两边的误差可能是由于求解精确解时级数展开项数不足造成的。
2. 观察图(2.4)，比较图(2,5)和图(2,6)。发现随时间步增加，误差在增加。
3. 观察图(2.4)、图(2,6)和图(27)。发现与例1相同结论，误差在边界附近较小，在区域中部误差较大。同时，随时间步增加，误差并不是无限制增长。显然，中部误差有下降趋势。
4. 标准误差为0.0110，满足一般工程要求。

例题3：矩形域的拉普拉斯方程边值问题
(3,1)
边界条件
(3,2)
将求解区域离散为有限个差分点，采用五点差分格式，
(3,3)
其中为方向的步长，为方向的步长。

图(3,1) 求解区域及边界条件 图(3,2) 五点差分格式示意

对内点使用格式，边界点值由边界条件直接获得，构成线性方程组。该线性方程组每个方程最多含有5个未知数且对角占优，可使用超松弛法求解，本次求解松弛因子取1.6。
方程精确解为
(3,4)
取，，计算结果如下。

图(3.3) 拉普拉斯方程解

为探究各差分点的误差关系，求出各点误差。同时由于精确解较多点趋于零，使用相对误差会出现不规则。故使用绝对误差，定义如下，
(3,5)
以及探究总体误差情况，计算标准误差
(3,6)
其中为方向节点数，为方向节点数
各差分点绝对误差分布如下

图(3.4) 各节点绝对误差分布

观察图(3,4)得，在边界附近，误差较小。在精确值趋于零的节点，误差较小。计算得标准误差，满足一般工程需求。

例题4：四边简支平板，受均布力荷载，求平板挠度。
其控制方程为
(4,1)

图(4,1) 四边简支板示意
其中，与板厚、材料常数有关。
差分格式取点如图(4,2)所示，其中为方向的步长，为方向的步长。对0点式各项使用下列公式代替，
(4,2)
(4,3)
图(4,2) 差分格式示意
(4,4)
将-代入得出差分格式
(4,5)
对于正则内点可直接使用格式建立方程。对于边界点，由于简支边挠度为零，已知，不需要建立方程。对于非正则内点（如6,3,7点在边界上，0点为非正则内点），要建立一个虚节点（如11号点），继续使用格式建立方程。在虚节点应使用简支变弯矩为零的边界条件建立方程。如
(4,6)
又因为，即边界处
(4,7)
为探究本方法误差情况，计算精确值进行比较。精确值如下
(4,8)
由于边界上挠度趋于零，不宜使用相对误差，故对各节点使用绝对误差，定义如下
(4,9)
其中为沿方向的节点编号。
同前例题计算标准误差
(4,10)
取的计算模型，荷载为。的计算网格，计算结果如下

图(4,3) 均布力作用下四边简支板挠度

图(4,4) 各节点绝对误差

如图(4,3)板受四边的约束，边界附近的挠度较小。在四角上，挠度受边界约束更明显。板中部挠度最大，约为2.5厘米。观察图(4,4)边界上误差最小。最大误差并不在区域中心部位。从边界向区域内部，误差先变大后缩小，大约经过10个节点，误差达到最大值。在四角上误差传播较慢。计算得标准误差，满足一般工程需求。
由上面四个例题可以看出，使用差分格式，由于受边界条件的约束，在边界上误差较小。边界向计算区域内部，误差逐渐变大。若选用的差分格式收敛且稳定，误差有上限。

附件：
changweifenfagncheng.f90 例题1计算程序
huitu4_1.m 例题1结果绘图代码
ReChuanDaoFangCheng.f90 例题2计算程序
huitu4_2.m 例题2结果绘图代码
resultgif.gif 例题2结果动态图
LaplaceEquation.f90 例题3计算程序
huitu4_3.m 例题3结果绘图代码
BanQiao.f90 例题4精确解计算程序
banqiao.m 例题4差分解计算及结果绘图程序