误差和有效数字 基本概念

避免危害

数值分析: 将复杂的计算式转换为指令集中定义过的+、- (×、÷)输入到计算机中, 得到近似解

1.2 误差的背景介绍

1.2.1 来源与分类

来源	分类
从实际问题中抽象出数学模型	模型误差
通过测量得到模型中参数的值	观测误差
求近似解	方法误差(截断误差)
机器字长有限	舍入误差

1.2.2 误差与有效数字

绝对误差

绝对误差: $e^* = x^* - x, 其中x为精确值, x^*为x的近似值$

• 精确值$x^*$通常不可求得, 我们只能求出近似值

绝对误差限: $|e^*|的上限即为\varepsilon^*, 称为\textcolor{#66ccff}{绝对误差限}$

• 工程上常记为$x = x^* \pm e^*$, 如$\int_0^1e^{-x^2}dx=0.743\pm0.006$

• 注意:

  • e* 理论上讲是唯一确定的，可能取正，也可能取负。
  • e*> 0 不唯一，当然 e* 越小越具有参考价值。

相对误差

相对误差: $e_r^* = \frac{e^*}x ≈ \frac{e^*}{x^*}\\$

相对误差限: $\varepsilon_r^* (= \frac{\varepsilon^*}{x}) = \frac{\varepsilon^*}{|x^*|}\\$

有效数字

看不懂的有效数字定义: 用科学记数法, 记$近似值x^*=a_1.a_2a_3...a_n*10^m(其中a_1\neq 0)$, 若$(|e^*|=)\ |x^*-x|\leq \textcolor{red}{0.5}*10^{m-n+1}$(即$a_n$的截取按照四舍五入规则), 则称$x^*$为n位有效数字, 精确到$10^{m-n+1}$

• $m+1$: $x^*$小数点前的位数 $n$: 有效位数 $m-n+1$: 小数点后位数 (经过四舍五入的结果)
• 有效数字即为绝对误差限

有效数字与相对误差的关系

有效数字→相对误差限

• 已知$x^*$有n位有效数字, 其相对误差限为 $\begin{aligned}\varepsilon_r^*&=\left|\frac{\varepsilon^*}{x^*}\right|=\frac{0.5\times10^{-(n-1)}\times10^m}{a_1.a_2\cdots a_n\times10^m}=\frac{10^{-(n-1)}}{2\times a_1.a_2\cdots}\leq\frac1{2a_1}\times10^{-(n-1)}\end{aligned}$

相对误差限→有效数字

• 已知 x* 的相对误差限可写为$\varepsilon_r*=\frac1{2(a_1+1)}\times10^{-(n-1)}\\$

• 则: $\begin{aligned}\mid x-x^*\mid&\leq\varepsilon_r^*\cdot\mid x^*\mid=\frac{10^{-(n-1)}}{2(a_1+1)}\times a_1.a_2\cdots\times10^m\\&<\frac{10^{-(n-1)}}{2(a_1+1)}\cdot(a_1+1)\times10^m=0.5\times10^{m-n+1}\end{aligned}$

• $(a_1+1)$是为了通过不等式消除$0.a_2a_3...$

例: 为使$\pi^*$的相对误差小于0.001%,至少应取几位有效数字？

• 需求: 从有效数字确定相对误差限

• 设$\pi^*$取n位有效数字, 则其相对误差限为$\varepsilon_r^*\leq\frac1{2a_1}\times10^{-n+1}$

• 要保证相对误差小于0.001%, 只要保证$\varepsilon_r^*\leq\frac1{2a_1}\times10^{-n+1}<0.001\%$

• 已知 $\pi的a_1 = 3$，则从以上不等式可解得 $n > 6 - log6$，即$n \geq 6$，应取$\pi^* = 3.14159$

数值运算中的误差限

误差传递公式

误差: $\varepsilon = x^*-x$, 相对 误差$\varepsilon_r=\frac \varepsilon{x^*}\\$, 函数值相对误差$\frac{f(x^*)-f(x)}{f(x^*)}\\$

相对误差限比值(条件数): $C_p=|\frac{f(x^*)-f(x)}{f(x^*)}|/|\varepsilon_r|\approx|\frac{xf'(x^*)}{f(x^*)}| \\$​​

误差限公式

• 经过$f(x)$计算, 相对误差会被放大

加法: $y^*=x_1^*+x_2^*$

• $\varepsilon(y^*)\leq\varepsilon(x_1^*)+\varepsilon(x_2^*)$ 误差限直接相加

乘法: $y^*=x_1^*\cdot x_2^*$

• $\varepsilon(y^*)\leq|x_2^*|\varepsilon(x_1^*)+|x_1^*|\varepsilon(x_2^*)$
• 会受$|x_1^*|, |x_2^*|$影响

除法: $y^*=\frac{x_1^*}{x_2^*}\\$

• $\varepsilon(y^*)\leq\frac{|x_2^*|\varepsilon(x_1^*)+|x_1^*|\varepsilon(x_2^*)}{\left|x_2^*\right|^2}\\$

• $|x_2^*|\downarrow\ \Rightarrow\ \varepsilon(y^*)↑↑$

• 要尽量避免除法的出现, 防止误差暴涨

其他注意事项

1. 避免小分母: 分母过小会造成浮点溢出
2. 避免相近二数相减
   • 会导致有效数字减少
   • 几种经验性避免方法：
     • $\sqrt{x+\varepsilon}-\sqrt{x}=\frac\varepsilon{\sqrt{x+\varepsilon}+\sqrt{x}};\quad\ln(x+\varepsilon)-\ln x=\ln(1+\frac\varepsilon x);\\$
     • 当 | x | << 1 时：$\begin{aligned}&1-\cos x=2\sin^2\frac x2;\\&e^x-1=x\left(1+\frac12x+\frac16x^2+...\right)\end{aligned}$
3. 避免大数吃小数
   • 计算机浮点运算, 指数对齐导致的基数部分丢失
4. 尽量减少运算次数
   • 使用秦九昭算法

第2章 插值法

插值多项式的唯一性

拉格朗日插值公式, 误差, 余项

差商, 均差, 均差的性质

牛顿插值多项式的形式, 误差

• 2.3.4 不用

埃尔米特插值

• 例题, 三个点值+一个点的导数的计算
• 均差表绘制方式

分段低次插值: 知道概念即可

?

三次样条插值 不用

什么是插值法

当精确函数 y = f(x) 非常复杂或未知时，在一系列节点 $x_0 … x_n$处测得函数值$y_0 = f(x_0), … yn _= f(x_n)$，由此构造一个简单易算的近似函数 $P(x) ≈ f(x)$，满足条件$P(x_i) = f(x_i) (i = 0, … n)$​。这里的 P(x) 称为f(x) 的插值函数。

最常用的插值函数是多项式

插值多项式

为了使插值函数更方便在计算机上运算,一般插值函数都使用代数多项式和有理函数

代数插值多项式的存在唯一性

设函数$y=f(x)$在区间[a,b]上的代数插值多项式为$P(x)=a_0+a_1x+a_2x^2+...+a_nx^n, P_n(x_i)=y_i, i=0, 1,2,...,n$

即多项式$P_n(x)$的系数$a_0, a_1, ...,a_n$满足线性方程组

$\begin{cases}a_0+a_1x_0+a_2x_0^2+...+a_nx_0^n=y_0 \\ a_0+a_1x_1+a_2x_1^2+...+a_nx_1^n=y_1\\......\\ a_0+a_1x_n+a_2x_n^2+...+a_nx_n^n=y_n\end{cases}$

上述方程组的系数行列式为n+1阶的Vandermond行列式 ($x_0, x_1...$是已知量)

$V=\left | \begin{matrix} 1 & x_0 & ... & x_0^n \\ 1 & x_1 & ... & x_1^n \\ ... & ... & ... & ... \\ 1 & x_n & ... & x_n^n \\ \end{matrix} \right | =\prod_{i=0}^{n-1}\prod_{j=i+1}^n(x_j-x_i) \neq 0 (x_i\neq x_j)$

定理: 满足$P(x_i)=y_i, i=0,1,...,n$ , 次数不超过n的插值多项式一定唯一存在

• 若多项式次数≠n, 则插值多项式不唯一 例如$P(x)=L_n(x)+p(x)\prod_{i=0}^n(x-x_i)$也是一个插值多项式，其中p(x)可以是任意多项式。

2.2 拉格朗日插值

2.2.1 线性插值&抛物线插值

插值的目的是求出n次多项式$P(x)=a_0+a_1x+a_2x^2+...+a_nx^n\ 使得\ P_n(x_i)=y_i, i=0,1,...,n$

线性插值

即n=1, 已知$y_0=f(x_0), y_1=f(x_1)$ 求$P_1(x)=a_0+a_1x$满足$P_1(x_0)=y_0, P_2(x_1)=y_1$

由几何意义易得

$\begin{aligned} \boldsymbol{P_{1}}(\boldsymbol{x})& =y_0+\frac{y_1-y_0}{x_1-x_0}(x-x_0) \\ &={\left[\frac{x-x_1}{x_0-x_1}\right]}y_0+{\left[\frac{x-x_0}{x_1-x_0}\right]}y_1=\sum_{i=0}^{1}\textcolor{red}{l_i(x)}y_i \end{aligned}$ $\begin{aligned} l_0(x)=\frac{x-x_1}{x_0-x_1}\\l_1(x)=\frac{x-x_0}{x_1-x_0}\end{aligned}$

抛物线插值

即n=2, … 略!

2.2.2 拉格朗日插值多项式

定义—n次插值基函数

若n次多项式$l_j(x)(j=0,1,...,n)$在$n+1$个节点$x_0<x_1<...<x_n$上

满足条件$l_j(x_k)=\delta_{jk}=\begin{cases}1,&k=j\\0,&k\neq j\end{cases}\quad j,k=0,1,...,n$,

则称这n+1个n次多项式$l_0(x), l_1(x),...,l_n(x)$为插值节点$x_0,x_1,...,x_n$​上的n次插值基函数

与上面的例子类似, 可以推导得出n次插值基函数为$\begin{aligned} l_k(x)& =\frac{(x-x_0)(x-x_1)\cdots(x-x_{k-1})(x-x_{k+1})\cdots(x-x_n)}{(x_k-x_0)(x_k-x_1)\cdots(x_k-x_{k-1})(x_k-x_{k+1})\cdots(x_k-x_n)} \\ &=\prod_{i=0}^n\frac{\left(x-x_i\right)}{\left(x_k-x_i\right)}\quad\quad k=0,1,2,\cdots,n \\ \end{aligned}$​

• (分子只少了一项$x-x_k$, 分母对应的少一项$x_k-x_k$)

所以$\color{blue}L_n(x)=\sum_{k=0\\}^ny_kl_k(x)$, 由$l_k(x)$的定义可知$L_n(x_j)=\sum_{k=0\\}^n y_kl_k(x_j)=y_j,j=0,1,...,n$ 满足插值多项式的定义

将形如上述$L_n(x)$的插值多项式称为拉格朗日(Lagrange)插值多项式

拉格朗日插值多项式基函数的简化写法

记$\omega_{n+1}(x)=(x-x_0)(x-x_1)...(x-x_n)$, 可得$\omega_{n+1}^{\prime}(x_k)=(x_k-x_0)(x_k-x_1)\cdot \cdot\cdot (x_k-x_{k-1})(x_k-x_{k+1})\cdot \cdot \cdot (x_k-x_n)$ (导数乘法公式)

得$\begin{aligned} l_k(x)& \begin{aligned}&=\frac{(x-x_0)(x-x_1)\cdots(x-x_{k-1})(x-x_{k+1})\cdots(x-x_n)}{(x_k-x_0)(x_k-x_1)\cdots(x_k-x_{k-1})(x_k-x_{k+1})\cdots(x_k-x_n)}\end{aligned}=\textcolor{blue}{\frac{\omega_{n+1}(x)}{\omega_{n+1}^{\prime}(x_k)(x-x_k)}}\quad\quad\quad k=0,1,2,\cdots,n \end{aligned}$

公式总结:

$$\textcolor{blue}{L_n(x)}=\sum_{k=0\\}^ny_kl_k(x): Lagrange插值多项式\\ \textcolor{blue}{l_k(x)}=\frac{\omega_{n+1}(x)}{\omega_{n+1}^{\prime}(x_k)(x-x_k)}: Lagrange插值基函数$$

2.2.3 插值余项&误差估计

设区间[a, b]上f(x)的差值多项式为n阶的$L_n(x)$, 令余项$R_n(x) = f(x) - L_n(x)$

显然, 在插值节点$x_i$上有$R_n(x_i) = f(x_i) - L_n(x_i) = 0$

因此, $R_n(x)$在[a, b]上有n+1个零点

设$R_n(x) = K(x)\omega_{n+1}(x);\quad \omega_{n+1}(x) = (x-x_0)(x-x_1)...(x-x_n)$

此时有$R_n(x) = K(x)\omega_{n+1}(x) = f(x) - L_n(x)$

….

得

拉格朗日型余项定理:

$$R_n(x) = K(x)\omega_{n+1}(x) =\frac{f^{(n+1)}(\xi)}{(n+1)!}\omega_{n+1}(x), 其中\omega_{n+1}(x)=\prod_{i=0}^{n}(x-x_i), \xi\in(a,b)$$

2.3 均差与牛顿插值多项式

2.3.1 插值多项式的逐次生成

Lagrange插值多项式的基函数为

当需要增加节点时, 所有基函数$l_i(x)$都需要重新计算

考虑使用一种逐次生成插值多项式的方法, 记为$P_n(x)$

$对0次插值,P_0(x)=f(x_0)$

对1次插值, $\begin{cases}P_1(x_0)=f(x_0)\\P_1(x_1)=f(x_1)\end{cases}$

• $P_1(x)=f(x_0)+\frac{f(x_1)-f(x_0)}{x_1-x_0}(x_1-x_0)(点斜式)=P_0(x)+a_1(x-x_0)$

对2次插值, $\begin{cases}P_2(x_0)=f(x_0)\\P_2(x_1)=f(x_1)\\P_2(x_1)=f(x_2)\end{cases}$

• $P_2(x)=P_1(x)+a_2(x-x_0)(x-x_1)$

可得, $P_n(x)=a_0+a_1(x-x_0)+a_2(x-x_0)(x-x_1)+...+a_n(x-x_0)(x-x_1)...(x-x_n)$

其中$a_0, a_1, ...,a_n可由P_n(x_i)=f(x_i)$计算得到

$$\begin{aligned} &P(x_0)=f_0=a_0 &a_0=f_0 \\ &P(x_1)=f_1=a_0+a_1(x_1-x_0)&a_1=\frac{f_1-f_0}{x_1-x_0} \\ &P(x_2)=f_2=a_0+a_1(x_2-x_0)+a_2(x_2-x_0)(x_2-x_1)& a_{2}=\frac{\frac{f_2-f_0}{x_2-x_0}-\frac{f_1-f_0}{x_1-x_0}}{x_{2}-x_{1}} \end{aligned}$$

这样计算下去会变得很麻烦, 定义均差来表示他们

2.3.2 均差及其性质

差商(均差, devided difference)的定义:

$$\begin{aligned} &一阶差商: f[x_0,x_i]=\frac{f(x_1)-f(x_0)}{x_1-x_0}\\ &二阶差商: f[x_0, x_1,x_2] =\frac{f[x_1,x_2]-f[x_0,x_1]}{x_2-x_0}\\&...\\ &k+1阶差商: f[x_0,...,x_{k+1}]=\frac{f[x_0,x_1,...,x_k]-f[x_1,...,x_k,x_{k+1}]}{x_0-x_{k+1}} \\ &\qquad\qquad\qquad\qquad\qquad\qquad\ =\frac{f[x_0,...,x_{k-1},x_k]-f[x_0,...,x_{k-1},x_{k+1}]}{x_k-x_{k+1}} \end{aligned}$$

• 是递归定义的

• k+1阶差商只需要任选两个k阶差商相减, 并没有固定选择的要求

  • 所以在上面选择了$f[x_0,...,x_{k-1},x_{k+1}]$作为减数

均差的性质 ※

①k阶均差-f(x)的线性表示

$$\begin{aligned} f[x_0,x_1,\cdots,x_{k-1},x_k]&=\sum_{i=0}^k\frac{f(x_i)}{(x_i-x_0)\cdots(x_i-x_{i-1})(x_i-x_{i+1})\cdots(x_i-x_k)} \\ &=\sum_{i=0}^k\frac{f(x_i)}{\omega_{k+1}^{\prime}(x_i)} \\ &\textbf{其中}\quad\omega_{k+1}(x)=\prod_{i=0}^k\left(x-x_i\right),\quad\omega_{k+1}^{\prime}(x_i)=\prod_{j=0}^k\left(x_i-x_j\right) \end{aligned}$$

②差商与x的顺序无关

• 如 $f[x_0,x_1,x_2]=f[x_0,x_2,x_1]=f[x_2,x_1,x_0]$

③k阶差商与k阶导

当$f^{(k)}(x)$在包含节点$x_0,x_1,\cdots,x_k$的区间存在时，在$x_0,x_1,\cdots,x_k$之间必存在一点$\xi$,使得$f[x_0,x_1,\cdots,x_k]=\frac{f^{(k)}(\xi)}{k!}\\$

2.3.3 牛顿插值多项式

$N_n(x)=\alpha_0+\alpha_1(x-x_0)+\alpha_2(x-x_0)(x-x_1)+....+\alpha_n(x-x_0)...(x-x_{n-1})$

$\left.\left\{\begin{array}{cccc|c}f(x)=f(x_0)+f[x,x_0](x-x_0) & \\ f[x,x_0]=f[x_0,x_1]+f[x,x_0,x_1](x-x_1) \\ ....................\\ f[x,x_0,...,x_{n-1}]=f[x_0,...,x_{n}]+f[x,x_0,...,x_{n}](x-x_{n})\quad& \end{array}\right.\right.$

把后一式带入前一式得

$$\begin{aligned} f(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)+...\\&+f[x_0,...,x_n](x-x_0)...(x-x_{n-1}) &牛顿插值多项式N_n(x)\\ &+f[x,x_0,...,x_n](x-x_0)...(x-x_{n-1})(x-x_n)&牛顿插值多项式余项R_n(x) \end{aligned}$$

由上面的式子可以得出, n阶插值多项式$P_n(x)=f[x_0,...,x_n]x^n+?x^{n-1}+...+?x+?1$, 最高项系数一定为n阶差商

↑可以用来证明性质1(结合拉格朗日插值公式)

具体计算方式: 画出差商表

2.4 埃尔米特插值

埃米特插值: 不仅要求函数值相等, 还要求若干阶导数相等

埃尔米特插值: 要求插值函数 P (x) 满足$P(x_i) = f (x_i), P’ (x_i) = f ’ (x_i),…, P^{(m)} (x_i) = f^{ (m)} (x_i)$

注意:

• N个条件可以确定N-1阶多项式
• 要求在1个节点$x_0$处直到m阶导数都重合的插值多项式即为Taylor多项式
  • $\begin{gathered} \begin{aligned}P(x)=f(x_0)+f'(x_0)(x-x_0)+...+\frac{f^{(m)}(x_0)}{m!}(x-x_0)^m\end{aligned} \\ \textbf{其余项为}\quad R(x)=f(x)-P(x)=\frac{f^{(m+1)}(\xi)}{(m+1)!}(x-x_0)^{(m+1)} \end{gathered}$
• 一般只考虑$f$与$f ’$的值。

Hermite插值计算出的多项式次数

给定n+1个点前面的插值函数最高只有n阶 (列出n+1个方程, 解得n个待定系数)

两点三次Hermite插值

考虑只有两个节点的插值问题

设$f(x)$在节点$x_0,x_1$处的函数值为$y_0,y_1$ 在节点 $x_0,x_1$处的的一阶导数值为 $y_0^{\prime},y_1^{\prime}$

两个节点最高可以用$2\times1+1=3$次$Hermite$多项式$H_{3}(x)$

$H_3(x)$应满足插值条件

$$H_3(x_0)=y_0\quad H_3(x_1)=y_1\\ H_3^{\prime}(x_0)=y_0^{\prime}\quad H_3'(x_1)=y_1'$$

用四个基函数表示$\color{blue}H_3(x) = y_0α_0(x) + y_1\alpha_1(x)+y_0'\beta_0(x)+y'\beta_1(x)$

可得

$$\begin{aligned}&\color{red}{{\alpha_0(x_0)=1}}\quad\color{red}{{\alpha_0(x_1)=0}}\quad\color{red}{{\alpha_0^{\prime}(x_0)=0}}\quad\color{red}{{\alpha_0^{\prime}(x_1)=0}}\\ &\alpha_1(x_0)=0\quad\alpha_1(x_1)=1\quad\alpha_1^{\prime}(x_0)=0\quad\alpha_1^{\prime}(x_1)=0\\ &\beta_0(x_0)=0\quad\beta_0(x_1)=0\quad\beta_0^{\prime}(x_0)=1\quad\beta_0^{\prime}(x_1)=0\\ &\beta_0(x_0)=0 \quad\beta_1(x_0)=0\quad\beta_1^{\prime}(x_1)=0\quad\beta_1^{\prime}(x_1)=1 \ \end{aligned} \Longrightarrow \begin{cases}对\alpha(x): 只有\textcolor{blue}{{a_i(x_i)=1}} ,其余为0\\ 对\beta(x): 只有\textcolor{blue}{\beta_i'(x_i)=1},其余为0\end{cases}$$

因为$H_3(x)$是3次多项式, 所以$\alpha_0(x)$最高也为3次, $\color{orange}\begin{cases}\alpha_0(x_0)=1&\alpha_0'(x_0)=0&①\\\alpha_0(x_1)=0 &\alpha_0'(x_1)=0&②\end{cases}$

• 在$x_1$处, $\color{orange}②\quad \alpha_0(x_1)=\alpha_0'(x_1)=0$, 即$x_1$为$\alpha_0$的二重零点, 函数值为0, 导数也为0 可设$\color{blue}\alpha_0(x)=(x-x_1)^2(\alpha x+b)$
• 在$x_0$处, $\color{orange}①\quad \alpha_0(x_0)=1\ \ \alpha_0'(x_0)=0$ , 代入得$a=-\frac2{\left(x_0-x_1\right)^3}\quad b=\frac1{\left(x_0-x_1\right)^2}+\frac{2x_0}{\left(x_0-x_1\right)^3}$

代入得

$$\begin{aligned} \alpha_0(x)& =\textcolor{blue}{(x-x_1)^2(ax+b)} \\ &=(x-x_1)^2\left(-\frac{2x}{\left(x_0-x_1\right)^3}+\frac1{\left(x_0-x_1\right)^2}+\frac{2x_0}{\left(x_0-x_1\right)^3}\right) \\ &=\frac{(x-x_1)^2}{\left(x_0-x_1\right)^2}\quad\left(1+\frac{2x_0}{x_0-x_1}-\frac{2x}{x_0-x_1}\right) \\ &=\left(1+2\frac{x-x_0}{x_1-x_0}\right)\left(\frac{x-x_1}{x_0-x_1}\right)^2=\color{blue} (1+2l_1(x))\cdot l_0^2(x)\\ \text{即}&\quad\alpha_0(x)=(1+2l_1(x))\cdot l_0^2(x)~=\left(1+2\frac{x-x_0}{x_1-x_0}\right)\!\left(\frac{x-x_1}{x_0-x_1}\right)^2 \end{aligned}$$

类似可得

$$\begin{aligned}\\\alpha_1(x)&=(1+2l_0(x))\cdot l_1^2(x)=\left(1+2\frac{x-x_1}{x_0-x_1}\right)\left(\frac{x-x_0}{x_1-x_0}\right)^2\\\beta_0(x)&=(x-x_0)\cdot l_0^2(x)~=\left(x-x_0\right)\left(\frac{x-x_1}{x_0-x_1}\right)^2\\\beta_1(x)&=(x-x_1)\cdot l_1^2(x)~=\left(x-x_1\right)\left(\frac{x-x_0}{x_1-x_0}\right)^2\end{aligned}$$

代入$\color{blue}H_3(x) = y_0α_0(x) + y_1\alpha_1(x)+y_0'\beta_0(x)+y'\beta_1(x)$ 得

$$\begin{aligned} \textcolor{blue}{H_{3}}(& x)=y_{0}\alpha_{0}(x)+y_{1}\alpha_{1}(x)+y_{0}^{\prime}\beta_{0}(x)+y_{1}^{\prime}\beta_{1}(x) \\ &=y_0(1+2l_1(x))\cdot l_0^2(x)+y_1(1+2l_0(x))\cdot l_1^2(x) \\ &\qquad+y_0^{\prime}(x-x_0)\cdot l_0^2(x)^2+y_1^{\prime}(x-x_1)\cdot l_1^2(x) \\ &=y_0\Bigg(1+2\frac{x-x_0}{x_1-x_0}\Bigg)\Bigg(\frac{x-x_1}{x_0-x_1}\Bigg)^2+y_1\Bigg(1+2\frac{x-x_1}{x_0-x_1}\Bigg)\Bigg(\frac{x-x_0}{x_1-x_0}\Bigg)^2 \\ &\qquad+y_0^{\prime}(x-x_0)\left(\frac{x-x_1}{x_0-x_1}\right)^2+y_1^{\prime}(x-x_1){\left(\frac{x-x_0}{x_1-x_0}\right)}^2 \end{aligned}$$

两点三次Hermite插值的余项

两点三次Hermite插值的误差为$R_3(x) = f(x)-H_3(x)$

有$$\begin{cases} R_3(X_i) =f(x_i)-H_3(x_i)=0 \ R_3^{\prime}(X_i) =f^{{\prime}(x_i)-H}{\prime}_3(x_i)=0 \end{cases} \qquad i=0,1$$

所以$x_0, x_1$均为$R_3(x)$的二重零点, 可设$R_3(x) = K(x)(x-x_0)^2(x-x_1)^2, K(x)待定$

设辅助函数$\varphi(t)=f(t)-H_3(t)-K(x)(t-x_0)^2(t-x_1)^2$

则有

$$\begin{cases} \varphi(x_i) \begin{aligned}&=f(x_i)-H_3(x_i)-K(x)(x_i-x_0)^2(x_i-x_1)^2=0\end{aligned} &i=0,1\\ \varphi(x) =f(x)-H_3(x)-K(x)(x-x_0)^2(x-x_1)^2=0 \end{cases}$$

至少有5个零点 ($\varphi(x_0), \varphi(x_1)四个, \varphi(x)至少一个$)

用4次Rolle定理, 即可得, $至少存在一点\xi \in[x_0, x_1], 使得\varphi^{(4)}(x)=0$, 即 $\varphi^{(4)}(\xi) = f^{(4)}(\xi)-4!K(x) = 0 \quad(H_3(x)三次,导没了)$

所以,两点三次Hermite插值的余项为

$$R_3(x) = \frac{f^{(4)}(\xi)}{4!}(x-x_0)^2(x-x_1)^2 \quad其中\xi \in[x_0, x_1]$$

高次Hermite插值

作为多项式插值,三次已是较高的次数，次数再高就有可能发生Runge现象

因此，对有n+1节点的插值问题，我们可以使用分段两点三次Hermite插值

三点+一导数 插值

将导数的值算作两个点相近的一阶差商, 使用牛顿插值计算

例：求一个次数不高于3的多项式$P_3(x)$，使其满足$P_3(0)=0,P_3(1)=1,P'_3(1)=3,P_3(2)=1$​ 。

$$P_3(x)=f(0)+f[0,1](x-0)+f[0,1,1](x-0)(x-1)+f[0,1,1,2](x-0)(x-1)(x-1)\\ P_3(x)=0+(x-0)+2(x-0)(x-1)-\frac52(x-0)(x-1)(x-1)=-\frac52x^3+7x^2-\frac72x$$

这里相当于设置了四个点0, 1, 1, 2, 其中一阶差商f[0,1], f[1,2]正常计算

对f[1, 1]

• 观察差商定义$f[x_0,x_i]=\frac{f(x_1)-f(x_0)}{x_1-x_0}$和导数定义$f'(x_0)=lim_{h→0}\frac {f(x_0+h)−f(x_0)}{h}$
• 可以发现$f[x_0, x_0]=f'(x_0)$

2.5 分段低次插值

Runge现象

$\textbf{在}[-5,5]\textbf{上考察}f(x)=\frac1{1+x^2}\textbf{ 的}L_n(x)\text{。 取 }x_i=-5+\frac{10}ni\quad(i=0,...,n)$

• n越大，端点附近抖动越大，称为Runge现象 ($L_n(x) \not\rightarrow f (x)$)

分段线性插值

在每个区间$[x_i, x_{i+1}]$​上，用1阶多项式 (直线) 逼近 f (x):

$$f(x)\approx P_1(x)=\frac{x-x_{i+1}}{x_i-x_{i+1}}y_i+\frac{x-x_i}{x_{i+1}-x_i}y_{i+1} \quad (for\ each x\in [x_i, x_{x+1}])$$

• 当取区间($max|x_{i+1}-x_i|$​)小时, 逼近f(x)
• 但是失去了光滑性

分段Hermite插值

$给定x_0, . . . , x_n; y_0, . . . , y_n; y_0^{\prime }, . . . , y_n^{\prime }$ 在$[ x_i, x_{i+ 1}] { 上 利 用 两 点 的 }y{及 }y$’构造3次Hermite函数

• 导数一般不易得到。

第三章 不需要了

第4章 数值积分

概论

• 求积节点 求积系数概念
• 待定系数的确定
• 求积公式余项证明p102公式 梯形公式余项

Newton-Cotes公式

• 一阶二阶
• 稳定性 为什么不用这个公式
• 偶阶xxx

复合求积公式

• 复合梯形, 辛普森 及其余项
• p108 例题

龙贝格求积公式

• 龙贝格算法 (注意什么时候停止计算

• p113 例6

4.5

• 自适应skip
• 高斯: 知道概念
• 高斯 lelangde公式

4.1 数值积分概论

4.1.1 数值积分基本思想

对于积分 $I(f)=\int_a^bf(x)dx$

如果知道$f(x)$的原函数$F(x)$, 则由牛顿-莱布尼茨公式有

$$\int_a^bf(x)dx = F(x)|_a^b = F(b)-F(a)$$

但是在工程和科研中, 常出现以下问题

1. $f(x)$的解析式不存在, 只给出了$f(x)$的数值
2. $f(x)$的原函数$F(x)$求不出来, 如$F(x)$不是初等函数
3. $f(x)$表达式复杂, $F(x)$很难求出

简单的积分近似计算方式:

$$梯形公式: T=\frac{b-a}2[f(b)+f(a)]\\ 中矩形公式: R=(b-a)f(\frac{b+a}2)\\ 机械求积: \textcolor{blue}{\int_a^bf(x)dx\approx \sum_{k=0}^nA_kf(x_k)}, 其中x_k称为求积节点, A_k称为求积系数$$

• 求积节点一般是给定的, 我们的目标就是确定求积系数$A_k$

4.1.2 插值多项式计算积分

积分的近似计算方法很多,但为方便起见,最常用的一种方法是利用插值多项式来构造数值求积公式,具体步骤如下:

在积分区间$[a, b]$上取n+1个插值节点$a\leq x_0< x_1...< x_n\leq b$

$f(x)$的n次插值多项式: $L_n(x) = \sum_{k=0}^nf(x_k)l_k(x), l_k(x)为插值基函数$

用$L_n(x)$作为$f(x)$的近似, 此时积分的计算为

$$\begin{gathered} \int_a^bf(x)dx \approx\int_a^bL_n(x)dx=\int_a^b\sum_{k=0}^nf(x_k)l_k(x)dx \\ =\sum_{k=0}^nf(x_k)\int_a^bl_k(x)dx \end{gathered}$$

设求积系数$A_k=\int_a^bl_k(x)dx$则$f(x)$的积分$I(f)=\int_a^bf(x)dx\approx\sum_{k=0}^nA_kf(x_k)=I_n(f)$

4.1.3 代数精度

定义 若求积公式$\int_a^bf(x)dx\approx\sum_{k=0}^nA_kf(x_k)$

• 对任意次数不超过$m$次 的代数多项式 $P_i( x) ( i\leq m)$都准确成立，即$\int_a^bP_i(x)dx=\sum_{k=0}^nA_kP_i(x_k)\quad i=0,1,\cdots,m$
• 但对$m+1$次多项式却不能准确成立，即只要$\int_a^bx^{m+1}dx\neq\sum_{k=0}^nA_kx_k^{m+1}$​​
• 则称该求积公式具有m次的代数精度

例: 求梯形公式的代数精度

例：对于$[a,b]$上1次插值，有 $L_1(x)=\frac{x-b}{a-b}f(a)+\frac{x-a}{b-a}f(b)$$\Longrightarrow A_1=A_2=\frac{b-a}2\quad\Longrightarrow\int_a^bf(x)dx\approx\frac{b-a}2[f(a)+f(b)]$ 考察其代数精度。

解：逐次检查公式是否精确成立

代入0次的代数多项式$P_0=1:\int_a^b1dx=b-a=\frac{b-a}2[1+1]$

代入1次$P_1=x:\int_a^bxdx=\frac{b^2-a^2}2=\frac{b-a}2[a+b]$

代入2次$P_{2}= x^{2}: \int _{a}^{b}x^{2}dx= \frac {b^{3}- a^{3}}3\neq \frac {b- a}2[ a^{2}+ b^{2}]$

4.2 Newton-Cotes数值求积分

4.2.1 NewTon-Cotes公式

牛顿–柯特斯公式

Newton-Cotes公式是指等距节点下使用Lagrange插值多项式建立的数值求积公式

设函数$f(x)\in C[a,b]$, ($C[a,b]: 在[a,b]上连续的函数集合$)

将积分区间$[a,b]$分割成n等分, 各个节点为$x_k = a + kh, h = \frac{b-a}{n}$

使用Lagrange插值得: $L_n(x)=\sum_{k=0}^nf(x_k)l_k(x)\quad R_n(x)=\frac{f^{(n+1)}(\xi)}{(n+1)!}\omega_{n+1}(x)\\$

• 其中 $l_k(x)=\prod_{\substack{0\leq j\leq n\\j\neq k}}\frac{x-x_j}{x_k-x_j}\quad\xi\in[a,b]\quad\omega_{n+1}(x)=\prod_{i=0}^n(x-x_i)$

此时积分准确值$I = \int_a^bf(x)dx = \int_a^b[L(x)+R(x)]dx = \int_a^bf(x_k)l(x_k)dx + \int_a^bR(x)dx = \sum_{k=0}^n A_kf(x_k)+\int_a^bR(x)dx\\$

• 其中 $A_k=\int_a^bl_k(x)dx=\int_a^b\prod_{0\leq j\leq n \and j\neq k}\frac{x-x_j}{x_k-x_j}dx\\$

令$I_n(f)=\sum_{k=0}^nA_kf(x_k)\\R(I_n)=\int_a^bR_n(x)dx$得, $I(f) = I_n(f)+R(I_n), I(f)\approx I_n(f)$

• n阶Newton-Cotes求积公式: $I_n(f)=\sum_{k=0}^nA_kf(x_k)$
• Newton-Cotes公式的余项: $R(I_n)=\int_a^bR_n(x)dx \quad 其中, R_n(x)=\frac{f^{(n+1)}(\xi)}{(n+1)!}\omega_{n+1}(x)\\$

A_k的计算与Cotes系数

$A_k=\int_a^bl_k(x)dx=\int_a^b\prod_{0\leq j\leq n \and j\neq k}\frac{x-x_j}{x_k-x_j}dx\\$

$$\begin{aligned} &设x=a+th, 由x\in[a, b], 得t\in[0, n]\\ A_{k} &=\int_a^b\prod_{\substack{0\leq j\leq n\\j\neq k}}\frac{x-x_j}{x_k-x_j}dx=\int_0^n\Bigg(\prod_{\substack{0\leq j\leq n\\j\neq k}}\frac{(t-j)h}{(k-j)h}\Bigg)\cdot h\cdot dt\\ & =\frac{h\cdot(-1)^{n-k}}{k!\cdot(n-k)!}\int_0^n\prod_{\begin{matrix}0\leq j\leq n\\j\neq k\end{matrix}}(t-j)dt\\ & =(b-a)\cdot\frac{(-1)^{n-k}}{n\cdot k!\cdot(n-k)!}\int_0^n\prod_{0\leq j\leq n}(t-j)dt\\ \end{aligned}$$

$A_k=\begin{pmatrix}b-a\end{pmatrix}\cdot C_k^{(n)}$

$\therefore I_n(f)=\sum_{k=0}^nA_kf(x_k)=\begin{pmatrix}b-a\end{pmatrix}\cdot \sum_{k=0}^nC_k^{(n)}f(x_k)\\$

• $C_k^{(n)}$称为Cotes系数
• Cotes 系数仅取决于 n和k，可查表得到。与 f (x) 及区间[a, b]均无关。

4.2.2 低阶Newton-Cotes公式及其余项

在Newton-Cotes公式中,n=1,2,4时的公式是最常用也最重要三个公式,称为低阶公式

梯形公式

梯形求积公式

$$\begin{aligned} &n=1,x_0=a,x_1=b,h=b-a\\ &\text{Cotes系数为} \begin{cases}C_0^{(1)}=-\int_0^1(t-1)dt=\frac12\\ C_1^{(1)}=\int_0^1tdt=\frac12\end{cases} \\&I_1(f)=(b-a)\sum_{k=0}^1C_k^{(1)}f(x_k)=\frac{b-a}2[f(x_0)+f(x_1)] \\ &\text{即 }I_1(f)=\frac{b-a}2[f(a)+f(b)]\end{aligned}$$

梯形求积公式(两点公式): $T=I_1(f)=\frac{(b-a)}2[f(a)+f(b)]\\$

梯形公式余项

$R_n(x)=\frac{f^{(n+1)}(\xi)}{(n+1)!}\omega_{n+1}(x)$

$R(T)=R(I_1)=\int_a^bR_1(x)dx$

$\begin{aligned} R(T)& =\int_a^b\frac{f''(\xi_x)}{2}(x-a)(x-b)dx \\ &=\frac{f^{\prime\prime}(\eta)}2\int_a^b(x-a)(x-b)dx \\ &=-\frac{f^{\prime\prime}(\eta)}2\frac{(b-a)^3}6 \\ &=-\frac{(b-a)^3}{12}f^{\prime\prime}(\eta) \end{aligned}$

$\therefore\ \mid R(T)\mid\leq\frac{(b-a)^3}{12}M_2\quad M_2=\max_{x\in[a,b]}\mid f^{\prime\prime}(x)\mid$

• 梯形(trapezia)公式具有1次代数精度

Simpson公式

$$\begin{aligned} &\text{取}n=2,\text{则}x_0=a,x_1=\frac{b+a}2,x_2=b,h=\frac{b-a}2& \\ &\text{Cotes系数为}\quad \begin{cases}C_0^{(2)}=\frac{1}{4}\int_0^2(t-1)(t-2)dt=\frac{1}{6}\\ C_1^{(2)}=\frac{-1}{2}\int_0^2t(t-2)dt=\frac{4}{6} \\ C_2^{(2)}=\frac{1}{4}\int_0^2(t-1)tdt=\frac{1}{6} \end{cases}\\ &\text{求积公式为} \quad I_2=(b-a)\sum_{k=0}^2C_k^{(2)}f(x_k)\\ \end{aligned}$$

Simpson求积公式(三点公式, 抛物线公式): $\begin{aligned}S=I_2(f)&=[b-a](\frac16f(x_0)+\frac46f(x_1)+\frac16f(x_2))\\&=\frac{b-a}6[f(a)+4f(\frac{a+b}2)+f(b)]\end{aligned}$

Simpson公式的余项: $R(S)=R(I_2)=\int_a^bR_2(x)dx =-\frac{b-a}{180}\Big(\frac{b-a}{2}\Big)^4f^{(4)}(\eta)\\$

• Simpson公式具有3次代数精度

	公式	余项$R[f]$	步长h	余项$R[f]$ (用h表示)	代数精度
n=1(梯形公式)	$\frac{(b-a)}2[f(a)+f(b)]\\$	$-\frac{(b-a)^3}{12}f^{\prime\prime}(\xi)\\$	$b-a$	$-\frac{(b-a)^3}{12}f^{\prime\prime}(\xi), \xi\in[a,b]\\$	1
n=2(Simpson公式)	$\frac{b-a}6[f(a)+4f(\frac{a+b}2)+f(b)]\\$	$-\frac{b-a}{180}(\frac{b-a}2)^4f^{(4)}(\xi)\\$	$\frac{b-a}2\\$	$-\frac1{90}h^5f^{(4)}(\xi), \xi\in[a,b]\\$	3
n = 3: Simpson’s 3/8-Rule			$\frac{b-a}3\\$	$-\frac3{80}h^5f^{(5)}(\xi)\\$	3
n = 4: Cotes Rule		$-\frac{2(b-a)}{945}(\frac{b-a}4)^6f^{(6)}(\eta)\\$	$\frac{b-a}4\\$	$-\frac8{945}h^7f^{(6)}(\xi)\\$	5

定理：当n为偶数时，Newton-Cotes公式至少具有n+1次代数精度

证明：当n为偶数时， Newton-Cotes公式对$x^{n+1}$的余项为0

$$\begin{aligned} &\begin{cases} R_n(x)=f(x)-P_n(x)=\frac{f^{(n+1)}(\xi)}{(n+1)!}(x-x_0)\cdots(x-x_n)=\prod_{j=0}^n(x-x_j) \\\\ R(I_n)=\int_a^bR_n(x)dx=\int_a^b\prod_{j=0}^n(x-x_j)dx\\ \end{cases} \\ &令x=a+th\quad \int_a^bR_n(x)dx=\int_0^nh\prod_{j=0}^n(a+th-a-jh)dt=h^{n+2}\int_0^n\prod_{j=0}^n(t-j)dt \\ &\text{令}t=\mu+\frac n2 \quad\int_0^n\prod_{j=0}^n(t-j)dt=\int_{-\frac n2}^{\frac n2}\prod_{j=0}^n(\mu+\frac n2-j)d\mu=\int_{-\frac n2}^{\frac n2}\prod_{j=-\frac n2}^{\frac n2}(\mu-j)d\mu \\ &被积函数\prod_{j=-\frac n2}^{\frac n2}(\mu-j)是奇函数 \\ &令g(\mu)= \prod_{-\frac n2}^{\frac n2}(\mu-j)=(\mu+\frac n2)(\mu+\frac n2-1)\cdots\mu\cdots(\mu-\frac n2+1)(\mu-\frac n2) \\ &g(-\mu)= (-\mu+\frac{n}{2})(-\mu+\frac{n}{2}-1)\cdots(-\mu)\cdots(-\mu-\frac{n}{2}+1)(-\mu-\frac{n}{2}) \\ &g(-\mu)= \begin{pmatrix}-1\end{pmatrix}^{n+1}g(\mu)=-g(\mu) \\ &\therefore \int_{-\frac n2}^{\frac n2}\prod_{j=-\frac n2}^{\frac n2}(\mu-j)d\mu = 0\\ &\therefore R(I_n)=\int_a^bR_n(x)dx=\int_0^n\prod_{j=0}^n(t-j)dt = \int_{-\frac n2}^{\frac n2}\prod_{j=-\frac n2}^{\frac n2}(\mu-j)d\mu = 0 \end{aligned}$$

4.2.3 Newton-Cotes公式的稳定性

对Cotes系数: $C_k^{(n)}=\frac{(-1)^{n-k}}{n\cdot k!\cdot(n-k)!}\int_0^n\prod_{0\leq j\leq n}(t-j)dt\\$

值只和积分区间$[a,b]$的阶段$x_j$划分有关, 与函数无关

因此用Newton-Cotes公式计算积分的舍入误差主要由函数值$f(x_k)$的计算引起

考虑$f(x_k)$的舍入误差对公式的影响

记精确值为$f(x_k)$,近似值为$\tilde{f}(x_k)$, 误差$\varepsilon_k=f(x_k)-\tilde{f}(x_k)$

则积分精确值和近似值误差为: $I_n-\tilde{I}_n=(b-a)\sum_{k=0}^nC_k^{(n)}[f(x_k)-\tilde{f}(x_k)]\\$

$$\begin{aligned}&I_n-\bar{I}_n&&=(b-a)\sum_{k=0}^nC_k^{(n)}\varepsilon_k\\ &\left|I_n-\bar{I}_n\right|&&\leq(b-a)\sum_{k=0}^n\left|C_k^{(n)}\right| |\varepsilon_k|\\&&&\leq(b-a)\varepsilon\sum_{k=0}^n\left|C_k^{(n)}\right|\end{aligned}$$

• $\text{性质}:\sum_{k=0}^nC_k^{(n)}=1$
• $\varepsilon=\max\{|\varepsilon_k|\}$

$\text{若 }\forall k\leq n\text{,}C_k^{(n)}>0\text{,有}\left|I_n-\vec{I}_n\right|\leq(b-a)\varepsilon$

结论

Newton-Cotes公式的舍入误差只是函数值误差的**(b-a)倍**

即$\forall k\leq n\text{ ,C}_k^{(n)}>0$时, Newton-Cotes公式是稳定的

• 稳定: 指误差是否会在计算过程中显著增长
• 实际上, 当n<8时, 公式都是稳定的

若$C_k^{(n)}$有正有负，有$(b-a)\varepsilon\sum_{k=0}^n\left|C_k^{(n)}\right|\geq(b-a)\varepsilon\sum_{k=0}^nC_k^{(n)}=(b-a)\varepsilon$

此时,公式的稳定性将无法保证

• 因此,在实际应用中一般不使用高阶Newton-Cotes公式, 而是采用低阶复合求积法

4.3 复合求积公式

当积分区间$[a,b]$长度较大, 节点数n+1固定时, 使用Newton-Cotes的误差较大

• 考虑提高n, 过高次数的插值, 公式的舍入误差很难控制
• Newton-Cotes公式在$n\geq 8$时不稳定

为了提高公式的精度,又使算法简单易行,往往使用复合方法

• 将$[a,b]$分为若干小区间, 在小区间内使用低阶Newton-Cotes, 最后累加

高次插值有Runge 现象，故采用分段低次插值 → 分段低次合成的 Newton-Cotes 复合求积公式。

4.3.1 复合梯形公式

分割: $h = \frac{b-a}{n}, x_k=a+hk(k=0,1,\cdots,n)$

在第k个区间($[x_k,x_{k+1}]$)上用梯形公式, 然后累加:

$$\int_{x_k}^{x_{k+1}}f(x)dx\approx\frac{x_{k+1}-x_k}2[f(x_k)+f(x_{k+1})],\quad k=0,...,n-1\quad \Longrightarrow \\\int_a^bf(x)dx\approx\sum_{k=0}^{n-1}\frac h2[f(x_k)+f(x_{k+1})]=\frac h2{\left[f(a)+2\sum_{k=1}^{n-1}f(x_k)+f(b)\right]}=\color{blue}{T_n}$$

$$\textcolor{blue}{T_n}=\frac h2[f(a)+2\sum_{k=1}^{n-1}f(x_k)+f(b)]$$

• 头, 末尾($a = x_0, b = x_n$)加1次, 其他节点被加两次

4.3.2 复合Simpson公式

复合Simpson公式:$h=\frac{b-a}n,x_k=a+kh\quad(k=0,...,n) \\ \int_{x_k}^{x_{k+1}}f(x)dx\approx\frac h6[f(x_k)+4f(x_{k+\frac12})+f(x_{k+1})]$

• 因为simpson最少需要三个插值点, 在$x_k和x_{k+1}$之间插入一个点$x_{k+\frac12}$
  • n=2→n=3, 并不会增加过多的误差神奇的分数下标
• $x_k$被加两次→$2f(x_k)$, $x_{k+\frac{1}{2}}$不会重复加→$4f(x_{k+\frac12})$

得$\int_a^bf(x)dx\approx\frac h6[f(a)+4\sum_{k=0}^{n-1}f(x_{k+\frac12})+2\sum_{k=1}^{n-1}f(x_k)+f(b)]=\color{blue}{S_n}\\$

例: 使用各种复合求积公式计算定积分 $I=\int_0^1\frac{\sin x}xdx\\$

为简单起见,依次使用8阶复合梯形公式、4阶复合Simpson公式和2阶复合Cotes公式

可得各节点的值如表

Trapz	Simpson	Cotes	$x_i$	$f(x_i)$
$x_0$	$x_0$	$x_0$	0	1
$x_1$	$x_{0+\frac{1}{2}}$	$x_{0+\frac{1}{4}}$	0.125	0.99739787
$x_2$	$x_1$	$x_{0+\frac{1}{2}}$	0.25	0.98961584
$x_3$	$x_{1+\frac{1}{2}}$	$x_{0+\frac{3}{4}}$	0.375	0.97672674
$x_4$	$x_2$	$x_1$	0.5	0.95885108
$x_5$	$x_{2+\frac{1}{2}}$	$x_{1+\frac{1}{4}}$	0.625	0.93615564
$x_6$	$x_3$	$x_{1+\frac{1}{2}}$	0.75	0.90885168
$x_7$	$x_{3+\frac{1}{2}}$	$x_{1+\frac{3}{4}}$	0.875	0.87719257
$x_8$	$x_4$	$x_2$	1	0.841470981

$\begin{aligned} &T_8=\frac{1}{16}[f(0)+2\sum_{k=1}^7f(x_k)+f(1)]=0.94569086 \\ &S_4=\frac{1}{24}[f(0)+4\sum_{k=0}^3f(x_{k+\frac12})+2\sum_{k=1}^3f(x_k)+f(1)]=0.94608331 \\ &C_{2}=\frac1{180}[7f(0)+\sum_{k=0}^1[32f(x_{k+\frac14})+12f(x_{k+\frac24})+32f(x_{k+\frac34})]+14\sum_{k=1}^1f(x_k)+7f(1)]=0.94608307 \end{aligned}$

$\begin{aligned} &T_8 = 0.94569086 \quad 精度最低\\&S_4=0.94608331 \quad 精度次高\\&C_2=0.94608307\quad 精度最高\\ &精确值I =0.946083070367183\end{aligned}$

• 但是误差和收敛速度还需要考虑

4.3.3 复合求积公式的余项和收敛的阶

余项

	求积公式	单纯求积公式的余项	复合求积公式的余项
$T_n$	$\frac h2[f(a)+2\sum_{k=1}^{n-1}f(x_k)+f(b)]\\$	$-\frac{(b-a)}{12}(b-a)^2f''(\eta)\\$	$-\frac{(b-a)}{12}(h)^2f''(\eta)\\$
$S_n$	$\frac h6[f(a)+4\sum_{k=0}^{n-1}f(x_{k+\frac12})+2\sum_{k=1}^{n-1}f(x_k)+f(b)]\\$	$-\frac{b-a}{180}(\frac{b-a}2)^4f^{(4)}(\eta)\\$	$-\frac{b-a}{180}(\frac{h}2)^4f^{(4)}(\eta)\\$
$C_2$		$-\frac{2(b-a)}{945}(\frac{b-a}4)^6f^{(6)}(\eta)\\$	$-\frac{2(b-a)}{945}(\frac{h}4)^6f^{(6)}(\eta)\\$

• 复合求积公式的误差就是将小区间(长度为步长h)的误差逐个累加

收敛阶

定义: 若一个积分公式的误差满足 $\lim_{h\to0}\frac{R[f]}{h^p}=C<\infty, 且C\neq 0\\$ 则称该公式是$p$​ 阶收敛的。

• $T_n\thicksim O\left(h^2\right),S_n\thicksim O\left(h^4\right),C_n\thicksim O\left(h^6\right)$

例: 计算$π = \int_0^1\frac{4}{1+x^2}dx\\$

• 
• 

例2: 给定精度$\varepsilon$, 如何取n ( 要求$|I-T_n|< \varepsilon$, 如何判断n=?)

$|R[f]| = |-\frac{h^2}{12}(b-a)f''(\xi)|\leq |-\frac{h^2}{12}(b-a)|M_2$

上例中若要求$|I-T_n|<4\times10^{-6},\text{则}|R_n[f]|\leq\left|-\frac{h^2}{12}(b-a)\right|M_2=\frac{2h^2}3<4\times10^{-6}$

$\Longrightarrow h< 0. 00244949$ 即：取 $n=409$ 通常采取将区间不断对分的方法，即取$n=2^k$

上例中$2^k\geq 409\Rightarrow k= 9\textbf{ 时 , }T_{512}= 3. 14159202$

4.4 龙贝格求积公式

复合梯形公式: $\textcolor{blue}{T_n}=\frac {b-a}{2n}[f(a)+2\sum_{k=1}^{n-1}f(x_k)+f(b)]\\$

在复合梯形公式的基础上, 将[a, b]分隔成2n等份, 且$h = \frac{(b-a)}{n}$不变

即多出来下标为$k+\frac12, k=0,1,...,n-1$的点

$$\begin{aligned} T_{2n} &= \frac{(b-a)}{4n}[f(a)+2\sum_{k=1}^{n-1}f(x_k)+2\sum_{k=0}^{n-1}f(x_{k+\frac12})+f(b)]\\ &= \frac{h}{4}[f(a)+2\sum_{k=1}^{n-1}f(x_k)+f(b)]+\frac{h}{2}\sum_{k=0}^{n-1}f(x_{k+\frac12})\\ &= \textcolor{blue}{\frac{1}{2}T_n + \frac{h}{2}\sum_{k=0}^{n-1}f(x_{k+\frac12})}\\ &= \frac{1}{2}T_n + \frac{b-a}{2n}\sum_{k=0}^{n-1}f(a+(k+\frac12)h) \end{aligned}$$

• 在已经计算出$T_n$基础上, 如果精度不够, 就可以加n个点递归计算, 提高精度

新步长$h_k$

$$\begin{aligned} &记h_k = \frac{(b-a)}{2^k}, n=2^{k-1}\\ & 当节点数n=1时, h_0=b-a, n=2时, h_1=\frac{b-a}2\\ & T_1 = \frac{b-a}2[f(a)+f(b)] \quad(基本的梯形公式)\\ & T_2=\frac{1}{2}T_1+\frac{b-a}{2}f(a+\frac{1}{2}h)=T_0(1)=\frac{1}{2}T_0(0)+h_1f(a+h_1)\\ & 当n=2^k时, 记T_n=T_{2^k}=T_0(k) \quad {k=0,1,2,\dots} \end{aligned}$$

递推公式

$$\begin{aligned}T_0(0)&=\frac{b-a}2[f(a)+f(b)]\\T_0(k)&=\frac12T_0(k-1)+h_k\sum_{j=0}^{2^{k-1}-1}f(a+(2j+1)h_k)\quad k=1,2,\cdots\end{aligned}$$

上式称为递推的梯形公式

外推加速公式

复合梯形公式余项: $R_T=-\frac{(b-a)}{12}(h)^2f''(\eta)\\$

$$\begin{aligned} &由余项得\frac{R_{T2n}}{R_{Tn}} =\frac{I-T_{2n}}{I-T_n}\approx \frac14, 展开得I\approx \frac43T_{2n}-\frac13T_n\\\\ I&\approx \frac43T_{2n}-\frac13T_n\\ &=\frac46T_n+\frac23\sum f(x_{k+\frac12})-\frac13T_n\\ &=\frac13T_n + \frac23h\sum f(x_{k+\frac12})\\ &=\frac13 \frac{h}2[f(a)+2\sum f(x_k)+f(b)] + \frac23h\sum f(x_{k+\frac12})\\ &=\textcolor{red}{\frac h6[f(x_k)+4f(x_{k+\frac12})+f(x_{k+1})]=S_n} 即复合Simpson公式\\ \end{aligned}$$

• 可知, $S_n$的精度>

龙贝格公式(算法)

用上前面的新步长: $n=2^{k-1}, T_n=T_0(k-1), T_{2n}=T_0(k), I\approx \frac43T_0(k)-\frac13T_0(k-1)=S_{n}$

• 既然有$T_0(k)$, 那么肯定有$T_1(k)$

$T_1(k-1)=\frac43T_0(k)-\frac13T_0(k-1)=S_n=S_{2^{k-1}}\\$

同理, 对Simpson公式$S_n=S_{2^{k-1}}=T_1(k-1);\\S_{2n}=S_{2^{k}}=T_1(k)$

…

加速公式总结为: $T_m(k-1)=\frac1{4^m-1}[4^mT_{m-1}(k)-T_{m-1}(k-1)]\\$

• Romberg算法的代数精度为m的两倍
• Romberg算法的收敛阶高达m+1的两倍

4.5 高斯求积

只讲概念: 什么是高斯求积公式, 什么是高斯求积公式的的节点

构造具有2$n+1$次代数精度的求积公式

$$\int_a^b\rho(x)f(x)dx\approx\sum_{k=0}^nA_kf(x_k)$$

将节点 $x_0 … x_n$ 以及系数 $A_0 … A_n$ 都作为待定系数。令 $f (x) = 1, x, x^2, …, x^{2n+1}$​ 代入可求解，得到的公式具有2n+1 次代数精度。这样的节点称为Gauss点，公式称为Gauss 型求积公式。

高斯-勒让德求积公式

Legendre 多项式族 : 定 义 在 [ - 1, 1] 上 , $\rho ( x) \equiv 1$$P_k(x)=\frac1{2^kk!}\frac{d^k}{dx^k}(x^2-1)^k\\$ 满足：$(P_k,P_l)=\begin{cases}0&k\neq l\\\frac2{2k+1}&k=l&\end{cases}$

由$P_0= 1, P_1= x$有 递 推$( k+ 1) P_{k+ 1}= ( 2k+ 1) xP_k- kP_{k- 1}$

以$P_n+ 1$的根为节点的求积公式称为高斯-勒让德公式

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

5.1

• 5.1.2 skip

5.2

note: 5.2是分解法的推导过程, 5.3是公式(但是公式背不了一点)

• 高斯消去法
• 矩阵三角分解(LU分解) (推导过程)
• 列主元消去

5.3

• 矩阵三角分解的公式

  • 做题用算法/公式无所谓
  • p153/3.2 3.3 公式不一定要背 会做就行

• skip: 选主元三角分解法; 平方根分解法; 追赶法;

5.4

• 矩阵范数–各种算子范数计算
• 定理19 核半径

5.5

• 什么是病态良态, 计算
• 条件数, 计算

end

5.2 高斯消去法

5.2.1 消元与回代计算

对于线性方程组$\begin{cases}a_{11}x_1+a_{12}x_2+...+a_{1m}x_m = b_1\\...\\a_{n1}x_1+a_{n2}x_2+...+a_{nm}x_m = b_n \end{cases}$ 表示为Ax=b, A=\matrix{a_{11} &a_{12}&...&a_{1m}\\a_{21}&a_{22}&...&a_{2m}\\...\\a_{n1}&a_{n_2}&...&a_{nm}}

若$detA \neq 0$(行列式)对其增广矩阵施行行初等变换

$$\overline{A}=(A,b)\stackrel{\text{记}}{=}(A^{(1)},b^{(1)})=\begin{pmatrix}a_{11}^{(1)}&a_{12}^{(1)}&\cdots&a_{1n}^{(1)}&b_1^{(1)}\\a_{21}^{(1)}&a_{22}^{(1)}&\cdots&a_{2n}^{(1)}&b_2^{(1)}\\\vdots&\vdots&&\vdots&\vdots\\a_{n1}^{(1)}&a_{n2}^{(1)}&\cdots&a_{nn}^{(1)}&b_n^{(1)}\end{pmatrix}$$

经过n-1步消除列元素得

$$(A^{(1)},b^{(1)})\longrightarrow(A^{(n)},b^{(n)})=\begin{pmatrix}a_{11}^{(1)}&a_{12}^{(1)}&\cdots&a_{1n}^{(1)}&b_{1}^{(1)}\\&a_{22}^{(2)}&\cdots&a_{2n}^{(2)}&b_{2}^{(2)}\\&&\ddots&\vdots&\vdots\\&&&a_{nn}^{(n)}&b_{n}^{(n)}\end{pmatrix}$$

可知$a_{ii}^{( i) }\neq 0$ $i= 1, 2, \cdots , n$

因此，上三角形方程组$A^{(n)}x=b^{(n)}$有唯一解

因此可得线性方程组$Ax=b$的解：

$$\begin{cases}x_n=\frac{b_n}{a_{nn}^{(n)}}\\x_i=\frac{b_i^{(i)}-\sum_{j=i+1}^na_{ij}^{(i)}x_j}{a_{ii}^{(i)}},i=n-1,n-2,\cdots,2,1 \end{cases}$$

定理: 若A的所有顺序主子式均不为0，则高斯消元无需换行即可进行到底，得到唯一解。

增广矩阵通常用于判断矩阵的有解的情况，比如说秩（A)<秩（A|B) 方程无解； 秩（A)=秩（A|B) =n方程有唯一解； 秩（A)=秩（A|B) <n方程有无穷多解； 秩（A)>秩（A|B)不可能。

