圲 坃坏坎坔坅坎坔坓 圲圮圴圮圱 坅坵坬坥坲方法圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮 圴圱 圲圮圴圮圲 龙格坼库塔方法圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮 圴圲 圲圮圴圮圳 两点边值问题圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮 圴

Transcription

1 Contents 0 绪论 7 1 计算机数和误差 11 圱圮圱 计算机数及其表示 圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮 圱圱 圱圮圲 舍入误差对计算的影响圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮 圱圲 圱圮圳 减法的计算 圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮 圱圳 2 物理学中的常用数值方法 15 圲圮圱 数值积分方法圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮 圱圵 圲圮圱圮圱 梯形法和辛普生方法圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮 圱圵 圲圮圱圮圲 反常积分的计算圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮 圱圷 圲圮圱圮圳 高斯积分方法圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮 圱圹 圲圮圱圮圴 哥西主值的计算圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮 圲圲 圲圮圱圮圵 急速振荡函数的积分圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮 圲圴 圲圮圱圮圶 高维积分的计算圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮 圲圶 圲圮圱圮圷 固体热容量的计算 圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮 圲圶 圲圮圲 方程的求根和最优化方法 圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮 圲圷 圲圮圲圮圱 二分法圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮 圲圷 圲圮圲圮圲 牛顿法和割线法圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮 圲圸 圲圮圲圮圳 一维最优化圬黄金分割法圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮 圲圹 圲圮圲圮圴 高维最优化圬共扼梯度方法圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮 圳圱 圲圮圲圮圵 约束最优化圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮 圳圵 圲圮圲圮圶 最小二乘法及曲线拟合圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮 圳圶 圲圮圳 函数的求值 圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮 圳圸 圲圮圳圮圱 多项式的求值和级数求和 圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮 圳圸 圲圮圳圮圲 连分式的求值圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮 圴地 圲圮圴 常微分方程的数值解法圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮 圴圱

2 圲 坃坏坎坔坅坎坔坓 圲圮圴圮圱 坅坵坬坥坲方法圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮 圴圱 圲圮圴圮圲 龙格坼库塔方法圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮 圴圲 圲圮圴圮圳 两点边值问题圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮 圴圶 圲圮圵 插值和逼近 圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮 圵地 圲圮圵圮圱 多项式插值圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮 圵地 圲圮圵圮圲 分式有理函数插值和外推 圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮 圵圱 圲圮圵圮圳 三次样条插值圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮 圵圲 圲圮圵圮圴 列表函数的积分圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮 圵圴 圲圮圵圮圵 坐坡坤圓坥插值与外推 圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮 圵圴 圲圮圵圮圶 发散级数的坃坥坳圓坡坲坯求和及其推广圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮 圵圸 圲圮圵圮圷 坂坯坲坥坬求和圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮 圶地 圲圮圵圮圸 坂圓坥坺坩坥坲逼近圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮 圶圱 圲圮圵圮圹 多元函数的插值及逼近圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮 圶圳 圲圮圶 快速付里叶变换圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮 圶圵 圲圮圶圮圱 付里叶变换圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮 圶圵 圲圮圶圮圲 离散付里叶变换圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮 圶圸 圲圮圶圮圳 快速付里叶变换圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮 圷圱 圲圮圷 附录场特殊函数的计算 圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮 圷圴 圲圮圷圮圱 误差函数的计算圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮 圷圶 圲圮圷圮圲 指数积分的计算圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮 圷圷 圲圮圷圮圳 整数阶贝塞尔函数及虚宗量贝塞尔函数的计算圮圮圮圮圮圮圮圮圮圮圮圮 圷圸 圲圮圷圮圴 球谐函数的计算圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮 圸圷 圲圮圷圮圵 椭圆积分的计算圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮 圸圹 圲圮圸 附录场高斯积分的节点和权重 圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮 圹圱 3 数值线性代数 97 圳圮圱 主元消去法解线性代数方程组圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮 圹圷 圳圮圲 坌坕分解法圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮 圹圸 圳圮圳 三对角方程组的求解圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圱地圱 圳圮圴 实对称矩阵的本征值和本征向量圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圱地圱 圳圮圴圮圱 坈坯坵坳坥坨坯坬坤坥坲方法 圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圱地圲 圳圮圴圮圲 坑坌算法圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圱地圴 圳圮圵 负本征值方法和递推方法 圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圱地圶

3 坃坏坎坔坅坎坔坓 圳 4 电磁场的计算 109 圴圮圱 坍坡坸坷坥坬坬方程圬边值问题圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圱地圹 圴圮圲 差分和差商 圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圱圱圱 圴圮圳 二维坐坯坳坳坩坯坮方程的五点差分格式圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圱圱圲 圴圮圴 差分方程的求解圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圱圱圵 圴圮圴圮圱 简单迭代法坼坊坡坣坯坢坩方法 圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圱圱圶 圴圮圴圮圲 逐次迭代法坼均坡坵坳坳坻坓坥坩坤坡坬迭代方法圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圱圱圶 圴圮圴圮圳 逐次超松弛迭代法在坓坵坣坣坥坳坳坩坶坥坏坶坥坲圭坒坥坬坡坸坡坴坩坯坮圩圮圮圮圮圮圮圮圮圮圮圮圱圱圶 圴圮圵 变分方法复习圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圱圱圸 圴圮圶 有限元方法的理论基础圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圱圲圱 圴圮圷 用有限元方法求解二维坌坡坰坬坡坣坥方程 圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圱圲圳 5 Monte Carlo 方法 129 圵圮圱 随机变量及其分布 圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圱圲圹 圵圮圲 赝随机数的产生圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圱圳圱 圵圮圳 用坍坯坮坴坥坃坡坲坬坯方法计算积分 圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圱圳圳 圵圮圴 自相似性与分形圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圱圳圶 圵圮圵 扩散限制粘结的计算机模拟圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圱圳圹 圵圮圶 高分子的模拟和无规行走 圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圱圳圹 6 逾渗和统计物理问题 143 圶圮圱 逾渗简介 圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圱圴圳 圶圮圲 逾渗的计算机模拟和分析 圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圱圴圶 圶圮圳 正则系综和统计模型圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圱圴圷 圶圮圳圮圱 坉坳坩坮坧模型圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圱圴圸 圶圮圳圮圲 坌坥坮坡坲坤圭坊坯坮坥坳模型 圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圱圴圹 圶圮圴 统计物理的坍坯坮坴坥坃坡坲坬坯模拟圬坍坥坴坲坯坰坬坩坳方法圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圱圴圹 圶圮圵 坉坳坩坮坧模型的坍坯坮坴坥坃坡坲坬坯模拟圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圱圵圲 圶圮圶 经典统计物理复习 圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圱圵圳 圶圮圶圮圱 坔坨坥坭坩坣坲坯坣坡坮坯坮坩坣坡坬坥坮坳坥坭坢坬坥在坎坖坅圩圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圱圵圴 圶圮圶圮圲 坔坨坥坣坡坮坯坮坩坣坡坬坥坮坳坥坭坢坬坥在坎坖坔圩 圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圱圵圴 圶圮圶圮圳 坔坨坥坎坐坔坥坮坳坥坭坢坬坥在坎坐坔圩圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圱圵圵 圶圮圶圮圴 坔坨坥坧坲坡坮坤坣坡坮坯坮坩坣坡坬坥坮坳坥坭坢坬坥在圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圱圵圵 圶圮圶圮圵 坔坨坥坲坭坯坤坹坮坡坭坩坣坡坬坲坥坬坡坴坩坯坮坳 圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圱圵圶 圶圮圶圮圶 坃坯坲坲坥坬坡坴坩坯坮坦坵坮坣坴坩坯坮坳圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圱圵圶 圶圮圶圮圷 坔坩坭坥坣坯坲坲坥坬坡坴坩坯坮坦坵坮坣坴坩坯坮坡坮坤坴坲坡坮坳坰坯坲坴坣坯坥圎坣坩坥坮坴坳圮圮圮圮圮圮圮圮圮圱圵圸

4 圴 坃坏坎坔坅坎坔坓 圶圮圷圶圮圸圶圮圹圶圮圱地 液体模型的坍坯坮坴坥坃坡坲坬坯模拟圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圱圵圹圶圮圷圮圱坍坯坮坴坥坃坡坲坬坯坓坩坭坵坬坡坴坩坯坮坯坦坔坨坥坌坥坮坡坲坤圭坊坯坮坥坳坍坯坤坥坬圮圮圮圮圮圮圮圮圱圶地 圶圮圷圮圲自由能计算圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圱圶圲 圶圮圷圮圳王圭坌坡坮坤坡坵方法及其它圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圱圶圴分子动力学方法圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圱圶圵 圶圮圸圮圱 均坥坮坥坲坡坬坰坲坯坣坥坤坵坲坥坯坦坍坄在坎坖坅坥坮坳坥坭坢坬坥圩圮圮圮圮圮圮圮圮圮圮圮圮圮圱圶圵 圶圮圸圮圲 坓坩坭坵坬坡坴坩坯坮坯坦坌坥坮坮坡坲坤圭坊坯坮坥坳坬坩坱坵坩坤坳 圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圱圶圹 圶圮圸圮圳 坓坩坭坵坬坡坴坩坯坮坯坦坈坡坲坤圭坳坰坨坥坲坥坳坹坳坴坥坭坳 圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圱圷地 坓坩坭坵坬坡坴坩坯坮坯坦坌坡坮坧坥坶坩坮坤坹坮坡坭坩坣坳圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圱圷圱 坌坯坮坧坲坡坮坧坥坦坯坲坣坥坡坮坤坅坷坡坬坤坳坵坭坭坡坴坩坯坮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圱圷圲 7 原子结构的计算 173 圷圮圱 原子结构问题圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圱圷圳 圷圮圲 变分法圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圱圷圴 圷圮圳 坖坩坲坩坡坬定理和坈坥坬坬坭坡坮坮圭坆坥坹坮坭坡坮定理 圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圱圷圶 圷圮圴 轨道近似和坈坡坲坴坲坥坥坻坆坯坣坫方程圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圱圷圹 圷圮圵 统计近似和 X α 方程 圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圱圸圱 圷圮圶 径向方程的求解方法圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圱圸圴 圷圮圷 计算程序 圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圱圸圸 圷圮圸 坓坴坵坲坭圭坌坩坯坵坶坩坬坬坥问题圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圮圲圱地

5 Chapter 0 绪论 不论是理论物理圬还是实验物理都离不开数值计算圬例如圬量子力学中简谐振子的波函数可以用厄米多项式来表示圬但为了得到波函数的数值圬仍然需要作数值计算圻实验上测量到的数据要进行数据处理圬也涉及到数据拟合圬分析并计算误差等一系列数值计算圮大家在过去的学习中已经知道圬牛顿力学方程只有二体问题是可解的在自由粒子或小振动等特殊问题除外圩圬三体以上的问题曾折磨了全世界的许多优秀的数学家和理论物理学家圬但最终也未得到解析解圻麦克斯韦方程组只有在一些非常特殊的条件下才能求得解析解圻量子力学的薛定格方程圬除了氢原子和简谐振子外圬还没有一个真实的物理问题可以找到解析解圻统计物理告诉我们圬只要求得了系统的配分函数圬则一切物理量都可通过求微分得到圬然而圬除了可数的几个二维模型外圬没有一个实际的物理系统的配分函数被求出来过圮一些大规模的物理实验圬耗资在数亿美元圬如果没有充分的论证就去做的话圬其浪费是十分巨大的圮对于这样一类问题, 大规模的数值计算往往可以发挥重要作用圮大规模数值计算也提供了理论物理和实验物理之间的桥梁圮几乎所有的实际物理问题都无法得到精确解圬发展各种理论模型和近似方法就成为理论物理学的重要方面圮理论模型的正确性最终需要经受实践的检验圬但对于一个确定的物理模型圬其对应的实验体系往往很难获得圬为了检验针对一个确定模型的近似方法圬最好能够得到这个模型的精确解圬而数值模拟方法正好能够提供这样的结果圮大量的实验体系往往包含多种因素圬这些因素有些可以排除圬而有些非常难以排除圮如果一个理论的出发点包含了所有的实验因素圬则这样的理论几乎肯定得不到解析结果圮但是圬数值模拟方法一般总能够把各种因素包含进去圬从而能够更好的模拟实验环境圮需要指出的是圬长期以来圬我们总有一些误解圬一讲到计算圬似乎就需要超级计算机圮事实上圬现在使用的笔记本计算机圬各种微型计算机的计算能力往往比圱地年前的超级计算机的计算能力更强圮用一台目前似乎已经被人看不起的坉坂坍圭坐坃圯坘坔微机在坉坂坍圭坐坃是基于坩坮坴坥坬芯片的第一代个人计算机圬大约圱圹圸地年推出圬其坃坐坕是坩坮坴坥坬圭圸地圸圶芯片圬有圶圴坫内存圬

6 圶坃坈坁坐坔坅坒地圮绪论两个圵吋软驱圬机器带有坍坓圭坄坏坓圱圮地和坃坐圯坍两个操作系统, 坉坂坍圭坐坃圯坘坔是坉坂坍圭坐坃的扩展, 增加了一个圱地坍的硬盘, 并且把内存增加到圲圵圶坋圩圬就可以很容易地解决三体问题圬四体问题圬圮圮圮圮所有元素周期表上的原子的能级和波函数也可以相当精确地计算出来圮如果你拥有一台目前市场上一般配置的微机圬你就可以模拟一些二维和三维模型系统的配分函数圬如坌坥坮坡坲坤圭坊坯坮坥坳流体圬三维坉坳坩坮坧模型等圬可以计算元素固体和简单化合物晶体的电子结构圬声子谱等圮进一步圬你如果拥有或能够使用到大型计算机圬或利用微机构成计算集群, 你就可以作真实流体的流体力学计算圬可以计算复杂边界条件下的三维电磁场圬可以计算三维真实系统的配分函数圬圮圮圮圮同样你也能够在大型物理实验开始前作很多模拟圬对实验的各种情况有一个清楚的概念圮本课程的目的在于对物理中的数值计算进行一些入门介绍圬使大家在学完本课程后圬在组织一些较大规模的计算时心中有数圬少走弯路圮课程分两部分圬一部分是基本算法圬主要介绍计算机数的特点圬舍入误差圬计算稳定性等基本概念和一些常用的数值计算方法圻第二部分则结合物理问题圬讲一圬二个专题圬使大家学到一些组织较大规模计算的步骤和技巧圮学习本课程要求学过四大力学在经典力学圬电动力学圬热力学和统计物理学及量子力学圩和数学物理方法圬掌握至少一门计算机高级语言在如坃圬坐坁坓坃坁坌圬坆坏坒坔坒坁坎圬坊坡坶坡等圩圮计算是一门实践性很强的课程圬上机实习是本课程的一个有机组成部分圮 关于本讲义的一点说明场这是笔者在圱圹圹圳年为开设 计算物理 课程准备的讲义圮圱圹圹圲年圬上海交通大学应用物理系代理系主任张仲渊教授找到我, 告诉我决定在本科生高年级开设 计算物理 课程圬并安排由笔者开课圮在接到这一任务后圬笔者翻阅了当时国内出版的几本 计算物理 教材圬发现都过于专门圬不适于做为本科高年级的教材圬而更适合于某一专业的研究生教材圮为此圬我只好试图自己写一本讲义圬以应急需圮基本算法部分主要参考了块圮坈圮坐坲坥坳坳圬坓圮坁圮坔坥坵坫坯坬坳坫坹圬块圮坔圮坖坥坴坴坥坲坬坩坮坧坡坮坤坂圮坐圮坆坬坡坮坮坥坲坹所著坎坵坭坥坲坩坣坡坬坒坥坣坩坰坥坳圬坔坨坥坁坲坴坯坦坓坣坩坥坮坴坩圌坣坃坯坭坰坵坴坩坮坧一书, 当时, 这本书出版时间还不算太长, 但已经成为大多数经常从事数值计算的物理学工作者的手头必备书, 我自己当时也买了一本, 书价是我的两个月的工资 这本书一方面讲解算法, 另一方面对每一个算法给出了透明, 简洁的源程序 尽管这些程序的效率, 安全性方面比一些成熟的软件包要差, 但由于其简洁清楚, 特别适合于修改后放到自己的程序中去 讲义的内容曾在上海交通大学应用物理系讲授多次圬各个年级的同学积极参与了本课程的学习并对课程的讲法提出了很多非常宝贵的意见和建议圬这些建议在讲义的成文中对笔者有很大帮助圬在此表示感谢圮圲地地地年后圬物理系决定把计算物理课程改为双语课程圬这本讲义就算完成了其历史使

