第二章 有限差分法初步

《第二章 有限差分法初步》由会员分享，可在线阅读，更多相关《第二章 有限差分法初步（61页珍藏版）》请在金锄头文库上搜索。

1、一 差商与微商 第二章 有限差分法初步 1有限差分法基本概念 i 有限差分的数学基础是用差商代替微商 有如下两种数学形式 i 微商 导数 的定义 是连续函数 则它的导数为 若 2 1 式 2 1 右边 是有限的差商 与 都不为零 而式 2 1 左边 是 当 趋于零时极限情形下的差商 称之微商 在 没有到达零之前 只是 的近似 趋于 的过程认为是近似向精确过渡 用 代替 就是精确向近似过渡 两者的差值 表示差商代替微商的偏差 ii 偏差 Taylor级数展开 2 2 稍加整理后可写成 2 3 可见与只能是近似相等 偏差为 iii 微商与差商的几何意义 图2 1差商与微商的比较 图2 1表示了差商

2、与微商之间的关系 应当指出 用不同方法得到的差商去代替微商 它们带来的偏差是不同的 向右 前 差商 iV 差商的几种表示 2 4 向左 后 差商 2 5 中心差商 取向右差商与向左差商的平均值 2 6 偏差分析 将Taylor级数写成 2 7 Taylor级数还可写成 2 8 由式 2 7 可得 由式 2 8 可得 2 9 2 10 2 9 2 10 得到 2 11 比较式 2 9 2 10 2 11 可看到 用不同的差商形式去代替微商 所带来的偏差是不同的 这些偏差都是截去了Taylor级数展开式中的高阶项而引起的 常称 截断误差 用向右差商与向左差商代替微商 其截断误差为与 同量级的小量

3、同量级的小量 中心差商的截断误差小于向右差商或向左差商 而用中心差商代替微商 其截断误差是与 讨论 上述一阶差商一般仍是x的函数 对它们还可以求差商 这种一阶差商的差商称为二阶差商 它是二阶微商的近似 常用向右差商的向左差商来近似二阶微商 即 V 二阶差商 根据式 2 7 2 8 可得 2 12 由式 2 12 知 二阶差商的截断误差也为与同阶的小量 结论 由于用差商代替微商必然带来截断误差 相应地用差分方程代替微分方程也必然带来截断误差 这是有限差分法固有的 因此 在应用有限差分法进行数值解时 必须对差分的构成及其对方程造成的误差引起注意 2从微分形式出发的差分格式 图2 2给出了一个简单边

4、界值问题 图2 2矩形区域离散化 问题是求图2 2所示的边值问题的解 其数学表达如下 方程 2 13 边界条件 2 14 2 15 2 16 2 17 式 2 13 2 17 构成定解问题 在问题的提法已经明白之后 差分格式的构成 下面分别予以说明 vi 构成差分格式 ii 建立区域内差分方程 iii 边界条件的差分形式 i 区域离散法 1 区域离散化所谓离散化 就是把几何上连续的区域用一系列网格线把它划分开 一般说来 网格形式应视几何区域的不同而不同 对于矩形区域而言 用矩形的网格 如图2 2 用五条水平网线与五条垂直网线把矩形区域离散掉 网线与网线的交点称之为 节点 节点与节点的距离称之为

5、步长 x方向的步长表示为 y方向的步长表示为 节点编号 为便于计算 需对节点逐个编号 常用 i j 表示节点位置 其中 i j是与网线相对应的正整数 i j的排列 可有不同的方式 习惯上 与x y轴相一致 i由左而右逐个增长 j由下而上逐个增长 但也有 考虑到与矩阵的格式相一致 i表示行数 由上而下逐个增长 j表示列数 由左而右逐个增长 这种从上到下 从左到右的编排与一般书写习惯也是一致的 因此 在计算机上算题也常被采用 在本课程中 大都采用与坐标相一致的编排方法 在区域内的节点称 内节点 在边界上的节点称 边界节点 图2 2所示边界是规则的 则节点或在区域内 或正好落在边界上 步长或可以是不

