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

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

• 函数的逼近（多项式插值）
• 线性方程组的数值解（Gauss消元法，复矩阵的QR分解，QR方法和最小二乘）
• Fourier变换（离散/快速，三角多项式）
• 数值微分
• 数值积分（Newton-Cotes公式，Gauss公式和正交基）
• 非线性方程数值解
• 常微分方程的数值解（显式Euler法，隐式Euler法，改进显式Euler法和Runge-Kutta法）

函数的逼近

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

1. Newton差商表示：计算差商$[x_0,\ldots,x_n]f$，然后以这些差商作为系数，乘以$(x-x_0)\cdots(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）。考虑某个小区间$\tau$，定义其中的变化为$d_\tau=\max |u_n(M_\tau)-u_{n-1}(M_\tau)|$，其中$M_\tau$是区间$\tau$的中间点。假设这些变化的最大值是$D$，那么对于给定$\alpha>0$，可以把所有满足$d_\tau\geq\alpha D$的区间中点都收集起来，和原来就有的区间端点一起，共同构成新的支撑点集合。

线性方程组的数值解

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

当然，手算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,\ldots,(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)$.

2. 快速Fourier变换

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

• 没有用到求和中单位根的性质;
• 指数函数的周期性;

于是, 为了减少算法的复杂性, 我们希望将这些性质用上, 最终将算法复杂度减少到$O(N\log 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$是奇数, 那么可以将首项忽略, 其余项数都是$\xi^{-k}_N$​的倍数. 重新提取出来即可, 和偶数项只差一次加法和一次乘法.

例. 假设我们给定了8个采样, 需要计算其Fourier系数$F[0],\ldots, F[7]$​​, 那么我们只需要知道
$$
E(0),\ldots,E(3), O(0),\ldots,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_0$处的取值，我们可以直接取$x_0$为第一个插值点，然后后面的微分计算就会更简便：因为Newton插值多项式中每一项都包含$x-x_0$，求导在$x_0$处计算更为简便。这时，误差是$\Vert f^{(n+1)}\Vert_\infty\cdot |I|^n$.

数值积分

给定函数$f$和权重函数$w$，希望通过函数$f$在某些点上的取值，逼近积分$\int wf$。首先定义某个求积公式的误差阶数为：最大的保证多项式积分不变的多项式阶数。比如某个求积公式对次数$\leq 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)$$

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

• Newton-Cotes公式：等距节点，$w(x)=1$
• Gauss公式：希望能寻找合适的节点，使得其构造的积分公式误差最小。给定插值点数量n，这n个点应该选取成“n次正交多项式”的基底，其中这里正交的含义是和$w$有关的，因为内积是$\int (fg)(x) w(x)\mathrm{d}x$.
  • 如果$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(n\arccos x).$$

数值求解非线性方程

• 二分法：无需赘述，注意误差为$\epsilon$的计算次数；
• Newton法：线性化+迭代，先选择初值，然后用线性化逼近$f$，$x_{n+1}$就是线性化f后的线性方程的解。该方法依赖于初值的选取：如果初值的选择不好，可能方法不收敛。为此，我们为了保证收敛性，需要首先
  1. 通过二分法找到一个解存在的区间；
  2. 选择这个区间，使得区间长度充分小，并且在这个区间上，$f’$的下界和$f”$的上界也有控制；

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

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

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

1. 找到一个区间，使得零点在这个区间内；
2. 算出这个大区间的M；
3. 通过二分法缩小大区间，使得$4\delta M<1$, 其中$3\delta$是小区间的区间半径；
4. 任选小区间内的点作为初值点，进行牛顿法。

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

数值求解ODE

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

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

1. 显式欧拉法：$\Phi(x,y,h)=f(x,y)$;
2. 隐式欧拉法：\[\Phi(x,y,h)=f(x+h,y(x+h))=f(x+h,y_{\mathrm{next}})\]
3. Crank-Nicolson方法：\[\Phi(x,y,h)=\frac 1 2 (f(x,y)+f(x+h,y(x+h)) \] 或者\[\Phi(x,y,h)=\frac 1 2 (f(x,y)+f(x+h, x+hf(x,y) ) \] 其中，第一个是隐式的，第二个是显式的；
4. 改进欧拉法（2阶）：令$k_1=f(x,y)$，$k_2=f(x+h/2, y+1/2 hk_1))$，然后定义\[\Phi(x,y,h)=k_2\] 就是先用显式欧拉法走半步，再用走的半步走到的点对应的$f$作为“速度”，从原地走一步；
5. Runge-Kutta法：事实上是改进欧拉法的推广——走两次推广为走s次，定义\[k_1=f(x,y),\]\[k_s=f(x+\mu_s h,y+h\sum_{j=1}^{s-1} \mu_{s,j} k_j)\]然后定义\[y_{\mathrm{next}} = y+h\sum_{s=1}^r \alpha_s k_s\]并且取参数$\sum \alpha_s=1$来保证一致性.

---