高等数值算法与应用

Size: px

Start display at page:

Download "高等数值算法与应用"

型圣南门
5 years ago
Views:

1 - 偏微分方程 (chap11) - 喻文健 1

2 基本概念与典型问题有限差分法及泊松方程规则区域泊松方程的快速解法数值稳定性 L 型区域的波方程 2

3 Partial differential equation Three kinds of problems 3

4 基本概念用多元函数的偏微分 ( 偏导数 ) 描述的方程自变量 : 时间 (t), 空间位置 (x, y, z) 二阶线性偏微分方程 ( 含, 也记 u xx ) 各种物理现象例 : 电场中的电势 ; 弹性力学中点的位移含两个自变量的二阶线性 PDE(partial differential equation) 一般形式 : 可根据判别式分类 1. 椭圆型方程, 例 : 泊松方程 u xx +u yy = f(x,y) 2. 抛物型方程, 例 : 热扩散方程 u t u xx =f(x,t) 3. 双曲型方程, 例 : 波动方程 u tt u xx =f(x,t) 这三类方程的定义可推广到含更多自变量的问题椭圆型 ( 与时间无关的稳态 ); 其他两种 ( 与时间相关的过程 ) 4

5 泊松方程 (Poisson equation) ( u) 静电场高斯定理 E ds dv E 电场强度线积分与路线无关, 定义标量电势 u, 拉普拉斯 (Laplace) 算子 : ( ) S 二维泊松方程 : 特例为 Laplace 方程 : u 0 热方程 ( 一维无热源热传导问题 ) 热传导定律 T T Qi 1, i ka Q ka x x V ( 一维 ) ( 二维 ) x T T i i 1 Q i 1, i 能量守恒定律 Q T ( Qi 1, i Qi, i1) t c A x T c A x t T 0 0 E u T ( x,0) u ( x) T 0 t c x 2 T k T 2 TL L x 5

6 热方程 ( 续 ) 一维无热源热传导方程 (u 为温度 ): 热传导方程 ( 考虑热源 ): 波方程 u 为电压 u( x, t) I( x, t) = L x t I( x, t) u( x, t) = C x t 2 2 u( x, t) u( x, t) = LC 2 2 x t 假设系数 LC=1, 得需给初始条件, 例如, 初始条件一维均匀无损传输线 (transmission line)... 假设 Laplace 算子前系数为 1... ( 像波一样, 传输线上的信号会向前传递反射 ) 6

7 偏微分方程的求解多个自变量使问题复杂多个空间自变量使问题的定义域形状复杂对二维空间问题, 考虑 (x,y) 平面内的有界区域, 区域边界上需指定边界条件, 例如 u, 或者 u 的偏导数值已知例 : 拉普拉斯方程 : Δu = 0 一维空间问题, 定义域 axb, 解 u(x) 是连接边界值的线性函数二维空间问题, 定义域 x, y[a, b], 解 u(x, y) 可能很复杂数值求解方法有限差分有限元边界元, 等等求定义域内离散点上的 u 的近似值 7

8 Finite difference method 2-D Poisson equation 8

9 基本思想用差分 ( 差商 ) 代替方程中的偏导数对时间空间定义域进行离散化重点是空间离散, 例如采用均匀网格差分公式拉普拉斯算子 Δu 即二阶中心差分近似二维情况离散 Laplace 算子五点差分格式 2 d u( x) 2 x 也记为 : ( 设单方向有 m 个内部网格点 ) x 据此可联立线性方程组解 Poisson 方程 9

10 差分公式热方程波方程中还有对时间的一阶二阶偏导数一阶偏导可用向前差分 ( 欧拉法 ) 代替, 得到显格式解法热方程问题 ( 设为时间步长 ) t 利用初始条件, 可逐个计算后续时间点的值波方程问题 : 本章考虑的初始条件 h u(x, t) t (x,t ) h u(x, t) x 沿时间轴可逐步算出所有点的函数值! (x,t ) x 10

