数学物理方程第九章 数学物理中的近似解法.pdf

第 9 章 数学物理中的近似解法 我们在前面几章讨论数学物理问题时 求得的都是解析解 或说是精确解 这些解法一 般只对较简单的问题和较规则的区域才有效 但实际问题往往很复杂 这时要解出精确解就 很困难 有时甚至不可能 另一方面 在建立数学模型时 我们已作了很多近似 所以求出 的精确解也知识推导出的数学问题的精确解 并非真正实际问题的精确解 因此 我们有必 要研究近似解法 只要使所求得的近似解与精确解之间的误差在规定的范围内 则仍能满足 实际的需要 电子计算机技术的飞速发展 为求得近似解提供了极大的可能 近似解一般可 分为近似解析解和近似数值解 近似解析解是指对解析解在规定误差范围内的近似 近似数 值解是指解在讨论的区域上 足够多 的点处的近似值 9 1 解析近似解 9 1 1 正则摄动法求解非线性偏微分方程 由于大多数非线性方程无法求得其解析解 这就使我们不得不借助于解析近似法来求其近似 解析解本节将介绍一种求非线性方程近似解析解的正则摄动法 摄动法的解题思想是 如果非线性方程的非线性项是高阶小量 则作为初步近似 可将其略 去 从而非线性问题便化为了线性问题 求解该线性问题 并将所得的解作为非线性问题的 零级近似 再把原非线性定解问题的解看作它的零级近似解与一个待求的含有小参量的解的 和 代入原定解问题 略去更高阶小量 得到关于摄动解的线性定解问题 求解该线性问题 并将求得的解作为原定解问题的高一级的近似 一级近似 仿此步骤进行下去便可得到各级 近似解 这种求近似解析解的方法被称为摄动法或小参数法 而以上这种带有小参量的定解 问题被称为摄动问题 我们用一个具体的例子来介绍正则摄支法解题的步骤 例 1 试在单位圆内求解定解问题 cos 10 1 0 1r ru r u r u u 9 1 1 解 这是一个含有小参量的非线性问题 当0 时 显然 式 9 1 1 变为了线性定 解问题 cos 1 0 1r ru ru 9 1 2 由分离变量法即可求得其解 视为零级近似 并记作 0 u 即 cos 0 ru 9 1 3 又设 110 cosuruuu 9 1 4 将式 9 1 4 代入式 9 1 1 并略去 2 o项 得 0 2sin 2 1 11 1 r u ru 9 1 5 这是泊松方程的狄氏问题 我们采用下述方法来求解 由式 8 1 5 右边函数的形式 我们 猜测方程的解应为 2sin 1 rfru 式中 rf为待定函数 将之代入式 9 1 5 得 7 1 9 0 0 1 6 1 9 2 1 4 32 ff rrfrf rrfr 这是欧拉方程 我们可求得一个特解为 3 10 1 rrf 则式 9 1 6 的通解为 32 2 2 1 10 1 rrcrcrf 由边界条件式 9 1 7 求得 0 10 1 21 cc 即 32 10 1 10 1 rrrf 则 2sin 1 10 2 1 r r u 故定解问题的一级近似为 2sin 10 cos 23 10 rrruu 一般地 可令 L 2 2 10 rururuu 代入式 9 1 1 中 比较 的同次幂系数 有 cos 0 10 0 r u u 0 11 00 1 r u u r u u M 0 12 010 2 r u uu r u u 这是一系列泊松方程的狄氏问题 依次求解每一个定解问题 0 2 1 1 110 rk kk u kuuufuLL 从而可得到定解问题式 9 1 1 的和级近似解 用以上的摄动法求定解问题时 若小参量 是方程中最高阶偏导数的系数 则求解会遇 到困难 因为此时略去含 的高阶项后 会降低方程的阶数或改变方程的类型 因此所得的 线性的近似方程的解将不可能满足原来定解问题的边界条件 因而上面所叙述的正则摄动法 失效 需考虑所谓的奇异摄动法 这已超出了本课程的内容 在此不再赘述 有兴趣的读者可查阅相关书籍 9 1 2 积分方程的近似解 在研究实际问题时 我们经常遇到积分号内出现未知函数的方程 称之为积分方程 例 如 b a xfdtttxkx 若fk 为已知函数 为参数 为未知函数 则它就是一个积分方程 称之为第二类弗 雷德霍姆方程 其中 txk称为核函数 大部分积分方程难以求出其精确解 下面介绍如 何求出方程的解析近似解 在第二类弗雷德霍姆方程 b a xfdtttxkx 9 1 8 中 若核函数具有形式 n i ii tbxatxk 1 其中 tbxa ii 都是线性无关的实函数 组 则称其为退化核 对于退化核的积分方程 可将之化为代数线性方程组求解 事实上 若 txk为退化核 则式 9 1 8 为 b a i n i i xfdtttbxax 1 9 1 9 令 b a ii dtttbZ 则式 9 1 9 化为 n i ii Zxaxfx 1 将其代入积分式 9 1 8 得 n i n k kk b a iii dtZtaxftbZxa 11 0 因 2 1 nixaiL 线性无关 故 2 1 0 1 nidtZtatftbZ b a n k kkii L 记 b a ikki b a ii adttatbfdttftb 则有 2 1 1 nifZaZ n k ikiki L 9 1 10 式 9 1 10 是一个关于 2 1 niZiL 的线性方程组 易见 若式 9 1 9 有连续解 x 则 i Z必满足式 9 1 10 反之 若式 9 1 10 有解 i Z 代入方程 n i ii Zxaxfx 1 即得式 9 1 9 的连续解 x 利用退化核方程的求解方法 可求任意核积分方程的近似解 把已知的核 txk用与 它相近的退化核 txk代替 把相应的具有退化核方程的解 x 作为原积分方程的近似 解 一般可取 txk的泰勒 Taylor 级数的有限项 txk 用退化核逼近任意核方法求出 的是近似解析解 我们以具有例子说明此方法 例 2 求方程 1 0 sin cos1 xdttxtxx 的近似解 解 用 Taylor 级数表示核函数xtxtxkcos1 则 L 4 2 1 4523 txtx xtxk 取 2 1 23t x xtxk 考虑方程 求解此具有退化核的方程 得 3 21 1 sin xcxcxx 式中 dtttCdttC 2 1 1 0 2 2 1 0 1 由 的表达式 得 1cos 2 1 11sin 12 1 24 1 sin 2 1 1cos1 4 1 2 1 1 sin 21 1 0 5 2 32 1 2 2 21 1 0 3 211 CC dttCttCttC CC dttCtCtC 解得 111cos51sin12 43 4 3 8 1cos 6 19 1 sin 43 24 21 CC xxCxCxsin 1 3 21 就是所求积分方程的近似解 9 2 数学物理方程的差分解法 我们以二维拉普拉斯方程的定解问题为例 介绍一种数值近似方法 差分解法 其基 本思想是用函数的差商代替微商 从而把微分运算化成代数运算 求出解在定解区域中足够 的点上的近似值 由此得出近似解 由微分法知道 函数 xf的导数的增量与自变量增量的比值在自变量增量趋于零的极 限 即 x xxfxxf x xxfxf x xfxxf xf x xx 2 lim lim lim 0 00 我们把函数的增量f 与自变量x 的比值称为差商 因此 当x 很小时 可用下列三 种差商中的一种代替微商 x xfxxf xf 9 2 1 x xxfxf xf 9 2 2 x xxfxxf xf 2 9 2 3 我们称式 9 2 1 9 2 3 中的三个差商分别为 xf在x处的向前 向后和中心差商 由 于它们都是一阶导数的近似 所以称之为一阶差商 当x 固定时 一阶差商也可看成x的 函数 再考虑二阶或更高阶差商 它们是二阶微商或更高阶微商的近似 比如用一阶向前差 商的向后差商 或一阶向后差商向前进差商 来近似二阶微商 求得 x x xxfxf x xfxxf xf 2 2 x xxfxfxxf 对于一个多元函数 也可用差商近似偏导数 例如对函数 yxuu 可用对x 或 对y 的差商来近似代替对x 或对y 的偏导数 22 2 22 2 2 2 y yyxuyxuyyxu y u x yxxuyxuyxxu x u y yxuyyxu y u x yxuyxxu x u 因此 可用一个关于差商的方程 称为差分方程 来近似代替一个偏微分方程 对于二维拉普拉斯方程 0 2 2 2 2 y u x u 其差分方程为 0 2 2 2 2 y yyxuyxuyyxu x yxxuyxuyxxu 9 2 4 这种近似的替代所产生的误差可估计为 4 4 4 4 3 24 4 32 4 1 3 1 2 1 4 1 3 1 2 1 2 xfxxfxxf xxfxxfxxf xxfxxf xxfxxfxfxxf 式中 1 0 2121 xxxx 于是我们得到 2 2 2 2 2 2 2 2 2 2 x y yxu y yyxuyxuyyxu x x yxu x yxxuyxuyxxu 因此 用差分方程 9 2 4 代替二维拉普拉斯方程 其截断误差 这种误差是由于用 Taylor 级数无穷项而造成的 为 22 yx 的数量级 则阶段误差为 2 h少数量级 考虑定解问题 26 9 25 9 0 2 2 2 2 yxfu yx y u x u 式中 表示平面有界区域 的边界 我们作平行于坐标轴的两族直线 jhyy ihxx i i 0 0 这两族直线将区域 分割成许多小方格 如图 9 1 所示 这些小方格称为网格 小方格 的边长h称为步长 以 h 表示由一些正方形网格的边所联成的封闭折线 如图 9 1 中粗线 所示 选取 h 应尽可能与原边界 接近 h 所包围的所有网格的全体记作 h 显然 h 是 一个与区域 近似的区域 点 ji yx称为网格的特点 所谓差分解法就是以 h 代替 以 h 代替 然后在 h 内部 不包括边界 h 上的 所有节点上求出定解问题式 9 2 5 式 9 2 6 解的近似值 要求到这个目的 必须做两件事 一是将原微分方程式 9 2 5 化成差分方程 即建 立解在 h 内任一节点 ji yx上近似值所满足的方程 二是将原定解问题中的边界条件式 9 2 6 确定解在 h 上各个节点处的近似值 在 h 的任一节点 ji yx处 分别用二阶差商 2 2 2 2 y yyxuyxuyyxu x yxxuyxuyxxu jijiji jijiji 代替 2 2 2 2 y u x u 用 ji U 表示u在 ji yx处的近似值 则对应于二维拉普拉斯的差分方程为 04 1 11 1 jijijijiji UUUUU 或者写成 4 1 11 1 jijijiji ji UUUU U 9 2 7 式 9 2 7 表明 任一内节点 ji yx上解的值等于其周围四个相邻节点上解的近似值的算 术平均值 9 2 这正是调和函数平均值公式的体现 图 9 1 图 9 2 再看边界条件的差分近似 由于边值 yxf定义在 上 在 h 上并未定义 所以要转 移到 h 上 对于任一边界节点 hji yx 取 ji yx 使 ji yx与 ji yx的距离最短 如图 9 3 令 jiji yxfU 9 2 8 图 9 3 这可作为差分方程式 9 2 7 的边界条件 这样 式 9 2 7 式 9 2 8 组成了边值问题式 9 2 6 式 9 2 6 的差分形式 hji jijijijiji jiji yxyxfU UUUUU 4 1 1 11 1 这是一个关于 ji U 的线性代数方程组 可以证明它的解是存在唯一的 我们可以利用消元法直接求解方程组式 9 2 7 由于步长h很小 导致 hh 上节 点很多 直接消去法往往难以奏效 所以通常都采用迭代法求解 这里 我们介绍两种迭代 法 同步迭代法和异步迭代法 1 同步迭代法 首先 任意给定网格区域 h 内节点 ji yx上一组数值 0 ji U作为解的零阶近似 代 入差分方程式 9 2 7 的右端 得到 4 1 0 1 0 1 0 1 0 1 1 jijijijiji UUUUU 将 1 ji U作为解的一次近似 右端四个值当中若涉及到边界点上的值 均用相应的已知 值 ji yxf代入 一般来讲 在已得到解的第k次近似 k ji U后 由公式 4 1 1 1 1 1 1 k ji k ji k ji k ji k ji UUUUU 9 2 9 得到解的第1 k次近似 这样就得到一个近似解的序列 2 1 0 kU k ji 可以证明 不论 零次近似 0 ji U如何选项 当 k时 此序列必收敛于差分方程式 9 2 7 的解 因此当 k相当大时 k ji U就给出所要求的近似值 通常在计算一定步数后 相邻两次迭代解 1 k ji U和 k ji U之间的误差小于某规定正数 即 1 max k ji k ji ji UU 可结束迭代过程 而取 k ji U作为所求近似解 2 异步迭代法 一般说来 同步迭代法的收敛速度是比较慢的 为了加快迭代程序的收敛性 常常采 用异步迭代法 所谓异步迭代法就是在计算第1 k次近似值 1 k ji U时 如果图 9 3 所示的四 个相邻结点中有些结点处的第1 k次近似值已经求得 就用这些值代替式 9 2 9 右端原 来的第k次近似值 在使用异步迭代法时 必须将网格区域的结点按一定的顺序进行排序 并逐个进行迭代 通常是按结点的自然顺序进行迭代 即在每一横排上从左到右依次进行迭 代 等这一所有结点全部做完了之后 再紧接着对上一排的所有接点用同一顺序进行迭代 如图 9 4 显然 在求结点 ji yx处的第1 k次近似值 1 k ji U时 其周围四个相邻结点 中有两个结点 1ji yx 及 1 ji yx处还只有第k次近似值 因此异步迭代法的相应迭代公 式为 4 1 1 1 1 1 1 k ji k ji k ji k ji k ji UUUUU 9 2 10 与同步迭代法类似 当此式右端涉及到边界结点上的值时 均用边界条件式 9 2 8 中所给 的已知值代入 由于在异步迭代法中有一半是用了迭代的新值 所以可以预计异步迭代法的收敛速度比 同步迭代法的收敛速度要快一倍左右 因此在求解拉普拉斯方程的定解问题时 异步迭代法 是一个常被采用的方法 图 9 4 图 9 5 例3 求边界为3 0 4 0 yyxx 边界条件如图 9 5 所示的拉普拉斯方程的狄 式问题的近似解 解 取1 h 由于结点个数不多 因此零次近似可取为边界的平均值 即 4582 0 12707 02290 4 1 0 ji U 然后采用异步迭代公式 4 1 1 1 1 1 1 k ji k ji k ji k ji k ji UUUUU 进行计算 9095 0 7071 024582 02 4 1 1 1 1 U 7056 0 124582 09059 0 4 1 1 1 2 U 4677 0 7071 04582 007056 0 4 1 1 1 3 U 8410 0 4582 09059 002 4 1 1 2 1 U 5012 0 07056 04582 07410 0 4 1 1 2 2 U 2422 0 4677 0005012 0 4 1 1 2 3 U 这样一次次地迭代下去 一直算到 16 ji U 可以发现这些值与 15 ji U相比 小数点后面三位 数字都相同 所以就可取 16 ji U作为原定解问题的近似解 2916 04616 0 8667 0 4168 0 7405 0 0836 1 16 2 3 16 2 2 16 2 1 16 1 3 16 1 2 16 1 1 UU UU UU 上面我们仅就第一类边界条件讨论了拉普拉斯方程的解法 如果边界条件不是第一类的条件 而是第二类或第三类边界条件 那么怎样用差分方法来解呢 当然在网格内节点处前面所列 的差分方程式 9 2 7 仍适用 现在的问题是如何处理网格边界节点上解的值 由于第二或第三 类边界条件中含有未知函数的法向导数 所以在用差分方法求解时 必须根据所给的边界条件 在边界节点 hji xxp 处列出相应的差分方程 和方程组式 9 2 7 一起联立求解 我们以第三类边界条件 yxfu n u 9 2 11 为例 说明如何列出边界节点上的差分方程 其中f 是 上的已知函数 n u 0 表示u 的外法向导数 图 9 6 对 h 上任一节点p 自p向边界 作法线 此法线与 相交于 p点 图 9 6 我们就以 p处的外法向作为在p点的n 设如此规定的法向矢量n与两坐标轴的夹角分别为 则 sincos cos cos y u x u y u x u n u p 9 2 12 将偏导数 x u 及 y u 分别用差商 h PUQU 及 h PURU 来代替 并以 pf代替 pf 以 p 代替 p 则由式 9 2 12 得 sin cos pfPUp h PURU h PUQU 或 sincos sin cos Ph PhfRUQU PU 9 2 13 式中 PU QU等分别表示解在QP 等点处的近似值 式 9 2 13 就是与边界条件式 8 2 11 相对应的差分方程 到此为止 我们已经解决了在正方形网格的情况下 怎样用差分方法求拉普拉斯方程的定解问 题 但是 读者应该明白 网格不一定非得用正方形的 为了使 h 更接近于边界曲线 我们可 以用矩形网格 平行四边形网格 正六边形网格等 即使采用正方形网格 所得的差分方程 也不一定只有式 9 2 7 这一种形式 只要用不同的近似式代替u 就会引出不同形式的 差分方程 此外 差分方程的解法也有多种 上面所叙述的迭代法只不过是其中的一种方法 关于这些问题 读者可以参阅其他书籍 9 3 积分方程的数值积分法 考虑第二类弗雷德霍姆方程 b a xfdtttxkx 9 3 1 假定 txk xf具有一阶连续导数 2 1 njxx j 得 2 1 njxfdtttxkx b a jjj b a n i iijij xxxkAdtttxk 1 2 1 1 njxfxxxkAx i n i iijij 9 3 2 式 9 3 2 为一线性方程组 解此方程组得 x 在各 j x处的近似值 ba上节点 j x取得越 多 1 jjj xxx越小 则误差越小 对不同的公式 jj xA 的取法不同 如对于矩形公式 取 2 1 1 nj n ab hhAhjax jj 而对梯形公式 取 1 2 1 121 n ab hhAA h AAhjax nnj 对Simpson公式 则取 bnhaxhaxax n 2 1221 hAAA h AA nn 3 4 3 242121 n ab hhAAA n 2 3 2 1253 用数值积分法求出的是近似数值解 对其他线性积分方程可用类似方法求数值解 例 4 求 x tx dt tx 2 sin1625 cos 2 38 6 的近似解 解 易见 xx 将12 等分 并记 j j 6 由对称性 有 jj 6 即只要 在 2 0 中确定 j 即可 用矩形公式得方程组 919 135 031 015 0 1318 034 132 016 0 2116 032 034 118 0 2515 031 035 019 1 3210 3210 3210 3210 解之得953 0 726 4 274 12 074 16 3210 实际上 方程的精解为 xx2cos 17 128 2 17 即970 0 2 735 4 3 265 12 6 03 16 0 由此可见误差不超过 2 附录 探讨定解问题的适定性 能量积分法 附录 探讨定解问题的适定性 能量积分法 在本书中 我们讨论了定解问题的求解 至于解的性质基本未讨论 定解问题的适定性 包含三个方面的内容 解的存在性 唯一性 稳定性 我们利用所谓的能量积分法 证明一 维波动方程定解问题是唯一的 也是稳定的 首先 我们导出长为l的弦振动时的动能和位能的表达式 然后再给出总能量 动能的 总和 所满足的不等式 由物理学知识 若质点的质量是m 在时刻t其速度是v 则它在t时刻的动能是 2 2 1 mv 现在考虑弦上的元素ds 当弦作微小横向振动时 dxds 它在时刻t速度为 t u 所以ds 在时刻t所具有的动能近似地为dxut2 2 1 其中 是弦的密度 一般来说 是x的函数 整个弦在t时刻的动能为 t tdx uU 0 2 2 1 1 再看弦在t时刻的位能 或称势能 所谓位能就是使弦变形时所做的功 假设弦不受 外力作用 则使弦变形的力只有张力 反抗张力所做的功就是弦的位能的增量 当弦的振幅 很小时 它的张力可以看作是一个常向量 其大小记作T 张力在位移方向的分量近似于 x Tu 在这个力的作用下 弦变形了 其位移的增量为dxuduu x 弦上元素ds在t时 刻的位能近似为dxTux 2 2 1 整个弦在t时刻的位能为 dxTuV l x 0 2 2 1 2 当然 如果除了张力T以外 在t时刻弦在位移方向还受到密度为 txf的外力作用 这时 位能应为 l x dxfuTuV 0 2 2 1 3 将式 1 和式 2 或式 3 相加 即得弦在t时刻的总能量 l xt dxTuutE 0 22 2 1 4 l xt dxfuTuutE 0 22 2 1 2 1 2 1 5 如果 是常数 并不计常数因子 式 4 可以表示为 l xt dxtxuatxutE 0 222 6 式中 T a 2 或者更简单地写成 l xt dxtxutxutE 0 22 7 我们把式 6 或式 7 称为能量积分 或称之为u的能量模 有时将其正平方根称为u的 能量模 现在来考虑初值问题 0 0 0 2 2 2 2 2 xxu xxu txtxf x u a t u tt t 8 设u是式 8 的解 古典解 为了建立能量不等式 如图 1 我们过点 00 tx作特征线 00 ttaxx 它们与x轴相交于 0 00 atx 及 0 00 atx 这两条特征线与x轴 所围成的区域记作K 任取 0 0t 令 0 tKKI 侧边为特征线的梯形 0000 taxtaxtK I 区间 图 1 在式 8 中的波动方程两端同乘以 t u 并在 K上积分 得 KK dxdt t u fdxdt x u a t u t u 2 2 2 2 2 9 先计算上式左端的积分 由于 2 2 2 2 1 t u tt u t u 2 2 2 2 1 x u tx u t u x x u t u xx u t u xx u t u 代入式 9 可得 KK dxdt t u fdxdt x u t u x a x u a t u t u 2 1 2222 利用格林公式得 KK dxdt t u fdx x u a t u dt x u t u a 2 1 2222 10 式中 K 表示 K的边界 它由上底 下底 0 及两侧 1 与 2 所组成 把式 10 左端记作J 则 21 0 2 1 2 1 2 1 2222 222222 U dx x u a t u dt t u x u a dxadx x u a t u J x 11 将式 11 右端第三项记成 1 J 现在来估计 1 J 在 1 上 adtdx 在 2 上 adtdx 故 0 2 2 2 1 2 1 0 2 0 2 2222 2222 1 00 00 2 1 dt x u a t ua dt x u a t ua dx x u a t u adt t u x u a dxa x u a t u t u x u aJ ttaxx ttaxx 12 由式 10 式 11 与式 12 可得 K x dxdt t u fdxadx x u a t u 2 0 222222 13 利用代数不等式 22 2baab 可得 KK dxdtf t u dxdt t u f 2 22 从而 KK x KK x dxdtfdxdt x u a t u dxa dxdtfdxdt t u dxa dx x u a t u 2222222 22222 222 0 0 14 令 0 222 0 222 222 00 00 t dxdt x u a t u dxdt x u a t u dxdt x u a t u G ttax ttax K 由式 14 可知 G满足微分不等式 FG d dG 15 式中 K x dxdtfdxaF 2222 0 为了从式 15 中解出 G 用 e乘其两端 得 FeGe d d 在 0 上积分 得 1 0 eFdttFeGe t 故 1 FeG 16 将式 16 代入式 15 的右端 得 Fe d dG 再由 G的表达式可得能量不等式 2222222 0 K x dxdtfdxaedx x u a t u 17 注 1 上面我们得到了一维波动方程的能量积分及能量不等式 完全类似地可以得到弹 性膜或弹性体振动时的动能及位能的表达式为 动能 dxuU t 2 2 1 位能 dxu T V 2 2 式中 表示弹性物体所占的空间区域 x表示二维或三维欧氏空间的向量 21 xxx 或 321 xxxx u 表示u的梯度 例如 在三维空间内 k x u j x u i x u u 321 故 2 3 1 2 i i x u u 此外 用推导式 17 的同样方法也可以得到高维波动方程的能量不等式 所不同的是过空 间内一点向下作特征锥面 例如 对二维波动方程 过点 000 tyx作特征锥面 2 0 2 00 yyxxtta 而不是作特征线 注 2 不等式式 16 称为 Gronwall 不等式 它说明只要非负函数满足微分不等式式 15 即可得到式 16 更一般的情形是 如果非负函数 G在 0 上连续可微 0 0 G 且对 0 T 它满足 FcG d dG 则 1 1 Fe c G c 证明的方法和推导式 16 完全相同 只要用 c e 代替 e 由于一维波动方程及其初始条件都是线性的 要证明式 8 的解是唯一的 只要证明 齐次问题 0 0 0 0 0 2 2 2 2 2 xu xu tx x u a t u