【笔记】计算方法
课程《计算方法》笔记
参考书目：《数值分析原理》（吴勃英，王德明，丁效华，李道华）
非线性方程和方程组的数值解法
非线性方程：
我们这里只考虑求实根、单根，首先需要满足零点存在定理：\(f(x)\) 在 \([a, b]\) 上连续，且 \(f(a)f(b)<0\)，并且我们假设 \((a, b)\) 上只有一个根 \(\alpha\)
二分法
初始有 \([a _0, b _0]\)，每次取 \([a _n, b _n]\) 中点 \(x _0\) 判断 \(f(a _n)f(x _0)\) 符号，得到 \([a _{n+1}, b _{n+1}]\)
直到恰好为零点，或者与根的差小于 \(\epsilon\)：
迭代法
原式 \(f(x)=0\)，总能转化成 \(x=φ(x)\) 的形式，称 \(φ(x)\) 为“迭代函数”，比如 \(x+f(x),\ x+kf(x),\ x+g(x)f(x)\)
取初值 \(x _0\)，迭代 \(x _{n} = φ(x _{n-1})\)，得数列 \(x _n\)
若数列收敛 \(x _n \to α\)，即满足 \(α=φ(α)\)，则此为原式的解
定理：收敛的充分条件
若 \(φ(x)\in [a, b], x\in[a, b]\)，且满足
则迭代会收敛于唯一根 \(\alpha\)
存在性证明：连续函数 \(g(x)=x-φ(x)\)，继承原方程的符号有 \(g(a)\le 0, g(b)\ge 0\)，则存在 \(\alpha\in[a, b]\) 为根
唯一性证明：反证，假设 \(α _1 ≠ α _2\)，由于 \(α _1-α _2=φ(α _1)-φ(α _2)=φ'(ξ)(α _1-α _2)\)，同取绝对值，代入 \(φ'≤L\)，得 \(L≥1\)，矛盾！
收敛性证明：考虑
\(|\alpha-x _{n+1}| = |φ(\alpha)-φ(x _n)|=|φ'(ξ)||\alpha-x _n|\le L|\alpha-x _n|\le L^2 |\alpha-x _{n-1}|\le \dotsb\)，则 \(\lim |\alpha-x _n|=0\)
或者：考虑级数 \(x _0, x _1-x _0, …, x _n-x _{n-1},…\)，级数收敛等价于前 \(n\) 项和 \(S _n=x _n\) 有限
可以考虑充分条件，证明绝对收敛，即绝对值之和 \(|x _0|+|x _1-x _0|+…\) 有限；用拉氏中值\[x _{n+1}-x _n=φ(x _n)-φ(x _{n-1})=φ'(ξ)(x _n-x _{n-1}) \]两边同取绝对值，再代入 \(|φ'(x)|≤L\)，得 \(|x _{n+1}-x _n|\le L|x _n-x _{n-1}|\)，代回求和式：\(|x _0|+|x _1-x _0|+…\le |x _0| + |x _1-x _0|(1+L+L^2+…)\)，故 \(L<1\) 可令其有限
牛顿迭代法（切线法）
取 \((x _0, f(x _0))\)，切线方程 \(y-f(x _n)=f'(x _n)(x-x _n)\) 代入 \(y=0\) 迭代，即
它还有一个含义，考虑二阶泰勒展开 \(0=f(x)=f(x _k)+f'(x _k)(x-x _k)+o\)，丢掉余项，即得该递推式
它形如迭代法的 \(φ(x)=x+gf,g=-1/f'\)，所以收敛条件可以为：\(|φ'| = |\frac{f * f''}{(f')^2}| ≤L<1\)
定理：收敛的一个充分条件
若 \(f\) 在有根区间 \([a, b]\) 二阶导存在，且 \(f(a)f(b)<0, f'(x)\neq 0, f''(x)\) 不变号，初值 \(x _0\) 有 \(f(x _0)f''(x _0)>0\)，则 \(x _n\) 收敛于唯一根
例：设 \(a>0\)，试建立 \(1/a\) 的牛顿迭代式，要求不含有除法运算，且初值满足 \(0<x _0<2/a\) 时收敛
解：
\(f(x) = x-\frac{1}{a} = 0,f'=1\)，这是个线性式子，迭代就..寄
\(f(x)=a-\frac{1}{x}=0,f'=\frac{1}{x^2}\)，迭代式 \(x _{k+1}=x _k(2-ax _k)\) 可行
现考虑 \(r _k=1-ax _k\)，则 \(r _{k+1} = r _k^2 = \dotsb=r _0^{2^k}\)，由于 \(0<x _0<2/a\) 有 \(-1<r _0<1\)，故 \(r _k\to 0\)，收敛得证
延伸：矩阵求逆，可以考虑 \(X _{k+1}=X _k(2E-AX) \to A^{-1}\)
收敛阶
定义：若有 \(x _k \to α\)，定义 \(ε _k = α - x _k\)，若有下式成立：
则称此迭代是 \(p\) 阶收敛，易见 \(p\) 越大收敛所需次数越少
其中，\(p=1\) 称为线性收敛，\(p=2\) 称为平方收敛
考察式子：\(ε _{k+1}=φ(α)-φ(x _k)=φ'(ξ)ε _k\)
于是有 \(\lim \frac{ε _{k+1}}{ε _k}=\lim φ'=φ'(α)\)，可见若 \(φ'(α)≠0\) 则线性收敛，\(φ'(α)=0\) 则至少可以是平方收敛。事实上有定理：
定理 设 \(φ(x)\) 在 \(x=φ(x)\) 的根 \(\alpha\) 邻域内有充分多阶导数，则 \(p\) 阶收敛充要条件是：
充分性证明：泰勒展开，代入，得：
\[x _{i+1}=φ(x _i)=φ(\alpha) + \frac{1}{p!}φ^{(p)}(\xi)(x _i-\alpha)^p,\quad \xi\in (\alpha, x _i) \]代入定义式，即
\[\lim \frac{|\epsilon _{i+1}|}{|\epsilon _i|^p} = \lim \frac{|x _{i+1}-\alpha|}{|x _i-\alpha|^p} = \frac{1}{p!}|φ^{(p)}(\alpha)|\neq 0 \]
而牛顿法 \(x _{k+1} = x _k - \frac{f(x _k)}{f'(x _k)}\) 是平方收敛的，因为：
收敛效率
刻画计算效率，称 \(EI = p^{\frac{1}{\theta}}\) 为效率指数，其中 \(\theta\) 每次迭代的计算量，\(p\) 为收敛阶；\(EI\) 越大，则计算效率越高
重根的迭代法 - 牛顿迭代法和修正
\(\alpha\) 是 \(f(x)\) 的 \(r\) 重根，又即 \(f(x)=(x-\alpha)^r g(x)\)
继续考虑牛顿迭代法，\(f'(x)=r(x-\alpha)^{r-1}g(x)+(x-\alpha)^r g'(x)\)，求导下去会发现 \(f(\alpha)=f'(\alpha)=\dotsb=f^{(r-1)}(\alpha)=0\)，之前关于导数的推导不能用了，但是我们还是考虑牛顿迭代公式：
且 \(r>1, |φ'(\alpha)|<1\)，故为一阶收敛
可以考虑修正，若已知 \(r\)，递推式为：
则其是二阶收敛的
证明：见 (1.5.3)，总之经典定义+泰勒展开
若 \(r\) 未知，考虑事实：\(f'\) 的重根数比 \(f\) 少一，则设 \(u(x)=f(x)/f'(x)\)，则 \(u(x)\) 具有单根 \(\alpha\)，故递推式为：
例题 1.3 求 \((\sin x-\frac{x}{2})^2=0\) 的正根（0 处重根）
用上面三个方法都行；提一嘴，判断根数，可以考虑求导后代入，到几阶导不为零就有几个根
多点迭代法 - 两点迭代（割线法）
用割线替代切线，迭代公式为
考虑收敛阶，有结论：\(p=\frac{1+\sqrt{5}}{2}\)（有一些条件，见书 定理 1.6）
非线性方程组 - 拟牛顿法
非线性方程组 \({\bf F(x)=0, F(x)} = (f _1({\bf x}), \dots, f _n({\bf x}))^T, {\bf x}\in {\mathbb R}^n\)，我们一样考虑泰勒展开取前两项 \(f(\mathbf{x})=f\left(\mathbf{x} _{k}\right)+\left[\nabla f\left(\mathbf{x} _{k}\right)\right]^{T}\left(\mathbf{x}-\mathbf{x} _{k}\right)+o^{n}\)，得到递推式：
其中
线性方程组的数值解法
研究 \(n\) 阶线性方程组
主要有几类方法：直接法、迭代法、梯度法（方程解转化成最小值问题）
直接法：方程转化成 \({\bf G x = d}\)，\(\bf G\) 为简单矩阵（比如上三角）
迭代法：方程转化成 \(\bf x={\rm G}(x) = Bx + d\)，迭代
迭代法
Jacobi 迭代
设 \(\bf A\) 非奇异，且 \(a _{ii}\ne 0\)，则可将方程写为：
视为迭代式，即
其为 Jacobi 迭代法
从矩阵来看，可将 \(\bf A\) 分解成 \(\bf A = D + L + U\)，分别为 \(\bf A\) 的对角线元素，和对角线为零的下三角、上三角阵
则有：
Gauss-Seidel 迭代
对 Jacobi 迭代 \((1)\) 稍加修改，总是代入新算出的分量：
其为 Gauss-Seidel 迭代法
考虑矩阵，可以写成：
超松弛迭代法（SOR 方法）
为了加快收敛速度，变化 Gauss-Seidel 迭代公式 \((3)\)：
记分量误差 \(r^{(k+1)} _i = x _i^{(k+1)} - x _i^{(k)} = \frac{1}{a _{ii}}(b _i - \sum _{j=1}^{i-1} a _{ij}x _j^{(k+1)}- \sum _{j=i}^{n} a _{ij}x _j^{(k)})\)，引入松弛因子 \(\omega\)，则新的迭代公式为：
矩阵形式：
向量范数
记空间 \(\mathbb R\) 或 \(\mathbb C\) 为 \(\mathbb K\)，定义：
向量范数公理
任意向量 \(\bf x\in \mathbb K^n\) 对应一个非负实数 \(\bf \lVert x\rVert\)；对任意 \(\bf x, y\in \mathbb{K}^n\)，满足：
- 非负性：\(\bf \lVert x\rVert\ge 0,\quad \lVert x\rVert=0\Leftrightarrow x=0\)
- 齐次性：\(\bf \lVert \alpha x\rVert = | \alpha| \lVert x\rVert\)
- 三角不等式：\(\bf \lVert x+y\rVert\le \bf \lVert x\rVert + \bf \lVert y\rVert\)
则称 \(\bf \lVert x\rVert\) 为 \(\bf x\) 的范数
常见范数有：
- 1 范数 \({\bf \lVert x\rVert} _1 = \sum |x _i|\)
- 2 范数 \({\bf \lVert x\rVert} _2 = \sqrt{\sum |x _i|^2}\)
- \(∞\) 范数 \({\bf \lVert x\rVert} _∞ = \max |x _i|\)
- 一般的 \(p\) 范数 \({\bf \lVert x\rVert} _p = (\sum|x _i|^p)^{1/p}\)
范数的等价性定理：
\(\forall p,q\)，\(p\) 范数和 \(q\) 范数是等价的，含义为：
矩阵范数
矩阵范数公理
任意矩阵 \(\bf A\in \mathbb K^{n\times n}\) 对应一个非负实数 \(\bf \lVert A\rVert\)；对任意 \(\bf A, B\in \mathbb{K}^{n\times n}\)，满足：
- 非负性：\(\bf \lVert A\rVert\ge 0,\quad \lVert A\rVert=0\Leftrightarrow A=0\)
- 齐次性：\(\bf \lVert \alpha A\rVert = | \alpha| \lVert A\rVert\)
- 三角不等式：\(\bf \lVert A+B\rVert\le \bf \lVert A\rVert + \bf \lVert B\rVert\)
则称 \(\bf \lVert A\rVert\) 为 \(\bf A\) 的范数
定义 矩阵范数与向量范数相容，若满足：\(\bf {\lVert Ax\rVert}\le {\lVert A\rVert}{\lVert x\rVert}\)，则称此矩阵范数与向量范数相容
由此我们继续定义：
定义 矩阵的范数
称 \(\bf \lVert A \rVert = \max _{\lVert x \rVert \ne 0}\frac{\lVert Ax \rVert}{\lVert x \rVert} = \max _{\lVert x \rVert = 1}\lVert Ax \rVert\) 为矩阵 \(\bf A\) 的范数，也称其为向量范数产生的从属范数；易见，从属范数一定与所给定的向量范数相容
可以证明该定义符合上述三个公理（用左定义式去证公理，而右式由 \(\bf \max _{\lVert x \rVert \ne 0}\frac{\lVert Ax \rVert}{\lVert x \rVert} = \max _{\lVert x \rVert \ne 0}\lVert A\frac{x}{\lVert x \rVert}\rVert\) 可得）
对应地，定义 \(p\) 范数为 \({\bf \lVert A \rVert} _p = \max _{{\bf \lVert x \rVert} _p \ne 0}\frac{{\bf \lVert Ax \rVert} _p}{{\bf \lVert x \rVert} _p}\)
在该定义下，可以得到若干常见范数：
- 1 范数/“列范数” \({\bf \lVert A\rVert} _1 = \max _{1\le j\le n}\sum _{i=1}^{n}|a _{ij}|\)
- 2 范数 \({\bf \lVert A\rVert} _2 =\sqrt{\rho({\bf A^T A})} = \sqrt{\max _i |\lambda _i|}\)，\(\lambda _i\) 为对称阵 \(\bf A^T A\) 的第 \(i\) 个特征值
- \(∞\) 范数/“行范数” \({\bf \lVert A\rVert} _\infty = \max _{1\le i\le n}\sum _{j=1}^{n}|a _{ij}|\)
以 \(∞\) 范数 为例，证明思路：
\[\begin{aligned} \bf Ax &= \begin{bmatrix} \sum _j a _{1j}x _j & \sum _j a _{2j}x _j & \dotsb & \sum _j a _{nj}x _j \end{bmatrix}^T \\ { \lVert {\bf A}\rVert} _\infty &= \max _{\lVert {\bf x} \rVert _\infty = 1}\lVert {\bf Ax} \rVert _\infty = \max _{\lVert {\bf x} \rVert _\infty = 1} \max _i |\sum _j a _{ij}x _j| \\ & \le \max _{\lVert {\bf x} \rVert _\infty = 1} \max _i \sum _j |a _{ij}x _j| = \max _{\lVert {\bf x} \rVert _\infty = 1} \max _i \sum _j |a _{ij}| |x _j| \\ & \le \max _i \sum _j |a _{ij}| \quad\quad;|x _j|\le \lVert {\bf x} \rVert _\infty = 1\\ {\bf \lVert A\rVert} _\infty &= \max _{\lVert {\bf x} \rVert _\infty = 1}\lVert {\bf Ax} \rVert _\infty \\ & \ge \lVert {\bf Ax}^{(0)} \rVert _\infty = \max _{i}|\sum _j a _{ij}x^{(0)} _j| \quad\quad; \lVert {\bf x}^{(0)} \rVert _\infty = 1 \\ & = |\sum _j a _{kj}x^{(0)} _j| \quad\quad; k=\arg \max _{i} |\sum _j a _{ij}x^{(0)} _j|\\ &= \sum _j |a _{kj}| \quad\quad; 取 x _j^{(0)}=1 \text{ if } a _{kj}\ge 0 \text{ else } -1 \\ &= \max _i \sum _j |a _{ij}| \quad\quad; 好像 {\bf x}^{(0)} 和 k 互相定义？怪耶 \end{aligned}\]
收敛性判别
收敛记作 \({\bf x}^{(k)}\to {\bf x}^*\)，即 \(\lim _{k\to ∞} x _i^{(k)}= x _i^*,\quad i=1,\dots,n\)
我们也可以将向量收敛用范数刻画，具体来说：
定理：收敛性判别，充要条件
以下三个命题等价：
- 迭代法 \(\bf x={\rm G}(x) = Bx + d\) 收敛
- 谱半径 \(\rho({\bf B})< 1\)
- 存在一种范数使得 \({\lVert \bf B\rVert}< 1\)
不严谨说明：若收敛，有 \(\bf x^* = Bx^ * + d\)，代入递推式：
可见 \({\lVert \bf B\rVert}< 1\) 时，\({\bf x}^{(k+1)} - {\bf x}^{*} \to 0\)
迭代收敛的其他判别方法
引理：若 \({\bf B} \in {\mathbb K}^{n\times n}\) 满足 \({\lVert \bf B\rVert}< 1\)，则 \(\bf I-B\) 可逆
证明：反证，\(\bf I-B\) 不可逆，则 \(\bf (I-B)x=0\) 有非零解 \(\bf x^ *\)，则 \(\bf x^ * = B x^ *,\lVert x^ *\rVert = \lVert B x^ *\rVert\le \lVert B\rVert\lVert x^ *\rVert\implies \lVert B\rVert\ge 1\)，矛盾！
定义 对角占优矩阵
称方阵 \(\bf A\) 为对角占优矩阵，若 \(\bf A\) 满足：\(|a _{ii}|> \sum _{j\ne i}|a _{ij}|,\quad i=1,\dots,n\)
定理：对角占优矩阵 \(\bf A\) 可逆
证明：\(\bf A=L+U+D\)，由于对角占优，将式子展开可证得 \(\lVert\bf I-D^{-1}A\rVert _\infty<1\)，又由上文引理得 \(\bf I-(I-D^{-1}A)=D^{-1}A\) 可逆，故 \(\bf A\) 可逆
定理：对角占优 \(\implies\) Jacobi, Gauss-Seidel 迭代法均收敛
证明 i：Jacobi 迭代式中 \(\bf B = I-D^{-1}A\)，可证得 \(\lVert\bf I-D^{-1}A\rVert _\infty<1\)
证明 ii：Gauss-Seidel 迭代式中 \(\bf B = -(L+D)^{-1}U,\quad \lVert B \rVert _{\infty}=\max _{\lVert x \rVert _\infty =1}\lVert Bx \rVert _{\infty}=\max _{\lVert x \rVert _\infty =1}\lVert y \rVert _{\infty}\)
\(\bf y=Bx = -(L+D)^{-1}U,\quad (L+D)y = -Ux,\quad y = -D^{-1}Ly-D^{-1}Ux\)
\(({\bf Ly}) _i = a _{i1}y _1 + a _{i2}y _2 + \dotsb + a _{i (i-1)}y _{i-1}=\sum _{j=1}^{i-1}a _{ij}y _j, \quad ({\bf D^{-1}Ly}) _i = \frac{1}{a _{ii}}\sum _{j=1}^{i-1}a _{ij}y _j,\quad ({\bf D^{-1}Ux}) _i = \frac{1}{a _{ii}}\sum _{j=i+1}^{n}a _{ij}x _j\)
\(y _i = -\frac{1}{a _{ii}}\sum _{j=1}^{i-1}a _{ij}y _j - \frac{1}{a _{ii}}\sum _{j=i+1}^{n}a _{ij}x _j\)
\(\lVert y \rVert _\infty = \max _i \big|y _i\big| = \big|y _k\big| = \big|\frac{1}{a _{kk}}\sum _{j=1}^{k-1}a _{kj}y _j + \frac{1}{a _{kk}}\sum _{j=k+1}^{n}a _{kj}x _j\big|\le \big|\frac{1}{a _{kk}}\big|\sum _{j=1}^{k-1}\big|a _{kj}\big|\cdot\big|y _j\big| + \big|\frac{1}{a _{kk}}\big|\sum _{j=k+1}^{n}\big|a _{kj}\big|\cdot\big|x _j\big|\le \big|\frac{1}{a _{kk}}\big|\sum _{j=1}^{k-1}\big|a _{kj}\big|\cdot{\bf \lVert y\rVert} _{\infty} + \big|\frac{1}{a _{kk}}\big|\sum _{j=k+1}^{n}\big|a _{kj}\big|;\quad {\bf \lVert x \rVert} _\infty =1, |x _j|<1\)
\(\big( 1- \big|\frac{1}{a _{kk}}\big|\sum _{j=1}^{k-1}\big|a _{kj}\big|\big){\bf \lVert y\rVert} _{\infty}\le \big|\frac{1}{a _{kk}}\big|\sum _{j=k+1}^{n}\big|a _{kj}\big|\)，且由对角占优得 \(1- \big|\frac{1}{a _{kk}}\big|\sum _{j=1}^{k-1}\big|a _{kj}\big|> \big|\frac{1}{a _{kk}}\big|\sum _{j=k+1}^{n}\big|a _{kj}\big|\)
故 \({\bf \lVert y\rVert} _{\infty}< 1\)，故收敛
考虑 SOR，\(\bf B = (D+\omega L)^{-1}(D-\omega D-\omega U)\)
若收敛，至少需满足 \(0<\omega< 2\)，才有可能收敛
例题：若有 \({\bf x} = [1, -1],\quad x _1+2 x _2=-1, 3x _1 + x _2 = 2\)，求迭代式
解：\({\bf A} = \begin{bmatrix} 1 & 2 \\ 3 & 1 \end{bmatrix}\)。若为 Jacobi，\({\bf B = I-D^{-1}A} = \begin{bmatrix} 0 & -2 \\ -3 & 0 \end{bmatrix}\)，求 \(\rho({\bf B})\)，计算特征值：\({\bf |\lambda I-B|=0}, \lambda^2=6, \rho({\bf B}) = \sqrt{6}>1\)，不收敛！同理 Gauss-Seidel 也不收敛
直接交换行，方程等价于 \({\bf A} = \begin{bmatrix} 3 & 1 \\ 1 & 2 \end{bmatrix}\)，此时对角占优，Jacobi, Gauss-Seidel 均收敛
所以说：先看看能不能 \(\bf A\) 直接交换行得到对角占优，实在不行，再看看 \(\rho<1\) 是否成立
Gauss 消元法
行选主元素 Gauss 消元法
方程组：
第一步，假设 \(a _{11}\ne 0\)，令第一行 \(x _1\) 系数为 \(1\)，分别代入后面方程，得
第二步，假设 \(a _{22}^{(1)}\ne 0\)，令第二行 \(x _2\) 系数为 \(1\)，分别代入后面方程，得
以此类推，最后得：
从最后一个方程倒回去递推，称之为（行选主元素）Gauss 消元过程，复杂度 \(\text{O}(\frac{1}{3}n^3)\)
可做条件： \(a _{kk}^{(k-1)}\ne 0,1\le k\le n\)，等价于所有顺序主子式 \(D _k\ne 0\)
因为每次进行消元只进行行变换和除法；而行变换不改变行列式值，对行除一个数使行列式除同一个数
考察 \((a _{ij}),1\le i,j\le k\) 的方阵，消元时其行列式值 \(D _k\) 先后除以 \(a _{11},a _{22}^{(1)},\dotsb,a _{kk}^{(k-1)}\)，最后变为行列式等于 \(1\) 的上三角阵，所以 \(D _k = a _{11}a _{22}^{(1)}\dotsb a _{kk}^{(k-1)}\)，故上述等价关系成立
列选主元素 Gauss 消元法
若 \(a _{kk}^{(k-1)}\) 数量级极小时，按行会出现精度问题，比如：
\(10^{-4}x _1+x _2=1,x _1+x _2=2 \implies x _1+10^4 x _2 = 10^4,(1-10^4)x _2=2-10^4\)，若精度不够导致 \(10^4-1\approx 10^4, 10^4-2\approx 10^4\)，则由最后的方程组递推得 \(x _2=1, x _1=0\)，是错解！
列选主元素做法，也就是第 \(k\) 次时，选择使除数最大的行 \(r _k=\arg \max _{k\le i\le n} |a _{ik}^{(k-1)}|\)，将其与 \(k\) 行做交换，然后消元
回到上面的例子，交换一二行，变成 \(x _1+x _2=2,10^{-4}x _1+x _2=1\implies x _1+x _2=2,(1-10^{-4})x _2=1-2\times 10^{-4}\)，最后解得 \(x _2=1,x _1=1\)
此外还有全选主元的 Gauss 消元法，不过较为少用
三角分解法
若 \(\bf A\) 任意顺序主子式均非零，则从矩阵变换来看，每次对 \(\bf A\) 做消元，相当于左乘下三角矩阵，最后得到一个上三角矩阵 \(\bf L _{n}L _{n-1}\dotsb L _{1}A=U\implies A = LU\)，故 \(\bf A\) 为下三角和上三角的乘积
此时由 \(\bf LUx=b\)，记 \(\bf Ux=y\)，则等价于解方程组 \(\bf Ly=b,Ux=y\)，称之为 三角分解法
对于不一定满足顺序主子式均非零，但是 \(\bf A\) 非奇异的，可以通过行交换排列 \(\bf P\) 得到 \(\bf PA=LU\)，不过这里不考虑了
Crout 克劳特分解方法
对于高斯消元过程，\(\bf L _{n}L _{n-1}\dotsb L _{1}A=U\)，\(\bf U\) 为单位上三角矩阵
- 初始计算：
- \(\bf L\) 第一列：\(l _{i1}=a _{i1},\quad i=1,\dotsb,n\)
- \(\bf U\) 第一行：\(u _{1j}=a _{1j}/l _{11},\quad j=2,\dotsb,n\)
- 递推：
- \(\bf L\) 第 \(r\) 列：\(l _{ir}=a _{ir}-\sum _{k=1}^{r-1}l _{ik}u _{kr},\quad i=r,\dotsb,n\)
- \(\bf U\) 第 \(r\) 行：\(u _{rj}=(a _{rj}-\sum _{k=1}^{r-1}l _{rk}u _{kj})/l _{rr},\quad i=r+1,\dotsb,n\)
- 求解 \(\bf Ly=b\)：
- \(y _1 = b _1 / l _{11}\)
- \(y _i = (b _{i}-\sum _{j=1}^{i-1}l _{ij}y _j)/l _{ii},\quad i=2,\dotsb,n\)
- 求解 \(\bf Ux=y\)：
- \(x _n = y _n\)
- \(x _i = y _{i}-\sum _{j=i+1}^{n}u _{ij}x _j,\quad i=n-1,\dotsb,1\)
Doolittle 杜利特尔分解法
\(\bf A\) 分解成单位下三角和上三角乘积：
解法类同 Crout
矩阵条件数与误差分析
直观感受，比如两个方程组：
\(2x _1 + 6x _2 = 8,2x _1 + 6.00001x _2 = 8\implies x _1 = 4, x _2 = 0\)
\(2x _1 + 6x _2 = 8,2x _1 + 6.00001x _2 = 8.00001\implies x _1 = 1, x _2 = 1\)
现实生活中得到的矩阵不一定是准确的（可能导致不稳定方程组-“病态方程组”），我们来考虑误差导致近似解 \(\bf \tilde{x} = x+\delta x\) 的精度变化
以下默认精确解满足：\(\bf Ax = b\)，且 \(\bf A\) 可逆
- \(\bf A\) 精确，\(\bf b\) 有小的扰动：
\(\bf A(x+\delta x) = b + \delta b\implies \delta x = A^{-1}\delta b\implies \lVert \delta x\rVert = \lVert A^{-1} \delta b\rVert\le \lVert A^{-1}\rVert \lVert \delta b\rVert\)
且 \(\bf \lVert b\rVert = \lVert Ax\rVert \le \lVert A\rVert \lVert x\rVert\implies \lVert x\rVert\ge \lVert b\rVert / \lVert A\rVert\)
合并不等式，得\[\bf \frac{\lVert \delta x\rVert}{\lVert x\rVert}\le \lVert A\rVert \lVert A^{-1}\rVert \frac{\lVert \delta b\rVert}{\lVert b\rVert} \]称 \(\lVert A\rVert \lVert A^{-1}\rVert\) 为 \(\bf A\) 的条件数，记作 \(\bf Cond(A)\) 或 \(k({\bf A})\)
可见，若条件数过大（比如好几个数量级），就有可能误差较大 - \(\bf A\) 有小的扰动，\(\bf b\) 精确
可以证明，\(\bf A + \delta A\) 非奇异的充分条件是 \(\bf \lVert A^{-1}\delta A \rVert<1\) 或 \(\bf \lVert A^{-1}\rVert \lVert\delta A \rVert<1\)，下面假设这两式成立
\(\bf (A + \delta A)(x + \delta x) = b\implies \delta x = -A^{-1}\delta A(x +\delta x)\)
\(\bf \implies \lVert \delta x \rVert\le \lVert A^{-1} \rVert \lVert \delta A \rVert \lVert x \rVert + \lVert A^{-1} \rVert \lVert \delta A \rVert \lVert \delta x \rVert\)
\(\bf \implies (1-\lVert A^{-1} \rVert \lVert \delta A \rVert)\lVert \delta x \rVert\le \lVert A^{-1} \rVert \lVert \delta A \rVert \lVert x \rVert\)
故\[\bf \frac{\lVert \delta x \rVert}{\lVert x \rVert}\le \frac{\lVert A^{-1} \rVert\lVert \delta A \rVert}{1-\lVert A^{-1} \rVert \lVert \delta A \rVert} = \frac{Cond(A)\frac{\lVert \delta A \rVert}{\lVert A \rVert}}{1-Cond(A)\frac{\lVert \delta A \rVert}{\lVert A \rVert}} \]误差也和 \(\bf Cond(A)\) 有关
插值法与数值逼近
给出未知函数 \(f(x)\) 在 \(n+1\) 个插值节点 \(x _0, x _1, \dotsb,x _n\) 上的函数值，寻求一个插值函数 \(y(x)\) 以逼近 \(f(x)\)，使得其在插值节点处的值与 \(f(x)\) 相等
同时，称误差函数 \(E(x) = f(x) - y(x)\) 为“插值余项”；关于插值余项可能会用到的定理：
罗尔定理：\(f(a)=f(b)\implies \exist \xi\in(a,b),f'(\xi)=0\)
插值法 - 多项式插值
形如 \(y = a _m x^m + a _{m-1} x^{m-1} + \dotsb + a _1 x + a _0\)，取 \(m=n\)，转化成 \(n+1\) 个 \(n+1\) 元方程组：\(\bf Ax = b\)，即：
\(\bf A\) 可逆，等价于 \(|{\bf A}|\ne 0\)，其为范德蒙行列式，有 \(|{\bf A}| = \prod _{i\ne j}(x _i - x _j)\ne 0\)，即 \(x _0, \dotsb, x _n\) 互不相等
结论：插值节点互不相等时，插值多项式存在且唯一
Lagrange 插值
先考虑简单情况
一点插值：\((x _0, f(x _0))\)，插值零次多项式 \(L _0(x) = f(x _0)\)
两点插值，插值一次多项式（线性插值）\(L _1(x) = \frac{x-x _1}{x _0 - x _1}f(x _0) + \frac{x-x _0}{x _1 - x _0}f(x _1)\)
观察，现考虑形式
且满足 \(l _j(x _k) = [j=k]\)，那么就是合法的插值
显然我们可以令
记多项式 \(P _{n+1}(x) = (x-x _0)\dotsb(x-x _n)\)，那么拉式插值为：
现考虑插值余项，以一次为例，有：
\(E(x)=f(x)-L _1(x)=C _x(x-x _0)(x-x _1),\quad x _0<x _1,C _x\) 为某个待求的关于 \(x\) 的函数
引入关于 \(t\) 的中间函数 \(g(t)=f(t)-L _1(t)-C _x(t-x _0)(t-x _1)\)（注意到 \(g(t)\) 式中 \(C _x\) 仍是关于 \(x\) 的函数，也就是关于 \(t\) 的常函数，便于我们之后求解 \(C _x\)），固定常数 \(x\in (x _0,x _1)\)，有 \(g(x _0)=g(x _1)=g(x)= 0\)，应用两次罗尔定理可得 \(\exist \xi\in (x _0, x _1), g''(\xi)=0\)；再考虑 \(g''=f''-0-C _x 2!\)，代入 \(\xi\)，得 \(0=g''(\xi)=f''(\xi)-C _x 2!\)，即 \(C _x=\frac{f''(\xi)}{2!}\)
对于 \(n\) 次同理，应用 \(n+1\) 次罗尔定理，得：
称之为拉格朗日插值余项
我们可以用这个式子，说明一下拉式插值函数在多项式线性空间里的意义——基函数表示
代入具体函数：
- 若 \(f(x)=1\)，则有 \(1-\sum _j l _j(x)=0\)
- 若 \(f(x)=x^k,1\le k\le n\)，则有 \(x^k-\sum _j x _j^k l _j(x)=0, x _j\) 是已知插值节点
现考虑次数 \(\le n\) 的多项式全体，构成了一个线性空间，且 \(1,x,\dotsb, x^n\) 为其一组基
由 1,2 式 \(1=\sum _j l _j(x),x ^k = \sum _j x _j^k l _j(x)\)，可见 \(l _0(x),l _1(x),\dotsb, l _n(x)\) 也是一组基——其实这也就是我们上面提到的 \(L _n(x) = f(x _0)l _0(x) + f(x _1) l _1(x) + \dotsb + f(x _n) l _n(x)\)，所以我们也称 \(l _0(x),l _1(x),\dotsb, l _n(x)\) 为 “插值基函数”
Newton 插值
我们希望把 \(n+1\) 点、\(n\) 次插值多项式写成升幂形式：\(y=a _0+a _1(x-x _0)+a _2(x-x _0)(x-x _1)+\dotsb+a _{n}(x-x _0)(x-x _1)\dotsb(x-x _{n-1})\)
代入插值条件 \(y(x _j) = f(x _j)=f _j\)，不断递推即可得 \(a _0=f _0,a _1=(f _1-f _0)/(x _1-x _0),a _2=\frac{(f _2-f _0)/(x _2-x _0)-(f _1-f _0)/(x _1-x _0)}{x _2-x _1},\dotsb\)
我们引入 “差商”：
称 \(f[x _0\ x _1] = \frac{f(x _1)-f(x _0)}{x _1-x _0}\) 为一阶差商
称 \(f[x _0\ x _1\ x _2] = \frac{f[x _0\ x _2]-f[x _0\ x _1]}{x _2-x _0}\) 为二阶差商
以此类推，称
为 \(n\) 阶差商
定理：差商关于点对称，即与 \(x _0,x _1,\dotsb,x _n\) 顺序无关
证明：比较拉式插值和牛顿插值的首项系数，得
\(f[x _0\ x _1\dotsb x _n] = \sum _{j=0}^{n}\frac{f(x _j)}{P' _{n+1}(x _j)}\)
故我们的 \(n+1\) 点、\(n\) 次牛顿插值为：
从基的意义来看，\(1, (x-x _0), (x-x _0)(x-x _1),\dotsb, (x-x _0)\dotsb(x-x _{n-1})\) 为一组基函数
现考虑插值余项，对于某点 \((\bar{x},f(\bar{x}))\)，引入 \(n+1\) 次插值 \(f(\bar x)=N _{n+1}(\bar{x})=N _n(\bar{x})+f[x _0 \dotsb x _n\ \bar{x}] (x-x _0)\dotsb(x-x _n)\)，得余项为：
Hermite 插值，埃尔米特插值
不仅要求插值点处值相等，还要求一阶导相等——相切，即插值条件为：
取 \(2n+1\) 次多项式 \(H(x)\)，称之为（完全的）Hermite 插值多项式（部分插值点处一阶导相等的称为不完全的），我们先考虑完全的
仿照拉式插值形式，设：
满足 \(h _j(x _k) = [j=k],\quad h' _j(x _k) = 0,\quad \bar{h} _j (x _k)=0,\quad \bar{h}' _j(x _k) = [j=k]\)
则称 \(h _j(x),\bar{h} _j(x)\) 为 Hermite 插值基函数
沿用拉插记号 \(P _{n+1}(x) = (x-x _0)\dotsb(x-x _n),\quad l _j(x) = \frac{P _{n+1}(x)}{(x-x _j)P' _{n+1}(x _j)}\)，且 \(l _j(x _k) = [j=k]\) 具有一重根 \(x _0,x _1,\dots,x _{j-1},x _{j+1},\dots,x _n\)；考察各个零点的重数
\(\bar{h} _j(x)\) 具有二重根 \(x _0,x _1,\dots,x _{j-1},x _{j+1},\dots,x _n\) 和一重根 \(x _j\)，故可以令 \(\bar{h} _j(x) = (x-x _j) l _j^2(x)\)
\(h _j(x)\) 具有二重根 \(x _0,x _1,\dots,x _{j-1},x _{j+1},\dots,x _n\)，且要满足 \(h _j(x _j) = 1, h' _j(x _j) = 0\)，待定系数该部分乘式为 \((a _0 + a _1(x-x _j))\)，代入解得 \(h _j(x) = [1-2(x-x _j)l' _j(x _j)]l _j^2(x)\)
故有：
插值余项 \(f(x)-H(x) = C _x (x-x _0)^2(x-x _1)^2\dotsb (x-x _n)^2\)
设 \(g(t) = f(t)-H(t) - C _x (t-x _0)^2(t-x _1)^2\dotsb (t-x _n)^2\)
用 \(2n+2\) 次罗尔定理，得 \(0 = \big(\frac{{\rm d}g}{{\rm d}t}\big)^{2n+2} (\xi) = f^{(2n+2)}(\xi) - 0 - C _x (2n+2)!\)
故 \(C _x = f^{(2n+2)}(\xi) / (2n+2)!\)，故插值余项为：
对于不完全的 Hermite 插值，依然是考虑根的重数，就不写了。
Runge 现象，龙格
考虑无限次可微函数，若取很高的插值多项式系数，可能会发现在插值区间外，函数发生剧烈震荡甚至发散，以至于无法拟合插值区间外的函数部分；故实际应用中，很少使用高次多项式
函数逼近 - 最佳平方逼近
有时甚至观测得到的 \(f(x _i)\) 本身也带有误差，我们也不强求 \(y\) 在离散点处相等，而是选择最小化误差，使得 \(f\) 和 \(y\) 尽量逼近，称 \(y\) 为拟合函数
比如，最小二乘拟合：最小化 \(\sum _{i=0}^{n} (y(x _i) - f(x _i))^2\)，或者加权和：\(\sum _{i=0}^{n} \rho _i (y(x _i) - f(x _i))^2\)
正交多项式
定义 权函数 \(\rho(x)\) 在有限或无限区间 \([a,b]\) 上，称 \(\rho(x)\) 为 \([a,b]\) 上的权函数，若：
- \(\rho(x)\ge 0,x\in [a, b]\)
- \(k=0,1,2,\dotsb,\int _a^b \rho(x)x ^k\) 都存在
- 对于非负的 \(f(x)\in C[a,b]\)，若 \(\int _a^b f\rho {\rm d}x = 0\)，则 \(f(x)\equiv 0\)
定义 带权正交 \(f(x),g(x)\in C[a,b]\)，称 \(f(x),g(x)\) 在 \([a,b]\) 上带权 \(\rho\) 正交，记为 \(f\bot g\)，若两者内积为零：
定义 正交函数族 若函数序列 \(\{φ _j \} _0^{\infty}\) 在 \([a,b]\) 上带权 \(\rho\) 两两正交，则称其为 \([a,b]\) 上带权 \(\rho\) 的正交函数集
给定权函数 \(\rho\)，可以通过正交化方法，构造其多项式序列 \(\{φ _n(x) \} _0^{\infty}\)：
或者用其递推式计算，先算出 \(φ _0, φ _1\)，然后有：
根据不同的权函数和区间，我们有常见正交多项式序列：
Legendre 勒让德多项式，区间 \([-1, 1]\)，权函数 \(\rho(x)=1\)
性质：奇偶性 \(P _n(-x) = (-1)^n P _n(x)\)，即 \(n\) 偶数偶函数，奇数基函数
Chebyshev 切比雪夫多项式，区间 \([-1, 1]\)，权函数 \(\rho(x) = \frac{1}{\sqrt{1-x^2}}\)
性质：奇偶性 \(T _n(-x) = (-1)^n T _n(x)\)
第二类 Chebyshev 多项式，区间 \([-1, 1]\)，权函数 \(\rho(x) = \sqrt{1-x^2}\)
性质：奇偶性 \(s _n(-x) = (-1)^n s _n(x)\)
Laguerre 拉格朗日多项式，区间 \([0, \infty)\)，权函数 \(\rho(x) = e^{-x}\)
Hermite 埃尔米特多项式，区间 \((-\infty, \infty)\)，权函数 \(\rho(x) = e^{-x^2}\)
定理 函数列线性无关的充要条件 设 \(\{φ _j \} _0^{n}\sub C[a,b]\)，其线性无关，等价于行列式非零：
函数的最佳平方逼近
线性子空间：\(φ _0,φ _1,\dots,φ _n\) 为 \([a,b]\) 上 \(n+1\) 个线性无关函数，用 \(\Phi = \text{Span}\{φ _0,φ _1,\dots,φ _n \}\) 表示由这些函数张成的线性子空间，即任意 \(φ\in \Phi,φ = \sum _{j=0}^{n}a _j φ _j\)
连续形式
已知 \(f(x)\in C[a, b]\)，我们寻求一个 \(φ ^*\) 逼近 \(f\)，使得最小化：
这里 2 范数的平方，其实就是各个点的平方和；权重 \(\rho(x)\) 刻画各个点的重要程度
求解 \(φ ^*\)，即求解系数 \(a _0,\dotsb,a _n\)，等价于使函数 \(F[a _0,\dotsb,a _n] = \int _a^b \rho(x)[\sum a _jφ _j - f(x)]^2 {\rm d}x\) 取极值：
由于已有 \(φ _0,φ _1,\dots,φ _n\) 线性无关，则这就是一个有唯一解的线性方程组，称之为法方程
离散形式 - 最小二乘逼近
现实情况已知的是 \(f(x)\) 的离散点信息 \((x _i, f _i),i=0,1,\dots,m(m>n)\)，问题变成最小化：
和连续情况同理，求 \(F[a _0,\dotsb,a _n]\) 的极值，得到法方程：
此处的内积为离散求和形式：
称得到的 \(φ ^*=\sum a^* _jφ _j\) 为 “最小二乘逼近函数”
模型选择
如何选择 \(φ\) 的具体函数形式？根据数据猜想、根据背景推导、多尝试
此外，对于猜想的非线性函数模型 \(\phi\)，我们往往选择变换成若干函数的线性组合形式 \(φ=F[\phi]=\sum a _jφ _j\)，从而选择 \(\{φ _j \} _0^{n}\) 为基底做逼近
数值积分
考察积分式 \(\int _a^b f(x) {\rm d}x\) 的近似计算
朴素算法为 Newton-Leibnitz 公式：\(\int _a^b f(x) {\rm d}x = F(x)| _a^b, F'(x) = f(x)\)
但是实际应用中会遇到很多问题，比如只知道 \(f\) 的采样数据表；就算知道 \(f\) 表达式、但是原函数 \(F\) 经常不能以有限形式表达；就算 \(F(x)\) 能以有限形式表达、但是过于复杂、反而导致精度不够等等...
数值积分思想：对于线性泛函 \(I(f) = \int _a^b f(x) {\rm d}x\)，设想用另一个线性泛函 \(Q(f)\) 来逼近：
一般的我们只考虑 \(m=0\) 的情形，即
称之为机械求积公式；有点类似用矩形条近似
现在我们只需考虑如何确定 \(H _j,x _j\)
定义 代数精度：称一个机械求积公式具有 \(r\) 次代数精度，若其对于所有次数不超过 \(r\) 的多项式均能 精确成立 \(E(f)=0\)；
这个条件等价于对 \(f(x)=1,x,\dotsb, x ^r\) 都精确成立，也就是：
由定义，反过来求某个求积公式的精度时，依次代入 \(f=1,x,x^2,\dotsb\)，直到某个次数不成立为止
现在观察条件式，它有 \(2n+2\) 个未知量 \(H _j, x _j\)，\(n+1\) 个等式；未知数太多了，不过，我们有定理，从而不用太关注/直接固定 \(x _j\) 的取值：
定理 对于任意给定的 \(n+1\) 个互异节点 \(a\le x _0< x _1 < \dotsb < x _n \le b\)，总存在一组求积系数 \(H _0, H _1, \dotsb, H _n\)，使得求积公式至少有 \(n\) 次精度
证明 考虑 Lagrange 插值多项式：\(L _n(x) = \sum _{j=0}^{n} l _j(x) f(x _j)\)，取 \(H _j = \int _a^b l _j(x) {\rm d} x\)，那么 \(E(f)=I(f)-Q(f) = \int _a^b f(x){\rm d} x - \sum _{j=0}^n \big(\int _a^b l _j(x) {\rm d} x\big) f(x _j)=\int _a^b [f(x)-L _n(x)]{\rm d} x\)，由拉插余项，又等于 \(\frac{1}{(n+1)!}\int _a^b f^{(n+1)}(\xi)P _{n+1}(x){\rm d}x\)
由于 \(f(x)\in M _n\) 最多为 \(n\) 次，\(f^{(n+1)}=0\)，于是 \(E = 0\)；得证
例题
对于给出带未知系数的求积公式，使得代数精度尽量高，求解系数
解法：由代数精度定义式得到系数的方程；注意还需要验证对于更高的精度是否满足，直到不满足为止
等距节点的求积公式 - Newton-Cotes 公式
将区间 \([a,b]\) 划分成 \(n\) 等分，即 \(x _i = a + i \frac{b-a}{n}\)，代入上述定理中的拉式插值系数：\(H _j = \int _a^b l _j(x) {\rm d} x\)（\(l _j(x)\) 为多项式，对任意 \(n\) 可以直接求得其积分），并定义 Cotes 系数 \(C _j = \frac{H _j}{b-a}\)，得（闭型）Newton-Cotes 公式：
此外，由于 \(f(x)\equiv 1\) 时应当精确成立，故 \(\sum C _j = 1\) 恒成立
常见公式：
- \(Q _1(f) = \frac{b-a}{2}[f(a)+f(b)]\) - 梯形公式
用过 \(f _a,f _b\) 的一次函数的面积拟合 - \(Q _2(f) = \frac{b-a}{6}[f(a)+4f(\frac{a+b}{2})+f(b)]\) - Simpson 公式
用过 \(f _a, f _{\frac{a+b}{2}}, f _b\) 的二次函数的面积拟合 - \(Q _4(f)\) - Cotes 公式
- 当 \(n\) 很大时，逼近效果也会受到拉插的 Rouge 现象影响；具体来说，从 \(n=8\) 开始，系数 \(C _j\) 会有正有负，逼近效果变差
开型 Newton-Cotes 公式：对于等距点 \(a=x _0,x _1,\dotsb, x _{n-1}, x _n=b\)，只将 \(x _1,\dotsb, x _{n-1}\) 作为求积节点，但是除了 \(n=2\) 有用，其他不比闭型的好，一般不考虑
\(n=2,\int _a^b f{\rm d}x \approx (b-a)f(\frac{a+b}{2})\) - 中点公式
Newton-Cotes 的数值稳定性
考虑计算过程的机器舍入造成误差，精确值 \(f\)，计算值 \(\tilde{f}\)，舍入误差 \(\epsilon = \tilde{f} - f\)
总的误差为：\(\eta = (b-a)\sum C _j\epsilon _j\)
记 \(\epsilon = \max |\epsilon _j|\)
若 \(C _j\) 均为正（\(n\le 7\) 或 \(n=9\)），\(\sum |C _j| = \sum C _j=1\)，则 \(|\eta|\le (b-a)\epsilon\sum |C _j|=(b-a)\epsilon\)，此时是稳定的
若 \(C _j\) 有正有负时，\(\sum |C _j|>1\) 会很大，无法保证稳定性
Newton-Cotes 的代数精度
之前提到，拉插能使得精度至少为 \(n\)；而再引入等距取点的 Newton-Cotes 公式能使得精度为：
定理 对于给定 \(n+1\) 个互异节点的 Newton-Cotes 公式；\(n\) 为奇数时有 \(n\) 次精度；\(n\) 为偶数则有 \(n+1\) 次精度
某些计算后，得各阶 Newton-Cotes 公式的余项 \(E(f)\)，见书 (4.2.21)
复化 Newton-Cotes 公式
实际计算时，我们无法通过增大次数以提高精度，我们可以将区间 \([a,b]\) 划分成若干子区间（等分，而且通常不断地对半分），子区间内用低阶 Newton-Cotes 求积公式
- 复化梯形公式，假设 \(n=2^k\) 等分，\(h=\frac{b-a}{2^k},x _i = a+ih\)，下同
\(T _n = \frac{h}{2}\sum _{i=0}^{n-1}[f(x _i) + f(x _{i+1})] = \frac{h}{2}[f(a)+f(b) + 2\sum _{i=0}^{n-1}f(x _i) ]\)
若对精度不满意，继续对半分（但是仍为梯形的 2 次精度）：
\(T _{2n}=\frac{h}{4}\sum _{i=0}^{n-1}[f(x _i) +2f(x _{i+1/2})+ f(x _{i+1})] = \frac{h}{4}\sum _{i=0}^{n-1}[f(x _i)+ f(x _{i+1})] + \frac{h}{2}\sum _{i=0}^{n-1}f(x _{i+1/2})=\frac{1}{2}T _n + \frac{h}{2}\sum _{i=0}^{n-1}f(x _{i+1/2})=\frac{1}{2}T _n+\frac{1}{2}U _n\)，可以利用已算好的 \(T _n\) 小小递推一下
余项为 \(E _n = -\frac{b-a}{12}h^2 f''(\eta),\quad \eta\in (a,b)\) - 复化 Simpson 公式：
\(S _n = \frac{h}{6}\sum _{i=0}^{n-1}[f(x _i) + 4f(x _{i+1/2})+f(x _{i+1})] = \frac{h}{6}[f(a)+f(b) + 2\sum _{i=1}^{n-1}f(x _i) + 4\sum _{i=0}^{n-1}f(x _{i+1/2})] = \frac{1}{3}T _n + \frac{2}{3} U _n\)
会发现 \(S _n = \frac{4 T _{2n}- T _n}{4-1}\)，低次精度（梯形的 2 次）的组合能够达到更高阶（抛物线的 3 次）！
Romberg 积分法
回顾通过对低阶精度的组合，复化梯形公式有 \(T _{2n} = \frac{T _n + U _n}{2}\)
现考虑将积分区间分成 \(n = 2^k\) 等分，\(T _{0,k}\) 表示将区间 \(2^k\) 等分后使用复化梯形公式得到的数值积分，递推式就是上一段那个，递推初值为：
求出 \(T _{0,1}\) 后，用如下公式可以求出 Simpson 公式的值：
同理，对区间 \([a,b]\) 进行 \(2^i\) 等分，相应计算出复化梯形公式 \(T _{0,i}\)，于是有：
最后递推可得到 \(T _{i,0}\)，判断其和 \(T _{i-1,0}\) 之差是否在精度误差之内，若满足 \(|T _{i,0}- T _{i-1,0}|\le ε\)，即可退出
Gauss 求积公式
...
常微分方程的数值解法
初值问题概念
一阶常微分方程初值问题，即求解函数 \(y\)：
数值解法，为寻求解 \(y(x)\) 在一系列等距离散点 \(a=x _0\le x _1 \le \dotsb\le x _n\le \dotsb\) 上的近似值 \(y _0=\eta, y _1, \dotsb, y _n, \dotsb\)，且后面的点由前面点的信息递推而来，称此数值解法为 “差分方法”
下文默认步长恒为 \(h\)，即 \(x _n = x _0 + nh\)
下文默认真解为 \(y(x)\)
根据递推式使用的信息，分为 单步法（Euler、Runge-Kutta、Taylor）、多步法（线性多步法）；此外还需要考虑是否收敛、误差累积等问题
离散化方法
就是差分方法的推导思路
几乎所有差分方法可以通过三种离散化方法之一得到：Taylor 展开、化导数为差商、数值积分方法。
Taylor 展开方法
在 \(x _n\) 点展开：\(y(x _{n+1}) = y(x _n) + hy'(x _n) + \frac{h^2}{2!}y''(\xi _n)\)
\(h\) 足够小时，略去余项（需要注意余项可能累积），得到：
\(y(x _{n+1}) \approx y(x _n) + hy'(x _n) = y(x _n) + hf(x _n, y(x _n))\)
利用其，可以得到递推式，称为 Euler（折线）法：
化导数为差商
对于充分小的 \(h\)：\(\frac{y(x _{n+1}) -y(x _n)}{h}\approx y'(x _n) = f(x _n, y(x _n))\)
得到 \(y(x _{n+1}) \approx y(x _n) + hf(x _n, y(x _n))\)，也能推出 Euler 法
数值积分方法
在 \([x _n, x _{n+1}]\) 上对 \(y'(x) = f(x, y(x))\) 积分：\(y(x _{n+1}) - y(x _{n}) = \int _{x _n}^{x _{n+1}}f(x, y(x)){\rm d}x\approx hf(x _n, y(x _n))\)，若 \(h\) 足够小
也能推出 Euler 法
Euler 折线法
推导见上，递推式为
几何意义显然：每次沿着切线方向飞出去
但是误差积累严重（尤其一阶导不变号时，总是积累同一个符号的误差）
线性多步法
思想：通过之前若干节点和其一阶导来逼近，形式化地说：
若 \(b _{-1} \ne 0\)，则右侧含有 \(y _{n+1}\)，此为“隐式方法”，需要解方程，但是更精确更多采用；反之称“显式方法”
线性多步法的逼近准则
准确成立：称上述递推式对 \(y(x)\) 准确成立，若：
若对于任意 \(y(x)\in {\bf M} ^r\)，均 准确成立；而当 \(y\) 为 \(r+1\) 次多项式时则 不 准确成立，则称此线性多步法是 \(r\) 阶的
例如，Euler 法是一阶的（只对直线恒准确成立）
局部截断误差
引入误差概念，对于某一步 \(n\to n+1\)，计算 \(y(x _{n+1})\) 和代入 \(y (x _{n-i}), y' (x _{n-i})\)（而不是自己算出的 \(y _{n-i}, y' _{n-i}\)，即假设之前每一步计算都是准确的，为了使问题研究简单，这是模型误差）的递推式的差 \(T _n\)，称为局部截断误差
我们再将所有 \(y (x _{n-i}), y' (x _{n-i})\) 在 \(x _n\) 点泰勒展开，整理得：
其中 \(C _q\) 都是关于 \(a _i, b _i\) 的式子，具体形式不表
回到准确成立，也就是局部误差 \(T _n = 0\)，有定理：
定理 \(r\) 阶准确成立的充要条件
\(r\) 阶准确成立，充要条件为 \(C _0 = C _1 = \dotsb = C _r = 0, C _{r+1}\ne 0\)
理解：对于阶数小等于 \(r\) 的多项式 \(y(x)\)，\(q> r\) 时有 \(y^{(q)}(x)\equiv 0\)；而对于 \(q\le r\) 部分有 \(C _q = 0\)，故 \(T _n = 0\) 恒成立；对于 \(r+1\) 阶多项式，\(T _n = C _{r+1}h^{r+1}y^{(r+1)}(x _n)\ne 0\)，称此 \(C _{r+1}\) 为误差常数项
定义 相容：称满足条件 \(C _0 = C _1 = 0\) 的线性多步法是相容的
求解递推方程
递推方程有 \(2p+3\) 个未知数；若要求精度为 \(r\)，则 \(y(x)\in {\bf M}^r\) 由基 \(1,x,x^2,\dotsb, x^r\) 构成，即需要对这 \(r+1\) 个基均准确成立，故有 \(r+1\) 个方程，可见：\(r+1\le 2p+3\implies r\le 2p+2\)
接下来解方程：\(C _0 = C _1 = \dotsb = C _r = 0\) 即可
for exam：设定 \(p\) 和某些系数，求解精度最高阶方法
例 1，\(p=0\) 时，使精度达到最高阶的方法
解，\(r\le 2p+2 = 2\)，\(C _0 = C _1 = C _2 = 0\)，解得 \(a _0, b _{-1,0}\)，为一步、二阶、隐式法：
\(y _{n+1} = y _n + \frac{h}{2} (y' _{n+1} + y' _{n})\)，误差常数 \(C _3 = -\frac{1}{12}\)
称之为梯形方法例 2，\(p=1\) 时，以 \(a _1\) 为自由参数，使精度达到最高阶方法
解，\(r\le 2p+2-1=3\)，\(C _0 = C _1 = C _2 = C _3 = 0\)，解得 \(a _0, b _{-1,0,1}\)，为两步、至少三阶、隐式法例 3，例 2 中 \(a _1\) 为何值时达到最高阶
解，其实就是让 \(r=2p+2=4\)，引入 \(C _4=0\) 即可；顺便判断下 \(C _5\ne 0\)
得 \(y _{n+1} = y _{n-1} + \frac{h}{3} (y' _{n+1} + 4y' _{n} + y' _{n-1})\)
称之为 Simpson 方法
线性多步法-收敛性
回顾：微分方程 \(y(a) = \eta,\frac{{\rm d}y}{{\rm d}x} = f(x, y),\quad a\le x\le b\)
回顾：线性多步法 \(y _{n+1} = \sum _{i=0}^{p} a _i y _{n-i} + h \sum _{i=-1}^{p} b _i y' _{n-i},\quad n \ge p\)
定义 收敛性
称一个线性多步法是收敛的，若线性多步法初始条件 \(y _k = y _k(h),k = 0,1,\dotsb,p\)，满足 \(\lim _{h\to 0} y _k(h) = \eta\)（\(h\to 0\) 且 \(p\) 有限，其实就是线性递推的初值为 \(\eta\) 的表述）
且若递推方程的解 \(y _{n}, n\to \infty\) 对任意固定的 \(x\in [a, b]\)，都满足 \(\lim _{h\to 0, n\to \infty} y _n = y(x), \color{deeppink}{nh = x - a}\)，则称此线性多步法是收敛的
试验方程
为了判断某个线性多步法是否收敛，引入
称此为试验方程，实际上解为 \({\color{royalblue}{y(x)}} = \eta e ^{\lambda(x-a)}\)；
而将此方程用线性多步法求解时，设 \(y _n = r ^n\)，会得到 \((1-h\lambda b _{-1})r ^{p+1} = \sum _{i=0}^{p}(a _i + h\lambda b _i)r ^{p-i}\) 称为此线性多步法的特征方程，它具有 \(p+1\) 个根 \(r _0(h\lambda),r _1(h\lambda),\dotsb,r _p(h\lambda)\)，而通解为 \({\color{deeppink}{y _n}}= \sum _{i=0}^{p} d _i[r _i (h\lambda)]^n\)，其中 \(d _i\) 为任意常数
验证一个线性递推式是否收敛，我们可以研究此通解 \(\color{deeppink}{y _n}\) 和实际解 \(\color{royalblue}{y(x)}\) 是否满足收敛的定义
Runge-Kutta 法
考虑单步法——泰勒展开的高阶形式 \(y _{n+1} = y _n + h y' _n + \frac{h^2}{2!} y'' _n +\dotsb + \frac{h^r}{r!}y ^{(r)} _n\)
其中 \(y' = f\)，那么对于高阶导 \(y^{(l)}\) 或者说 \(f^{(l-1)}\)，如果将其一并以 \(f\) 表示呢？
考虑引入某个 \((x+a,y+b)\) 处对二元函数 \(f\) 的泰勒展开：\(f(x+a,y+b) \approx f(x, y) + (a\frac{\partial}{\partial x} + b\frac{\partial}{\partial y})f + \frac{1}{2!}(a\frac{\partial}{\partial x} + b\frac{\partial}{\partial y}) ^2 f + \dotsb + \frac{1}{r!}(a\frac{\partial}{\partial x} + b\frac{\partial}{\partial y})^r f + \dotsb\)
可见，我们可以以 \(a,b\) 为待解的参数，表达对应的高阶导数；这大概就是 RK 法递推形式的来源
引入 RK 法，一般形式为：
其中 \(\alpha, b, c\) 均为常数，\(c _1 = 0, \alpha _{1j} = 0, j=1,2,\dotsb, s-1\)
希望这个递推式准确成立，引入局部截断误差，定义为 \(T _n = y(x _{n+1}) - y _{n+1}\)；如果希望 RK 法能精度足够高，也就是希望能达到某个 \(r\) 阶的准确成立，那就是让对应截断误差为零
二级 RK 法
考虑 \(s = 2\)
\(K _1=f\) 显然含义为近似一阶导
对 \(K _2\) 展开：
代回递推式，得 \({\color{deeppink}{y _{n+1}}} = y _n + h(b _1 + b _2)f + h^2 b _2 (c _2 f' _x + \alpha _{21} f f' _y) +\dotsb + O(h ^4)\)
而对于 \(y(x _{n+1})\) 在 \(x _n\) 展开，并将 \(y ^{(l)}\) 用 \(f^{(l-1)}\) 表示，得 \({\color{deeppink}{y(x _{n+1})}} = y(x _n) + hf + \dotsb + O(h^4)\)
回到截断误差（假定 \(y _n = y (x _n)\)）有
希望准确成立——截断误差为零，解方程：\(b _1 + b _2 = 1, b _2 c _2 = \frac{1}{2}, \alpha _{21}b _2 = \frac{1}{2}\)，而 \(h^3\) 系数里带常数无法恒零，故只解此方程，最后二级 RK 法截断误差：\(T _n = O(h^3)\)，最多达到 2 阶准确成立
上述方程有无穷解，考虑 \(c _2\) 为参数（含义为对 \(K _2 = f(x _n + c _2 h, y _n + \dotsb)\) 中，\(x _n\le x _n + c _2 h \le x _{n+1}\) 的位置选择），求解其余参数。
其中 \(c _2 = \frac{2}{3}\) 时误差取最小，称之为 Heun 法；此外还有 \(c _2 = \frac{1}{2}\) 中点方法、\(c _2=1\) 改进 Euler 法。
四级 RK 法 - Runge-Kutta 法
虽然可以增大极数以增大准确成立的阶数，但是级数 \(s\) 增大时，阶数增长速度减慢（并不是跟着线性增长的）；应用中最常用的是四级 RK 法，也就是 Runge-Kutta 法；经典形式为
