第五讲——显式差分和隐式差分ppt课件

回顾 有限差分法基础差分格式差分方程边界条件的处理相容性 稳定性和收敛性 回顾 1 有限差分法的相容性 稳定性和收敛性 相容性 针对差分格式而言 在时间步长和空间步长趋近于零的情况下 如果差分格式的截断误差 差分格式与原有偏微分方程之差 的模趋近于零 则该差分格式与原偏微分方程是相容的 或称该差分方程与原偏微分方程具有相容性 稳定性 stability 如果偏微分方程的严格解析解有界 差分格式给出的解也有界 称该差分格式是稳定的 如果差分格式给出的解是无界的 则称该差分格式是不稳定的 稳定性反映了差分格式在计算中控制误差传递的能力 收敛性 convergence 如果当时间和空间步长趋于零时 FDE解趋于PDE解 称该差分格式是收敛的 如果则称该差分格式是收敛的 收敛性描述的是当差分网格无限细化时 差分方程的解是否具有无限逼近偏微分方程的解的能力 Lax等价定理 Laxequivalencetheorem 如果逼近一个给定问题的差分格式是相容的 那么该差分格式的收敛性与稳定性互为充分必要条件 相容性是比较容易满足的 在此基础上 如果满足了稳定性条件 差分格式的收敛性就自动满足 2 5有限差分法实例 forj 2 n 1fori 2 m 1 a j 1 m i j 1 m i 1 1 a j 1 m i j 1 m i 1 1 a j 1 m i j m i 1 a j 1 m i j 2 m i 1 a j 1 m i j 1 m i 4 endend 内部节点 边界节点 A矩阵非零系数减少 同时引入第一类边界 方程右端项B向量出现非零元素 局部节点编号 总体节点编号 组建A和B矩阵 求解线性方程组得到X Matlab2Dclear clc figure color w a zeros 135 135 fori 1 135a i i 1 end fori 1 7a 15 i 1 15 i 2 0 25 a 15 i 1 15 i 16 0 25 a 15 i 1 15 i 14 0 25 endfori 1 7a 15 i 15 15 i 14 0 25 a 15 i 15 15 i 30 0 25 a 15 i 15 15 i 0 25 Enda 1 2 0 25 a 1 16 0 25 a 121 122 0 25 a 121 106 0 25 a 135 134 0 25 a 135 120 0 25 a 15 14 0 25 a 15 30 0 25 fori 2 14a i i 1 0 25 a i i 1 0 25 a i i 15 0 25 endfori 122 134a i i 1 0 25 a i i 1 0 25 a i i 15 0 25 endfori 1 7forj 2 14 a 15 i j 15 i j 1 0 25 a 15 i j 15 i j 1 0 25 a 15 i j 15 i j 15 0 25 a 15 i j 15 i j 15 0 25 endend b a 1 c zeros 135 1 fori 121 135c i 1 25 endd b c s zeros 11 17 fori 2 16s 11 i 100 endfori 1 9forj 1 15 s i 1 j 1 d 15 i 1 j 1 endendsubplot 1 2 1 mesh s axis 0 17 0 11 0 100 subplot 1 2 2 contour s 32 2 5应用实例 南加州一次未来大地震的强地面运动的数值模拟 盆地效应 Cui 2013 Cui 2013 Cui 2013 Cui 2013 总结 1 有限差分方法给出的数值解的精度取决于所用的差分形式 向前 向后 中心 2 偏微分方程的显式有限差分格式通常是有条件稳定的 为了保证得到精确的数值解 最关键的是需要根据稳定性条件选取正确的空间和时间步长 显式与隐式差分格式 主讲人 胡才博中国科学院大学地球科学学院中国科学院计算地球动力学重点实验室 显式差分格式 explicitdifferencescheme 差分方法中可逐层逐点分别求解的格式 特点1 不联立解方程 2 时间步长和空间步长的选择受限制 通常要求时间步长足够小 隐式差分格式 implicitdifferencescheme 特点时间步长和空间步长的选择不受限制 需要联立解方程组 显式和隐式 求解问题与时间相关 例子 1 显式差分格式 左端 n 1时刻的值 右端 n时刻的值 特点 结构简洁 直接求解 求解速度快 但是 时间步长需满足 显式差分格式才能得到稳定的数值解 否则 数值解将会不稳定而振荡 显示差分格式示意图 2 隐式差分格式 时间一阶精度空间二阶精度 隐式有限差分格式 Crank Nicolson隐式差分格式 Crank Nicolson隐式差分格式 Forward TimeCentral Spacemethod Backward TimeCentral Spacemethod Crank Nicolson隐式差分格式 一般差分格式 求解区域 边界条件 初始条件 一种隐式差分格式的程序实现 A sparse nx nx fori 2 nx 1A i i 1 s A i i 1 2 s A i i 1 s end A 1 1 1 A nx nx 1 rhs zeros nx 1 rhs 2 nx 1 Told 2 nx 1 rhs 1 Tleft rhs nx Tright 内部节点 边界节点 载荷项 内部 边界 边界条件 初始条件 Crank Nicolson隐式差分格式的程序实现 A sparse nx nx fori 2 nx 1A i i 1 s A i i 2 2 s A i i 1 s end A 1 1 1 A nx nx 1 内部节点 边界节点 B sparse nx nx fori 2 nx 1B i i 1 s B i i 2 2 s B i i 1 s end B 1 1 1 B nx nx 1 内部节点 边界节点 例子 牛顿冷却定律 温度高于周围环境的物体向周围媒质传递热量逐渐冷却时所遵循的规律 当物体表面与周围存在温度差时 单位时间从单位面积散失的热量与温度差成正比 Tair 一阶常微分方程的数值解 首先对时间和温度进行离散 利用向前差分形式 得到以下的显式差分格式 Tcap 利用向前差分格式 现在改用向后差分形式进行近似 得到隐式差分格式 可以验证 当时间步长趋近于零时 以上近似解趋于解析解 因此 该格式收敛 稳定性条件 T单调减小的条件 显式差分格式 隐式差分格式 当dt 1 25 tau 0 7时 显式差分格式不稳定 结果振荡 隐式差分格式稳定 结果不精确 隐式差分格式 无条件稳定 重点考察差分格式的收敛性 当dt 1 tau 0 7时 显式差分格式不稳定 结果振荡 隐式差分格式稳定 结果不精确 当dt 0 5 tau 0 7时 显式差分格式稳定 隐式差分格式稳定 结果不精确 两者都不精确 当dt 0 1 tau 0 7时 显式差分格式稳定 隐式差分格式稳定 结果都比较精确 当dt 0 01 tau 0 7时 显式差分格式稳定 隐式差分格式稳定 结果都相当精确 当dt和tau都大于零时 该式无条件满足 因此混合差分格式无条件稳定 xt nt 1 nt dt plot xt T e b xt T i g xt T m m xt T a r holdon set gca DataAspectRatio max xt min xt max T e min T e 311 xlabel Time s Fontname timesnewroman FontSize 14 ylabel Temperature Fontname timesnewroman FontSize 14 title dt 0 01tau 0 7 Malab 1Dclear clc figure color w t0 1 initialtemperaturetau 0 7 timeconstantdt 0 01 timeintervalt total 10 nt round t total dt totaltimestepsT e 1 t0 T i 1 t0 T m 1 t0 fori 1 nt xt i i 1 dt T e i 1 T e i 1 dt tau explicitT i i 1 T i i 1 dt tau implicitT m i 1 T m i 1 dt 2 tau 1 dt 2 tau mixT a i t0 exp xt i tau analyticalresultsend 不同差分格式的matlab程序 混合差分格式精度最高 不同差分格式计算结果对比 混合差分格式精度最高 不同差分格式计算结果对比 混合差分格式精度最高 不同差分格式计算结果对比 混合差分格式精度最高 不同差分格式计算结果对比 混合差分格式精度最高 不同差分格式计算结果对比 显式差分格式 1 对步长有要求 2 无需解方程 二阶精度 Matlab2Dclear clc figure color w lx 17 ly 11 v1 zeros ly lx forj 2 lx 1v1 ly j 100 end v2 v1 maxt 1 t 0 k 0 while maxt 1e 6 k k 1maxt 0 fori 2 ly 1forj 2 lx 1 v2 i j v1 i j 1 v1 i 1 j v1 i 1 j v1 i j 1 4 t abs v2 i j v1 i j if t maxt maxt t endendendv1 v2 end subplot 1 2 1 mesh v1 axis 0 17 0 11 0 100 subplot 1 2 2 contour v1 32 迭代解法 K 迭代步数 K 419 总结 显式格式算法简单 易于编程 可以从给定的初始条件开始 在时间上逐层前进求解 一些与时间有关的偏微分方程的求解 需要用到隐式差分格式 在时间上计算数值解的传播时 需要求解线性方程组 通常在计算的每一个时间步 需要求矩阵的逆矩阵 因此 隐式格式算法相对于显式格式更复杂 编程更困难 显式格式通常比隐式格式的稳定性差 如果时间步长取得过大 可能会给出物理上不正确的结果 某些隐式格式的优点是其无条件稳定性 因此时间步长可以取得大一些 但是并不能保证精度很高 可以利用显式 隐式混合格式 如 Crank Nicholsonscheme 它们无条件稳定 而且精度高于相应的显式和隐式格式 总结有限差分法的主要内容包括 1 根据问题的特点将定解区域作网格剖分 在所有网格节点上用有限差分格式对导数求近似 对函数 初始和边界条件求近似 把原方程离散化为代数方程组 即有限差分方程组 2 有限差分模型形态的理论分析 以保证计算过程可行及计算结果正确 解的相容性 对于一个微分方程建立的各种差分格式 为了有实用意义 一个基本要求是它们能够任意逼近微分方程 这就是相容性要求 解的稳定性 因为差分格式的计算过程是逐层推进的 在计算第m 1层的近似值时要用到第m层的近似值 直到与初始值有关 前面各层若有舍入误差 必然影响到后面各层的值 如果误差的影响越来越大 以致差分格式的精确解的面貌完全被掩盖 这种格式是不稳定的 相反如果误差的传播是可以控制的 就认为格式是稳定的 只有在这种情形下 差分格式在实际计算中的近似解才可能任意逼近差分方程的精确解 解的收敛性 一个差分格式是否有用 最终要看差分方程的精确解能否任意逼近微分方程的解 这就是收敛性的概念 总结 3 差分格式 对于有限差分格式 从格式的精度来划分 有一阶格式 二阶格式和高阶格式 考虑时间因子的影响 差分格式还可以分为显式 隐式 显隐交替格式等 目前常见的差分格式 主要是上述几种形式的组合 不同的组合构成不同的差分格式 通过对时间和空间这几种不同差分格式的组合 可以组合成不同的差分计算格式 显式格式算法简单 易于编程 可以从给定的初始条件开始 在时间上 逐层前进求解 隐式格式需要联立求解线性方程组 总结 4 差分格式的构造 构造差分的方法有多种形式 目前主要采用的是泰勒级数展开方法 其基本的差分表达式主要有三种形式 一阶向前