110 微分方程的数值解法与程序实现双曲型偏微分方程是对波动和振动现象进行描述的一种微分方程, 常常涉及到更为复杂的情况, 如激波粘性等, 它在海洋气象等许多流体力学的实际问题中都有重要应用. 由于双曲型方程本身描述的是一种波动或振动现象, 需要体现波的传播性质, 而且方程的解对于初值具有局部

Size: px

Start display at page:

Download "110 微分方程的数值解法与程序实现双曲型偏微分方程是对波动和振动现象进行描述的一种微分方程, 常常涉及到更为复杂的情况, 如激波粘性等, 它在海洋气象等许多流体力学的实际问题中都有重要应用. 由于双曲型方程本身描述的是一种波动或振动现象, 需要体现波的传播性质, 而且方程的解对于初值具有局部"

雾元
5 years ago
Views:

1 微分方程的数值解法与程序实现双曲型偏微分方程是对波动和振动现象进行描述的一种微分方程, 常常涉及到更为复杂的情况, 如激波粘性等, 它在海洋气象等许多流体力学的实际问题中都有重要应用. 由于双曲型方程本身描述的是一种波动或振动现象, 需要体现波的传播性质, 而且方程的解对于初值具有局部依赖性, 这些特征都是抛物型方程所没有的. 本章仅对较为简单的双曲型方程进行讨论, 让读者对这一类方程的特性有所了解, 掌握其基本的理论分析方法和数值格式设计. 第一节一阶双曲型方程的若干差分方法本小节考察一阶双曲型方程中最简单的线性模型, 其定解条件可以仅有初始条件, 也可以初边值条件都有. 此处先讨论一阶对流方程, 其半无界的初值问题为 xt (, ) xt (, ) + =, < x<+, t > (4-) ux (,) = ϕ( x), < x<+ 其中, 为非零常数. 在对式 (4-) 进行数值离散求解之前, 需要对连续型的原方程有较深入的了解, 以便分析和建立理论基础并获得后续数值格式的设计. 一精确解所具有的波的传播性质及对初值的局部依赖性可以直接验证 uxt (, ) = ϕ( x t) (4-) du d 是式 (4-) 精确的解析解. 事实上, 若将 x 看作 t 的函数, 则由全微分公式 u x = + 知, 当 dt dt dx dt = du 时就有 dt =, 这说明存在一族特征线 : x = t+ ξ,ξ 为任意常数 (4-3) du 使得在这样的特征线上就有 dt =, 也就是 u 值为常数, 即 uxt ( ( ), t) = ux ( (),) = u( ξ,) = ϕξ ( ). 换言之, 要获得在 x t 平面上的任意一点 ( x, t ) 处的函数值 ux (, t ), 只要把 ( x, t ) 沿特征线投影到 x 轴上得到投影点 ( ξ,), 其中, ξ = x t, 则初始波形 ϕ ( x) 在这一点的值 ϕ( ξ ) 就等于 ux (, t ) ( 见图 4-). 可见, 任意一点 ( x, t ) 处的函数值局部依赖于 ϕ ( x) 在 x 轴上的投影点 ( ξ,) 处的初值. 特征线 x = t+ ξ 的方向代表了波的传播方向, 当 > 时波向右传播, 当 < 时波向左传播, 但波形不变, 波速为. 双曲型方程的解的一个主要特征就是反映了波的传播, 见式 (4-). 对流方程 (4-) 表示的是一个单向传播的波, 因而也称作单向波方程. 以下对一阶双曲型方程半无界初值问题式 (4-) 进行各种差分格式的设计并对它们作进一步的分析. 第一步, 进行网格剖分, 即对 x t 上半平面进行矩形网格剖分, 分别取等距空间步长和时间步长为,, 得到网格节点, Z, =,,, 其中 t =. 第二步, 将原方程弱化, 使之仅在离散的节点处成立, 即

2 第四章双曲型偏微分方程的有限差分法 + =, Z, > t ( x, t) x ( x, t) ux (, t) = ϕ( x), Z (4-4) 第三步, 处理偏导数. 对偏导数用不同的差商近似将建立不同的差分格式. 下面进行具体的讨论. 图 4- > 时的特征线二迎风格式情形 : 关于时间和空间的一阶偏导数都利用一阶向前差商近似, 则有 ux (, t+ ) ux (, t) ux ( +, t) ux (, t), 然后在式 (4-4) 中将精确解 ux (, t ) 用数值解 u 替换并忽略高阶小项, 可得以下数值格式 : 易见此格式的局部截断误差为 O( + ). 若记则式 (4-5) 可以简化为 + u u u+ u + =, Z, u = ϕ( x), Z + + ϕ (4-5) r = (4-6) u = ( + r) u r u, u = ( x ), Z,, (4-7) 数值格式的式 (4-7) 稳定性的讨论不再适宜用直接的矩阵特征值分析方法了, 这主要是因为虽然式 (4-7) 可以写成 u + = Au 的形式, 但此时矩阵 A 原则上是无穷阶的且不是正规矩阵. 所以在这一章我们将统一采用傅里叶分析方法利用增长因子或增长矩阵来讨论稳定性条件, 这本质上是将初始波形函数 ϕ ( x) 看成一系列由各种不同频率的正弦波的叠加, 然后考察在波的传播过程中振幅是否会放大. iω 在实际操作中, 由于方程的线性性, 可以让 ϕ ( x) 是只包含一个波, 即简谐波, 也就是设 ϕ ( x) = e x, 其中 ω 为任意给定的常数, 表示简谐波的频率. 这时在第个时间层上的精确解 uxt (, ) 也相应地只有 iωx 一项记为 uxt (, ) = v( ω)e, 其中 i 为虚数单位, 即满足 i =. 同样, 对离散的数值解设有 u iωx = v ( ω)e (4-8) 式 (4-8) 相当于对离散的 u 进行变量分离. 对数值格式稳定性的考察现在就转化为对振幅 v ( ω ) 是否 + 会放大进行讨论. 如果 v ( ω) = G( ω, ) v ( ω), 则 G( ω, ) 就称为对应于频率为 ω 的增长因子. 这时利

