数值计算方法与算法

教学目标 掌握常用的数值计算方法 掌握计算方法的数学原理 学会选择恰当的计算方法 学会使用流行的计算软件

教学计划 时间：7:50-8:35 地点：2106 2-27 绪论 3-06 多项式插值 3-13 3-20 分段和样条插值 3-27 时间：7:50-8:35　　地点：2106 2-27 绪论 3-06 多项式插值 3-13 3-20 分段和样条插值 3-27 数值微分 4-03 数值积分 4-10 4-17 最小二乘法 4-24 方程求根 5-01 五一放假 5-08 5-15 求解线性方程组 5-22 5-29 计算矩阵特征值 6-05 6-12 微分方程数值解 6-19 6-26 复习 7-04 期末考试 7-13 报送成绩

教材及参考 数值计算方法与算法，张韵华、奚梅成、陈效群编，科学出版社，2006。 科学计算导论，Michael T. Heath著，张威、贺华、冷爱萍译，清华大学出版社， 2005。 数值计算方法解题指导，张韵华编，科学出版社，2003。 网络教学资源 http://bb.ustc.edu.cn

联系方式 教师：王新茂 xinmao@ustc.edu.cn 理化大楼16#016，3607565 助教：杨 荣 助教：杨　荣 rongyang@mail.ustc.edu.cn 136-15607693

第0章 绪论

什么是数值计算方法？ 什么是“好的”数值计算方法？ 实际问题 数学问题 近似解 误差小 ─ 误差分析 耗时少 ─ 复杂度分析 抗干扰 ─ 稳定性分析 数学建模 数值计算 实际问题　　　　数学问题　　　　近似解

误差的类型 绝对误差＝真实值－近似值 相对误差＝绝对误差／真实值 误差的来源 原始误差、截断误差、舍入误差 输入 计算 输出 真实值 近似值

一些例子： 计算地球的体积 计算 如何减小计算误差？ 选择好的算法、提高计算精度 范数的定义 满足非负性,齐次性,三角不等式的实函数

常用的向量范数 常用的矩阵范数 矩阵的谱半径 例：计算矩阵 的范数和谱半径。 例：范数在误差估计中的应用

作业一、145页习题6第1,2题. 作业二、利用公式 编写两个计算ex 的C程序(一个用单精度, 另一个用双精度).令x=±1,±5,±10,±15,±20, 比较它们和库函数exp(x)之间的运行时间和计算误差. 思考如何减小误差？

第1章 插值

函数逼近 用未知函数f(x)的值构造近似函数φ(x)。要求误差小、形式简单、容易计算。 常用的函数逼近方法 插值：φ(xi)=yi, i=0,1,…,n. 拟合：||φ(x)-f(x)||尽可能小 通常取 φ(x) = a0φ0(x) + … + anφn(x)，其中 {φi(x)}为一组基函数。

多项式插值 给定平面上n+1个插值点(xi,yi), 构造n次多项式φ(x), 满足φ(xi)=yi, i=0,1,…,n.

单项式插值

Lagrange 插值

Newton 插值

k阶差商 差商表 1 2 … n ...

差商的性质 以x0,…,xn为节点的n次插值多项式φ(x)的首项系数等于f[x0,…,xn]。 证明：分别以x0,…,xn-1和x1,…,xn为节点构造n-1次插值多项式φ1(x)和 φ2(x)，则有 对n用归纳法。 f[x0,…,xn]与x0,…,xn的顺序无关。

误差估计： 证明：设 ，则 　 有n+2个零点。 根据中值定理，存在 于是 。

Hermite插值 给定平面上n+1个插值点(xi,yi,mi), 构造2n+1次多项式φ(x), 满足φ(xi)=yi, φ’(xi)=mi, i=0,1,…,n.

单项式 基函数 Lagrange 基函数

误差估计： 证明：设 ，则 　 有2n+3个零点。根据中值定理，存在 于是 。

