数值分析（一）

基础知识

误差类型及传播

误差类型	定义
绝对误差	$e(x^{\ast })=x^{\ast }-x$
绝对误差限	$\mid e(x^{\ast })\mid \le \epsilon (x^{\ast })$
相对误差	$e_r(x^) = \cfrac{x^ - x}{x} $
相对误差限	$\mid e_r(x^{\ast })\mid \le {\epsilon }_r(x^{\ast })$

• 设$x$的近似值$x^{\ast }$有如下标准形式
  $$x^{\ast }=\pm 10^m\times 0.a_1{a}_2\cdots a_n\cdots a_p$$
  其中$m$为整数， { a i } ∈ { 0 , 1 , 2 , ⋯   , 9 } \{a_i\} \in \{0,1,2,\cdots,9\} {ai​}∈{0,1,2,⋯,9}且$a_1\ne 0,p\ge n$。如果$n$满足
  $$\mid e^{\ast }\mid =\mid x^{\ast }-x\mid \le \frac{1}{2}\times 10^{m-n}$$
  的最大正整数，则有效数字与相对误差的关系

  • 若$x^{\ast }$具有$n$位有效数字，则(有省略，不要用)
    $$\mid e_r^{\ast }\mid \le \frac{1}{2a_1}\times 10^{1-n}$$

  • 若
    $$\mid e_r^{\ast }\mid \le \frac{1}{2(a_1+1)}\times 10^{1-n}$$
    则$x^{\ast }$至少具有$n$位有效数字

  • 近似数的绝对误差限是末尾数位的一半

• 数值运算的误差估计

  • 一元函数误差估计
    $$\begin{array}{rl}& e(y^{\ast })\approx f^{\prime }(x^{\ast })e(x^{\ast })\\ & e_r(y^{\ast })\approx \frac{x^{\ast }{f}^{\prime }(x^{\ast })}{f(x^{\ast })}{e}_r(x^{\ast })\end{array}$$

  • 多元函数误差估计
    $$\begin{array}{rl}& e(y^{\ast })\approx \sum _{i=1}^n{f}_i^{\prime }(x_1^{\ast },\cdots ,x_n^{\ast })e(x_i^{\ast })\\ & e_r(y^{\ast })\approx \sum _{i=1}^n\frac{x_i^{\ast }}{y^{\ast }}{f}_i^{\prime }(x_1^{\ast },\cdots ,x_n^{\ast })e_r(x_i^{\ast })\end{array}$$

  • 算术运算的误差估计
    $$\begin{array}{rl}& \epsilon (x_1^{\ast }\pm x_2^{\ast })\approx \epsilon (x_1^{\ast })+\epsilon (x_2^{\ast })\\ & \epsilon (x_1^{\ast }{x}_2^{\ast })\approx \mid x_2^{\ast }\mid \epsilon (x_1^{\ast })+\mid x_1\mid \epsilon (x_2^{\ast })\\ & \epsilon \left(\frac{x_1^{\ast }}{x_2^{\ast }}\right)\approx \frac{\mid x_2^{\ast }\mid \epsilon (x_1^{\ast })+\mid x_1\mid \epsilon (x_2^{\ast })}{\mid x_2^{\ast }{\mid }^2}\end{array}$$

    ε r ( x 1 ∗ ± x 2 ∗ ) ≤ max ⁡ { ε r ( x 1 ∗ ) , ε r ( x 2 ∗ ) } ε r ( x 1 ∗ x 2 ∗ ) ≈ ε r ( x 1 ∗ ) + ε r ( x 2 ∗ ) ε r ( x 1 ∗ x 2 ∗ ) ≈ ε r ( x 1 ∗ ) + ε r ( x 2 ∗ ) \begin{aligned} &\varepsilon_r(x_1^* \pm x_2^*) \leq \max\{\varepsilon_r(x_1^*),\varepsilon_r(x_2^*)\}\\ &\varepsilon_r(x_1^* x_2^*) \approx \varepsilon_r(x_1^*) + \varepsilon_r(x_2^*) \\ &\varepsilon_r\Big(\cfrac{x_1^*}{ x_2^*}\Big) \approx \varepsilon_r(x_1^*)+ \varepsilon_r(x_2^*) \end{aligned} ​εr​(x1∗​±x2∗​)≤max{εr​(x1∗​),εr​(x2∗​)}εr​(x1∗​x2∗​)≈εr​(x1∗​)+εr​(x2∗​)εr​(x2∗​x1∗​​)≈εr​(x1∗​)+εr​(x2∗​)​