7 圷命圮在教学过程中圬我也发现这一讲义对一些问题的讲解过于简洁圬而内容的深度也不太适合本科生的学习圮这一段时间圬国际上出现了几本比较合适的教材圬所以圬当时认为这一讲义不需要继续修订圬也就放在一边了圮随后圬研究生的课程做了一些改革圬增加了圳圶学时的数值计算圬任课老师感觉这个讲义的内容作为研究生数值计算的课程比较合适圬遂采用这本讲义的电子版教学圮圲地圱圳年圱地月, 致远学院的叶曦副院长和负责物理班教学的郑杭教授找我, 希望我为致远学院的同学开设一门物理中数值计算的课程 注意到致远学院已经有至少两门相关课程, 这门课程应该更多强调物理方面 为了这个课程, 利用寒假圬我修订了讲义圬主要是对原来的一些表述做了修改和补充, 同时尽可能强调物理方面的内容圮鉴于目前已经存在大量优秀的数学软件圬如坍坡坴坬坡坢圬坍坡坰坬坥等圮大量的小规模的计算往往可以通过这些软件方便进行, 因此, 此次修订时时也适当提到这方面的内容 讲义中提到的一些程序将放在课程主页上圬这些程序中, 有些来自前面提到的坎坵坭坥坲坩坣坡坬坒坥坣坩坰坥坳, 有些是我曾经多次使用过的圬也有些是专为此讲义而准备的圬由于时间关系圬没有经过仔细调试和考验圬因此在用于实际问题时圬一定要经过仔细试算圬以免造成不必要的损失圬同时圬我也愿在此声明圬笔者对因使用此讲义中的程序而导致的任何直接的和间接的损失将不负担任何责任圮所有程序都以坆坏坒坔坒坁坎圷圷写出圬鉴于坃语言在实际工作中使用的越来越广泛圬同时也是很多同学喜欢的语言圬部分程序已经翻译为坃语言圮由于程序是用坆坏坒坔坒坁坎写的圬翻译中自然留有坆坏坒坔坒坁坎的痕迹圮在上次讲义完成到这次修订之间, 面向对象的程序设计逐步走向主流, 其代表性语言是坣圫圫和坊坡坶坡, 建议有志于从事计算物理的同学注意这一点圮

8 圸坃坈坁坐坔坅坒地圮绪论

9 Chapter 1 计算机数和误差 在这一章中圬我们将给出电子数字计算机中数的表示圬计算规则圬舍入误差及算法稳定性 等概念圬在处理上圬并不追求数学上的严密圬而是用一些启发式的推导和例子来说明主要 概念圬需要进一步钻研的同学圬可以参看有关计算数学的书籍圮 1.1 计算机数及其表示 我们在初中的数学上就已经知道圬数有整数圬实数和复数圬整数可以取地圬 ± 圱圬 ± 圲圬... 圬 ±n 圬... 圬 ± 圬实数则是整个数轴上的连续统圮数学分析就是建立在实数系上的圮实数的运算满足加法和乘法的交换律圬接合律等一般的运算律圮而计算机所能表达的数的全体圬则是一个离散的圬有限的集合圮通常计算机整数的取值范围为 N 到 N 圱圬 N 的数值随计算机不同而异圬如坉坂坍圭坐坃机的 N 圳圲圷圶圸圮实数也有一个取值范围圬在坉坂坍圭坐坃上圬实数的 绝对值最大约为圱地圬最小约为圱地圬大于最大值的数将造成上溢圬小于最小值的数将造成下溢圮在在很多计算机系统中圬下溢作为地来处理圩圮计算机实数在其范围内的分布是非常不均匀的圬在靠近地的地方较密圬离开地越远越稀疏圮计算机实数一般不满足实数的运算律圮计算机中的实数通常用规格化的二进制浮点数来表示圬例如圱圱. 地地圱圬地. 地圱圱地地圱 和地. 圱圱地地圱分别表示为地. 圱圱地地圱 圲圬地. 圱圱地地圱 圲和地. 圱圱地地圱 圲圮由于计算机只有有限的字长圬所以只能表示有限个实数圬例如圬对于一个字长为五的计算机圬将无法区分地. 圱圱地地圱圬地. 圱圱地地圱地圱地圱圮既然计算机的实数是一个有限的集合圬那么初始数据和中间结果都有可能不在这个集合中圬于是必须用计算机实数去近似表示相应的实数圬从而引进舍入误差圮如果按照四舍五入的原则去近似表示圬则由计算机表示一个实数时引进的舍入误差不大于最后一位数的一半圬但很多计算机采用了只舍不入的方法圬则由计算机数表示一个实数时引进的舍入误差最大可为最后一个二进制位圮

10 圱地坃坈坁坐坔坅坒圱圮计算机数和误差 1.2 舍入误差对计算的影响 在数值计算过程中圬一般每一步都有舍入误差圬舍入误差的积累有可能使整个计算失败圮 例如圬考虑下述问题圬计算积分 利用分部积分法得到或 1 E n 1 0 x n e x 1 dx, n 圱, 圲,..., 圲地. x n e x 1 dx x n e x n x n 1 e x 1 dx, E n 圱 n E n 1 容易算出 E 1 圱 /e 圮利用上述递推公式圬在坉坂坍圭坐坃微机上用双精度计算结果如下场 E 1 地. 圳圶圷圸圷圹圴圴圱圱圷圱圴圴 E 3 地. 圲地圷圲圷圶圶圴圷地圲圸圶圵 E 5 地. 圱圴圵圵圳圲圹圴地圵圷圳地圸 E 7 地. 圱圱圲圳圸圳圵地圴地圶圹圳圶 E 9 地. 地圹圱圶圱圲圲圹圲圹圹圴圱圷圱 E 11 地. 地圷圷圳圵圲圲圲圹圳圵圸圷圸地 E 13 地. 地圶圶圹圴圷圷圷圹圹圶圹圷圲圳 E 15 地. 地圵圹地圳圳圷圹圳圶圴圱圹地圲 E 17 地. 地圵圷圱圹圱圸圷地圵圹圷圳地圸 E 19 圱. 圵圵圹圶圱圹圷圴圴圲圷圹圲 E 2 地. 圲圶圴圲圴圱圱圱圷圶圵圷圱圲 E 4 地. 圱圷地圸圹圳圴圱圱圸圸圵圳圸 E 6 地. 圱圲圶圸地圲圳圵圶圵圶圱圵圲 E 8 地. 圱地地圹圳圱圹圶圷圴圴圵地圹 E 10 地. 地圸圳圸圷圷地圷地地圵圸圲圹圳 E 12 地. 地圷圱圷圷圳圲圴圷圶圹圴圶圳圷 E 14 地. 地圶圲圷圳圱地圸地圴圲圳圸圷圳 E 16 地. 地圵圵圴圵圹圳地圱圷圲圹圵圷地 E 18. 地圲圹圴圵圳圶圷地圷圵圱圵圳圶 E 20 圳地. 圱圹圲圳圹圴圸圸圵圵圸圴 被积函数在积分区间坛地, 圱坝的内部总是正的圬因此积分值不可能为负圬更进一步圬由于 E n 1 0 x n e x 1 dx 1 0 x n dx 圱 n 圫圱 所以圬必有 E 20 圱 / 圲圱圮显然圬上述计算结果是不正确的圬而问题就出在舍入误差上面圮我 们对此作一个简单的分析圮双精度数可表示圱圴位圱地进制数圬所以 E 1 的误差约为 ɛ 圱地 15 在也可能大于此数圬这里仅作一数量级的估计圩圮在递推过程中圬这个误差将被传播圬放大圬 最后导致计算结果的错误圮误差的传播规律为 E 2 圱 圲在 E 1 圫 ɛ 圩圱 圲 E 1 圲 ɛ E 3 圱 圳在圱 圲 E 1 圲 ɛ 圩圱 圳在圱 圲 E 1 圩圫圳圡 ɛ

11 圱圮圳圮减法的计算圱圱 计算到 E 15 时的误差为圱圵圡 ɛ 圱地 3 圬而计算到 E 18 时的误差已为圱圸圡 ɛ 圶. 地圬已经远大 于 E 18 的上限圱 / 圱圹圮 现在我们考虑把递推关系改写成 E n 1 圱 E n, n, 圴, 圳, 圲. n 利用关系圱 E n n 圫圱取 E 30 地. 地开始计算圬此时误差为 ɛ 圱 / 圳圱 地. 地圳圲圬作一次递推后得到 E 29 圬其误差减小 到圱圱 地. 地地圱圱圳地圳圱在计算 E 20 时圬误差已减小到圲地圡 16 圳 圱地圳圱圡准确到圱圴位有效数字圬初始误差的影响已完全消除了圮而每一步的舍入误差, 都在其后 的迭代中减小和消除 下表是计算结果圮 E 30 地. 地 E 29 地. 圳圳圳圳圳圳圳圳圳圳圳圳圳圳 E 28 地. 地圳圳圳圳圳圳圳圳圳圳圳圳圳圳 E 26 地. 地圳圵圷圵圸圳圷圷圴圲圵地圴圴 E 24 地. 地圳圸圵圱圶圵圵圱圳圴圹圸圸圵 E 22 地. 地圴圱圷圳圶圴圴圳地圲圷圸地圸 E 20 地. 地圴圵圵圴圴圸圸圴地圷圵圸圱圸 E 18 地. 地圵地圱圱圹圸圵圴圹圵圸地圹圴 E 16 地. 地圵圵圷圱圹圳圴圵圹圳圱圲圳圶 E 14 地. 地圶圲圷圳圲圱圶圳圹圴圱圳圸地 E 12 地. 地圷圱圷圷圳圲圵圳圶圴圸地圳地 E 10 地. 地圸圳圸圷圷地圷地圱地圳圳圹圴 E 8 地. 圱地地圹圳圱圹圶圷圴圴圵圵圹 E 6 地. 圱圲圶圸地圲圳圵圶圵圶圱圵圳 E 4 地. 圱圷地圸圹圳圴圱圱圸圸圵圳圸 E 27 地. 地圳圴圵圲圳圸地圹圵圲圳圸圱地 E 25 地. 地圳圷地圸圶圲圱圶圲圵圲圸圸圳 E 23 地. 地圴地地圶圱圸圱地圳圶地圴圲圱 E 21 地. 地圴圳圵圵圷圴圳圴圴地圷圸圲圷 E 19 地. 地圴圷圷圲圲圷圵圵圷圹圶圲地圹 E 17 地. 地圵圲圷圷圱圱圱圹圱圶圸圹圹圵 E 15 地. 地圵圹地圱圷圵圴地圸圷圹圲圹圸 E 13 地. 地圶圶圹圴圷圷地圲圵圷圵圶圱圶 E 11 地. 地圷圷圳圵圲圲圲圸圸圶圲圶圶圴 E 9 地. 地圹圱圶圱圲圲圹圲圹圸圹圶圶圱 E 7 地. 圱圱圲圳圸圳圵地圴地圶圹圳地 E 5 地. 圱圴圵圵圳圲圹圴地圵圷圳地圸 E 3 地. 圲地圷圲圷圶圶圴圷地圲圸圶圵 E 2 地. 圲圶圴圲圴圱圱圱圷圶圵圷圱圲 E 1 地. 圳圶圷圸圷圹圴圴圱圱圷圱圴圴由上例可以看出舍入误差对计算的影响圬如果我们在计算中不去分析算法圬即使编出了 正确的程序圬也可能得出错误的圬有时甚至是荒谬的结果圮 1.3 减法的计算 也许圬大家认为减法是一种非常简单的计算圬但是在用计算机做减法时圬如不小心圬就很

12 圱圲坃坈坁坐坔坅坒圱圮计算机数和误差 容易出错圬问题就出在舍入误差上面圮我们考虑下面的例子场对于大的 x 圬计算 x 圫圱 x 在圱圮圳圮圱圩 取 x 圱圶圲圸圳圵. 地圬在一个具有圶位十进制有效数字的计算机上圬计算结果为地. 地地圱圬而准确 到圶位十进制有效数字的精确结果为地. 圱圲圳圹地圷 圱地 2 圬计算中损失了圵位有效数字圮对于更 大的 x 圬则可能根本得不到正确的结果圮但是圬如果把计算公式改为 在 x 圫圱 x 圫圱圫 x x 圩 x 圫圱圫 x 圱 x 圫圱圫 x 在圱圮圳圮圲圩 就可以得到较好的结果圮又如对于大的 N 圬计算 N+1 N dx 圱圫 x 2 坡坲坣坴坡坮在 N 圫圱圩 坡坲坣坴坡坮在 N 圩 时圬可以使用等价的公式场 ( ) 圱坡坲坣坴坡坮在 N 圫圱圩 坡坲坣坴坡坮在 N 圩坡坲坣坴坡坮圱圫 N 在 N 圫圱圩 同理圬对于小的 ɛ ( 坳坩坮在 x 圫 ɛ 圩 坳坩坮在 x 圩圲坣坯坳 x 圫 ɛ ) 坳坩坮圲 ( ɛ ) 圲 上面分析的是一些简单的例子圬在实际问题的计算中圬常常会碰到一些十分复杂的问 题圬即使有经验的计算研究人员也经常在减法计算中犯错误圮 练习 : 分析下面的表达式, 给出对任何 x( 只要在所给函数定义域内 ) 都能得出正确结果的计算 公式. 坣坯坳 2 x 坳坩坮 2 x 坬坮在 x 圩 圱 e 2x e x 圱 v 2 圱 圱圫 x 2 圱 x2 坳坩坮在 x 圩 x 2 x 坴坧在 x 圩

13 Chapter 2 物理学中的常用数值方法 在这一章和随后的几章中圬我们将介绍几种常用算法圬这些内容应该在数学课中讨论圬但 目前的数学课程并不包含这些属于计算数学课程的内容圬所以我们先作一些介绍圮 2.1 数值积分方法 积分分为定积分和不定积分两种圬我们在这里只讨论定积分圬不定积分的计算可作为变 上限的定积分来计算圮考虑定积分 它代表 f 在 x 圩曲线下的面积圮 I b a f 在 x 圩 dx 在圲圮圱圮圱圩 梯形法和辛普生方法 在微积分中圬我们知道积分是由求和的极限定义的圬这就是第一种计算积分的方法坼梯 形法圮把区间坛 a, b 坝分成一些子区间坛 a x 0 < x 1 < x 2 < < x n b 坝圬则在圲圮圱圮圱圩的积分 I 可 近似为场 I 圱 n 1 在 x i+1 x i 圩在 f 在 x i 圩圫 f 在 x i+1 圩圩圲 i=0 在圲圮圱圮圲圩 在实际计算中圬常采用等间距分法圬即 x i+1 x i h b a n 为一常数圬上面公式可简化为场 I h n 1 在 f 在 x i 圩圫 f 在 x i+1 圩圩圲 i=0 为了节省工作量圬实际使用的计算公式为场 n 1 I h f 在 x i 圩圫 i=1 h 在 f 在 x 0 圩圫 f 在 x n 圩圩圲 在圲圮圱圮圳圩

14 圱圴坃坈坁坐坔坅坒圲圮物理学中的常用数值方法 当子区间越来越小圬 h 地时, 求和的结果将趋于积分值 由于我们事先并不知道积分 的结果, 所以我们也就不知道求和结果和积分结果之间的差别 为此, 我们需要估计 计算误差, 一个直观的估计是比较不同大小的子区间的结果之差 为了在每次子区间 变小时能够用到前次的计算结果, 我们可以从 n 圱开始计算圬每次把区间数增加一倍圬 则前一次计算的函数结果仍然有用圬通过比较连续两次的计算结果之差是否小于给定的 误差限来结束计算圮 梯形方法实际上把每一个子区间上的函数近似为线性函数 αx 圫 β 圬然后计算其积分 值圬再把结果加起来圮在作为一个练习, 请自行证明圩圬因此圬我们说梯形方法是具有线性 代数精度的方法圮下面考虑具有二次代数精度的方法圬坼辛普生方法圮先考虑把积分区 间坛 a, b 坝分为二个子区间圬分点为在 a, a 圫 h, a 圫圲 h b 圩圬 h b a 2 圬在区间上把被积函数用二次函数来近似圬即取 f 在 x 圩 α 在 x a 圩 2 圫 β 在 x a 圩圫 γ 圬其中 α, β, γ 可利用分点上的函数 值 f 在 a 圩, f 在 a 圫 h 圩, f 在 b 圩表示为场 α f 在 b 圩 圲 f 在 a 圫 h 圩圫 f 在 a 圩圲 h 2 β f 在 b 圩圫圴 f 在 a 圫 h 圩 圳 f 在 a 圩圲 h γ f 在 a 圩 简单地计算得到场 b a 在 α 在 x a 圩 2 圫 β 在 x a 圩圫 γ 圩 dx h 在 f 在 b 圩圫圴 f 在 a 圫 h 圩圫 f 在 a 圩圩圳 在圲圮圱圮圴圩 这就是著名的辛普生求积公式圮显然圬它具有二次代数精度圮 实际计算中圬把被积区间坛 a, b 坝先分为 N 个子区间圬然后在每个子区间中用辛普生求积 公式计算圬并把最后结果加起来圮其计算公式为场 b a f 在 x 圩 dx h 在 f 在 a 圩圫圴 f 在 a 圫 h 圩圫圲 f 在 a 圫圲 h 圩圫圴 f 在 a 圫圳 h 圩圫 圫 f 在 b 圩圩圳 在圲圮圱圮圵圩 这是一个在实际计算中常用的方法圬它对于一般的有限区间的积分都能给出较好的正确结果圮辛普生方法有各种不同的实现圬我们在这里给出简单的一种圬取每个子区间均相等圬从二个子区间开始圬每次子区间数目加倍圬这样可以在每次计算时使用前次的计算值圮这种实现的缺点是在整个积分区间均匀取点圬对于在部分区间变化剧烈而在其余地方变化平缓的被积函数计算效率不高圮这一实现可以充分利用它与梯形法的关系圮假 k k 定 n 圲圬用 S k 表示用圲个子区间梯形公式求得的结果圮则式在圲圮圱圮圵圩可写为场

