数值计算方法

第一章——数值计算中的误差

误差来源

观测误差、模型误差、截断误差（方法误差）、舍入误差

绝对误差、相对误差和有效数字

设准确值为$x$，计算近似值为$x^{\ast }$
绝对误差：$\epsilon =|x-x^{\ast }|$
相对误差：${\epsilon }_r=\frac{|x-x^{\ast }|}{|x|}\approx \frac{|x-x^{\ast }|}{|x^{\ast }|}$
有效数字：若$|x-x^{\ast }|\le 0.5\times 10^{-n}$，则称$x^{\ast }$精确到小数点后n位，并把第一个非零位到这一位的所有数字称为有效数字

误差传播

加减：绝对误差累计
乘除：相对误差累计

数值计算时注意的问题

1. 相近数相减
2. 大数吃小数（相差大的数加减）
3. 分母绝对值太小

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

解决问题

小型线性方程组求解。

消去法

• 高斯（Gauss）消去法：只做行加
• 列主元素法：先通过换行使得主元最大，再用高斯消去法求解
• 全主元素法：先通过换列和换行使得主元最大，再高斯消去法求解

分解法

方法	格式	要求
杜利特尔（Doolittle）分解	$A=LU$	$A$对称正定
平方根法（Cholesky）	$A=L{L}^{\top }$	$A$对称正定
改进平分根法	$A=LD{L}^{\top }$	$A$对称正定
追赶法	$A=LU$	三对角矩阵、严格对角占优

误差估计

向量范数

• 1范数：$||x|{|}_1={\sum }_{i=1}^n|x_i|$
• 2范数：$||x|{|}_2=({\sum }_{i=1}^n{x}_i^2{)}^{\frac{1}{2}}$
• 无穷范数：$||x|{|}_{\infty }={\max }_{1\le i\le n}|x_i|$
• p范数：$||x|{|}_p=({\sum }_{i=1}^n{x}_i^p{)}^{\frac{1}{p}}$

矩阵范数

• 1（列）范数：$||A|{|}_1=ma{x}_{1\le j\le n}({\sum }_{i=0}^n|a_{ij}|)$
• 无穷（行）范数：$||A|{|}_{\infty }=ma{x}_{1\le i\le n}({\sum }_{j=0}^n|a_{ij}|)$
• 2范数：$||A|{|}_2=\sqrt{{\lambda }_{max}(A^{\top }A)}$
• F范数：$||A|{|}_F=({\sum }_{i=0}^n{\sum }_{j=0}^n{a}_{ij}^2{)}^{\frac{1}{2}}$

谱半径$\rho (A)=max({\lambda }_i)$, $\rho (A)\le 范数$

条件数

$$cond(A)=||A||\cdot ||A^{-1}||$$
条件越大数大，矩阵越病态，系数矩阵或右侧常数项有微小扰动，会造成巨大的误差。

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

解决问题

利用迭代法求解大型方程组。

雅可比（Jacobi）迭代法

