分片近似法求解数学物理方程反问题.pdf

1 分片近似法求解数学物理方程反问题分片近似法求解数学物理方程反问题 童海滨 1 1河海大学水资源环境学院 210098 E mail lbthb 摘摘 要 要 指出了基于优化理论的数理方程反演方法的局限性 分析了直接法的优点和缺陷 从弥补直接法的缺陷出发 借鉴有限元的思想并利用双差准则的优势 提出了数理方程反演 问题的一种新的直接解法 分片拟合法 论述了该方法的基本原理和需要解决的关键问 题 近似函数的收敛性和系数确定过程中的计算量以及相应的对策 采用分片拟合以 降低计算量 采用双差准则以防止近似函数的发散 随后列举了一个常微分方程反演的例子 以说明该方法的构造和求解过程 与基于优化理论的反演方法相比 本文提出的方法不需进 行微分方程的正演过程 数值正演或解析正演 不需经过参数的试算过程 即可直接由观 测数据一次求出所需反演的参数 数值实验的结果表明了本文所提算法的精确性和有效性 关键词 关键词 数学物理方程 反演 优化 分片近似法 双差准则 1 1 引言 引言 若在两个问题中 一个问题的表述或处理涉及到或包含了有关另一个问题的全部或部分 的知识 我们称其中一个为正问题 另一个反问题 1 微分方程是数学里的一个经典分支 在自然科学的各个领域中存在着各种各样的微分方 程 对微分方程中的计算问题也分为两类 1 已知方程的初始条件和 或 边界条件以及 各种参数 求方程的解 数值的或解析的 2 已知方程的解 求方程的初始条件和 或 边界条件以及方程中的未知参数 各种应用领域的 预测 预报 等问题常常可以归结为 微分方程的正演问题 而 参数估计 和 模型识别 问题常可归结为微分方程的反演问题 而更多的例子则同时涉及这两个过程 如 控制 领域的某些问题 依据常规思路 微分方程的反演问题可以看成一个优化问题 但换一种思路 微分方程 的反演问题也可以看成一个方程求解问题 依据前一种思路 首先要面对正演问题 无论是 数值的 还是解析的 其次要面对优化问题 在很多情况下 这意味着复杂的解析计算和 或 数值计算量 以及对最终结果是否是最优值的怀疑 而第二种思路 如果不追求过于 完美的解析式和 或 过于精确 实际上达不到 或不需要 的数值解 就给新方法的开辟 拆除了人为的障碍 2 基于优化理论的微分方程反演方法极其局限性 2 基于优化理论的微分方程反演方法极其局限性 基于优化理论的反演方法 算法流程如图 1 的局限性包含在如下几个方面 1 正问 题的解析解一般不存在 而数值解的计算量又较大 2 如果需要反演的参数较多 则又需 解决高维优化问题 而高维优化 除非十分特殊的问题 目前还没有有效的求解方法 或者 说 在有限的计算量下 是不可能求得最优解的 3 解的不适定性 2 2 图 1 基于优化理论的数理方程反演方法的算法流程图 既然有那么多的缺陷 是否有更简便一些的方法呐 据前人的研究 还有所谓直接法 即根据实测点据的 x y 值通过数值近似计算直接代入原微分方程 解出所求参数 但此种 方法对观测数据提出了较高的要求 如果因缺少数据 而进行线性插值来补缺的话 则又会 给计算带来严重的后果 但在作者看来 这些缺陷是可以采取某些措施加以弥补的 如果从实用的角度来看 与 其花费巨大的计算量寻求一个不知是否适定的解 倒不如寻求一个计算量小且接近真值 或 尚算实用 的解 另外 由直接解法求参数 当观测数据少时 也并非只有 插值 一条路可走 即使插 值 如果处理得当 则亦不会 对计算带来严重的后果 正如有限元法通过分片插值克服 了传统的 Ritz Galerkin 方法在构造子空间时所遇到的困难 并且可以使得形成 Ritz Galerkin 方程的系数矩阵时计算量大为减少以及使系数矩阵变为一个带装稀疏矩阵 从而减少了计算 量和存储量 3 因而本文的思路仍然是直接法 并且主要的目的是构造对计算不会带来 严 重后果 并且 适定 的直接法 3 分片 段 近似法求解数学物理方程反问题 3 分片 段 近似法求解数学物理方程反问题 3 1 分片近似法的基本思想分片近似法的基本思想 由直接法 算法流程见图 2 求解数学物理反问题的最关键之处即在于如果近似解 y yh x 中包含大量的待定参数 并且逼近方法是采用的连续高次插值或拟合 则插值或拟合本 身可能是不收敛的 或者即使收敛 需要求解一个计算量庞大的 甚至病态的 代数方程组 才能求出近似解中所包含的待定参数 但正如有限元法通过分片插值解决了传统 Ritz Galerkin 法所固有的数值计算量大的难 题一样 在直接法求反问题时 仍然可以采用这种 分而治之 的策略 只不过有限元法是 面向正问题 而本方法是面对的反问题而已 另外插值或拟合的收敛性问题 可以采用特殊的优度准则加以处理 本文采用最小二乘 试探参数 微分方程 解析解 数值解 实 测 点 据 x 值 实 测 点 据 x 值 计算 y 值 实测 y 值 误差 计算 y 值 实测 y 值 误差 调整参数 调整参数 3 法和双差准则联合实用的方法 由 观测 数据拟合出所需的 最佳 多项式函数 4 图 2 直接法求解数理方程反问题的算法流程图 3 2 一个算例一个算例 3 2 1 正问题正问题 为了说明这种分片近似法求解逆问题的具体操作步骤 这里举一个常微分方 程的例子 正问题如下 已知一阶微分方程的初值问题 00 yy kxpkyk dx dy x 为常数 1 易知其解析解为 x kxxk dxexpekyy 0 0 2 现不妨设 k 1 p x sin x y0 0 则相应的初值问题的解析解解为 x xx dxexey 0 sin 3 3 2 2 反演问题反演问题 1 现在假定已知条件为 00 yy kxpkyk dx dy x 为常数 4 k 1 y xi yi i 1 n 试求 p x 3 2 3 算法流程算法流程 1 由观测数据 xi y xi i 1 n 作 y x 的逼近函数 yh x 本例采用多项式拟合 关 于多项式的最佳次数的选择 文 3 给出了用最小二乘法和双差准则分别确定多项式 的最佳系数和最佳拟合次数的算法流程 该算法防止了多项式拟合 插值在经验点处 拟合精度较高 而在非经验点处又容易发散的缺点 具有较高的精度和较好的稳健 性 实测数据 x 实测数据 y 插值 拟合 近似解 y yh x 微分方程 L y f 参数 4 2 把拟合出的 y yh x 和 k 1 代入微分方程 4 式中 直接解出函数 p x 的近似 函数 p ph x 因为 yh x 是多项式函数 所以求导运算也非常简单 3 把求出的 ph x 和 p x 在本例中即指 sin x 在区间 x1 xn 上作比较 或计算两 函数之间的距离 估计结果的精度 3 2 4 数值实验数值实验 按前面的算法在 maple10 上计算结果如下 分别取 xi i L i 1 10 以微分方程正演 的解析解计算出的 y xi 值作为相应的 观测 数据 yi 当 L 0 1 时 所得 观测 数据如 表 1 所示 按最小二乘法和双差准则根据上述 xi yi i 1 10 拟合所得的最佳多项式 yh x 的次数和系数见表 2 函数曲线见图 3 反演所得 p x 的近似解为多项式函数 ph x 各 项的次数和系数见表 3 表1 观测 数据 xi yi i 1 10 L 0 1 xi 0 1 0 2 0 3 0 4 0 5 0 6 0 7 0 8 0 9 1 yi 0 004833 334702 0 018666 75302 0 040500 96911 0 069338 69714 0 104186 1812 0 144059 2473 0 187980 4019 0 234989 1728 0 284143 3005 0 334524 0601 表 2 y x 的近似函数 yh x 各项所对应的系数和次数 L 0 1 xi x0 x1 x2 x3 ai 0 000141185314799867324 0 00167191325041083350 0 494800105936482104 0 161831842252137814 表 3 反演所得 p x 的多项式函数近似解 ph x 各项所对应的系数和次数 L 0 1 xi x0 x1 x2 x3 ai 0 001530727935 0 9912721250 0 0093045790 0 161831842252137814 为直观的比较 ph x 和 p x 在本例即 sin x 的接近情况 下面给出这两个函数在 区间 0 1 1 上的函数图像 如图 4 其中红色线条代表 sin x 绿色线条代表 ph x 图 3 函数 y yh x 的图像 L 0 1 图 4 ph x 和 sin x 的吻合情况 L 0 1 当 L 0 7 时 数值实验的各项结果如下 表 4 观测 数据 xi yi i 1 10 L 0 7 xi 0 7 1 4 2 1 2 8 3 5 4 2 4 9 5 6 6 3 7 0 yi 0 187980 4019 0 531039 7754 0 745255 9500 0 669010 2770 0 307935 4216 0 18315 96871 0 58075 91995 0 70156 73261 0 49060 42156 0 04800188683 5 表 5 y x 的近似函数 yh x 各项所对应的系数和次数 L 0 7 xi x0 x1 x2 x3 x4 x5 x6 ai 0 19222936472 8831870 0 59189338672 0155457 1507950 0 49480010593 6482104 0 07227188163 70973572 0 16183184225 2137814 0 00003344977 23159651932 表 6 反演所得 p x 的多项式函数近似解 ph x 各项所对应的系数和次数 L 0 7 xi x0 x1 x2 x3 x4 x5 x6 ai 0 3996640220 1 698521861 0 319221514 0 19905551930 0560973729 2 0 0034356003 79 0 000033449772315 9651932 图 5 ph x 和 sin x 的吻合情况 L 0 7 图 6 ph x 和 sin x 的吻合情况 L 0 3 因篇幅所限 下面仅给出当 L 0 3 0 5 1 0 1 2 1 6 1 7 2 0 时 ph x 和 sin x 曲线的 对比情况 图 6 12 图 7 ph x 和 sin x 的吻合情况 L 0 5 图 8 ph x 和 sin x 的吻合情况 L 1 0 图 9 ph x 和 sin x 的吻合情况 L 1 2 图 10 ph x 和 sin x 的吻合情况 L 1 6 6 图 11 ph x 和 sin x 的吻合情况 L 1 7 图 12 ph x 和 sin x 的吻合情况 L 2 0 由数值实验结果 我们看到在观测数据十分 稀疏 例如 L 1 时 时 反演出的 ph x 和 p x 仍然十分接近 这一方面表明 直接法 在运用适当时 还是有效的 另一方 面 也表明 双差准则 在防止 插值 或 拟合 的发散方面 有较高的性能 当数据量特别大时 例如观测数据有 1000 个 那么高次插值或高次拟合 将面临较大 的计算量 但通过 分而治之 的办法 比如十个相邻的数据一组 在每组里面仍然可以采 用本文的算法处理之 从而避免一次处理成千上万个数据所带来的计算量大和 病态矩阵 等问题 3 2 5 反演问题反演问题 2 在 3 2 2 中所提到的反演问题是已知观测数据 xi yi 和 k 求 p x 如果现在换一下已 知条件 即已知 p x 和观测数据 xi yi 求 k 基本思路和 3 2 2 节一样 先根据观测数据 xi yi 拟合出 y x 的近似函数 yh x 然 后把 yh x 和 p x 代入微分方程 即可得 xyxp dxxyd k h h h 本来 k 是常数 但通过这种方法反演出的 k 的近似值 kh却是一个函数 把 xi yi 代 入上式 即可得 n 个 k 的估计值 可以取其均值作为最终估计值 满足 2 1 kk n i i h 最小 也可以建立目标函数 xyxp dxxyd kk h h 1n xxx 求使目标函数在 x1 xn 上取得最小值的 k 值作为 k 的最终取值 其中 表示某种范数 在本例中 由 xi yi 计算 yh x 的方法和 3 2 2 中是一样的 过程雷同 不在赘述 4 总结 总结 分别叙述了求解物理方程反问题的两类方法优化法和直接法各自的优缺点 在分析了通 常导致直接法缺陷的原因后 借鉴有限元法的思想和双差准则的优点 提出了求解数理方程 逆问题的一种新的直接法 分片近似法 随后以一个简单的常微分方程的反演问题为例 说明了该算法的基本思路和主要计算流程 和基于优化的反演方法相比 它无须多次 试算 和 调整 参数 一次 即可达到最终目的 数值实验结果表明了该算法在 观测 数据 十分 稀疏 的情况下依然保持了较高的精度 并附带验证了双差准则的优越性和实用性 7 致谢致谢 本文所述问题的原型来源于岩土工程领域 相关信息由河海大学土木工程学院硕士生 葛国昌提供 在此特致感谢 参考文献参考文献 1 Keller J B Inverse Problems Am Math Mon 83 107 118 1976 2 肖庭延 于慎根 王彦飞 反问题的数值解法 北京 科学出版社 2003 8 13 3 陈玉田 偏微分方程数值解法 南京 河海大学出版社 1999 231 232 4 童 海 滨 双 差 准 则 及 其 在 多 项 式 拟 合 和 水 文 频 率 计 算 中 的 应 用 Partitioned Approximating Method of Solving Indirect Problems of Mathematical Physics Equation TONG Hai bin1 1College of Water Resource and Environment Hohai Univ Nanjing 210098 China 1 Abstract The limitation of inversion method that based on optimization theory is pointed out The advantage and disadvantage of direct method is analized In oder to cover the shortage of the direct methods a new method partitioned approximating method of solving the indirect problems of mathematical physics equation is proposed after borrowing ideas from finite element method and making use of double difference criterion s merits The priciples of the method and the key problems approximating function s convergence and computation time of dertermining the coefficient that need to solve is discussed The corresponding countermeasure is presented adopt partitioned approximation to avoid the computing dificulty and adopt double difference criter