高等数值算法与应用 ( 十二 ) Advanced Numerical Algorithms & Applications 计算机科学与技术系喻文健
Today Summary of Lecture 11 Introduction to Partial Differential Equation Preliminaries of Partial Differential Equation Time-Dependent p Problems Time-Independent Problems Matlab Topics Wenjian Yu 2
Lecture 11 ODE-BVP Boundary Value Problem Separated, linear boundary condition 解的存在 唯一性 线性 BVP 问题解的存在 唯一性 Mode solution: 充分必要条件 : 敏感性 稳定性 Linear BVP 一般地, 不一定成立 非奇异 解的不稳定性受边界条件限制 线性 BVP 问题可推导出条件数, 表明解的敏感程度 Numerical methods Shooting Multiple shooting Finite difference 多次 time integration; 结合非线性方程求解方法 分多个子区间, 增加问题维度, 提高稳定性 差分近似微分 得到离散点处的近似解 Wenjian Yu 3
Lecture 11 ODE-BVP Numerical methods (cont d) 函数空间逼近 确定系数 选基函数 Collocation: 在若干配置点处满足 ODE 和边界条件 最小二乘 : 最小化剩余函数的范数 加权余量法 Galerkin: 基函数的选择 : 谱方法, FEM 相关的特征值问题, 分部积分变换降低求导阶数 global support localized support 将 ODE 离散成代数方程, 转化为代数方程特征值问题 Wenjian Yu 4
Preliminaries of Partial Differential Equation Wenjian Yu 5
两方面的复杂性 我们就一些简单问题讨论基本概念和方法 Wenjian Yu 6
Continuous phenomena modeled d by PDE Maxwell s equations in electro-magnetics 静电场的高斯定律 εe ds = ρdv εe = ρ ( ) S 安培定律 法拉弟电磁感应定律 V H = ε E t + J ( 电流 变化的电场产生磁场 ) ( μh ds) S H E dl = = μ t E l t 磁场是无源场 ( μ H ) = 0 均为 x, y, z, t 的函数 x x Ex ρ Ey = y ε E z z Hx x t E x H y E = μ y y t E z H z z t 其他 : 流体的 Navier-Stokes 弹性力学中 linear elasticity 量子力学中 Schrodinger s 广义相对论中 Einstein s Wenjian Yu 7
类似于 ODE 边值问题 类似于 ODE 初值问题 记号说明 : 求解未知函数 u 在定义域范围内满足给定 PDE 同时满足初始 和边界条件 Wenjian Yu 8
对流方程, 或单向波方程 纯初值问题 验证? Wenjian Yu 9
定义在 (x, t) 二维坐标系统上的曲面 ( 或一系列曲线 ) 等值线, 等高线 x ct = x 0 c =1 Initial function u 0 is propagated to the right (or left) with velocity c Wenjian Yu 10
等值线 刻画依赖关系 特征线 理论上很重要, 确定问题的可解条件 例如, 定义域为 0 x 1, t 0 的常系数对流方程 阴影区域的解由初值 u 0 决定, 还需加边界条件确定其他区域的解 Wenjian Yu 11
( 阶数为最高阶偏导数 ) ( 考虑两个自变量 ) 双曲 判别式 通过变量代换转化为 u u + = 0 tt xx 抛物 u u + = 0 椭圆 u + u + = 0 回忆平面二次曲线方程 : 2 2 ax + bxy + cy + dx + ey + f = 0 通过变量代换 2 b 2 4ac b 2 d b ( x + y ) + y + ( x + y ) + α y + β = 0 2 2a 4a a 2a t xx xx yy Wenjian Yu 12
双曲 特征线的概念对前两种仍适用 对流, 波动 抛物 耗散 扩散 椭圆 更多与物理问题的联系, 见 p. 387 Wenjian Yu 13
Example - Heat equation ( 抛物型 ) 假设不从空气散热, 求时间 t 的温度分布 Δx Ti T i + 1 Q i +1, i T 0 α T T (0, x) = f( x) 热传导定律 0 Qi+ 1, i ΔT Q T = ka = ka 2 Δt Δx t x 能量守恒定律 k T x = cρ 2 Qi+ 1, i Qi, i 1 ΔT Q T cρ AΔ x = 0 ( ) dx= cρ A dx Δt Δt x t t 两端的边界条件 T t L β x T L 给定了初值 以及 Wenjian Yu 14
Example - Laplace equation ( 椭圆型 ) 静电场中的高斯定理 : E εe ds = ρ dv εe = ρ ( ) S V 其中 S 是封闭曲面 由于静电场 E 是保守场, 即沿任意闭合曲线积分为 0, 则可引入标量电位 u, 使得 : E = u 2 2 2 u u u ρ u = + + = 2 2 2 x y z ε 注意 Laplace 算符的不同用法和意义 2 ρ 是空间点 (x,y,z) 处的电荷密度 2D example: Poisson 方程, 若右端为零则称 Laplace 方程 uxx + uyy = 0 给定区域的边界条件后, 可求电势 u 的分布 Wenjian Yu 15
总结 --PDE 基本概念 PDE 问题求解的复杂性 变量多导致问题类型多 : 初值 边值 混合 高维的复杂定义域形状 PDE 问题描述各种自然现象 仅讨论含两独立变量的 PDE 纯边值问题, 不含时间变量 初边值问题, 含时间变量 单向波方程 -- 特征线的概念 如何加边界条件使问题可解 PDE 的阶数 二阶 PDE 分类 ( 双曲 抛物 椭圆 ) 二阶 PDE 举例 -- 热方程 静电场方程 Wenjian Yu 16
Time-Dependent Problems Wenjian Yu 17
抛物型 (heat) 双曲型 (wave) 数值求解方法可分为两类 : 只对空间进行离散的半离散方法 ( 有限差分 collocation 方法 ) 全离散方法 ( 显格式有限差分 隐格式有限差分 ) Wenjian Yu 18
ODE 初值问题 Wenjian Yu 19
因此, 得到齐次线性 ODE = y Ay Wenjian Yu 20
x 轴上的离散点, 对应的时变曲线 u x Wenjian Yu 21 t
关于 Gerischgorin 圆盘定理 : 设矩阵 A = ( a ij ) n n, 定义复数平面的 n 个圆盘区域 D = z a n a i ii ij j = 1, j i 则矩阵 A 的特征值 λi UD i 选择解 ODE 算法时要考虑这种 stiffness Wenjian Yu 22
对给定的时间 t, u(t, x) 变成关于 x 的函数, 再将其表示为一组基函数的线性组合 Wenjian Yu 23
并非通常讨论 ODE 问题时用的显格式 Mα () t = cnα() t Wenjian Yu 24
Wenjian Yu 25
并非直接得到 u 的近似值, 而得到近似函数 而不是平行于时间轴 仍有 stiff 问题 Wenjian Yu 26
并非仅在空间上离散 有限差分近似替代求导 求一系列离散点上的近似解 Wenjian Yu 27
代数方程系统 线性或非线性 从初值函数, 沿时间轴一步步得到后续的解 沿时间轴 推进 的公式可能是显格式或隐格式 Wenjian Yu 28
已知 ( 定解 ) 条件 向前差分中心差分 t k+1 时间点的函数值 Wenjian Yu 29
蜡纸标记 t t k+1 时间点的函数值如何通过前面时间点的值来计算? 通常用 stencil 图直观表示 x 不同于 ODE 中的定义, 但也反映当前计算步的误差 Wenjian Yu 30
Local truncation ti error 定义 : 将准确解带入差分方程, 计算方程的残差 ut ( k+ 1, xi) ut ( k, xi) ut ( k, xi+ 1) 2 ut ( k, xi) + ut ( k, xi 1) c = 2 Δt ( Δx )? ut ( k+ 1, xi) ut ( k, xi) ut ( k, xi) = + O ( Δ t ) Δt t ut ( 2 k, xi+ 1) 2 ut ( k, xi) + ut ( k, xi 1) ut ( k, xi ) 2 = + O(( Δx) ) 2 2 ( Δx) x 2 相减 u u = c 2 t x ut ( k+ 1, x i) ut ( k, x i) ut ( k, x i+ 1) 2 ut ( k, x i) + ut ( k, x i 1) c = O ( Δ t ) + O (( Δ x ) 2 ) Δt ( Δx) 2 Wenjian Yu 31
由于有对 t 的二阶偏导, 增加一个初始条件 局部截断误差 O t O x 2 2 (( Δ ) ) + (( Δ ) ) Wenjian Yu 32
t x 用向前差分近似 Wenjian Yu 33
Semidiscrete 法中, 时间步长选择交给 ODE solver 去做 时间上用向前差分格式 u t = cu xx 完全等价于 MOL 得到的方程用前向 Euler 解 2 Δt λ Wenjian Yu 34
时间上用向后差分格式 完全等价于 MOL 得到的方程用向后 Euler 解 无条件稳定 t x Wenjian Yu 35
Wenjian Yu 36
时间离散 : 梯形公式 k + 1 k ui ui = Δ t k k k 1 ui+ 1 2ui + ui 1 [ c 2 2 ( Δx) k+ 1 k+ 1 k+ 1 ui+ 1 2 ui + ui 1 + c ] 2 ( Δx) 二阶准确度 : 2 阶截断误差 t Wenjian Yu 37 x
例 : 半离散方法, 向后 Euler = y Ay y = y + Δt Ay k+1 k k+1 适度的 ( I - Δt A) y k+1 = y k 稀疏矩阵 二维空间 : 五对角 ( 二阶差分公式中涉及 5 个点 ) 三维空间 : 七对角 ( 二阶差分公式中涉及 7 个点 ) Wenjian Yu 38
收敛性 : 时间 空间步长 0, 则数值解 准确解 一致性, 相容性 : 局部截断误差 0 稳定性 : 同解 ODE- IVP 问题方法的稳定性 充分必要条件 Wenjian Yu 39
分析截断误差, 阶数 p 1 即保证相容 稳定性分析比较困难 y = Ay k k 0 y = ( I + Δt A) y -44 c Eig ( A ) 0 ( Δx) 2 4c Δt 1 - Eig ( I + Δ t A ) 2 ( Δx) 1 依赖区域 显格式 双曲方程 Wenjian Yu 40
依赖区域 特征线 差分格式 显格式双曲型 PDE 的依赖域 显格式! 有限差分格式的依赖域 Wenjian Yu 41
基本解函数 ψ ( x + ct ) ψ ( x ct) 两个特征线 PDE 依赖域 Wenjian Yu 42
Stencil: Slop: Δ t > Δx 1 c Slop: 阴影 : 差分格式的依赖域 CFL 条件 : 差分格式依 Δ t 1 赖域 > PDE 依赖域则稳 < Δxx c 定 Wenjian Yu 43
总结 时变初边值问题的求解 半离散方法 ( 离散空间变量 ) 有限差分法 (Method of Line) Collocation method( 谱方法 有限元 ) 全离散方法 ( 沿时间轴推进 ) 热方程 ( 抛物型 ) 波动方程 ( 双曲型 ) 局部截断误差 稳定性 差分格式的 stencil 图 隐格式 时间后向欧拉 梯形公式 差分方法的收敛性 : 一致性 稳定性 (Lax 等价性定理 ) 稳定性分析 矩阵法 依赖域分析 ( 针对显格式双曲型方程的 CFL 条件 ) Wenjian Yu 44
Time-Independent Problems Wenjian Yu 45
椭圆型 稳态问题 和时间无关 模型问题 第一类第二类 混和边界条件 Wenjian Yu 46
并非按时间一步一步求 通过解代数方程直接得到所有点的近似解, 纯边值 问题 Wenjian Yu 47
为了准确性, 实际步长可能很小 Wenjian Yu 48
通过解代数方程直接得到所有点的近似解 ; 纯边值 问题 并非按时间一步一步求步求 二阶中心差分 稀疏矩阵解法 Wenjian Yu 49
基函数是 B 样条 高维定义域空间的离散 四面体六面体 Wenjian Yu 50
二维 : 双线性 双三次插值 三维 : 三线性 三三次 含未知量增加 矩阵仍然稀疏 Wenjian Yu 51
总结 -- 时不变边值问题的求解 椭圆型方程的种类 -- 亥姆荷兹方程 泊松方程 拉普拉斯方程 三种边界条件 --Dirichlet Neumann mixed 有限差分解法 基本步骤 一个 Laplace 方程的例子 有限元解法 基本步骤 高维情况带来的一些复杂情况 Wenjian Yu 52
Matlab topics Matlab commands for PDE 一维空间的抛物型, 椭圆型方程 pdepe (Solve initial-boundary value problems for parabolic-elliptic PDEs in 1-D) Syntax: sol =pdepe(m,pdefun,icfun,bcfun,xmesh,tspan) m 表示定义域类型 : 长条形 (0), 圆柱对称 (1), 球对称 (2) Pde 方程 : pdefun: [c,f,s] = pdefun(x, t, u, dudx) bcfun: Demo u u, t 0, u(0, x) sinx t xx c 为对角阵, 其余向量 = = 热传导问题 : PDE Toolbox: 2-D FEM 方法 pde_1 Wenjian Yu 53