2 目 录 2.4.1 Euler 方 法................................ 39 2.4.2 龙 格 库 塔 方 法............................ 40 2.4.3 两 点 边 值 问 题...........................

目 录 第 零 章 绪 论 5 第 一 章 计 算 机 数 和 误 差 9 1.1 计 算 机 数 及 其 表 示................................ 9 1.2 舍 入 误 差 对 计 算 的 影 响............................. 10 1.3 减 法 的 计 算.................................... 12 第 二 章 物 理 学 中 的 常 用 数 值 方 法 13 2.1 数 值 积 分 方 法.................................. 13 2.1.1 梯 形 法 和 辛 普 生 方 法.......................... 13 2.1.2 反 常 积 分 的 计 算............................. 15 2.1.3 高 斯 积 分 方 法.............................. 17 2.1.4 哥 西 主 值 的 计 算............................. 20 2.1.5 急 速 振 荡 函 数 的 积 分.......................... 22 2.1.6 高 维 积 分 的 计 算............................. 24 2.1.7 固 体 热 容 量 的 计 算........................... 24 2.2 方 程 的 求 根 和 最 优 化 方 法............................ 25 2.2.1 二 分 法.................................. 25 2.2.2 牛 顿 法 和 割 线 法............................. 26 2.2.3 一 维 最 优 化, 黄 金 分 割 法........................ 27 2.2.4 高 维 最 优 化, 共 扼 梯 度 方 法....................... 29 2.2.5 约 束 最 优 化............................... 33 2.2.6 最 小 二 乘 法 及 曲 线 拟 合......................... 34 2.3 函 数 的 求 值.................................... 36 2.3.1 多 项 式 的 求 值 和 级 数 求 和....................... 36 2.3.2 连 分 式 的 求 值.............................. 38 2.4 常 微 分 方 程 的 数 值 解 法............................. 39

2 目 录 2.4.1 Euler 方 法................................ 39 2.4.2 龙 格 库 塔 方 法............................ 40 2.4.3 两 点 边 值 问 题.............................. 44 2.5 插 值 和 逼 近.................................... 47 2.5.1 多 项 式 插 值............................... 48 2.5.2 分 式 有 理 函 数 插 值 和 外 推....................... 49 2.5.3 三 次 样 条 插 值.............................. 50 2.5.4 列 表 函 数 的 积 分............................. 51 2.5.5 Padé 插 值 与 外 推............................. 52 2.5.6 发 散 级 数 的 Cesáro 求 和 及 其 推 广................... 56 2.5.7 Borel 求 和................................ 58 2.5.8 Bézier 逼 近................................ 59 2.5.9 多 元 函 数 的 插 值 及 逼 近......................... 60 2.6 快 速 付 里 叶 变 换................................. 63 2.6.1 付 里 叶 变 换............................... 63 2.6.2 离 散 付 里 叶 变 换............................. 65 2.6.3 快 速 付 里 叶 变 换............................. 68 2.7 附 录 : 特 殊 函 数 的 计 算.............................. 71 2.7.1 误 差 函 数 的 计 算............................. 73 2.7.2 指 数 积 分 的 计 算............................. 74 2.7.3 整 数 阶 贝 塞 尔 函 数 及 虚 宗 量 贝 塞 尔 函 数 的 计 算........... 75 2.7.4 球 谐 函 数 的 计 算............................. 84 2.7.5 椭 圆 积 分 的 计 算............................. 86 2.8 附 录 : 高 斯 积 分 的 节 点 和 权 重.......................... 88 第 三 章 数 值 线 性 代 数 93 3.1 主 元 消 去 法 解 线 性 代 数 方 程 组......................... 93 3.2 LU 分 解 法.................................... 94 3.3 三 对 角 方 程 组 的 求 解.............................. 97 3.4 实 对 称 矩 阵 的 本 征 值 和 本 征 向 量........................ 98 3.4.1 Householder 方 法............................ 98 3.4.2 QL 算 法................................. 100 第 四 章 电 磁 场 的 计 算 103 4.1 Maxwell 方 程, 边 值 问 题............................. 103 4.2 差 分 和 差 商.................................... 105