15 圲圮圱圮数值积分方法圱圵 b a f 在 x 圩 dx 圱坛圴 h 在 f 在 a 圫 h 圩圫 f 在 a 圫圲 h 圩圫 圫 f 在 b h 圩圫圳 圱 f 在 a 圩圫圲圱 圱 f 在 b 圩圩圲 圱 圲 h 在 f 在 a 圫圲 h 圩圫 f 在 a 圫圴 h 圩圫 圫 f 在 b 圲 h 圩圫 f 在 a 圩圫 f 在 b 圩圩圲圲圫 h 在 f 在 a 圩圫 f 在 b 圩圩 圲 h 在 f 在 a 圩圫 f 在 b 圩圩圫 h 在 f 在 a 圩圫 f 在 b 圩圩坝圱在圴 S k S k 1 圩在圲圮圱圮圶圩圳计算定积分还有很多方法圬如龙贝格方法圬坎坥坷坴坯坮圭坃坯坴坥坳方法等圬我们将不在此讨论圬 有兴趣的同学可参看有关计算数学的书籍圮坑坵坡坤坰坡坣坫 ( 坮坥坴坬坩坢圮坯坲坧 ) 中包含了若干非常 稳定的求积分的程序, 建议在实际工作中使用 在坍坡坴坬坡坢中, 梯形法, 辛普森均作为标准函数提供, 实际应用时, 只要写一个 被积函数的坭程序, 调用积分程序即可 常用的几个求积分的函数是坱坵坡坤圬坱坵坡坤坬圬 坱坵坡坤坧坫等 其中坱坵坡坤使用自适应的辛普森算法, 坱坵坡坤坬应用龙贝格加速来改善计算效率 和提高精度 例如, 计算积分 1 先定义函数坦坵坮坣并存入坦坵坮坣圮坭 funtion y=fun(x) end y=sin(x.^2)./(1.+x.^4); 然后在坍坡坴坬坡坢的命令窗口运行 0 坳坩坮在 x 2 圩圱圫 x 4 dx >> format long >> quad(@fun,0,1) 结果为 : ans = 反常积分的计算 反常积分分为两类圬一类是积分区间有限圬在积分区间内被积函数有奇点圬另一类是积分区间为无限圮对于二者兼而有之的积分圬可以分为两个或多个积分来处理圬使每个积分为上述两类之一圮

16 圱圶坃坈坁坐坔坅坒圲圮物理学中的常用数值方法 积分区间内含有奇点的积分 对于积分区间内含有奇点的积分圬根据奇点的性质圬可采用下面的方法处理圮可去奇点圬设 x 0 为一可去奇点圬且知道在 x 0 点被积函数的极限圬则只要在计算函数 sin(x) 时圬在 x 0 的一个邻域内用极限取代函数值即可圮例如圬在 x 地有可去奇点圬但 x sin(x) 当 x 地时圬 圱 x2 圫 圬在计算函数时圬可使用下面的语句场 x 6 IF(ABS(X).LT.1.0E-4) THEN F=1.0-X*X/6.0 ELSE F=SIN(X)/X ENDIF 极限方法圬有时我们并不知道奇点处函数的极限表达式圬有时甚至不知道奇点是否可 去圬此时可用极限方法来逼近圮例如圬若已知 x 0 为奇点圬对于积分 b x 0 f 在 x 圩 dx 坬坩坭 r x 0 b r f 在 x 圩 dx 定义一个收敛于 x 0 的序列 b > r 1 > r 2 > 圬例如可取 r n x 0 圫圲 n 圬则上述积分可写为 一个求和场 b b r1 r2 r3 f 在 x 圩 dx f 在 x 圩 dx 圫 f 在 x 圩 dx 圫 f 在 x 圩 dx 圫 f 在 x 圩 dx x 0 r 1 r 2 r 3 r 4 等式中每一项都是正常积分圬当 r n r n+1 f 在 x 圩 dx ɛ 时圬终止计算圮如果上式不收敛圬则很可 能原积分不存在圮 消除奇点圬有些奇点可以通过变量替换或对被积函数进行变换而消去圮我们看几个例 子场 1 作变量替换 x t 2 圬 dx 圲 t dt 圬则有场 或把积分写为场 1 0 坣坯坳在 x 圩 x dx 坣坯坳在 x 圩 x dx 坣坯坳在 x 圩 x dx 圲 圱 1 dx 圫 x 0 第二项积分中的奇点已成为可去奇点圮 I 在 x 圩 1 0 坣坯坳在 t 2 圩 dt 坣坯坳在 x 圩 圱 1 dx 圲圫 x 0 x 0 e t 圱 t 在 t 圱的邻域圬被积函数与 e 1 / 在圱 t 圩很相像圬因此变为场 I 在 x 圩 x 0 e 1 圱 t dt 圫 x 0 e t e 1 圱 t 所有的奇异性都在第一项圬而这一项可解析处理圮 圱 dt 坬坮 圱 x 圫 e x 0 坣坯坳在 x 圩 圱 dx x e t e 1 dt 圱 t 除上面例子演示的方法外圬还可用分部积分法消除奇点圬这里不再一一举例了圮

17 圲圮圱圮数值积分方法圱圷 积分区间为无限的积分 下面我们讨论积分区间为无限的情形圮为具体起见圬只讨论形为 0 f 在 x 圩 dx 的积分圬其它形式的无穷区间的积分均可化为上面的形式圮 变量替换法场通过变量替换圬可把无穷区间的积分化为有限区间上的积分圬如 令 x 坬坮在 t 圩圬 dx dt t 圬则有场 0 f 在 x 圩 dx 1 0 f 在 坬坮在 t 圩 1 dt t 0 g 在 t 圩 dt t 在圲圮圱圮圷圩 g(t) 这里 g 在 t 圩 f 在 坬坮在 t 圩圩圬若在 t 地的邻域内有界圬则上式的积分成为一个正常积分圬否 t 则就是前面讨论过的有限区间上的反常积分圬可以用前面的方法来处理圮上述变换对于 当 x 时圬 f 在 x 圩以 e kx 的形式衰减的函数圬计算结果十分令人满意圮其它常用的变换还 有 x t 1 圬 dx dt 1 t (1 t) 2 以及 t 坴坡坮坨 x 等圮 极限法场由定义 0 f 在 x 圩 dx f ( ) t 圱圱 t 在圱 t 圩 dt 2 r f 在 x 圩 dx 坬坩坭 f 在 x 圩 dx r 0 在圲圮圱圮圸圩 定义一个趋于 的序列地 < r 1 < r 2 < 圬例如可取 r n 圲 n 圬则上述积分可写为一个求 和场 0 f 在 x 圩 dx r1 0 f 在 x 圩 dx 圫 r2 r 1 f 在 x 圩 dx 圫 r3 r 2 f 在 x 圩 dx 圫 在圲圮圱圮圹圩 式中每一项都是正常积分圬当 r n+1 r n f 在 x 圩 dx ɛ 时圬终止计算圮如果上式不收敛圬则很可 能原积分不存在圮如果能找到尾项 f 在 x 圩 dx 的一个较好的渐近表达式圬则上述积分只需 R 进行到某一足够大的 r n 处圬再加上尾项圮 还有一些其它的有效方法圬这里不再一一列举圬我们下面就要讲到的高斯积分方法提 供了另一种处理反常积分的工具圮 高斯积分方法 前面的积分公式在梯形公式及辛普生公式圩的代数精度都较低圬一个自然的问题是圬能否构造出高代数精度的公式圬回答是肯定的圬这就是高斯积分方法圮我们将不讨论高斯积分方法的理论基础圬而只是给出有关的计算公式圬对理论基础有兴趣的同学可以参看有关的计算数学的书籍圮

18 圱圸坃坈坁坐坔坅坒圲圮物理学中的常用数值方法 有限区间的高斯型积分公式场 高斯坼勒让德积分公式场 1 1 f 在 x 圩 dx n w i f 在 x i 圩 此处圬 x i 为 n 阶勒让德多项式的根圬通常称为积分节点圬 w i 称为积分权重圮 等权重切贝雪夫积分公式场 这里 x i 为积分节点圬权重为圲 /n 圮 第一类高斯坼切贝雪夫积分公式场 f 在 x 圩 圱 x 2 dx 第二类高斯坼切贝雪夫积分公式场 1 1 f 在 x 圩圱 x 2 dx π n 圫圱 f 在 x 圩 dx n i=1 π n i=1 圲 n n f i=1 n f 在 x i 圩 i=1 ( ( 坣坯坳 )) 在圲 i 圫圱圩 π 圲 n ( ) ( ( )) 2 iπ iπ 坳坩坮 f 坣坯坳 n 圫圱 n 圫圱 如果积分区间不是坛 圱, 圱坝圬而是坛 a, b 坝圬可利用变换 ( ) ( ) b a b 圫 a x t 圫圲圲 在圲圮圱圮圱地圩在圲圮圱圮圱圱圩在圲圮圱圮圱圲圩在圲圮圱圮圱圳圩在圲圮圱圮圱圴圩 将积分区间变为坛 圱, 圱坝圬然后利用上述公式进行计算圮在实际计算时圬常把积分区间分为几个小区间圬在每一个小区间内使用高斯积分公式圬并把结果加起来圮练习 : 试编写上述四种高斯积分公式的程序, 并利用所编写的程序计算积分 1 1 x 圫圱. 圵 dx, 1 1 圱圲圫 x dx, 1 0 dx x 2 0 e os2 (x) dx, ( 积分所需的节点及权重在后面的表中给出 ). 1 0 坳坩坮在 x 圩 x dx 无限区间的高斯型积分公式场 高斯坼拉盖尔积分公式场 0 e x f 在 x 圩 dx n w i f 在 x i 圩 i=1 在圲圮圱圮圱圵圩

19 圲圮圱圮数值积分方法圱圹 高斯坼厄米积分公式场 e x2 f 在 x 圩 dx n w i f 在 x i 圩 i=1 在圲圮圱圮圱圶圩 这里 x i 分别为 n 阶拉盖尔多项式和厄米多项式的根圮其数值和对应的权重的数值在附 录的表中给出圮 练习场 试编写上述二种高斯积分公式的程序, 并利用所编程序计算积分 0 e x 圱圫 x 4 dx, e x2 0 坣坯坳在 x 圩 dx 各阶高斯积分公式的取点并不重和, 当增加分点时, 所有的函数值都要重 算 圱圹圶圵年坋坲坯坮坲坯坤提出了结合高斯方法和加速收敛技巧的算法 这一算法也可以 处理有一定奇异性的被积函数, 但需要把奇点放在积分限上 在坍坡坴坬坡坢中已实现这一 算法, 其调用命令是坱坵坡坤坧坫圮 例如, 计算积分 定义函数坦坵坮坳坩坮圬并存入坦坵坮坳坩坮圮坭 0 坳坩坮 x x e x dx funtion y=funsin(x) end y=sin(x).*exp(-x)./x; 在坍坡坴坬坡坢命令窗口键入如下命令 format long quadgk(@funsin,0,inf) ans = 这一积分的精确结果是 π 4 地. 圷圸圵圳圹圸圱圶圳圳圹圷圴圴圸

20 圲地坃坈坁坐坔坅坒圲圮物理学中的常用数值方法 哥西主值的计算 在物理计算中圬哥西主值的计算占了很重要的位置圬如物质的介电函数的实部和虚部之 间满足著名的坋坲坡坭坥坲坳圭坋坲坯坮坩坧关系场 ɛ 在 ω 圩 ɛ 1 在 ω 圩圫 iɛ 2 在 ω 圩 圲 ɛ 1 在 ω 圩圱圫 π P ɛ 2 在 ω 圩 圲 ω π P 0 0 sɛ 2 在 s 圩 s 2 ω ds 2 ɛ 1 在 s 圩 s 2 ω ds 2 在圲圮圱圮圱圷圩 这个关系是一个严格的关系圬它来源于因果关系圬是非常普遍的圮在物理上圬 ɛ 2 在 ω 圩的测量 和计算都比较容易圬因此圬一般不直接测量或计算 ɛ 1 在 ω 圩圬而是利用上述关系由 ɛ 2 在 ω 圩的数 据进行计算圮 上述坋坲坡坭坥坲坳圭坋坲坯坮坩坧关系中出现的坐就代表哥西主值积分圬设 a < < b 圬 f 在 x 圩在 x 的邻域内无界圬则哥西主值积分定义为场 P b a f 在 x 圩 dx 坬坩坭 r 0 [ r a f 在 x 圩 dx 圫 不失一般性圬我们取 地并将积分取成 a a f 在 x 圩 dx 的形式圮令 则 g 在 x 圩 圱坛 f 在 x 圩 f 在 x 圩坝, 圲 由于 g 在 x 圩是奇函数而 h 在 x 圩是偶函数圬于是场 因此圬 P b a r a r 圲 a a r f 在 x 圩 dx 圫 g 在 x 圩 dx 圫 h 在 x 圩 dx f 在 x 圩 dx 圲坬坩坭 r 0 a r h 在 x 圩 f 在 x 圩 g 在 x 圩圫 h 在 x 圩 a r a r f 在 x 圩 dx g 在 x 圩 dx 圫 h 在 x 圩 dx 圲 a 0 r a b h 在 x 圩 dx +r ] f 在 x 圩 dx 圱坛 f 在 x 圩圫 f 在 x 圩坝圲 h 在 x 圩 dx 圫 a 0 a r h 在 x 圩 dx 坛 f 在 x 圩圫 f 在 x 圩坝 dx 在圲圮圱圮圱圸圩 这样我们就把哥西主值积分的计算变为在 x 地处有奇点的一个反常积分的计算圮 h 在 x 圩也 可能在 x 地处无奇点圬此时只要计算一个常规积分圮例如前面的坋坲坡坭坥坲坳圭坋坲坯坮坩坧关系中 积分的计算可化为场 圲 π P 0 sɛ 2 在 s 圩 s 2 ω 2 ds 圱 圲 π P ɛ 2 在 s 圩 s ω ds 圱 π 0 ɛ 2 在 ω 圫 x 圩 ɛ 2 在 ω x 圩 dx x

21 圲圮圱圮数值积分方法圲圱 圲 ω π P 0 ɛ 1 在 s 圩 s 2 ω 2 ds 圱 圲 π P ɛ 1 在 s 圩 s ω ds 圱 π 0 在得到上述关系时圬我们利用了 ɛ 1 为偶函数而 ɛ 2 为奇函数的事实圮 ɛ 1 在 ω 圫 x 圩 ɛ 1 在 ω x 圩 dx x 当哥西积分区间内的无界邻域不止一个时圬可以把积分区间分为几段圬使得每一段中 只包含一个无界邻域圬然后在每一段用上述方法计算圮 另一种计算主值积分的方法是直接从定义出发圬设 f 在 x 圩在积分区间内无奇点圬则积分 P b a f 在 x 圩 x x 0 dx x0 f 在 x 圩 dx 圫 a x x 0 I 圫坬坩坭 I ɛ, ɛ 0 b x 0 + f 在 x 圩 x x 0 dx 圫 P x0 + x 0 I 为第一行右方第一和第二项之和圬是正常积分圬可用前述任一方法计算圮 I ɛ, x0 ɛ x 0 把 f 在 x 圩在 x 0 点展开圬并逐项积分圬得到 取极限 ɛ 地圬得到 f 在 x 圩 x x 0 dx 圫 x0 + x 0 +ɛ f 在 x 圩 x x 0 dx I ɛ, 圲 f 在 x 0 圩在圁 ɛ 圩圫 f (3) 在 x 0 圩在圁 3 ɛ 3 圩圫 I 0,δ 圲 这样圬原主值积分为 I 0, 圫 I 圮 练习场 f 在 x 圩 x x 0 dx 在圲圮圱圮圱圹圩 f (2n+1) 2n+1 圁在 x 0 圩. 在圲圮圱圮圲地圩在圲 n 圫圱圩圡在圲 n 圫圱圩 n=0 试用上述方法计算主值积分 1 e x P, 0 x 地. 圵 对于不同的圁, 在 (2.1.20) 中取两项, 比较计算结果与精确结果. ( 到 16 位的精 确结果为 地. 圶圱圵地圱圸圱圴圶圲圸地圵圶圷圴 ) 另外一种常用的计算主值积分的方法是基于如下公式 坬坩坭 η 0 + 圱 x iη P 圱圫 iπδ 在 x 圩 x 在圲圮圱圮圲圱圩 式中 η 地 + 表示 η 从大于地的方向趋于地圮这一公式当然是在积分的意义下理解 利用上 述结果, 就可以得到计算主值积分的一个方法 : 例如要计算 P b a F 在 x 圩 dx b 取一个小量 η, 计算 F 在 x iη 圩 dx 圬其实部给出 P b F 在 x 圩 dx 的依赖于 η 的近似值, 通过 a a 减小 η, 可以得到精度要求下收敛的结果 例如, 对于前面的练习题, 分别取不同的 η 圬计算结果如下

