计算物理学习提纲

计算物理

课程内容

• 科学计算基本原则
• 插值与函数拟合
• 数值微积分
• 非线性方程求根与方程组
• 矩阵、特征值与特征向量
• 常微分方程数值解
• 偏微分方程初边值问题
• Monte Carlo 方法

科学计算基本流程

1. 建立描述实际系统的模型，能解析求解时优先分析其解析性质。
2. 设计算法，将连续问题离散化为计算机可处理的形式。
3. 编写程序并运行，得到数值结果。
4. 检验结果，并与实验结果对比以修正模型或算法。

计算机适合处理离散系统，以及四则运算、逻辑判断和迭代计算；使用时必须注意其精度有限。

基本原则：数值模型允许合理近似，但结果必须可用。

数值误差

误差来源

• 模型误差：理论模型与实际系统之间的差异
• 测量误差：输入数据由测量引入的误差
• 截断误差：有限项展开或有限次迭代对无限过程的近似
• 舍入误差：有限精度数值表示与运算产生的误差

绝对误差与相对误差

设真值为 $x^{\ast }$ ，近似值为 $x$ 。

绝对误差：

Copy LaTeX$$|x-x^{\ast }|.$$

相对误差：

Copy LaTeX$$\frac{|x-x^{\ast }|}{|x^{\ast }|}.$$

真值通常不可得，实际估计时可用可计算的近似值代替分母中的真值。

一阶误差传播

若

Copy LaTeX$$y=f(x_1,\dots ,x_n),$$

则在误差较小时，可用一阶 Taylor 展开估计输出误差：

Copy LaTeX$$y-y^{\ast }\approx \sum _{i=1}^n\frac{\partial f(x_1^{\ast },\dots ,x_n^{\ast })}{\partial x_i^{\ast }}(x_i-x_i^{\ast }).$$

判断方法：偏导数的绝对值越大，相应输入误差对结果的放大越明显。

稳定计算方法

选择稳定的递推方向

对

Copy LaTeX$$I_n=e^{-1}{\int }_0^1{x}^n{e}^{x}dx,$$

分部积分得到

Copy LaTeX$$I_n=1-n{I}_{n-1}.$$

正向递推会把上一步误差乘以 $n$ ，误差可能迅速放大，甚至得到与 $I_n>0$ 相矛盾的结果。

改写为

Copy LaTeX$$I_{n-1}=\frac{1}{n}(1-I_n)$$

并反向递推时，误差每步缩小为原来的 $1/n$ ，数值上更稳定。

高阶初值可利用

Copy LaTeX$$\frac{e^{-1}}{n+1}<I_n<\frac{1}{n+1}$$

进行估计，再向低阶回推。

要点：代数等价的公式未必具有相同的数值稳定性，应选择抑制误差传播的计算方向。

秦九韶算法

多项式

Copy LaTeX$$P(x)=a_n{x}^n+\cdots +a_{1}x+a_0$$

可写成嵌套形式：

Copy LaTeX$$P(x)=\{\cdots [(a_{n}x+a_{n-1})x+a_{n-2}]\cdots \}x+a_0.$$

秦九韶算法只需 $n$ 次乘法和 $n$ 次加法，避免分别计算各次幂。

避免相消误差

两个大小相近的数相减时，有效数字容易大量损失。例如

Copy LaTeX$$\sqrt{x+1}-\sqrt{x}$$

在 $x$ 很大时应有理化为

Copy LaTeX$$\frac{1}{\sqrt{x+1}+\sqrt{x}}.$$

判断方法：若表达式包含两个近似相等的大数相减，优先寻找因式分解、有理化或其他等价形式。

避免以过小的数作除数

除数的绝对值过小时，输入或舍入误差会被显著放大。计算前应检查分母是否接近零，并评估问题的条件性。

割圆术与稳定实现

面积近似思路

对半径为 $r$ 的圆，用内接正多边形逼近。设当前边数为 $m$ 、边长为 $l$ ，半边长与中心到边的距离分别为

Copy LaTeX$$x=\frac{l}{2},h=\sqrt{r^2-x^2}.$$

每个等腰三角形的面积为 $xh$ ，故正多边形面积为

Copy LaTeX$$S_m=mxh.$$