3 微分方程的数值解法与程序实现 + + 用递推关系可得 v ( ω) = G ( ω, ) v ( ω), 即 v ( ω) = G ( ω, ) v ( ω), 从而差分格式的稳定性只要考察 G ( ω, ) 是否关于所有的时间层一致有界, 也就是稳定性可以用增长因子的大小来刻画, 即有下面的定理. 定理差分格式 v ( ω) = G( ω, ) v ( ω) 稳定的充要条件是存在不依赖于频率 ω 的正数 M,,, 使得当空间步长 < 时间步长 < 且 T( 有界 ) 时, 有 G ( ω, ) M, 其中, 为非负整数. + 定理 4.. 差分格式 v ( ω) = G( ω, ) v ( ω) 稳定的充要条件是存在不依赖于频率 ω 的常数 C >, 使得 G( ω, ) + C (4-9) 证明 : 定理 4.. 的证明与定理 3.. 的证明相仿, 故略. 式 (4-9) 称为 Von Neumnn 条件. 现在研究差分格式的式 (4-7) 的增长因子. 易见, 利用式 (4-8) 得 ω ( ) + iωx iωx iωx iωx + i v ( ω)e = ( + r) v ( ω)e r v ( ω)e = ( + r) re v ( ω)e iω 从而 G( ω, ) = + r re = + r( cos ω) irsinω. 注意到 ( ) G( ω, ) = + r( cos ω) + r sin ω = + 4 r( r+ )sin ω 故由定理 4.. 知数值格式稳定等价于 + 4 rr ( + )sin + C, 从而获得式 (4-7) 的稳定性条件 : r <, 即 < 且 (4-) 情形 : 关于时间和空间的一阶偏导数分别利用一阶向前差商和一阶向后差商近似, 则有 ux (, t+ ) ux (, t) ux (, t) ux (, t), 然后在式 (4-4) 中将精确解 ux (, t ) 用数值解 u 替换并忽略高阶小项, 可得以下数值格式 : ω + u u u u + =, Z, u = ϕ( x), Z 易见, 此格式的局部截断误差仍为 O( + ) 且式 (4-) 可以简化为 (4-) + ϕ 仿照情形, 不难得到式 (4-) 的增长因子为此时, ( ) u = ru + ( r) u, u = ( x ), Z, (4-) iω G( ω, ) = r+ r e = r( cos ω) irsinω (, ) ( cos ) sin 4 ( )sin ω G ω = r ω + r ω= r r, 故数值格式稳定等价于 ω 4 r( r)sin + C, 从而获得式 (4-) 的稳定性条件 : < r, 即 > 且 (4-3) 综上, 由稳定性条件式 (4-) 和式 (4-3), 有以下迎风格式 :

4 第四章双曲型偏微分方程的有限差分法 3 < 时, 取数值格式 + u u u+ u + =, Z, u = ϕ( x), Z (4-4) > 时, 取数值格式且稳定性条件为 + u u u u + =, Z, u = ϕ( x), Z (4-5) (4-6) 这样我们可以根据原方程中系数的符号来选取恰当的步长及合适的数值格式. 式 (4-4) 和式 (4-5) 说明迎风格式实际上是在双曲型方程离散的过程中将关于空间的偏导数用在特征方向一侧的单边差商来代替, 体现了原方程中波的传播方向, 它们都是一阶格式. 事实上, 方程 (4-) 含有未知函数关于空间的一阶偏导数项, 也就是对流项, 尽管在数学理论上对这个一阶偏导数进行离散是没有什么特殊困难的, 但在物理过程看却不是这样, 因为对流作用带有强烈的方向性, 所以对流项的离散是否合适直接影响数值格式的性能, 这也就说明了迎风格式之所以有效是因为使用了定向的单边差商. 三一个完全不稳定的差分格式情形 3: 前面讨论了关于时间和空间的一阶偏导数均用一阶差商近似的情况, 接下来容易想到可以对空间的偏导数采用二阶中心差商来近似, 从而有 ux (, t+ ) ux (, t) ux ( +, t) ux (, t), + u u u+ u 然后在式 (4-4) 中将精确解 ux (, t ) 用数值解 u 替换并忽略高阶小项, 可得以下数值格式 : + =, Z, u = ϕ( x), Z 易见此格式的局部截断误差为 O( + ) 且式 (4-7) 可以简化为 (4-7) + r r u = u + u u+, u = ϕ ( x ), Z, (4-8) 利用分解式 (4-8) 可以得到式 (4-8) 的增长因子为 r iω iω G( ω, ) = + (e e ) = irsinω 显然, 对任何 sinω 都有 G( ω, ) = + r sin ω >, 即不存在不依赖于时间步长空间步长的常数 C >, 使得 Von Neumnn 条件式 (4-9) 成立, 从而数值格式的式 (4-8) 完全不稳定. 四蛙跳 (Lepfrog) 格式情形 4: 对情形 3 进行改进, 容易想到对时间和空间的偏导数都采用二阶中心差商来近似, 从而有 ux (, t+ ) ux (, t ) ux ( +, t) ux (, t),