11 差分算符 Δ h ( 离散的拉普拉斯算子 ) 一维情况下, u(x) 对应的离散值组成的向量 u a u 1 u 2 b 对各节点依次写, 排成一列得算符 Δ h 对应一个三对角矩阵 A u 矩阵 A 稀疏负定 ; 记为 h 2 Δ h 则泊松方程变为 Au = b, 其中向量 b 的值由 f(x) 值构成, 考虑边界值的条件 b: =h 2 f ; b 1 : = b 1 u(a); b m : = b m u(b); 热 / 波方程也可用 Δ h 的矩阵形式 ; 显示公式不需解线性方程组对二维泊松方程, 类似地可得差分算符对应的矩阵应考虑复杂区域的离散 ; 矩阵的形成与点的编号有关 11

12 复杂二维区域上的泊松方程一个小的 L 型区域 h 2 (Δ h u) 1 = u 2 +u 3 +u x 1,y m +u x 0,y m 1 4u 1 h 2 (Δ h u) 8 = u 6 +u 7 +u 9 +u 13 4u 8 算符 h 2 Δ h 对应矩阵 A 对称稀疏矩阵 ( 负定 ), 每行非零元不超过 5 个求解 Au = b 同样, b 包含了 f(x, y) 值及边界值的贡献 y A x

13 用 Matlab 解泊松方程 numgrid 生成几种二维区域的网格编号矩阵 : L=numgrid('L', 7) A=-delsq(L) 针对网格编号 L, 生成 h 2 Δ h ( 五点差分格式 ) 对应的矩阵 A 返回线性编号根据 f(x, y) 及边界值 ( 都乘以 h 2 ) 生成向量 b, 解方程 : u=a\b 另一种生成网格编号 L 的方法 xv = [ ]; yv = [ ]; [x, y] = meshgrid(-1:h:1); [in, on] = inregion(x, y, xv, yv); p = find(in-on); 严格内部点 L = zeros(size(x)); L(p) = 1: length(p); 进行编号区域顶点的坐标值 ( 逆时针顺序 ) 生成网格点 ( 含边界点 ) 的坐标值得到区域上 (in)/ 边界上 (on) 网格点注 : inregion 是 NCM 程序

14 生成网格编号的程序演示与 Δ h 有关的 Matlab 命令一旦有了网格编号矩阵 L, delsq 命令生成 h 2 Δ h 对应的矩阵反过来, 若已知网格上函数值 u, del2 命令计算各网格点上的 Δ h u ( 再乘以 h 2 /4) xv = [ ]; yv = [ ]; [x, y] = meshgrid(-1:1/3:1); y= flipud(y); y 方向向下 [in, on] = inregion(x, y, xv, yv); p = find(in-on); L = zeros(size(x)); L(p) = 1: length(p); 例 : h = 1/20; (Laplace 算子理论值 ) [x,y] = meshgrid(-1:h:1); u = x.^2 + y.^2; d 的维度与 x, y 一样, 值 h 2 4 Δ h, d = (4/h^2) * del2(u); 对此例, 等于 4. 14

15 Matrix eigen-decomposition Fast solver for 2-D region FFT based fast solver 15

16 一维问题中心差分格式 f(x) 回忆习题 10.4, 矩阵 f u h 1 的特征值 f 1 N 记为矩阵 TN k 易知 TN 的特征值为 k 2 2cosk, k, k 1, 2,, N N 1 单位特征向量为 z k, 其元素 2 jk T T zk( j) sin( ) N ZZ N 1 N 1 a u 1 u 2 N b

17 二维问题 T NN = T 2 N N u h f j=4 j=3 j=2 j=1 y h u 1,3 j=0 i=0 i=1 i=2 i=3 i=4 变量按列分成 N 组 T NN T -2I I I N N N N I N I T -2I N N N 矩阵阶数 : N 2 4, 直接解法 ( N ) h x

18 - 将矩阵方程写成另一种形式 2 4u j k u j 1 k u j + 1 k u j k 1 u j k + 1 h f j k 1 jk, N 设 u j,k 组成一个 NN 的矩阵 U ( 每行对应于一列网格点 ) - ( ) 2u j k u j -1 k u j + 1 k T N U j k - ( ) 2u j k u j k -1 u j k + 1 U T N j k k=4 k=3 k=2 k=1 y h u 1,3 k=0 j=0 j=1 j=2 j=3 j=4 h x T U U T h 2 F N F 为 f jk N 组成的矩阵 T N