目 录 3 4.3 二 维 Possion 方 程 的 五 点 差 分 格 式....................... 106 4.4 差 分 方 程 的 求 解................................. 109 4.4.1 简 单 迭 代 法 Jaobi 方 法....................... 110 4.4.2 逐 次 迭 代 法 Gauss Seidal 迭 代 方 法................. 110 4.4.3 逐 次 超 松 弛 迭 代 法 (Suessive Over-Relaxation)........... 110 4.5 变 分 方 法 复 习.................................. 112 4.6 有 限 元 方 法 的 理 论 基 础............................. 115 4.7 用 有 限 元 方 法 求 解 二 维 Laplae 方 程...................... 118 第 五 章 Monte Carlo 方 法 125 5.1 随 机 变 量 及 其 分 布................................ 125 5.2 赝 随 机 数 的 产 生................................. 127 5.3 用 Monte Carlo 方 法 计 算 积 分.......................... 129 5.4 自 相 似 性 与 分 形................................. 132 5.5 扩 散 限 制 粘 结 的 计 算 机 模 拟.......................... 134 5.6 高 分 子 的 模 拟 和 无 规 行 走............................ 135 第 六 章 逾 渗 和 统 计 物 理 问 题 137 6.1 逾 渗 简 介..................................... 137 6.2 逾 渗 的 计 算 机 模 拟 和 分 析............................ 140 6.3 正 则 系 综 和 统 计 模 型.............................. 141 6.3.1 Ising 模 型................................ 142 6.3.2 Lenard-Jones 模 型............................ 143 6.4 统 计 物 理 的 Monte Carlo 模 拟, Metroplis 方 法................. 143 6.5 Ising 模 型 的 Monte Carlo 模 拟......................... 146 6.6 经 典 统 计 物 理 复 习................................ 147 6.6.1 The miro anonial ensemble(nve).................. 148 6.6.2 The anonial ensemble(nvt)..................... 148 6.6.3 The NPT ensemble(npt)........................ 149 6.6.4 The grand anonial ensemble(..................... 149 6.6.5 Thermodynamial relations....................... 150 6.6.6 Correlation funtions........................... 150 6.6.7 Time orrelation funtion and transport oeffiients.......... 152 6.7 液 体 模 型 的 Monte Carlo 模 拟.......................... 153 6.7.1 Monte Carlo Simulation of The Lenard -Jones Model.......... 154 6.7.2 自 由 能 计 算............................... 156

4 目 录 6.7.3 王 -Landau 方 法 及 其 它......................... 157 6.8 分 子 动 力 学 方 法................................. 159 6.8.1 General proedure of MD (NVE ensemble).............. 159 6.8.2 Simulation of Lennard-Jones liquids.................. 163 6.8.3 Simulation of Hard-sphere systems................... 163 6.9 Simulation of Langevin dynamis........................ 164 6.10 Long range fore and Ewald summation..................... 166 第 七 章 原 子 结 构 的 计 算 167 7.1 原 子 结 构 问 题.................................. 167 7.2 变 分 法...................................... 168 7.3 Viral 定 理 和 Hellmann-Feynman 定 理..................... 170 7.4 轨 道 近 似 和 Hartree Fok 方 程......................... 173 7.5 统 计 近 似 和 X α 方 程............................... 175 7.6 径 向 方 程 的 求 解 方 法.............................. 178 7.7 计 算 程 序..................................... 181

第 零 章 绪 论 不 论 是 理 论 物 理, 还 是 实 验 物 理 都 离 不 开 数 值 计 算, 例 如, 量 子 力 学 中 简 谐 振 子 的 波 函 数 可 以 用 厄 米 多 项 式 来 表 示, 但 为 了 得 到 波 函 数 的 数 值, 仍 然 需 要 作 数 值 计 算 ; 实 验 上 测 量 到 的 数 据 要 进 行 数 据 处 理, 也 涉 及 到 数 据 拟 合, 分 析 并 计 算 误 差 等 一 系 列 数 值 计 算. 大 家 在 过 去 的 学 习 中 已 经 知 道, 牛 顿 力 学 方 程 只 有 二 体 问 题 是 可 解 的 ( 自 由 粒 子 或 小 振 动 等 特 殊 问 题 除 外 ), 三 体 以 上 的 问 题 曾 折 磨 了 全 世 界 的 许 多 优 秀 的 数 学 家 和 理 论 物 理 学 家, 但 最 终 也 未 得 到 解 析 解 ; 麦 克 斯 韦 方 程 组 只 有 在 一 些 非 常 特 殊 的 条 件 下 才 能 求 得 解 析 解 ; 量 子 力 学 的 薛 定 格 方 程, 除 了 氢 原 子 和 简 谐 振 子 外, 还 没 有 一 个 真 实 的 物 理 问 题 可 以 找 到 解 析 解 ; 统 计 物 理 告 诉 我 们, 只 要 求 得 了 系 统 的 配 分 函 数, 则 一 切 物 理 量 都 可 通 过 求 微 分 得 到, 然 而, 除 了 可 数 的 几 个 二 维 模 型 外, 没 有 一 个 实 际 的 物 理 系 统 的 配 分 函 数 被 求 出 来 过. 一 些 大 规 模 的 物 理 实 验, 耗 资 在 数 亿 美 元, 如 果 没 有 充 分 的 论 证 就 去 做 的 话, 其 浪 费 是 十 分 巨 大 的. 对 于 这 样 一 类 问 题, 大 规 模 的 数 值 计 算 往 往 可 以 发 挥 重 要 作 用. 大 规 模 数 值 计 算 也 提 供 了 理 论 物 理 和 实 验 物 理 之 间 的 桥 梁. 几 乎 所 有 的 实 际 物 理 问 题 都 无 法 得 到 精 确 解, 发 展 各 种 理 论 模 型 和 近 似 方 法 就 成 为 理 论 物 理 学 的 重 要 方 面. 理 论 模 型 的 正 确 性 最 终 需 要 经 受 实 践 的 检 验, 但 对 于 一 个 确 定 的 物 理 模 型, 其 对 应 的 实 验 体 系 往 往 很 难 获 得, 为 了 检 验 针 对 一 个 确 定 模 型 的 近 似 方 法, 最 好 能 够 得 到 这 个 模 型 的 精 确 解, 而 数 值 模 拟 方 法 正 好 能 够 提 供 这 样 的 结 果. 大 量 的 实 验 体 系 往 往 包 含 多 种 因 素, 这 些 因 素 有 些 可 以 排 除, 而 有 些 非 常 难 以 排 除. 如 果 一 个 理 论 的 出 发 点 包 含 了 所 有 的 实 验 因 素, 则 这 样 的 理 论 几 乎 肯 定 得 不 到 解 析 结 果. 但 是, 数 值 模 拟 方 法 一 般 总 能 够 把 各 种 因 素 包 含 进 去, 从 而 能 够 更 好 的 模 拟 实 验 环 境. 需 要 指 出 的 是, 长 期 以 来, 我 们 总 有 一 些 误 解, 一 讲 到 计 算, 似 乎 就 需 要 超 级 计 算 机. 事 实 上, 现 在 使 用 的 笔 记 本 计 算 机, 各 种 微 型 计 算 机 的 计 算 能 力 往 往 比 10 年 前 的 超 级 计 算 机 的 计 算 能 力 更 强. 用 一 台 目 前 似 乎 已 经 被 人 看 不 起 的 IBM-PC/XT 微 机 (IBM-PC 是 基 于 intel 芯 片 的 第 一 代 个 人 计 算 机, 大 约 1980 年 推 出, 其 CPU 是 intel-8086 芯 片, 有 64k 内 存, 两 个 5 吋 软 驱, 机 器 带 有 MS-DOS 1.0 和 CP/M 两 个 操 作 系 统,IBM-PC/XT 是 IBM-PC 的 扩 展, 增 加 了 一 个 10M 的 硬 盘, 并 且 把 内 存 增 加 到 256K), 就 可 以 很 容 易 地 解 决 三 体 问 题, 四 体 问 题,. 所 有 元 素 周 期 表 上 的 原 子 的 能 级 和 波 函 数 也 可 以 相 当 精 确

6 第 零 章 绪 论 地 计 算 出 来. 如 果 你 拥 有 一 台 目 前 市 场 上 一 般 配 置 的 微 机, 你 就 可 以 模 拟 一 些 二 维 和 三 维 模 型 系 统 的 配 分 函 数, 如 Lenard-Jones 流 体, 三 维 Ising 模 型 等, 可 以 计 算 元 素 固 体 和 简 单 化 合 物 晶 体 的 电 子 结 构, 声 子 谱 等. 进 一 步, 你 如 果 拥 有 或 能 够 使 用 到 大 型 计 算 机, 或 利 用 微 机 构 成 计 算 集 群, 你 就 可 以 作 真 实 流 体 的 流 体 力 学 计 算, 可 以 计 算 复 杂 边 界 条 件 下 的 三 维 电 磁 场, 可 以 计 算 三 维 真 实 系 统 的 配 分 函 数,. 同 样 你 也 能 够 在 大 型 物 理 实 验 开 始 前 作 很 多 模 拟, 对 实 验 的 各 种 情 况 有 一 个 清 楚 的 概 念. 本 课 程 的 目 的 在 于 对 物 理 中 的 数 值 计 算 进 行 一 些 入 门 介 绍, 使 大 家 在 学 完 本 课 程 后, 在 组 织 一 些 较 大 规 模 的 计 算 时 心 中 有 数, 少 走 弯 路. 课 程 分 两 部 分, 一 部 分 是 基 本 算 法, 主 要 介 绍 计 算 机 数 的 特 点, 舍 入 误 差, 计 算 稳 定 性 等 基 本 概 念 和 一 些 常 用 的 数 值 计 算 方 法 ; 第 二 部 分 则 结 合 物 理 问 题, 讲 一, 二 个 专 题, 使 大 家 学 到 一 些 组 织 较 大 规 模 计 算 的 步 骤 和 技 巧. 学 习 本 课 程 要 求 学 过 四 大 力 学 ( 经 典 力 学, 电 动 力 学, 热 力 学 和 统 计 物 理 学 及 量 子 力 学 ) 和 数 学 物 理 方 法, 掌 握 至 少 一 门 计 算 机 高 级 语 言 ( 如 C, PASCAL, FORTRAN, Java 等 ). 计 算 是 一 门 实 践 性 很 强 的 课 程, 上 机 实 习 是 本 课 程 的 一 个 有 机 组 成 部 分. 关 于 本 讲 义 的 一 点 说 明 : 这 是 笔 者 在 1993 年 为 开 设 计 算 物 理 课 程 准 备 的 讲 义. 1992 年, 上 海 交 通 大 学 应 用 物 理 系 代 理 系 主 任 张 仲 渊 教 授 找 到 我, 告 诉 我 决 定 在 本 科 生 高 年 级 开 设 计 算 物 理 课 程, 并 安 排 由 笔 者 开 课. 在 接 到 这 一 任 务 后, 笔 者 翻 阅 了 当 时 国 内 出 版 的 几 本 计 算 物 理 教 材, 发 现 都 过 于 专 门, 不 适 于 做 为 本 科 高 年 级 的 教 材, 而 更 适 合 于 某 一 专 业 的 研 究 生 教 材. 为 此, 我 只 好 试 图 自 己 写 一 本 讲 义, 以 应 急 需. 基 本 算 法 部 分 主 要 参 考 了 W. H. Press, S. A. Teukolsky, W. T. Vetterling and B. P. Flannery 所 著 Numerial Reipes, The Art of Sientifi Computing 一 书, 当 时, 这 本 书 出 版 时 间 还 不 算 太 长, 但 已 经 成 为 大 多 数 经 常 从 事 数 值 计 算 的 物 理 学 工 作 者 的 手 头 必 备 书, 我 自 己 当 时 也 买 了 一 本, 书 价 是 我 的 两 个 月 的 工 资 这 本 书 一 方 面 讲 解 算 法, 另 一 方 面 对 每 一 个 算 法 给 出 了 透 明, 简 洁 的 源 程 序 尽 管 这 些 程 序 的 效 率, 安 全 性 方 面 比 一 些 成 熟 的 软 件 包 要 差, 但 由 于 其 简 洁 清 楚, 特 别 适 合 于 修 改 后 放 到 自 己 的 程 序 中 去 讲 义 的 内 容 曾 在 上 海 交 通 大 学 应 用 物 理 系 讲 授 多 次, 各 个 年 级 的 同 学 积 极 参 与 了 本 课 程 的 学 习 并 对 课 程 的 讲 法 提 出 了 很 多 非 常 宝 贵 的 意 见 和 建 议, 这 些 建 议 在 讲 义 的 成 文 中 对 笔 者 有 很 大 帮 助, 在 此 表 示 感 谢. 2000 年 后, 物 理 系 决 定 把 计 算 物 理 课 程 改 为 双 语 课 程, 这 本 讲 义 就 算 完 成 了 其 历 史 使 命. 在 教 学 过 程 中, 我 也 发 现 这 一 讲 义 对 一 些 问 题 的 讲 解 过 于 简 洁, 而 内 容 的 深 度 也 不 太 适 合 本 科 生 的 学 习. 这 一 段 时 间, 国 际 上 出 现 了 几 本 比 较 合 适 的 教 材, 所 以,

7 当 时 认 为 这 一 讲 义 不 需 要 继 续 修 订, 也 就 放 在 一 边 了. 随 后, 研 究 生 的 课 程 做 了 一 些 改 革, 增 加 了 36 学 时 的 数 值 计 算, 任 课 老 师 感 觉 这 个 讲 义 的 内 容 作 为 研 究 生 数 值 计 算 的 课 程 比 较 合 适, 遂 采 用 这 本 讲 义 的 电 子 版 教 学. 2013 年 10 月, 致 远 学 院 的 叶 曦 副 院 长 和 负 责 物 理 班 教 学 的 郑 杭 教 授 找 我, 希 望 我 为 致 远 学 院 的 同 学 开 设 一 门 物 理 中 数 值 计 算 的 课 程 注 意 到 致 远 学 院 已 经 有 至 少 两 门 相 关 课 程, 这 门 课 程 应 该 更 多 强 调 物 理 方 面 为 了 这 个 课 程, 利 用 寒 假, 我 修 订 了 讲 义, 主 要 是 对 原 来 的 一 些 表 述 做 了 修 改 和 补 充, 同 时 尽 可 能 强 调 物 理 方 面 的 内 容. 鉴 于 目 前 已 经 存 在 大 量 优 秀 的 数 学 软 件, 如 Matlab, Maple 等. 大 量 的 小 规 模 的 计 算 往 往 可 以 通 过 这 些 软 件 方 便 进 行, 因 此, 此 次 修 订 时 时 也 适 当 提 到 这 方 面 的 内 容 讲 义 中 提 到 的 一 些 程 序 将 放 在 课 程 主 页 上, 这 些 程 序 中, 有 些 来 自 前 面 提 到 的 Numerial Reipes, 有 些 是 我 曾 经 多 次 使 用 过 的, 也 有 些 是 专 为 此 讲 义 而 准 备 的, 由 于 时 间 关 系, 没 有 经 过 仔 细 调 试 和 考 验, 因 此 在 用 于 实 际 问 题 时, 一 定 要 经 过 仔 细 试 算, 以 免 造 成 不 必 要 的 损 失, 同 时, 我 也 愿 在 此 声 明, 笔 者 对 因 使 用 此 讲 义 中 的 程 序 而 导 致 的 任 何 直 接 的 和 间 接 的 损 失 将 不 负 担 任 何 责 任. 所 有 程 序 都 以 FORTRAN 77 写 出, 鉴 于 C 语 言 在 实 际 工 作 中 使 用 的 越 来 越 广 泛, 同 时 也 是 很 多 同 学 喜 欢 的 语 言, 部 分 程 序 已 经 翻 译 为 C 语 言. 由 于 程 序 是 用 FORTRAN 写 的, 翻 译 中 自 然 留 有 FORTRAN 的 痕 迹. 在 上 次 讲 义 完 成 到 这 次 修 订 之 间, 面 向 对 象 的 程 序 设 计 逐 步 走 向 主 流, 其 代 表 性 语 言 是 ++ 和 Java, 建 议 有 志 于 从 事 计 算 物 理 的 同 学 注 意 这 一 点.

8 第 零 章 绪 论

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

10 第 一 章 计 算 机 数 和 误 差 1.2 舍 入 误 差 对 计 算 的 影 响 在 数 值 计 算 过 程 中, 一 般 每 一 步 都 有 舍 入 误 差, 舍 入 误 差 的 积 累 有 可 能 使 整 个 计 算 失 败. 例 如, 考 虑 下 述 问 题, 计 算 积 分 利 用 分 部 积 分 法 得 到 或 1 E n = 1 0 x n e x 1 dx, n = 1, 2,..., 20. x n e x 1 dx = x n e x 1 1 0 n 1 0 0 x n 1 e x 1 dx, E n = 1 n E n 1 容 易 算 出 E 1 = 1/e. 利 用 上 述 递 推 公 式, 在 IBM-PC 微 机 上 用 双 精 度 计 算 结 果 如 下 : E 1 = 0.36787944117144 E 2 = 0.26424111765712 E 3 = 0.20727664702865 E 4 = 0.17089341188538 E 5 = 0.14553294057308 E 6 = 0.12680235656152 E 7 = 0.11238350406936 E 8 = 0.10093196744509 E 9 = 0.091612292994171 E 10 = 0.083877070058293 E 11 = 0.077352229358780 E 12 = 0.071773247694637 E 13 = 0.066947779969723 E 14 = 0.062731080423873 E 15 = 0.059033793641902 E 16 = 0.055459301729570 E 17 = 0.057191870597308 E 18 =.029453670751536 E 19 = 1.5596197442792 E 20 = 30.192394885584 被 积 函 数 在 积 分 区 间 [0, 1] 的 内 部 总 是 正 的, 因 此 积 分 值 不 可 能 为 负, 更 进 一 步, 由 于 E n = 1 0 x n e x 1 dx 1 0 x n dx = 1 n + 1 所 以, 必 有 E 20 1/21. 显 然, 上 述 计 算 结 果 是 不 正 确 的, 而 问 题 就 出 在 舍 入 误 差 上 面. 我 们 对 此 作 一 个 简 单 的 分 析. 双 精 度 数 可 表 示 14 位 10 进 制 数, 所 以 E 1 的 误 差 约 为 ϵ = 10 15 ( 也 可 能 大 于 此 数, 这 里 仅 作 一 数 量 级 的 估 计 ). 在 递 推 过 程 中, 这 个 误 差 将 被 传 播, 放 大, 最 后 导 致 计 算 结 果 的 错 误. 误 差 的 传 播 规 律 为 E 2 = 1 2(E 1 + ϵ) = 1 2E 1 2ϵ E 3 = 1 3(1 2E 1 2ϵ) = 1 3(1 2E 1 ) + 3!ϵ

1.2 舍 入 误 差 对 计 算 的 影 响 11 计 算 到 E 15 时 的 误 差 为 15!ϵ 10 3, 而 计 算 到 E 18 时 的 误 差 已 为 18!ϵ 6.0, 已 经 远 大 于 E 18 的 上 限 1/19. 现 在 我 们 考 虑 把 递 推 关 系 改 写 成 E n 1 = 1 E n, n =, 4, 3, 2. n 利 用 关 系 E n 1 n + 1 取 E 30 = 0.0 开 始 计 算, 此 时 误 差 为 ϵ 1/31 0.032, 作 一 次 递 推 后 得 到 E 29, 其 误 差 减 小 到 1 1 30 31 0.0011 在 计 算 E 20 时, 误 差 已 减 小 到 20! 31! 3 10 16 准 确 到 14 位 有 效 数 字, 初 始 误 差 的 影 响 已 完 全 消 除 了. 而 每 一 步 的 舍 入 误 差, 都 在 其 后 的 迭 代 中 减 小 和 消 除 下 表 是 计 算 结 果. E 30 = 0.0 E 29 = 0.33333333333333 E 28 = 0.033333333333333 E 27 = 0.034523809523810 E 26 = 0.035758377425044 E 25 = 0.037086216252883 E 24 = 0.038516551349885 E 23 = 0.040061810360421 E 22 = 0.041736443027808 E 21 = 0.043557434407827 E 20 = 0.045544884075818 E 19 = 0.047722755796209 E 18 = 0.050119854958094 E 17 = 0.052771119168995 E 16 = 0.055719345931236 E 15 = 0.059017540879298 E 14 = 0.062732163941380 E 13 = 0.066947702575616 E 12 = 0.071773253648030 E 11 = 0.077352228862664 E 10 = 0.083877070103394 E 9 = 0.091612292989661 E 8 = 0.10093196744559 E 7 = 0.11238350406930 E 6 = 0.12680235656153 E 5 = 0.14553294057308 E 4 = 0.17089341188538 E 3 = 0.20727664702865 E 2 = 0.26424111765712 E 1 = 0.36787944117144 由 上 例 可 以 看 出 舍 入 误 差 对 计 算 的 影 响, 如 果 我 们 在 计 算 中 不 去 分 析 算 法, 即 使 编 出 了 正 确 的 程 序, 也 可 能 得 出 错 误 的, 有 时 甚 至 是 荒 谬 的 结 果.

12 第 一 章 计 算 机 数 和 误 差 1.3 减 法 的 计 算 也 许, 大 家 认 为 减 法 是 一 种 非 常 简 单 的 计 算, 但 是 在 用 计 算 机 做 减 法 时, 如 不 小 心, 就 很 容 易 出 错, 问 题 就 出 在 舍 入 误 差 上 面. 我 们 考 虑 下 面 的 例 子 : 对 于 大 的 x, 计 算 x + 1 x (1.1) 取 x = 162835.0, 在 一 个 具 有 6 位 十 进 制 有 效 数 字 的 计 算 机 上, 计 算 结 果 为 0.001, 而 准 确 到 6 位 十 进 制 有 效 数 字 的 精 确 结 果 为 0.123907 10 2, 计 算 中 损 失 了 5 位 有 效 数 字. 对 于 更 大 的 x, 则 可 能 根 本 得 不 到 正 确 的 结 果. 但 是, 如 果 把 计 算 公 式 改 为 ( x + 1 x + 1 + x 1 x) = (1.2) x + 1 + x x + 1 + x 就 可 以 得 到 较 好 的 结 果. 又 如 对 于 大 的 N, 计 算 N+1 N dx = artan(n + 1) artan(n) 1 + x2 时, 可 以 使 用 等 价 的 公 式 : ( ) 1 artan(n + 1) artan(n) = artan 1 + N(N + 1) 同 理, 对 于 小 的 ϵ ( sin (x + ϵ) sin(x) = 2 os x + ϵ ) ( ϵ sin 2 2) 上 面 分 析 的 是 一 些 简 单 的 例 子, 在 实 际 问 题 的 计 算 中, 常 常 会 碰 到 一 些 十 分 复 杂 的 问 题, 即 使 有 经 验 的 计 算 研 究 人 员 也 经 常 在 减 法 计 算 中 犯 错 误. 练 习 : 分 析 下 面 的 表 达 式, 给 出 对 任 何 x( 只 要 在 所 给 函 数 定 义 域 内 ) 都 能 得 出 正 确 结 果 的 计 算 公 式. os 2 x sin 2 x 1 + x2 1 x 2 ln(x) 1 e 2x e x 1 v2 1 x2 sin(x) x tg(x)

第 二 章 物 理 学 中 的 常 用 数 值 方 法 在 这 一 章 和 随 后 的 几 章 中, 我 们 将 介 绍 几 种 常 用 算 法, 这 些 内 容 应 该 在 数 学 课 中 讨 论, 但 目 前 的 数 学 课 程 并 不 包 含 这 些 属 于 计 算 数 学 课 程 的 内 容, 所 以 我 们 先 作 一 些 介 绍. 2.1 数 值 积 分 方 法 积 分 分 为 定 积 分 和 不 定 积 分 两 种, 我 们 在 这 里 只 讨 论 定 积 分, 不 定 积 分 的 计 算 可 作 为 变 上 限 的 定 积 分 来 计 算. 考 虑 定 积 分 I = 它 代 表 f(x) 曲 线 下 的 面 积. b a f(x)dx (2.1) 2.1.1 梯 形 法 和 辛 普 生 方 法 在 微 积 分 中, 我 们 知 道 积 分 是 由 求 和 的 极 限 定 义 的, 这 就 是 第 一 种 计 算 积 分 的 方 法 梯 形 法. 把 区 间 [a, b] 分 成 一 些 子 区 间 [a = x 0 < x 1 < x 2 < < x n = b], 则 (2.1) 的 积 分 I 可 近 似 为 : I = 1 n 1 (x i+1 x i )(f(x i ) + f(x i+1 )) (2.2) 2 i=0 在 实 际 计 算 中, 常 采 用 等 间 距 分 法, 即 x i+1 x i = h = b a n 为 一 常 数, 上 面 公 式 可 简 化 为 : I = h n 1 (f(x i ) + f(x i+1 )) 2 i=0 为 了 节 省 工 作 量, 实 际 使 用 的 计 算 公 式 为 : n 1 I = h f(x i ) + h 2 (f(x 0) + f(x n )) (2.3) i=1 当 子 区 间 越 来 越 小, h 0 时, 求 和 的 结 果 将 趋 于 积 分 值 由 于 我 们 事 先 并 不 知 道 积 分 的 结 果, 所 以 我 们 也 就 不 知 道 求 和 结 果 和 积 分 结 果 之 间 的 差 别 为 此, 我 们 需 要 估 计 计

14 第 二 章 物 理 学 中 的 常 用 数 值 方 法 算 误 差, 一 个 直 观 的 估 计 是 比 较 不 同 大 小 的 子 区 间 的 结 果 之 差 为 了 在 每 次 子 区 间 变 小 时 能 够 用 到 前 次 的 计 算 结 果, 我 们 可 以 从 n = 1 开 始 计 算, 每 次 把 区 间 数 增 加 一 倍, 则 前 一 次 计 算 的 函 数 结 果 仍 然 有 用, 通 过 比 较 连 续 两 次 的 计 算 结 果 之 差 是 否 小 于 给 定 的 误 差 限 来 结 束 计 算. 梯 形 方 法 实 际 上 把 每 一 个 子 区 间 上 的 函 数 近 似 为 线 性 函 数 αx + β, 然 后 计 算 其 积 分 值, 再 把 结 果 加 起 来.( 作 为 一 个 练 习, 请 自 行 证 明 ), 因 此, 我 们 说 梯 形 方 法 是 具 有 线 性 代 数 精 度 的 方 法. 下 面 考 虑 具 有 二 次 代 数 精 度 的 方 法, 辛 普 生 方 法. 先 考 虑 把 积 分 区 间 [a, b] 分 为 二 个 子 区 间, 分 点 为 (a, a + h, a + 2h = b), h = b a 2, 在 区 间 上 把 被 积 函 数 用 二 次 函 数 来 近 似, 即 取 f(x) α(x a) 2 + β(x a) + γ, 其 中 α, β, γ 可 利 用 分 点 上 的 函 数 值 f(a), f(a + h), f(b) 表 示 为 : 简 单 地 计 算 得 到 : b a α = β = f(b) 2f(a + h) + f(a) 2h 2 f(b) + 4f(a + h) 3f(a) 2h γ = f(a) (α(x a) 2 + β(x a) + γ) dx = h (f(b) + 4f(a + h) + f(a)) (2.4) 3 这 就 是 著 名 的 辛 普 生 求 积 公 式. 显 然, 它 具 有 二 次 代 数 精 度. 实 际 计 算 中, 把 被 积 区 间 [a, b] 先 分 为 N 个 子 区 间, 然 后 在 每 个 子 区 间 中 用 辛 普 生 求 积 公 式 计 算, 并 把 最 后 结 果 加 起 来. 其 计 算 公 式 为 : b a f(x)dx h (f(a) + 4f(a + h) + 2f(a + 2h) + 4f(a + 3h) + + f(b)) (2.5) 3 这 是 一 个 在 实 际 计 算 中 常 用 的 方 法, 它 对 于 一 般 的 有 限 区 间 的 积 分 都 能 给 出 较 好 的 正 确 结 果. 辛 普 生 方 法 有 各 种 不 同 的 实 现, 我 们 在 这 里 给 出 简 单 的 一 种, 取 每 个 子 区 间 均 相 等, 从 二 个 子 区 间 开 始, 每 次 子 区 间 数 目 加 倍, 这 样 可 以 在 每 次 计 算 时 使 用 前 次 的 计 算 值. 这 种 实 现 的 缺 点 是 在 整 个 积 分 区 间 均 匀 取 点, 对 于 在 部 分 区 间 变 化 剧 烈 而 在 其 余 地 方 变 化 平 缓 的 被 积 函 数 计 算 效 率 不 高. 这 一 实 现 可 以 充 分 利 用 它 与 梯 形 法 的 关 系. 假 定 n = 2 k, 用 S k 表 示 用 2 k 个 子 区 间 梯 形 公 式 求 得 的 结 果. 则 式 (2.5) 可 写 为 : b a f(x)dx 1 3 [4h(f(a + h) + f(a + 2h) + + f(b h) + 1 2 f(a) + 1 2 f(b)) 2h(f(a + 2h) + f(a + 4h) + + f(b 2h) + 1 2 f(a) + 1 2 f(b)) + h(f(a) + f(b)) 2h(f(a) + f(b)) + h(f(a) + f(b))] = 1 3 (4S k S k 1 ) (2.6)

2.1 数 值 积 分 方 法 15 计 算 定 积 分 还 有 很 多 方 法, 如 龙 贝 格 方 法, Newton-Cotes 方 法 等, 我 们 将 不 在 此 讨 论, 有 兴 趣 的 同 学 可 参 看 有 关 计 算 数 学 的 书 籍. Quadpak(netlib.org) 中 包 含 了 若 干 非 常 稳 定 的 求 积 分 的 程 序, 建 议 在 实 际 工 作 中 使 用 在 Matlab 中, 梯 形 法, 辛 普 森 均 作 为 标 准 函 数 提 供, 实 际 应 用 时, 只 要 写 一 个 被 积 函 数 的 m 程 序, 调 用 积 分 程 序 即 可 常 用 的 几 个 求 积 分 的 函 数 是 quad, quadl, quadgk 等 其 中 quad 使 用 自 适 应 的 辛 普 森 算 法,quadl 应 用 龙 贝 格 加 速 来 改 善 计 算 效 率 和 提 高 精 度 例 如, 计 算 积 分 1 先 定 义 函 数 fun 并 存 入 fun.m funtion y=fun(x) end y=sin(x.^2)./(1.+x.^4); 然 后 在 Matlab 的 命 令 窗 口 运 行 0 sin(x 2 ) 1 + x 4 dx >> format long >> quad(@fun,0,1) 结 果 为 : ans = 0.229253354700751 2.1.2 反 常 积 分 的 计 算 反 常 积 分 分 为 两 类, 一 类 是 积 分 区 间 有 限, 在 积 分 区 间 内 被 积 函 数 有 奇 点, 另 一 类 是 积 分 区 间 为 无 限. 对 于 二 者 兼 而 有 之 的 积 分, 可 以 分 为 两 个 或 多 个 积 分 来 处 理, 使 每 个 积 分 为 上 述 两 类 之 一. 积 分 区 间 内 含 有 奇 点 的 积 分 对 于 积 分 区 间 内 含 有 奇 点 的 积 分, 根 据 奇 点 的 性 质, 可 采 用 下 面 的 方 法 处 理. 可 去 奇 点, 设 x 0 为 一 可 去 奇 点, 且 知 道 在 x 0 点 被 积 函 数 的 极 限, 则 只 要 在 计 算 函 数 时, 在 x 0 的 一 个 邻 域 内 用 极 限 取 代 函 数 值 即 可. 例 如, sin(x) 在 x = 0 有 可 去 奇 点, 但 当 x x 0 时, sin(x) x 1 x2 6 +, 在 计 算 函 数 时, 可 使 用 下 面 的 语 句 :

16 第 二 章 物 理 学 中 的 常 用 数 值 方 法 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 = lim r x 0 b r f(x)dx 定 义 一 个 收 敛 于 x 0 的 序 列 b > r 1 > r 2 >, 例 如 可 取 r n = x 0 + 2 n, 则 上 述 积 分 可 写 为 一 个 求 和 : b b r1 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 ϵ 时, 终 止 计 算. 如 果 上 式 不 收 敛, 则 很 可 能 原 积 分 不 存 在. 消 除 奇 点, 有 些 奇 点 可 以 通 过 变 量 替 换 或 对 被 积 函 数 进 行 变 换 而 消 去. 我 们 看 几 个 r2 r3 例 子 : 1 作 变 量 替 换 x = t 2, dx = 2t dt, 则 有 : 或 把 积 分 写 为 : 1 0 os(x) x dx = 1 1 第 二 项 积 分 中 的 奇 点 已 成 为 可 去 奇 点. 0 0 0 os(x) x dx os(x) x dx = 2 1 0 os(t 2 )dt 1 1 os(x) 1 1 dx + dx = 2 + x 0 x 0 I(x) = x 0 e t 1 t 在 t = 1 的 邻 域, 被 积 函 数 与 e 1 /(1 t) 很 相 像, 因 此 变 为 : I(x) = x 0 e 1 x 1 t dt + 0 所 有 的 奇 异 性 都 在 第 一 项, 而 这 一 项 可 解 析 处 理. os(x) 1 dx x e t e 1 dt = 1 x 1 t e ln 1 x + e t e 1 dt 1 t 除 上 面 例 子 演 示 的 方 法 外, 还 可 用 分 部 积 分 法 消 除 奇 点, 这 里 不 再 一 一 举 例 了. 0

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

18 第 二 章 物 理 学 中 的 常 用 数 值 方 法 有 限 区 间 的 高 斯 型 积 分 公 式 : 高 斯 勒 让 德 积 分 公 式 : 1 1 f(x)dx = n w i f(x i ) (2.10) 此 处, x i 为 n 阶 勒 让 德 多 项 式 的 根, 通 常 称 为 积 分 节 点, w i 称 为 积 分 权 重. 等 权 重 切 贝 雪 夫 积 分 公 式 : 这 里 x i 为 积 分 节 点, 权 重 为 2/n. 1 第 一 类 高 斯 切 贝 雪 夫 积 分 公 式 : 1 1 1 f(x)dx = 2 n f(x) 1 x 2 dx = π n n f i=1 i=1 n f(x i ) (2.11) i=1 ( os ( )) (2i + 1)π 2n (2.12) 第 二 类 高 斯 切 贝 雪 夫 积 分 公 式 : 1 1 f(x) 1 x 2 dx = π n + 1 n ( ) ( ( )) iπ iπ sin 2 f os n + 1 n + 1 i=1 如 果 积 分 区 间 不 是 [ 1, 1], 而 是 [a, b], 可 利 用 变 换 ( ) ( ) b a b + a x = t + 2 2 (2.13) (2.14) 将 积 分 区 间 变 为 [ 1, 1], 然 后 利 用 上 述 公 式 进 行 计 算. 在 实 际 计 算 时, 常 把 积 分 区 间 分 为 几 个 小 区 间, 在 每 一 个 小 区 间 内 使 用 高 斯 积 分 公 式, 并 把 结 果 加 起 来. 练 习 : 试 编 写 上 述 四 种 高 斯 积 分 公 式 的 程 序, 并 利 用 所 编 写 的 程 序 计 算 积 分 1 1 x + 1.5dx, 1 1 1 2 + x dx, 1 0 dx x 2 0 e os2 (x) dx, ( 积 分 所 需 的 节 点 及 权 重 在 后 面 的 表 中 给 出 ). 1 0 sin(x) x dx 无 限 区 间 的 高 斯 型 积 分 公 式 : 高 斯 拉 盖 尔 积 分 公 式 : 0 e x f(x)dx = n w i f(x i ) (2.15) i=1

2.1 数 值 积 分 方 法 19 高 斯 厄 米 积 分 公 式 : e x2 f(x)dx = n w i f(x i ) (2.16) i=1 这 里 x i 分 别 为 n 阶 拉 盖 尔 多 项 式 和 厄 米 多 项 式 的 根. 其 数 值 和 对 应 的 权 重 的 数 值 在 附 录 的 表 中 给 出. 练 习 : 试 编 写 上 述 二 种 高 斯 积 分 公 式 的 程 序, 并 利 用 所 编 程 序 计 算 积 分 0 e x 1 + x 4 dx, 0 e x2 os(x)dx 各 阶 高 斯 积 分 公 式 的 取 点 并 不 重 和, 当 增 加 分 点 时, 所 有 的 函 数 值 都 要 重 算 1965 年 Kronrod 提 出 了 结 合 高 斯 方 法 和 加 速 收 敛 技 巧 的 算 法 这 一 算 法 也 可 以 处 理 有 一 定 奇 异 性 的 被 积 函 数, 但 需 要 把 奇 点 放 在 积 分 限 上 在 Matlab 中 已 实 现 这 一 算 法, 其 调 用 命 令 是 quadgk. 例 如, 计 算 积 分 定 义 函 数 funsin, 并 存 入 funsin.m 0 sin x x e x dx funtion y=funsin(x) end y=sin(x).*exp(-x)./x; 在 Matlab 命 令 窗 口 键 入 如 下 命 令 format long quadgk(@funsin,0,inf) ans = 0.785398163397455 这 一 积 分 的 精 确 结 果 是 π 4 = 0.785398163397448

20 第 二 章 物 理 学 中 的 常 用 数 值 方 法 2.1.4 哥 西 主 值 的 计 算 在 物 理 计 算 中, 哥 西 主 值 的 计 算 占 了 很 重 要 的 位 置, 如 物 质 的 介 电 函 数 的 实 部 和 虚 部 之 间 满 足 著 名 的 Kramers-Kronig 关 系 : ϵ(ω) = ϵ 1 (ω) + iϵ 2 (ω) ϵ 1 (ω) = 1 + 2 π P sϵ 2 (s) s 2 ω ds 2 ϵ 2 (ω) = 2ω π P 0 0 ϵ 1 (s) ds (2.17) s 2 ω2 这 个 关 系 是 一 个 严 格 的 关 系, 它 来 源 于 因 果 关 系, 是 非 常 普 遍 的. 在 物 理 上, ϵ 2 (ω) 的 测 量 和 计 算 都 比 较 容 易, 因 此, 一 般 不 直 接 测 量 或 计 算 ϵ 1 (ω), 而 是 利 用 上 述 关 系 由 ϵ 2 (ω) 的 数 据 进 行 计 算. 上 述 Kramers-Kronig 关 系 中 出 现 的 P 就 代 表 哥 西 主 值 积 分, 设 a < < b, f(x) 在 x = 的 邻 域 内 无 界, 则 哥 西 主 值 积 分 定 义 为 : P b a f(x)dx = lim r 0 [ r a f(x)dx + a 不 失 一 般 性, 我 们 取 = 0 并 将 积 分 取 成 f(x)dx 的 形 式. 令 a b +r ] f(x)dx (2.18) g(x) = 1 [f(x) f( x)], 2 h(x) = 1 [f(x) + f( x)] 2 则 f(x) = g(x) + h(x) 由 于 g(x) 是 奇 函 数 而 h(x) 是 偶 函 数, 于 是 : 因 此, P b a = = 2 r a r a a r f(x)dx = 2 lim f(x)dx + g(x)dx + h(x)dx r 0 a r a r a r f(x)dx g(x)dx + h(x)dx = 2 a 0 r a h(x)dx = h(x)dx + a 0 a r h(x)dx [f(x) + f( x)]dx 这 样 我 们 就 把 哥 西 主 值 积 分 的 计 算 变 为 在 x = 0 处 有 奇 点 的 一 个 反 常 积 分 的 计 算. h(x) 也 可 能 在 x = 0 处 无 奇 点, 此 时 只 要 计 算 一 个 常 规 积 分. 例 如 前 面 的 Kramers-Kronig 关 系 中 积 分 的 计 算 可 化 为 : 2 π P 0 sϵ 2 (s) s 2 ω ds = 1 2 2π P ϵ 2 (s) s ω ds = 1 π 0 ϵ 2 (ω + x) ϵ 2 (ω x) dx x

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

22 第 二 章 物 理 学 中 的 常 用 数 值 方 法 η 积 分 值 0.01 0.596196166026032 + 1.884174006156568i 0.001 0.613114998941270 + 1.903350575685698i 0.0001 0.614827622300549 + 1.905260181139026i 0.00001 0.614999091791223 + 1.905451057233541i 2.1.5 急 速 振 荡 函 数 的 积 分 换 在 各 种 变 换 特 别 是 付 里 叶 变 换 中, 经 常 会 遇 到 急 速 振 荡 函 数 的 积 分, 例 如, 付 里 叶 变 b 又 如 付 里 叶 贝 塞 尔 变 换, a f(x) os(nx)dx, 1 这 里 0 < γ 1 < γ 2 < 是 贝 塞 尔 函 数 的 根. 一 般 我 们 可 把 这 类 函 数 写 为 : I(t) = b a 0 b a f(x)xj n (γ m x)dx f(x) sin(nx)dx f(x)k(x, t)dx, < a < b < (2.22) 式 中 K(x, t) 为 一 振 荡 核 而 f(x) 为 非 振 荡 部 分. 由 于 振 荡 型 函 数 在 积 分 区 间 内 多 次 取 几 乎 相 等 的 正 值 和 负 值, 如 果 用 通 常 的 方 法 计 算, 则 由 于 正 负 数 相 加, 造 成 有 效 数 字 的 损 失, 几 乎 得 不 到 正 确 的 结 果. 下 面 先 给 出 最 常 遇 到 的 振 荡 函 数 的 积 分 付 里 叶 系 数 的 一 种 计 算 方 法, 然 后 再 讨 论 较 为 一 般 的 方 法. 付 里 叶 系 数 的 计 算 : 付 里 叶 系 数 的 复 数 形 式 为 : b 对 上 式 分 部 积 分, 假 定 f(x) 存 在 直 到 n 阶 的 导 数, 可 以 得 到 : b a n 1 f(x)e isx dx = e isa k=0 a f(x)e isx dx (2.23) ( ) k+1 i n 1 f (k) (a) e isb s 当 n 趋 于 无 穷 大 时, 余 项 R n 趋 于 0, 而 求 和 项 就 给 出 了 积 分 的 近 似 值. 练 习 : 写 出 (2.24) 的 实 部 和 虚 部, 即 k=0 ( ) k+1 i f (k) (b) + R n (2.24) s b a f(x) os sxdx 和 b a f(x) sin sxdx 的 积 分 公 式.

2.1 数 值 积 分 方 法 23 上 述 公 式 需 要 计 算 f(x) 的 直 到 n 1 阶 的 导 数, 这 对 于 较 为 复 杂 的 函 数 来 说, 曾 是 一 个 很 重 的 负 担. 不 过 目 前 计 算 机 代 数 的 发 展, 计 算 导 数 已 经 没 有 多 少 困 难. 为 了 看 出 这 一 方 法 的 效 力, 我 们 来 看 一 个 例 子. 计 算 I = 2π 0 x os(x) sin 30xdx 用 公 式 (2.24) 的 虚 部, 取 n = 3, 由 f(x) = x os(x), a = 0, b = 2π, f(0) = 0, f(2π) = 2π, f (0) = 0, f (2π) = 2π, 代 入 得 到 : 2π I = x os(x) sin 30xdx = 1 ( ) 3 1 30 2π + ( 2π) = 0.20967222 30 0 而 该 积 分 的 精 确 值 为 0.20967247. 基 于 Euler-Malaurin 公 式 的 梯 形 法 : 可 以 证 明 ( 见 王 竹 溪, 郭 敦 仁 特 殊 函 数 概 论 第 一 章 ), 函 数 f(x) 的 积 分 可 以 写 成 a+mh [ ] f(a) f(a + mh) f(x)dx = h + f(a + h) + + f(a + (m 1)h) + a 2 2 n ( 1) k B k h 2k [ + f 2k 1 (a + mh) F 2k 1 (a) ] + R n (2k)! 其 中 k=1 R n = θ ( 1)n+1 B n+1 h [ (2n+1) f (2n+1) (a + mh) f (2n+1) (a) ] (2n + 2)! (2.25) 这 里 0 θ 1,B n 是 Bernoulli 数 如 果 f(x) 是 周 期 函 数,f(a + mh) = f(a), 且 f(x) 存 在 直 到 (2n + 1) 阶 的 导 数, 则 显 然 有 f (2k+1) (a) = f (2k+1) (a + mh), 于 是, 利 用 梯 形 法 可 以 得 到 精 确 的 结 果 作 为 例 子, 取 不 同 的 点, 试 用 梯 形 法 计 算 如 下 积 分, 并 于 精 确 结 果 比 较 1 Bessel 函 数 的 积 分 表 示 之 一 是 0 1 1 + 1 2 sin 10πxdx, 精 确 值 = 1.15470054 J n (t) = 1 π π 计 算 n = 2, t = 1, 3, 7, 10 时 Bessel 函 数 的 值 0 os (t sin x nx)dx 由 Euler-Malaurin 公 式 可 知, 如 果 一 个 函 数 满 足 f (a) = f (a + mh), f (3) (a) = f (3) (a+mh),, 则 梯 形 公 式 也 可 以 得 到 很 精 确 的 结 果 例 如 函 数 f(x) = exp(x 2 (1 x 2 )), 在 0 和 1 满 足 上 述 关 系, 从 而, 可 以 用 梯 形 公 式 计 算 1 0 exp(x 2 (1 x 2 ))dx 这 个 积 分 的 精 确 结 果 是 1.03414105, 请 分 别 用 梯 形 法 和 辛 普 森 法 计 算 并 比 较 其 结 果

24 第 二 章 物 理 学 中 的 常 用 数 值 方 法 子 区 间 法, 在 一 般 情 况 下, 如 果 K(x, t) = 0 的 根 容 易 求 出, 且 记 为 a x 1 < x 2 < < x n b, 则 可 用 普 通 的 积 分 规 则 计 算 每 一 个 子 积 分 xi+1 x i, 并 把 结 果 加 起 来. 用 这 种 方 法 得 到 的 一 般 是 一 个 交 错 级 数, 且 相 邻 项 的 绝 对 值 比 较 接 近, 其 求 和 可 用 Euler 变 换 计 算. (Euler 变 换 将 在 后 面 讨 论 ) 在 讲 了 样 条 插 值 之 后, 我 们 再 介 绍 一 种 基 于 样 条 插 值 的 计 算 积 分 (2.23) 的 方 法. 更 多 的 方 法 可 以 参 看 (P. J. Davis and P. Rabinowitz,Methods of Numerial Integration, Aademi Press, 1984) 2.1.6 高 维 积 分 的 计 算 高 维 积 分 的 计 算 原 则 上 可 以 通 过 化 成 累 次 积 分 来 进 行, 例 如, 对 于 积 分 I = g(x, y, z) dxdydz (2.26) 其 中 Ω 为 积 分 区 域, 总 可 以 化 为 I = b dx Ω d(x) dy f(x,y) a (x) e(x,y) g(x, y, z) (2.27) 这 相 当 于 依 次 计 算 三 个 一 重 积 分, 每 一 积 分 以 其 前 一 积 分 作 为 被 积 函 数. 稍 加 考 虑 就 可 发 现, 这 一 过 程 的 计 算 量 随 积 分 重 数 的 增 加 而 变 得 无 法 完 成. 如 每 重 积 分 取 10 个 点 ( 这 是 一 个 非 常 小 的 数 字 ), 则 n 重 积 分 就 要 计 算 10 n 个 函 数 值. 因 此, 这 一 方 法 只 能 用 于 积 分 重 数 较 低 的 积 分. 高 维 积 分 还 可 以 用 高 斯 方 法 计 算, 这 只 要 把 每 重 积 分 用 高 斯 积 分 公 式 代 入 即 可. 对 于 n > 3 的 高 维 积 分, 计 算 量 总 是 相 当 大 的, 如 果 你 对 计 算 精 度 要 求 不 高, 那 么 可 以 利 用 后 面 将 要 讲 到 的 Monte Carlo 方 法. 另 外, 如 果 可 能, 请 尽 可 能 地 解 析 积 出 内 层 积 分, 以 减 小 计 算 量. 2.1.7 固 体 热 容 量 的 计 算 在 德 拜 近 似 下, 固 体 的 热 容 量 由 下 式 给 出 ( ) 3 T Θd /T C V = 9Nk Θ D 0 ξ 4 e ξ (e ξ 1) 2 dξ 式 中 Θ D 为 德 拜 温 度. 为 了 求 得 热 容 量, 需 要 计 算 上 式 中 的 积 分, 这 一 积 分 只 有 当 Θ d /T 1 或 Θ d /T 1 时 才 能 解 析 积 出, 一 般 情 况 下 只 能 进 行 数 字 计 算. 考 虑 f(x) = x 0 ξ 4 e ξ (e ξ 1) 2 dξ

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

26 第 二 章 物 理 学 中 的 常 用 数 值 方 法 2.2.2 牛 顿 法 和 割 线 法 在 所 有 一 元 函 数 求 根 方 法 中, 牛 顿 法 也 许 是 最 值 得 推 荐 的 一 个. 其 想 法 是 这 样 的, 假 定 f(x) = 0 的 根 为 x, 而 我 们 有 一 个 对 x 的 猜 测 值 x n, f(x) 可 在 x n 附 近 展 开 为 泰 勒 级 数 : f(x) = f(x n ) + f (x n )(x x n ) + 略 去 高 次 项, 令 上 式 为 0, 可 解 出 x, 将 此 值 作 为 x 的 一 个 更 好 的 近 似, 记 为 x n+1, 则 x n+1 = x n f(x n) f (x n ). (2.28) 这 样, 给 定 一 个 x 的 初 始 值 x 0, 利 用 下 面 的 迭 代 公 式 反 复 迭 代 : x n+1 = x n f(x n) f (x n ) (2.29) 如 果 迭 代 收 敛, 则 必 收 敛 于 f(x) = 0 的 一 个 根 x. 这 一 迭 代 公 式 (2.29) 就 是 著 名 的 牛 顿 迭 代 公 式. 可 以 证 明, 如 果 迭 代 到 第 n 步 时 的 误 差 为 ϵ n = x n x, 则 第 n + 1 步 的 误 差 为 f (x ) ϵ n+1 = ϵ 2 n, 也 就 是 说, 牛 顿 法 是 平 方 收 敛 的. 证 明 如 下 : 由 迭 代 公 式 (2.29), 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 ϵ n + 1 f (x ) 2 f (x ) ϵ2 n f ) (x ) f (x ) ϵ2 n = x + 1 f (x ) 2 f (x ) ϵ2 n 二 分 法 是 线 性 收 敛 的, 请 证 明 之. 但 是, 另 一 方 面, 对 于 任 一 初 始 值, 牛 顿 法 并 不 总 是 收 敛, 事 实 上, 牛 顿 法 的 收 敛 范 围 是 较 小 的, 因 此, 在 实 际 使 用 时, 必 须 找 一 个 较 好 的 x 0, 使 x 0 尽 可 能 地 靠 近 x, 这 可 以 把 二 分 法 与 牛 顿 法 结 合 起 来 以 完 成. 试 计 算 10. 令 f(x) = x 2 10, f(x) = 0 的 根 就 是 所 求 结 果 取 初 始 值 x 0 = 3, f (x) = 2x. 迭 代 结 果 如 下 : 3. 3.166666666666667 3.162280701754386 3.162277660169842 3.162277660168380 3.162277660168379

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