误差定性分析及数值运算中的若干原则

病态问题与条件数

• 病态问题：对一个数值问题，若由于问题本身而使计算结果绝对误差或相对误差的绝对值很大，这种问题称为绝对误差意义或相对意义下的病态问题

• 条件数

  条件数较大，则表示自变量的微小扰动将引起函数值较大的变化

  • 绝对意义下计算函数值$f(x)$的条件数：$Con{d}_a(f)=\mid f^{\prime }(x)\mid$
  • 相对意义下计算函数值$f(x)$的条件数：$Con{d}_r(f)=\mid x\frac{f^{\prime }(x)}{f(x)}\mid$

算法的数值稳定性

一个算法，如果初始数据微小的误差仅使最终结果产生微小的误差，或在运算过程中舍入误差在一定条件下能够得到控制，则称该算法数值问题，否则不稳定

数值运算中若干原则

• 避免大数吃掉小数
• 避免用绝对值很小的数作除数
• 避免两个相近数相减
• 减小计算工作量

非线性方程求根

• 如果有$\alpha$使$f(\alpha )=0$，则称$\alpha$为方程$f(x)=0$的根或函数$f(x)$的零点，如果函数$f(x)$可分解为
  $$f(x)=(x-\alpha {)}^{m}g(x)$$
  且$g(\alpha )\ne 0$，$m$为正整数，则当$m=1$时，称$\alpha$为$f(x)=0$的单根或$f(x)$的单零点；当$m\ge 2$时，称$\alpha$为$f(x)=0$的$m$重根或$f(x)$的$m$重零点
• 若$f(x)$有$m$阶连续导数，则$\alpha$是$f(x)$的$m$重零点的充要条件为
  $$f(\alpha )=0,f^{\prime }(\alpha )=0,\cdots ,f^{(m-1)}(\alpha )=0,f^{(m)}(\alpha )\ne 0$$

二分法

设$\alpha$为方程$f(x)=0$在区间$[a,b]$内的唯一根，其中$f(x)$为满足$f(a)f(b)<0$的连续函数，则二分法计算过程中第$k$个区间$[a_k,b_k]$的中点$x_k$满足不等式
$$\mid x_k-\alpha \mid \le \frac{1}{2}(b_k-a_k)=\frac{1}{2^k}(b-a)$$
且二分法产生的序列 { x k } k = 1 ∞ \{x_k\}_{k=1}^\infty {xk​}k=1∞​收敛于$\alpha$

不动点迭代的基本理论

将$f(x)=0$改写为同解方程组$x=\phi (x)$的形式。如果$f(\alpha )=0$，则$\alpha =\phi (\alpha )$，称$\alpha$为函数$\phi (x)$的不动点，选取初始近似值$x_0$，并构造迭代格式
$$x_{k+1}=\phi (x_k)$$
可产生序列 { x k } k = 0 ∞ \{x_k\}_{k=0}^\infty {xk​}k=0∞​，上述迭代式称为不动点迭代，$\phi (x)$称为迭代函数。