19 基于特征值分解的解法 T U U T h 2 F N N T T 将 T ZZ 代入, 左乘 Z 右乘 N Z 解出 T T T T T 2 Z ( ZZ ) UZ Z U( ZZ ) Z Z ( h F) Z Z T UZ Z T UZ h 2 Z T FZ 2 T T U U h F, 设 U Z UZ, F Z FZ u jk 2 hf 1 jk, N j jk k 算法 : 1) 2) For all j and k, 3) F U ZFZ ZU Z u jk 计算量 2 hf j jk k (Z 矩阵对称 ) 3 ( N )

20 基于 FFT 的快速解法在上一个算法中, 2 jk Zjk sin, N 1 N 1 1 jk, N 说明矩阵 Z 为 DFT 矩阵的 (2:N+1, 2:N+1) 子矩阵的虚部的倍算法 : 1) F ZFZ 2) For all j and k, 3) U ZU Z u jk 2 (2 2) 2N 2 e i N 考虑 2N+2 阶的旋转因子, 对应 DFT 变换矩阵为 : jk jk N 1 N 1 jk F i, 2N2 2N2 cos sin jk Z( ZF ) F2 N 2 T T 若干矩阵向量乘 2 hf j jk k 计算 Zx 0 jk, 2N 1 2 N 1 用零扩充 x 向量为 2N+2 维 ; FFT(x); 结果的 N 分量虚部计算量 2 ( N log2n ) 2 N 1

21 基于 FFT 的快速解法将计算复杂度由 O(N 4 ) 降为 O(N 2 logn) 适合二三维规则区域结构, 均匀离散网格区域的边界条件也有一定的要求 : 各个面为一种边界条件

22 Numerical stability 22

23 热方程的显式有限差分法 h u(x, t) 一维空间, M 对角元为 1-2, 行元素之和 1 二维空间, M 对角元为 1-4, 行元素之和 1 若 M 对角元 0, M = 1 沿时间轴逐步算出函数值 t (x,t ) 设时刻 t k, 各空间离散点上 u 值组成向量 u (k), x t k +δ 的为 u (k+1), 则得向量方程 : u (k+1) = u (k) + δδ h u (k) = u (k) + δ A σ = δ, 其中 M=I + σa h 2 u(k) = Mu (k) h 2 ( 设 1- -1, 即 2) u (k+1) = Mu (k) u (k+1) M u (k) = u (k) 稳定! 若 M 有对角元 <0, 则 M >1, 误差逐步放大不稳定! A 对一般的二维离散三维问题, 可类似推导! 23

24 热方程的显式有限差分法稳定条件是对角元 1-2, 或一维 : 1/2, 二维 : 1/4 即得到对时间步长的限制 : δ h 2 2, 或 δ h 2 4 波方程的显式有限差分法 t 显格式算法的时间步不能太大! u (k+1) = Mu (k) u (k 1), M = 2I + σa (x,t ) 类似前面分析, M 的对角元为 2-2 ( 一维 ), 或 2-4 ( 二维 ) 若对角元 0, 则 M = 2, 可证明算法稳定稳定性对时间步长的要求 : δ h, 或 δ h 2 δ 2 h 2 x 24

25 有关稳定性的更多讨论为保证显式有限差分的计算稳定, 不能取大的时间步长热方程的稳定性条件较波方程更苛刻 :δ h 2 4 vs. δ h 2 这些稳定性条件也叫 CFL 条件 (by Courant, Friedrichs, Lewy), 他们 1928 年的论文从依赖域的角度阐述了这些条件若采用隐式差分公式 ( 回忆 ODE), 每个时间步内求解线性方程组, 但时间步长可以取得较大 ( 稳定性好 ) 泊松方程 : u = 1 各种二维区域 ([-1, 1] 范围内 ), 三种方程 : 热方程 : u = u+1 t 可选择差分网格的 h, 及 σ = δ h 2 波方程 : 2 u 对热方程与波方程, 观察稳定性 = u t 2 (1/4) (1/2) ( 热方程的解不会稳定住 ) 演示程序 pdegui 25

26 Wave equation for the L region 26

27 零边值波方程的分离变量法 2 u(x, t) t 2 = u(x, t), 考察的情况可验证它满足第 2 个初始条件代入微分方程得 v x 满足的方程为若 u 满足全 0 的边界条件, 则 v 也满足全 0 的边界条件非零解函数 v x 为特征函数, 对应称为特征值 ( 也待定一维情况 v (x)+λv(x)= 0, [a, b] 上的解为 : v k (x)=sin(kπ x a ) b a ) 对应 λ k = kπ b a 2 ; 若定义域为 [0,π], v k (x)=sin(kx) 特征函数的线性组合仍满足 PDE, 根据第 1 个初始条件, 解组合系数, 变为求 u 0 (x) 的正弦级数 (k 为任意整数 ) ( 定义域 [0,π]) 27