5.2.2 矩阵的三角分解

Doolittle分解法和Crout分解法

LU分解法 : 用矩阵描述高斯消去法的过程

矩阵A的LU分解法定义: 不带行交换的Gauss 消去法的消元过程,产生一个单位下三角矩阵L和一个上三角矩阵U,即

$$A = LU\\ \begin{gathered}L=\begin{pmatrix}1\\m_{21}&1\\m_{31}&m_{32}&1\\\vdots&\vdots&\vdots&\ddots\\m_{n-1,1}&m_{n-1,2}&m_{n-1,2}&\cdots&1\\m_{n,1}&m_{n,2}&m_{n,3}&\cdots&m_{n,n-1}&1\end{pmatrix}\\U=A^{(n)}=\begin{pmatrix}a_{11}^{(1)}&a_{12}^{(1)}&\cdots&a_{1n}^{(1)}\\&a_{22}^{(2)}&\cdots&a_{2n}^{(2)}\\&&\ddots&\vdots\\&&&a_{nn}^{(n)}\end{pmatrix}\end{gathered}$$

• 其中U的第一行$a_{1i}^{(1)}$​等于矩阵A的第一行
• L的第一列$m_{i1}=\frac{a_{i1}}{a^{(1)}_{11}}\\$
• 其他元素通过矩阵乘法 (解方程)计算 出来