不断将边数加倍时，$S_m$ 从下方趋近圆面积 $\pi r^2$ ；当 $r=1$ 时即趋近 $\pi$ 。

手写代码中的问题

• x、S 未声明，会导致编译错误。
• 计算了 S 却没有 printf，程序运行后看不到近似结果。
• sqrt((r-h)*(r-h)+x*x) 中的 r-h 会在 $h\approx r$ 时发生相消，导致有效数字丢失。
• i 表示边数但名称含义不清；边数不断翻倍，还需留意整数溢出。
• 应检查 scanf 是否成功，并限制迭代次数为正数。

稳定化来自

Copy LaTeX$$r-h=\frac{r^2-h^2}{r+h}=\frac{x^2}{r+h}.$$

于是新边长可计算为

Copy LaTeX$$l_{\mathrm{new}}=\sqrt{(r-h{)}^2+x^2}=\sqrt{\frac{2r{x}^2}{r+h}}.$$

当 $r=1$ 时：

Copy LaTeX$$l_{\mathrm{new}}=\sqrt{\frac{2x^2}{1+h}}.$$

一个可运行且避免该相消的版本：

#include <math.h>
#include <stdio.h>

int main(void)
{
    const double radius = 1.0;
    double side = 1.0;
    long long sides = 6;
    int iterations;

    if (scanf("%d", &iterations) != 1 || iterations <= 0)
        return 1;

    for (int k = 0; k < iterations; ++k) {
        const double half_side = 0.5 * side;
        const double height = sqrt(radius * radius
                                   - half_side * half_side);
        const double area = sides * half_side * height;

        printf("%lld %.17g\n", sides, area);

        side = sqrt(2.0 * radius * half_side * half_side
                    / (radius + height));
        sides *= 2;
    }

    return 0;
}

多项式插值

给定 $n+1$ 个互异节点

Copy LaTeX$$a<x_0<\cdots <x_n<b,y_i=f(x_i),$$

存在唯一的次数不超过 $n$ 的多项式 $P_n$ 满足

Copy LaTeX$$P_n(x_i)=y_i(i=0,\dots ,n).$$

唯一性的两种判断方法：

• 系数方程的 Vandermonde 行列式为 $\prod _{0\le j<i\le n}(x_i-x_j)\ne 0$ 。
• 若两个插值多项式都满足条件，它们的差次数不超过 $n$ 却有 $n+1$ 个不同零点，只能恒为零。

单项式组 $\{1,x,\dots ,x^n\}$ 是 ${\mathcal{P}}_n$ 的基；没有另行指定内积时，不能称它为正交基。

Lagrange 插值

基函数

Copy LaTeX$$l_i(x)=\prod _{\begin{array}{c}j=0\\ j\ne i\end{array}}^n\frac{x-x_j}{x_i-x_j}$$

满足 $l_i(x_j)={\delta }_{ij}$ 。插值多项式为

Copy LaTeX$$L_n(x)=\sum _{i=0}^n{y}_i{l}_i(x).$$

一次插值的常用形式：

Copy LaTeX$$L_1(x)=y_0+\frac{y_1-y_0}{x_1-x_0}(x-x_0).$$

优点：公式对称、可以直接写出。缺点：新增节点后通常需要重新计算全部基函数，递推性差。

差商

一阶和高阶差商递推定义为

Copy LaTeX$$f[x_0,x_1]=\frac{f(x_1)-f(x_0)}{x_1-x_0},$$

Copy LaTeX$$f[x_0,\dots ,x_k]=\frac{f[x_1,\dots ,x_k]-f[x_0,\dots ,x_{k-1}]}{x_k-x_0}.$$

常用性质：

• 对节点顺序对称
• 对函数满足线性
• 次数低于 $k$ 的多项式，其 $k$ 阶差商为零
• 次数等于 $k$ 的多项式，其 $k$ 阶差商等于最高次项系数
• 若 $f\in C^k$ ，存在节点区间内的 $\xi$ 使 $f[x_0,\dots ,x_k]=f^{(k)}(\xi )/k!$

Newton 插值

Copy LaTeX$$\begin{array}{rl}{N}_n(x)& =f[x_0]+(x-x_0)f[x_0,x_1]\\ & +(x-x_0)(x-x_1)f[x_0,x_1,x_2]+\cdots \\ & +\prod _{j=0}^{n-1}(x-x_j)f[x_0,\dots ,x_n].\end{array}$$