迭代公式：
$$\{\begin{array}{rl}{x}_1^{(k+1)}& =(b_1-a_{12}{x}_2^{(k)}-a_{13}{x}_3^{(k)})/a_{11}\\ x_2^{(k+1)}& =(b_2-a_{21}{x}_1^{(k)}-a_{22}{x}_2^{(k)})/a_{22}\\ x_3^{(k+1)}& =(b_3-a_{31}{x}_1^{(k)}-a_{32}{x}_2^{(k)})/a_{33}\end{array}$$
迭代矩阵：
$$\begin{gathered} A=D-(L+U) \\ M=D^{-1}(L+U) \end{gathered}$$

高斯赛德尔（Gauss-Sedel）迭代法

迭代公式：
$$\{\begin{array}{rl}{x}_1^{(k+1)}& =(b_1-a_{12}{x}_2^{(k)}-a_{13}{x}_3^{(k)})/a_{11}\\ x_2^{(k+1)}& =(b_2-a_{21}{x}_1^{(k+1)}-a_{22}{x}_2^{(k)})/a_{22}\\ x_3^{(k+1)}& =(b_3-a_{31}{x}_1^{(k+1)}-a_{32}{x}_2^{(k+1)})/a_{33}\end{array}$$
迭代矩阵：
$$\begin{gathered} A=(D-L)+U \\ M=(D-L{)}^{-1}U \end{gathered}$$

超松弛法

迭代公式：
$$\{\begin{array}{rl}{x}_1^{(k+1)}& =(1-w)x_1^{(k)}+w(b_1-a_{12}{x}_2^{(k)}-a_{13}{x}_3^{(k)})/a_{11}\\ x_2^{(k+1)}& =(1-w)x_2^{(k)}+w(b_2-a_{21}{x}_1^{(k+1)}-a_{22}{x}_2^{(k)})/a_{22}\\ x_3^{(k+1)}& =(1-w)x_3^{(k)}+w(b_3-a_{31}{x}_1^{(k+1)}-a_{32}{x}_2^{(k+1)})/a_{33}\end{array}$$

收敛条件

根据系数矩阵A判断（充分条件）

• 雅可比迭代法：A严格对角占优 $\cup$ A弱对角占优且不可约
• 高斯赛德尔：A严格对角占优 $\cup$ A弱对角占优且不可约 $\cup$ A对称正定
• 超松弛：A严格对角占优$(0<w<1)$ $\cup$ A对称正定$(0<w<2)$

根据迭代矩阵判断（充要条件）

谱半径$\rho (M)<1$是迭代收敛的充要条件
推论：若迭代矩阵范数$||M||<1$，则迭代收敛

误差估计

估计迭代次数：
$$||x^{(k)}-x^{\ast }||\le \frac{||M|{|}^k}{1-||M||}||x^{(1)}-x^{(0)}||$$
停止条件：
$$||x^{(k)}-x^{\ast }||\le \frac{||M||}{1-||M||}||x^{(k)}-x^{(k-1)}||$$

第四章——矩阵的特征值与特征向量的数值计算

解决问题

快速求解矩阵的某一特征值和特征向量。

幂法

用于迭代求解矩阵A的最大特征值和特征向量。

核心思想

设矩阵$A_{n\times n}$的特征向量为$u_1,u_2,...,u_n$且其对应的特征值满足$|{\lambda }_1|>|{\lambda }_2|\ge |{\lambda }_3|...\ge |{\lambda }_n|$。任取初始化向量x，则x可以用A的特征向量表示：
$$x=a_1{u}_1+a_2{u}_2+...+a_n{u}_n$$
将x不断左乘A，得得到：
$$\begin{array}{rl}{A}^{k}x& =A^k{a}_1{u}_1+A^k{a}_2{u}_2+...+A^k{a}_n{u}_n\\ & ={\lambda }_1^k{a}_1{u}_1+{\lambda }_2^k{a}_2{u}_2+...+{\lambda }_n^k{a}_n{u}_n\\ & ={\lambda }_1^k(a_1{u}_1+\epsilon )\\ & \approx {\lambda }_1^k{a}_1{u}_1\end{array}$$
所以，将x不断左乘A，可以将其转换成最大的特征向量且每次增大的倍数就为对应的特性值。

具体算法

1. 输入矩阵A，误差$ϵ$，最大迭代次数N，初始化向量x，初始特征值$\lambda$
2. 归一化初始向量：$y=\frac{x}{max(x)}$
3. 迭代：$x=Ay$，$\beta =max(x)$，$y=\frac{x}{\beta }$
4. 判断终止条件：是否满足最大迭代次数 $\cup$ $|\beta -\lambda |<ϵ$，不满足 ⇒ 3；满足 ⇒ 5
5. 特征值为$\beta$，特征向量为$y$

反幂法

用于求解最小的特征值和特征向量 ⇒ 求接近${\lambda }^{\ast }$的特征值和特征向量

具体算法

1. 输入矩阵A，误差$ϵ$，最大迭代次数N，初始化向量x，目标特征值${\lambda }^{\ast }$，初始特征值为$\lambda =0$
2. 设$B=A-{\lambda }^{\ast }I$，求$B$的LU分解
3. 归一化初始向量：$y=\frac{x}{max(x)}$
4. 迭代：用求解$LUx=y$代替$x=B^{-1}y$，$\beta =max(x)$，$y=\frac{x}{\beta }$
5. 判断终止条件：是否满足最大迭代次数 $\cup$ $|\beta -\lambda |<ϵ$，不满足 ⇒ 4；满足 ⇒ 6
6. 特征值为${\lambda }^{\ast }+\frac{1}{\beta }$，特征向量为$y$

第五章——插值法

解决问题

给出离散数值点，用多项式函数做局部拟合。

拉格朗日（Lagrange）插值

结构规整，由方程组导出的函数通解形式：
$$L_n(x)=\sum _{i=0}^n\frac{{\prod }_{j=0,j\ne i}^n(x-x_j)}{{\prod }_{j=0,j\ne i}^n(x_i-x_j)}\times y_i$$
n阶式需要n+1个节点，例如2阶式子如下：
$$L_2(x)=\frac{(x-x_1)(x-x_2)}{(x_0-x_1)(x_0-x_2)}{y}_0+\frac{(x-x_0)(x-x_2)}{(x_1-x_0)(x_1-x_2)}{y}_1+\frac{(x-x_0)(x-x_1)}{(x_2-x_0)(x_2-x_1)}{y}_2$$

牛顿（Newton）插值法

可以逐次添加，可以使用先前的计算结果：
$$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]+...+(x-x_0)(x-x_1)...(x-x_{n-1})f[x_0,x_1,...,x_n]$$
n阶式需要n+1个节点，通常在求解时会构造差商表进行求解。

埃尔米特（Hermite）插值法

利用导数信息做插值，常用降解法求解：

1. 利用牛顿法或拉格朗日法求解普通插值
2. 设一般式：$H(x)=N_n(x)+f(x)(x-x_0)(x-x_1)...(x-x_n)$，在代入导数信息求解

误差估计

对于n+1个节点，可以构造n阶多项式，其误差为（上述三种方法都适用）：
$$R_n(x)=f(x)-{\varphi }_n(x)=\frac{f^{(n+1)}(\xi )}{(n+1)!}\prod _{i=0}^n(x-x_i)$$

第六章——函数逼近

解决问题

给出一些列点，拟合函数表达式的参数

多项式最小二乘法

记n个点的m阶多项式拟合误差为：$F(a_0,a_1...a_m)={\sum }_{i=0}^n(y_i-{\sum }_{k=0}^m{a}_k{x}_i^k{)}^2$
令各个$a_i$求偏导为0：$\frac{dF}{d{a}_p}={\sum }_{i=0}^n[2(y_i-{\sum }_{k=0}^m{a}_k{x}_i^k)x_i^p]=0$
由此可以得到m+1阶方程组：${\sum }_{i=0}^n{\sum }_{k=0}^m{a}_k{x}_i^{k+p}={\sum }_{i=0}^n{y}_i{x}_i^p$,当m=2时，方程简记为：
$$\left[\begin{array}{ccc}{S}_n& S_x& S_{x^2}\\ S_x& S_{x^2}& S_{x^3}\\ S_{x^2}& S_{x^3}& S_{x^4}\end{array}\right]\left[\begin{array}{c}{a}_0\\ a_1\\ a_2\end{array}\right]=\left[\begin{array}{c}{S}_y\\ S_{yx}\\ S_{y{x}^2}\end{array}\right]$$
其中, $S_x={\sum }_{i=0}^n{x}_i$ ，其它表达式同理

补充一些变换

求$y=b{e}^{ax}$时，可以先两边取对数得到$\ln y=\ln b+ax$，用$Y=a_0+a_{1}x$拟合
求$y=a\frac{x}{x+b}$时，先去倒数得到$\frac{1}{y}=\frac{b}{ax}+\frac{1}{a}$, 用$Y=a_0+a_{1}X$拟合，再代会x, y

第七章——数值积分和数值微分

解决问题

利用数值计算的方法求某点处的微分和某个区间内的导数

数值微分

名称	公式	误差
一阶向前差商	$f^{\prime }(x)=\frac{f(x+h)-f(x)}{h}$	$\frac{h}{2}{f}^{\prime \prime }(\xi )$
一阶向后差商	$f^{\prime }(x)=\frac{f(x)-f(x-h)}{h}$	$\frac{h}{2}{f}^{\prime \prime }(\xi )$
一阶中心差商	$f^{\prime }(x)=\frac{f(x+h)-f(x-h)}{2h}$	$\frac{h^2}{6}{f}^{\prime \prime \prime }(\xi )$
二阶差商	$f^{\prime }(x)=\frac{f(x+h)-2f(x)+f(x-h)}{h^2}$	$\frac{h^2}{12}{f}^{(4)}(\xi )$
三点式1	$f^{\prime }(x_0)=\frac{-3f(x0)+4f(x_1)+f(x_2)}{2h}$	$\frac{h^2}{3}{f}^{(3)}(\xi )$
三点式2	$f^{\prime }(x_1)=\frac{f(x_2)-f(x_1)}{2h}$	$\frac{h^2}{6}{f}^{(3)}(\xi )$
三点式3	$f^{\prime }(x_2)=\frac{f(x_0)-4f(x_1)+3f(x_2)}{2h}$	$\frac{h^2}{3}{f}^{(3)}(\xi )$

牛顿-科特斯（Newton-Cotes）积分

按比例构造函数值的线性组合求出积分。

阶级	别名	比例	误差
1	梯形公式	1 1	$\frac{(b-a{)}^3}{12}{f}^{\prime \prime }(\xi )$
2	辛普森(Simpsen)公式	1 4 1	$\frac{(b-a{)}^5}{2880}{f}^{(4)}(\xi )$

衡量准确性

利用代数精度衡量准确性：某式子能准备积分多少阶代数，则具有多少代数精度。

复化求积公式

通过减小积分区间来减小积分的误差。

名称	公式	误差
梯形公式	$\frac{h}{2}[f(a)+f(b)+2{\sum }_{i=1}^{n-1}f(x_i)]$	$\frac{(b-a)}{12}{h}^2{f}^{\prime \prime }(\xi )$
辛普森公式	$\frac{h}{3}[f(a)+f(b)+2\sum f(x_{偶})+4\sum f(x_{奇})]$	$\frac{b-a}{180}{h}^4{f}^{(4)}(\xi )$

逐次半分法

合理缩小积分粒度，求得更准确的结果。
逐次半分梯形：
$$T_{2^{m+1}}=\frac{1}{2}{T}_{2^m}+h_{m+1}\sum f(x_{新})$$

龙贝格（Romberg）积分法

在逐次半分法的基础上增加外推法（4倍）获得更准确的结果。

高斯（Gauss）型求积公式

具有2n-1阶代数精确度的求积式称为高斯型求积公式。

第八章——非线性方程的解法

解决的问题

$$f(x)=0$$

对分区间法

若f(x)在[a, b]上连续，且有$f(a)f(b)<0$，则可以依照零点定理通过每次取中点的方法迭代缩小范围。

迭代法

将原方程转换成$x_{n+1}=\varphi (x_n)$的形式迭代求解。注意不是所有的格式都收敛，可以根据下面定理判断。

判断收敛性

压缩印象原理

1. 对于$x\in [a,b]$，有$\varphi (x)\in [a,b]$
2. 对于$\forall x,y\in [a,b]$，$|\varphi (x)-\varphi (y)|\le L(L<1)$

则$\varphi (x)$在[a, b]中收敛到唯一解$x^{\ast }$且误差满足$|x^{\ast }-x_n|\le \frac{L^n}{1-L}|x_1-x_0|$

推论

1. 对于$x\in [a,b]$，有$\varphi (x)\in [a,b]$
2. 对于$\forall x,y\in [a,b]$，$\varphi (x{)}^{\prime }\le L(L<1)$

则$\varphi (x)$在[a, b]中收敛到唯一解$x^{\ast }$且误差满足$|x^{\ast }-x_n|\le \frac{L^n}{1-L}|x_1-x_0|$

局部收敛定理

若$\varphi (x)$在$x^{\ast }$的某邻域$O(x^{\ast },\theta )$中连续可微，且$|\varphi (x{)}^{\prime }|<1$，则$\varphi (x)$在局部收敛。

收敛速度

若${\lim }_{n\to \infty }\frac{|x_{n+1}-x^{\ast }|}{|x_n-x^{\ast }{|}^r}=a$，则是r阶收敛，可以利用洛必达法则对$\varphi (x)$反复求导，求收敛阶。

牛顿法

设$f(x)=0$, 提供一种构造方法：
$$x_{n+1}=x_n-\frac{f(x)}{f(x{)}^{\prime }}$$

牛顿弦截法

用差商代替导数，得到迭代公式:
$$x_{n+1}=x_n-\frac{f(x_n)}{f(x_n)-f(x_{n-1})}(x_n-x_{n-1})$$

第九章——常微分方程的数值解法

解决的问题

$$\{\begin{array}{rl}\frac{dy}{dx}& =f(x,y)\\ y(a)& =y_0\end{array}$$

欧拉（Euler）法

将该问题直接作为一个积分问题，从$x=a$开始按定长$h$积分：
$$y_{n+1}=y_n+h\times f(x_n,y_n)$$

例题

改进欧拉法

先利用Euler法求得$y_{n+1}$的值，由$(x_{n+1},y_{n+1})$再做一次Euler进行修正：
$$\{\begin{array}{rl}{y}_p& =y_n+h\times f(x_n,y_n)\\ y_q& =y_n+h\times f(x_{n+1},y_q)\\ y_{n+1}& =\frac{y_p+y_q}{2}\end{array}$$

例题