数值分析笔记

此篇完全是为了通过考试. 我不喜欢数值分析, 但是不得不学.

总的来说，这门课似乎包括以下几个板块

• 函数的逼近（多项式插值）

• 线性方程组的数值解（Gauss消元法，复矩阵的QR分解，QR方法和最小二乘）

• Fourier变换（离散/快速，三角多项式）

• 数值微分

• 数值积分（Newton-Cotes公式，Gauss公式和正交基）

• 非线性方程数值解

• 常微分方程的数值解（显式Euler法，隐式Euler法，改进显式Euler法和Runge-Kutta法）

函数的逼近

函数的逼近通过插值多项式来实现，给定点集及其上的函数值，输出是插值多项式，有两种表示插值多项式的办法

1. Newton差商表示：计算差商$\[x_0,\ldots,x_n\]f$，然后以这些差商作为系数，乘以(x − x₀)⋯(x − x_n − 1);

2. Langrange表示：考虑Lagrange多项式l_i，然后直接用f(x_i)l_i(x)的求和作为插值多项式。

逼近问题的好坏由收敛速度来刻画，收敛速度的定义是：当点集的最大区间长度h趋近于0时，绝对误差的收敛阶数h^p。因此，如果希望得到收敛阶数，可以绘制loglog图看图像的斜率是多少。

另外还可以考虑在每个区间内，用样条函数来逼近。区间和区间之间的光滑性是给定的m，区间内部都是某个次数的、用等距点来插值的多项式。用得比较多的有两个情形：常数样条函数、线性样条函数。线性样条函数的基本单位是b_i，它是在第i个区间内从0到1，在第i+1个区间内从1到0，在其它区间上为0的分段线性函数。利用这个基本单位函数，一般函数f的线性样条函数是$$f_\mathrm{int}=\sum_{i=0}^N f(x_i)b_i(x).$$

另外，用连续线性样条函数来插值，其误差还和节点的选取有关：如果节点选取得比较好，那么误差可能比较小。那么我们可以归纳地逐渐增加节点，使得每一步增加节点时，误差都变化比较小，最终在第N步时停止，这种办法叫做自适应细化（adaptive refinement）。考虑某个小区间τ，定义其中的变化为d_τ = max |u_n(M_τ) − u_n − 1(M_τ)|，其中M_τ是区间τ的中间点。假设这些变化的最大值是D，那么对于给定α > 0，可以把所有满足d_τ ≥ αD的区间中点都收集起来，和原来就有的区间端点一起，共同构成新的支撑点集合。

线性方程组的数值解

考虑复数域内的线性方程组Ax = b，希望将A作QR分解，然后就可以通过左乘$Q\*$来抵消Q，得到方程Rx = y，最后从后往前逐渐求解。这里Q是一个酉矩阵，而一般的消元得到的矩阵是下三角矩阵。因此，我们考虑Householder变换：逐步计算$P_j=I-2ww^\*$，其中P_jA的第一列只有第一个元素不为0. 计算次数是4/3n³ + O(n²)

当然，手算QR分解依旧可以用Gram-Schmidt正交化的思路，逐个计算。这里注意复数矩阵和（列）向量的内积是自己和共轭转置的乘积。

快速Fourier变换

1. 离散傅里叶变换

首先, 我们考虑一般的离散傅里叶变换: 一般地, 复变函数$f:\C\rightarrow \C$的傅里叶变换为
$$ F(t):=\int_\R f(s)\exp(-is\cdot t)\mathrm{d}s $$
其反变换为
$$ f(t)=\frac{1}{2\pi} \int_\R F(s)\exp(is\cdot t)\mathrm{d}t $$
给定间隔T_s, 如果我们接受Dirac函数$\delta:\R\rightarrow \R$的存在性, 那么f的离散化 (采样函数) 可以用Dirac函数表示为
$$ f_s(t)=\sum_{n=-\infty}^\infty f(t)\delta(t-nT_s) $$
于是f_s的Fourier变换 (代入定义公式) 就是
$$ F_s(t)= \sum_{n=-\infty}^\infty f(nT_s)\exp(-inT_s\cdot t). $$
然而, 在现实生活中, 如果我们并不知道f本身的表达式, 我们没有办法精确计算出无限多个f(nT_s). 所以我们的思路是: 在有限多个等距的点0, T_s, 2T_s, …, (N − 1)T_s上进行采样, 从而得到离散的Fourier变换
$$ F\[k\]=\sum_{n=0}^{N-1} f(nT_s)\exp(-\frac{2k\pi i}{N}\cdot n),\quad 0\leq k\leq N-1. $$
这个算法的复杂度是O(N²).

2. 快速Fourier变换

从离散Fourier变换中, 我们事实上

• 没有用到求和中单位根的性质;

• 指数函数的周期性;