6、变的常量 即等步长 也可以在区域内的不同处是不同的 即变步长 如果区域内各处的温度梯度变化很大 则在温度变化剧烈处 网格布得密些 在温度变化不剧烈处 网格布得疏些 至于网格布置多少 步长取多大为宜 要根据具体问题 兼顾到计算的精确度与计算的工作量等因素而定 步长 从物理方面对区域离散化可作这样的理解 即认为区域内离散的每个节点 都集中着它周围区域 尺度为步长 的热容 或者说 区域内连续分布的热容都被分别地集中到离散的节点上去了 这样 节点的温度代表着它周围区域的某种平均温度 一系列离散的节点温度值代表着连续区域内的温度分布 区域离散化物理理解 节点 i j 处的温度表示成 2 差分方程代替微分

7、方程在上节我们已对有限差分法的数学基础作了简要的介绍 说明了如何用差商代替微商 以及由此带来的误差 这里介绍用差商代替微商的办法来处理导热方程 2 13 得到相应的差分方程 方程 2 13 对区域内各个点都成立的 当然对任意一个内节点 i j 也成立 2 13 在 i j 处存在二阶偏微商与这些二阶偏微商所对应的差商可表示成 i 2 3 4 j 2 3 4 2 18 2 19 i 2 3 4 j 2 3 4 其中 与表示相应的二阶差商与二阶偏微商的差别为与的数量级 将式 2 18 与 2 19 代入方程 2 13 得 2 20 式 2 20 中去掉项 得到 2 21 i 2 3 4 j 2 3

8、 4 式 2 21 被称为对应于方程 2 13 的差分方程 方程 2 21 被改写成 2 22 若 则式 2 22 又被改写成 或 2 23 物理意义 一点 i j 处的温度是它周围4点温度的平均值 由于差分方程 2 21 是从式 2 20 中去掉项得来的 称去掉的项为差分方程 2 21 的截断误差 当与趋于零时 差分方程的截断误差也趋于零 即差分方程逼近微分方程 我们称这种逼近的差分方程与相应的微分方程为 相容 相容性问题 3 边界条件的差分形式 对流换热边界条件 2 14 用T对x的向前差商代替式 2 14 中的T对x的一阶偏微商 使式 2 14 变成为如下差分形式 这里介绍用差商代替微商

9、的办法把定解问题中的各种边界条件表示成差分的形式 或 2 24 i 1 j 2 3 4 热流边界条件 2 15 用T对y的向前差商代替式 2 15 中T对y的一阶偏微商 使式 2 15 变成为如下差分形式 i 1 2 3 4 5 j 1 2 25 或 绝热边界条件 变成为 i 5 j 2 3 4 2 26 2 16 给定温度边界条件 i 1 2 3 4 5 j 5 2 27 至此 我们对全部节点 包括内节点与边界节点 都用差分形式代替了原来的函数形式 对于内节点上差分形式 我们通称差分方程 因为内节点上温度都是未知的 对于边界节点的差分形式 在边界节点的温度为未知量时 它是差分方程 而对于边界

10、节点为给定的温度时 得到的就不是差分方程了 但在实际应用中 人们往往习惯地把由内节点与边界节点建立起来的差分形式 都统称为差分方程 笼统地讲 一个节点对应一个差分方程 在边界节点的处理方面还有几点需要强调 i 每一边界节点只应属于一种边界条件 如图2 2中 i 1 j 5的节点只属于边界条件式 2 16 iii 若边界节点不正好落在区域的边界上 则需对它们进行特殊处理 ii 对应不同边界的差分方程式 2 24 2 25 2 26 都是用一阶向前差商代替一阶微商得到的 也即它们的截断误差为 或 量级 与内节点差分方程的截断误差相比 低了一个量级 这一点也是从微分形式出发建立差分格式的弱点 4 差

11、分格式的构成 由于式 2 13 2 17 所表示的方程式与边界条件都是线性的 由此而得到的内节点与边界节点的差分方程也都是线性代数方程 由全部节点的差分方程构成一个线性代数方程组 在这个方程组里 方程式的个数等于节点的个数 针对前面讨论的例子 我们可以看到 这个有25个节点所构成的方程组 只需要5个式子 即式 2 21 2 24 2 27 就可表示 这里每个式子都表示几个节点方程 也就是说 在组成方程组时 不必把每个节点方程都写出来 而只要写出几个规格化了的方程就可以了 因此 人们把规格化了的 由内节点与边界节点全部差分方程所构成的线性代数方程组 称之为 差分格式 一般地说 差分格式被写成如下