28 第 二 章 物 理 学 中 的 常 用 数 值 方 法 的 最 优 化 问 题 等 价. 这 里 w i, i = 1, 2,, 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 + 1.618(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) (2.32) 这 里 w 是 一 个 比 例 因 子, 我 们 现 在 来 确 定 它. 由 于 要 求 每 一 步 都 满 足 上 述 关 系, 显 然 有 x b = w( b) = w(( a) (b a)) = w( a) w 2 ( a) = (w w 2 )( a) (2.33) 另 一 方 面, 由 (2.32) 可 得 (b a) + ( x) = ( a) (x b) = 2w( a) (2.34) 即 由 (2.33) 和 (2.35), 我 们 得 到 : x b = (1 2w)( a) (2.35) w 2 3w + 1 = 0 (2.36) 由 此 解 得, w = 3 5 0.38197, 1 w 0.61803 2 1 w 通 常 被 称 为 黄 金 分 割 数 或 黄 金 数. 一 般 在 开 始 计 算 时, a, b, 三 者 并 不 满 足 黄 金 分 割 关 系, 但 几 次 迭 代 后, 这 一 关 系 便 可 满 足, 并 在 其 后 一 直 保 持.

2.2 方 程 的 求 根 和 最 优 化 方 法 29 2.2.4 高 维 最 优 化, 共 扼 梯 度 方 法 这 一 节 给 出 多 维 问 题 的 最 优 化 方 法. 对 于 一 个 含 有 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 ) = min t F (x i + tp i ) (2.37) 则 x i+1 = x i + t i p i 便 为 所 求. 因 此, 高 维 最 优 化 的 的 不 同 方 法 实 际 上 是 选 择 不 同 的 一 维 寻 找 方 向. 一 个 直 观 的 推 断 就 是 在 每 一 点 选 择 该 点 目 标 函 数 下 降 最 快 的 方 向 p i = F (x i ) (2.38) 作 为 寻 找 方 向, 这 一 方 法 称 为 最 陡 下 降 法. 直 观 的 想 象 可 能 会 认 为 最 陡 下 降 方 向 是 一 种 理 想 的 寻 找 方 向, 但 事 实 并 非 如 此, 最 陡 下 降 方 向 只 表 示 在 x i 附 近 的 下 降 性 质, 而 对 整 个 最 优 化 过 程 来 说, 我 们 寻 找 的 是 全 局 的 最 优 方 向. 作 为 一 个 例 子, 考 虑 如 下 函 数 的 最 优 化 : F = x 2 + 10y 2 (2.39) 显 然, 这 一 函 数 的 极 小 点 是 (0, 0). 函 数 (2.39) 的 负 梯 度 为 g = (2x, 20y) T (2.40) 如 果 从 点 (2, 1) 出 发, 利 用 最 陡 下 降 法 计 算 所 得 逐 x 如 下 : (2, 1), (1.792829, 0.035857), (0.461013, 0.230507), (0.413259, 0.008265), (0.106267, 0.053133), (0.095259, 0.001905), (0.024495, 0.012248), (0.021958, 0.000439), (0.005646, 0.002823), (0.005061, 0.000101), (0.001302, 0.000651), (0.001167, 0.000023), (0.000300, 0.000150), (0.000269, 0.000005), (0.000069, 0.000035), (0.000062 0.000001). 显 然, 收 敛 是 很 慢 的. 最 陡 下 降 方 向 并 不 是 理 想 的 下 降 方 向, 为 了 改 进 最 优 化 的 收 敛 速 度, 需 要 寻 求 其 他 方 向. 为 了 判 断 一 个 寻 找 方 向 的 好 坏, 应 该 有 一 个 标 准, 注 意 到 一 般 函 数 在 最 小 点 附 近 近 似 于 二 次 函 数, 因 此 如 果 一 个 方 法 对 二 次 函 数 比 较 有 效, 则 可 望 其 对 一 般 函 数 也 比 较 有 效, 如 果 一 个 方 法 对 二 次 函 数 的 效 果 都 不 好 的 话, 很 难 期 望 它 对 一 般 函 数 会 有 好 的 效 果.