• 全局收敛性
  对于方程$x=\phi (x)$，若${\phi }^{\prime }(x)$在$[a,b]$上连续，且

  • 封闭性：$x\in [a,b],\phi (x)\in [a,b]$
  • 压缩性：$\exists L\in (0,1),x\in [a,b],\mid {\phi }^{\prime }(x)\mid \le L<1$

  $⟹$

  • $\phi (x)$在$[a,b]$上有唯一的不动点$\alpha$
  • $\forall x_0\in [a,b],x_{k+1}=\phi (x_k)$收敛，且${\lim }_{k\to \infty }{x}_k=\alpha$
  • 序列 { x k } k = 0 ∞ \{x_k\}_{k=0}^\infty {xk​}k=0∞​有误差估计式
    $$\mid x_k-\alpha \mid \le \frac{L}{1-L}\mid x_k-x_{k-1}\mid \text{\hspace{1em}(误差后验估计式)}$$
    $$\mid x_k-\alpha \mid \le \frac{L^k}{1-L}\mid x_1-x_0\mid \text{\hspace{1em}(误差先验估计式)}$$
    $$\underset{k\to \infty }{\lim }\frac{x_k-\alpha }{x_{k-1}-\alpha }={\phi }^{\prime }(\alpha )\text{\hspace{1em}(误差渐进表达式)}$$
    设函数$\phi (x)$在区间$[a,b]$内有不动点$\alpha$，且当$x\in [a,b]$时，$\mid {\phi }^{\prime }(x)\mid \ge 1$，则对$\forall x_0\in [a,b],x_0\ne \alpha$,迭代格式$x_{k+1}=\phi (x_k)$发散

• 局部收敛性
  设函数$\phi (x)$有不动点$\alpha$，如果存在$\alpha$的一个邻域$N(\alpha )=[\alpha -\sigma ,\alpha +\sigma ],\sigma >0$，使迭代格式$x_{k+1}=\phi (x_k)$对$\forall x_0\in N(\alpha )$均收敛，则称该迭代格式具有局部收敛性
  设函数$\phi (x)$有不动点$\alpha$，${\phi }^{\prime }(x)$在$\alpha$连续，则

  • $\mid {\phi }^{\prime }(\alpha )\mid <1\Rightarrow$迭代格式收敛
  • $\mid {\phi }^{\prime }(\alpha )\mid >1\Rightarrow$迭代格式发散
    设序列 { x k } k = 0 ∞ \{x_k\}_{k=0}^\infty {xk​}k=0∞​收敛到$\alpha$。令$k$次的迭代误差$e_k=x_k-\alpha ,x_k\ne \alpha$。如果存在实数$p\ge 1$及非零常数$c>0$，使得
    $$\underset{k\to \infty }{\lim }\frac{\mid e_{k+1}\mid }{\mid e_k{\mid }^p}=c$$
    则称序列 { x k } k = 0 ∞ \{x_k\}_{k=0}^\infty {xk​}k=0∞​为$p$阶收敛
  • 若$\phi (x)$在不动点$\alpha$的领域内有$p$阶连续导数，则当初值$x_0$取得充分靠近$\alpha$时，迭代格式$x_{k+1}=\phi (x_k)$为$p$阶收敛的充要条件是
    $$\phi (\alpha )=\alpha ,{\phi }^{\prime }(\alpha )={\phi }^{\prime \prime }(\alpha )=\cdots ={\phi }^{(p-1)}(\alpha )=0,{\phi }^{(p)}(\alpha )\ne 0$$
    且有
    $$\underset{k\to \infty }{\lim }\frac{e_{k+1}}{e_k^p}=\frac{1}{p!}{\phi }^{(p)}(\alpha )\ne 0$$
    如果$p=1$时，要求$0<\mid {\phi }^{\prime }(\alpha )\mid <1$

不动点迭代的加速

$$\begin{array}{rl}& x_{k+1}-\alpha =\phi (x_k)-\phi (\alpha )={\phi }^{\prime }({\xi }_k)(x_k-\alpha ),{\xi }_k\in [\alpha ,x_k]\\ & x_{k+2}-\alpha =\phi (x_{k+1})-\phi (\alpha )={\phi }^{\prime }({\xi }_{k+1})(x_{k+1}-\alpha ),{\xi }_{k+1}\in [\alpha ,x_{k+1}]\end{array}$$
令${\phi }^{\prime }({\xi }_k)\approx {\phi }^{\prime }({\xi }_{k+1})$,则有
$$\frac{x_{k+2}-\alpha }{x_{k+1}-\alpha }\approx \frac{x_{k+1}-\alpha }{x_k-\alpha }$$
解出
$$x_{k+1}=\psi (x_k)=\alpha \approx x_k-\frac{(x_{k+1}-x_k{)}^2}{x_{k+2}-2x_{k+1}+x_k}=x_k-\frac{[\phi (x_k)-x_k{]}^2}{\phi (\phi (x_k))-2\phi (x_k)+x_k}$$
称其为$\mathrm{Steffensen}$迭代格式