Runge现象：并非插值点取得越多越好。 解决办法：分段插值

三次样条插值 给定平面上n+1个插值点(xi,yi), 构造分段三次多项式φ(x), 满足φ(xi)=yi, φ’(x)可微，φ”(x)连续。

作业一、习题1第2,4,6,8,10,12,14,16题。 作业二、在半圆 上随机选取10个点, 构造插值多项式, 画出函数图像, 并比较3种插值方法的计算误差。 作业三、思考3种插值系数之间的关系。 比较3种插值方法的优缺点和应用范围。

第2章 数值微分和数值积分

数值微分 差商法 向前差商 向后差商 中心差商

插值法 在 x 附近取点(xi,f(xi))构造插值多项式φ。 f’(x)≈φ’(x)

例：用向后差商公式计算f’’(0.2), f’’(0.4)。 例：用中心差商公式计算f’(xi)。 例：用向后差商公式计算f’’(0.2), f’’(0.4)。 x 0.0 0.1 0.2 0.3 0.4 f(x) 0.818731 0.904837 1.000000 1.105171 1.221403 f’(x) x 0.0 0.1 0.2 0.3 0.4 f(x) 1.7 1.5 1.6 2.0 1.9 f’(x) f”(x)

例：设xi=x0+i*h, i=1,...,n。计算φ’(xk)。 解：

误差估计 前后差商 中心差商 插值微分

数值积分 插值法

若积分公式对任意m次多项式都取等号，则称积分公式具有至少m阶的代数精度。 插值型积分公式的代数精度≥n。 当积分节点 x0,...,xn 给定时， 代数精度≥n的积分公式唯一。

例：设xi=a+i*h, i=0,...,n, h=(b-a)/n。 计算Newton-Cotes积分 解：

特别，当n=1,2时，积分公式分别称为 梯形公式 Simpson公式 n a1 a2 a3 a4 a5 1 ½ 2 1/6 4/6 3 1/8 3/8 4 7/90 32/90 12/90

误差估计 特别，梯形公式和Simpson公式的误差为 代数精度＝1 代数精度＝3

复化数值积分

梯形公式 Simpson公式

Richardson外推法 我们要计算 假设 则 有比　　 和　　 更高的精度。 误差估计

Romberg积分公式 等分的梯形公式， 瑕积分 重积分 Gauss-Legendre积分

定理：假设 满足 则插值积分公式具有2n+1阶的代数精度。 证明： 课本21页性质1.3：若f(x)为m次多项式， 则f[x0,...,xn,x]为m-n-1次多项式。

求多项式空间在内积 下的标准正交基。 解法1：对任意基作Gram-Schmidt正交化。 解法2：对任意度量方阵作相合对角化。 解法3：求解正交关系的线性方程组。 解法4：Legendre多项式

作业一、习题2第1~11题。 作业二、使用各种方法计算 的值，其中 0<x<1，要求误差<1e-9。

第3章 曲线拟合的最小二乘法

曲线拟合 对区间 I 上的连续函数 f ， 构造特定类型的函数 φ 使 φ≈f 。 对离散数据序列 (xi, yi), i=1,2,…,m， 构造特定类型的函数 φ 使 φ(xi)≈yi 。 最小二乘法 求 φ 使 最小。

多项式拟合 其中 是标准正交基， 。 求 使 最小。

奇异值分解 Moore-Penrose 广义逆 矛盾方程组 　　　　的解

其他类型的离散数据拟合

作业二、对下列数据，用最小二乘法求形如 的经验公式。 作业一、习题3第1~5题。 作业二、对下列数据，用最小二乘法求形如 的经验公式。 x 0.1 0.2 0.3 0.4 0.5 y 0.9059 0.5632 0.3835 0.4244 0.6730 0.6 0.7 0.8 0.9 1.0 1.0490 1.4315 1.6973 1.7608 1.6016

第4章 非线性方程求根