30 第 二 章 物 理 学 中 的 常 用 数 值 方 法 计 算 实 践 也 证 实 了 这 一 论 断. 下 面, 我 们 就 从 二 次 函 数 寻 求 一 种 比 较 有 效 的 下 降 方 向. 为 此, 我 们 先 给 出 几 个 数 学 定 理 ; 定 理 一 : 若 N 维 欧 几 里 德 空 间 中 的 向 量 q 与 N 个 线 性 独 立 向 量 p 1, p 2,, p N 都 正 交, 则 向 量 q 必 为 0. 这 一 定 理 的 结 论 是 显 然 的. 定 义 : A 共 轭 向 量. 设 A 是 一 个 N N 的 对 称 正 定 矩 阵, p, q 是 两 个 N 维 向 量, 若, p T Aq = 0 (2.41) 则 称 向 量 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 = 0 对 i = 1, 2,, N, 用 p T i A 左 乘 上 式 得 α i p T i Ap i = 0 因 为 p i 0, A 正 定, 所 以 p T i Ap i > 0 从 而 必 有 α i = 0, i = 1, 2,, N 因 此, 向 量 组 p 1, p 2,, p N 线 性 独 立. 有 了 以 上 数 学 准 备, 我 们 来 考 虑 二 次 函 数 的 最 优 化 问 题, 我 们 的 目 的 是 找 一 种 下 降 方 向, 能 够 在 有 限 步 达 到 二 次 函 数 的 极 小 点. 考 虑 二 次 N 元 函 数 F (x) = a + b T x + 1 2 xt Ax (2.42) 其 中 A 为 一 正 定 矩 阵. 函 数 (2.44) 的 梯 度 为 G(x) = b + Ax (2.43) 记 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. 现 在 有 了 两 个