5 4 微分方程的数值解法与程序实现然后在式 (4-4) 中将精确解 ux (, t ) 用数值解 u 替换并忽略高阶小项, 可得以下格式 : + u u u+ u + =, Z, u = ϕ( x), Z 易见, 此格式的局部截断误差为 O( + ) 且式 (4-9) 可以简化为三层蛙跳格式 : + + ϕ u + r( u u ) u =, u = ( x ), Z, (4-9) 注意, 当数值格式不是两层格式时, 前面的利用增长因子的大小来判定稳定性的定理不再适用, 而必须相应地修改. 事实上, 仿照第三章第三节中三层理查德森格式的改写思路, 通过引入中间变量 v = u, 则式 (4-9) 第一式可以写成 : + u = r( u+ u ) + v + v = u u 若记 w =, 就有 v + r r w = w + w + w + (4-) v ( ω) iωx 再设对离散的数值解有分解式 w = e, 那么式 (4-) 就成为 v ( ω) + v ( ω) iωx r i i ( ) ω r ω v ω iωx e = e + + e e + v ( ω) v ( ω) 从而得到增长矩阵为 r iω r iω irsinω G( ω, ) = e + + e = (4-) 于是, 关于差分格式的稳定性在多层情况下可以用增长矩阵来描述. 由于增长因子 ( 两层格式 ) 和增长矩阵 ( 三层及以上格式 ) 在不同的数值格式背景下可以区分, 所以此处都采用同一个记号. 定理 4..3 数值格式 w + = G( ω, ) w 稳定的充分必要条件是存在不依赖于频率 ω 的正数 M,,, 使得当空间步长 < 时间步长 < 且 T ( 有界 ) 时, G ( ω, ) M, 其中, 为非负整数. 定理 4..4 差分格式 w + = G( ω, ) w 稳定的必要条件是存在不依赖于频率 ω 的常数 C >, 使得 ρ G( ω, ) + C (4-) 式 (4-) 也称为 Von Neumnn 条件. ( ) 定理 4..5 当增长矩阵 G( ω, ) 为正规矩阵时, 差分格式 w + = G( ω, ) w 稳定的充要条件是存在不依赖于频率 ω 的常数 ρ G( ω, ) + C C >, 使得 ( ). 以上定理的证明分别类似于定理 3.. 定理 3.. 和定理 3..3 的证明, 故略. 此外, 定理 4..

6 第四章双曲型偏微分方程的有限差分法 5 是定理 4..3 的特例, 相当于定理 4..3 中的增长矩阵为一阶的数, 即增长因子 ; 定理 4.. 也是定理 4..5 的特例. 另外, 再补充几个 Von Neumnn 条件作为差分格式稳定的充要条件的定理 []. 定理 4..6 若记增长矩阵 G( ω, ) = G ( θ), 其中 θ = ω. 对任意给定的 θ R, 下列条件之一成立 : () G ( θ ) 有 p 个不同的特征值, 其中 p 为矩阵 G( ω, ) 的阶数 ; () G ( μ ) ( θ) = γ I, μ =,, ( s, s, 且 G ) ( θ ) 有 p 个不同的特征值 ; μ (3) ρ( G ( θ )) <. 则 Von Neumnn 条件是差分格式稳定的充要条件. 下面利用增长矩阵来分析蛙跳格式的稳定性. 注意到式 (4-) 中增长矩阵 G( ω, ) 的特征值为 λ, = irsinω± r sin ω. 显然, 当 r 时 λ, =, 从而 ρ( G( ω, )) = 也就是 Von Neumnn 条件满足, 但此时由于 G( ω, ) 不是正规矩阵, 所以 r 仅仅是蛙跳格式稳定的必要条件而非充分条件. 注意当 r < 时,Von Neumnn 条件满足且 G( ω, ) = G ( θ) 有两个不同的特征值, 从而由定理 4..6 知, r < π i 时差分格式稳定. 而当 r = 时, 特别地取 θ = ω =, 这时 G( ω, ) =, 且可以直接计算出 n n 3 i 4 5 4i n + i G ( ω, ) =, G ( ω, ) = 并归纳出 G ( ω, ) =, n, 由此 i 4i 3 n n i i = n n+ 知矩阵的行范数 G ( ω, ) = mx gi = + 并不关于时间层一致有界, 其中 g i 为 n G ( ω, ) 的元素, 也就是说 r = 时蛙跳格式不稳定. 综上, 蛙跳格式的稳定性条件为 r <. 蛙跳格式是一个三层格式, 不能自启动进行计算, 需要用一个二阶的格式 ( 如后面将会详细介绍的 Lx-Wendroff 格式 ) 算出第一层信息, 然后再用蛙跳格式进行后面时间层的计算. 五 Lx-Friedrics 格式情形 5: 在情形 3 中修改关于时间的一阶偏导数, 将 ux (, t ) 用其左右相邻两节点的算术平均来 ux (, t ) ( ux (, t) ux (, t)) u 近似, 就是取, 关于空间的一阶偏导数仍用二阶中心 ( x, t) ux ( +, t) ux (, t) 差商, 即, 然后在式 (4-4) 中将精确解 ux (, t ) 用数值解 u 替换并 ( x, ) t 忽略高阶小项, 可得 Lx-Friedrics 格式 : ( + ) + u u + u u+ u + =, Z, (4-3) u = ϕ( x), Z 易见, 其局部截断误差为 O( + ) + O. 式 (4-) 可以改写成 + u = ( + r) u + ( r) u+, u = ϕ( x), Z,