问题 求f(x)=0在区间[a,b]内的实根 求f(x)=0在x0附近的一个实根 求f(x)=0在x0附近的一个复根 求多项式f(x)=0的所有复根 求非线性方程组的根 方法 用近似函数φ(x)的根逼近f(x)的根。

二分法 已知f(a)f(b)<0，设c=(a+b)/2。若f(a)f(c)<0 则根在[a,c]内；若f(a)f(c)>0则根在[b,c]内。 当|f(c)|<ε或|b-a|<ε时，输出c。 迭代步数：O(log2ε)

当| φ’(x)|≤L<1时，|xk+1-α|≤L|xk-α|。 当|xn+1-xn|<ε时，输出 xn。 不动点 当| φ’(x)|≤L<1时，|xk+1-α|≤L|xk-α|。 当|xn+1-xn|<ε时，输出 xn。 迭代步数：O(logLε) Lipschitz 常数 线性收敛

当|f(xk)|<ε或|xk+1-xk|<ε时，输出xk+1。 迭代步数：O(loglogε) Newton法（一阶Taylor展开） 二次收敛 当|f(xk)|<ε或|xk+1-xk|<ε时，输出xk+1。 迭代步数：O(loglogε)

Newton法（p重根情形）

用Newton迭代法求 f(z)=z3−2z+2 的根。当初值分别位于红、蓝、绿色区域时，迭代收敛到三个根。当初值位于黑色区域时，迭代陷入死循环0→1→0。图片引自John Hubbard, Dierk Schleicher, Scott Sutherland, How to find all roots of complex polynomials, Inventiones mathematicae 146, 1-33 (2001).

弦截法 （线性插值） 当|f(xk)|<ε或|xk+1-xk|<ε时，输出xk+1。 迭代步数：O(loglogε)

弦截法的收敛速度

Newton法解非线性方程组 求 　的所有复根 等价于求 x1,…,xn 使 f(t)=(t-x1)…(t-xn)。

其他求根方法 Brent （反插值 x=φ(y)） Halley （二阶Taylor展开） Muller （二次插值） 有理插值……

作业一、习题4第1、3、5、7、8题。 作业二、比较各种求根方法的优缺点。

第5章 解线性方程组的直接法

问题：求解n元线性方程组AX=B。 　　　方法？速度？精度 ？存储？ 下三角方程组 上三角方程组 n(n-1)/2次加减法，n(n+1)/2次乘除法。

Gauss消元法解一般方程组 (2n3+3n2-5n)/6次加减法， (n3+3n2-n)/3次乘除法。

追赶法解三对角方程组 3n-3次加减法，5n-4次乘除法。

线性方程组解的精度 矩阵条件数

Gauss消元法的实质是LU分解 存在性？A的顺序主子式≠0。 唯一性？L1U1=L2U2  L1-1L2=U1U2-1 对角 精确度？A-1b的相对误差≈(L,U,b)的相对误差×cond(L)×cond(U)。

Dolittle分解 Courant分解

全/列/行主元分解 LDLT分解、Cholesky分解

QR分解

SVD分解 Givens旋转 Householder反射

作业一、习题5第1--8题(1)。 作业二、计算100阶Hilbert矩阵H的逆矩阵A以及AH-I的元素平方和。

第6章 解线性方程组的迭代法

Jacobi迭代 Gauss-Seidel迭代

迭代法解线性方程组 A X = B A Xk+1 – B = C (A Xk – B) C 称为 Conditioner，满足 ρ (C)<1　或　||C||<1 通常取 C = I – A Ã-1，其中 Ã≈A，于是 Xk+1 = Xk – Ã-1 (A Xk – B)

Jacobi迭代：Ã = D 定理：A行对角优、或A列对角优  Jacobi迭代收敛。 Gauss-Seidel迭代：Ã = D + L 定理：A行对角优、或A列对角优、或A正定  Gauss-Seidel迭代收敛。