于是, 为了减少算法的复杂性, 我们希望将这些性质用上, 最终将算法复杂度减少到O(Nlog N). 注意到离散Fourier变换中, 对于每一个k, Fourier系数$F\[k\]$都是一个N项的指数求和. 首先, 我们固定k, 记
$$ f\[n\]:=f(nT_s),\quad 0\leq n\leq N-1. $$
假设我们已经归纳地有了计算: 给定间隔T_s和采样值$(F\[0\],F\[1\],\ldots,F\[N-1\])$ 的”黑箱”
$$ \mathrm{FFT}_k(T_s; (F\[0\],F\[1\],\ldots,F\[N-1\]) ) \in \C, \quad 0\leq k\leq N-1. $$

1. 如果N​是偶数, 那么可以将求和分拆成奇数项求和与偶数项求和

$$\begin{aligned} F\[k\] &=\sum^{N/2-1}_{m=0} f\[2m\]\exp(-\frac{2k\pi i}{N}\cdot 2m) +\sum^{N/2-1}_{m=0} f\[2m+1\]\exp(-\frac{2k\pi i}{N}\cdot (2m+1))\\\\ &=\mathrm{FFT}_k(2T_s; x\[0\],x\[2\],\ldots,x\[N-2\])+\xi_N^{-k}\mathrm{FFT}_k(2T_s; x\[1\],x\[3\],\ldots, x\[N-1\])\\\\ &=:O\[k\]+\xi_N^{-k} E\[k\]. \end{aligned}$$

现在改变k: 注意到$O\[k\]=O\[k+\frac{N}{2}\]$, $E\[k\]=E\[k+\frac{N}{2}\]$, 于是由N次单位根的性质
$$F\[k+\frac{N}{2}\]=O(k)-\xi_N^k E\[k\].$$

2. 如果N是奇数, 那么可以将首项忽略, 其余项数都是ξ_N^−k​的倍数. 重新提取出来即可, 和偶数项只差一次加法和一次乘法.

例. 假设我们给定了8个采样, 需要计算其Fourier系数$F\[0\],\ldots, F\[7\]$​​, 那么我们只需要知道
E(0), …, E(3), O(0), …, O(3)
共计8个”四采样点FFT” 的计算, 再有2N = 16次加法和乘法组合. 四采样点FFT又可以变成4个二采样点+8次加法和乘法组合. 而每个二采样点都需要计算2次, 于是每个四采样点计算16次,

---

现在来计算上面这个递归算法的复杂度: 假设N个点的复杂度是T(N), 那么
$$ T(N)\sim2T(\frac{N}{2})+N\leq \cdots\sim N\sum\frac 1 n\sim N\log N. $$

数值微分

计算函数的数值微分可以用插值多项式的微分来实现。特别地，如果我们要计算微分在x₀处的取值，我们可以直接取x₀为第一个插值点，然后后面的微分计算就会更简便：因为Newton插值多项式中每一项都包含x − x₀，求导在x₀处计算更为简便。这时，误差是‖f^(n + 1)‖_∞ ⋅ |I|ⁿ.

数值积分

给定函数f和权重函数w，希望通过函数f在某些点上的取值，逼近积分∫wf。首先定义某个求积公式的误差阶数为：最大的保证多项式积分不变的多项式阶数。比如某个求积公式对次数 ≤ 3的多项式都精确，但对多项式次数=4时不成立，那么这个求积公式就是3次的。一般而言，求积公式具有下面的结构

$$I(f)=\sum_{j=1}^n \lambda_j f(x_j).$$

有一些常见的求积公式列在下面

1. （Riemann公式，0阶）f(x_j)(x_j + 1 − x_j)

2. （Midpoint公式，1阶）$$f(\frac{x_j+x_{j+1}}{2})(x_{j+1}-x_j)$$

3. （Trapezoidal公式，1阶）$$\frac{f(x_j)+f(x_{j+1})}{2}\cdot (x_{j+1}-x_j)$$

4. （Simpson公式，3阶）$$\frac 1 6 (f(x_j)+4f(\frac{x_j+x_{j+1}}{2}) +f(x_{j+1}))(x_{j+1}-x_j)$$

事实上，如果节点取定，那么有一种唯一的方式选取参数λ_j，使得阶数最大。

• Newton-Cotes公式：等距节点，w(x) = 1