7 6 微分方程的数值解法与程序实现利用分解式 (4-8) 可以得到其的增长因子为 + r iω r iω G( ω, ) = e + e = cosω irsinω, 此时 G( ω, ) = cos ω+ r sin ω = + ( r )sin ω, 故数值格式稳定等价于 + ( r )sin ω + C, 从而获得 Lx-Friedrics 格式的稳定性条件 : r (4-4) 当 r 取定为常数时, 根据式 (4-3) 的局部截断误差知,Lx-Friedrics 格式是一阶格式. 六 Lx-Wendroff 格式下面介绍一个二阶格式, 通过泰勒公式及原方程变形而获得. 由原方程 (4-) 可知 : = 及再根据泰勒公式则有 u u u u u = = = = u 3 ux (, t ) = ux (, t) O( ) + ( x, ) t u 3 = ux (, t) + + O( ) ( x, t ) (4-5) 然后将上式中一阶偏导数和二阶偏导数都用中心差商来近似, 再将精确解 ux (, t ) 用数值解 u 替换并忽略高阶小项, 可得 Lx-Wendroff 格式 : + u+ u u+ u + u u = u +, Z, (4-6) u = ϕ( x), Z 3 其局部截断误差为 O( + + ). 式 (4-6) 可简记为容易得到其增长因子为 + rr ( + ) rr ( ) u = u + ( r ) u + u+, u = ϕ( x), Z, (4-7) rr ( + ) iω rr ( ) iω G( ω, ) e = + ( r ) + e = r ( cos ω) irsinω 4 ω 4ω 显然, G( ω, ) = + 4 r ( r )sin, 故数值格式稳定等价于 + 4 r ( r )sin + C, 从而获得式 (4-7) 的稳定性条件为 r (4-8) 七 Bem-Wrming 格式最后再介绍一个二阶的 Bem-Wrming 格式, 本质上它充分考虑了迎风格式的迎风特点, 同时借用 Lx-Wendroff 格式的设计思想提高了精度.

8 先讨论 < 第四章双曲型偏微分方程的有限差分法 7 的情况. 显然式 (4-5) 仍然成立. 在式 (4-5) 中取 ux ( +, t) ux (, t) u 高阶项, 即取 = + O ( ), 此外, 也取迎风的 ( x, t) 为迎风的形式且兼顾 u ux ( +, t) ux ( +, t) + ux (, t) + O ( ), 把上面两式代入式 (4-5) 得 ux ( +, t) ux (, t) u u 3 ux (, t+ ) = ux (, t) + + O( + ) ( x, t) ( x, t) ( + ) u 3 = ux (, t) ( ux ( +, t) ux (, t) ) + + O( + ) ( x, t) = ux (, t) ( ux ( +, t) ux (, t) ) + ( + ) 3 ( ux ( +, t) ux ( +, t) + ux (, t) ) + O( + + ) 这样, 将精确解 ux (, t ) 用数值解 u 替换并忽略高阶小项, 可得以下 < 时的 Bem-Wrming 格式 : + r( + r) u = u r( u+ u) + ( u+ u+ + u), Z, u = ϕ( x), Z 3 = (4-9) 差分格式 (4-9) 的局部截断误差均为 O( + + ). 再利用分解式 (4-8), 代入式 (4-9) 可得增长因子为 iω rr ( + ) i ω iω G( ω, ) = r(e ) + (e e + ) 经过整理化简可得 G( ω, ) = r( + r) ( + r)(cosω ) +, 故数值格式稳定等价于 + r( + r) ( + r)(cosω ) + C, 从而得稳定性条件为 r, 再联合一开始的假设 <, 最终可得式 (4-9) 的稳定性条件为 r < (4-3) 用同样的思路, 可以得到 > 时的 Bem-Wrming 格式 : + r( r) u = u r( u u ) ( u u + u ), Z,, u = ϕ( x), Z (4-3) iω r( r) iω i ω 其增长因子为 G( ω, ) = r( e ) ( e + e ), 经过整理化简可得 G( ω, ) = r( r) ( r )(cosω ) +, 类似可得使数值格式 (4-3) 稳定的条件为 < r (4-3) 八隐格式的设计前面介绍的若干数值格式都是显格式, 从而必须附带稳定性条件成立才能保证数值解最终收敛到

9 8 微分方程的数值解法与程序实现精确解. 而事实上隐格式通常稳定性较好, 所以可以考虑设计隐格式来求解. 比如 > 时, 修改情形 3 那个完全不稳定的显格式为隐格式, 即关于时间和空间的一阶偏导数分别利用一阶向后差商和二阶中心差商近似, 就有 ux (, t) ux (, t ) ux ( +, t) ux (, t), 然后在式 (4-4) 中将精确解 ux (, t ) 用数值解 u 替换并忽略高阶小项, 可得以下数值格式 : u u u+ u + =, Z, u = ϕ( x), Z (4-33) 此格式的局部截断误差仍为 O( + ) 且式 (4-33) 可以简化为 r r u + u + u+ = u, u = ϕ( x), Z, 或者 r + + r + u + u + u+ = u, u = ϕ( x), Z, (4-34) 利用分解式 (4-8), 可得式 (4-34) 的增长因子为 G( ω, ) = = + (e e ) + r iω iω irsin 由于对任意 r 均有 G( ω, ) =,Von Neumnn 条件满足, 由定理 4.. 知差分格式的 + r sin ω 式 (4-33) 无条件稳定. 除了上面可以修改情形 3 为隐格式, 其他的迎风显格式等都可以这样操作, 这里就不一一列举了, 读者可以自己进行格式的设计和分析. 九 Cournt-Friedrics-Lewy 条件在进行收敛性分析之前, 我们还需要针对双曲型方程波的传播特征, 研究差分格式离散解的依赖区域. 从而获得差分格式收敛的一个必要条件 Cournt-Friedrics-Lewy 条件, 简称 CFL 条件, 或者 Cournt 条件. CFL 条件本质上说的是原方程解的依赖区间必须包含于差分格式解的依赖区间. 为此, 我们举例说明. 首先以在 > 时式 (4-5) 的迎风格式为例. 前面已经分析过, 原微分方程 (4-) 的精确解 uxt (, ) 在网格节点 Px (, t ) 处的值只依赖于过这一点的特征线在 x 轴上的投影点 P 处的初值, 从而 P 是 uxt (, ) 在点 Px (, t ) 处的依赖区域 ( 一个点 ). 而差分格式的式 (4-5) 的解在点 Px (, t ) 处的值 u 依赖于前一时间层上的 u, u, 而 u, u 又分别依赖于更前一层上的 u, u 和 u, u, 依次类推, 可得 u 依赖于初始时间层上的 u, u,, u, 这样, 数值解 u 的依赖区域是区间 [ x, x ], 或者说数值格式的式 (4-5) 的解 u 的依赖区域为 [ x, x ]. 再考察数值解与精确解依赖区域的关系可知, 在点 Px (, t ) 处, 如果精确解的依赖区域 P 在区间 [ x, x ] 之外, 那么用数值格式的式 (4-5) 算出来的解 u 就与原问题式 (4-) 的解毫无关系, 因此这个数值格式的解就不可能收敛到原方程的解. 这样说来, 数值解收敛到精确解的一个必要条件就是 P [ x, x], 这也就意味着 x x t x, 即 > 且, 也就是式 (4-5) ω