2.2 方 程 的 求 根 和 最 优 化 方 法 31 方 向,g 1 和 g 0, 最 陡 下 降 法 选 取 新 的 方 向 为 g 1. 我 们 保 留 原 来 方 向 的 一 点 信 息, 选 新 的 寻 找 方 向 为 这 两 个 方 向 的 组 合 p 1 = g 1 α 0 g 0 (2.44) 并 要 求 p 1 与 p 0 为 A 共 轭, 则 应 有 ( g 1 α 0 g 0 ) T Ap 0 = 0 (2.45) 由 于 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 (2.46) 因 此, 方 程 (2.45) 成 为 ( g 1 α 0 g 0 ) T (g 1 g 0 ) = 0 (2.47) 由 此 解 得 这 里 假 定 g 0 0, 否 则, 若 g 0 = 0, 则 x 0 就 是 所 求 的 极 值 点. α 0 = gt 1 g 1 g T 0 g 0 (2.48) 把 α 0 代 入 (2.44), 便 得 到 p 1. 沿 p 1 求 出 极 小 点 x 2 并 算 出 g 2. 由 于 p 0 和 p 1 为 A 共 轭, 因 此 由 (2.46) g T 0 (g 2 g 1 ) = 0 考 虑 到 g 0 与 g 1 正 交, g T 0 g 1 = 0, 由 上 式 可 知 g 0 与 g 2 也 正 交. 因 x 2 是 沿 p 1 方 向 的 极 小 点, 故 g 2 p 1, 也 就 是 g T 2 ( g 1 α 0 g 0 ) = 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 (2.49) 要 求 p 2 分 别 与 p 0 和 p 1 为 A 共 轭, 也 就 是 要 求 如 下 方 程 成 立 ( g 2 α 1 g 1 β 0 g 0 ) T (g 1 g 0 ) = 0 (2.50) ( g 2 α 1 g 1 β 0 g 0 ) T (g 2 g 1 ) = 0 (2.51) (2.52)