定理: 若A的所有顺序主子式均不为0，则 A 的 LU 分解唯一（其中 L 为单位下三角阵）。

对于杜立特分解: 固定i: $j = i, i+1, ...,n$有$a_{ij}=\sum_{k=1}^{i-1}l_{ik}u_{kj}+u_{ij}, (l_{ij}=1)$

得 $u_{ij} =a_{ij}- \sum_{k=1}^{i-1}l_{ik}u_{kj}$

5.2.3 列主元消去法

在使用高斯消去法时, 每次确定行时, 加上交换行的操作

换行: 选择绝对值最大的列元素,换到上面

5.3 矩阵三角分解法

5.3.1 直接分解法

直接法是将原方程组化为一个或若干个三角形方程组的方法，共有若干种

$$\begin{aligned} A=&\begin{pmatrix}a_{11}&a_{12}&\cdots&a_{1n}\\a_{21}&a_{22}&\cdots&a_{2n}\\\vdots&\vdots&\vdots&\vdots\\a_{n1}&a_{n2}&\cdots&a_{nn}\end{pmatrix} & x=\begin{pmatrix}x_1\\x_2\\\vdots\\x_n\end{pmatrix}\quad & b=\begin{pmatrix}b_1\\b_2\\\vdots\\b_n\end{pmatrix}\\ &\quad\quad\quad\text{系数矩阵}&\ \ \text{未知量向量}\quad&\text{常数项} \end{aligned}$$