28 二维 L 型区域分离变量法 ( 求空域上特征函数 v), 可推广到二维情况 L 型区域虽然简单, 但波动方程无解析解, 且在凹角处有奇异性 ( 梯度无界, 薄膜产生裂缝 ) 分离变量法中的特征函数 v 1 x 成为 Matlab 的 logo L 型区域波问题的实际背景 : 压住 1/4 在风中的毛巾 ; L 形手鼓 ; 脊形微波波导 ( 边缘弯曲效果 ) 波导 - 同轴线适配器电场等位线 28

29 求解 L 型区域的波方程 2 u(x, t) t 2 = u(x, t), 解为 u x,t = a k cos ( λ k t)v k (x) 特征函数仍然满足方程 v x + λv x = 0 可用有限差分法转化为求解矩阵特征值 A h2 v +λv =0 A h2 v = λv m=200; h = 1/m; A=delsq(numgrid('L', 2*m+1))/h^2; [V, D]= eigs(a, 6, 0); 得到特征值与离散的 {v k x }, 实际上可取部分特征值 ( 较小的 ) 与特征函数来近似 u x,t 利用初始条件, u 0 x = a k v k (x), 再用线性最小二乘法可求系数 {a k } k k 29

求解 L 型区域的波方程 m=200; h = 1/m; A=delsq(numgrid('L', 2*m+1))/h^2; [V, D]= eigs(a, 6, 0); 网格点数目达 119201, 常规的 eig 无法解 ( 内存不够 ) 所以用 eigs. 结果为画出第 1 个特征值对应的特征函数 0.02 n= 2*m+1; 0.

30 求解 L 型区域的波方程 m=200; h = 1/m; A=delsq(numgrid('L', 2*m+1))/h^2; [V, D]= eigs(a, 6, 0); 网格点数目达 , 常规的 eig 无法解 ( 内存不够 ) 所以用 eigs. 结果为画出第 1 个特征值对应的特征函数 0.02 n= 2*m+1; Grid= numgrid('l', n); 0.01 pos= find(grid); data(pos)= V(Grid(pos), 6); 0 20 surf(-data, 'LineStyle', 'none'); 40 详见 mymembrane.m 最接近 0 的 6 个, 按从大到小排列 Matlab 的 logo 源于它 30

求解 L 型区域的波方程上述计算采用的是数值方法 : lambda=eigs(a, 6, 0) 只有 3, 4 位准确的有效数字, 误差大是由于奇异点的存在详见课本准确值 membranetx 是 Matlab 程序 membrane 的简化版 ( 可输出特征值 ) membrane 用解析加数值方法计算 L 型区域波方程的特征函数,

31 求解 L 型区域的波方程上述计算采用的是数值方法 : lambda=eigs(a, 6, 0) 只有 3, 4 位准确的有效数字, 误差大是由于奇异点的存在详见课本准确值 membranetx 是 Matlab 程序 membrane 的简化版 ( 可输出特征值 ) membrane 用解析加数值方法计算 L 型区域波方程的特征函数, 使用了分数阶贝塞耳函数做基函数做最小二乘拟合, 比较复杂 ; 但准确度高 ( 含边缘弯曲效果的处理 ) 解析加数值计算更准确 for k=1:6, [L, lambda(k)]=membranetx(k); end; 前 3 个特征函数常作 Matlab 的 logo membrane(1) membrane(2) membrane(3) 31

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),

32 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) sin x t xx PDE Toolbox: 2-D FEM 方法一维空间的抛物型, 椭圆型方程 c 为对角阵, 其余向量热传导问题 : pde_1.m 32

Microsoft PowerPoint - new第十二讲PDE.ppt [兼容模式]

Microsoft PowerPoint - new第十二讲PDE.ppt [兼容模式] 高等数值算法与应用 ( 十二 ) Advanced Numerical Algorithms & Applications 计算机科学与技术系喻文健 Today Summary of Lecture 11 Introduction to Partial Differential Equation Preliminaries of Partial Differential Equation Time-Dependent