32 第 二 章 物 理 学 中 的 常 用 数 值 方 法 由 此 解 出 于 是 α 1 = gt 2 g 2 g T 1 g 1, β 0 = α 0 α 1. (2.53) p 2 = g 2 + α 1 p 1 (2.54) 现 假 定 已 求 得 p 0, p 1,, p j, 其 中 p 0 = g 0 p i = g i + α i 1 p i 1 α i 1 = g T i g i /g T i 1g i 1 i = 1, 2,, j 这 里 假 定 g 0, g 1,, g j 都 不 为 零 且 构 成 一 个 正 交 系. 沿 p j 求 出 x j+1 并 算 出 g j+1 = G(x j+1 ), 设 g j+1 0, 则 仿 前 可 证, 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 (2.55) 并 要 求 p j+1 与 p 0, p 1,, p j 为 A 共 轭, 即 ( g j+1 β 0 g 0 β 1 g 1 β j 1 g j 1 α j g j ) T (g j+1 g j ) = 0 (2.56) 得 到 α j = g T j+1g j+1 /g T j 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

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

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

2.2 方 程 的 求 根 和 最 优 化 方 法 35 ϕ i (x) } 的 一 个 线 性 组 合, y = i a i ϕ i (x) (2.64) 然 后 用 曲 线 拟 合 方 法 求 出 系 数 a i. 函 数 族 { ϕ i (x) } 通 常 从 物 理 考 虑 来 选 取, 一 种 常 用 的 方 法 是 选 取 { 1, x, x 2,, x m 1 } 用 于 拟 合 m 个 参 数 ; 也 可 选 取 三 角 函 数 { 1, sin x, os x, sin 2x, os 2x,, sin 2mx, os 2mx} 用 于 拟 合 2m + 1 个 参 数. 这 种 方 法 称 为 线 性 拟 合 方 法, 即 y 线 性 地 依 赖 于 参 数 a i, 对 于 有 些 物 理 问 题, y 和 x 的 函 数 关 系 中 非 线 性 地 依 赖 于 参 数 a i, 例 如 y = a 1 e a2x + a 3 e a4x sin(a 5 x + a 6 ) 就 是 一 个 具 有 6 个 参 数 的 非 线 性 拟 合 问 题. 曲 线 拟 合 通 常 采 用 最 小 二 乘 法, 如 果 x 的 测 量 误 差 可 以 忽 略, 观 测 点 i 的 残 差 δ i 定 义 为 在 该 点 y 的 观 测 值 与 理 论 值 之 差, 即 δ i = y i f(x i ; a 1, a 2,, a m ) (2.65) 若 观 测 值 y i 的 方 差 为 σi 2, 则 定 义 其 权 重 因 子 为 记 χ 2 为 观 测 点 残 差 的 加 权 平 方 和, χ 2 = N w i δi 2 = i=1 w i = 1 σ 2 i N w i (y i f(x i ; a 1, a 2,, a m )) 2 (2.66) i=1 曲 线 拟 合 的 最 小 二 乘 法 就 是 以 χ 2 为 目 标 函 数 的 优 化 问 题, 即 寻 找 一 组 参 数 a i, (i = 1, 2,, m) 使 χ 2 为 最 小. 在 观 测 点 的 方 差 未 知 的 情 况 下, 可 以 取 w i = 1, 也 可 以 在 每 一 观 测 点 多 测 量 几 次, 以 其 平 均 值 y i 作 为 该 点 的 观 测 值, 而 以 y 2 i y i 2 作 为 其 方 差. 在 曲 线 拟 合 问 题 中, 最 困 难 的 是 选 取 合 适 的 经 验 公 式, 这 需 要 结 合 物 理 考 虑 和 对 测 量 数 据 的 观 察, 也 依 赖 于 对 各 种 函 数 及 其 图 像 的 明 确 概 念. 例 题 在 某 一 能 谱 测 量 中, 得 到 下 表 所 示 数 据, 试 给 出 经 验 公 式 并 作 出 曲 线 拟 合. x 0.0 0.1 0.2 0.3 0.4 0.5 0.6 y 2.117 2.634 3.459 4.671 6.600 9.978 15.893 x 0.7 0.8 0.9 1.0 1.1 1.2 1.3 y 25.108 30.979 25.768 17.363 12.031 10.454 10.280 x 1.4 1.5 1.6 1.7 1.8 1.9 y 11.012 11.138 9.745 7.576 5.598 4.165 把 表 中 数 据 点 到 坐 标 纸 上 ( 见 图 2.1 中 圈 点 ), 可 以 看 出 数 据 有 二 个 峰, 为 此, 我 们 试 用 Lorentz 形 曲 线 去 拟 合, 即 取 经 验 公 式 为 y = a 1 (x a 2 ) 2 + a 2 3 a 4 + (x a 5 ) 2 + a 2 6