直接三角分解法

$A=LU \Rightarrow LUx=b \Rightarrow \begin{cases}Ly=b\\Ux=y\end{cases}$

杜立特分解法

$$A = LU\\ \begin{gathered}L=\begin{pmatrix}1\\m_{21}&1\\m_{31}&m_{32}&1\\\vdots&\vdots&\vdots&\ddots\\m_{n-1,1}&m_{n-1,2}&m_{n-1,2}&\cdots&1\\m_{n,1}&m_{n,2}&m_{n,3}&\cdots&m_{n,n-1}&1\end{pmatrix}\\U=A^{(n)}=\begin{pmatrix}a_{11}^{(1)}&a_{12}^{(1)}&\cdots&a_{1n}^{(1)}\\&a_{22}^{(2)}&\cdots&a_{2n}^{(2)}\\&&\ddots&\vdots\\&&&a_{nn}^{(n)}\end{pmatrix}\end{gathered}$$

• 其中U的第一行$a_{1i}^{(1)}$​等于矩阵A的第一行
• L的第一列$m_{i1}=\frac{a_{i1}}{a^{(1)}_{11}}\\$
• 其他元素通过矩阵乘法 (解方程)计算 出来

5.4 向量和矩阵范数

5.4.1 向量和矩阵的范数

向量范数