10 第四章双曲型偏微分方程的有限差分法 9 的稳定性条件式 (4-3). 换言之, 对迎风格式的式 (4-5) 而言, 其 CFL 条件与稳定性条件一致. CFL 条件是差分格式收敛的一个必要条件而不是充分条件. 这个可以从情形 3 那个完全不稳定的格式的式 (4-7) 得到说明. 对于式 (4-7) 而言, 易见其数值解在网格节点 Px (, t ) 处的值依赖于前一时间层上的 u, u, u+, 而 u, u, u+ 又分别依赖于更前一层上的 u, u, u 和 u, u, u+ 及 u, u+, u+, 即 u, u, u+ 依赖于 u, u, u, u+, u+, 依次类推, 可得 u 依赖于初始时间层上的 u, u +,, u, u, u+,, u+, u+, 于是数值解 u 的依赖区域是区间 [ x, x+ ]. 这样格式的式 (4-7) 的数值解收敛的必要 CFL 条件就是 P [ x, x+ ], 即 x x t x+, 即, 也就是 r. 而前文已经分析知式 (4-7) 完全不稳定, 即对任意 r 数值格式不稳定从而不收敛, 这就说明 CFL 条件只是差分格式收敛的必要条件而非充要条件. 最后再说明一下对于不同的数值格式,CFL 条件一般是不同的, 如上面分析的式 (4-5) 和式 (4-7) 的 CFL 条件是不同的. 十数值算例例 4.. 求解一阶对流方程初值问题 + =, < x<+, t > t x ux (,) = ϕ( x), < x<+, x 其中, 初值 ϕ( x) = 在 x = 处间断. 取空间步长和时间步长分别为 =., =.5, 给出, x > 时刻 t =.5 时 x [,] 区间内数值解的图像., 解 : 易见上述初值问题的精确解为 uxt (, ) = ϕ( x t) =, 易得 r =.5. 精确解见图 4-. x t 且对应上述空间时间步长的选取 x > t 图 4- 例 4.. 的精确解图像此处分别采用迎风格式 Lx-Friedrics 格式 Lx-Wendroff 格式 Bem-Wrming 格式进行计算, 对应程序分别为 Egc4_sec_.c Egc4_sec_.c Egc4_sec_3.c Egc4_sec_4.c. 为了保证在解的依赖区间的影响下能获得 x [,] 区间内的数值解, 在迎风格式 Lx-Friedrics 格式和

11 微分方程的数值解法与程序实现 Lx-Wendroff 格式中取子区域 x [,] 进行计算, 在 Bem-Wrming 格式中取子区域 x [,] 进行计算. 为了借助 MATLAB 画出直观图, 我们在程序设计时采用了文件指针, 以方便数据的存储与读取. 各种算法的数值解如图 4-3 所示. () 迎风格式 (b) Lx-Friedrics 格式 (c) Lx-Wendroff 格式 (d) Bem-Wrming 格式图 4-3 数值解的图像十一推广对于空间区域有界的一阶方程的初值问题, 要根据的正负号给出相应的边界条件才能使问题适定, 即唯一可解. 具体来说, 当 > 时, 波向右无限传播, 所以需要给出左边界条件 ; 当 < 时波向左无限传播, 需要给出右边界条件. 例如, 考察以下初边值问题 : xt (, ) xt (, ) + =, x x xb, t >, < 为常数 (4-35) ux (,) = ϕ( x), x x x b ux ( b, t) = ψ ( t), t> 其中, x, x b为已知常数, ϕ( x), ψ ( t) 为已知函数. 显然与前面半无界初值问题式 (4-) 一样, 式 (4-35) 也有一族特征线, 其方程为 x = t+ ξ,ξ 为任意常数, 且沿着特征线函数 uxt (, ) 为常值. 在这些特征线中有一条临界线 L, 如图 4-4 所示, 这条特征线将求解区域分成两部分 Ω 和 Ω. 由于式 (4-35) 表示的波函数是向左传播的, 所以在 Ω 内的点 ( x, t ) 处的函数 u 的值依赖于右边