12、的形式 2 28 其中 n是节点数 也即方程个数 每个方程对应一个节点 和bi i 1 2 n j 1 2 n 都是常数 方程组 2 28 可被进一步写成矩阵的形式 2 29 其中 如果在方程组中去掉其中已知温度节点的那些方程 由此构成的线性代数方程组中 方程的个数等于温度未知的节点个数 也即方程组的未知数 这样的方程组也是差分格式 也可写成式 2 29 综上所述 用有限差分法对式 2 13 2 17 所组成的边值问题的数值处理 最终归结成求解线性代数方程组 2 29 方程组的解即各节点的温度 如果整个区域的节点足够多 那么 离散节点的温度分布就近似代替了区域内的连续温度分布 为便于讨论各种差

13、分格式的优缺点 最好把方程组 2 29 中系数矩阵具体地写出来 但当我们着手书写由式 2 21 2 24 2 27 组成的代数方程组时 发现它所占的版面太大 造成印刷的困难 所以 为便于书写 采用图2 3所示的网格 并假定 得到由全部节点组成的线性代数方程组为 2 30 表示成矩阵形式 2 31 式 2 30 或 2 31 构成差分格式 若在方程组中去掉已知温度节点所对应的方程 即第1至第3个方程 则式 2 31 被改写成 2 32 将式 2 32 与式 2 31 进行对照 即可得到矩阵 的各个元素 这里特别提醒读者注意 在式 2 32 与 2 29 中温度T下角码的对应关系 在讨论二维稳定导

14、热问题时 人们习惯用二个角码 i j 来表示节点位置及对应的曾度 3二维扩散方程的差分解法 考虑平面上粒子扩散的情况 如图3 1所示 图3 1平面扩散 扩散是物质的输运过程 例如气体扩散 半导体扩散 半导体中掺杂 热处理中的渗碳等等 在实际问题中常常遇到 硅片 砷 在窗口有一些粒子向内扩散 例如砷在硅单晶中的扩散 某一t时刻 考虑平面内任意一个闭合回路 根据粒子数守恒 回路内粒子数变化有 其中是杂质 砷 的面密度积分是对回路内面积进行 J是粒子流 将线积分变成面积分 因此 2 33 简化假定 2 34 其中D是扩散系数 简化假定表示引起粒子流的原因是扩散结果 忽略了温度和压力的不均匀引起的粒子

15、流 同时把D看成常数 与面密度无关 是线性扩散问题 因此 2 33 式变为 2 35 这是D为常数时的二维扩散方程 它是一个二阶抛物型偏微分方程 根据初值条件可以对它进行求解 现在我们用差分法对方程 2 35 式进行数值计算 设时间步长为和y步长都是h 将平面分成个网络 网格线的交点称为结点 于是 x ih i 0 1 2 3 N y jh j 0 1 2 3 M 对结点 i j k时刻有 2 36 1 2 36 2 2 36 3 其中第一式用了t的向前差商 第二和第三式是关于x和y的二阶中心差商 即 2 36 式代入 2 35 式 得到方程 2 37 这是一组线性代数方程组 我们已经将偏微分

16、方程化为线性代数方程组 如在k时刻的值分布已算得 可根据方程 2 37 式求得时刻的面密度值 因此当时的值给定后 可依次求得以后各时刻的值 这种差分格式有时也称显式差分法 差分法的收敛条件是 2 38 我们给出此方法的误差是 对和而言为 对而言为 两者合起来误差可写为 其中O代表误差 这是与泰勒级数展开比较得到的截断误差 下面考虑计算流图和程序 二维平面如图3 1所示 M1 M2是窗口 在扩散过程中 窗口处的浓度始终保持一恒定值 我们取的约化值为1 在O M1和M1 M上覆有掩膜使流量为零 在边界处 即流量为零 上下边界处密度 初始时内部全部为零 于是 我们有边界条件 2 39 i 0 j M1 M1 1 M2 窗口 右边界 上下边界 其中k为时间的标号 我们用c表示面密度 用T表示时间步长 用H表示x y两个方向的距离步长h D为扩散系数 假定我们总共计算时间为MT 根据方程 2 37 和初值及边界条件 2 39 式 可画出流图 相应的源程序为 DIMENSIONC1 51 71 C2 51 71 dataN M M1 M2 D H T MT OPEN 6 RILE PRINTER