定义1 对于$n$维向量空间 $R^n$中任意一个向量 $x$, 若存在唯一一个实数 $\|x\|\in R$与$x$对应，且满足

• (正定性)$\|x\|\geq0$,且$\forall x\in R^n,\|x\|=0\Leftrightarrow x=0;$
• (齐次性)$\|\alpha x\| = | \alpha | \cdot \| x\|$, $\forall x\in R^n, \alpha \in R;$
• (三角不等式)$\|x+y\|\leq\|x\|+\|y\|,\forall x,y\in R^n.$
• 则称 $\|x\|$为向量$x$的范数

在向量空间$R^n(C^n)$中, $x = (x_1, x_2, ..., x_n)^T$, 常用的向量x的范数

• x的2-范数或欧氏范数: $\left\|x\right\|_2 =(\begin{array}{c}\left|x_1\right|^2+\left|x_2\right|^2+\cdots+\left|x_n\right|^2\end{array})^{1/2}$

• x的1-范数: $\left\|x\right\|_{1}=\left|x_1\right|+\left|x_2\right|+\cdots+\left|x_n\right|$

• x的∞范数(最大范数): $\left\|x\right\|_\infty =\max_{1\leq i\leq n}\left|x_i\right|\\$

• x的p-范数: $||x||_p = (|x_1|^p+|x_2|^p+...|x_n|^p)^{1/p}$

  • $\begin{gathered} \max_{1\leq i\leq n}\left|x_i\right| \leq\left(\left|x_1\right|^p+\left|x_2\right|^p+\cdots+\left|x_n\right|^p\right)^{1/p}\leq\left(n\max_{1\leq i\leq n}\left|x_i\right|^p\right)^{1/p} =n^{1/p}\max_{1\leq i\leq n}\left|x_i\right|\to\max_{1\leq i\leq n}\left|x_i\right|\left(p\to\infty\right) \end{gathered}$