• $\psi (\alpha )=\alpha \Rightarrow \phi (\alpha )=\alpha$；反之$\varphi(\alpha) = \alpha $，且$\varphi^\prime( x )$在 $\alpha$的邻域内连续，${\phi }^{\prime }(\alpha )\ne 1$,则$\psi (\alpha )=\alpha$
• 若$\phi (x)$在其不动点$\alpha$的邻域内有二阶连续导数，且${\phi }^{\prime }(\alpha )=A\ne 1,\ne 0$，则$\mathrm{Steffensen}$迭代局部收敛
• 若${\phi }^{\prime \prime }(\alpha )\ne 0$，则$\mathrm{Steffensen}$迭代平方收敛；若${\phi }^{\prime \prime }(\alpha )=0$，则$\mathrm{Steffensen}$迭代超平方收敛【可将线性收敛或不收敛的式子平方收敛】

Newton迭代

为了构造二阶收敛速度的迭代格式，令$\phi (x)=x+h(x)f(x)$，通过满足二阶收敛的条件求解设定函数。或者将原函数$f(x)$进行泰勒展开，得到二阶收敛格式
$$\phi (x)=x-\frac{f(x)}{f^{\prime }(x)}$$

• 局部收敛性
  设$\alpha$是方程$f(x)=0$的根，$f^{\prime }(x)$在$\alpha$的邻域内连续，则$\text{Newton}$迭代局部收敛。且求单根时，$\text{Newton}$迭代至少二阶收敛；而求重根时，$\text{Newton}$迭代只有一阶收敛
  当$f(x)=(x-\alpha {)}^{m}g(x),g(\alpha )\ne 0$时，有
  $${\phi }^{\prime }(\alpha )=1-\frac{1}{m}$$

• 全局收敛性判定$\text{I}$
  给定$f(x)=0$，如果$f^{\prime \prime}(x) $在$[a,b]$上连续，且

  • $f(a)f(b)<0$
  • $\forall x\in [a,b],f^{\prime }(x)\ne 0,f^{\prime \prime }(x)\ne 0$
  • 取$x_0\in [a,b],f(x_0)f^{\prime \prime }(x_0)>0$
    则$f(x)=0$在$[a,b]$内有唯一的根$\alpha$，且由初值$x_0$按$\mathrm{Newton}$迭代格式产生的序列收敛于$\alpha$，并有
    $$\underset{k\to \infty }{\lim }\frac{x_{k+1}-\alpha }{(x_k-\alpha {)}^2}=\frac{f^{\prime \prime }(\alpha )}{2f^{\prime }(\alpha )}$$
    即迭代二阶收敛

• 全局收敛性判定$\text{II}$
  给定$f(x)=0$，如果$f^{\prime \prime }(x)$在$[a,b]$上连续，且

  • $f(a)f(b)<0$
  • $\forall x\in [a,b],f^{\prime }(x)\ne 0,f^{\prime \prime }(x)\ne 0$
  • $\mid \frac{f(a)}{f^{\prime }(a)}\mid \le b-a\mid ,\frac{f(b)}{f^{\prime }(b)}\mid \le b-a$

  则$\forall x_0\in [a,b]$，由$\text{Newton}$迭代格式产生的序列 { x k } k = 0 ∞ \{x_k\}_{k=0}^\infty {xk​}k=0∞​收敛于$f(x)=0$的根$\alpha$

Newton迭代的变形

求重根的修正Newton法