Newton 形式和 Lagrange 形式表示同一个唯一插值多项式。Newton 形式的优势是递推性：加入新节点 $x_{n+1}$ 时，只需增加一项

Copy LaTeX$$\prod _{j=0}^n(x-x_j)f[x_0,\dots ,x_{n+1}].$$

差商表与原地计算

差商表按阶递推：先由函数值计算全部一阶差商，再依次计算二阶、三阶差商。Newton 多项式使用

Copy LaTeX$$f[x_0],f[x_0,x_1],\dots ,f[x_0,\dots ,x_n]$$

这一组系数。实现时可以使用二维差商表；也可以用一维数组原地更新，覆盖已经不再需要的低阶数据。

Runge 现象

在等距节点上使用高次多项式插值时，增加节点不一定提高整体精度，区间端点附近可能出现越来越强的振荡。例如

Copy LaTeX$$f(x)=\frac{1}{1+25x^2}$$

的高次等距节点插值。

常用处理：

• 使用分段低次插值
• 改用合适的非等距节点
• 使用三次样条等分段方法

Hermite 插值

Hermite 插值不仅匹配节点处的函数值，还匹配指定阶数的导数。若总共有 $m+1$ 个独立插值条件，可唯一确定一个次数不超过 $m$ 的多项式。

重节点差商

差商对节点连续。节点重合时，重节点差商由导数给出：

Copy LaTeX$$f[\underset{k+1\text{ 个}}{\underbrace{x_0,\dots ,x_0}}]=\frac{f^{(k)}(x_0)}{k!}.$$

Taylor 多项式可看作所有节点都重合于 $x_0$ 的 Hermite 插值多项式。

三点三次 Hermite 插值

若给定三个函数值：

Copy LaTeX$$p(x_i)=f(x_i)(i=0,1,2),$$

以及中间节点的一阶导数 $p^{\prime }(x_1)=f^{\prime }(x_1)$ ，可写成

Copy LaTeX$$\begin{array}{rl}p(x)& =f(x_0)+f[x_0,x_1](x-x_0)\\ & +f[x_0,x_1,x_2](x-x_0)(x-x_1)\\ & +A(x-x_0)(x-x_1)(x-x_2),\end{array}$$

再由导数条件确定 $A$ 。余项为

Copy LaTeX$$R(x)=\frac{f^{(4)}({\xi }_x)}{4!}(x-x_0)(x-x_1{)}^2(x-x_2).$$

两点三次 Hermite 插值

在区间 $[x_0,x_1]$ 上给定

Copy LaTeX$$y_i=f(x_i),m_i=f^{\prime }(x_i)(i=0,1),$$

则

Copy LaTeX$$\begin{array}{rl}{H}_3(x)& =y_0(1+2\frac{x-x_0}{x_1-x_0}){\left(\frac{x-x_1}{x_0-x_1}\right)}^2\\ & +y_1(1+2\frac{x-x_1}{x_0-x_1}){\left(\frac{x-x_0}{x_1-x_0}\right)}^2\\ & +m_0(x-x_0){\left(\frac{x-x_1}{x_0-x_1}\right)}^2\\ & +m_1(x-x_1){\left(\frac{x-x_0}{x_1-x_0}\right)}^2.\end{array}$$

其中函数值基函数在对应节点取 $1$ 、在其余节点取 $0$ ，且端点导数均为 $0$ ；导数基函数满足

Copy LaTeX$${\beta }_j(x_i)=0,{\beta }_j^{\prime }(x_i)={\delta }_{ij}.$$

余项为

Copy LaTeX$$R_3(x)=\frac{f^{(4)}({\xi }_x)}{4!}(x-x_0{)}^2(x-x_1{)}^2.$$

分段插值

把 $x_0<x_1<\cdots <x_n$ 划分成相邻小区间，并在每段使用低次多项式，可避免用单个高次多项式覆盖整个区间。

• 只有函数值：可用分段线性插值。
• 同时给出端点函数值和导数：可在每段使用两点三次 Hermite 插值。
• 若缺少导数而又要求更高阶光滑性：可使用三次样条，由整体连续性和边界条件确定导数或系数。