12 第四章双曲型偏微分方程的有限差分法界条件, 而在 Ω 内的点 ( x, t ) 处的函数 u 的值则依赖于初始条件. 为了得到整个求解区域内任意一点 ( x, t ) 处的精确解 uxt (, ), 分三步走. 图 4-4 精确解值与初边值条件的依赖关系 ( 向左传播 ) 首先, 在 x t 平面内容易得出临界线 L 的方程为 : t = ( x x b ) 其次, 考察 Ω 内的任意一点 A( x, t ) 处的函数值 ux (, t ). 注意到过点 A 的特征线方程为 t = ( x x ) + t, 且易知沿着这条特征线, 点 A 在右边界上的投影点为 A ( xb, α ), 也就是这条特征线与右边界的交点, 易知 α = ( x ) b x + t. 于是, A 点处的 u 值就是 A 点处的 u 值, 即有 u u ( ) ( x b x ) t = = ψα = ψ A A +. 再由 A 点的任意性知, Ω 内的任意一点 ( x, t ) 处的 u 值为 uxt (, ) = ψ ( xb x) + t, ( xt, ) Ω 最后, 再考察 Ω 内的任意一点 Bx (, t ) 处的函数值 ux (, t ). 注意到过点 B 的特征线方程为 t = ( x x ) + t, 且易知沿着这条特征线, 点 B 在 x 轴上的投影点为 B ( β,), 也就是这条特征线与 x 轴的交点, 易知 β = x t. 于是, B 点处的 u 值就是 B 点处的 u 值, 即有 u = u = ϕβ ( ) = ϕ( x t ). 再由 B 点的任意性知, Ω 内的任意一点 ( x, t ) 处的 u 值为 B B uxt (, ) = ϕ( x t), ( xt, ) Ω 综上, 一阶双曲型方程初边值问题式 (4-35) 的精确解为 ψ ( xb x) + t, t > ( x xb) uxt (, ) = ϕ( x t), < t ( x xb ) 类似地, 考察以下初边值问题 : xt (, ) xt (, ) + =, x x xb, t >, > 为常数 ux (,) = ϕ( x), x x x b ux (, t) = ψ ( t), t>

13 微分方程的数值解法与程序实现仿前面分析可知 ( 课后习题 ), 上述问题的精确解为 ψ ( x x) + t, t > ( x x) uxt (, ) = ϕ( x t), < t ( x x ) 关于上述问题的数值格式的设计可以由本节前面介绍的半无界初值问题的若干数值方法加以推广, 此处不再展开, 具体问题可通过课后习题得以深入学习. 第二节二阶双曲型方程的显式差分法本节主要研究下面的二阶双曲型方程的初边值问题 : uxt (, ) uxt (, ) = f( x, t), < x<, < t T (4-36) ux (,) = ϕ( x), ( x,) = ψ( x), x u(,) t = α(), t u(,) t = β(), t < t T 其中,u 表示一个与时间 t 和位置 x 有关的待求波函数, ϕ( x), ψ( x), α( t), β ( t) 及方程右端项函数 f ( xt, ) 都是已知函数,, T 是非零常数. 式 (4-36) 中第二式是初始条件, 第三式是边界条件. 式 (4-36) 中若 f( x, t), 则对应的初边值问题与一阶对流方程有类似的属性, 只是它有两条特征线, 精确解 uxt (, ) 在点 ( x, t ) 处的值的依赖区域为 [ x t, x+ t], 从而也有相应的 CFL 条件, 关于这些理论分析将不再展开. 本节重点处理较为一般的非齐次的二阶双曲型方程初边值问题式 (4-36), 通常无法用解析的办法求得其精确解, 所以我们需要进行数值计算求得原方程的近似解, 也就是所谓的数值解. 一三层显差分格式的建立有了以往的经验, 按照正常的步骤建立差分格式. 第一步, 区域剖分. 由于原问题的求解区域是矩形域, 所以我们将在矩形域网格点上近似求解. 首先, 对原始问题的求解区域 x, t T 进行网格剖分, 为计算简单, 确切地说是为了编程方便, 通常取为等距剖分, 即取 x = ( =,,, m), t = ( i =,,, n), 这里 = / m, = T / n 分别称为空间步长和时间步长. 第二步, 对原方程进行弱化, 将原来很强的处处成立的微分方程弱化为仅在某些点成立的方程, 这里设方程仅在网格节点成立, 即得节点处离散方程 : u u = f( x,,,, t ) < < m < n ( x, t) ( x, t) ux (, t) = ϕ( x), ( x, t) = ψ( x), m, ux (, t) = ( t), ux (, t) β( t), n. α m = (4-37) 第三步, 将各节点处的偏导数用差商来近似. 为此容易想到方程中的二阶偏导数都取成二阶中心差商近似, 即

14 第四章双曲型偏微分方程的有限差分法 3 ux (, t ) ux (, t) + ux (, t ), u + t ux (, t) ux (, t) + ux (, t ), u + x 误差为 O( ), 误差为 O ( ). 此外, 初始条件中的一阶偏导数可以用向前差商近似, 即 ux (, t) ux (, t), 误差为 O( ). 将上面各式代入式 (4-37) 第一式, 然后用各节点 ( x, t ) 处的数值解 u 代替精确解 ux (, t ), 并忽略高阶小项, 得以下差分格式 : + u u + u u+ u + u = f( x, ),,, t m n u u u = ϕ( x), = ψ( x), m, u = α( t), um = β( t), n. 为简单起见, 记网比 r =, 整理上述格式后得到以下三层显格式 : + u = ru + ( r) u + ru+ u + f( x, t), m, n, u = ϕ( x), m; u = u + ψ( x), m, u = α( t), um = β( t), n. 易见, 此格式的局部截断误差为 O( + ). 二显格式的稳定性收敛性分析同样, 只需要对齐次方程零边界情形讨论. 此时, 式 (4-38) 第一式为 + = ( ) u ru r u ru u (4-38) 同三层的蛙跳格式一样, 通过引入变量 v u =, 可将三层格式写成两层格式 : u = ru + r u + ru v + v = u + ( ) + + u r u ( r) u r u + 也就是, = v v v v + 记 w ( ω) = = e ( ω) u v v v iωx, 则由上式得 + v ( ω) iωx r i ( ) i ( ) ω r r ω v ω e = e + + e e + v ( ω) v ( ω) iωx