当$\alpha$为$f(x)=0$的$m\ge 0$重根时，将$\text{Newton}$迭代法修正为
$$x_{k+1}=x_k-m\frac{f(x_k)}{f^{\prime }(x_k)}\text{\hspace{1em}(修正的Newton迭代格式)}$$
该迭代方程局部收敛，且至少二阶收敛。
由于$m$一般事先不知道，因此取
$$\mu (x)=\frac{f(x)}{f^{\prime }(x)}$$
当$\alpha$为$f(x)=0$的$m$重根时，$\alpha$是$\mu (x)$的单重零点，则$\text{Newton}$修正为
$$\phi (x)=x-\frac{\mu (x)}{{\mu }^{\prime }(x)}=x-\frac{f(x)f^{\prime }(x)}{(f^{\prime }(x){)}^2-f(x)f^{\prime \prime }(x)}$$
可知其局部收敛，且至少二阶收敛
零点的重数$m$可通过以下方式估计
由误差渐进表达式和牛顿法迭代式与$m$的关系可得
$$\underset{k\to \infty }{\lim }\frac{x_{k+1}-\alpha }{x_k-\alpha }=1-\frac{1}{m}\Rightarrow m\approx \frac{\alpha -x_k}{x_{k+1}-x_k}$$

Newton下山法

由于$\text{Newton}$迭代法对初值的选择要求比较苛刻，为了放宽初值的选择范围，将$\mathrm{Newton}$法的计算结果${\overline{x}}_{k+1}=x_k-\frac{f(x_k)}{f^{\prime }(x_k)}$与前一步的近似值$x_k$进行加权后得到新的近似，即
$$\begin{array}{rl}{x}_{k+1}& =(1-{\lambda }_s)x_k+{\lambda }_s{\overline{x}}_{k+1}\\ & =(1-{\lambda }_s)x_k+{\lambda }_s[x_k-\frac{f(x_k)}{f^{\prime }(x_k)}]\\ & =x_k-{\lambda }_s\frac{f(x_k)}{f^{\prime }(x_k)},s=0,1,\cdots \end{array}$$
${\lambda }_s$称为下山因子,通过选择适当下山因子，能使得产生的迭代序列 { x k } k = 0 ∞ \{x_k\}_{k=0}^\infty {xk​}k=0∞​满足
$$\mid f(x_{k+1})\mid <\mid f(x_k)\mid$$
则该序列是一个下山序列。
$\mathrm{Newton}$下山法的具体算法如下：
计算开始前需要先给出根的绝对误差限${\epsilon }_x$，下山因子的下界${\epsilon }_{{\lambda }_s}$，$\mid f(x)\mid$的下界${\epsilon }_f$，适当的增量$\delta$以及最大迭代次数$N$

