本篇为应付考试而整理的复习笔记，因而只是整理我个人认为考前需要速通的点，因而不一定全面以及完备。
误差分析
两类基本误差
绝对误差：\(e_{\bar{x}}=x-\bar{x}\)
相对误差：\(r_{\bar{x}}=\frac{e_{\bar{x}}}x\)
由于 \(x\) 往往未知，在 \(r_{\bar{x}}\) 很小的时候，可以用 \(r^*_{\bar{x}}=\frac{e_{\bar{x}}}{\bar{x}}\) 近似
舍入误差与有效数字
记十进制数 \(a\) 某个数码 \(a_s\) 的位数为 \(k(a_s)\)，当 \(a_s\) 在个位时取 0，十位时取 1，小数点后第一位时取 -1，以此类推
舍入误差：将 a 四舍五入至 \(k\) 位小数，那么 \(|e_{\bar{a}}|\leq\frac12\times10^{-k}\)
有效数字：设 \(\bar{a}\) 为 \(a\) 的近似值，从 \(\bar{a}\) 中首位不为零的数码起，称数码 \(a_s\) 为有效数字，若满足： \[ |e_{\bar{a}}|\leq\frac12\times10^{k(a_s)} \] （因而对于四舍五入得到的近似值，从首位不为零的数码起，所有数码都是有效数字）
利用有效数字，我们有如下相对误差估计：
设 \(\bar{a}\) 有 \(k\) 位有效数字，其中第一位有效数字为 \(a_s\)，那么： \[ |r^*_{\bar{a}}|\leq\frac1{2a_s}\times10^{-(k-1)} \] 反之，若有上述不等式成立，则 \(\bar{a}\) 至少有 \(k\) 位有效数字
（这说明有效数字越多，相对误差越小）
运算中的误差传播
考察二元函数 \(z=f(x,y)\) 产生的误差：\(e_{\bar{z}}=f(x,y)-f(\bar{x},\bar{y})\)
当 \(|e_{\bar{x}}|,|e_{\bar{y}}|\) 很小，且 \(f\) 可微时，有近似公式： \[ e_{\bar{z}}\approx \frac{\partial f}{\partial x}(\bar{x},\bar{y})\:e_{\bar{x}}+\frac{\partial f}{\partial y}(\bar{x},\bar{y})\:e_{\bar{y}}\\ r_{\bar{z}}\approx \frac{\bar{x}}{\bar{z}}\frac{\partial f}{\partial x}(\bar{x},\bar{y})\:r_{\bar{x}}+\frac{\bar{y}}{\bar{z}}\frac{\partial f}{\partial y}(\bar{x},\bar{y})\:r_{\bar{y}} \] （对于 \(r^*\) 亦有类似的近似）
由此导出：
\(f(x,y)=x\pm y\) \[ e_{\bar{x}\pm \bar{y}}=e_{\bar{x}}\pm e_{\bar{y}}\\ r_{\bar{x}\pm \bar{y}}=\frac{\bar{x}}{\bar{x}\pm \bar{y}}r_{\bar{x}}\pm\frac{\bar{y}}{\bar{x}\pm \bar{y}}r_{\bar{y}} \] 可见 \(x\approx y\) 时作减法会产生较大的相对误差，这被称为相减相消，应尽量避免
\(f(x,y)=xy\) \[ e_{\bar{x}\bar{y}}=\bar{y}e_{\bar{x}}+\bar{x}e_{\bar{y}}\\ r_{\bar{x}\bar{y}}=r_{\bar{x}}+r_{\bar{y}} \]
\(f(x,y)=x/y\) \[ e_{\bar{x}/\bar{y}}=\frac{\bar{y}e_{\bar{x}}-\bar{x}e_{\bar{y}}}{\bar{y}^2}\\ r_{\bar{x}/\bar{y}}=r_{\bar{x}}-r_{\bar{y}} \] 可见 \(y\to0\) 时会产生较大的绝对误差，应尽量避免
计算机浮点系统
称 \(F(p,t,L,U)\) 为一个计算机浮点系统，其中：
- \(p\)：基数（进制）
- \(t\)：精度（最大小数位数）
- \(L,U\)：指数的取值范围 \([-L,U]\)（整数）
属于浮点系统 \(F\) 的数 \(x\) 具有如下形式： \[ x=\pm p^J\sum_{k=1}^td_kp^{-k}=\pm (0.d_1d_2...d_t)_p\times p^J\\ J\in \mathbb{Z}\cap[-L,U],\:d_k\in\{0,1,2,...,p-1\} \] 称 \((0.d_1d_2...d_t)_p\) 为 \(x\) 的尾数，\(J\) 为 \(x\) 的指数（阶）
称 \(x\) 为规格化浮点数，若 \(d_1\neq0\)
UFL 与 OFL
- 称 \(F\) 中的最小正数为 \(UFL=p^{-t}\times p^{-L}=p^{-(L+t)}\)
- 称 \(F\) 中的最大正数为 \(OFL=(1-p^{-t})\times p^U\)
机器精度
\(F\) 自然不能表示任意实数，称 \(F\) 中用于表示一个在 UFL 以及 OFL 范围内的实数 \(x\) 的最大的相对误差为机器精度 \(\varepsilon_{mach}\)
- 若采取截断法：\(\varepsilon_{mach}=p^{1-t}\)
- 若采取最近舍入法（四舍六入五取双）：\(\varepsilon_{mach}=\frac12p^{1-t}\)
一般地：\(0<UFL<\varepsilon_{mach}<OFL\)
机器算术运算
以 \(\oplus\:\ominus\:\otimes\) 等符号代表机器算术运算符号， 可以假定机器运算的定义为： \[ x\odot y=fl(fl(x)\cdot fl(y)) \] 其中 \(\cdot\) 泛指上述符号，\(fl(x)\) 为 \(F\) 中最接近 \(x\) 的数，也就是 \(x\) 在 \(F\) 中的表示
乘除运算是直接进行的，只需要分别处理底数部分和指数部分即可；
而对于加减运算，则需要对阶操作，即首先将阶数较小的数调整为与阶数较大的数同阶数（相应地，底数修改小数点位置），这个过程中若底数超出了精度 \(t\)，则需要截断或最近舍入；完成加减运算之后，仍然需要作截断或最近舍入
特别地，对于减法，我们有精度丢失定理：设 \(x>y\) 均为正的规格化二进制浮点数，且 \(2^{-q}\leq1-\frac yx\leq2^{-p}\) 则 \(x\ominus y\) 丢失至多 \(q\) 个，至少 \(p\) 个有效二进制位
数值稳定与病态问题
称一个算法是数值稳定的，若在计算过程中，误差不会增长，否则称为数值不稳定的
例：在计算递推式时，式中出现了形如 \(x_{n+1}=ax_n,\:a>1\) 的情况，那么误差就会被 \(a\) 不断放大，导致数值不稳定
称一个问题是病态的，若输入的微小变化能够导致输出的巨大误差
例：求解线性方程组时系数矩阵接近奇异，如 \(Hilbert\) 矩阵
插值法
插值问题是指给定一组离散数据 \(\{x_k,y_k\},k=0,1,...,n\)，我们试图寻求一函数形如： \[ \varphi(x)=a_0\varphi_0(x)+a_1\varphi_1(x)+...+a_n\varphi_n(x) \] 使得 \(\varphi(x_k)=y_k,\forall k=0,1,...,n\)
\(Lagrange\) 插值
记 \(l_i(x)=\prod_{j=0,j\neq i}^n\frac{x-x_j}{x_i-x_j}\)，则 \(l_i(x_j)=\delta_{ij}\)，称为 \(Lagrange\) 基本多项式
则 \(Lagrange\) 插值多项式为：\(L_n(x)=\sum_{i=0}^nf(x_i)l_i(x)\)，容易证明这是存在唯一的
方便起见，若记 \(\omega_{n+1}(x)=\prod_{i=0}^n(x-x_i)\)，那么 \(\omega_{n+1}'(x_i)=\prod_{j=0,j\neq i}^n(x_i-x_j)\)
于是 \(l_i(x)=\frac{\omega_{n+1}(x)}{(x-x_i)\omega_{n+1}'(x_i)}\)
余项
设离散数据 \(\{x_k,y_k\},k=0,1,...,n\) 是由函数 \(f\) 给出，也就是 \(y_i=f(x_i)\)，称： \[ r_n(x)=f(x)-L_n(x) \] 为 \(Lagrange\) 插值多项式的余项
假设 \(f\) 在包含 \(\{x_n\}\) 的区间 \([a,b]\) 上具有 \(n\) 阶连续导数，且在 \((a,b)\) 内存在 \(n+1\) 阶有界导数，
那么 \(\forall x \in [a,b],\:\exists \xi\in(a,b)\) 使得： \[
r_n(x)=\frac{f^{n+1}(\xi)}{(n+1)!}\omega_{n+1}(x)
\] 其证明需要简单了解，大致如下：
因为 \(x_0,x_1,...,x_n\) 都是 \(r_n(x)\) 的一重零点，故可以设 \(r_n(x)=\omega_{n+1}(x)K(x)\)，为确定 \(K(x)\)，引入辅助函数： \[
F(t)=f(t)-L_n(t)-K(x)\omega_{n+1}(t)
\] 则 \(t = x_0,x_1,...,x_n,x\) 时均有 \(F(t)=0\)，反复应用 \(Rolle\) 定理就有 \(0=F^{(n+1)}(\xi)=f^{(n+1)}(\xi)-K(x)(n+1)!\)
从而不难解出 \(K(x)\) （这说明 \(\xi\) 与 \(x_0,x_1,...,x_n,x\) 有关）
通过考察余项，我们可以估计误差，特别地，对于线性插值，有： \[ |r_1(x)|\leq\frac{(b-a)^2}{8}|f''(\xi)| \] 对于二次插值，有： \[ |r_2(x)|\leq \frac{\sqrt3}{27}(b-a)^3|f'''(\xi)| \]
\(Runge\) 现象
一般而言，\(n\) 越大，\(L_n(x)\) 在端点附近的抖动也越大，这被称为 \(Runge\) 现象，因而高次插值法一般不能用于外插
\(Newton\) 插值
均差
称 \(f[x_i,x_j]=\frac{f(x_j)-f(x_i)}{x_j-x_i}\) 为 \(f\) 关于基点 \(x_i,x_j\) 的一阶均差
称 \(f[x_i,x_j,x_k]=\frac{f[x_j,x_k]-f[x_i,x_j]}{x_k-x_i}\) 为 \(f\) 关于基点 \(x_i,x_j,x_k\) 的二阶均差
一般地，\(f[x_0,x_1,...,x_n]=\frac{f[x_1,x_2,...,x_n]-f[x_0,x_1,...,x_{n-1}]}{x_n-x_0}\)
方便起见，不妨规定 \(f[x_i]=f(x_i)\)
用归纳法可以证明： \[ f[x_0,x_1,...,x_n]=\sum_{i=0}^n\frac{f(x_i)}{\prod_{j=0,j\neq i}^n(x_i-x_j)} \]
差商表
常用差商表来计算各个均差：
\[ \left[ \begin{array}{} x_0 & f(x_0)\\ x_1 & f(x_1) & f[x_0,x_1]\\ x_2 & f(x_2) & f[x_1,x_2] & f[x_0,x_1,x_2]\\ ... & ... & ... & ... & ...\\ x_n & f(x_n) & f[x_{n-1},x_n] & f[x_{n-2},x_{n-1},x_n] & ... & f[x_0,...,x_n]\\ \end{array} \right] \]
\(Newton\) 插值多项式
\(Lagrange\) 插值的缺陷是每新增一个样本点都要全部重新算一遍，\(Newton\) 插值优化了这一点，其关键是将多项式考虑为如下形式： \[ N_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-1}) \] 从而不难看出 \(a_k=f[x_0,x_1,...,x_k]\)
若新增一插值基点 \(\bar{x}\)，那么 \(f(\bar{x})-N_n(\bar{x})=N_{n+1}(\bar{x})-N_n(\bar{x})=f[x_0,x_1,...,x_n,\bar{x}]\omega_{n+1}(\bar{x})\)
据插值多项式的唯一性，比较上一节中的余项就有： \[
f[x_0,x_1,...,x_n,x]=\frac{f^{(n+1)}(\xi)}{(n+1)!}
\]
有限差
下面我们考虑样本点满足 \(x_k-x_{k-1}=h,\forall k=1,2,..,n\) 的简单情形，方便起见，引入有限差的概念，分别称： \[
\Delta f(x)=f(x+h)-f(x)\\
\nabla f(x)=f(x)-f(x-h)\\
\delta f(x)=f(x+\frac h2)-f(x-\frac h2)
\] 为函数 \(f(x)\) 在点 \(x\) 的一阶向前差分，一阶向后差分以及一阶中心差分，易见它们都是线性算子
下面以向前差分为例定义高阶情形： \[
\Delta^2f(x)=\Delta(\Delta f(x))=\Delta f(x+h)-\Delta f(x)\\
\Delta^nf(x)=\Delta(\Delta^{n-1} f(x))=\Delta^{n-1} f(x+h)-\Delta^{n-1} f(x)\\
\] 类似地，方便起见我们规定 \(\Delta^0f(x)=f(x)\)，用归纳法可以证明： \[
\Delta^n f(x)=\sum_{i=0}^n(-1)^iC_n^if(x+(n-i)h)\\
\nabla^n f(x)=\sum_{i=0}^n(-1)^iC_n^if(x-ih)\\
\delta^n f(x)=\sum_{i=0}^n(-1)^iC_n^if(x+(\frac n2-i)h)
\] 反过来，我们也有： \[
f(x+nh)=\sum_{i=0}C_n^i\Delta^if(x)\\
f(x-nh)=\sum_{i=0}^n(-1)^iC_n^i\nabla^if(x)\\
f(x+nh)=\sum_{i=0}^nC_n^i\delta^if(x+\frac i2h)
\] 最后，对于满足 \(x_k-x_{k-1}=h,\forall k=1,2,..,n\) 的样本点，我们有： \[
f[x_0,x_1,...,x_n]=\frac{\Delta^nf(x_0)}{n!h^n}=\frac{\nabla^nf(x_n)}{n!h^n}
\] 对于中心差分则需要分样本点的奇偶个数作讨论，原理类似
差分表
可以用差商表类似的方法构造差分表以计算差分，以向前差分为例：
\[ \left[ \begin{array}{} x_i=x_0+ih & \Delta & \Delta^2 & \Delta^3 & ... & \Delta^{n-1} & \Delta^n\\ x_0 & \Delta f(x_0) & \Delta^2 f(x_0) & \Delta^3 f(x_0) & ... &\Delta^{n-1} f(x_0)& \Delta^n f(x_0)\\ x_1 & \Delta f(x_1) & \Delta^2 f(x_1) & \Delta^3 f(x_1) & ... &\Delta^{n-1} f(x_1)&\\ x_2 & \Delta f(x_2) & \Delta^2 f(x_2) & \Delta^3 f(x_2) & ...\\ ... & ... & ... & ... & ... \\ x_{n-1} & \Delta f(x_{n-1}) & \Delta^2 f(x_{n-1}) \\ x_n & \Delta f(x_n) \end{array} \right] \]
\(Newton\) 前插和后插公式
考虑样本点满足 \(x_k-x_{k-1}=h,\forall k=1,2,..,n\) 的情形，利用均差和向前差分的关系，作变换：\(x=x_0+sh\)
从而 \(Newton\) 插值公式可以整理为： \[
N_n(x)=N_n(x_0+sh)=\sum_{k=0}^n\left(\begin{array}{}s\\k\end{array}\right)\Delta^kf(x_0)\\
\left(\begin{array}{}s\\k\end{array}\right)=\frac{s(s-1)...(s-k+1)}{k!}
\] 即 \(Newton\) 前插公式；利用向后差分，作变换：\(x=x_n+th\)，同理有： \[
N_n(x)=N_n(x_n+th)=\sum_{k=0}^n\left(\begin{array}{}-t\\k\end{array}\right)(-1)^k\nabla^kf(x_n)\\
\left(\begin{array}{}-t\\k\end{array}\right)(-1)^k=\frac{t(t+1)...(t+k-1)))}{k!}
\] 即 \(Newton\) 后插公式
Hermite 插值
在 \(Hermite\) 插值中我们不仅已知样本点的值，还知道样本点的前 \(n\) 阶导数值，因而我们希望寻求一插值多项式使得其导数值也与我们的已知数据符合，下面讨论 \(n=1\) 的简单情形：
设已知 \(f(x)\) 在插值基点 \(x_0,x_1,...,x_n\) 处的函数值分别为 \(y_0,y_1,...,y_n\) 以及一阶导数值分别为 \(y'_0,y'_1,...,y'_n\)
我们试求一多项式 \(H_{2n+1}\in R[x]_{2n+2}\) 使得： \[
\left\{ \begin{aligned}
H_{2n+1}(x_i)=y_i\\
H'_{2n+1}(x_i)=y'_i
\end{aligned}\right.
\] 不妨设： \[
H_{2n+1}=\sum_{i=0}^ny_iA_i(x)+\sum_{i=0}^ny'_iB_i(x)\\
A_i(x_j)=\delta_{ij},\:A'_i(x_j)=0\\
B_i(x_j)=0,\:B'_i(x_j)=\delta_{ij}
\] 可以通过分析 \(A_i,\:B_i\) 各自零点的重数设出其形式进而求解，最后可以得到： \[
B_i(x)=(x-x_i)l_i^2(x)\\
A_i(x)=\left(1-2(x-x_i)\sum_{j=0,j\neq i}^n\frac1{x_i-x_j}\right)l_i^2(x)
\] 与 \(Lagrange\) 插值类似，可以证明 \(Hermite\) 插值的唯一性，以及其余项可表示为： \[
r(x)=\frac{f^{(2n+2)}(\xi)}{(2n+2)!}\omega_{n+1}^2(x)
\]
分段插值
为了解决高次多项式插值带来的 \(Runge\) 现象，可以使用分段插值，也就是将插值区间细分为若干子区间，并在这些子区间内使用较低次数的多项式进行插值的办法，其中最简单的分段插值就是分段线性插值
分段插值误差界
直观上，采用分段插值时，子区间分的越细，那么分段插值的误差界应该越小
下面我们考虑将插值区间分为等距基点 \(x_k-x_{k-1}=h,\forall k=1,2,..,n\) 的情形，由之前对余项的讨论，有：
分段线性插值： \[ |e_{f_1(x)}|\leq\frac{1}8Mh^2\\ M=\mathop{max}\limits_{x\in[a,b]}|f''(x)| \]
分段抛物线插值（即每一个区间含三个基点采用二次多项式插值）： \[ |e_{f_2(x)}|\leq\frac{\sqrt{3}}{27} M h^3\\ M=\mathop{max}\limits_{x\in[a,b]}|f'''(x)| \]
三次样条插值
直接作分段插值的一个严重的问题是得到的结果往往不光滑，因而我们可以考虑在每一区间 \([x_i,x_{i+1}]\) 上采用三次多项式插值，并要求每一段的拼接点的导数以及二阶导数相等，这种插值办法被称为三次样条插值
直观上，我们需要求解 \(4n\) 个待定系数，然而假设的条件只能给出 \(4n-2\) 个方程，因而我们需要附加两个条件，常采用边界条件：
- 自然边界条件：取插值函数在 \(x_0, x_n\) 处的二阶导数为 \(0\)，这样得到的插值称为自然三次样条插值
- 完备边界条件：假设我们已知待插值函数 \(f\) 在端点处的导数值，取插值函数在 \(x_0, x_n\) 处的导数与之相等即可，这样得到的插值称为完备三次样条插值
下面简要讨论具体的计算，给定 \([a,b]\) 的一组基点： \[ a=x_1<x_2<...<x_{n+1}=b \] 记 \([x_i,x_i+1]\) 上的三次插值多项式为 \(S_i(x),\:i=1,2,...,n\)，已知条件为： \[ S_{i-1}(x_i)=S_i(x_i)=f(x_i),\:i=2,...,n\\ S'_{i-1}(x_i)=S'_i(x_i),\:i=2,...,n\\ S''_{i-1}(x_i)=S''_i(x_i),\:i=2,...,n\\ S_1(x_1)=f(x_1),\:S_n(x_{n+1})=f(x_{n+1}) \] 记 \(S''_i(x_i)=m_i,\:h_i=x_{i+1}-x_i\)，因为 \(deg\:S''_i=1\)，结合 \(S''_{i-1}(x_i)=S''_i(x_i)\) 就知道其必形如： \[ S''_i(x)=m_i\frac{x_{i+1}-x}{h_i}+m_{i+1}\frac{x-x_i}{h_i} \] 对其积分得： \[ S'_i(x)=-m_i\frac{(x_{i+1}-x)^2}{2h_i}+m_{i+1}\frac{(x-x_i)^2}{2h_i}+A_i\\ S_i(x)=m_i\frac{(x_{i+1}-x)^3}{6h_i}+m_{i+1}\frac{(x-x_i)^3}{6h_i}+A_i(x-x_i)+B_i \] 结合 \(S_i(x_i)=f(x_i),\:S_i(x_{i+1})=f(x_{i+1})\) 可以解得： \[ B_i=f(x_i)-m_i\frac{h_i^2}6\\ A_i=f[x_i,x_{i+1}]-\frac{h_i}6(m_{i+1}-m_i) \] 再由 \(S'_{i-1}(x_i)=S'_i(x_i)\) 就得到： \[ h_{i-1}m_{i-1}+2(h_{i-1}+h_i)m_i+h_im_{i+1}=6(f[x_i,x_{i+1}]-f[x_{i-1},x_i]) \] 也就是： \[ \left(\begin{array}{} h_1 &2(h_1+h_2) & h_2 & 0 & ... & 0 \\ 0 & h_2 &2(h_2+h_3) & h_3 & ... & 0 \\ ... & ... & ... & ... & ... & ... \\ 0 & ... & h_{n-2} &2(h_{n-2}+h_{n-1}) & h_{n-1} & 0\\ 0 & 0 & ... & h_{n-1} &2(h_{n-1}+h_{n}) & h_{n} \end{array}\right) \left(\begin{array}{} m_1\\ m_2\\ ...\\ m_n\\ m_{n+1} \end{array}\right)=6 \left(\begin{array}{} d_2\\ d_3\\ ...\\ d_{n-1}\\ d_n \end{array}\right)\\\\ d_i=f[x_i,x_{i+1}]-f[x_{i-1},x_i] \] 前已述该方程没有唯一解，需要补上一定的边界条件，即：
自然边界条件：\(m_1=m_{n+1}=0\)
完备边界条件：\(S'_1(a)=f'(a),\:S'_n(b)=f'(b)\)，进而： \[ 2h_1m_1+h_1m_2=6(f[x_1,x_2]-f'(a))\\ h_nm_n+2h_nm_{n+1}=6(f'(b)-f[x_n,x_{n+1}]) \]
补上这些边界条件即有唯一解
基样条插值
注意到 \(Lagrange\) 插值多项式总可以写成 \(Lagrange\) 基本多项式 \(l_i(x)\) 的线性组合，亦可寻找一基三次样条系来表示三次样条插值函数
常见的办法是构造一在区间 \([-2,2]\) 上的分段三次多项式 \(\varphi(x)\)，使得其在 \(x=-1,0,1\) 上面的取值以及导数值，二阶导数值较好，以及在 \(x=-2,2\) 上的函数值，导数值，二阶导数值均为 \(0\)，从而可以通过将其根据实际的基点位置作平移和伸缩变换得到一系列基函数
函数逼近
正交函数系
称 \(W(x)\) 为 \([a,b]\) 上的权函数，若满足：
- 在 \((a,b)\) 中，\(W(x)\geq 0\)，且只有有限个零点
- \(\int_a^bx^iW(x)dx,i=1,2,...\) 存在且有限
则对于非负连续函数 \(f\)，若 \(\int_a^bf(x)W(x)dx=0\)，则在 \([a,b]\) 上有 \(f(x)\equiv 0\)
给定 \(f,g\in C[a,b]\)，称 \((f,g)=\int_a^bW(x)f(x)g(x)dx\) 为函数 \(f,g\) 在 \([a,b]\) 上关于权 \(W(x)\) 的内积
（满足内积的各种性质，亦可以由此定义范数）
我们常用下面三类范数： \[ \Vert f\Vert_\infty=\mathop{max}_\limits{x\in[a,b]}|f(x)W(x)|\\ \Vert f\Vert_1=\int_a^b|f(x)W(x)|dx\\ \Vert f\Vert_2=\sqrt{(f,f)} \]
称函数系 \(\{f_n(x)\}\) 为 \([a,b]\) 上关于权函数 \(W(x)\) 的正交函数系，若 \((f_i,f_j)=A_i\delta_{ij}\)
线性无关
可以与线性代数理论类似考虑一族函数的线性无关性，称函数系 \(\{f_n(x)\}\) 在 \([a,b]\) 上线性无关，若： \[ \sum_{i=0}^na_if_i=0 \Rightarrow a_i=0,\forall i=0,1,...,n \] 可以证明这等价于： \[ \left| \begin{array}{} (f_0,f_0) & (f_0,f_1) &... &(f_0,f_n)\\ (f_1,f_0) & (f_1,f_1) &... &(f_1,f_n)\\ ... & ... & ... & ... \\ (f_n,f_0) & (f_n,f_1) &... &(f_n,f_n) \end{array} \right|\neq0 \] 该行列式被称为 \(Gram\) 行列式
可以用 \(Schmidt\) 正交化将一族线性无关函数化为正交函数系，即： \[ \varphi_0=f_0\\ \varphi_n=f_n-\sum_{i=0}^{n-1}\frac{(f_n,f_i)}{(f_i,f_i)}f_i \]
直交多项式
取线性无关的多项式序列 \(\{1,x,...,x^n\}\)，称由此通过 \(Schmidt\) 正交化得到的正交函数系为直交多项式 \(\{p_n(x)\}\)
可以证明有递推关系： \[
\begin{aligned}
&p_{k+1}(x)=(x-\alpha_k)p_k(x)-\beta_kp_{k-1}(x)\\
&p_0=1,\:p_{-1}=0\\
&其中:\\
&\alpha_k=\frac{(xp_k(x),p_k(x))}{(p_k(x),p_k(x))}\\
&\beta_k=\left\{
\begin{align}
&\frac{(p_k(x),p_k(x))}{(p_{k-1}(x),p_{k-1}(x))},k\geq1\\
&0,k=0
\end{align}
\right.\\
\end{aligned}
\]
\(Chebyshev\) 多项式
取权函数 \(W(x)=\frac1{\sqrt{1-x^2}}\)，区间 \([a,b]=[-1,1]\)，由此得到 \(Chebyshev\) 多项式，可表示为： \[ T_n(x)=cos(n\:arccos\:x) \] 显然，取 \(x=cos\theta\)，据三角公式就有： \[ T_{n\pm1}(x)=cos(n\pm1)\theta=cos(n\theta)cos\theta\mp sin(n\theta)sin\theta \] 两式相加就能得到递推关系： \[ T_{n+1}(x)+T_{n-1}(x)=2xT_n(x) \] 利用 \(Euler\) 公式，有：\(cos\:n\theta=\frac{e^{in\theta}+e^{-in\theta}}2\)，由此亦可直接证明： \[ T_n(x)=\frac12\left((x+\sqrt{x^2-1})^n+(x-\sqrt{x^2-1})^n\right) \] #### 近似最佳一致逼近
固定 \(z>1\)，记 \(\theta_n \subseteq R[x]_{n+1}\) 为满足条件 \(\forall p_n \in \theta_n\) ，都有 \(p_n(z)=1\) 的集合，那么： \[ \left\Vert\frac{T_n(x)}{T_n(z)}\right\Vert_\infty=\mathop{inf}_\limits{p_n\in \theta_n}\Vert p_n \Vert_\infty \] 类似地，记 \(\bar{T}_n(x)\) 为 \(T_n(x)\) 的首一化多项式，也就是：\(\bar{T}_n(x)=\frac1{2^{n-1}}T_n(x)\)，那么任给 \(n\) 次首一多项式 \(p_n(x)\)，都有： \[ \Vert p_n\Vert_\infty\geq\Vert\bar{T}_n\Vert_\infty=\frac1{2^{n-1}}\:\:\:\:(W(x)=1) \] 当且仅当 \(p_n=\bar{T}_n\) 时取等，这说明要想使得对定义在 \([-1,1]\) 上的函数 \(f\) 作 \(n\) 次多项式插值得到的余项较小，可以将插值基点选取为 \(T_{n+1}(x)\) 的零点，即： \[ x_i=cos\frac{(2i+1)\pi}{2(n+1)},i=0,1,...,n \] 对于一般的定义在 \([a,b]\) 上的函数 \(f\)，作平移和伸缩变换再类似处理即可，也就是： \[ [a,b]\to[-1,1]:x\to\frac2{b-a}(x-\frac{b-a}2)\\ \frac2{b-a}(x_i-\frac{b-a}2)=cos\frac{(2i+1)\pi}{2(n+1)},i=0,1,...,n \]
由此可以得到一个近似的最佳一致逼近多项式，关于最佳一致逼近的具体讨论在后续部分
第二类 \(Chebyshev\) 多项式
与第一类的唯一差别是权函数取为 \(W(x)=\sqrt{1-x^2}\)，将其记为 \(\{U_n(x)\}\)，可以表示为： \[ U_n(x)=\frac{sin\left((n+1)arccos\:x\right)}{\sqrt{1-x^2}} \] 其与第一类有关系式： \[ T_{n+1}'(x)=(n+1)U_n(x) \]
\(Legendre\) 多项式
取最简单的权函数 \(W(x)=1\)，区间 \([a,b]=[-1,1]\)，得 \(Legendre\) 多项式，可以表示为： \[ p_n(x)=\frac1{2^nn!}\frac{d^n}{dx^n}(x^2-1)^n \] 其用处在这门课主要体现在高斯型求积公式中，此处暂且不谈
最佳一致逼近
取 \(W(x)=1\)，设 \(f\in C[a,b]\)，若有一列函数 \(\{f_n\}\) 使得： \[ \lim_{n\to\infty}\Vert f-f_n\Vert_\infty=0 \] 称 \(\{f_n\}\) 为 \(f\) 的一致逼近函数，不难看出也就是在 \([a,b]\) 上一致收敛于 \(f\)
我们在数学分析的课程上已经知道 \(Weierstrass\) 给出的关于多项式一致逼近的相关结果，即 \(\forall \varepsilon>0,\exists p_n(x)\) 使得： \[ \Vert f-p_n\Vert_\infty<\varepsilon \] 然而我们知道多项式拟合的次数太高并不好，下面我们固定 \(n\)，考虑在次数不超过 \(n\) 的多项式集合 \(H_n\) 中寻求 \(p_n\) 使得： \[ \Vert f-p_n\Vert_\infty=\mathop{inf}_\limits{q_n\in H_n}\Vert f-q_n\Vert_\infty \] 这就是最佳一致逼近多项式问题，关于这个问题，\(Chebyshev\) 给出了以下结果：
设 \(a\leq x_1<x_2<...<x_n\leq b\)，称 \(\{x_1,x_2,...,x_n\}\) 为 \(f\) 在 \([a,b]\) 上的交错点组，若满足： \[ |f(x_i)|=\mathop{max}_\limits{x\in[a,b]}|f(x)|,\forall i=1,2,...,n\\ f(x_i)=-f(x_{i+1}),\forall i=1,2,...,n-1 \]
（\(Chebyshev\) 定理）设 \(p\in H_n\)，则 \(p\) 是 \(f\) 的最佳一致逼近多项式当且仅当 \(f(x)-p(x)\) 在 \([a,b]\) 上存在 \(n+2\) 个点组成的交错点组
上述 \(p(x)\) 在 \(H_n\) 中是唯一的
若 \(f^{(n+1)}(x)\) 在 \([a,b]\) 中不变号，则 \(a,b\) 属于交错点组
最佳平方逼近
与最佳一致逼近类似，只是考虑的范数由 \(\infty\) 范数改为了 \(2\) 范数，我们考察一般情形：
设 \(f\in C[a,b],\:\{\varphi_n\}\) 为 \([a,b]\) 上的线性无关函数族，我们希望寻求 \(a_0,a_1,...,a_n\) 使得： \[ \Vert f-\sum_{k=0}^na_k\varphi_k\Vert_2 \] 达到极小；考察其平方不影响结果，从而可以使用多元函数求极值的办法给出结论： \[ F(a_0,a_1,...,a_n)=\int_a^b\left(f(x)-\sum_{k=0}^na_k\varphi_k(x)\right)^2W(x)dx\\ \frac{\partial F}{\partial a_j}=2\int_a^b\left(f(x)-\sum_{k=0}^na_k\varphi_k(x)\right)(-\varphi_j(x))W(x)dx=0 \] 也就是： \[ \left( \begin{array}{} (\varphi_0,\varphi_0) & (\varphi_0,\varphi_1) &... &(\varphi_0,\varphi_n)\\ (\varphi_1,\varphi_0) & (\varphi_1,\varphi_1) &... &(\varphi_1,\varphi_n)\\ ... & ... & ... & ... \\ (\varphi_n,\varphi_0) & (\varphi_n,\varphi_1) &... &(\varphi_n,\varphi_n) \end{array} \right) \left( \begin{array}{} a_0\\ a_1\\ ...\\ a_n \end{array} \right)=\left( \begin{array}{} (f,\varphi_0)\\ (f,\varphi_1)\\ ...\\ (f,\varphi_n) \end{array} \right) \] 这被称为法方程，据 \(\{\varphi_n\}\) 线性无关知其系数矩阵行列式不为 \(0\)，从而存在唯一解
特别地，若 \(\{\varphi_n\}\) 为正交函数系，则系数矩阵为对角矩阵，进而直接有： \[ a_k=\frac{(f,\varphi_k)}{(\varphi_k,\varphi_k)} \] 进一步，如果让 \(n\to\infty\)，这就给出了 \(f\) 的广义 \(Fourier\) 级数
若正交函数系取为 \(Chebyshev\) 多项式 \(\{T_n(x)\}\)，由此就得到 \(f\) 在 \([-1,1]\) 上的 \(Chebyshev\) 级数： \[ f(x)=\frac12a_0+\sum_{n=1}^\infty a_nT_n(x)\\ a_k=\frac{(T_k,f)}{(T_k,T_k)}=\frac2\pi\int_{-1}^1T_k(x)f(x)\frac{dx}{\sqrt{1-x^2}} \]
离散情形
与连续情形类似，若已知 \(f\) 在给定点集 \(\{x_1,x_2,...,x_m\}\) 上的值，记： \[ r_i=\sum_{j=0}^na_j\varphi_j(x_i)-f(x_i) \] 为残量，亦可寻求 \(a_0,a_1,...,a_n\) 使得： \[ \sum_{i=0}^mr_i^2 \] 达到极小，这样的问题被称为离散的最佳平方逼近问题，也就是线性最小二乘拟合问题，在后续部分讨论
在等距点集情形，我们有著名的 DFT 以及 FFT，此前已作过其浅析，详见：https://dasasdhba.github.io/study-FFT/
数据拟合
数据拟合问题是指寻求一函数 \(f(x)\) 来拟合一组离散数据 \(\{(x_i,y_i)\}_{i=1}^m\)，与插值法不同，我们采用均方误差（mean-square error MSE）来衡量拟合的好坏，也就是考虑： \[ \mathcal{L}(f)=\frac1m\sum_{i=1}^m(f(x_i)-y_i)^2 \] 一般地，给定函数族 \(\mathcal{F}\)，我们试图求解 \(f\in\mathcal{F}\) 使得： \[ \mathcal{L}(f)=\mathop{inf}_\limits{g\in\mathcal{F}}\mathcal{L}(g) \]
线性拟合
下面介绍线性拟合的最小二乘法，这与最佳平方逼近是类似的
设 \(f\in C[a,b],\:\{\varphi_n\}\) 为 \([a,b]\) 上的线性无关函数族，我们希望寻求 \(a_0,a_1,...,a_n\) 使得： \[ \varphi(x)=\sum_{j=0}^na_j\varphi_j\\ \mathcal{L}(\varphi)=\frac1m\sum_{i=1}^m(\varphi(x_i)-y_i)^2 \] 达到极小；类似地，可以使用多元微积分求极值的理论： \[ F(a_0,...,a_n)=\sum_{i=1}^m(\sum_{j=0}^na_j\varphi_j(x_i)-y_i)^2\\ \frac{\partial F}{\partial a_k}=2\sum_{i=1}^m\varphi_k(x_i)(\sum_{j=0}^na_j\varphi_j(x_i)-y_i)=0 \] 记 \(Y=(y_1,...,y_m)^T,\:\Phi_i=(\varphi_i(x_1),...,\varphi_i(x_m))^T\)，取 \(\mathbb{R}^m\) 这个 \(Euclid\) 空间的通常内积 \((X,Y)=X^TY\)，则上述方程可以改写为： \[ \left( \begin{array}{} (\Phi_0,\Phi_0) & (\Phi_0,\Phi_1) &... &(\Phi_0,\Phi_n)\\ (\Phi_1,\Phi_0) & (\Phi_1,\Phi_1) &... &(\Phi_1,\Phi_n)\\ ... & ... & ... & ... \\ (\Phi_n,\Phi_0) & (\Phi_n,\Phi_1) &... &(\Phi_n,\Phi_n) \end{array} \right) \left( \begin{array}{} a_0\\ a_1\\ ...\\ a_n \end{array} \right)=\left( \begin{array}{} (Y,\Phi_0)\\ (Y,\Phi_1)\\ ...\\ (Y,\Phi_n) \end{array} \right) \] 由此即可解出 \(a_0,...,a_n\)；特别地，若取 \(n=1, \varphi_0(x)=1,\varphi_1(x)=x\)，得到： \[ \begin{aligned} &\left\{\begin{align} &a_1=\frac{(\Phi_1,Y)-m\bar{x}\bar{y}}{(\Phi_1,\Phi_1)-m\bar{x}\bar{y}}=\frac{\sum_{i=1}^m(x_i-\bar{x})(y_i-\bar{y})}{\sum_{i=1}^m(x_i-\bar{x})^2}\\ &a_0=\bar{y}-a_1\bar{x}\\ \end{align}\right.\\ &其中:\\ &\bar{x}=\sum_{i=1}^mx_i,\:\bar{y}=\sum_{i=1}^my_i \end{aligned} \] 即中学阶段已经介绍过的最小二乘法求解回归直线方程
过拟合现象
若取 \(\{\varphi_n\}=\{x^n\}\)，直观上 \(n\) 越大误差越小，然而与 \(Runge\) 现象类似，过高的拟合次数会产生过拟合现象
非线性拟合
非线性拟合是考虑用不能被一族函数的线性组合所表出的函数作拟合，如： \[ f(x)=\frac a{1+be^{cx}} \] 我们的目标是寻求待定系数 \(a,b,c\) 使得 \(\mathcal{L}(f)\) 达到极小
线性化求近似解
若有一向量值函数 \(F:\mathbb{R}^2 \to \mathbb{R}^2\) 使得 \(\{F(x_i,y_i)\}_{i=1}^m\) 可以用一族函数的线性组合来拟合，则可以使用线性拟合的办法求近似解
（化曲为直）
举例，若离散数据符合非线性函数： \[ y=\frac1{1+ae^{bx}} \] 这可以改写为： \[ ln\left(\frac 1y-1\right)=bx+ln\:a \] 这就把问题转换为了线性拟合，注意到由此得到的是： \[ \bar{\mathcal{L}}(a,b)=\frac1m\sum_{i=0}^m\left|bx_i+lna-ln\left(\frac 1y_i-1\right)\right|^2 \] 达到极小，故这是原非线性拟合问题的一个近似解
梯度下降法
考虑使用非线性函数 \(f(x;a_1,...,a_n)\) 作非线性拟合，其中 \(a_1,...,a_n\) 为待定系数，我们试图寻求待定系数使得： \[ F(a_1,...,a_n)=\mathcal{L}(f(x;a_1,...,a_n)) \] 达到极小；据多元微积分理论，\(\forall a \in \mathbb{R}^n\)，\(F\) 在 \(a\) 处沿着梯度方向变化最快，其中梯度是指： \[ \nabla F(a)=\left(\frac{\partial F}{\partial a_1},...,\frac{\partial F}{\partial a_n}\right)(a) \] （方便起见，这里没有加上转置符号）
因而可以通过从某一初值 \(a\) 开始，沿着梯度方向不断迭代的办法近似求解 \(F\) 的极小值点，这被称为梯度下降法，算法如下： \[ \begin{aligned} &a=(b_1,...,b_n)\:为初值\\ &\eta=10^{-4}\:为步长，又称学习率\\ &\varepsilon=10^{-8}\:为终止条件\\ \\ &while\;\Vert\nabla F(a)\Vert_2 \geq \varepsilon:\\ &\;\;\;\;a=a-\eta\:\nabla F(a) \end{aligned} \] 其中 \(\eta\) 的选取至关重要，太大没法收敛，太小收敛太慢；另一方面，这种算法显然只能得到局部极小值，不能保证其为全局极小值点，因而适当的选取初值也较为重要，可以通过线性化先求解一近似解作为初值
数值微积分
假设需要处理的函数十分复杂，又或者我们只知道一些离散点 \(\{x_n,y_n\}\) 的数据，此时需要考察求解微分以及积分的近似值的方法。
数值微分
数值微分比较容易，因为只需要近似求解一个极限，下面记需要考察的函数为 \(f\)
差商法
据定义，可直接取 \(\frac{f(x+h)-f(x)}h\) 作为近似值，即向前差商法；同理还有向后差商法和中心差商法。
以向前差商法为例考察误差，据 \(Taylor\) 展开，有： \[ f(x+h)=f(x)+hf'(x)+\frac{h^2}{2!}f''(\xi) \] 由此不难看出误差为 \(O(h)\)
插值法
取插值函数 \(L(x)\simeq f(x)\)，并以插值函数的导数作为近似的方法
之前已讨论过插值函数的误差，从而不难估计此法的误差
数值积分
假设我们只知道离散点 \(\{x_n,f(x_n)\}\) 的数据，我们尝试考察 \(f\) 在 \([a,b]\) 上关于权函数 \(W(x)\) 的积分： \[ I(f)=\int_a^bf(x)W(x)dx \] 称基于这些数据给出的公式： \[ I_n(f)=\sum_{i=0}^na_if(x_i) \] 为数值积分公式，其中 \(a_i\) 为待定系数。为衡量公式的精确度，称 \(I_n(f)\) 有 k 阶代数精度，若： \[ I_n(x^i)=I(x_i),i=0,1,...,k,\:I_n(x^{k+1})\neq I(x^{k+1}) \]
插值法
与数值微分不同，若想直接使用积分定义作近似，作 \(Riemann\) 和难以估计误差，而以 \(Lebesgue\) 积分的观点来考察简单函数的逼近则更加困难，因而数值积分基本上只有插值法。以 \(Lagrange\) 插值为例： \[ I_n(f)=\int_a^bL_n(x)W(x)dx=\sum_{i=0}^n\left(\int_a^bl_i(x)W(x)dx\right)f(x_i) \] 据 \(Lagrange\) 插值的性质，不难看出此法的代数精度至少为 \(n\)，误差由插值函数的误差作积分给出
函数的数值积分
与之前讨论的仅知道离散点数据的情形不同，以下假设我们需要考察的被积函数 \(f\) 为已知函数，积分区间为 \([a,b]\)
尽管仍然只能采取插值法，然而我们可以自由地选取插值基点，由此可以深入研究各种方法
Newton-Cotes 求积公式
我们考察 \(W(x)=1\)，插值基点为 \(n+1\) 个等距基点（步长 \(h=\frac{b-a}n\)）的情形，此时有： \[
a_i=\int_a^bl_i(x)=\frac {(-1)^{n+1-i}\;h}{(i-1)!(n+1-i)!} \int_0^n\prod_{j=0,j\neq i-1}^n(t-j)dt
\] 若 \(n\) 为偶数，且 \(f(x)\) 在 \([a,b]\) 上有 \(n+2\) 阶连续导数，那么误差： \[
E_n(f)=\frac{h^{n+3}f^{(n+2)}(\eta)}{(n+2)!}\int_0^nt\prod_{j=0}^n(t-j)dt
\] 若 \(n\) 为奇数，且 \(f(x)\) 在 \([a,b]\) 上有 \(n+1\) 阶连续导数，那么误差： \[
E_n(f)=\frac{h^{n+2}f^{(n+1)}(\eta)}{(n+1)!}\int_0^nt\prod_{j=0}^n(t-j)dt
\] 还需要考虑舍入误差可能带来的数值不稳定性，有估计： \(|I_n(f)-I_n^*(f)|\leq\varepsilon\sum_{i=0}^n|a_i|\)，其中 \(\varepsilon\) 为最大的舍入误差
取 \(f(x)=1\)，此时若 \(n\geq8\)，那么存在 \(a_i<0\)，这会使得数值不稳定
通常我们只采用 \(n=1\) 或 \(n=2\) 的情形（原因见下节），这分别对应梯形公式与 \(Simpson\) 公式
梯形公式：\(I_1(f)=\frac{b-a}2(f(a)+f(b))\)
直观上，这是以梯形面积作为近似，因而叫梯形公式
\(Simpson\) 公式：\(I_2(f)=\frac{b-a}6\left(f(a)+4f(\frac{a+b}2)+f(b)\right)\)
直观上，这是以两端以及中点所作抛物线的面积作为近似，因而也叫抛物线公式
复合求积公式
前已述不能考虑通过提高 \(Newton-Cotes\) 求积公式的阶数来减少误差，而另一方面，积分具有区间可加性，这启发我们将 \([a,b]\) 分成若干小段，并对每一段分别使用梯形公式或 \(Simpson\) 公式再求和，这就是复合求积公式。
下面考虑将 \([a,b]\) 分为 \(n\) 段，进而步长为 \(h=\frac{b-a}n\)，为了衡量复合求积公式的精确度，称公式 \(p\) 阶收敛，若： \[ \lim_{h\to0}\frac{I(f)-I_n(f)}{h^p}=C<\infty,\:C\neq0 \] 下面分别考虑使用梯形公式和 \(Simpson\) 公式作复合求积：
复合梯形公式：\(T_n(f)=\frac h2(f(a)+f(b)+2\sum_{i=1}^{n-1}f(a+ih))\)
误差：\(E_n(f)=-\frac{h^2(b-a)}{12}f''(\xi)\)，因而是 \(2\) 阶收敛的
复合 \(Simpson\) 公式（为统一 \(h\) 的值，假定 \(n\) 为偶数，将 \([a,b]\) 分为 \(m=\frac n2\) 段）： \[ S_n(f)=\frac h3\left(f(a)+f(b)+4\sum_{i=1}^mf(a+(2i-1)h)+2\sum_{i=1}^{m-1}f(a+2ih)\right) \] 误差：\(E_m(f)=-\frac{h^4(b-a)}{180}f^{(4)}(\xi)\)，因而是 \(4\) 阶收敛的
可见，使用复合求积公式，只要通过控制 \(h\)，就可以有效地控制误差，因而我们并不需要使用高阶的 \(Newton-Cotes\) 求积公式
区间逐次分半法
书接上回，使用复合求积公式时，若想控制误差界为给定的 \(\varepsilon\)，由其误差估计式，往往需要考察 \(f\) 的高阶导数，这会带来麻烦。
为此，常使用一种被称为区间逐次分半的办法——不断将区间二分，直到达到精度要求
由于我们不希望直接考察误差估计式，常常采用以下办法近似（以复合梯形公式为例）：
我们知道： \[
E_n(f)=I(f)-T_n(f)=-\frac{h^2(b-a)^3}{12n^2}f''(\xi_1)\\
E_{2n}(f)=I(f)-T_{2n}(f)=-\frac{h^2(b-a)^3}{12\times4n^2}f''(\xi_2)
\] 如果我们认为：\(f''(\xi_1)\simeq f''(\xi_2)\)，那么 \(E_n(f)\simeq 4E_{2n}(f)\)，由此可以推出： \[
E_{2n}(f)\simeq \frac13(T_{2n}(f)-T_n(f))
\] 这样，我们就回避了直接考察 \(f''\)
另一方面，若记 \(n=2^{m-1}\) ，即将 \([a,b]\) 作 \(2^{m-1}\) 等分，记由此得到的梯形值为 \(T_{m,1}\)（下标与 \(T_n\) 作区分，无其他含义），那么： \[ T_{m,1}=\frac12\left(T_{m-1,1}+h_{m-1}\sum_{k=1}^{2^{m-2}}f(a+(k-\frac12)h_{m-1})\right)\\ h_m=\frac{b-a}{2^{m-1}} \] 即每一次半分时存在递推关系，降低了计算量；事实上，后面会介绍 \(Romberg\) 积分法，想法与此十分相似
自适应积分法
在区间逐次分半法中，每次我们都将区间的等分数量翻倍，尽管有递推式，但这样做仍然会导致速度较慢。
然而事实上，我们可以将需要控制的误差界 \(\varepsilon\) 在等分的 \(n\) 个区间上分别考虑，
只要每个区间上的误差都不超过 \(\frac\varepsilon n\)，那么也能达到同样的目的
这样，我们只需要对那些误差超过 \(\frac\varepsilon n\) 的区间应用逐次分半法，这就是自适应积分法。
不难看出，这种办法很适合递归实现：
1 | function y = self_adaptive_integral(f,a,b,eps) |
其中 numeric_integral() 是某个数值积分公式，numeric_integral_err() 往往类似于区间逐次分半法中的误差估计方法
需要指出的是，由于分半的过程中，各个区间需要处理的误差界也在减半，因而此法容易陷入无限递归，应当加上递归层数的限制
\(Euler-Maclaurin\) 公式
后面会稍微用到，但用得不多，故只给出相关结果，证明略
记 \(B_n\) 为 \(Bernoulli\) 数，这是由：\(\frac{x}{e^x-1}=\sum_{n=0}^\infty \frac{B_n}{n!}x^n\) 给出的，显式公式为： \[ B_n=\sum_{k=0}^n\frac1{k+1}\sum_{r=0}^kC_k^r(-1)^rr^n \] 设 \(a,b\in\mathbb{N},\:f\in C^\infty[a,b]\)，那么有 \(Euler-Maclaurin\) 公式： \[ \sum_{n=a}^bf(n)-\int_a^bf(x)dx=\frac12(f(a)+f(b))+\sum_{n=1}^\infty \frac{B_{2n}}{2n!}(f^{(2n-1)}(b)-f^{(2n-1)}(a)) \] （简单来说就是利用 \(Bernoulli\) 多项式的性质反复分部积分得到的）
一般地，若 \(f\in C^k[a,b]\)，那么同样有： \[ \sum_{n=a}^bf(n)-\int_a^bf(x)dx=\frac12(f(a)+f(b))+\sum_{n=2}^{k-1} \frac{B_{n}}{n!}(f^{(n-1)}(b)-f^{(n-1)}(a))+R_k(f) \] 为了描述余项，需要引入 \(Bernoulli\) 多项式，类似地，这由：\(\frac{te^{tx}}{e^t-1}=\sum_{n=0}^\infty B_n(x)\frac{t^n}{n!}\) 给出，有显式公式： \[ B_n(x)=(B+x)^n \] 其中 \(B^k=B_k\) 为 \(Bernoulli\) 数（只是个方便的记号，\(B\) 并不是一个数）
于是，记 \(\{x\}=x-[x]\) 为 \(x\) 的小数部分，那么余项就可以写为： \[ R_k(f)=(-1)^{k-1}\int_a^b\frac{B_k(\{x\})}{k!}f^{(k)}(x)dx \] 对于 \(a,b\in\mathbb{R}\) 的情形，可以固定 \(n\in\mathbb{N}\)，由 \(x\to \frac{b-a}nx+a\) 将 \([0,n]\to[a,b]\)，再使用上述公式
如果记 \(h=\frac{b-a}n\)，整理可以得到： \[ \sum_{i=0}^nhf(a+ih)-\int_a^bf(x)dx=\frac h2(f(a)+f(b))+\sum_{m=2}^{k-1} \frac{B_{m}}{m!}h^{m}(f^{(m-1)}(b)-f^{(m-1)}(a))+R_k(f)\\ R_k(f)=(-1)^{k-1}h^k\int_a^b\frac{B_k(\{\frac{x-a}h\})}{k!}f^{(k)}(x)dx \]
注意， \(\sum_{i=0}^nhf(a+ih)-\frac h2(f(a)+f(b))\) 恰为复合梯形公式！
故 \(Euler-Maclaurin\) 公式事实上给出了复合梯形公式的一个极为精确的误差估计
\(Romberg\) 积分法
沿用区间逐次分半法的相关记号，将 \([a,b]\) 作 \(2^{m-1}\) 等分，记由此得到的梯形值为 \(T_{m,1}\)
\(Romberg\) 积分法实质上是利用 \(Euler-Maclaurin\) 公式对此做了改进，从而大大加快了收敛速度
事实上，据 \(Euler-Maclaurin\) 公式，有： \[
I(f)-T_{m,1}=a_2h^2+a_4h^4+O(h^6)
\] 其中 \(a_2,a_4\) 只与 \(f\) 在 \(a,b\) 处的相关导数值有关，若进一步考察 \(T_{m-1,1}\)，有： \[
I(f)-T_{m-1,1}=2^2a_2h^2+2^4a_4h^4+O(h^6)
\] 由此，我们容易消去误差中的 \(h^2\) 项，从而直接将误差缩小至 \(O(h^4)\)！
我们记： \[
T_{m,2}=\frac{4T_{m,1}-T_{m-1,1}}{3}
\] （不难验证，这恰为复合 \(Simpson\) 公式）
仿此，可以继续考虑： \[ I(f)-T_{m,2}=b_4h^4+b_6h^6+O(h^8)\\ I(f)-T_{m-1,2}=2^4b_4h^4+2^6b_6h^6+O(h^8)\\ \Rightarrow T_{m,3}=\frac{16T_{m,2}-T_{m-1,2}}{15} \] 一般地，有： \[ T_{m,j}=\frac{4^{j-1}T_{m,j-1}-T_{m-1,j-1}}{4^{j-1}-1} \] 再结合区间逐次分半法中提到的递推式： \[ T_{m,1}=\frac12\left(T_{m-1,1}+h_{m-1}\sum_{k=1}^{2^{m-2}}f(a+(k-\frac12)h_{m-1})\right)\\ T_{1,1}=\frac{h_1}2(f(a)+f(b)),\;h_m=\frac{b-a}{2^{m-1}} \] 这就给出了 \(Romberg\) 积分法，计算顺序呈三角式，也就是： \[ \begin{aligned} &T_{1,1}\\ &T_{2,1} &T_{2,2}\\ &T_{3,1} &T_{3,2} &&T_{3,3}\\ &... \end{aligned} \]
\(Gauss\) 求积公式
到目前为止，我们的讨论大都基于 \(Newton-Cotes\) 求积公式——一个选取等距基点的方法
一个自然的问题是，若选取一些特殊的基点作插值，是否可以达到更好的效果？这就是 \(Gauss\) 求积公式要讨论的问题
之前已经讨论过，\(n+1\) 个插值基点 \(a\leq x_0<x_1<...<x_n\leq b\) 给出的求积公式： \[ I_n(f)=\int_a^bL_n(x)W(x)dx=\sum_{i=0}^n\left(\int_a^bl_i(x)W(x)dx\right)f(x_i) \] 的代数精度至少为 \(n\)，如果记 \(\omega_{n+1}(x)=\prod_{i=0}^n(x-x_i)\)，则上式又可以写为： \[ I_n(f)=\sum_{i=0}^n\left(\int_a^b\frac{\omega_{n+1}(x)}{(x-x_i)\omega'_{n+1}(x)}W(x_i)dx\right)f(x_i) \] 可以证明，上述求积公式的代数精度至多是 \(2n+1\)，且达到 \(2n+1\) 的充分必要条件为：
\(\omega_{n+1}(x)\) 与不高于 \(n\) 次的任何多项式都关于区间 \([a,b]\) 上的权函数 \(W(x)\) 直交
事实上，如果 \(p(x)\) 为一个不高于 \(2n+1\) 次的多项式，由带余除法： \[ p(x)=q(x)\omega_{n+1}(x)+r(x) \] 再结合直交的条件就不难看出结论成立
一般地，我们称选取了 \(n+1\) 个插值基点的求积公式为 \(Gauss\) 型求积公式，若其代数精度为 \(2n+1\)
下面我们简单考虑其误差，因为代数精度为 \(2n+1\)，可选 \(Hermite\) 插值逼近 \(f\) 来作误差估计，由此可以证明离散误差为： \[ E_n(f)=Cf^{(2n+2)}(\xi)\\ C=\frac{d^2}{(2n+2)!}\int_a^bp^2_{n+1}(x)dx \] 其中，\(p_{n+1}(x)\) 为 \([a,b]\) 上关于权函数 \(W(x)\) 的 \(n+1\) 次直交多项式，\(d\) 为 \(p_{n+1}(x)\) 首项系数的倒数
下面我们讨论几个特殊情形：
\([a,b]=[-1,1],\;W(x)=1\)，我们知道此时直交多项式为 \(Legendre\) 多项式，取 \(n=1\)，有 \(P_2(x)=\frac12(3x^2-1)\)
因而只要取 \(P_2(x)\) 的两个根 \(x=\pm\frac1{\sqrt{3}}\) 为插值基点，就得到代数精度为 \(3\) 的插值公式 \[ \int_{-1}^1f(x)dx\simeq f(-\frac1{\sqrt{3}})+f(\frac1{\sqrt{3}}) \]
称为 \(Gauss-Legendre\) 两点公式，离散误差 \(E_1(f)=\frac1{135}f^{(4)}(\xi)\)
\([a,b]=[-1,1],\;W(x)=\frac1{\sqrt{1-x^2}}\)，我们知道此时直交多项式为 \(Chebyshev\) 多项式 \(T_{n+1}(x)=cos((n+1)arccos\:x)\) 因而只要取插值基点为 \(x_i=cos\left(\frac{2i-1}{2(n+1)}\pi\right),\;i=1,2,...,n+1\)，就得到 \(Gauss-Chebyshev\) 求积公式
利用 \(Chebyshev\) 多项式的相关性质，化简后的最终结果为： \[ \int_{-1}^1f(x)\frac{dx}{\sqrt{1-x^2}}\simeq\frac{\pi}{n+1}\sum_{k=1}^{n+1}f(x_i) \]
对于一般的情形，我们可以仿照复合求积公式的想法，首先将 \([a,b]\) 等分为若干小区间，然后在每个小区间上通过积分变换将其变换为 \([-1,1]\) 上的积分，应用上述公式即可
重积分与广义积分
重积分只要处理为累次积分，就可以导出类似的积分公式（当然，概率论上有一个著名的重积分计算方法，不过此处不作介绍）
广义积分通常有两种办法处理，一是使用 \(Gauss\) 型求积公式，只需要寻找相关的直交多项式；二是截断法，我们知道广义积分收敛的必要条件是在奇点附近的积分趋于 \(0\)，从求解近似解的角度来讲，直接截去即可
ODE 初值问题数值解
给定一个存在唯一解的 ODE 初值问题：
\[ \left\{\begin{aligned} &y'=f(t,y),\;a\leq t \leq b\\ &y(a)=\eta \end{aligned}\right. \]
称在一系列离散点 \(t_0,t_1,...,,t_n\) 处的近似值 \(y_0,y_1,...,y_n\) 为该问题的数值解（称为离散变量法）
通常采用递推的方法，假定我们已知 \(y_0,y_1,...,y_{k-1}\)，我们设计算法：\(y_k\simeq \varphi(y_0,y_1,...,y_k)\)
注意，算法两边可能同时出现 \(y_k\)，这种情况我们称方法是隐式的（即需要解方程），否则是显式的
通常，我们有三种思路：
差商代替倒数，即：\(y'(t_k)\simeq \frac{y(t_{k+1})-y(t_k)}{h_k},\:h_k=t_{k+1}-t_k\)，于是： \[ y(t_{k+1})\simeq y(t_k)+h_kf(t_k,y(t_k)) \] 这被称为 \(Euler\) 法
\(Taylor\) 级数法：取等距离散点，步长 \(h=t_k-t_{k-1}\)，视 \(y(t+h)\) 为关于 \(h\) 的函数，考察 \(h=0\) 处的 \(Taylor\) 展开： \[ \begin{aligned} &y(t+h)=y(t)+hy'(t)+\frac12h^2y''(t)+\frac1{3!}h^3y'''(t)+...+\frac1{p!}h^py^{(p)}(t)+\frac1{(p+1)!}h^{p+1}y^{(p+1)}(\xi)\\ &=y(t)+hf(t,y)+\frac12h^2\frac d{dt}f(t,y)+\frac1{3!}h^3\frac {d^2}{dt^2}f(t,y)+...+\frac1{(p+1)!}h^{p+1}\frac{d^{p}}{dt^{p}}f(\xi,y(\xi)) \end{aligned} \] 截去余项就得到递推公式，不难看出 \(p=1\) 时这就是 \(Euler\) 法
数值积分法：\(y(t_{k+1})-y(t_k)=\int_{t_k}^{t_{k+1}}f(t,y)dt\)，将右边替换为数值积分公式即可
为了研究误差，可以假定 \(y_0,y_1,...,y_{k-1}\) 是精确的，称这种情况下算法 \(\varphi(y_0,y_1,...,y_k)\) 的误差为局部离散误差
记 \(h_k=t_{k+1}-t_{k}\)，我们称算法为 \(p\) 阶方法，若 \(y_{k+1}\) 的局部离散误差为 \(O^*(h_k^{p+1})\)
从整体上来看，我们称 \(y_n\) 与其精确解之间的误差为整体离散误差
单步法
单步法是说递推算法形如：\(y_k\simeq \varphi(y_{k-1},y_k)\)，即每次递推只用到上一步的数据
对于等距、显式的情形，我们通常写为： \[ y_{k+1}\simeq y_k+h\Phi(t_k,y_k,h) \] 我们首先讨论一些具体的方法，然后再从理论上系统研究误差问题
\(Euler\) 法及其改进
前已述 \(Euler\) 法的原理，在等距情形，也就是 \(\Phi(t_k,y_k,h)=f(t_k,y_k)\)
因为 \(Euler\) 法可以由 \(Taylor\) 级数法取 \(p=1\) 得到，这就不难看出其局部离散误差为 \(O^*(h^2)\)，因而是一阶方法
此法精度不高，故需改进，如果我们采用数值积分法中的梯形公式，可以得到一隐式方法： \[
y_{k+1}\simeq y_k+\frac h2(f(t_k,y_k)+f(t_{k+1},y_{k+1}))
\] 据梯形公式的误差估计，我们知道此法为二阶方法，然而此法作为隐式方法，一般需要通过迭代求近似解
迭代求解需要一个初值，我们可以使用 \(Euler\) 法提供，具体来说，也就是： \[
\begin{aligned}
&y_{k+1,0}\simeq y_k+hf(t_k,y_k)\\
&y_{k+1,m}\simeq y_k+\frac h2(f(t_k,y_k)+f(t_{k+1},y_{k+1,m-1}))
\end{aligned}
\] 这个过程也被称作 PECE，即 Predict-Evaluate-Correct-Evaluate，迭代 \(k\) 次的模式记为 \(P(EC)^kE\)
简单来说，就是先用显式方法作预测，再将该预测值用于隐式方法作校正（并反复迭代）的过程。
（注：如果对预测值和校正值的离散误差有较好的把握，可以通过 \(Romberg\) 积分法类似的想法去进一步修正误差）
当 \(h\) 足够小时，迭代一两次就已经能够得到较好的结果，由此就得到了改进的 \(Euler\) 方法： \[ y_{k+1}\simeq y_k+\frac h2(f(t_k,y_k)+f(t_{k+1},y_k+hf(t_k,y_k))) \]
### \(Runge-Kutta\) 方法
若直接采用 \(Taylor\) 法，则有： \[ \Phi(t_k,y_k,h)=f(t_k,y_k)+\frac12h\frac d{dt}f(t_k,y_k)+\frac1{3!}h^2\frac {d^2}{dt^2}f(t_k,y_k)+...+\frac1{p!}h^{p-1}\frac{d^{p-1}}{dt^{p-1}}f(t_k,y_k) \] 此法看似容易，然而计算 \(\frac{d^p}{dt^p}f(t,y(t))\) 绝非易事，\(Runge-Kutta\) 给出了一个想法，绕过了求解偏导数的麻烦
以 \(p=2\) 为例（称为二阶 \(Runge-Kutta\) 方法），此时： \[
\Phi(t,y,h)=f(t,y(t))+\frac12 h \frac d{dt}f(t,y(t))
\] 其中： \[
\frac d{dt}f(t,y(t))=f_t'(t,y)+f_y'(t,y)\cdot y'(t)=f'_t(t,y)+f'_y(t,y)\cdot f(t,y)
\] 我们用待定系数法求一 \(\Phi\) 的近似解，如下： \[
\begin{aligned}
&\Psi(t,y,h)=c_1K_1+c_2K_2\\
&K_1=f(t,y)\\
&K_2=f(t+a_1h,y+a_2hK_1)
\end{aligned}
\] 用 \(Taylor\) 展开 \(K_2\)（视为 \(h\) 的函数），得到： \[
K_2=f(t,y)+a_1hf_t'(t,y)+a_2hf(t,y)f_y'(t,y)+O(h^2)
\] 于是： \[
\Psi(t,y,h)=(c_1+c_2)f(t,y)+c_2h(a_1f_t'(t,y)+a_2f_y'(t,y)f(t,y))+O(h^2)
\] 我们知道： \[
y(t+h)=y(t)+h\Phi(t,y,h)+O(h^3)
\] 要想使得上式将 \(\Phi\) 替换为 \(\Psi\) 仍然成立，这只要求： \[
(c_1+c_2)f(t,y)+c_2h(a_1f_t'(t,y)+a_2f_y'(t,y)f(t,y))=\Phi(x,y,h)
\] 比较各项系数，得方程组： \[
\left\{\begin{aligned}
&c_1+c_2=1\\
&c_2a_1=c_2a_2=\frac12\\
\end{aligned}\right.
\] 可自由选取一个变量的值，从而解出该方程，这样就导出了二阶 \(Runge-Kutta\) 方法
特别地，取 \(a_1=1\)，得到 \(a_2=1,\;c_1=c_2=\frac12\)，由此恰好得到改进的 \(Euler\) 方法
至今最常用的是一个四阶 \(Runge-Kutta\) 方法： \[ \begin{aligned} &y_{k+1}=y_k+\frac h6(K_1+2K_2+2K_3+K_4)\\ &K_1=f(t_k,y_k)\\ &K_2=f(t_k+\frac12h,y_k+\frac12hK_1)\\ &K_3=f(t_k+\frac12h,y_k+\frac12hK_2)\\ &K_4=f(t_k+h,y_k+hK_3) \end{aligned} \] 其局部离散误差为 \(O^*(h^5)\)
理论误差分析
我们从理论角度研究一些一般性的误差问题
相容性
对于等距情形，我们当然希望近似解关于 \(h\) 能够收敛到精确解，而这有一个必要条件： \[ \lim_{h\to0}\left(\frac{y(t+h)-y(t)}{h}-\Phi(t,y,h)\right)=0 \] 如果 \(\Phi\) 是连续函数，这个事情实际上是在说：\(y'(t)=\Phi(t,y,0)\)，这给出如下定义：
称单步法与初值问题是相容的，若：\(\Phi(t,y,0)=f(t,y)\)
如果只知道一个单步法是相容的，实际上并不能给出太多信息，这只能保证这个方法至少为一阶方法
收敛性
首先我们给出收敛性的定义：任给 \(t\in[a,b]\)，固定 \(t_n=t\)，若有 \(\lim_{h\to0}y_n=y(t)\)，则称方法是收敛的
对于单步法，我们有如下定理：
若 \(\Phi(t,y,h)\) 对 \(y\) 满足 \(Lipschitz\) 条件，\(\forall t,h\)，则单步法收敛的充分必要条件是相容条件成立
\(Lipshitz\) 条件就是说，\(\forall y_1,y_2\in\mathbb{R}\)，存在 \(L>0\;s.t.\) \[
|\Phi(t,y_1,h)-\Phi(t,y_2,h)|\leq L|y_1-y_2|,\;\forall t,h
\] 进一步，这给出了整体离散误差的一个估计，我们记 \(R(t,h)\) 为局部离散误差，对于 \(p\) 阶方法来说，应该有： \[
|R(t,h)|\leq Mh^{p+1}
\] 于是，只要 \(Lipschitz\) 条件满足，那么整体离散误差 \(\varepsilon_n=y(t_n)-y_n\) 就满足： \[
|\varepsilon_n|\leq e^{L(b-a)}|\varepsilon_0|+h^p\frac ML(e^{L(b-a)}-1)
\] 如果有 \(\varepsilon_0=0\)，那么这个 \(p\) 阶单步法的整体离散误差就为 \(O^*(h^p)\)
稳定性
稳定性是指数值稳定性，在误差分析一章中有提到；
在 ODE 初值问题中，我们可以通过考察对初值作扰动而产生的变化来考察稳定性，具体来说：
若存在 \(h_0,C>0\)，使得对任意的初值 \(y_0,\bar{y}_0\)，通过步长小于 \(h_0\) 的单步法得到的相应的解 \(y_n,\bar{y}_n\)，都满足： \[ |y_n-\bar{y}_n|\leq C|y_0-\bar{y}_0| \] 就称单步法是稳定的；事实上，只要 \(\Phi(t,y,h)\) 对 \(y\) 满足 \(Lipschitz\) 条件，那么单步法就是稳定的
而在实际应用中，讨论某个固定的 \(h\) 的情形更有意义，为此引入绝对稳定性：
称单步法是绝对稳定的，若由初值 \(y_k,\bar{y}_k\) 满足： \(\bar{y}_k=y_k+\delta\) 分别给出的 \(y_{k+1},\bar{y}_{k+1}\) 都有：\(\vert \bar{y}_{k+1}-y_{k+1}\vert<\vert\delta\vert\)
通常，我们只讨论初值问题 \(y'=\mu y\) 的绝对稳定性，称 \((\alpha,\beta)\) 为绝对稳定区间，若任给 \(\mu h\in(\alpha,\beta)\) 单步法都绝对稳定
举例：考虑 \(Euler\) 法的绝对稳定区间，这由： \[ |\bar{y}_{k+1}-y_{k+1}|=|1+\mu h||\delta|<|\delta| \] 不难看出 \(\mu h\in(-2,0)\)
多步法
与单步法相对的，多步法是说递推算法形如：\(y_k\simeq \varphi(y_{k-},...,y_k),\;m>1\)
我们仍然考虑一些具体方法，然后理论分析误差
\(Adams\) 方法
使用数值积分法 \(y(t_{k+1})-y(t_k)=\int_{t_k}^{t_{k+1}}f(t,y)dt\) 的问题在于右边的积分公式一般不能显式地给出
\(Adams\) 方法指出，可以用 \(y_{k-m},...,y_k\) 这些已知点作 \(f(t,y)\) 的插值多项式 \(p_m(t)\)，于是就得到近似公式： \[
y_{k+1}\simeq y_k+\int_{t_k}^{t_{k+1}}p_m(t)dt
\] 可以证明，\(m\) 步显式 \(Adams\) 方法至少为 \(m\) 阶方法
亦可以将 \(y_{k+1}\) 也作为插值点，由此得到方法被称为隐式 \(Adams\) 方法
可以证明 ，\(m\) 步隐式 \(Adams\) 方法至少为 \((m+1)\) 阶方法
（注：\(m=1\) 时，隐式 \(Adams\) 公式恰为梯形公式）
由于是隐式方法，通常也需要配合 PECE 模式使用，这里不详细叙述
还有一被称为 \(Milne\) 方法的多步法与 \(Adams\) 方法类似，只是将考虑的积分区间扩大： \(y(t_{n+1})-y(t_{n-p})=\int_{t_{n-p}}^{t_{n+1}}f(t,y)dt\)
然后仍然考虑使用插值多项式代替 \(f(t,y)\)，故也不详细叙述
\(Hamming\) 方法
与 \(Runge-Kutta\) 方法的想法类似，在等距情形，我们可以使用待定系数结合 \(Taylor\) 展开的方法来构造多步法，令： \[
L(y(t),h)=\sum_{j=0}^k(\alpha_jy(t+jh)-h\beta_jy'(t+jh))
\] 视为 \(h\) 的函数，考虑 \(h=0\) 处的 \(Taylor\) 展开，我们知道其必然形如： \[
L(y(t),h)=c_0y(t)+c_1hy'(t)+...+c_ph^py^{(p)}(t)+...
\] 其中 \(c_j\) 与 \(y(t)\) 无关，故我们分别令 \(y(t)=1,t,t^2,...\) 并取 \(t=0\)，就不难得到 \(c_j\) 的值： \[
\begin{aligned}
&c_0=\alpha_0+\alpha_1+...+\alpha_k\\
&c_1=\alpha_1+2\alpha_2+...+k\alpha_k-(\beta_0+\beta_1+...+\beta_k)\\
&...\\
&c_p=\frac1{p!}(\alpha_1+2^p\alpha_2+...+k^p\alpha_k)-\frac1{(p-1)!}(\beta_1+2^{p-1}\beta_2+...+k^{p-1}\beta_k)\\
&...
\end{aligned}
\] 要想使得 \(L(y(t),h)\) 导出一个至少 \(p\) 阶的线性多步法，这就要求：\(c_0=c_1=...=c_p=0\)
解出 \(\alpha_j,\beta_j\) 即可
注：\(Hamming\) 方法似乎只是 \(k=3,\;p=4,\;\alpha_k=1,\;\alpha_{k-1}=\beta_0=0\) 时导出的特例
理论误差分析
与单步法对应小节类似，我们考察形如： \[ \sum_{j=0}^ka_jy_{m+j}=h \sum_{j=0}^k\beta_j f(t_{m+j},y_{m+j}) \] 这样的线性 \(k\) 步法，我们记： \[ \rho(\lambda)=\alpha_k\lambda^k+\alpha_{k-1}\lambda^{k-1}+...+\alpha_1\lambda+\alpha_0\\ \sigma(\lambda)=\beta_k\lambda^k+\beta_{k-1}\lambda^{k-1}+...+\beta_1\lambda+\beta_0\\ \] 显然，这两个多项式与线性 \(k\) 步法一一对应，称 \(\rho(\lambda)\) 为线性 \(k\) 步法的特征多项式
仿照单步法的讨论，有如下结果：
相容性：称线性 \(k\) 步法是相容的，若它至少是一阶方法
可以证明，这等价于：\(\rho(1)=0,\;\rho'(1)=\sigma(1)\)
收敛性：称线性 \(k\) 步法是收敛的，若任给 \(t\in[a,b]\)，固定 \(t_n=t\)，都有 \(\lim_{h\to0}y_n=y(t)\)
对于线性 \(k\) 步法而言，难以考察收敛性，但是有等价关系：收敛 \(\Leftrightarrow\) 相容且稳定
稳定性：称线性 \(k\) 步法是稳定的，若存在 \(h_0,C>0\)，使得对任意的初值 \(y_0,\bar{y}_0\)
通过步长小于 \(h_0\) 的线性 \(k\) 步法得到的相应的解 \(y_n,\bar{y}_n\)，都满足： \[ |y_n-\bar{y}_n|\leq C\mathop{max}_\limits{0\leq j \leq k-1}|y_j-\bar{y}_j| \]可以证明，这等价于 \(\rho(\lambda)\) 满足一个所谓的特征根条件，即：
\(\rho(\lambda)\) 的所有根都在单位圆中，且在单位圆上的根只能是单重根
强稳定性：上述稳定性一般又称弱稳定性（误差可能会震荡），特征根条件称为弱根条件，若进一步满足：
\(\rho(\lambda)\) 的全部根除了 \(\lambda =1\) 外均在单位圆内（强根条件），则称满足强稳定性
绝对稳定性：仍然只讨论 \(y'=\mu y\)，此时线性 \(k\) 步法可以看成 \(k\) 阶常系数递推数列，有特征方程： \[ \rho(\lambda)-\mu h\sigma(\lambda) = 0 \] 我们称线性 \(k\) 步法是绝对稳定的，若上述方程的根 \(\lambda_r\) 都满足 \(|\lambda_r|<1\)
称 \((\alpha,\beta)\) 为绝对稳定区间，若任给 \(\mu h\in(\alpha,\beta)\) 多步法都绝对稳定可以证明，当 \(h\) 接近 \(0\) 时，想要达成绝对稳定必须要有 \(\mu h<0\)
高阶 ODE 和 ODE 方程组
对于一阶 ODE 方程组，只要简单地记 \(Y=(y_1,...,y_m),\;F=(f_1,...,f_m),\;\eta=(\eta_1,...,\eta_m)\)，那么初值问题就形如： \[ \left\{\begin{aligned} &Y'=F(t,Y),\;a\leq t \leq b\\ &Y(a)=\eta \end{aligned}\right. \] 之前讨论的方法对此是完全适用的，而对于高阶 ODE： \[ \left\{\begin{aligned} &y^{(m)}=f(t,y,y',...,y^{(m-1)}),\;a\leq t \leq b\\ &y(a)=\eta_1,\;y'(a)=\eta_2,...,y^{(m-1)}(a)=\eta_m \end{aligned}\right. \] 只要令：\(y_1=y,\;y_2=y',...,y_m=y^{(m-1)}\)，就化为了方程组： \[ \left\{\begin{aligned} &y_1'=y_2,\;y_1(a)=\eta_1\\ &y_2'=y_3,\;y_2(a)=\eta_2\\ &...\\ &y_m'=f(t,y_1,y_2,...,y_m),\;y_m(a)=\eta_m\\ \end{aligned}\right. \] 故无需多言
后记
- 这门课其实还讲了点 M-P 神经网络，然而感觉讲的太快太浅，偏科普性质（问就是感觉自己学了个寂寞），故没有整理至此。
- 没有整理差分方程一节，个人认为本质上跟 \(k\) 阶递推数列相关理论没什么区别，反正都是线代的活。
- 封面是我有一次把梯度下降方向搞反了搞成了“梯度上升”的图像捏，挺壮观的是不是（（（