15 4 微分方程的数值解法与程序实现即 w + = G( ω, ) w, 从而得到增长矩阵为其中, r = >. 易见, 增长矩阵的特征值为 ω 4r sin G( ω, ) = ω ω ω λ, = rsin ± 4rsin r sin (4-39) ω ω ω 显然, r 时, λ, = rsin ± i 4rsin r sin 且 λ,, 也就是显差分格式 (4-38) 满足 Von Neumnn 条件, 但这只是式 (4-38) 稳定的必要条件. 当 r < 且 ω pπ (p 为整数 ) 时, 则 λ, λ 是增长矩阵的两个不同的特征值, 从而由定理 4..6 知, 此时差分格式稳定. 但如果 r < 且 + ω = pπ (p 为整数 ), 则 G( ω, ) =, 通过直接计算可知, G ( ω, ) =, ( ) 从而 G ( ω, ) ( ), 也就是 G ( ω, ) 关于不是一致有界的, 从而格式不稳定. 以上说明, 当 r < 时, 并不是对所有的频率 ω 显格式都稳定, 也就是 r < 时三层显格式不稳定. 此外, 当 r = 时, 若取 ω = π, 则 G( ω, ) =, 也会出现 G ( ω, ) ( ), 从而数值格式不稳定. 至于 r > 的情况, 取 ω = π, 则 λ = r 4 r( r ) <, 显然不满足 Von Neumnn 条件, 从而数值格式不稳定. 综上, 将得出显格式绝对不稳定的结论. 但实际上, 这个结论是错误的. 上述稳定性分析也是不正确的. 原因在于本章第一节给出的判断稳定性的充分条件都是针对一阶方程而言的, 而这次研究的对象却是二阶双曲型方程, 关于二阶方程差分格式的稳定性条件会更宽松一些. 事实上, 我们需要一个弱一点的稳定性概念阶稳定. 定义 4.. 增长矩阵 G( ω, ) 关于范数是阶稳定的, 如果存在非负常数 M, M, 对任意频率 ω 都有 G ( ω, ) M+ M, 其中代表时间层, n. 定理 4.. [] 若增长矩阵 G( ω, ) 关于范数是阶稳定的, 且第个时间层节点处误差均为 O( + ), 则差分格式 w + = G( ω, ) w 收敛. 证明 : 由于方程 (4-36) 中的二阶偏导数均用中心差商来近似, 就有 ux (, t+ ) ux (, t) + ux (, t ) ux ( +, t) u( x, t) + u( x, t ) = f( x, ) ( ) t + O + 也就是, ux (, t ) = rux (, t) + ux (, t) + ( rux ) (, t) ux (, t ) + f( x, t) + R ( ) + + 其中, R = O( + ). 将式 (4-38) 第一式减去上式可得误差 e = u u 满足方程 + ( ) + e = re + r e + re e R. 引入中间变量 ε = e, 则上式可改写为 e = re + r e + re R + ε = e + ( ) + ε

16 第四章双曲型偏微分方程的有限差分法 5 从而可得 + e r e ( r) e r e + R = ε ε + ε ε + e v ( ω) iωx 再令 e w = =, ε v ( ω) R iωx = ξ e 就得 + v ( ω) v ( ω) ξ = G( ω, ) + + v ( ω) v( ω) 不妨记为 W + = G( ω, ) W + R (4-4) T T 其中, W = ( v ( ω), v ( ω)), R = ( ξ,). 易见, R = O( + ). 利用式 (4-4) 进行递推可得 + l l W = G ( ω, ) W + G R. 若增长矩阵 G( ω, ) 关于范数是阶稳定的, 则 l= + G ω + M+ l M l l= W (, ) W ( ) R G ( ω, ) W + ( M + n M ) R l= l * + W ( M M ) ( M nm ) nc O( ) 其中, C * l 为各个 R 中的项 O( + ) 中关于精确解 uxt (, ) 高阶导数的最大值. 注意到 ε = e = u u =, 所以 W 为 O( + ) 时, 即 e = O( + ) 时, T T = ( v( ω), v( ω)) = ( v( ω),) W v ω = e = O + ( ) ( ), 且当第个时间层节点处误差均故 + * M+ M + M+ nm nc O + W [( ) ( ) ] ( ) (4-4) 可见当, 时, 由于 n = T 有界, 上式右端方括号内的表达式有界, 从而就有 e +, 说明 u u, 差分格式二阶收敛. 证毕. + W, 也定理 4.. 说明在稳定性条件从增长矩阵的范数一致有界 ( G ( ω, ) M ) 放宽到增长矩阵次幂的范数关于线性增长 ( G ( ω, ) M+ M) 的情况下, 如果初始计算 ( 指的是第个时间层各节点处 ) 的误差较小 ( 第个时间层各节点的值是准确的 ), 小到 O( + ), 则差分格式是二阶收敛的. 所以, 对于二阶双曲型方程初边值问题三层显格式 (4-38) 的稳定性条件可以重新进行讨论. 下面就来重新考察增长矩阵式 (4-39) 的范数. 事实上, 前面已经分析出当 r < 且 ω pπ (p 为整数 ) 时, 差分格式稳定. 现在, 如果 r < 且 ω = pπ, 则 G( ω, ) = 且 + G ( ω, ) = ( ) 从而 G ( ω, ) G ( ω, ) = +, 其中, 是矩阵的行范数, 故 G( ω, ) 关于范数是阶稳