1. 选取初值$x_0$
2. 令下山因子${\lambda }_s=1$
3. 利用$\mathrm{Newton}$下山公式由$x_k$计算$x_{k+1}$，并计算$f(x_{k+1})$
4. 若$k>N$，则转到$6$，否则执行$5$
5. 判断$\mid f(x_{k+1})\mid <\mid f(x_k)\mid$是否成立
   1. 当$\mid f(x_{k+1})\mid <\mid f(x_K)\mid$成立，分情况执行
      a. 若$\mid x_{k+1}-x_k\mid <{\epsilon }_x$，则取$\alpha \approx x_{k+1}$，输出结果
      b. 若$\mid x_{k+1}-x_k\mid \ge {\epsilon }_x$，则将$x_{k+1}$作为新的$x_k$，$k$增加$1$，转向$2$
   2. 当$\mid f(x_{k+1}) \mid \geq \mid f(x_k) \mid $成立，分情况执行a.若$\mid f(x_{k+1}) \mid < \varepsilon_f$，则取$\alpha \approx x_{k+1}$，输出结果b.若$\mid f(x_{k+1} \mid \geq \varepsilon_f$，分情况执行
      • 若${\lambda }_s>{\epsilon }_{{\lambda }_s}$，则以$\frac{{\lambda }_s}{2}$代替${\lambda }_s$，转向$3$
      • 若${\lambda }_s\le {\epsilon }_{{\lambda }_s}$，则以$x_{k+1}+\sigma$代替$x_k$，$k$增加$1$，转向$2$
6. 已迭代$N$次，仍未达到要求，下山失败，此时需重新另选初值$x_0$，再转向$2$后开始新的计算

弦割法

$$\begin{array}{rl}& x_{k+1}=x_k-\frac{f(x_k)}{f^{\prime }(x_k)}\\ & f^{\prime }(x)\approx \frac{f(x_k)-f(x_{k-1})}{x_k-x_{k-1}}\\ \Rightarrow & x_{k+1}=x_k-\frac{f(x_k)}{f(x_k)-f(x_{k-1})}(x_k-x_{k-1})\end{array}$$

解线性代数方程组的直接法

Gauss消元法

Gauss顺序消元法

KaTeX parse error: Undefined control sequence: \mbox at position 30: …ongrightarrow_{\̲m̲b̲o̲x̲{初等行变换，依次消去各列的主…

• $\text{Gauss}$顺序消元法能进行下去$⟺$ 前$n-1$阶顺序主子式不为零

Gauss主元消元法

• 列主元消元法
• 全主元消元法

矩阵三角分解法

直接三角分解法

$$A=LU\{\begin{array}{ll}\text{L is Unit lower triangular matrix}& \text{Doolittle}\\ \text{U is Unit upper triangular matrix}& \text{Crout}\end{array}$$

前$n-1$阶顺序主子式不为零，$\text{Doolittle,Crout}$分解唯一
$$A=LU,Ax=b⟹(LU)x=b⟹\{\begin{array}{ll}Ly=b& \text{lower triangular matrix}\\ Ux=y& \text{upper triangular matrix}\end{array}$$

$$\begin{array}{rl}& \{\begin{array}{l}{y}_1=b_1,\\ y_k=b_k-{\sum }_{m=1}^{k-1}{l}_{km}{y}_m(k=2,3,\cdots ,n).\end{array}\\ & \{\begin{array}{l}{x}_n=\frac{y_n}{u_{nn}}\\ x_k=\frac{1}{u_{kk}}\left(y_k-{\sum }_{m=k+1}^n{u}_{km}{x}_m\right)(k=n-1,n-2,\cdots ,1)\end{array}\end{array}$$

$$\{\begin{array}{l}{u}_{kj}=a_{kj}-{\sum }_{m=1}^{k-1}{l}_{km}{u}_{mj}\Rightarrow a_{kj}(j=k,k+1,\cdots ,n),\\ l_{ik}=\frac{1}{u_{kk}}\left(a_{ik}-{\sum }_{m=1}^{k-1}{l}_{im}{u}_{mk}\right)\Rightarrow a_{ik}(i=k+1,k+2,\cdots ,n)\end{array}$$

列主元三角分解法

• 选第一列主元后，分解第一行第一列，以及$y_1$
• 求解第二列分解后的元素的临时大小，求最大行交换
• 求第2行第2列，$y_2$，重复

平方根法

$A$是$n$阶对称正定矩阵$\Longrightarrow $
$$A=L{L}^T\text{\hspace{1em}(Cholesky)}$$
$L$是非奇异下三角矩阵（如果规定对角元为正，分解唯一）

方程组的性态与误差分析

向量和矩阵范数见矩阵论
$$Cond(A)=\parallel A^{-1}\parallel \parallel A\parallel$$
降低方程组系数矩阵条件数的方法
$$Ax=b\to {\Lambda }_{1}A{\Lambda }_{2}y={\Lambda }_{1}b,x={\Lambda }_{2}y$$
使得
$$Cond({\Lambda }_{1}A{\Lambda }_2)=min$$
${\Lambda }_1,{\Lambda }_2$简单的设置方法有
$${\Lambda }_2=I,{\Lambda }_1=diag({\frac{1}{s}}_1,{\frac{1}{s}}_2,\cdots ,{\frac{1}{s}}_n),s_i=\underset{1\le j\le n}{\max }\mid a_{ij}\mid$$