定义 向量序列$\{\vec{x}^{(k)}\}$收敛于向量$\vec{x}^*$是指对每一个$1\leq i\leq n$, 都有$\lim_k\to\infty x_i^{(k)}=x_i^*$。

• 可以理解为$\|\vec{x}^{(k)}-\vec{x}^*\|_\infty\to0$

矩阵范数

定义: $R^{m\times m}$空间的矩阵范数$||\cdot||$ 对任意$A, B\in R^{m\times m}$ 满足以下条件

• 正定性: $\|A\|\geq0;\quad\|A\|=0\Leftrightarrow A=0$
• 齐次性: $||alpha A||=|\alpha|\cdot||A||$
• 三角不等式 : $||A+B||\leq ||A||+||B||$
• 相容 当 m = n 时, $|| AB || < || A || · || B ||$

算子范数:

$$\begin{aligned} &\textbf{由向量范数}\parallel\cdot\parallel_p\textbf{导出关于矩阵 }A\in R^{n\times n}\text{ 的 }p\text{ 范数:}\\&\parallel A\parallel_p=\max_{{\vec{x}}\neq0}\frac{\parallel A\vec{x}\parallel_p}{\parallel\vec{x}\parallel_p}=\max_{\parallel\tilde{x}\parallel_p=1}\parallel A\vec{x}\parallel_p\quad\text{则} \begin{cases} \parallel AB\parallel_p\leq\parallel A\parallel_p\parallel B\parallel_p\\ \parallel A\vec{x}\parallel_p\leq\parallel A\parallel_p\parallel\vec{x}\parallel_p \end{cases} \end{aligned}$$

特别有:

• 行和范数: $\parallel A\parallel_\infty=\max_{1\leq i\leq n}\sum_{j=1}^n\mid a_{ij}\mid\\$ (最大的 行之和)
• 列和范数: $\parallel A\parallel_1=\max_{1\leq j\leq n}\sum_{i=1}^n\mid a_{ij}\mid\\$ (最大的 列之和)
• 2-范数: $\parallel A\parallel_2=\sqrt{\lambda_{\max}(A^TA)}$ ($A^TA$矩阵的最大特征值)

谱半径

定义 矩阵A的谱半径记为$\rho(A)=\max_{1\leq i\leq n}|\lambda_i|$,其中$\lambda_i$​ 为A的特征根

5.5 线性方程组的误差分析

$||A||\cdot ||A^{-1}||$是关键的误差放大因子, 称为A的条件数, 即为$cond(A)$

• 越大, A越病态, 难以求得准确解

根据算子范数的不同, 条件数也不同

$$\begin{aligned} &cond(A)_1 =\left\|A\right\|_1\cdot\left\|A^{-1}\right\|_1 \\ &cond\left(A\right)_\infty =\left\|A\right\|_\infty\cdot\left\|A^{-1}\right\|_\infty \\ &cond\left(A\right)_2 =\left\|A\right\|_2\cdot\left\|A^{-1}\right\|_2=\sqrt{\lambda_{\max}\left(A^TA\right)}\sqrt{\frac{1}{\lambda_{\min}\left(A^TA\right)}} =\sqrt{\frac{\lambda_{\max}(A^TA)}{\lambda_{\min}(A^TA)}} \end{aligned}$$

条件数的性质

• A可逆, 则$cond(A)_p \geq 1$
• A可逆, $a\in R$,则$cond(\alpha A)=cond(A)$
• A正交, 则$cond(A)_2=1$
• A可逆, R正交, 则$cond(RA)_2=cond(AR)_2=cond(A)_2$

注：一般判断矩阵是否病态，并不计算$A^{-1}$，而由经验得出。

• 行列式很大或很小（如某些行、列近似相关）；
• 元素间相差大数量级，且无规则；
• 主元消去过程中出现小主元；
• 特征值相差大数量级。

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

6.1

• 什么是迭代 发散
• 序列极限
• 一些定理结论 不证明

6.2

• 雅可比 高斯-塞德尔迭代
• 定理9(1) ; (2), 定理10不用
• 超松弛skip

end

6.1 引言

在用直接法解线性方程组时要对系数矩阵不断变换

如果方程组的阶数很高，则运算量将会很大, 并且大量占用计算机资源

因此对线性方程组$Ax=b$, 要求找寻更经济、适用的数值解法

设$A\in R^{n\times n},b\in\mathbb{R}^n,x\in\mathbb{R}^n$

可以将线性方程组变换为$x=Bx+f$, 其中$B\in R^n\times n,f\in R^n,x\in R^n$

• B称为迭代矩阵, f为常数项
• 显然上面两式同解，我们称两个方程组等价

对第二个线性方程组，采用以下步骤：

取初始向量 $\vec x^{(0)}$,代入，可得$x^{(1)}=B\vec x^{(0)}+f$

依此类推

$$\begin{gathered} \vec x^{(2)}=B\vec x^{(1)}+f \\ \begin{array}{c}\bullet\\\bullet\end{array} \\ \vec x^{(k+1)}=B\vec x^{(k)}+f \\ \begin{pmatrix}k=0,1,2,\cdots\end{pmatrix} \end{gathered}$$

这种方式就称为迭代法 ,以上过程称为迭代过程

迭代法产生一个序列$\{\vec x^{(k)}\}_0^\infty$

如果其极限存在，即$\lim_{k\to\infty}\vec x^{(k)}=\vec x^*$, 则称迭代法收敛，否则称为发散

迭代法收敛的充分条件

设有线性方程组$x=Bx+f$ 以及一阶定常迭代法$x^{(k+1)}=Bx^{(k)}+f$

如果有迭代矩阵B的某种算子范数$||B||=q<1$, 则

• 迭代法收敛,即对任取$x^{(0)}$均有$\lim_{k\to\infty}\vec x^{(k)}=\vec x^*$, 且$x^{*}=B\vec x^{*}+f$

• $||x^*-x^{(k)}||\leq q^k||x^*-x^{(0)}||$

• $||x^*-x^{(k)}||\leq \frac{q}{1-q}||x^*-x^{(0)}||$

• $||x^*-x^{(k)}||\leq \frac{q^k}{1-q}||x^{(1)}-x^{(0)}||$

6.2 雅可比迭代法与高斯-塞德尔迭代法

雅可比迭代法

设线性方程组的一般形式为

$$\begin{cases}a_{11}x_1+a_{12}x_2+\cdots+a_{1n}x_n=b_1\\a_{21}x_1+a_{22}x_2+\cdots+a_{2n}x_n=b_2\\\cdots\cdots\cdots\cdots\\a_{n1}x_1+a_{n2}x_2+\cdots+a_{nn}x_n=b_n\end{cases}$$

设$a_{ii}\neq0\quad(i=1,2,\cdots,n),$则可从上式解出 $x_i$

$$如从第一行得: x_1=\frac1{a_{11}}[b_1-(a_{12}x_2+\cdots+a_{1n}x_n)]\\ x_2 = \frac1{a_{22}}[b_2-(a_{21}x_1+\cdots+a_{2n}x_n)]\\... \\ 此时原线性方程组化为: \\ \begin{cases}&x_1=\frac1{a_{11}}(b_1-\sum_{j=1}^na_{1j}x_j)=x_1+\frac1{a_{11}}(b_1-\sum_{j=1}^na_{1j}x_j)\\&x_2=\frac1{a_{22}}(b_2-\sum_{j=1}^na_{2j}x_j)=x_2+\frac1{a_{22}}(b_2-\sum_{j=1}^na_{2j}x_j)\\&\cdots\cdots\cdots\cdots\\&x_i=\frac1{a_{ii}}(b_i-\sum_{j=1}^na_{ij}x_j)=x_i+\frac1{a_{ii}}(b_i-\sum_{j=1}^na_{ij}x_j)\\&\cdots\cdots\cdots\cdots\\&x_n=\frac1{a_{nn}}(b_n-\sum_{j=1}^na_{nj}x_j)=x_n+\frac1{a_{nn}}(b_n-\sum_{j=1}^na_{nj}x_j)\end{cases}$$

令$D = diag(a_{11}, a_{22}, ...,a_{nn})$, 为方程组系数矩阵A的对角线

$$\begin{aligned}L=&\begin{pmatrix}0&0&\cdots&0\\-a_{21}&0&\cdots&0\\\vdots&\ddots&\ddots&\vdots\\-a_{n1}&-a_{n2}&\cdots&0\end{pmatrix}&A\text{的下三角部分的负矩阵} \\\\ U=&\begin{pmatrix}0&-a_{12}&\cdots&-a_{14}\\0&0&\cdots&-a_{24}\\\vdots&\ddots&\ddots&\vdots\\0&0&\cdots&0\end{pmatrix}&A\text{的上三角部分的负矩阵} \\\end{aligned}$$

$$A=D-L-U\\\begin{aligned} Ax=b&\Leftrightarrow(D-L-U)x = b \\ &\Leftrightarrow Dx=(L+U)x+b\\&\Leftrightarrow{x}=\underbrace{D^{-1}(L+U)}_{\boldsymbol{B_J}}{x}+\underbrace{D^{-1}\vec{b}}_{\vec{\boldsymbol{f}}}\end{aligned}$$

故迭代过程化为$x^{(k+1)}=D^{-1}(L+U)x^{(k)}+D^{-1}b$,

令$B_J=D^{-1}(L+U),f=D^{-1}b$,于是$x^{(k+1)}=B_Jx^{(k)}+f\ (k=0,1,2,\cdots)$

等价线性方程组为$x=B_Jx+f\longleftrightarrow Ax=b$, 即

$$\vec{x}^{(k+1)}=B_J\vec{x}^{(k)}+D^{-1}\vec{b}, B_J = D^{-1}(L+U)$$

• 称上式为解线性方程组的Jacobi迭代法(J法), 矩阵$B_J$为迭代法的迭代矩阵

Gauss-Seidel迭代法

雅可比迭代法中描述的$\vec x^{(k)}, \vec x^{(k+1)}$都是向量, 每轮迭代时还要从$\vec x^{(k)}=(x^{(k)}_1, x^{(k)}_2, ..., x^{(k)}_n)^T$的$x^{(k)}_i$中逐个计算$x^{(k+1)}_i$

• 我们知道$x^{(k+1)}_i$的误差会比$x^{(k)}_i$更小
• 在每轮迭代时, 对$x^{(k+1)}_i$, 它前面的$x^{(k+1)}_j(j=1, 2, ...,i-1)$都计算出来了, 可以将它用上

$$\begin{aligned} &x_1^{(k+1)} =\frac1{a_{11}}(-a_{12}x_2^{(k)}-a_{13}x_3^{(k)}-a_{14}x_4^{(k)}-\cdots-a_{1n}x_n^{(k)}+b_1) \\ &x_2^{(k+1)}= \frac1{a_{22}}(\textcolor{red}{-a_{21}x_1^{(k+1)}}-a_{23}x_3^{(k)}-a_{24}x_4^{(k)}-\cdots-a_{2n}x_n^{(k)}+b_2) \\ &x_3^{(k+1)}=\frac1{a_{33}}(\textcolor{red}{-a_{31}x_1^{(k+1)}-a_{32}x_2^{(k+1)}}-a_{34}x_4^{(k)}-\cdots-a_{3n}x_n^{(k)}+b_3) \\ &\cdots\\ &x_n^{(k+1)} =\frac1{a_{nn}}(\textcolor{red}{-a_{n1}x_1^{(k+1)}-a_{n2}x_2^{(k+1)}-a_{n3}x_3^{(k+1)}-\cdots-a_{nn-1}x_{n-1}^{(k+1)}}+b_n) \end{aligned}$$

• 红色部分用了$x^{(k+1)}_i$来更精确的进行计算

转换为矩阵形式