17 6 微分方程的数值解法与程序实现 ω 4sin μ 定的. 再考察 r = 的情况, G( ω, ) = : =, 易见 μ. μ = 的情况刚 ( ) ( + ) ( ) 才已经讨论过. μ = 时, G( ω, ) = 且类似地有 G ( ω, ) =, + + ( ) ( ) ( ) 所以 G ( ω, ) G ( ω, ) = +, 故 G( ω, ) 关于范数是阶稳定的. 现在只需要讨论 μ < 的情形. 此时, G( ω, ) 有一对共轭复根, 从而由定理 4..6 知数值格式稳定. 相比较之前得到的错误结论显格式绝对不稳定, 现在知道错误的真正原因在于我们之前给出的稳定性判断条件太强了. 最后, 分析一下第个时间层节点处的误差. 注意到第个时间层节点处我们采用的是式 (4-38) 中 u = u +ψ ( x ), 而利用泰勒公式精确值应为 u u = + u u + ( x, tξ ), 其中, tξ ( t, t ). 注意到原方程的初始条件 u( x, t ) = ϕ ( x) = u 和 u ( x, t ) ( x = ψ ), 故 e = u u( x, t ) = O( ), 所以此时的式 (4-4) 只能达到所以当, 时, 由于 n = T + * W M+ M + M+ nm nc O + [( ) ( ) ] ( ) + 有界, 上式方括号内表达式有界, 从而 W, 也就有 e +, 说明 u u, 差分格式收敛, 但关于时间只有一阶收敛, 关于空间是二阶收敛的. 综上, 二阶双曲型方程初边值问题式 (4-36) 的三层显差分格式的式 (4-38) 的稳定性条件是 r, 且关于时间一阶收敛, 关于空间二阶收敛. [,3,4] 另外说明一下, 有不少教材将二阶双曲型方程降阶为一阶双曲型方程组, 然后对一阶双曲型方程组进行差分格式设计, 其设计思路不清晰, 得到的格式也不常规, 虽然最终的数值格式与我们直接通过二阶差商近似得到的显格式等价, 但编者认为, 这些教材这样做的目的只是为了说明在上述操作过程下可以得到正确的显格式的稳定性条件, 而这些其实都不是本质的, 真正本质的是二阶方程的稳定性条件比一阶方程的弱. 三改进的三层显格式数值格式 (4-38) 精度比较低的原因是在对初始条件进行一阶偏导数的近似时, 用了精度较低的一阶向前差商近似, 导致 e 只能达到 O( ). 为了提高精度, 可以考虑关于时间的一阶偏导数用中心差商近似, 即取 ux (, t ) ux (, t ), 误差为 O( ) 从而得数值格式 : + u = ru + ( r) u + ru+ u + f( x, t), m, n, u = ϕ( x), m; u = u + ψ( x), m, (4-4) u = α( t), um = β( t), n. 但这样做又会引入新的虚拟越界点 u. 有了以往的经验, 就应该知道下面的操作了, 即认为式 (4-4)

18 中第一式对 = 也成立, 即第四章双曲型偏微分方程的有限差分法 7 + u = ru + ( r) u + ru u + f, m 从而解出 u = ru + ( r) u + ru+ u + f( x, t), m, 这样, 在式 (4-4) 的第二式中可以消去 u, 因此可以得到 + ) u = ( ru + ( r) u + ru + f( x, t + ψ( x ))/, m (4-43) 从而得到以下改进的三层显格式 : + u = ru + ( r) u + ru+ u + f( x, t), m, n, u = ϕ( x ), i m, u = ( ru + ( r) u + ru+ + f( x, t) + ψ( x) )/, m, u = α( t), um = β( t), n. (4-44) 显然, 此格式的局部截断误差为 O( + ). 此外, 注意到由式 (4-43) 泰勒公式和原方程 u u = f( x, t), 可知此时 e = u u( x, t ) = O( ), 即式 (4-4) 成立. 这样在稳定性条件 r 成立的前提下, 式 (4-44) 关于时间和空间都是二阶收敛的. 四数值算例例 4.. 用一阶显格式 (4-38) 计算双曲型方程初边值问题 : uxt (, ) uxt (, ) t = e sin x, < x< π, < t, ux (,) = sin x, ( x,) = sin x, x π, u(, t) =, u(π, t) =, < t. t 已知其精确解为 uxt (, ) = e sinx. 分别取步长为 π =, = 和 =, π =, 给出在节点 5 iπ 4,, i =,,9 处的数值解及误差. 5 解 : 程序见 Egc4_sec_.c, 计算结果列表如下 ( 表 4-). 表 4- 一阶显格式的计算结果 ( x, t ), 步长下 u 误差 u (, ) u x t, 步长下 u 误差 u (, ) u x t (π/,.8) e e 3 (π/,.8) e e 3 (3π/,.8) e e 3 (4π/,.8) e e 3 (5π/,.8) e e 3 (6π/,.8) e e 3 (7π/,.8) e e 3 (8π/,.8) e e 3 (9π/,.8) e e 3

没有幻灯片标题

没有幻灯片标题第九章常微分方程数值解法 Euler 方法 Ruge-Kutta 法 3 单步法的绝对稳定性 4 线性多步法 5 一阶方程组与高阶方程的初值问题 -- 常微分方程数值解法必要性在工程和科学技术的实际问题中, 常需要求解微分方程只有简单的和典型的微分方程可以求出解析解, 而在实际问题中的微分方程往往无法求出解析解 y xy 如微分方程初值问题 y(0 0, 其解析解 ( 精确解为 : x t y(