松弛迭代： Ã = w-1D + L 定理：松弛迭代收敛  0<w<2 定理：A正定且0<w<2  松弛迭代收敛 Newton迭代求 A-1： Xk+1 = 2 Xk – Xk A Xk

作业一、145页习题3、6、7。 作业二、用迭代法计算10阶Hilbert矩阵H的逆矩阵A以及AH-I的元素平方和。

第7章 计算矩阵的特征值 和特征向量

问题1：求复方阵的模最大(最小)特征值。 方法：幂法、反幂法 问题2：求复方阵的所有特征值。 方法：QR 迭代 问题3：求Hermite方阵的所有特征值。 方法：Jacobi 方法

幂法 当 A 只有一个模最大的特征值λmax ，并且 x0 与λmax 的特征向量 amax 不正交时 当 A 的模最大的特征值都相同时，以上迭代仍然收敛。

当 A 的模最大的特征值各不相同时，可以选取数 s 使 A - s I 的模最大的特征值只有一个。 当 A 恰有 m 个模最大的特征值时，有 R 的特征值就是 A 的模最大的特征值。

反幂法 当 A 只有一个模最小的特征值λmin ，并且 x0 与λmin 的特征向量 amin 不正交时 计算 A - s I 的模最小的特征值　等价于　计算 A 的最接近 s 的特征值。

QR 迭代 利用 QR 分解，酉相似 A 为上三角。 QR 迭代的本质是幂法 当 时，QR 迭代收敛。 可对 A - s I 作 QR 分解，加速收敛。

Jacobi 方法 通过 Givens 旋转，逐渐减小非对角元。本质是 2 阶 Hermite 方阵的酉相似。 Jacobi 方法具有 2 阶收敛速度。

复矩阵的奇异值分解 A = UΣV 一般方法 AH A = VHΣ2 V 或 A AH = UΣ2 UH QR 迭代 Jacobi 方法 计算 2 阶方阵的 SVD

作业一、167页习题1(3)、2(2)、3(4)。 作业二、计算10阶Hilbert矩阵的正交相似标准形以及过渡矩阵。

第8章 常微分方程数值解

问题：求解一阶常微分方程的初值问题： 解法：化微分方程为积分方程

Euler折线法 向前Euler公式 向后Euler公式 Picard迭代 中心Euler公式 梯形公式 改进的 Euler公式

Runge-Kutta方法 选取 {xi, cij} 使 yr 有最高精度 p，即 r 1, 2, 3, 4 5, 6, 7 8, 9 10, 11, ... p p = r p = r-1 p = r-2 p ≤ r-2

Runge-Kutta方法的误差估计 设　　　满足Lipschitz条件 设　　满足 初值误差 截断误差 整体误差

线性多步法 (*) 其中 φ(t) 为 f(t,y(t)) 的 q 次插值多项式。 当 xn,…,xn-q 为插值节点时，(*) 称为显式Adams公式。 当 xn+1,…,xn+1-q 为插值节点时，(*) 称为隐式Adams公式。

一阶常微分方程组的初值问题： 解法：同一阶常微分方程的初值问题。

高阶常微分方程的初值问题 解法：令 化方程为

单步法 (*) 收敛性 稳定性 将(*)应用于方程y’=y，得 yn+1=E(h)yn。 当|E(h)|<1时，称(*)绝对稳定的。 称{复数 h : |E(h)|<1}为绝对稳定区域。称{实数 h : |E(h)|<1}为绝对稳定区间。

作业、191-192页习题1--6。

插值 拟合 微积分 非线性方程 线性方程组 特征值问题 常微分方程 插值公式、截断误差、计算复杂度、差商、Runge现象 最小二乘法，SVD分解 微分公式、积分公式、截断误差、外推法、代数精度 解法、计算复杂度、收敛阶数 直接法、迭代法、误差、计算复杂度、收敛条件、LU分解、范数、条件数 算法、计算复杂度、收敛条件、QR分解 解法、误差估计、稳定性、显式、隐式

结 束