$$\begin{aligned} &\vec{x}^{(k+1)}=D^{-1}(L\vec{x}^{(k+1)}+U\vec{x}^{(k)})+D^{-1}\vec{b} \\ &\Leftrightarrow\quad(D-L)\vec{x}^{(k+1)}=U\vec{x}^{(k)}+\vec{b} \\ &\Leftrightarrow\vec{x}^{(k+1)}= \underbrace{{(D-L)^{-1}U}}_{B}\ \vec{x}^{(k)} +(\underbrace{D-L)^{-1}}_{\boldsymbol{\overline{f}}}\vec{b} \end{aligned}$$

• B: Gauss-Seidel 迭代阵

Jacobi迭代法和Gauss-Seidel迭代法统称为简单迭代法

---

例题:

用Jacobi迭代法求解方程组,误差不超过1e-4

取初值$x^{(0)} = [ 0 0 0] ^T$,使用Jacobi迭代法 $x^{(k+1)}=B_Jx^{(k)}+f\quad(k=0,1,2,\cdots n,\cdots)$

用Guass-Seidel迭代法, 选择相同初值$x^{(0)} = [ 0 0 0] ^T$

通过迭代,至第7步得到满足精度的解x7

• 可以看出,Gauss-Seidel迭代法的收敛速度比Jacobi迭代法要高

迭代法的收敛性

设解线性方程组的迭代格式: $\vec{x}^{(k+1)}=B\vec{x}^{(k)}+f$ (第k次迭代)

而方程组精确解$x^*$, 有$x^* = Bx^*+f$

两式相减得: $x^{(k+1)}-x^*=Bx^{(k)}-Bx^*=B\left(x^{(k)}-x^*\right)$

令$\varepsilon = x^{(k)}-x^*, k=0,1,\cdots$, 则$\varepsilon^{(k+1)}=B\varepsilon^{(k)}=B^2\varepsilon^{(k-1)}=\cdots=B^{k+1}\varepsilon^{(0)}$

• $\varepsilon^{(0)}=x^{(0)}-x^*$, 是一个非零常数向量, ($x^{(0)}$是选取的初值)
  • 当然, 如果$\varepsilon^{(0)}=0$, 说明恰好选到了精确值, 是特殊情况

因此迭代法收敛的充要条件为

$$\lim_{k\to\infty}\varepsilon^{(k+1)}=\lim_{k\to\infty}\left(x^{(k+1)}-x^{\star}\right)=0\\\text{可转变为}\lim_{k\to\infty}B^{k+1}=0$$

定理1: 迭代格式收敛的充要条件为 $\lim_{k\to\infty}B^{k}=0\\$

定理2: 迭代格式收敛的充要条件为谱半径$\rho(B)<1$

• 根据矩阵与其Jordan标准形及特征值的关系,可知 $\lim_{k\to\infty}B^{k+1} \Leftrightarrow B的所有特征值的绝对值小于1$

又因为矩阵的谱半径不超过其任一种算子范数,即$\rho(B)<||B||_v$ , 可得

定理 (充分条件) 若存在一个矩阵范数使得$\parallel B\parallel=q<1$, 则迭代收敛，且有下列误差估计：

$$①\varepsilon^{(k)} = \parallel\vec{x}^*-\vec{x}^{(k)}\parallel\leq\frac q{1-q}\parallel\vec{x}^{(k)}-\vec{x}^{(k-1)}\parallel\\ ②\parallel\vec{x}-\vec{x}^{(k)}\parallel\leq\frac{q^k}{1-q}\parallel\vec{x}^{(1)}-\vec{x}^{(0)}\parallel$$

证明:

• &\begin{aligned}①\vec{x}^*-\vec{x}^{(k)}&=B(\vec{x}^*-\vec{x}^{(k-1)})\\&=B(\vec{x}^*-\overline{x}^{(k)}+\vec{x}^{(k)}-\vec{x}^{(k-1)})\\\Rightarrow\parallel\vec{x}^*-\vec{x}^{(k)}\parallel&\leq q(\parallel\vec{x}^*-\vec{x}^{(k)}\parallel+\parallel\vec{x}^{(k)}-\vec{x}^{(k-1)}\parallel)\quad\checkmark\end{aligned}

• $\begin{aligned}& ②\vec{x}^{(k)}-\vec{x}^{(k-1)}=B(\vec{x}^{(k-1)}-\vec{x}^{(k-2)})=...=B^{k-1}(\vec{x}^{(1)}-\vec{x}^{(0)})\\ &\Rightarrow\parallel\vec{x}^{(k)}-\vec{x}^{(k-1)}\parallel\leq q^{k-1}\parallel\vec{x}^{(1)}-\vec{x}^{(0)}\parallel\end{aligned}$

例:

• (1) 求Jacobi法的迭代矩阵

  • $\begin{aligned}B_J=D^{-1}\big(L+U\big)&=\begin{pmatrix}1&0&0\\0&1&0\\0&0&1\end{pmatrix}\cdot\begin{pmatrix}0&-2&2\\-1&0&-1\\-2&-2&0\end{pmatrix}=\begin{pmatrix}0&-2&2\\-1&0&-1\\-2&-2&0\end{pmatrix}\end{aligned}$
  • 显然, $B_J$的常用算子范数$||B_J||>1$, 不能用范数, 只能用特征值来判断
  • $\begin{aligned}\det(\begin{array}{c}\lambda I-B_{J}\end{array})&=\det\begin{pmatrix}\lambda&2&-2\\1&\lambda&1\\2&2&\lambda\end{pmatrix}=\lambda^{3}=0\end{aligned}$
  • 所以$\quad\lambda=0\quad\rho(B_{J})=\max(\mid\lambda\mid)=0<1$, Jacobi迭代法收敛

• (2) 求Gauss-Seidel法的迭代矩阵

  • $B_G=(D-L)^{-1}U=\begin{pmatrix}1&0&0\\1&1&0\\2&2&1\end{pmatrix}^{-1}\cdot\begin{pmatrix}0&-2&2\\0&0&-1\\0&0&0\end{pmatrix}=\begin{pmatrix}0&-2&2\\0&2&-3\\0&0&2\end{pmatrix}$
  • $\begin{aligned}&\lambda=0, \lambda=2\\&\rho(B_{G})=\max(|\lambda|)=2>1\end{aligned}$ , 所以Gauss-Seidel迭代法发散

• Gauss-Seidel迭代法收敛速度快, 但是有可能发散, 并不一定优于Jacobi迭代法

另外,给出系数矩阵对角占优线性方程组的一个结论

定理: 若线性方程组$Ax=b$的系数矩阵$A$为严格对角占优矩阵，则Jacobi法和$G-S$法均收敛

• 矩阵严格对角占优: $\mid a_{ii}\mid>\sum_{j\neq i}\mid a_{ij}\mid\quad i=1,2,3,\cdots,n\\$ 对角线上的值比行上的其他值加起来都大

  • 可得 $\frac1{\mid a_{ii}\mid}\sum_{j\neq i}\mid a_{ij}\mid<1\quad i=1,2,3,\cdots,n\\$

• 对于Jacobi迭代法,其迭代矩阵为$B_J = D^{-1}(L+U)B_J=-\begin{pmatrix}0&\frac{a_{12}}{a_{11}}&\cdots&\frac{a_{1n}}{a_{11}}\\\frac{a_{21}}{a_{22}}&0&\cdots&\frac{a_{2n}}{a_{22}}\\\vdots&\vdots&\ddots&\vdots\\\frac{a_{n1}}{a_{nn}}&\frac{a_{n2}}{a_{nn}}&\cdots&0\end{pmatrix}$

  • $\left\|B_J\right\|_\infty=\max_i\frac1{\mid a_{ii}\mid}\sum_{j\neq i}\mid a_{ij}\mid<1\\$, 所以Jacobi迭代法收敛

• 对于G—S迭代法,其迭代矩阵为$B_G=(D-L)^{-1}U$, $B_G$的形式不像$B_J$那样好确定, 不用范数, 使用特征值判断

  • $$\begin{aligned} B_G\text{的特征值}\lambda\text{满足 }\det(\lambda I-B_G)&=0 \\ \text{即}\det[\lambda I-(D-L)^{-1}U]& =0 \\ \text{从而 }\det(D-L)^{-1}\cdot\det[\lambda(D-L)-U]& =0 \\ \text{因此 }\det[\lambda(D-L)-U]& =0 \\ 由|a_{ii}|>\sum_{j\neq i}|a_{ij}| 得 |\lambda|\cdot|a_{ii}|>|\lambda|\cdot\sum_{j=1}^{i-1}|a_{ij}|+|\lambda|\cdot\sum_{j=i+1}^n|a_{ij}| \\ \qquad=\mid\lambda\mid\cdot\sum_{j=1}^{i-1}\mid a_{ij}\mid+\sum_{j=i+1}^n\mid a_{ij}\mid+(\mid\lambda\mid-1)\cdot\sum_{j=i+1}^n\mid a_{ij}\mid \end{aligned}$$

第7章 非线性方程与方程组的数值解法

7,1

• 二分法和误差

7.2

• 不动点概念和公式
• 不动点收敛性
• 局部收敛
• 收敛的阶

7.3 skip

7.4 牛顿

• 简化牛顿 牛顿下山skip
• 重根情形

7.1 方程求根与二分法

7.1.1 多项式基础

略, 高数

7.1.2 二分法

零点定理: 若$f \in C[a, b]$，且 f (a) · f (b) < 0，则 f 在 (a, b) 上必有一根。

算法描述

设[a,b]为单根区间

取中点$x_0=\frac12(a+b)$

• 若$f(x_0)=0$, $x_0$记为$[a,b]$中的根

• 若$f(a)\cdot f(x_0)<0$, 则$[a, x_0]$为有根区间, 令$a_1=a, b_1=x_0$

• 若$f(x_0)\cdot f(b)<0$, 则$[x_0, b]$为有根区间, 令$a_1=x_0, b_1=b$

一轮操作后, 有根区间$[a,b]$缩小为一半$[a_1, b_1]$

循环, 继续取$[a_1, b_1]$的中点$x_1=\frac12(a_1+b_1)$, 得到新区间$[a_2, b_2]$

以此类推, 有$x_n = \frac12(a_n+b_n)$

对于每个小区间都有$(b_n-a_n)=\frac1{2^{n+1}}(b-a)\\|x_n-x_{n-1}|=\frac{1}{2^{n+1}}(b-a)$

确定适当的n, 可以得到任意要求的精度 ($|x_{k+1}-x_k|<\varepsilon_1$)

误差分析

对于第0步的$x_0=\frac12(a+b)$有误差$|x_0-x^*|\leq\frac{b-a}2$

第k步的$x_k$误差:$|x_k-x^*|\leq \frac{b-a}{2^{k+1}}$

对于给定的精度ε,可估计二分法所需的步数 k ：$\frac{b-a}{2^{k+1}}<\varepsilon\quad\Rightarrow\quad k>\frac{\left[\ln\left(b-a\right)-\ln\varepsilon\right]}{\ln2}-1\\$

特点

优点

• 简单;
• 对f (x) 要求不高(只要连续即可)

缺点

• 无法求复根及偶重根
• 收敛慢

注：用二分法求根，最好先给出$f(x)$ 草图以确定根的大概位置。或用搜索程序，将[a, b]分为若干小区间，对每一个满足$f (a_k)·f (b_k) < 0$的区间调用二分法程序，可找出区间[a, b]内的多个根，且不必要求$f (a)·f (b) < 0$​。

例

例 证明1-x-sin $x=0$在[0,1]内仅有一个根，使用二分法求误差不大于$\frac12\times10^{-4}$​的根需要对分多少次？

解:

• 设$f(x)=1-x-\sin x$,则$f(0)=1>0$,f$(1)=-\sin(1)<0$,且$f(x)$在上[0,1]上连续，故方程$f(x)=0$在[0,1]内至少有一根
• 又因为$f^\prime(x)=-1-\cos x<0,x\in[0,1]$,故$f(x)$在[0,1]上单调递减，因此$f(x)$在[0,1]上有且仅有一个根
• 使用二分法，使误差限$\left|x_k-x^*\right|\leq\frac1{2^{k+1}}(b-a)=\frac1{2^{k+1}}\leq\frac12\times10^{-4}$,解得 $2^{k}\geq10^{4},k\geq\frac{4\ln10}{\ln2}=13.2877$​

所以对分14次即可

7.2 迭代法及其收敛性

不动点迭代法

将非线性方程$f(x)=0$化为一个同解方程$x=\varphi(x)$, 且设$\varphi(x)$连续

任取初值$x_0$代入, 得$x_1=\varphi(x_0), x_2=\varphi(x_1)\cdots, x_k=\varphi(x_{k-1})$

称上式为求解非线性方程$x=\varphi(x)$的不动点迭代法

• 称$\varphi(x)$为迭代函数, $x_k$​为第k步迭代值
• 若存在$x^*,使得lim_{k→∞}x_k=x^*\\$, 则$x^*$为$\varphi$的不动点, 也就是方程f=0的根

发散与收敛

若存在一点$x^*$, 使得迭代序列$\{x_k\}_0^\infin$满足$lim_{k→\infin}x_k=x^*\\$则称迭代法收敛,否则称为发散