22 圲圲坃坈坁坐坔坅坒圲圮物理学中的常用数值方法 η 地. 地圱地. 地地圱地. 地地地圱地. 地地地地圱 积分值 地. 圵圹圶圱圹圶圱圶圶地圲圶地圳圲圫圱. 圸圸圴圱圷圴地地圶圱圵圶圵圶圸 i 地. 圶圱圳圱圱圴圹圹圸圹圴圱圲圷地圫圱. 圹地圳圳圵地圵圷圵圶圸圵圶圹圸 i 地. 圶圱圴圸圲圷圶圲圲圳地地圵圴圹圫圱. 圹地圵圲圶地圱圸圱圱圳圹地圲圶 i 地. 圶圱圴圹圹圹地圹圱圷圹圱圲圲圳圫圱. 圹地圵圴圵圱地圵圷圲圳圳圵圴圱 i 急速振荡函数的积分 在各种变换特别是付里叶变换中圬经常会遇到急速振荡函数的积分圬例如圬付里叶变换 b a f 在 x 圩坣坯坳在 nx 圩 dx, b a f 在 x 圩坳坩坮在 nx 圩 dx 又如付里叶坼贝塞尔变换圬 1 这里地 < γ 1 < γ 2 < 是贝塞尔函数的根圮一般我们可把这类函数写为场 I 在 t 圩 b a 0 f 在 x 圩 xj n 在 γ m x 圩 dx f 在 x 圩 K 在 x, t 圩 dx, < a < b < 在圲圮圱圮圲圲圩 式中 K 在 x, t 圩为一振荡核而 f 在 x 圩为非振荡部分圮由于振荡型函数在积分区间内多次取几乎相等的正值和负值圬如果用通常的方法计算圬则由于正负数相加圬造成有效数字的损失圬几乎得不到正确的结果圮下面先给出最常遇到的振荡函数的积分坼付里叶系数的一种计算方法圬然后再讨论较为一般的方法圮付里叶系数的计算 : 付里叶系数的复数形式为场 b a f 在 x 圩 e isx dx 对上式分部积分圬假定 f 在 x 圩存在直到 n 阶的导数圬可以得到场 在圲圮圱圮圲圳圩 b a n 1 f 在 x 圩 e isx dx e isa k=0 ( ) k+1 i n 1 f (k) 在 a 圩 e isb s k=0 ( ) k+1 i f (k) 在 b 圩圫 R n s 在圲圮圱圮圲圴圩 当 n 趋于无穷大时圬余项 R n 趋于地圬而求和项就给出了积分的近似值圮 练习 : 写出 (2.1.24) 的实部和虚部, 即 式. b a f 在 x 圩坣坯坳 sxdx 和 b a f 在 x 圩坳坩坮 sxdx 的积分公

23 圲圮圱圮数值积分方法圲圳 上述公式需要计算 f 在 x 圩的直到 n 圱阶的导数圬这对于较为复杂的函数来说圬曾是一个 很重的负担圮不过目前计算机代数的发展圬计算导数已经没有多少困难圮为了看出这一方 法的效力圬我们来看一个例子圮计算 I 2π 0 x 坣坯坳在 x 圩坳坩坮圳地 xdx 用公式在圲圮圱圮圲圴圩的虚部圬取 n 圳圬由 f 在 x 圩 x 坣坯坳在 x 圩圬 a 地圬 b 圲 π 圬 f 在地圩地圬 f 在圲 π 圩圲 π 圬 f 在地圩地圬 f 在圲 π 圩 圲 π 圬代入得到场 I 2π 0 而该积分的精确值为 地. 圲地圹圶圷圲圴圷圮 ( 圱 x 坣坯坳在 x 圩坳坩坮圳地 xdx 圲 π 圫圳地 圱圳地 ) 3 在 圲 π 圩 地. 圲地圹圶圷圲圲圲 基于 Euler-Malaurin 公式的梯形法 : 可以证明 ( 见王竹溪, 郭敦仁 特殊函数概 论 第一章 ), 函数 f 在 x 圩的积分可以写成 a+mh a 其中 [ f 在 a 圉 f 在 x 圩 dx h 圲 n 在 圱圩 k B k h 2k 圫在圲 k 圩圡 k=1 R n θ 圫 f 在 a 圫 h 圩圫 圫 f 在 a 圫在 m 圱圩 h 圩圫 [ ] f 2k 1 在 a 圫 mh 圩 F 2k 1 在 a 圩圫 R n ] f 在 a 圫 mh 圩圲 在 圱圩 n+1 B n+1 h [ ] (2n+1) f (2n+1) 在 a 圫 mh 圩 f (2n+1) 在 a 圩在圲 n 圫圲圩圡 在圲圮圱圮圲圵圩 这里地 θ 圱,B n 是坂坥坲坮坯坵坬坬坩数 如果 f 在 x 圩是周期函数,f 在 a 圫 mh 圩 f 在 a 圩圬且 f 在 x 圩存 在直到在圲 n 圫圱圩阶的导数, 则显然有 f (2k+1) 在 a 圩 f (2k+1) 在 a 圫 mh 圩圬于是, 利用梯形法可 以得到精确的结果 作为例子, 取不同的点, 试用梯形法计算如下积分, 并于精确结 果比较 1 圱 1 0 圱圫坳坩坮圱地 πx dx, 精确值圱. 圱圵圴圷地地圵圴 2 坂坥坳坳坥坬函数的积分表示之一是 J n 在 t 圩 计算 n 圲圬 t 圱圬圳圬圷圬圱地时坂坥坳坳坥坬函数的值 圱 π 坣坯坳圈 t 坳坩坮 x nx 圩 dx π 0 由坅坵坬坥坲圭坍坡坣坬坡坵坲坩坮公式可知, 如果一个函数满足 f 在 a 圩 f 在 a 圫 mh 圩圬 f (3) 在 a 圩 f (3) 在 a 圫 mh 圩圬 圬则梯形公式也可以得到很精确的结果 例如函数 f 在 x 圩坥坸坰在 x 2 在圱 x 2 圩圩圬在地和圱满足上述关系, 从而, 可以用梯形公式计算 1 0 坥坸坰在 x 2 在圱 x 2 圩圩 dx