36 第 二 章 物 理 学 中 的 常 用 数 值 方 法 40 30 20 10 0 0 0.5 1 1.5 2 图 2.1: 用 上 式 和 上 表 中 的 数 据 计 算 χ 2, 并 用 共 轭 梯 度 方 法 优 化 χ 2, 求 得 各 参 数 值 为 a 1 = 1.210 a 2 = 0.800 a 3 = 0.202 a 4 = 0.776 a 5 = 1.502 a 6 = 0.295 在 图 2.1 中, 我 们 用 实 线 画 出 了 拟 合 的 结 果, 虚 线 给 出 了 二 个 峰 的 图 像. 2.3 函 数 的 求 值 在 数 值 计 算 中, 经 常 要 求 一 些 函 数 的 值, 一 些 常 用 的 初 等 函 数 如 三 角 函 数, e 指 数 函 数, 双 曲 函 数, 对 数 函 数 等, 均 已 写 入 计 算 机 语 言 的 定 义 做 为 内 部 函 数, 在 使 用 时 只 需 调 用 即 可, 而 其 它 的 函 数 和 表 达 式 则 要 自 己 计 算, 本 章 将 给 出 一 些 在 物 理 计 算 中 经 常 遇 到 的 函 数 及 表 达 式 的 计 算 方 法 和 部 分 程 序. 一 般 说 来, 计 算 机 的 内 部 函 数 都 是 由 专 家 设 计 的, 其 效 率 和 精 度 都 很 好, 因 此, 如 果 计 算 机 语 言 已 经 提 供 了 所 需 的 函 数, 则 应 尽 量 使 用 这 些 函 数, 下 面 提 供 的 方 法 仅 仅 作 为 计 算 机 系 统 不 提 供 这 些 函 数 时 的 补 充. 2.3.1 多 项 式 的 求 值 和 级 数 求 和 所 谓 多 项 式 是 指 由 下 式 给 出 的 表 达 式 : p(x) = a 0 + a 1 x + a 2 x 2 + a 3 x 3 + + a n x n 为 了 计 算 p(x) 在 给 定 x 时 的 数 值, 我 们 用 (n), n = 0, 1, 2,, n 来 代 表 多 项 式 的 系 数, 则 可 以 写 出 计 算 语 句. 以 n = 4 为 例, 直 接 用 FORTRAN 语 言 写 出 的 算 式 为 :