• Gauss公式：希望能寻找合适的节点，使得其构造的积分公式误差最小。给定插值点数量n，这n个点应该选取成“n次正交多项式”的基底，其中这里正交的含义是和w有关的，因为内积是∫(fg)(x)w(x)dx.

  • 如果w = 1，正交多项式就是Legendre多项式$$\pi_n(x)=\sqrt{n+\frac 1 2+\frac{1}{2^n n!} \frac{\mathrm{d}^n}{\mathrm{d}x^n} ((x^2-1)^n).$$

  • 如果$w=\frac{1}{\sqrt{1-x^2}}$，正交多项式就是Chebyshev多项式Q_n(x) = cos (narccos x).

数值求解非线性方程

• 二分法：无需赘述，注意误差为ϵ的计算次数；

• Newton法：线性化+迭代，先选择初值，然后用线性化逼近f，x_n + 1就是线性化f后的线性方程的解。该方法依赖于初值的选取：如果初值的选择不好，可能方法不收敛。为此，我们为了保证收敛性，需要首先
  1. 通过二分法找到一个解存在的区间；
  2. 选择这个区间，使得区间长度充分小，并且在这个区间上，f^′的下界和f^″的上界也有控制；

具体而言，下面的定理和推论保证的：

然而，在实际应用中，我们不可能一开始就知道零点在哪，从而也没办法确定以零点为中心的邻域。为此我们只能扩大区间长度：

利用这个推论，在实际操作中就可以用牛顿法来确定零点：

1. 找到一个区间，使得零点在这个区间内；

2. 算出这个大区间的M；

3. 通过二分法缩小大区间，使得4δM < 1, 其中3δ是小区间的区间半径；

4. 任选小区间内的点作为初值点，进行牛顿法。

• 定点法：将f(x) = 0转换为适当的不动点方程x = φ(x)，然后用φ来迭代。事实上，牛顿法也是一种定点法，它将方程转化为了$x=x-\frac{f(x)}{f'(x)}$.

数值求解ODE

一般求解实数域上的一阶ODE。如果右侧的函数f是Lipschitz连续的，那么解存在且唯一，且连续依赖于初值x₀, y₀。有了这个理论保证，可以考虑单步法的数值求解。一个单步法由函数\[:[a,b]^+\]决定：该函数给出的迭代是\[y_{}=y+h(x,y,h).\]假设ODE的精确解是κ，那么单步法的局部离散误差定义为\[] 如果局部离散误差在h趋近于0时，能够一致收敛到0，那么就称这个单步法是一致的。这个局部离散误差一致收敛到0的阶数就是单步法的阶数。

根据选择函数Φ的不同，有如下几种不同的单步法

1. 显式欧拉法：Φ(x, y, h) = f(x, y);

2. 隐式欧拉法：\[(x,y,h)=f(x+h,y(x+h))=f(x+h,y_{})\]

3. Crank-Nicolson方法：\[(x,y,h)= 2 (f(x,y)+f(x+h,y(x+h)) \] 或者\[(x,y,h)= 2 (f(x,y)+f(x+h, x+hf(x,y) ) \] 其中，第一个是隐式的，第二个是显式的；

4. 改进欧拉法（2阶）：令k₁ = f(x, y)，k₂ = f(x + h/2, y + 1/2hk₁))，然后定义\[(x,y,h)=k_2\] 就是先用显式欧拉法走半步，再用走的半步走到的点对应的f作为“速度”，从原地走一步；

5. Runge-Kutta法：事实上是改进欧拉法的推广——走两次推广为走s次，定义\[k_1=f(x,y),\]\[k_s=f(x+s h,y+h{j=1}^{s-1} _{s,j} k_j)\]然后定义\[y_{} = y+h_{s=1}^r _s k_s\]并且取参数∑α_s = 1来保证一致性.

---

给定闭区间$I=\[a,b\]$, 定义在I上的函数f, 以及I的一个分割\[{x_0, x_1,, x_n},\]其中x₀, …, x_n是互异的. Lagrangian多项式所做到的, 就是给出一个n次的多项式, 使得该多项式在这些给定的点{x_i}处和f相等. 关于这样的多项式, 我们有如下唯一性定理:
如果记((x)=_i(x-x_i)), 那么(L_i(x)=).
任何领域的experts都会注意到, 很难直接用定义的方式来计算一个对象. 在实际应用中, 一般用归纳法计算Lagrangian插值多项式:

这样一来,
\[]
而对于系数$f\[x_0,\ldots,x_k\]$, 有如下Newton公式归纳
\[\]
从Newton公式知, $f\[x_0,\ldots,x_k\]$不依赖于x₀, …, x_k的顺序.

在应用中, 我们会采取无穷范数来进行函数的误差估计. 给定插值点Θ_n = {x₀, …, x_n}, 这些插值点确定函数ω_n(x) = (x − x₀)⋯(x − x_n), 和n次的牛顿多项式p_n. 那么函数f和p_n的误差满足
\[f(x)-p_n(x)=_n().\]
特别地, 数值上说,
\[,\]
其中C_n := ‖f^(n + 1)‖_∞.

如果函数的光滑性不好, 那么我们可以通过缩短区间长度来估计, 这样得到的分段函数称为样条函数 (spline function). 给定区间$I=\[a,b\]$的支撑点集
\[_n= {a=x_0<x_1<<x_n=b}, \]

定义(h_i=x_i-x_{i-1}, h=_{1in}h_i), 以及所有支撑点构成的集合为
\[:={ [x_{i-1},x_i] : i=1,2,, n}. \]
给定光滑度k, 我们定义所有样条函数构成的空间为
\[{}^{k,m} := {uC^k(I) | , u|_m }. \]

特别地, 如果k > m, 即全局光滑度大于多项式次数, 则样条函数空间就是全局的多项式函数空间, 即𝒮_𝒢^k, m = ℙ_m. 因此, 在下面的讨论中, 我们总是假设k ≤ m.