微分方程数值解第五章答案.pdf

第五章 双曲型方程的数值解 1 第五章 第五章 1 0 0 0 1 2 0 0 0 x uu u xx tx x 1 1 对初值问题对初值问题 2 试分别用左偏心格式 试分别用左偏心格式 LW 格式计算其数值解格式计算其数值解u k 1 2 3 4 取取 1 h k 解 矩形网格剖分区域 取空间步长h 时间步长 的矩形网格剖分 区 域 用 节 点表 示 坐 标 点0 1 2 j j k jk x tjh k 0 1 2 3 4 k 0 k j k j x u t u 1 左偏心格式 在 t 上用向前差商 x 上用向后 差商 得0 1 1 h uuuu k j k j k j k j 中国地质大学 北京 廉海荣编 1 因为2 1 h 整理得到 k j k j k j uuu 2 1 2 1 1 1 把已知条件离散成 则可以根据下一层求上一层的值得 到 1 2 3 4 下图中节点处值即为求出来的值 0 0 0 2 1 0 1 j j j 0 j u k u k u 第五章 双曲型方程的数值解 3 0 0 0 x u x 0 x 0 x u 0 t t 0 uu atT tx tT 11 11 0 2 kkkk jjjj uuuu a h 1 a 111 1111 222 kkkkkk jjjjjj uuuuuu a hh b 0 试分析它们的稳定性 试分析它们的稳定性 解 首先用矩形网格剖分区域 取空间步长 时间步长h 的矩形网格 剖 分 区 域 用 节 点 表 示 坐 标 点0 1 2 j j k jk x tjh k 0 1 2 T k a 对于初值问题的方程在处取值 然后用向后差分代替时间偏 j k 第五章 双曲型方程的数值解 4 导 然后用中心差商代替空间偏导为 舍去误差项 得到差分格式 1 11 0 2h kkkk jjjj uuuu a 然后用 k 1 层代替 k 层有差分格式 11 11 0 2h kkkk jjjj uuuu a 1 5 1 下面我们对这个差分格式的稳定性进行分析 格式等价于 111 11 1 2 kkkk jjj uar uu j u h 令代入上式 得 kkij j uv e 111 1 2 kijhkijhi hkijhi hkijh vear veeveev e 进一步简化得传播因子为 11 1 1sin 1 2 i hi h G iarh ar ee 显然 222 111 1sin1sin 1sin G iarhiarh a rh 对于 任意的 r 可知 G 0a 恒小于 1 VonNeumann 条件成立 故差 分格式无条件稳定 b 对于初值问题的方程在 处取值 然后再用向前差分代替时间 偏导 接着用中心差商代替空间偏导数 舍去误差项定义差分格式 j k 1 11 0 2h kkkk jjjj uuuu a 5 2 1 2 然后 5 1 和 5 2 式两边均乘以 做和有 中国地质大学 北京 廉海荣编 4 第五章 双曲型方程的数值解 5 111 1111 0 22h2h kkkkkk jjjjjj uuuuuu a 下面对这个差分格式的稳定性进行分析 差分格式等价于 111 1111 11 44 kkkkk jjjjj uar uuuar uu k j h 令代入上式 得 kkij j uv e 111 11 44 kijhkijhi hkijhi hkijhkijhi hkijhi h vear veeveev ear v eev ee 继续化简有 1 11 1sin1sin 22 kk iarh viarh v 1 1sin 2 1 1sin 2 iarh G iarh 从而得到传播因子为 显然 222 222 1 11 1sin 1sin1sin 2 42 1 11 1 1sin1sin 1sin 22 4 iarh a rhiarh G iarhiarh a rh VonNeumann 条件成立 故差分格式无条件稳定 0 0 0 0 0 0 0 uu a xtTx tx u xxx utttT 3 对变系数方程对变系数方程的边初值问题建 立左偏心格式及 的边初值问题建 立左偏心格式及 LF 格式 并研究它们的稳定性 格式 并研究它们的稳定性 中国地质大学 北京 廉海荣编 5 第五章 双曲型方程的数值解 6 中国地质大学 北京 廉海荣编 6 解 首先用矩形网格剖分区域 取空间步长 时间步长h 的矩形网格剖 分 区 域 用 节 点表 示 坐 标 点 0 1 2 j j k jk x tjh k 0 1 2 T k 下面建立左偏心格式 方程在处取值 关于时间和空间的一阶偏 导数均用向前差商代替 舍去误差项可定义差分格式 j k 1 1 0 kkkk jjjj j uuuu a h h r 以上等式可化为 令 k jj k jj k j uauau 1 1 r r1 对于初边值问题可定义差分方程组为 0 0 0 1 2 0 1 2 jj k k uj uk 令代入左偏心格式 有 kkij j uv e h e 1 1r r ijhkijhkijhijhk jj evaeva eev 由此可得传播因子为 显然 1r r ijh jj Gaa 2 22 2 1cos sin 1 2 2 1 cos jjj jj Ga ra rha r a ra rh h 即时 2 jj a ra r 2 1G 又因为 0 0 j CaC1 即稳定性条件为 1 1rC 下面建立 LF 格式 首先采用特征线法 如图所示 第五章 双曲型方程的数值解 7 A D BC P k tt 1 k tt 为了求方程的解在 P 点的函数值 利用特征线法可知 u Pu D 故只需 求解函数在D点的函数值 ADC三点时间坐标相等 故可以把此三点方程 的解看成空间的一元函数 可以利用 AC 两点函数值的线性插值来求 D 点的函数值 记 1 1 1 1 1 1 1 1 1 jk jk jk jk Pj kx tjh k Ajkxtjh k Bj kx tjh k Cjkxtjh k CDDA DkAk ADCA xxxx u xtu xtu xt xxxx 中国地质大学 北京 廉海荣编 7 对点线性插值有CA Ck 即 1 111 1 1 1 222 jj kkkk jjjjjj ha rha r uuua r ua r hh 1 k j u 也可以用差商代替偏导数建立 LF 格式 关于时间的一阶偏导数用向 前差商代替 关于空间的一阶偏导数用中心差商代替 舍去误差项可定义 中心差分格式 1 11 0 2 kkkk jjjj j uuuu a h 此等价于 1 11 2 kkkk jjjjj r uua uu 1 1 2 kkk jjj uuur h 1 用代替代入上式得 LF 格式 其中 第五章 双曲型方程的数值解 8 中国地质大学 北京 廉海荣编 8 h h 令代入 LF 格式 可得传播因子为 kkij j uv e cossin jj Ga rhia r 2 222 1 1 sin j Ga r 2 1G h 容易看出1 j a r 时易得 又因 为 所以格式稳定的充分条件是 1 1rC 0 0 j CaC 1 0 0 uu ax tx a 2 sin1 10 ih G 1 ra 2 hn 时有 下面用归纳法求 21 10 i G k G 第五章 双曲型方程的数值解 10 中国地质大学 北京 廉海荣编 10 i i 2 3 4 212132 101021 322143 211032 432154 321043 1 1 2 3 1 kk ii G i iii G ii ii G ii kki Gik kik k G k G 显然关于无界 那么矩阵族k不是一致有界 则当 差分格式不稳定 当1 ra0a aar的情况下是稳定的 第五章 双曲型方程的数值解 11 2 对每一层的最后一个点 即 N 1 号点亦利用左偏心格式求解 编程求得不同网格比 r 下得到的图形 分析 由于设边值为 0 又 a 0 那么根据特征性的性质随着 t 的增加函 数 u 的值应当趋近于 0 显然在 r 1 4 1 2 都是满足这一性质的 当 r 2 时 u 的值在后期发生严 畸变 这说明当 r 1 时差分格式并不稳定 实际上这里要着重讨论的是 r 1 时的情况 但很可惜在所举例子中这 时本应当表现出的不稳定性并不是十分突出 甚至可以说并没有体现 经 过分析 认为这可能是由于初 边值的取值不当造成的 另外由于第一层 及每层最后一个点的计算采用了左偏心格式 这也是可能导致错误的一个 中国地质大学 北京 廉海荣编 11 第五章 双曲型方程的数值解 12 中国地质大学 北京 廉海荣编 12 地方 当然 也不排除程序的编写错误 附程序 clear h 0 01 tao 0 03 fain 0 1 edge 0 a 1 xmin 0 xmax 1 tmin 0 tmax 1 r tao h N fix xmax h T tmax tao r0 1 4 这里特别设置r0 使满足左偏心格式稳定的条件 value 1 1 N 1 fain 赋初值 value 1 T 1 1 edge 赋边值 for j 2 N 1 value 2 j r0 a value 1 j 1 1 r0 a value 1 j 利用左偏心格式 设置第一层的数值 由于编程的方便 这里下标设置为2 end for k 2 T for j 2 N value k 1 j value k 1 j a r value k j 1 value k j 1 利 用 leap frog格式计算 end value k 1 N 1 r0 a value k N 1 r0 a value k N 1 利用左偏心 格式计算每层第N 1个点 end X Y meshgrid xmin h xmax tmin tao tmax mesh X Y value xlabel x ylabel t zlabel u x t hold on 第五章 双曲型方程的数值解 13 surf X Y value 22 2 0 uu a 22 a tx 5 试研究波动方程5 试研究波动方程为正常数的加权格式的稳定性 为正常数的加权格式的稳定性 2 22122 22 1 12 kkk tjxjxjxj a uuu h 1 k u 解 将二阶中心差分展开 有 2 11111 11 22 111 111 1 2 2 12 2 2 kkkkkk jjjjjj kkkkkk jjjjjj a uuuuuu h uuuuuu 1 k1k jj vu 令 那么格式等价于 122111 11 111 1 2 2 12 2 2 kkkkkk jjjjjj kkkkkk jjjjjj kk jj uuva ruuu uuuvvv vu 1 得 kkk jjj wuv 令 2222 111 1 222222 1 11 22222222 1 020 0000 0 12 0000 22 12 12 12 1000 kkk jjj kk jj kk jj a ra r www a ra ra r ww a ra ra ra r ww ijhkk j wq e 令 那么有 其中传播因子 1 k Gq k q 中国地质大学 北京 廉海荣编 13 第五章 双曲型方程的数值解 14 中国地质大学 北京 廉海荣编 14 22222222 2222 2222 2222 2222 2222 22 12 2 12 cos122cos 10 1 0 122cos 01 22 12 2 12 cos 122cos 122cos 1 0 122cos a ra rha ra rh G a ra rh a ra rh a ra rh a ra rh a ra rh 其特征值方程为 2222 2 2222 22 1 2 2 1 2 cos 10 122cos a ra rh G a ra rh 特征值按模不大于 1 的充要条件是 2222 2222 22 12 2 12 cos 12 122cos a ra rh bc a ra rh 成立 而此充分条件为 22222222 22222222 22 12 2 12 cos244cos 22 12 2 12 cos244cos a ra rha ra rh a ra rha ra rh 2222 22cosa ra rh 0 第一个不等式等价于 0 故对于任意的参数其都 成立 而第二式化简为 1 4 22 2 1cos 41 4a rh 当时 此不 等式成立 而当 1 0 4 22 2 1 cos 14 a r h 1 14 ar 时 若即 时 不等式成立 进一步可以证明此就是格式稳定的充分条件 第五章 双曲型方程的数值解 15 vw a tx wv a tx 6 对方程组对方程组 试建立 试建立 LF 格式 并研究格式的稳定性格式 并研究格式的稳定性 中国地质大学 北京 廉海荣编 15 解 用矩形网格剖分区域 取空间步长 时间步长h 的矩形网格剖分区 域 用节点表示坐标点0 1 2 j 0 1 2 k j k jk x tjh k 方程在处取值 关于时间的一阶偏导数用向前差商代替 关于空间 的一阶偏导数用中心差商代替 舍去误差项可定义中心差分格式 j k 1 11 1 11 2 2 kkkk jjjj j kkkk jjjj j vvww a h wwvv a h 此等价于 1 11 1 11 2 2 kkkk jjjj kkkk jjjj ar vvww ar wwvv 1111 11 22 kkkkkk jjjjjj vvvwww r h 用其中代替代入上式 得 LF 格式 1 1111 1 111 1 22 1 22 kkkkk jjjjj kkkk jjjj ar vvvww ar wwwvv 1 k j 1 k k j v mr wh 1 变为 设 第五章 双曲型方程的数值解 16 1 111 10011 011022 kkkk jjjj ar mmmm 1 k j m h 令 则有 kkij j mq e 1 10011 011022 cossin sincos kijhi hi hkijhi hi hkijh kijh ar qeeeq eeeq e hiarh q e iarhh 故其 LF 格式的传播因子为 cossin sincos L hiarh G iarhh 中国地质大学 北京 廉海荣编 16 于是 2222222 cossin1 1 sin L Gha rha rh 1ar 所以 当 1 L G LF 格式稳定 时 7 试研究二维波动方程的显格式的稳定条件 7 试研究二维波动方程的显格式的稳定条件 122 1 1 1 1 2 12 kkkkkkk i 1 jijijiji ji jij ururuuuuu 解 这是一个二维三层格式 利用传播因子法分析格式的稳定性 首先将 其化为等价的双层格式组 有 12 2 12 kk ux yrux y 2 第五章 双曲型方程的数值解 17 中国地质大学 北京 廉海荣编 17 8 考虑带阻尼的波动问题考虑带阻尼的波动问题 22 2 22 2 01 0 0 0 01 0 1 0 0 uuu aCxtT ttx u xxu xxx t ututtT 22221 2 kkk tjxjjj ua ruC uu r h 当时是稳定的 其中当时是稳定的 其中1ar 2 而当时 必然有0 此时 2 222222 1 2 12sin12sin 12 22 hh Ca rCa rC 中国地质大学 北京 廉海荣编 20 2 1 若 222 12sin 2 h Ca r 0 22 1 2 C a r 即 6 1C 此时有 2 222222 12sin12sin 21 22 hh GCa rCa rC 2 222222 12sin12sin 22 hh Ca rCa r 1 von Neumann 条件成立 2 2 若 222 12sin 2 h Ca r 0 22 1 2 C a r 即 2 222222 12sin12sin 21 22 hh GCa rCa rC 2 222222 12sin12sin 22 hh Ca rCa r 222 24sin2 2 h a rC 要使 von Neumann 条件成立 那么 第五章 双曲型方程的数值解 21 222 24sin1 2 h a r 7 联立 6 7 可得稳定的必要条件 Gr 5 若矩阵 1 中国地质大学 北京 廉海荣编 21 的两个互异特征值 2 那么相应地有两个线性无关 的归范特征向 e1 e2 11 22 1 1 e 21 T iAC iAC 21 22 22 2 1 e 21 T iAC iAC 21 而以 e1 e2 为列向量的行列式之模等于 则定理 4 2 2 可知 格式稳定 当时 取 1ar 综上所述 差分格式的稳定性条件是 ar 1 方法二 类似于抛物型多层差分格式稳定性讨论 将三层格式转换为双层格式组 221 2 kkkkkk tjxjjj ua ruC uu k j 展开得 11221 11 2 2 2 kkkkkkk jjjjjjj uuua ruuuC uu k j k j uv 1 令 有 第五章 双曲型方程的数值解 22 中国地质大学 北京 廉海荣编 22 kk1222222 11 1 12 C222 kkk jjj kk jj ua r uCa rua r u vu jj v k j k j k j u w v 写成矩阵形式 令 222222 1 11 120222100 01100000 kkkk jjj CCa ra ra r www j w 1 求显格式组的传播因子 二阶矩阵 kkij j wq e h 可得传播矩阵为 22 222 1 cos 1 1212 10 Ca rh G CC 3 求传播矩阵的特征值 判别 von Neumann 条件 22 2 222 1 cos 1 0 1212 Ca rh CC 特征方程 2 1 按模不大于 1 的充要条件是 特征根 22 222 1 cos 2 12 Ca rh C 当 时上式成立 此时 von Neumann 条件满足 1 ar 4 利用定理 4 2 2 给出稳定的充分条件 第五章 双曲型方程的数值解 23 中国地质大学 北京 廉海荣编 23 2 12 22 222 22 2 2 1 1 1 cos 12 1 1 cos 12 1 12 1 1 cos 12 1 12 0 1 1 arG ca rhc ca rhcc Aca rhc Bc B AAB E AAB 1 2 所以对于特征值 的特征向量e 1 同理对于特征值的特征向量 1 e 故由4 4节的稳定性判别定理知差分格式稳定 时 当 1ar 2N 221 1212 10 C GCC 2 221 0 1212 C CC 0C 时 1 2 1 时特征方程有特征根 特别地 21