24 圲圴坃坈坁坐坔坅坒圲圮物理学中的常用数值方法这个积分的精确结果是圱. 地圳圴圱圴圱地圵圬请分别用梯形法和辛普森法计算并比较其结果 子区间法圬在一般情况下圬如果 K 在 x, t 圩地的根容易求出圬且记为 a x 1 < x 2 < < x n b 圬则可用普通的积分规则计算每一个子积分 xi+1 x i 圬并把结果加起来圮用这种方法 得到的一般是一个交错级数圬且相邻项的绝对值比较接近圬其求和可用坅坵坬坥坲变换计算圮在坅坵坬坥坲变换将在后面讨论圩在讲了样条插值之后圬我们再介绍一种基于样条插值的计算积分在圲圮圱圮圲圳圩的方法圮更多的方法可以参看 ( 坐圮坊圮坄坡坶坩坳坡坮坤坐圮坒坡坢坩坮坯坷坩坴坺, 坍坥坴坨坯坤坳坯坦坎坵坭坥坲坩坣坡坬坉坮坴坥圭坧坲坡坴坩坯坮圬坁坣坡坤坥坭坩坣坐坲坥坳坳圬圱圹圸圴圩 高维积分的计算 高维积分的计算原则上可以通过化成累次积分来进行圬例如圬对于积分 I g 在 x, y, z 圩 dxdydz Ω 在圲圮圱圮圲圶圩 其中圊为积分区域圬总可以化为 I b dx d(x) dy f(x,y) a (x) e(x,y) g 在 x, y, z 圩 在圲圮圱圮圲圷圩 这相当于依次计算三个一重积分圬每一积分以其前一积分作为被积函数圮稍加考虑就可发现圬这一过程的计算量随积分重数的增加而变得无法完成圮如每重积分取圱地个点在这是 n 一个非常小的数字圩圬则 n 重积分就要计算圱地个函数值圮因此圬这一方法只能用于积分重数较低的积分圮高维积分还可以用高斯方法计算圬这只要把每重积分用高斯积分公式代入即可圮对于 n > 圳的高维积分圬计算量总是相当大的圬如果你对计算精度要求不高圬那么可以利用后面将要讲到的坍坯坮坴坥坃坡坲坬坯方法圮另外圬如果可能圬请尽可能地解析积出内层积分圬以减小计算量圮 固体热容量的计算 在德拜近似下圬固体的热容量由下式给出 ( ) 3 T Θd /T C V 圹 Nk 圂 D 0 ξ 4 e ξ 在 e ξ 圱圩 2 dξ 式中圂 D 为德拜温度圮为了求得热容量圬需要计算上式中的积分圬这一积分只有 当圂 d/t 圱或圂 d/t 圱时才能解析积出圬一般情况下只能进行数字计算圮考虑 f 在 x 圩 x 0 ξ 4 e ξ 在 e ξ 圱圩 2 dξ

25 圲圮圲圮方程的求根和最优化方法圲圵 利用我们学到的方法圬可以对不同的 x 求出其值圬并做出 f 在 x 圩的表或图形圮为了计算积 分圬注意到 ξ 地是被积函数的可去奇点圬因此圬应在该点把被积函数展开为 ξ 的级数圻 当 ξ 圬被积函数指数减小圮已知 f 在 圩 0 ξ 4 e ξ 在 e ξ 圱圩 2 dξ 圴 π 4 圱圵 我们可以根据 x 的值的大小分别用两种方法计算圬取一值 x 圬当 x < x 时圬用原式计算圬 当 x > x 时圬可以计算 f 在 圩 f 1 在 x 圩圬其中 f 1 在 x 圩 x ξ 4 e ξ 在 e ξ 圱圩 2 dξ x 可以选择的大一点圬这样 f 1 在 x 圩的数值比较小圬计算精度要求可以低一些圮 我们选择 x 圱地圬 x < x 时圬用辛普生方法计算 f 在 x 圩圬被积函数为 { ξ 2 ξ4 ξ 圫 地 ξ < 地. 地圵 ξ 4 e ξ (e x 1) 2 地. 地圵 ξ x 当 x > x 时圬计算 f 1 在 x 圩圬对积分作变量变换 ξ η 圫 x 则 f 1 在 x 圩 e x 0 在 η 圫 x 圩 4 e η 在圱 e (η+x) 圩 2 dη 这一积分可用高斯坼拉盖尔方法计算圮作为练习圬请同学完成余下的计算并作出 f 在 x 圩 x 的图形 2.2 方程的求根和最优化方法 方程的求根问题是物理计算中最常遇到的问题之一圬对于低次代数方程圬已经找到了解的公式圬但除了二次方程之外圬这些公式并无多少实用价值圬而对于即使如坣坯坳在 x 圩 x 地这样简单的方程圬也无法找到解析形式的解圮数值计算几乎是求解任何有实际意义的问题的唯一方法圮在这一节中圬我们将讨论几种求解 f 在 x 圩地的根的数字方法圮最优化方法除了本身的重要性之外圬也是一种求解多元方程的有效方法圬因此圬最优化方法的介绍也是本节的一个重要内容圮 二分法 二分法是求解方程的最安全圬也是概念上最简单的方法圬它可以做为后面将要讲到的方法失败后的最后保险圮假定我们知道在区间 a < x < b 之内只有 f 在 x 圩地的一个根圬则显然 f 在 a 圩与 f 在 b 圩异号圬把区间分半圬即令 在 b 圫 a 圩 / 圲圬计算 f 在 圩的值圬若 f 在 a 圩与 f 在 圩同号圬则根在坛, b 坝之间圬否则圬根在坛 a, 坝之间圬对根所在的区间重复上述过程直到达到所需精度圮

26 圲圶坃坈坁坐坔坅坒圲圮物理学中的常用数值方法 牛顿法和割线法 在所有一元函数求根方法中圬牛顿法也许是最值得推荐的一个圮其想法是这样的圬假 定 f 在 x 圩地的根为 x 圬而我们有一个对 x 的猜测值 x n 圬 f 在 x 圩可在 x n 附近展开为泰勒级 数场 f 在 x 圩 f 在 x n 圩圫 f 在 x n 圩在 x x n 圩圫 略去高次项圬令上式为地圬可解出 x 圬将此值作为 x 的一个更好的近似圬记为 x n+1 圬则 x n+1 x n f 在 x n 圩. 在圲圮圲圮圲圸圩 f 在 x n 圩 这样圬给定一个 x 的初始值 x 0 圬利用下面的迭代公式反复迭代场 x n+1 x n f 在 x n 圩 f 在 x n 圩 在圲圮圲圮圲圹圩 如果迭代收敛圬则必收敛于 f 在 x 圩地的一个根 x 圮这一迭代公式在圲圮圲圮圲圹圩就是著名的牛顿 迭代公式圮 可以证明圬如果迭代到第 n 步时的误差为 ɛ n x n x 圬则第 n 圫圱步的误差为 ɛ n+1 ɛ 2 f (x ) n 2f (x ) 圬也就是说圬牛顿法是平方收敛的圮证明如下场由迭代公式在圲圮圲圮圲圹圩圬 x n+1 x n f 在 x n 圩 f 在 x n 圩 x 圫 ɛ n ɛ nf 在 x 1 圩圫 2 ɛ2 nf 在 x 圩 f 在 x 圩圫 ɛ n f 在 x 圩 x 圫 ɛ n x 圫 ( ɛ n 圫 圱 f 在 x 圩 ɛ 2 圲 f 在 x n 圩 圱圲 f 在 x 圩 f 在 x 圩 ɛ 2 n f 在 x 圩 ɛ 2 f 在 x n 圩 二分法是线性收敛的, 请证明之圮但是圬另一方面圬对于任一初始值圬牛顿法并不总是收 敛圬事实上圬牛顿法的收敛范围是较小的圬因此圬在实际使用时圬必须找一个较好的 x 0 圬 使 x 0 尽可能地靠近 x 圬这可以把二分法与牛顿法结合起来以完成圮 试计算圱地. 令 f 在 x 圩 x 2 圱地, f 在 x 圩地的根就是所求结果 取初始值 x 0 圳, f 在 x 圩圲 x. 迭代结果如下 : ) 圳. 圳. 圱圶圶圶圶圶圶圶圶圶圶圶圶圶圷圳. 圱圶圲圲圸地圷地圱圷圵圴圳圸圶圳. 圱圶圲圲圷圷圶圶地圱圶圹圸圴圲圳. 圱圶圲圲圷圷圶圶地圱圶圸圳圸地圳. 圱圶圲圲圷圷圶圶地圱圶圸圳圷圹

27 圲圮圲圮方程的求根和最优化方法圲圷牛顿法的一个重要的缺点是需要计算函数的导数圬对于一些复杂的函数来说圬这是一件十分麻烦的工作圮为了即保持较高的收敛速度圬又避免麻烦的导数计算圬可以对牛顿法稍作变形圬用差分代替导数圬从而得到下面的割线法圮 x n+1 x n x n x n 1 f 在 x n 圩 f 在 x n 圩 f 在 x n 1 圩 在圲圮圲圮圳地圩 这一迭代方式需要二个初始值 x 0 圬 x 1 圬可以取为所猜测的解附近的二个点圮割线法的收敛 速度比牛顿法稍慢圬但省去了导数的计算圬对于求导数比较困难的函数特别有用圮作为一 个编程练习圬请写出这一算法的程序并作为你自己的库程序圮 练习 : 证明割线法的收敛阶为 ɛ n+1 ɛ n 另一个不需要计算导数的方法是所谓的坓坴坥坷坥坮坳坯坮方法圬算法如下场 x n+1 x n f 在 x n 圩 2 f 在 x n 圩 f 在 x n f 在 x n 圩圩 在圲圮圲圮圳圱圩 对于单零点圬这一算法的收敛阶为圲圬在请证明圩圬而且不需要计算导数圬是一个较好的算法圮 试计算圱地. 令 f 在 x 圩 x 2 圱地, f 在 x 圩地的根就是所求结果 取初始值 x 0 圳. 用 Stewenson 法 迭代结果如下 : 一维最优化, 黄金分割法 圳. 圳. 圱圴圲圸圵圷圱圴圲圸圵圷圱圴圳 圳. 圱圶圱圹圶圵圴圲圳圱圱圱圹圲地 圳. 圱圶圲圲圷圷圵圷圸圱圱圳圵圶圶 圳. 圱圶圲圲圷圷圶圶地圱圶圸圳圷圴 圳. 圱圶圲圲圷圷圶圶地圱圶圸圳圸地 这一节开始我们讨论最优化方法圬或计算一目标函数的极值的方法圮最优化方法在 物理学中有十分广泛的应用圬函数的求根问题也可化为一个最优化问题圬例如圬求函 数 f 在 x 圩地的问题就与求目标函数为 F 在 x 圩 f 在 x 圩 2 的极小值的问题等价圬而方程组 f i 在 x 1, x 2, x 3,, x n 圩地, i 圱, 圲, 圳,, n 的求解问题则与 F 在 x 1, x 2,, x n 圩 n i=1 2 w i f i 在 x 1, x 2,, x n 圩

28 圲圸坃坈坁坐坔坅坒圲圮物理学中的常用数值方法 的最优化问题等价圮这里 w i 圬 i 圱圬圲圬 圬 n 为一组权重参数圬适当的选则可使计算更为 容易进行圬但对根的位置没有影响圮 这一节我们先考虑一个自变量的目标函数的最优化问题圮假定我们有一目标函 数 F 在 x 圩圬我们的任务是寻找一个 x x 圬使得目标函数 F 在 x 圩最小圮这一工作可以分两步 来完成圬首先圬我们要找出目标函数的极小点所在的大致位置圬或者更确切一些圬我们要 找三个点 a < b < 圬使得目标函数 F 在 b 圩 < F 在 a 圩圬 F 在 b 圩 < F 在 圩圮第二步则在区间坛 a, 坝内寻 找 F 在 x 圩的最小值圮 第一步原则上是十分简单的圬假定有两点 a 圬 b 圬我们可以比较 F 在 a 圩和 F 在 b 圩以确定一个 函数减小的方向圬不失一般性圬假定 F 在 b 圩 < F 在 a 圩圬然后沿减小方向前进一步圬到达新的一 点 圬一般取 b 圫圱. 圶圱圸在 b a 圩圬如果 F 在 圩 < F 在 b 圩圬则令 a b 圬 b 圬继续前进圬寻找新 的 圻如果 F 在 圩 > F 在 b 圩圬则我们的目的已经达到圮 下面考虑第二步的计算圬即从三个点 a < b < 圬 F 在 b 圩 < F 在 a 圩圬 F 在 b 圩 < F 在 圩出发圬寻 找 F 在 x 圩的极小点 x min 圮这里介绍黄金分割法圮寻找极小点的过程圬实际上就是用一种方 法不断缩小区间坛 a, 坝圬并保证极小点一直位于区间之内圮这一过程可以这样来完成圬假定 已经作了 n 步圬得到三个点 a, b, 和对应的目标函数值 F 在 a 圩 > F 在 b 圩圬 F 在 b 圩 < F 在 圩圬不失一 般性圬假定 b > b a 圬在区间坛 b, 坝内选一点 x 圬如果 F 在 x 圩 < F 在 b 圩圬则作代换 a b 圬 b x 圬 新的缩小了的区间为坛 a, 坝在即原坛 b, 坝圬丢掉的区间为原坛 a, b 坝圩圻如果 F 在 x 圩 > F 在 b 圩圬则作代 换 x 圬新的缩小了的区间为坛 a, 坝在即原坛 a, x 坝圬丢掉的区间为原坛 x, 坝圩圮问题在于如何选 择 x 圬一个自然的想法是圬不论 F 在 x 圩 < F 在 b 圩或 F 在 x 圩 > F 在 b 圩圬丢掉的区间都相同圬即选择 x 圬 使得 b a x w 在 a 圩 在圲圮圲圮圳圲圩 这里 w 是一个比例因子圬我们现在来确定它圮由于要求每一步都满足上述关系圬显然有 x b w 在 b 圩 w 在在 a 圩 在 b a 圩圩 w 在 a 圩 w 2 在 a 圩在 w w 2 圩在 a 圩在圲圮圲圮圳圳圩 另一方面圬由在圲圮圲圮圳圲圩可得 即 在 b a 圩圫在 x 圩在 a 圩 在 x b 圩圲 w 在 a 圩 x b 在圱 圲 w 圩在 a 圩 在圲圮圲圮圳圴圩 在圲圮圲圮圳圵圩 由在圲圮圲圮圳圳圩和在圲圮圲圮圳圵圩圬我们得到场 w 2 圳 w 圫圱地 在圲圮圲圮圳圶圩 由此解得圬圳 圵 w 地. 圳圸圱圹圷, 圱 w 地. 圶圱圸地圳圲圱 w 通常被称为黄金分割数或黄金数圮一般在开始计算时圬 a 圬 b 圬 三者并不满足黄金分割关系圬但几次迭代后圬这一关系便可满足圬并在其后一直保持圮

29 圲圮圲圮方程的求根和最优化方法圲圹 高维最优化, 共扼梯度方法 这一节给出多维问题的最优化方法圮对于一个含有 n 个自变量的目标函数 F 在 x 1, x 2, x n 圩圬 我们的目的是要求出一组在 x 1, x 2,, x n 圩使目标函数达到极小圮为了讨论问题的方便圬 下面用 x 代表 n 维空间的一个点圬最优化问题也可以表述为圬寻找 n 维空间的一点 x 使 得目标函数 F 在 x 圩最小圮 高维最优化问题有很多求解方法圬其基本思路是相同的圬这就是从一点 x 0 出发圬按照 某种规定的方向 p 0 圬用一维最优化的方法寻找函数 F 在 x 圩的极小点 x 1 圬继续这一过程直到 达到最优点圮若第 i 步到达 x i 圬设选定方向 p i 圬则新的 x i+1 可如下求得圮实际上圬我们只要 找一参数 t i 圬使得 F 在 x i 圫 t i p i 圩坭坩坮 F 在 x i 圫 tp i 圩 t 在圲圮圲圮圳圷圩 则 x i+1 x i 圫 t i p i 便为所求圮因此圬高维最优化的的不同方法实际上是选择不同的一维寻 找方向圮一个直观的推断就是在每一点选择该点目标函数下降最快的方向 p i F 在 x i 圩 在圲圮圲圮圳圸圩 作为寻找方向圬这一方法称为最陡下降法圮直观的想象可能会认为最陡下降方向是一种理想的寻找方向圬但事实并非如此圬最陡下降方向只表示在 x i 附近的下降性质圬而对整个最优化过程来说圬我们寻找的是全局的最优方向圮作为一个例子圬考虑如下函数的最优化场 F x 2 圫圱地 y 2 在圲圮圲圮圳圹圩显然圬这一函数的极小点是在地, 地圩圮函数在圲圮圲圮圳圹圩的负梯度为 T g 在圲 x, 圲地 y 圩 在圲圮圲圮圴地圩 如果从点在圲, 圱圩出发圬利用最陡下降法计算所得逐 x 如下场 在圲, 圱圩圬在圱. 圷圹圲圸圲圹, 地. 地圳圵圸圵圷圩圬在地. 圴圶圱地圱圳, 地. 圲圳地圵地圷圩圬 在地. 圴圱圳圲圵圹, 地. 地地圸圲圶圵圩圬在地. 圱地圶圲圶圷, 地. 地圵圳圱圳圳圩圬在地. 地圹圵圲圵圹, 地. 地地圱圹地圵圩圬 在地. 地圲圴圴圹圵, 地. 地圱圲圲圴圸圩圬在地. 地圲圱圹圵圸, 地. 地地地圴圳圹圩圬在地. 地地圵圶圴圶, 地. 地地圲圸圲圳圩圬 在地. 地地圵地圶圱, 地. 地地地圱地圱圩圬在地. 地地圱圳地圲, 地. 地地地圶圵圱圩圬在地. 地地圱圱圶圷, 地. 地地地地圲圳圩圬 在地. 地地地圳地地, 地. 地地地圱圵地圩圬在地. 地地地圲圶圹, 地. 地地地地地圵圩圬在地. 地地地地圶圹, 地. 地地地地圳圵圩圬 在地. 地地地地圶圲 地. 地地地地地圱圩圮显然圬收敛是很慢的圮 最陡下降方向并不是理想的下降方向圬为了改进最优化的收敛速度圬需要寻求其他方 向圮为了判断一个寻找方向的好坏圬应该有一个标准圬注意到一般函数在最小点附近近似 于二次函数圬因此如果一个方法对二次函数比较有效圬则可望其对一般函数也比较有效圬 如果一个方法对二次函数的效果都不好的话圬很难期望它对一般函数会有好的效果圮计

30 圳地坃坈坁坐坔坅坒圲圮物理学中的常用数值方法算实践也证实了这一论断圮下面圬我们就从二次函数寻求一种比较有效的下降方向圮为此圬我们先给出几个数学定理圻定理一场若 N 维欧几里德空间中的向量 q 与 N 个线性独立向量 p 1, p 2,, p N 都正交圬则向量 q 必为地圮这一定理的结论是显然的圮定义场 A 共轭向量圮设 A 是一个 N N 的对称正定矩阵圬 p 圬 q 是两个 N 维向量圬若圬 p T Aq 地 在圲圮圲圮圴圱圩 则称向量 p 和向量 q 互为 A 共轭或互为 A 正交圮定理二场若 A 为 N N 对称正定矩阵圬 p 1 圬 p 2 圬 圬 p N 为 A 共轭的 N 维非零向量圬则此向量组必为线性独立圮证明场设向量组 p 1 圬 p 2 圬 圬 p N 之间存在线性关系 α 1 p 1 圫 α 2 p 2 圫 圫 α N P N 地对 i 圱, 圲,, N 圬用 p T i A 左乘上式得 α i p T i Ap i 地 因为 p i 地圬 A 正定圬所以 p T i Ap i > 地 从而必有 α i 地, i 圱, 圲,, N 因此圬向量组 p 1 圬 p 2 圬 圬 p N 线性独立圮 有了以上数学准备圬我们来考虑二次函数的最优化问题圬我们的目的是找一种下降方 向圬能够在有限步达到二次函数的极小点圮考虑二次 N 元函数 F 在 x 圩 a 圫 b T x 圫 其中 A 为一正定矩阵圮函数在圲圮圲圮圴圴圩的梯度为 圱 x T Ax 圲 在圲圮圲圮圴圲圩 G 在 x 圩 b 圫 Ax 在圲圮圲圮圴圳圩 记 x i 点的梯度为 g i G 在 x i 圩圬从 x 0 点出发圬对于第一个寻找方向圬除梯度外没有别的信息圬 我们就取 p 0 g 0 圬沿此方向找到极小点 x 1 并可求得这一点的梯度 g 1 圬因为 x 1 是沿 p 0 方 向的极小点圬这一点的梯度方向一定与 p 0 垂直,g 1 p 0 圬从而 g 1 g 0 圮现在有了两个方

31 圲圮圲圮方程的求根和最优化方法圳圱向,g 1 和 g 0 圬最陡下降法选取新的方向为 g 1 圮我们保留原来方向的一点信息, 选新的寻找方向为这两个方向的组合 p 1 g 1 α 0 g 0 在圲圮圲圮圴圴圩并要求 p 1 与 p 0 为 A 共轭圬则应有在 g 1 α 0 g 0 圩 T Ap 0 地在圲圮圲圮圴圵圩由于 x 1 x 0 圫 t 0 p 0 则 g 1 g 0 A 在 x 1 x 0 圩 t 0 Ap 0 一般圬若第 i 步的寻找方向为 p i 圬则同理可证 因此圬方程在圲圮圲圮圴圵圩成为 g i+1 g i A 在 x i+1 x i 圩 t i Ap i T 在 g 1 α 0 g 0 圩在 g 1 g 0 圩地 在圲圮圲圮圴圶圩 在圲圮圲圮圴圷圩 由此解得 α 0 g T 1 g 1 g T 0 g 0 这里假定 g 0 地圬否则圬若 g 0 地, 则 x 0 就是所求的极值点圮 在圲圮圲圮圴圸圩 把 α 0 代入在圲圮圲圮圴圴圩, 便得到 p 1 圮沿 p 1 求出极小点 x 2 并算出 g 2 圮由于 p 0 和 p 1 为 A 共轭圬因 此由在圲圮圲圮圴圶圩 g T 0 在 g 2 g 1 圩地 考虑到 g 0 与 g 1 正交圬 g T 0 g 1 地圬由上式可知 g 0 与 g 2 也正交圮因 x 2 是沿 p 1 方向的极小点, 故 g 2 p 1 圬也就是 g T 2 在 g 1 α 0 g 0 圩地 即 g 2 g 1 圮 g 0 圬 g 1 和 g 2 构成一个正交矢量组圬我们可以在这个正交组构成的三维空间寻 求与 p 0 圬 p 1 均为 A 共轭的寻找方向 p 2 圮令 p 2 g 2 α 1 g 1 β 0 g 0 在圲圮圲圮圴圹圩 要求 p 2 分别与 p 0 和 p 1 为 A 共轭圬也就是要求如下方程成立 T 在 g 2 α 1 g 1 β 0 g 0 圩在 g 1 g 0 圩地在圲圮圲圮圵地圩 T 在 g 2 α 1 g 1 β 0 g 0 圩在 g 2 g 1 圩地在圲圮圲圮圵圱圩 在圲圮圲圮圵圲圩

32 圳圲坃坈坁坐坔坅坒圲圮物理学中的常用数值方法 由此解出 于是 α 1 g T 2 g 2 g T 1 g 1, β 0 α 0 α 1. 在圲圮圲圮圵圳圩 p 2 g 2 圫 α 1 p 1 在圲圮圲圮圵圴圩 现假定已求得 p 0 圬 p 1 圬 圬 p j 圬其中 p 0 g 0 p i g i 圫 α i 1 p i 1 α i 1 gi T g i /gi 1g T i 1 i 圱, 圲,, j 这里假定 g 0 圬 g 1 圬 圬 g j 都不为零且构成一个正交系圮 沿 p j 求出 x j+1 并算出 g j+1 G 在 x j+1 圩圬设 g j+1 地圬则仿前可证圬 g 0 圬 g 1 圬 圬 g j+1 也构成正交系圮令 p j+1 g j+1 β 0 g 0 β 1 g 1 β j 1 g j 1 α j g j 在圲圮圲圮圵圵圩 并要求 p j+1 与 p 0 圬 p 1 圬 圬 p j 为 A 共轭圬即 T 在 g j+1 β 0 g 0 β 1 g 1 β j 1 g j 1 α j g j 圩在 g j+1 g j 圩地 在圲圮圲圮圵圶圩 得到 α j gj+1g T j+1 /gj T g j, β j 1 α j α j 1, β 1 α j α j 1 α 1, β 0 α j α j 1 α 1 α 0. 且 β 0 g 0 β 1 g 1 β j 1 g j 1 α j g j α j α 1 在 α 0 g 0 圫 g 1 圩 β 2 g 2 α j g j α j 在 α 1 p 1 g 2 圩 α j g j α j p j

33 圲圮圲圮方程的求根和最优化方法圳圳 即 p j+1 g j+1 圫 α j p j, α j g T j+1g j+1 /g T j g j 在圲圮圲圮圵圷圩 根据归纳法圬在在圲圮圲圮圵圷圩中令 j 地, 圱, 圲, 圬所得方向为 A 共轭方向圮上述方法最多只要 N 步就可找到函数在圲圮圲圮圴圴圩的极小点圬这是因为所有的 g j 互相正交圬即 gn T g j 地对所有 j 地, 圱, 圲,, N 圱都成立圬而我们的问题在 N 维空间圬故必有 g N 地圬亦即 x N 为函数在圲圮圲圮圴圴圩的极小点圮这就是所谓共轭梯度方法圬理论上圬对于二次函数圬只要 N 次迭代就可以达到极小点圬但实际计算时圬由于舍入误差的影响圬一般 N 次迭代并不能达到极小点圮对于一般的目标函数圬有限次在 N 次圩的迭代也一般不能达到极小点圬因此实际计算时圬从一初始点出发圬迭代 N 次后圬以所得结果为出发点圬重新开始计算圬直到收敛为止圮共轭梯度方法由于占用内存小圬方法简单圬编程方便圬目前在物理问题中应用比较多圮作为例子圬我们再考虑前面的例题圬对函数 x 2 圫圱地 y 2 T T 圬取 x 0 在圲, 圱圩圬则 g 0 在圴, 圲地圩圬 T T 取 p 0 g 0 圬求得 x 1 在圱. 圷圹圲圸圳, 地. 地圳圵圸圵圶圶圩圬 g 1 在圳. 圵圸圵圶圶, 地. 圷圱圷圱圳圱圩圬 α 0 T T 地. 地圳圲圱圴圲圳圬由此可构造出 p 1 在 圳. 圷圱圴圲圳, 地. 地圷圴圲圸圴圵圩圬沿 p 1 方向求得 x 2 在地. 地地地地地地, 地. 地地地地地地圩圬即已经达到极小点圮 习题 : 编制共轭梯度法的程序并求下面函数的极小值点及极小值, n [ F 在 x 1,, x n 圩圱圫圱地地在 x i x i 1 圩圫在圱 x i 1 圩 ] i=2 取 n 圱地, 选择初始点为 x i i/ 圱地, i 圱, 圲,, 圱地. 在圲圮圲圮圵圸圩 约束最优化 前面几节所讲的是无约束最优化问题圬在实际物理问题中圬还经常会遇到有约束最优化 问题圬其数学提法是圬给定一目标函数 F 在 x 1, x 2,, x n 圩圬在一定的约束下求解其极小值圬 约束一般写成一组函数圮即 { F 在 x 1, x 2,, x n 圩坍坩坮 G i 在 x 1, x 2,, x n 圩地 在 i 圱, 圲,, m 圩 在圲圮圲圮圵圹圩 式中 m < n 圬 m 个方程 G i 在 x 1, x 2,, x n 圩代表 m 个约束圮解决约束优化问题的通常做法 是构造一个新的函数圬把约束优化问题改为一无约束最优化问题圬然后利用前面的方法 求解圮通常的做法是圬定义 F F 圫 m λ i G 2 i i=1 在圲圮圲圮圶地圩

34 圳圴坃坈坁坐坔坅坒圲圮物理学中的常用数值方法 对于合适选择的 λ i 圬在 i 圱, 圲,, m 圩圬 F 的极小值就代表了由方程在圲圮圲圮圵圹圩定义的约束最 优化问题圮 例题在参数 x 1 圬 x 2 圬 圬 x n 的一个限制区域 D 内求解 F 的极小值圮为此圬我们把求解 区域扩展到整个区域圬如果 F 在 D 外无定义圬则应扩充其定义圮同时定义函数 G 为 { 地 {x} 坄 G i λ 在圲圮圲圮圶圱圩 ix 2 i {x} 坄 则通过计算 F 圫 G 的极小值就可得到问题的解圮在使用这种方法时应注意如果在 D 外 F 会 变小时圬 G 的选择应使得 F 圫 G 在 D 外总是较快增加的圮式在圲圮圲圮圶圱圩给出的仅仅是一种可 能的选择圬实际上 G 的选择可以有无穷多种圮 最小二乘法及曲线拟合 在物理实验中圬常常要测量一对物理量之间的关系圬在有些情况下圬一对物理量之间的函 数关系是已知的圬但包含若干个未知参数圬若用 x 和 y 代表这一对物理量圬则有 y f 在 x 圻 a 1, a 2,, a m 圩 在圲圮圲圮圶圲圩 式中 a 1 圬 a 2 圬 圬 a m 为一组未知参量圬 f 的函数形式是已知的圬在这种情况下圬曲线拟合的任务圬就是根据在 x, y 圩的 N 个观测点在 x 1, y 1 圩, 在 x 2, y 2 圩,, 在 x N, y N 圩来确定参数 a i 圬在 i 圱, 圲,, m 圩圬如果测量中没有误差圬则每一组在 x i, y i 圩都应准确地落在理论曲线上圬即应有 y i f 在 x i 圻 a 1, a 2,, a m 圩在 i 圱, 圲,, N 圩 在圲圮圲圮圶圳圩 为了确定诸 a i 的值圬只要选出在圲圮圲圮圶圳圩中的 m 个方程求解即可圬但实际测量总是有误差的圬在 x i, y i 圩并不准确落在理论曲线上圬在 N > m 的情况下圬方程组在圲圮圲圮圶圳圩成为矛盾方程组圬不能用解方程的方法求参数 a 1, a 2,, a m 圬只能用曲线拟合的方法根据带有误差的观测值点求得理论曲线参数的估计值圮在很多实际问题中圬 x 和 y 这两个物理量中总有一个量的测量精度比另一个高得多圬其测量误差可以忽略圮把可以忽略测量误差的量选作自变量 x 圬其观测值可看作是准确的圬对于每个 x i 圬对应的 y 的观测值 y i 是一个随机变量在见第圵章圩圬其误差可由其方差来估计圮物理中还经常遇到另外一种曲线拟合的任务场物理量 y 和 x 间的函数形式未知圬需要从观测点求出 y 和 x 的函数关系的一个经验公式圮一般的作法是选取某一函数族 { φ i 在 x 圩 } 的一个线性组合圬 y a i φ i 在 x 圩在圲圮圲圮圶圴圩 i

35 圲圮圲圮方程的求根和最优化方法圳圵 然后用曲线拟合方法求出系数 a i 圮函数族 { φ i 在 x 圩 } 通常从物理考虑来选取圬一种常用 的方法是选取 { 圱圬 x 圬 x 2 圬 圬 x m 1 } 用于拟合 m 个参数圻也可选取三角函数 { 圱圬坳坩坮 x 圬 坣坯坳 x 圬坳坩坮圲 x 圬坣坯坳圲 x 圬 圬坳坩坮圲 mx 圬坣坯坳圲 mx} 用于拟合圲 m 圫圱个参数圮这种方法称为线性 拟合方法圬即 y 线性地依赖于参数 a i 圬对于有些物理问题圬 y 和 x 的函数关系中非线性地依 赖于参数 a i 圬例如 y a 1 e a 2x 就是一个具有圶个参数的非线性拟合问题圮 圫 a 3 e a 4x 坳坩坮在 a 5 x 圫 a 6 圩 曲线拟合通常采用最小二乘法圬如果 x 的测量误差可以忽略圬观测点 i 的残差 δ i 定义 为在该点 y 的观测值与理论值之差圬即 δ i y i f 在 x i 圻 a 1, a 2,, a m 圩 在圲圮圲圮圶圵圩 若观测值 y i 的方差为 σ 2 i 圬则定义其权重因子为 w i 圱 σ 2 i 记 χ 2 为观测点残差的加权平方和圬 χ 2 N w i δi 2 i=1 N i=1 2 w i 在 y i f 在 x i 圻 a 1, a 2,, a m 圩圩 在圲圮圲圮圶圶圩 曲线拟合的最小二乘法就是以 χ 2 为目标函数的优化问题圬即寻找一组参数 a i 圬在 i 圱, 圲,, m 圩使 χ 2 为最小圮在观测点的方差未知的情况下圬可以取 w i 圱圬也可以在每一 观测点多测量几次圬以其平均值 y i 作为该点的观测值圬而以 yi 2 y i 2 作为其方差圮 在曲线拟合问题中圬最困难的是选取合适的经验公式圬这需要结合物理考虑和对测量 数据的观察圬也依赖于对各种函数及其图像的明确概念圮 例题在某一能谱测量中圬得到下表所示数据圬试给出经验公式并作出曲线拟合圮 坸 地圮地 地圮圱 地圮圲 地圮圳 地圮圴 地圮圵 地圮圶 坹 圲圮圱圱圷 圲圮圶圳圴 圳圮圴圵圹 圴圮圶圷圱 圶圮圶地地 圹圮圹圷圸 圱圵圮圸圹圳 坸 地圮圷 地圮圸 地圮圹 圱圮地 圱圮圱 圱圮圲 圱圮圳 坹 圲圵圮圱地圸 圳地圮圹圷圹 圲圵圮圷圶圸 圱圷圮圳圶圳 圱圲圮地圳圱 圱地圮圴圵圴 圱地圮圲圸地 坸 圱圮圴 圱圮圵 圱圮圶 圱圮圷 圱圮圸 圱圮圹 坹 圱圱圮地圱圲 圱圱圮圱圳圸 圹圮圷圴圵 圷圮圵圷圶 圵圮圵圹圸 圴圮圱圶圵 把表中数据点到坐标纸上在见图圲圮圱中圈点圩圬可以看出数据有二个峰圬为此圬我们试 用坌坯坲坥坮坴坺形曲线去拟合圬即取经验公式为 y a 1 在 x a 2 圩 2 圫 a 2 3 a 4 圫在 x a 5 圩 2 圫 a 2 6

36 圳圶 坃坈坁坐坔坅坒圲圮 物理学中的常用数值方法 用上式和上表中的数据计算 χ 2 圬并用共轭梯度方法优化 χ 2 圬求得各参数值为 a 1 圱. 圲圱地 a 2 地. 圸地地 a 3 地. 圲地圲 a 4 地. 圷圷圶 a 5 圱. 圵地圲 a 6 地. 圲圹圵 在图圲圮圱中圬我们用实线画出了拟合的结果圬虚线给出了二个峰的图像圮 坆坩坧坵坲坥圲圮圱场 2.3 函数的求值 在数值计算中圬经常要求一些函数的值圬一些常用的初等函数如三角函数圬 e 指数函数圬双曲函数圬对数函数等圬均已写入计算机语言的定义做为内部函数圬在使用时只需调用即可圬而其它的函数和表达式则要自己计算圬本章将给出一些在物理计算中经常遇到的函数及表达式的计算方法和部分程序圮一般说来圬计算机的内部函数都是由专家设计的圬其效率和精度都很好圬因此圬如果计算机语言已经提供了所需的函数圬则应尽量使用这些函数圬下面提供的方法仅仅作为计算机系统不提供这些函数时的补充圮 多项式的求值和级数求和 所谓多项式是指由下式给出的表达式场 p 在 x 圩 a 0 圫 a 1 x 圫 a 2 x 2 圫 a 3 x 3 圫 圫 a n x n

37 圲圮圳圮函数的求值圳圷 为了计算 p 在 x 圩在给定 x 时的数值圬我们用 在 n 圩圬 n 地, 圱, 圲,, n 来代表多项式的系数圬则 可以写出计算语句圮以 n 圴为例圬直接用坆坏坒坔坒坁坎语言写出的算式为场 p=(0)+(1)*x+(2)*x**2+(3)*x**3+(4)*x**4 这显然不是一个好的算法在为什么圿圩圮较好的算法是下面二者之一场 p=(0)+x*((1)+x*((2)+ x*((3)+x*(4)))) p=((((4)*x+(3))*x+(2))*x+(1))*x+(0) 当 n 较大时圬可采用下面的语句段场 P=C(N) DO 11 J=N-1,1,-1 P=P*X+C(J) 11 CONTINUE 在科研中圬还常常碰到级数求和的问题圬例如计算 S i=0 u i 我们当然不可能计算到无穷多项圬所以上述级数必需在某处截断圬如果级数收敛很慢圬就要计算很多项圬一方面占用大量的坃坐坕时间圬更重要的是当计算的项数非常多时圬舍入误差的积累有可能完全破坏整个计算圮为此必须使用加速收敛的技巧圮坁坩坴坫坥坮算法是常用的加速技巧之一圮它是对相邻的三个部分和进行处理圬得到一个更好的近似圮记部分和 为场 S n 则由 S n 1, S n, S n+1 可计算出 n i=0 u i 在圲圮圳圮圶圷圩 S n 2 在 S n+1 S n 圩 S n+1 S n+1 圲 S n 圫 S n 1 在圲圮圳圮圶圸圩 对 n 圱圬圲圬圳圬 圬先求出 S n 圬同时算出 S n 圬一般 S n 更快地收敛于 S 圮在实际使用时圬在圲圮圳圮圶圸圩最好用写出的式子计算圬一些代数上等价的式子可能导致较大的舍入误差圮 Euler 变换对于由下式定义的交错级数 在 圱圩 s u s 在圲圮圳圮圶圹圩 这里 u s > 地圮坅坵坬坥坲变换是一个强有力的算法圬它由下式给出在取 n 为偶数圩 s=0 在 圱圩 s u s u 0 u 1 圫 u 2 u n 1 圫 s=0 s=0 s 在 圱圩在圁 s u 圲 s+1 n 圩 在圲圮圳圮圷地圩

38 圳圸坃坈坁坐坔坅坒圲圮物理学中的常用数值方法 其中场 圁 u n u n+1 u n 圁 2 u n 圁在圁 u n 圩 u n+2 圲 u n+1 圫 u n 在圲圮圳圮圷圱圩 为向前差分圮式在圲圮圳圮圷地圩可推导如下场 定义位移算符 E 为场 u n+1 Eu n 由于在圱圫圁圩 u n u n 圫 u n+1 u n u n+1 Eu n 圬因此有算符等式 E 圱圫圁 而 u 0 u 1 圫 u 2 u 3 圫 在圱 E 圫 E 2 E 3 圫 圩 u 0 圱圲 圱 圱圫 2 u 0 s=0 s 在 圱圩圁 s u 圲 s+1 0 圱圱圫 E u 连分式的求值 所谓连分式圬是指由下式给出的表达式圬 f b 0 圫 b 1 圫 b 2 圫 a 1 a 2 a 3 a 4 b 3 圫 a 5 b 4 圫 b 5 圫 数学家们对它做过很多研究圬这是一个非常有趣的函数表示方法圮除了上面的表达式外圬人们还经常使用另一种记法圮 a 1 a 2 a 3 a 4 a 5 f b 0 圫 b 1 圫 b 2 圫 b 3 圫 b 4 圫 b 5 圫 这里圬系数 a n 和 b n 也可以是 x 的函数圮由于连分式只能从右向左求值圬因此在具体计算之前圬很难判断一个无限的连分式应该在何处截断圬一种直接的方法就是猜一个截断点圬求出一个值圬再猜一个更大一点的截断点计算圬通过比较不同截断处的结果来判断是否收敛圮这显然不是一个好的方法圮另一

39 圲圮圴圮常微分方程的数值解法圳圹 种在实践中常用的方法是用一个分式有理函数来近似连分式圬用 f n 表示计算到 a n 和 b n 的 连分式的值圬则场 这里 A n 圬 B n 由下面的递推公式给出场 练习 A 1 圱 f n A n B n B 1 地 在圲圮圳圮圷圲圩 A 0 b 0 B 0 圱在圲圮圳圮圷圳圩 A j b j A j 1 圫 a j A j 2 B j b j B j 1 圫 a j B j 2 j 圱, 圲,, n 试证明递推关系 (2.3.73) 并写出计算机程序 2.4 常微分方程的数值解法 常微分方程在物理学中具有重要意义圬牛顿第二定律就表示为一组二阶常微分方程圬第一章中提到的著名的三体问题就是由九个二阶常微分方程组成的方程组圻在量子力学中和电动力学中圬由薛定谔方程及坍坡坸坷坥坬坬方程分离变量得到的所谓坓坴坵坲坭坻坌坩坯坵坶坩坬坬坥本征值问题也是二阶常微分方程的求解问题圮在这一章中圬我们将给出几种数值求解常微分方程的常用方法圮 Euler 方法 对于常微分方程的初值问题 dy dx f 在 x, y 圩 y 在地圩 η 在圲圮圴圮圷圴圩 数值求解的第一步是把方程在圲圮圴圮圷圴圩变为一个差分方程圬为此圬让 x 取一组分立值 x n 圬 n 地圬圱圬圲圬 圬 x 0 地圬 x n+1 x n 圫 h 圬这里 h 称为积分的步长圬把 y 在 x 圩在 x n 处展开圬有 y 在 x n+1 圩 y 在 x n 圩圫 h dy 圫 O 在 h 2 圩在圲圮圴圮圷圵圩 dx x=xn 利用方程在圲圮圴圮圷圴圩圬便可得到 y 在 x n+1 圩 y 在 x n 圩圫 hf 在 x n, y n 圩圫 O 在 h 2 圩在圲圮圴圮圷圶圩这就是坅坵坬坥坲差分格式圬从 x 地出发圬利用这一方程一步一步地向前计算圬就可以得到任一 x n 处的函数值 y n y 在 x n 圩圮坅坵坬坥坲方法是一个显式方法圬即第 n 圫圱步的 y n+1 只依赖

40 圴地坃坈坁坐坔坅坒圲圮物理学中的常用数值方法于 x n 圬 y n 的值圬而且每一步只需计算一次 f 在 x n, y n 圩的值圬这是坅坵坬坥坲方法的优点圬但另一方面圬坅坵坬坥坲方法也存在严重的不足圬因为这一方法是一次的在每一步的误差项为 h 2 的数量级圬 圁 y 在 x N 圩 Nh 2 h 圩圬所以 h 必须取得足够小圬以保证有满意的精度圬但步长太小时圬积分到一预定的 x 则要做很多步的计算圬舍入误差的积累将会使计算结果严重失真圮由于这一原因圬在实际计算中几乎不使用坅坵坬坥坲方法圬而常常利用其推导上的简单和概念上非常清楚而用于教学圮 龙格 库塔方法 具有二阶以上精度的直接方法称为龙格坼库塔方法圮做为一种训练圬我们在此给出二阶 龙格坼库塔方法的推导圮下面先给出一种比较直接的推导圬然后再介绍一种易于推广到 高阶情形的较为系统的推导方法圬希望同学仔细体会这里介绍的方法和思路圬以求举一 反三圬为自己设计算法做准备圮 第一种推导方法场如果计算 x n+ 1 2 y n+ 1 2 x n 圫 h 2 处 y 的一阶导数圬可以得到 y n+1 y n h 这可用坔坡坹坬坯坲展开加以验证圮把上式稍加整理圬并代入 y n+ 1 2 圫 O 在 h 2 圩在圲圮圴圮圷圷圩 f 在 x n+ 1 圬 y 2 n+ 1 圩圬得到 2 y n+1 y n 圫 hf 在 x n+ 1, y 2 n+ 1 圩圫 O 在 h 3 圩在圲圮圴圮圷圸圩 2 上式中出现了半个分点处的值圬我们必须设法用分点 n 和 n 圫圱处的量来表示 n 圫 1 的量圬为此圬把 f 在 x n+ 1, y 2 n+ 1 圩和 f 在 x n+1, y n+1 圩在 n 点展开场 2 f 在 x n+ 1 2, y n+ 1 圩 f 在 x n, y n 圩圫 2 h df 圲 dx f 在 x n+1, y n+1 圩 f 在 x n, y n 圩圫 h df dx 把上面的第一式乘以圲减去第二式圬整理后得 注意到 f 在 x n+ 1 2, y n+ 1 圩 2 xn xn 圫 O 在 h 2 圩 2 处 圫 O 在 h 2 圩在圲圮圴圮圷圹圩 圱在 f 在 x n, y n 圩圫 f 在 x n+1, y n+1 圩圩圫 O 在 h 2 圩在圲圮圴圮圸地圩圲 y n+1 y n 圫 hf 在 x n, y n 圩圫 O 在 h 2 圩 及 f 在 x n+1, y n 圫 hf 在 x n, y n 圩圫 O 在 h 2 圩圩 f 在 x n 圫 h, y n 圫 hf 在 x n, y n 圩圩圫 O 在 h 2 圩 我们得到 y n+1 y n 圫 h 在 f 在 x n, y n 圩圫 f 在 x n 圫 h, y n 圫 hf 在 x n, y n 圩圩圫 O 在 h 3 圩在圲圮圴圮圸圱圩圲

41 圲圮圴圮常微分方程的数值解法圴圱这就是二阶龙格坼库塔格式的一种圬由它可以通过 n 点的解给出 n 圫圱点的解圮它显然具有二阶精度圮第二种推导方法场因为我们的目的是从 n 点的解给出 n 圫圱点的解并要求每一步的误差为 O 在 h 3 圩圬所以可以一般的把差分格式写为 y n+1 y n 圫 h 圁在 x n, y n 圻 h 圩 在圲圮圴圮圸圲圩 其中 圁在 x n, y n 圻 h 圩 y n 圫 h y n 圫圲 如果我们用另一函数圈在 x n, y n 圻 h 圩代替圁在 x n, y n 圻 h 圩圬且使得 则可得到 p 阶差分格式圮 对于 p 圲圬取 h 2 y n (3) 圫 圳圡 圈在 x n, y n 圻 h 圩 圁在 x n, y n 圻 h 圩 O 在 h p 圩 圈在 x, y 圻 h 圩 1 f 在 x, y 圩圫 2 f 在 x 圫 ha 2, y 圫 b 21 hf 在 x, y 圩圩 把上式对 h 展开 [ ] f f 圈在 x, y 圻 h 圩 1 f 在 x, y 圩圫 2 f 在 x, y 圩圫 h 2 a 2 圫 b 21 x y f 在 x, y 圩 另一方面场 圁在 x, y 圻 h 圩 f 在 x, y 圩圫 圫 O 在 h 2 圩在圲圮圴圮圸圳圩 [ ] 圱 f f h 圫圲 x y f 在 x, y 圩圫 O 在 h 2 圩在圲圮圴圮圸圴圩 比较式在圲圮圴圮圸圳圩和在圲圮圴圮圸圴圩圬只要在圲圮圴圮圸圳圩中的系数满足下述关系圬就能达到我们的要求圮 1 圫 2 圱 圱 2 a 2 在圲圮圴圮圸圵圩圲 圱 2 b 21 圲 方程在圲圮圴圮圸圵圩的解并不唯一圬也就是说圬有无穷多组系数可以满足要求圮为此圬令 2 α 圬 则 取 α 1 2 圬得到 1 1 圱 α a 2 b 21 圱圲 α 圱圲 α 圱, a 2 圱, b 21 圱圲 在圲圮圴圮圸圶圩

42 圴圲坃坈坁坐坔坅坒圲圮物理学中的常用数值方法 差分格式为 y n+1 y n 圫 这与第一种方法得到的结果相同圮若取 α 圱圬则 h 在 f 在 x n, y n 圩圫 f 在 x n 圫 h, y n 圫 hf 在 x n, y n 圩圩圲 在圲圮圴圮圸圷圩 1 地, a 2 圱, b 21 圲 圱圲 差分格式为 y n+1 y n 圫 h 在 f 在 x n 圫 圱圱 h, y n 圫 hf 在 x n, y n 圩圩圲圲 当 α 取不同值时圬我们还可得到其它形式的二阶差分格式圮 在圲圮圴圮圸圸圩 计算实践表明圬对于绝大多数简单问题圬四阶龙格坼库塔方法是最好的选择圮这一方 法的一种常用差分格式可用下面一组公式来表示场 k 1 hf 在 x n, y n 圩 ( ) h k 1 k 2 hf x n 圫, y n 圫圲圲 ( ) h k 2 k 3 hf x n 圫, y n 圫圲圲 k 4 hf 在 x n 圫 h, y n 圫 k 3 圩在圲圮圴圮圸圹圩 y n+1 y n 圫 k 1 圶 圫 k 2 圳 圫 k 3 圳 圫 k 4 圶 圫 O 在 h 5 圩在圲圮圴圮圹地圩 实际上圬四阶的龙格坼库塔方法的形式并不唯一圬这里给出的是常用的一种圮四阶龙 格坼库塔方法的推导比较繁琐圬但并不难圬其基本思路与前面推导二阶方法的第二种 方法相同圬块圮坃圮均坥坡坲的名著坎坵坭坥坲坩坣坡坬坩坮坩坴坩坡坬坖坡坬坵坥坐坲坯坢坬坥坭坳坩坮坏坲坤坩坮坡坲坹坄坩國坥坲坥坮坴坩坡坬 坅坱坵坡坴坩坯坮坳在坅坮坧坬坥坷坯坯坤坃坬坩國坳圬坎圮坊圮坐坲坥坮坴坩坣坥坻坈坡坬坬圬圱圹圷圱圩给出了一般龙格坼库塔方法的 推导圬同时也列出了其它四阶龙格坼库塔方法的形式圬有兴趣的同学可以参考圮 上面给出了用四阶龙格坼库塔方法积分一步的方法圬我们当然可以利用这一方法圬选 定一个步长 h 圬从初绐点开始一步一步地积分到所需的点圮那么圬步长如何选取呢圿选多 大才是合适的圬是否每一步都选相同的步长圿为此圬我们需要一个判定精度的方法圬同时圬 直观告诉我们圬在函数变化剧烈的地方圬步长应该取得小一点以保证精度圬而在函数变化 平缓的地方圬步长可以取的大一点以减小计算量圮对于每一步圬判定步长是否合适的一个 简单的方法就是把步长减半圬比较步长为圲 h 的一步的积分结果和步长为 h 的两步的计算 结果圬确切地说圬 y 在 x 圫圲 h 圩 y 1 圫 O 在 h 5 圩 y 在 x 圫 h 圫 h 圩 y 2 圫 O 在 h 5 圩在圲圮圴圮圹圱圩

43 圲圮圴圮常微分方程的数值解法圴圳 第二个方程是指用二步步长为 h 的积分结果圮二者之差 圁 y 2 y 1 在圲圮圴圮圹圲圩 显然可以作为误差的一个合适的量度圮由方程在圲圮圴圮圹圱圩可知圬圁正比于 h 5 圬因此圬如果步长为 h 1 时的误差为圁 1 圬而当步长为 h 0 时的误差为圁 0 圬则显然应有 h 0 h 1 1/5 圁 0 圁 1 在圲圮圴圮圹圳圩 如果我们令圁 0 为充许误差圬则方程在圲圮圴圮圹圳圩可以用以调整步长圬如果 圁 1 大于 圁 0 圬则 h 0 可作为新的步长重新计算圬否则圬如果 圁 1 小于 圁 0 圬则 h 0 可作为下一步步长的一个合理的推断圮 6 下面考虑如何选取圁 0 圬你可以把它取为某个小的数字圬例如圬圱地圬但略加思考就会发现这是不行的圬理由如下圬所求函数 y 的绝对值可以很大圬也可以很小圬因此圬绝对精度是没有意义的圻所以圬我们可以考虑取圁 0 ɛy 圬而 ɛ 为某一小的数字圻但另一方面圬有些函数是振荡的圬在某一极大和极小之间圬并且经常穿过地点圬对于这类问题圬圁 0 的一个合适的选则应为某一小量乘以 y 的极大值的绝对值圻为了照顾到上述两种情况圬我们可以取 圁 0 ɛ y sal 在圲圮圴圮圹圴圩 其中 y sal y 圫 h dy dx 上面是以一步为基础来考虑的圬在有些问题中圬我们可能要求解满足某种全局的精度圬 由于舍入误差的积累在在最坏的情况下圬舍入误差全部相加圬总体误差为 N 圁 0 L h 圁 0 圬其 中 L 为对 x 的积分的长度圩圬当 h 缩小时圬圁 0 也要相应的缩小圬以保证总的舍入误差固定圬 即选取 y 坳坣坡坬正比于 h 圬所以当 h 0 变化时圬 y sal 也相应的变化圬从而方程在圲圮圴圮圹圳圩中的指数 成为圱 / 圴圬另外圬我们只考虑了误差的主导项圬高阶项也应有一些贡献圮基于所有这些考 虑圬我们把方程在圲圮圴圮圹圳圩改写为 h 0 Sh 1 Sh 1 1/5 圁 0 圁 1 1/4 圁 0 圁 1 圁 0 圁 1 圁 0 < 圁 1 在圲圮圴圮圹圵圩 这里 S 为一考虑了高阶效应的安全系数圬指数的选择从最安全出发圮当放大步长时圬用小 指数圱 / 圵在少放一点圩圻当缩小步长时圬用大指数圱 / 圴在多缩一些圩圮圁 0 由在圲圮圴圮圹圴圩给出圮 为了简洁起见圬这里以一个方程为基础给出了四阶龙格坼库塔方法圬但不难把它们推 广到多个方程圮

44 圴圴坃坈坁坐坔坅坒圲圮物理学中的常用数值方法在结束本小节前圬我们指出圬虽然这里讲的是一阶常微分方程组的解法圬但实际上已经包括了高阶常微分方程圬因为任何一个高阶常微分方程均可化为一组等价的一阶常微分方程组圬如果有同学对这一表述还不清楚圬那么请立即复习 高等数学 的有关内容圮本节的方法基本上可以满足实际计算的需要圮但在实际工作中圬也会碰到一些特殊的问题圬比较重要的一类问题是所谓刚性微分方程求解问题圬关于这一点圬请参看块圮坃圮均坥坡坲的著作圮另外一些特殊问题将在后面结合具体物理问题进行处理圮 练习 : 单摆的方程为 x 坿圫 ω 2 坳坩坮 x 地 设初始位移为 x 地, 试用二阶龙格 - 库塔方法对不同的初始速度求解此方程, 并讨论所得结果 两点边值问题 所谓两点边值问题是指下面的定解问题 在点 x 1 圬有边界条件 dy i dx g i 在 x 圻 y 1, y 2,, y n 圩 i 圱, 圲,, n 在圲圮圴圮圹圶圩 B 1j 在 x 1 圻 y 1, y 2,, y n 圩地, j 圱, 圲,, n 1 在圲圮圴圮圹圷圩 在点 x 2 圬有边界条件 B 2j 在 x 2 圻 y 1, y 2,, y n 圩地, j 圱, 圲,, n n 1 在圲圮圴圮圹圸圩这样一类问题与前节所处理的问题的最大差别在于定解条件是在 x 的两个不同的点在端点圩处给出的圬而前节的定解条件是在一个点给出的圬因此圬前节的问题称为初值问题而本节将要解决的问题称为两点边值问题圮初值问题一般来说总是有解的圬但边值问题的解有时并不存在圬边值问题中常常包含一个参数圬只有对参数的某些值圬边值问题的解才存在圬这类问题称为微分方程的本征值问题圬使得方程有解的那些参数值称为本征值圮边值问题有各种不同的解法圬我们将在本节介绍打靶法圮并通过一个具体算例给出处理边值问题的一些技巧和方法圮在后面的章节中圬我们将结合具体物理问题圬再回到本节的问题圮这里要强调指出的是圬由于边值问题的特殊性圬试图像初值问题那样给出通用程序是几乎不可能的圬在实际工作中应根据具体的物理问题设计合适的方法圮打靶法的基本思路是从 x 1 点出发圬利用 n 1 个已知的边界条件并指定 n 2 个条件圬积分到 x 2 点圬看是否与 x 2 点的 n 2 个条件一致圬若否圬则调节在 x 1 指定的条件直到符合要求

45 圲圮圴圮常微分方程的数值解法圴圵 T 为止圮现在我们给出仔细的分析圬令 V 坛 V 1, V 2,, V n2 坝为一 n 2 维空间的矢量圬它给出了 x 1 点的 n 2 个指定的条件圬则方程在圲圮圴圮圹圷圩可以写为 y i 在 x 1 圩 y i 在 x 1 圻 V 1, V 2,, V n2 圩 i 圱, 圲,, n 在圲圮圴圮圹圹圩因此圬给定 V 圬我们就有了一个完整的初值问题圬利用前节的方法在区间坛 x 1, x 2 坝之间积分方程在圲圮圴圮圹圶圩圬便可得到一组 y i 在 x 2 圻 V 1, V 2,, V n2 圩圬在 x 2 圬定义一个 n 2 维的偏离矢量 F 圬 F 的取法并不唯一圬可视方便而定圬例如可取 F k B 2k 在 x 2 圻 y 1, y 2,, y n 圩 k 圱, 圲,, n 2 在圲圮圴圮圱地地圩如果 F 地圬则说明 V 的选则是正确的圮因此圬问题转化为求解方程 F 地圬这里有 n 2 个未知量 V k 圬 n 2 个方程 F k 地圮为了求解这一组方程圬我们假定 V 与真实的解相差不大圬从而可以把方程在圲圮圴圮圱地地圩线性化圬即令 这里 由方程在圲圮圴圮圱地圱圩求出 δv j 圬并令新的 V 为 n 2 F k 在 V δv 圩 A kj δv j A kj F k V j j=1 在圲圮圴圮圱地圱圩 在圲圮圴圮圱地圲圩 V new k V old k 圫 δv k 在圲圮圴圮圱地圳圩 通过反复迭代圬最终得到结果圮 现在圬让我们回到在圲圮圴圮圱地圲圩圬由于 F k 是由数值给出的圬因此 A kj 的计算也只能是数值 的圬最简单的方法是令 F k V j F k 在 V 1, V 2,, V j 圫 δv j, 圩 F k 在 V 1, V 2,, V j, 圩 δv j 在圲圮圴圮圱地圴圩 这样圬每完成一次迭代需要求解 n 2 圫圱次初值问题圬一次计算 F k 在 V 1, V 2,, V n2 圩圬 n 2 次计 算偏导数圬显然圬与初值问题相比圬边值问题的计算量是相当大的圮 V 的初始叠代值的选 择是一个困难的问题圬如果选择的不好圬则可导致计算的失败圬通常应根据对被求解问题 的了解程度尽量选择接近其解的量作为初始叠代值圮在后面的例子中圬将给出一个针对 具体问题的选法圮 为了帮助理解打靶法的应用圬这里给出一个例子圬即旋转椭球谐振子圮当在旋转椭球 坐标系中对某些偏微分方程如坌坡坰坬坡坣坥方程分离变量时圬将会导致下面的方程 [ ] ) d 在圱 x 2 ds 圩圫 (λ 2 x 2 m2 S 地 dx dx 圱 x 2 在圲圮圴圮圱地圵圩

46 圴圶坃坈坁坐坔坅坒圲圮物理学中的常用数值方法这里 m 为一整数圮这是一个坓坴坵坲坭坻坌坩坯坵坶坩坬坬坥本征值问题圬 λ 是方程的本征值圮方程在 x ± 圱具有奇点圬从坓坴坵坲坭坻坌坩坯坵坶坩坬坬坥本征值问题的一般理论可知在参见梁昆淼圬 数学物理方法 圬如果我们寻求在 x ± 圱为正则的解圬则只有对 λ 的某些特殊值才有可能圬这些特殊值称为本征值圬且所有的本征值均大于地圮当 地时圬方程化为球谐振子的方程圬其解为 Pn m 在 x 圩圬本征值为 λ mn n 在 n 圫圱圩圬 n m, m 圫圱, 圬当 地时圬我们记在圲圮圴圮圱地圵圩的本征值为 λ mn 在 圩圬对应的本征函数记为 S mn 在 x 圩圮 λ mn 在 圩和 S mn 在 x 圩的计算曾经是非常困难的问题圬其十分复杂的级数解圬递推关系等可在一些非常专门的书籍中查到圮在例如圬可参看坍圮坁坢坲坡坭坯坷坩坴坺坡坮坤坉坓坴坥坧坵坮坈坡坮坤坢坯坯坫坯坦坍坡坴坨坥坭坡坴坩坣坡坬坆坵坮坣坴坩坯坮坳在坄坯坶坥坲坐坵坢坬坩坣坡坴坩坯坮坳圬坎坥坷坙坯坲坫圬圱圹圶圸圩圩其数值解相对来说要容易得多圮首先圬我们作代换在从天上掉下来的圿圬请复习 数学物理方法 中常微分方程的有关内容圮圩 S 在圱 x 2 圩 m/2 y 则原方程化为 这里 在圱 x 2 d 2 y 圩 dx 圲在 m 圫圱圩 x dy 圫在 µ 2 x 2 圩 y 地在圲圮圴圮圱地圶圩 2 dx µ λ m 在 m 圫圱圩在圲圮圴圮圱地圷圩 方程在圲圮圴圮圱地圴圩和在圲圮圴圮圱地圶圩均在代换 x x 下不变圬因此 S 及 y 除了一个常数因子外也应在同一变换下不变圬由于解最终将归一化圬因此这一常数因子只能是 ± 圱圮方程在圲圮圴圮圱地圶圩也是一个坓坴坵坲坭坻坌坩坯坵坶坩坬坬坥方程圬其本征值 µ 地圬参照球谐振子的结果圬如果令 λ mn α mn 圫 n 在 n 圫圱圩圬且要求 α mn 地圬则由在圲圮圴圮圱地圷圩可知 n m 圬 n m 对应于最低本征值圮利用坓坴坵坲坭坻坌坩坯坵坶坩坬坬坥方程的从小到大排列的本征值与其对应的本征函数的零点的定理可知圬当 n m 为偶数时圬 y mn 在坛 圱, 圱坝之间有偶数个零点圬所以我们取 y 在 x 圩 y 在 x 圩圬否则圬取 y 在 x 圩 y 在 x 圩圬于是 y mn 在 x 圩在 圱圩 n m y mn 在 x 圩在圲圮圴圮圱地圸圩 在 x 圱时圬可以解得 在原点圬我们有 y 在 圱圩 µ 2 y 在 圱圩圲在 m 圫圱圩 y 在地圩地, n m 为奇数 在圲圮圴圮圱地圹圩 y 在地圩地, n m 为偶数在圲圮圴圮圱圱地圩 最后圬我们令 S mn 在 x 圱时与 P m n 有相同的极限行为圮在这相当于选则一个特定的归一 化条件圩 坬坩坭在圱 x 2 圩 m/2 S mn 在 x 圻 圩坬坩坭在圱 x 2 圩 m/2 Pn m 在 x 圩在圲圮圴圮圱圱圱圩 x 1 x 1

47 圲圮圴圮 常微分方程的数值解法 圴圷 现在圬令 y 1 y y 2 y y 3 µ 在圲圮圴圮圱圱圲圩 这样圬方程在圲圮圴圮圱地圶圩变为 y 1 y 2 y 2 圱 [ ] 圲 x 在 m 圫圱圩 y 圱 x 2 2 在 y 3 2 x 2 圩 y 1 y 3 地在圲圮圴圮圱圱圳圩我们选择区间坛 圱, 地坝来求解在为什么不选坛地, 圱坝圿圩圬边界条件为圬在 x 圱圬有 y 2 y 3 2 y 1 圲在 m 圫圱圩 y 1 坬坩坭在圱 x 2 圩 m/2 Pn m 在 x 圩 x 1 在 圱圩 m 在 n 圫 m 圩圡圲 m m 圡在 n m 圩圡 在圲圮圴圮圱圱圴圩 在 x 地圬有场 y 1 地, n m 为奇数 y 2 地, n m 为偶数在圲圮圴圮圱圱圵圩 在这一问题中圬未知量是本征值 µ mn 圬为了给出一个好的初值圬我们注意到在方 程在圲圮圴圮圱地圶圩中圬如果把 µ 2 x 2 作为本征值看待圬则为球谐振子方程圬其本征值已知 为 n 在 n 圫圱圩 m 在 m 圫圱圩圬因为地 x 2 圱圬因此圬把 x 2 代之以 1 2 得到的结果应是原问题的一个良好的近似圬为此圬 µ mn 的初始叠代值可选为 n 在 n 圫圱圩 m 在 m 圫圱圩圫 圮在这个例子中圬我们不仅示范了打靶法的使用圬同时也示范了如何把高阶方程化为一 阶方程以及如何处理本征值问题圬请同学仔细体会圬掌握这一方法并应用于解决实际问 题圮 练习 : 对于方程 y 圫 圱 x y 圫 ky 地 在边界条件 y 在地圩有限, y 在圱圩地下计算本征值 k. 提示, 为了得到 k 的叠代初值, 考虑 k 很大时, 方程成为 y 圫 ky 地, 其解 为 y 坣坯坳在 kx 圫 δ 圩, 若令 δ 地 (δ 应在区间坛地, π 2 坝内 ), 则 k 在 n 圫 1 2 圩 2 π 2.

48 圴圸坃坈坁坐坔坅坒圲圮物理学中的常用数值方法 2.5 插值和逼近 给定一个列表函数圬在 x 1, f 在 x 1 圩圩, 在 x 2, f 在 x 2 圩圩,, 在 x n, f 在 x n 圩圩圬要求寻找一个函数 φ 在 x 圩来近似代替 f 在 x 圩圬使得在在 x 1, x 2,, x n 圩圬 f 在 x i 圩 φ 在 x i 圩圬而对于其它的 x 圬 φ 在 x 圩是 f 在 x 圩的一个很好的近似圮这样的问题就称为插值问题圮如果读者稍为细心一点圬就会发现我们的最后一个要求是有问题的圬因为我们仅仅知道 f 在 x 圩在一些分立点 x i 上的值圬事先并不知道 f 在 x 圩在这些点以外如何变化圬因此除了要求 φ 在 x 圩在诸 x i 点上与 f 在 x 圩相等外圬无法提出更多的要求圮然而圬在物理上和数学上圬假如我们认为一组分立点的数值大致决定了 f 在 x 圩圬这就隐含着我们假定了 f 在 x 圩在任何两个给定点之间是光滑变化的圮因此圬我们的目的实际上是寻找一条光滑地通过 n 个给定点在 x 1, f 在 x 1 圩圩, 在 x 2, f 在 x 2 圩圩,, 在 x n, f 在 x n 圩圩的一条曲线圮 多项式插值 通过两点可以唯一地作一条直线圬通过三点可以唯一地作一条抛物线圬圮圮圮圮圮圮圮通过 n 个 数据点在 x 1, y 1 圩, 在 x 2, y 2 圩,, 在 x n, y n 圩可以唯一地决定一个 n 圱阶的多项式 P n 在 x 圩圬这一多 项式就称为 y f 在 x 圩的插值多项式圮它由著名的坌坡坧坲坡坮坧坥公式给出圮 P 在 x 圩 在 x x 2 圩在 x x 3 圩 在 x x n 圩在 x x 1 圩在 x x 3 圩 在 x x n 圩 y 1 圫 y 2 圫在 x 1 x 2 圩在 x 1 x 3 圩 在 x 1 x n 圩在 x 2 x 1 圩在 x 2 x 3 圩 在 x 2 x n 圩在 x x 1 圩在 x x 2 圩 在 x x n 1 圩 y n 在圲圮圵圮圱圱圶圩在 x n x 1 圩在 x n x 2 圩 在 x n x n 1 圩 它共有 n 项圬每一项为一 n 圱次多项式圬且第 i 项除在点 x i 为 y i f 在 x i 圩外圬在其余诸 x j 处均为地圮 直接对在圲圮圵圮圱圱圶圩编写程序虽然不是一个很好的算法圬但也是完全可行的圮读者可 作为一个编程练习来写出计算程序圮我们在这里将介绍一个更好的算法圬坎坥坶坩坬坬坥算 法圬这一算法的基本思想是从一个低次多项式出发圬逐次用递推方式加高并最 终得到结果圮下面给出这一算法的一个简短的说明圮对于给定的 x 圬令 P 0 11 是通 过在 x 1, y 1 圩点的唯一的地次多项式圬即 P 0 11 y 1 圬在同样定义 P 0 22, P 0 33, P 0 nn 圩圻定义 P 1 12 为 通过在 x 1, y 1 圩, 在 x 2, y 2 圩的唯一一次多项式圮同理有 P 1 23, P 1 34, P 1 n 1,n 圻 圬定义 P m i,i+m 为通 过在 x i, y i 圩, 在 x i+1, y i+1 圩, 在 x i+m, y i+m 圩的 m 次多项式圬则 P n 1 1n 的唯一 n 圱次多项式圬即为在圲圮圵圮圱圱圶圩的等价形式圮 即为通过在 x 1, y 1 圩, 在 x 2, y 2 圩,, 在 x n, y n 圩 所有这些 P 可以排成一个表圬 其坜祖代圢排在左边圬坜子代圢在右边并最终导致唯一的坜后代圢即为我们的解圮 如下式

49 圲圮圵圮插值和逼近圴圹所示在以 n 圵为例圩场 x 1 场 y 1 P11 0 P 1 12 x 2 场 y 2 P 0 22 P 2 13 P 1 23 P 3 14 x 3 场 y 3 P 0 33 P 2 24 P 4 15 在圲圮圵圮圱圱圷圩 x 5 场 y 5 P55 0 坜子代圢与坜父代圢之间满足如下递推关系场 P34 1 P25 3 x 4 场 y 4 P44 0 P35 2 P 1 45 P m i,i+m 在 x x i+m 圩 P m 1 i,i+m 1 圫在 x i x 圩 P m 1 i+1,i+m x i x i+m 在圲圮圵圮圱圱圸圩 练习证明上述递推关系直接用在圲圮圵圮圱圱圸圩计算是完全可行的圬但还可进一步改进圬一种改进的方法是在计算中做小量坜子代圢和坜父代圢之差的迭代圮即定义 C mi P m i,i+m P m 1 i,i+m 1 D mi P m i,i+m P m 1 i+1,i+m 在圲圮圵圮圱圱圹圩 从在圲圮圵圮圱圱圸圩可容易推出 C mi 和 D mi 所满足的递推关系为场 D m+1,i C m+1,i 在 x i+m+1 x 圩在 C m,i+1 D mi 圩 x i x i+m+1 在 x i x 圩在 C m,i+1 D mi 圩 x i x i+m+1 在圲圮圵圮圱圲地圩 在每一次递推中圬 C 和 D 是对插值的更高一次的修正圬最后的结果 P n 1 1n 系列从 y i 出发而最终到达坜家族圢终点的任一路径上的 C 及 D 圮 是任一 y i 加上一 分式有理函数插值和外推 有些函数在大部分函数圿圩不能很好地用多项式来近似圬我们在后面要举一个著名的例子圬但大部分函数都能用分式有理函数来近似圬这里边有一些深入的数学原理圬有兴趣的同学可以研读有关的数学书籍圮

50 圵地坃坈坁坐坔坅坒圲圮物理学中的常用数值方法 一个通过 m 圫圱个点在 x 0, y 0 圩, 在 x 1, y 1 圩,, 在 x m, y m 圩的分式有理函数可记为 R µ,ν 在 x 圩 P µ 在 x 圩 Q ν 在 x 圩 p 0 圫 p 1 x 圫 圫 p µ x µ 圱圫 q 1 x 圫 圫 q ν x ν 在圲圮圵圮圱圲圱圩 由于有 µ 圫 ν 圫圱个末知参数 p 及 q 圬因此圬必有 m 圫圱 µ 圫 ν 圫圱圮在做分式有理函数插值时圬还需要给定分子或分母的阶圮方程在圲圮圵圮圱圲圱圩两边乘以 Q ν 在 x 圩圬分别代入逐插值点圬得到一组线性代数方程组圬可由此解得 p i 和 q i 圬得到有理插值函数圮类似于多项式插值圬我们也可以构造一个递推算法圬这一算法称为坂坵坬坩坲坳坣坨坼坓坴坯坥坲算法圮这里不予介绍圬有兴趣的同学可参看在坊圮坓坴坯坥坲圬坒圮坂坵坬坩坲坳坣坨圬坉坮坴坲坯坤坵坣坴坩坯坮坴坯坎坵坭坥坲坩坣坡坬坁坮坡坬坹坳坹坳圬坓坰坲坩坮坧坥坲圭坖坥坲坬坡坧圬圱圹圹圲圬坃坨坡坰坴坥坲圲圩圮 三次样条插值 三次样条插值是目前应用最广泛的一种方法圬它可以保持插值函数具有二次连续导数圬 而且不会出现多项式插值中的一些问题圮三次样条插值问题可以表述如下圬对于给定的 n 个点在 x 1, y 1 圩, 在 x 2, y 2 圩,, 在 x n, y n 圩圬寻找一个分段三次多项式 S 在 x 圩 S 1 在 x 圩在 x 1 x x 2 圩 S 2 在 x 圩 在 x 2 x x 3 圩 S n 1 在 x 圩在 x n 1 x x n 圩 在圲圮圵圮圱圲圲圩 其中圬 S i 在 x 圩为一定义在区间坛 x i, x i+1 坝上的三次多项式并要求 S 在 x 圩在诸插值点上的一阶 及二阶导数连续圮 我们先来分析这一插值问题的解是否存在圬由于决定每个三次多项式需要圴个条件圬 所以总共需要圴在 n 圱圩个条件圬由每个多项式都要通过两个端点圬给出圲在 n 圱圩个条件圬 一阶导数和二阶导数连续给出圲在 n 圲圩个条件圬这样圬我们总共有圴在 n 圱圩 圲个条件圮 因此圬为了确定插值问题圬还需要补充二个条件圬由这两个条件的不同圬又有几种不同 的三次样条插值方法圬这里我们将介绍二种圬一种是再给定两个端点 x 1, x n 的一阶导数 值 y 1, y n 圬另一种是令两个端点上的二阶导数为地圬即令 y 1 y n 地圮 现在来求解这一问题圬由于二阶导数连续圬所以 y i, i 圲, 圳,, n 圱是存在的圬在任 一子区间坛 x i, x i+1 坝圬有 S i 在 x 圩 x x i x i+1 x i y i+1 圫 x i+1 x y i x i+1 x i 在圲圮圵圮圱圲圳圩

51 圲圮圵圮插值和逼近圵圱 对上式积分二次圬得到 S i 在 x 圩 3 圱在 x i+1 x 圩 y 圶 x i+1 x i 圫在 y i+1 圫在 y i i 圫 3 圱在 x x i 圩 y i+1 圶 x i+1 x i x x i 圩 x i+1 x i x i+1 x 圩 x i+1 x i 圱 y 2 i+1 在 x i+1 x i 圩圶 圱 y 2 i 在 x i+1 x i 圩圶 在圲圮圵圮圱圲圴圩 上式包含 n 个未知量 y i, i 圱, 圲, 圳,, n 圬为了决定这些未知量圬对上式微分一次圬得到 S i 圱在 x 圩 圲 在 x i+1 x 圩 2 x i+1 x i y i 圫 圫在 y i+1 在 y i 2 圱在在 x x i 圩 y 圲 x i+1 x i 圱 y 2 i+1 在 x i+1 x i 圩圶 圱 y 2 i 在 x i+1 x i 圩圶 i+1 圱圩 x i+1 x i 圱圩 x i+1 x i 在圲圮圵圮圱圲圵圩 由 S i 在 x 圩的连续性条件圬可得 n 圲个方程如下在对 i 圲, 圳,, n 圱圩 x i x i 1 y i 1 圫圶 x i+1 x i 1 y j 圫圳 x i+1 x i y i+1 圶 y i+1 y i x i+1 x i y i y i 1 x i x i 1 在圲圮圵圮圱圲圶圩 如前所述圬我们还需要圲个补充条件以决定 n 个未知量 y i 圬其选择方法是圬如果给定端 点的一阶导数圬则把 y 1 和 y n 代入方程在圲圮圵圮圱圲圵圩可得到两个关于 y 的方程圬与在圲圮圵圮圱圲圶圩一 起构成一个完整的线性代数方程组圻如果给定端点的二阶导数为地圬则在圲圮圵圮圱圲圶圩足以决 定 n 圲个未知量圮一旦求得了诸 y i 圬则对给定的 x 圬可利用在圲圮圵圮圱圲圴圩计算其对应的 y 值圮 下面我们看一个例子圬考虑下面的函数 y 圱圱圫圲圵 x 2 在圲圮圵圮圱圲圷圩 我们在坛 圱, 圱坝之间取圱圱个点圬分别用多项式圬分式有理函数和三次样条函数进行插值圬计算结果见图在圷圮圱圩圬其中实线为原函数圬点线为多项式插值的结果圬除了在插值点外圬与原函数的差别是相当大的圬点虚线和虚线分别为分式有理函数和三次样条插值的结果圬在图上几乎完全与准确值重合圮

52 圵圲坃坈坁坐坔坅坒圲圮物理学中的常用数值方法 坆坩坧坵坲坥圲圮圲场多项式圬分式有理函数和三次样条函数插值的比较 列表函数的积分 在物理上圬经常会遇到计算列表函数的积分圮实验测量的结果总是在一些分立点上圬大型数值计算的结果也只能是一些分立值圮样条插值方法可以用来计算此类积分圮对于列表函数圬可首先建立其样条插值圬求出由式在圲圮圵圮圱圲圲圩所指定的函数 S 在 x 圩圬由于诸 S i 在 x 圩都是三次多项式圬其积分可解析求出圬把每一分段上的积分值加起来就得到最终结果圮练习分别对于端点和内点计算 ( ) 的积分, 并最终给出样条插值方法计算列表函数的算法和程序 练习 对于函数 x 坣坯坳 x f 在 x 圩圱圫 x 3 2 (i), 用辛普生方法计算 f 在 x 圩 dx 1.5 (ii), 在区间坛圱, 圲坝上分 50 个子区间, 计算 f 在 x i 圩的值, 利用这些数值并用样条 函数方法 ( 见上题 ) 计算上述程分, 比较所得结果 Padé 插值与外推 在理论物理中圬微扰方法是进行实际计算的最重要的方法之一在在力学界称为摄动方法圩圬 粗略地讲圬微扰方法就是在计算某一物理量时圬利用这一物理问题中的某个小参量作级