如果将原方程表示为$\{\begin{array}{l}y=x\\y=\varphi(x)\end{array}$与原方程同解

例题: 用迭代法求解方程$2x^3-x-1=0$

• 将原方程化为等价方程$x=2x^3-1$

  • 取$x_0=0$, 得$\begin{aligned} &x_{0}=0 \\ &x_1=2x_0^3-1=-1 \\ &x_2=2x_1^3-1=-3 \\ &x_3=2x_2^3-1=-55 \end{aligned}$, 显然发散

• 如果将原方程化为等价方程$x=\sqrt[3]{\frac{x+1}2}$

  • $\begin{aligned}&x_0=0\\&x_1=\sqrt[3]{\frac{x_0+1}2}=\sqrt[3]{\frac12}\approx0.7937\\&x_2=\sqrt[3]{\frac{x_1+1}2}=\sqrt[3]{\frac{1.7937}2}\approx0.9644\end{aligned}$

  • 得$x2 = 0.9644; x3 = 0.9940;x4 = 0.9990;x5 = 0.9998;x6 = 1.0000;x7 = 1.0000$

  • 收敛, 原方程的解为$x=1.0000$

迭代法收敛定理

定理: 设迭代函数$\varphi(x)$在[a,b]上连续, 且满足

1. 当$x\in[a,b]$时,$\ a\leq\varphi(x)\leq b$
2. 存在一个整数L, 满足$0<L<1$且$\forall x\in [a,b]$, 有$|\varphi'(x)|\leq L$(上确界为L)

则有以下结论

1. 方程$x=\varphi(x)$在$[a,b]$内有唯一解$x^*$
2. 对于任意初值$x_0\in[a,b]$, 迭代法$x_{k+1}=\varphi(x_k)$均收敛于$x^*$
3. $\left|x_k-x^*\right|\leq\frac L{1-L}\left|x_k-x_{k-1}\right|\\$
4. $\left|x_k-x^*\right|\leq\frac{L^k}{1-L}\left|x_1-x_0\right|\\$
   • $\left|x_k-x^*\right|\leq\frac L{1-L}\left|x_k-x_{k-1}\right|\leq\frac{L^k}{1-L}\left|x_1-x_0\right|\\$​

定理证明:

结论1: 方程$x=\varphi(x)$在$[a,b]$内有唯一解$x^*$

1. ①构造了一个$(a,a)→(b,b)$的正方形空间, 对角线$y=x$; $f(b)=\varphi(a)-a \geq0; f(b)=\varphi(b)-b\leq0$
2. ②对于两个函数交点:$f(x)=\varphi(x)-x= 0$, $f'(x)=\varphi'(x)-1<0$, 单调递减
3. 零点定理得$\varphi(x)-x$在$(a,b)$有唯一零点$x^*$

结论2: 对于任意初值$x_0\in[a,b]$, 迭代法$x_{k+1}=\varphi(x_k)$均收敛于$x^*$

• 对$x_{k+1}=\varphi(x_k)$, 由微分中值定理和$|\varphi'(x)|<L$得
• $\begin{aligned}\left|x_{k}-x^{*}\right|& =\left|\varphi(x_{k-1})-\varphi(x_{k})\right| \\&\leq L|x_{k-1}-x^*|\leq\cdots\leq L^k|x_0-x^*| \end{aligned}$

结论3&4: $\left|x_k-x^*\right|\leq\frac L{1-L}\left|x_k-x_{k-1}\right|\\$ 和$\left|x_k-x^*\right|\leq\frac{L^k}{1-L}\left|x_1-x_0\right|\\$

• $\begin{aligned}\left|\begin{array}{c}x_{k+p}-x_k\end{array}\right|& \leq\left|\begin{array}{c}x_{k+p}-x_{k+p-1}\end{array}\right|+\left|\begin{array}{c}x_{k+p-1}-x_{k+p-2}\end{array}\right| +\cdots+\left|\begin{array}{c}x_{k+1}-x_k\end{array}\right| \\ &\leq\left(L^p+L^{p-1}+\cdots+L\right)|\left.x_k-x_{k-1}\right| \\&\leq\frac L{1-L}\Big|x_k-x_{k-1}\Big| (等比数列\frac{L-L^{p+1}}{1-L})\end{aligned}$

• $\begin{aligned}\left|x_k-x^*\right|&\leq\frac L{1-L}\Big|x_k-x_{k-1}\Big|\\&\leq\frac{L^2}{1-L}\Big|x_{k-1}-x_{k-2}\Big|\\&\cdots\cdots\cdots\\&\leq\frac{L^k}{1-L}\Big|x_1-x_0\Big|\end{aligned}$​

• 另一种方式证明收敛: 由于L<1,$\lim_{k\to\infty}(|x_k-x^*|)=0$

由定理可知, 只要迭代函数满足$|\varphi'(x)|\leq L<1$​, 迭代法收敛

• 注: 虽然收敛, 但不一定是唯一根

对于实际要求中的误差限$\varepsilon$, 只要$\frac L{1-L}|x_k-x_{k-1}|<\varepsilon\\$

• 所以当$\left|x_k-x_{k-1}\right|<\frac{1-L}L\varepsilon\approx\varepsilon\\$时, 即可停止迭代, 将$x_k$作为近似解

收敛的阶

由定理1的结论可以看出,L或$|\varphi'(x)|$在[a, b]上越小, 迭代法收敛越快

设$e_k=|x_k-x^*|$

定义: 若存在实数$p\geq1$和$c>0$满足$\lim_{k\to\infty}\frac{e_{k+1}}{e_k^p}=c\\$($e_{k+1}和e_k^p$同阶无穷小), 则称迭代法p阶收敛

• 当p=1称为线性收敛, p>1时称为超 线性收敛, p=2时称为平方收敛

收敛阶的计算

若$\varphi(x)$在精确解$x^*$处处可导, 泰勒展开得: $\varphi(x)=\varphi(x^*)+\varphi'(x^*)(x-x^*)+\frac{\varphi''(x^*)}{2!}(x-x^*)^2+...\\+\frac{\varphi^{(p-1)}(x^*)}{(p-1)!}(x-x^*)^{p-1}+\frac{\varphi^{(p)}(\xi)}{p!}(x-x^*)^p\\$

$\text{如果 }\varphi^{\prime}(x^*)=\varphi^{\prime\prime}(x^*)=\cdots=\varphi^{(p-1)}(x^*)=0\quad\text{ 而 }\varphi^{(p)}(x^*)\neq0\\\varphi(x)=\varphi(x^*)+\frac{\varphi^{(p)}(\xi)}{p!}(x-x^*)^p$

• $x_{k+1}=\varphi(x_k)=\varphi(x^*)+\frac{\varphi^{(p)}(\xi)}{p!}(x_k-x^*)^p\\$
• $x_{k+1}-x^*=\frac{\varphi^{(p)}(\xi)}{p!}(x_k-x^*)^p\\$
• $\frac{e_{k+1}}{e_k^p}=\frac{\left|x_{k+1}-x^*\right|}{\left|x_k-x^*\right|^p}\quad\to\left|\frac{\varphi^{(p)}(x^*)}{p!}\right|,(k\to\infty)\\$, 即迭代法$x_{k+1}=\varphi(x_k)$的收敛阶是p

定理: 如果迭代法迭代函数 $\varphi(x)$在根$x^*$附近满足：

• $(1)\varphi(x)$存在$p$​阶导数处连续；
• $(2)\varphi^{\prime}(x^*)=\varphi^{\prime\prime}(x^*)=\cdots=\varphi^{(p-1)}(x^*)=0$,而$\varphi^{(p)}{(x^*)}\neq0$

则迭代法$x_k+1=\varphi(x_k)$的收敛阶是$p$

7.4 牛顿法

原理：将非线性方程线性化 —— Taylor 展开 /Taylor’s expansion/

取$x_0=x^*$, 将f(x)在x0处做一阶泰勒展开,得

$f(x)=f(x_0)+f^{\prime}(x_0)(x-x_0)+\frac{f^{\prime\prime}(\xi)}{2!}(x-x_0)^2,\xi\text{ 在 }x_0\text{ 和 }x\text{ 之间。}$

将将 $(x* - x_0)^2$ 看成高阶小量，则有：

$0=f(x^*)\approx f(x_0)+f^{\prime}(x_0)(x^*-x_0)\quad\Rightarrow\quad x^*\approx(x_0-\frac{f(x_0)}{f^{\prime}(x_0)})\\$ (去掉了f''项, 所以约等于)

$x_{k+1}=x_k-\frac{f(x_k)}{f^{\prime}(x_k)}\\$:

推导过程

如果将非线性方程$f(x)=0$化为等价方程 $x=x-k(x)f(x)$ , 且$k(x)\neq0$​

求k(x)

令$\begin{aligned}\varphi(x)&=x-k(x)f(x)\\\varphi'(x)&=1-k'(x)f(x)-k(x)f'(x)\end{aligned}$

• 设$x^*$为f(x)=0的根, $|\varphi^{\prime}(x)|$在$x^*$​附近越小,则收敛速度越快 (f(x)对x的影响大)

• 如果$f'(x^*)\neq 0$, 令$\varphi'(x^*)=0$(收敛速度最大), 即$1-k'(x^*)f(x^*)-k(x^*)f'(x^*)=0$, $k(x^*)=\frac1{f'(x^*)}$

取$k(x)=\frac1{f'(x)}$, 则$\varphi(x)=x-\frac{f(x)}{f^{\prime}(x)}\quad x=x-\frac{f(x)}{f^{\prime}(x)}$

取初值$x_0$, 由上面k(x)构造迭代函数得

牛顿迭代法: $x_{k+1}=x_k-\frac{f(x_k)}{f'(x_k)}\quad(k=0,1,2,\cdots)\\$

• 具有局部收敛性, 只要$f'(x^*)\neq 0$, 牛顿迭代法至少平方收敛

例1: 用Newton迭代法求方程的根:

例: $\text{设}x^*\text{是方程}f(x)=0\text{的}m(\geq2)\text{重根},\text{证明牛顿法}x_{k+1}=x_k-\frac{f(x_k)}{f'(x_k)}\\$为线性收敛

• 因为$x^*$是方程$f(x)=0$的m重根, $\therefore f(x)=(x-x^*)^mg(x)$且$g(x^*)\neq 0, m\geq2$
• 得$f^{\prime}(x)=m(x-x^*)^{m-1}g(x)+(x-x^*)^mg^{\prime}(x)$
• $\begin{gathered} x_{k+1}=x_k-\frac{f(x_k)}{f^{\prime}(x_k)} =x_k-\frac{(x_k-x^*)^mg(x_k)}{m(x_k-x^*)^{m-1}g(x_k)+(x_k-x^*)^mg^{\prime}(x_k)} \\ =x_k-\frac{(x_k-x^*)g(x_k)}{mg(x_k)+(x_k-x^*)g^{\prime}(x_k)} \end{gathered}$​
• $x_{k+1}-x^*=(x_k-x^*)(1-\frac{g(x_k)}{mg(x_k)+(x_k-x^*)g^\prime(x_k)})\\$
• $\lim_{k\to\infty}\frac{x_{k+1}-x^*}{x_k-x^*}=\lim_{k\to\infty}(1-\frac{g(x_k)}{mg(x_k)+(x_k-x^*)g^{\prime}(x_k)})=1-\frac1m\\$
• $m\geq2\text{时,1}-\frac1m>0$ 由定义可知, 该迭代法对m(>=2)重根是线性收敛的

例: 设$f(a)=0$, 且$f'(a)\neq0$证明$x_{k+1}=x_k-\frac{f(x_k)}{f^{\prime}(x_k)}\\$, 至少是平方收敛

• 令$\varphi(x)=x-\frac{f(x)}{f'(x)}$
• 则$\varphi'(x)=1-\frac{[f'(x)]^2-f(x)f''(x)}{[f'(x)]^2}=\frac{f(x)f''(x)}{[f'(x)]^2}\\$
• 所以$\varphi^{\prime}(a)=0$,该迭代法至少平方收敛

收敛的充分条件

设$f\in C^2[a, b]$若

• $f(a)f(b)<0$ (区间端点异号, 有根)
• 在整个$[a,b]$区间上, $f''$不变号, 且$f'(x)\neq 0$ (根唯一)
• 选取$x_0\in[a,b]$使得$f(x_0)f''(x_0)>0$ (产生的序列单调有界, 保证收敛)

则牛顿迭代法产生的序列$\{x_k\}$收敛到$f(x)$在$[a,b]$区间内的唯一根

局部收敛性

设 $f \in C^2[a, b]$，若 x* 为 f (x) 在[a, b]上的单根，且 $f ’(x^*) \neq 0$，则存在 x* 的邻域$B_\delta (x^*)$使得任取初值$x_0\in B_\delta (x^*)$ ，Newton’s Method产生的序列$\{ x_k \}$收敛到x*，且满足$\lim_{k\to\infty}\frac{x^*-x_{k+1}}{\left(x^*-x_k\right)^2}=-\frac{f^{\prime\prime}(x^*)}{2f^{\prime}(x^*)}\\$​

例: 证明牛顿迭代法实际上是一种特殊的不动点迭代 其中$g(x)=x-\frac{f(x)}{f'(x)}$ 则$\left|g^{\prime}(x^*)\right|=\left|\frac{f^{\prime\prime}(x^*)f(x^*)}{f^{\prime2}(x^*)}\right|=0<1\quad\Rightarrow\quad\text{收敛}$