2.3 函 数 的 求 值 37 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 我 们 当 然 不 可 能 计 算 到 无 穷 多 项, 所 以 上 述 级 数 必 需 在 某 处 截 断, 如 果 级 数 收 敛 很 慢, 就 要 计 算 很 多 项, 一 方 面 占 用 大 量 的 CPU 时 间, 更 重 要 的 是 当 计 算 的 项 数 非 常 多 时, 舍 入 误 差 的 积 累 有 可 能 完 全 破 坏 整 个 计 算. 为 此 必 须 使 用 加 速 收 敛 的 技 巧. Aitken 算 法 是 常 用 的 加 速 技 巧 之 一. 它 是 对 相 邻 的 三 个 部 分 和 进 行 处 理, 得 到 一 个 更 好 的 近 似. 记 部 分 和 为 : S n = n u i (2.67) i=0 则 由 S n 1, S n, S n+1 可 计 算 出 S n = S n+1 (S n+1 S n ) 2 S n+1 2S n + S n 1 (2.68) 对 n = 1, 2, 3,, 先 求 出 S n, 同 时 算 出 S n, 一 般 S n 更 快 地 收 敛 于 S. 在 实 际 使 用 时, (2.68) 最 好 用 写 出 的 式 子 计 算, 一 些 代 数 上 等 价 的 式 子 可 能 导 致 较 大 的 舍 入 误 差. Euler 变 换 对 于 由 下 式 定 义 的 交 错 级 数 ( 1) s u s (2.69) 这 里 u s > 0. Euler 变 换 是 一 个 强 有 力 的 算 法, 它 由 下 式 给 出 ( 取 n 为 偶 数 ) 其 中 : s=0 ( 1) s u s = u 0 u 1 + u 2 u n 1 + s=0 s=0 ( 1) s 2 s+1 ( s u n ) (2.70) u n u n+1 u n 2 u n ( u n ) = u n+2 2u n+1 + u n (2.71)

38 第 二 章 物 理 学 中 的 常 用 数 值 方 法 为 向 前 差 分. 式 (2.70) 可 推 导 如 下 : 定 义 位 移 算 符 E 为 : u n+1 Eu n 由 于 (1 + )u n = u n + u n+1 u n = u n+1 = Eu n, 因 此 有 算 符 等 式 E = 1 + 而 u 0 u 1 + u 2 u 3 + = (1 E + E 2 E 3 + )u 0 = 1 1 + E u 0 = 1 1 u 2 1 + 0 = ( 1) s 2 s+1 s u 0 2 s=0 2.3.2 连 分 式 的 求 值 所 谓 连 分 式, 是 指 由 下 式 给 出 的 表 达 式, f = b 0 + b 1 + b 2 + a 1 b 3 + a 2 a 3 a 4 b 4 + a 5 b 5 + 数 学 家 们 对 它 做 过 很 多 研 究, 这 是 一 个 非 常 有 趣 的 函 数 表 示 方 法. 除 了 上 面 的 表 达 式 外, 人 们 还 经 常 使 用 另 一 种 记 法. 这 里, 系 数 a n 和 b n 也 可 以 是 x 的 函 数. a 2 a 3 a 4 a 5 f = b 0 + a 1 b 1 + b 2 + b 3 + b 4 + b 5 + 由 于 连 分 式 只 能 从 右 向 左 求 值, 因 此 在 具 体 计 算 之 前, 很 难 判 断 一 个 无 限 的 连 分 式 应 该 在 何 处 截 断, 一 种 直 接 的 方 法 就 是 猜 一 个 截 断 点, 求 出 一 个 值, 再 猜 一 个 更 大 一 点 的 截 断 点 计 算, 通 过 比 较 不 同 截 断 处 的 结 果 来 判 断 是 否 收 敛. 这 显 然 不 是 一 个 好 的 方 法. 另 一 种 在 实 践 中 常 用 的 方 法 是 用 一 个 分 式 有 理 函 数 来 近 似 连 分 式, 用 f n 表 示 计 算 到 a n 和 b n 的 连 分 式 的 值, 则 : 这 里 A n, B n 由 下 面 的 递 推 公 式 给 出 : A 1 1 B 1 0 f n = A n B n (2.72) A 0 b 0 B 0 1 (2.73) 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 = 1, 2,, n

2.4 常 微 分 方 程 的 数 值 解 法 39 练 习 试 证 明 递 推 关 系 (2.73) 并 写 出 计 算 机 程 序 2.4 常 微 分 方 程 的 数 值 解 法 常 微 分 方 程 在 物 理 学 中 具 有 重 要 意 义, 牛 顿 第 二 定 律 就 表 示 为 一 组 二 阶 常 微 分 方 程, 第 一 章 中 提 到 的 著 名 的 三 体 问 题 就 是 由 九 个 二 阶 常 微 分 方 程 组 成 的 方 程 组 ; 在 量 子 力 学 中 和 电 动 力 学 中, 由 薛 定 谔 方 程 及 Maxwell 方 程 分 离 变 量 得 到 的 所 谓 Sturm Liouville 本 征 值 问 题 也 是 二 阶 常 微 分 方 程 的 求 解 问 题. 在 这 一 章 中, 我 们 将 给 出 几 种 数 值 求 解 常 微 分 方 程 的 常 用 方 法. 2.4.1 Euler 方 法 对 于 常 微 分 方 程 的 初 值 问 题 dy dx = f(x, y) y(0) = η (2.74) 数 值 求 解 的 第 一 步 是 把 方 程 (2.74) 变 为 一 个 差 分 方 程, 为 此, 让 x 取 一 组 分 立 值 x n, n = 0, 1, 2,, x 0 = 0, x n+1 = x n + h, 这 里 h 称 为 积 分 的 步 长, 把 y(x) 在 x n 处 展 开, 有 y(x n+1 ) = y(x n ) + h dy dx + O(h 2 ) (2.75) x=xn 利 用 方 程 (2.74), 便 可 得 到 y(x n+1 ) = y(x n ) + hf(x n, y n ) + O(h 2 ) (2.76) 这 就 是 Euler 差 分 格 式, 从 x = 0 出 发, 利 用 这 一 方 程 一 步 一 步 地 向 前 计 算, 就 可 以 得 到 任 一 x n 处 的 函 数 值 y n y(x n ). Euler 方 法 是 一 个 显 式 方 法, 即 第 n + 1 步 的 y n+1 只 依 赖 于 x n, y n 的 值, 而 且 每 一 步 只 需 计 算 一 次 f(x n, y n ) 的 值, 这 是 Euler 方 法 的 优 点, 但 另 一 方 面, Euler 方 法 也 存 在 严 重 的 不 足, 因 为 这 一 方 法 是 一 次 的 ( 每 一 步 的 误 差 项 为 h 2 的 数 量 级, y(x N ) Nh 2 h), 所 以 h 必 须 取 得 足 够 小, 以 保 证 有 满 意 的 精 度, 但 步 长 太 小 时, 积 分 到 一 预 定 的 x 则 要 做 很 多 步 的 计 算, 舍 入 误 差 的 积 累 将 会 使 计 算 结 果 严 重 失 真. 由 于 这 一 原 因, 在 实 际 计 算 中 几 乎 不 使 用 Euler 方 法, 而 常 常 利 用 其 推 导 上 的 简 单 和 概 念 上 非 常 清 楚 而 用 于 教 学.