这学期选了一个数值分析。说老实话，是实在找不到特别有意思的课，稍微有点凑学分拔将军的意味。不过这课教的知识的确非常有用，所以最好还是写一下笔记，以备以后不时之需。这课用的教材是heath的scientific computing。这笔记基本上就是抄书。

数值分析笔记

观念比原理更重要，原理比方法更重要。

第一章

第一章是对于科学计算中重要概念的介绍，没有提到具体的算法。从中学到了一些基本的考虑问题的抓手。

  • 科学计算的核心就是误差和误差分析。误差大致有两类，一种是“软性”起源的，比如说解析函数需要做taylor展开，在我们需要的范围内留前几项，或者把非线性方程近似成线性来解，或者解稀疏矩阵，等等。另一种是“硬性”起源的——我这么形容的时候，想的是软件和硬件的区别——它指的是计算机进行一般运算时的精度问题。

    比如说，求导。对于一个小量h，$f^{\prime}(x) \approx \frac{f(x+h)-f(x)}{h}$就是不考虑二阶导的结果。考虑了呢？会多一项$h\cdot\frac{|f''(z)|}{2}\propto h$。这就是“软性”起源的truncation error。另一方面，实际计算$\frac{f(x+h)-f(x)}{h}$的时候，由于计算机的精度最多只能到$\epsilon>0$，这会导致这个计算多了$2\epsilon/h\propto h^{-1}$的误差，这是“硬性”起源的rounding error。

    有一些时候，计算过程中有许多步。我们担心误差会被步步放大。定量描述这种现象的概念是conditioning number $\frac{|\Delta y / y|}{|\Delta x / x|}$。还有forward/backward error的概念，分别指的是因变量的误差和其对应的自变量的误差，我们也可以把conditioning number理解为relative fwd error和relative bkwd error之比：$|(\Delta xf'(x)/f(x))/(\Delta x/x)|$，或$|xf'(x)/f(x)|$。

    计算机喜欢处理连续而不过分陡峭的数学问题。如果不连续，就斥之为ill-posed；连续但过分陡峭的，斥之为ill-conditioned。

  • 理解Rounding error，或计算机精度导致的误差，需要先理解计算机的浮点数机制。浮点数的表示公式如下：

    $x=\pm\left(d_{0}+\frac{d_{1}}{\beta}+\frac{d_{2}}{\beta^{2}}+\cdots+\frac{d_{p-1}}{\beta^{p-1}}\right) \beta^{E}$

    其中，beta是进制，E是$[L, U]$范围内的一个数，描述这个浮点数的数量级；p是精度，所以括号内的部分代表有效数字。一般浮点数都是要normalize的，意思就是说，除非x=0，否则$d_0$不能是0。这不难理解，我们只要考虑科学计数法，并与上面的公式作比较，就能明白。

    Remark: 我居然才意识到科学计数法里面的“科学”居然指的是科学计算。我小时候一直以为是表示大数方便才这么规定，现在才发现，原来是计算上的需求。

    会有一个underflow level $\mathrm{UFL}=\beta^{L}$ 和一个overflow level $\mathrm{OFL}=\beta^{U+1}\left(1-\beta^{-p}\right)$。

    机器处理那些精度过高的数字的时候会给round到一定精度去，这会导致一个误差。比如，如果是四舍五入，那么可能出现的（在有效数字部分的）最大误差为$\epsilon_{\mathrm{mach}}=\frac{1}{2} \beta^{1-p}$。

    浮点数运算会造成一定的误差。一种情况是，数量级不同的数相加，较小的数会被损失掉很多信息。例如，对于上述$\epsilon_{mach}$，有$(1+\epsilon)+\epsilon=1$，而$1+(\epsilon+\epsilon)>1$。一般的，对于加减乘除来说，有$\mathrm{fl}(x$ op $y)=(x$ op $y)(1+\delta)$（fl指float，op指operation）。其中$|\delta| \leq \epsilon_{\mathrm{mach}}$限制了这部分rounding error的程度。这很好理解，并且可以直接用于误差分析。

    另一种情况是，首位数字及数量级相同的数相减，会放大这些数本来自带的误差，因为误差的绝对大小不变的情况下，产生误差的数却小了很多。这叫cancellation error。

第二章

第二章粗浅介绍了一下解线性系统的事情。

  • 引入了矩阵的范数来讨论误差分析的事情，因为矩阵的范数有性质$\|\boldsymbol{A} \boldsymbol{x}\| \leq\|\boldsymbol{A}\| \cdot\|\boldsymbol{x}\|$。例如，对于线性系统Ax=b，如果b出了一些偏差，那么x的误差将会被控制在$\frac{\|\Delta \boldsymbol{x}\|}{\|\boldsymbol{x}\|} \leq\left\|\boldsymbol{A}^{-1}\right\| \cdot\|\Delta \boldsymbol{b}\| \frac{\|\boldsymbol{A}\|}{\|\boldsymbol{b}\|}$的范围，我们说$\|\boldsymbol{A}\| \cdot\left\|\boldsymbol{A}^{-1}\right\|$是系统的conditioning number。

  • 引入了残差的概念：对于有偏差的解$\hat{\boldsymbol{x}}$，残差是$r=b-A \hat{x}$，但更合理的应该是scale过的“相对”残差，即$\|r\| /(\|A\| \cdot\|\hat{x}\|)$。同样，$\frac{\|\Delta \boldsymbol{x}\|}{\|\hat{\boldsymbol{x}}\|} \leq \operatorname{cond}(\boldsymbol{A}) \frac{\|\boldsymbol{r}\|}{\|\boldsymbol{A}\| \cdot\|\hat{\boldsymbol{x}}\|}$，所以残差与误差之间还隔着一个condition。

  • 高斯消元法和L-U分解，这是复习大一的内容了。回忆一下，我们把一个矩阵乘以许多的下三角矩阵（$L^{-1}$），它们各处理一个列，最终把待分解的矩阵变成一个上三角矩阵（$U$）。

    在做三角矩阵的时候，我们如果想把矩阵的某一列中，某一元素以上/下的元素都变成0，常规的办法是用这一元素与那些变零的元素来消去，这意味着这一元素（pivot）不能是0。实际计算中，应当避免太小的pivot，因为它的倒数太大，会造成严重的rounding误差。

    这一方法（LU分解）需要的时间复杂度是$\mathcal{O}(n^{3})$，但分解完之后求解LUx=b的时间复杂度则是$\mathcal{O}(n^{2})$，比求逆矩阵快，而且精度高。分解完了之后，对于不同的b都可以解；对于稍微变动（变了一个秩一矩阵）的A也可以用$\mathcal{O}(n^{2})$的Sherman-Morrison公式来处理，仍然可以解得快。

    对于正定对称阵，可以用Cholesky分解，这是说$\boldsymbol{A}=\boldsymbol{L} \boldsymbol{L}^{T}$。这种办法的时间和空间复杂度都约只有普通LU分解的一半。

第三章

第三章介绍了最小二乘法。

  • 在解决$A x=b$问题时，如果A不是方阵而使得问题overdetermined了，就会成为一个最小二乘法问题$A x \cong b$。人们这时希望求解使得残差$r=b-A x$在2-norm之下最小的x。这种情况是很常见的，一个例子就是多项式拟合，因为通常数据点的数目远多于多项式的次数。

  • 可以由$\phi(\boldsymbol{x})=\|\boldsymbol{r}\|_{2}^{2}=\boldsymbol{r}^{T} \boldsymbol{r}=(\boldsymbol{b}-\boldsymbol{A} \boldsymbol{x})^{T}(\boldsymbol{b}-\boldsymbol{A} \boldsymbol{x})=\boldsymbol{b}^{T} \boldsymbol{b}-2 \boldsymbol{x}^{T} \boldsymbol{A}^{T} \boldsymbol{b}+\boldsymbol{x}^{T} \boldsymbol{A}^{T} \boldsymbol{A} \boldsymbol{x}$求导解出最优的$\boldsymbol{A}^{T} \boldsymbol{A} \boldsymbol{x}=\boldsymbol{A}^{T} \boldsymbol{b}$。另有一种几何的办法，即注意到最优的情形下事实上对应将b投影到Ax张成的空间，$A x=P b$，这时在两边同乘$\boldsymbol{A}^{T}$可得到相同结果。

  • 在condition的时候，按理说可以使用伪逆来计算conditioning number，但又注意到，如果b几乎垂直于Ax张成的平面，那么b的微扰会极大干扰x的求解。为将这一事实纳入考虑，可用$\frac{\|\boldsymbol{A} \boldsymbol{x}\|_{2}}{\|\boldsymbol{b}\|_{2}}=\frac{\|\boldsymbol{y}\|_{2}}{\|\boldsymbol{b}\|_{2}}=\cos (\theta)$去bound这一因素，经计算可得

    $\frac{\|\Delta \boldsymbol{x}\|_{2}}{\|\boldsymbol{x}\|_{2}} \leq\operatorname{cond}_{previous}(\boldsymbol{A}) \frac{1}{\cos (\theta)} \frac{\|\Delta \boldsymbol{b}\|_{2}}{\|\boldsymbol{b}\|_{2}}$

    对A的微扰造成的影响复杂些，但精神是一样的。

  • 求解$\boldsymbol{A}^{T} \boldsymbol{A} \boldsymbol{x}=\boldsymbol{A}^{T} \boldsymbol{b}$有多种办法。如A满秩，那么$A^TA$可做Cholesky分解。但数值上还有诸多tricky之处，例如矩阵乘法带来的cancellation等。故人们也去探索别的方法，而非仅从此式求解最小二乘问题。

    • 例如可将问题embed到更大的问题中，使得解线性问题时的pivot更加友好。如，注意到残差垂直于span(A)，故可以写下

      $\left[\begin{array}{ll}\boldsymbol{I} & \boldsymbol{A} \\ \boldsymbol{A}^{T} & \boldsymbol{O}\end{array}\right]\left[\begin{array}{l}\boldsymbol{r} \\ \boldsymbol{x}\end{array}\right]=\left[\begin{array}{l}\boldsymbol{b} \\ \boldsymbol{o}\end{array}\right]$

      当然还得properly scale r和x，于是实际上是$\left[\begin{array}{ll}\alpha \boldsymbol{I} & \boldsymbol{A} \\ \boldsymbol{A}^{T} & \boldsymbol{O}\end{array}\right]\left[\begin{array}{c}\boldsymbol{r} / \alpha \\ \boldsymbol{x}\end{array}\right]=\left[\begin{array}{l}\boldsymbol{b} \\ \boldsymbol{o}\end{array}\right]$的形式，实验上一般凭经验取$\alpha=\max _{i, j}\left|a_{i j}\right| / 1000$即可。

    • 另外还有一种办法。若我们的A可以变成上三角矩阵，那么可以将问题划分成$\left[\begin{array}{l}\boldsymbol{R} \\ \boldsymbol{O}\end{array}\right] \boldsymbol{x} \cong\left[\begin{array}{l}c_{1} \\ c_{2}\end{array}\right]$，这个时候实际上可以求得$\|\boldsymbol{r}\|_{2}^{2}=\left\|\boldsymbol{c}_{2}\right\|_{2}^{2}$。于是只需要考虑问题的上半部分就可以上一章的办法来解了。为达到这个情况所采用的的手段叫做QR分解。就是把A写成$\boldsymbol{A}=\boldsymbol{Q}\left[\begin{array}{l}\boldsymbol{R} \\ \boldsymbol{O}\end{array}\right]$的形式。

  • QR分解当然可以使用前述LU分解的类似手段。但是我们希望能够保A的2-范数，这个时候应当做一些正交变换。有许多种这样的变换。

    • householder（这是个人名）变换。首先，设想存在一个正交矩阵叫做householder矩阵，它作用在一个向量a上，得到的a’只有前n个元素非零。那么我们就可以用一连串的householder矩阵将A造成一个上三角矩阵。

      实际上，将变换的头尾连起来，若它们的长度相等，那么过此连线v的中点做一垂直平面，则可认为变换前后的向量是以此平面为镜面做的反射变换的结果。故我们算出a的模长，得到其连线v的方向$a-\alpha a'$。

      用$\boldsymbol{P}=\boldsymbol{v}\left(\boldsymbol{v}^{T} \boldsymbol{v}\right)^{-1} \boldsymbol{v}^{T}=\left(\boldsymbol{v} \boldsymbol{v}^{T}\right) /\left(\boldsymbol{v}^{T} \boldsymbol{v}\right)$可以将a投影到v的方向。故用$I-P$可投影到上述垂直平面，而$H=I-2 P=I-2\left(v v^{T}\right) /\left(v^{T} v\right)$则意味着投影到a’。H当然也确实是正交的，这就是householder矩阵。

    • givens（这也是个人名）旋转。非常简单，就是取两个分量，拉出来做一个平面旋转， 然后把其中一个转到0。例如

      $\left[\begin{array}{rrrrr}1 & 0 & 0 & 0 & 0 \\ 0 & c & 0 & s & 0 \\ 0 & 0 & 1 & 0 & 0 \\ 0 & -s & 0 & c & 0 \\ 0 & 0 & 0 & 0 & 1\end{array}\right]\left[\begin{array}{c}a_{1} \\ a_{2} \\ a_{3} \\ a_{4} \\ a_{5}\end{array}\right]=\left[\begin{array}{c}a_{1} \\ \alpha \\ a_{3} \\ 0 \\ a_{5}\end{array}\right]$

    • Gram-Schimidt正交化。这个也是很久以前学过的，把下一个向量投影到已经处理好的部分张成的空间，然后把这个投影砍掉，只留下垂直的部分，然后scale。

    • 如果A不满秩，我们在正交化的每一步应该用permutation矩阵把范数大的列挪到前面来。这叫column pivoting。这样最后就得到$Q^{T} A P=\left[\begin{array}{ll}R & S \\ O & O\end{array}\right]$，其中R是一个上三角的方阵，并且满秩。

  • SVD也是一个经典方法。

  • 不加说明地列举一下各个方法的性能。

    若A为m*n矩阵，那么在$m \approx n$时，householder QR不如最普通的办法（$\boldsymbol{A}^{T} \boldsymbol{A} \boldsymbol{x}=\boldsymbol{A}^{T} \boldsymbol{b}$）。$m \gg n$时，普通法更优。但householder法的误差基本永远更小。

    cholesky分解在$\operatorname{cond}(\boldsymbol{A}) \approx 1 / \sqrt{\epsilon_{\mathrm{mach}}}$时不可用，而householder法仅在$\operatorname{cond}(\boldsymbol{A}) \approx 1 / \epsilon_{\mathrm{mach}}$时才不可用。

    SVD是最稳定，同时也最慢的。

第四章

特征值和与之有关的性质刻画了线性算子的几何性质。这一章就来讲一讲特征值问题的数值方法。

  • 一些评论

    • 数学上，需要意识到，由于五次以上的方程没有根式解，大矩阵的特征值也很难去精确求解。不过，我们有收敛很快的近似方法。
    • 我们可以大致划定特征值的范围。Gershgorin定理说明本征值被限制在一系列圆盘里面，其圆心为对角元，而半径则是这一行其他元素的模长之和。

  • 问题的敏感性

    • 取方阵A，其有本征值$\lambda_i$，本征向量组X。对A做微扰，变为$A+E$，有本征值$\mu$，则（可以证明）：

      $\left|\mu-\lambda_{k}\right|\le\operatorname{cond}_{2}(\boldsymbol{X})\|\boldsymbol{E}\|_{2}$。

      因此，本征向量的性质决定conditioning的好坏。本征向量越是倾向于两两垂直，则conditioning就越好。正交矩阵也因此总是well-conditioned的。对于接近不满秩的矩阵，则本征向量接近于线性相关，而conditioning不好。

  • 问题的清洗

    进行预处理，尽量不改变本征值，而将问题转化为可解的形式。一种最普遍的办法是利用相似形$B=T^{-1} A T$，其中T满秩。则：

    $B y=\lambda y \quad \Rightarrow \quad T^{-1} A T y=\lambda y \quad \Rightarrow \quad A T y=\lambda T y$。

    例如，你可以（通过一些相似变换）把矩阵变成General Schur form

    $\boldsymbol{A}=\left[\begin{array}{cccc}\boldsymbol{A}_{11} & \boldsymbol{A}_{12} & \cdots & \boldsymbol{A}_{1 p} \\ & \boldsymbol{A}_{22} & \cdots & \boldsymbol{A}_{2 p} \\ & & \ddots & \vdots \\ & & & \boldsymbol{A}_{p p}\end{array}\right]$

    其中$\lambda(\boldsymbol{A})=\bigcup_{j=1}^{p} \lambda\left(\boldsymbol{A}_{j j}\right)$。

    这基本上是最坏的情况。如果好一些，A是实的，那么该form的对角元要么是实的本征值，要么是2*2的block，对应一组复的本征值。如果还要好一些（但不是同一个意义上的好），就是说，如果A是满秩的，那么对于相似变换$\boldsymbol{B}=\boldsymbol{T}^{-1} \boldsymbol{A T}$，可以让B是diag$(\lambda_i)$，T是本征向量。如果再更好，那么B和T的性质也可以更好。

    对于下面的一些办法，例如QR迭代，应该尽量把矩阵洗干净，弄成尽量tridiagonal（原矩阵对称的情况）或hessenberg（原矩阵不对称的情况）的形式。这样做，可以把每次QR的花销从$\mathcal{O}\left(n^{3}\right)$ 降到 $\mathcal{O}\left(n^{2}\right)$，并且减少QR的次数。总的来说，QR迭代是一个$\mathcal{O}\left(n^{3}\right)$的算法。

  • 计算方法。

    • Power iteration。任取向量$x_{0}=\sum_{j=1}^{n} \alpha_{j} v_{j}$，那么乘以A（并重新归一化）许多次，肯定会迭代到对应最大本征值的那个本征向量。当然，有可能会失败，比如$\alpha_0$恰巧为0，或者最大的本征值是简并的，等等。
      • 换成A的逆，这么做就可以知道最小的本征值。
      • 平移一下，弄成$\boldsymbol{A}-\sigma \boldsymbol{I}$（的逆），就可以知道离$\sigma$最近的本征值。（这就是inverse iteration强的地方）恰当的平移可以在任何情况下改进最大和次大的本征值的比例。

  • 瑞利商。取一个x，当成是近似的本征向量，可以像求解最小二乘问题那样求解本征值：$\boldsymbol{x} \lambda \cong \boldsymbol{A} \boldsymbol{x}$。结果是$\lambda=\frac{x^{T} A x}{x^{T} x}$，这被称为瑞利商。在用A的逆去迭代求解时，你可以先得到一个迭代的结果x，算出瑞利商$\sigma$，然后平移。

  • 算出一个本征值之后，就像算出多项式的根那样，你可以除以（x-a），然后去求解别的根。

    • 一种做法是，用Householder变换对第一个本征向量得出$H x_{1}=\alpha e_{1}$，然后以$H A H^{-1}=\left[\begin{array}{cc}\lambda_{1} & \boldsymbol{b}^{T} \\ \mathbf{0} & \boldsymbol{B}\end{array}\right]$得到B，B的本征值就是接下来要求的了。但它的本征向量和A不一样。求出B的本征向量$y_2$后，求出$\alpha=\frac{\boldsymbol{b}^{T} \boldsymbol{y}_{2}}{\lambda_{2}-\lambda_{1}}$，再求出A的本征向量$x_{2}=H^{-1}\left[\begin{array}{l}\alpha \\ y_{2}\end{array}\right]$。
    • 或者，取$\boldsymbol{u}_{1}^{T} \boldsymbol{x}_{1}=\lambda_{1}$，然后$\boldsymbol{A}-\boldsymbol{x}_{1} \boldsymbol{u}_{1}^{T}$。我觉得这样更自然。

  • 还有一些迭代法。比如，从$n \times p$的$\boldsymbol{X}_{0}$开始，一次迭代p个向量，$\boldsymbol{X}_{k}=\boldsymbol{A} \boldsymbol{X}_{k-1}$，你将会得到最大的p个本征值所对应的那些本征向量张成的空间。但是，你会总是需要归一化你的各个向量。一个简单的办法是每轮都做QR分解，然后从单位矩阵$X_{0}=I$开始，用

    $\begin{aligned} \hat{\boldsymbol{Q}}_{k} \boldsymbol{R}_{k} &=\boldsymbol{X}_{k-1} \\ \boldsymbol{X}_{k} &=\boldsymbol{A} \hat{\boldsymbol{Q}}_{k} \end{aligned}$

    来取代之前的过程。我们可以把这个认为是每轮更新Q，使得$\boldsymbol{A}_{k}=\hat{\boldsymbol{Q}}_{k}^{H} \boldsymbol{A} \hat{\boldsymbol{Q}}_{k}$是一个三角之类容易算本征值的形式（$H\equiv\dagger$）。同样，也可以在这个里面引入平移来加速整个计算。从$\boldsymbol{A}_{0}=\boldsymbol{A}$开始，用

    $\boldsymbol{Q}_{k} \boldsymbol{R}_{k}=\boldsymbol{A}_{k-1}-\sigma_{k} \boldsymbol{I}$
    $\boldsymbol{A}_{k}=\boldsymbol{R}_{k} \boldsymbol{Q}_{k}+\sigma_{k} \boldsymbol{I}$

    来迭代。这么做能行，是因为在迭代$\boldsymbol{A}^{k}=\hat{\boldsymbol{Q}}_{k} \hat{\boldsymbol{R}}_{k}$的同时，你实际上也迭代了$A^\dagger$的逆：$\left(\boldsymbol{A}^{-H}\right)^{k}=\hat{\boldsymbol{Q}}_{k} \hat{\boldsymbol{R}}_{k}^{-H}$（因为Q是厄密的）。实用上说，$\boldsymbol{A}_{k-1}$的右下角取一个2*2矩阵，它的本征值可以作为对平移距离$\sigma_k$（或对最大本征值）的估计。或者更粗略，直接取右下角的元素$a_{n n}^{(k-1)}$。

  • Krylov子空间方法。首先，将Power iteration的各步作为列向量排成一个矩阵：

    $\boldsymbol{K}_{k}=\left[\begin{array}{llllll}\boldsymbol{x}_{0} & \boldsymbol{x}_{1} & \cdots & \boldsymbol{x}_{k-1}\end{array}\right]=\left[\begin{array}{lllll}\boldsymbol{x}_{0} & \boldsymbol{A} \boldsymbol{x}_{0} & \cdots & \boldsymbol{A}^{k-1} \boldsymbol{x}_{0}\end{array}\right]$。

    该序列的元素张成Krylov子空间$\mathcal{K}_{k}=\operatorname{span}\left(\boldsymbol{K}_{k}\right)$。此时，我们发现

    $\begin{array}{llllll}\boldsymbol{A} \boldsymbol{K}_{n} & = & {\left[\begin{array}{llll}\boldsymbol{A} x_{0} & \cdots & \boldsymbol{A} \boldsymbol{x}_{n-2} & \boldsymbol{A} \boldsymbol{x}_{n-1}\end{array}\right]}\end{array}=\boldsymbol{K}_{n}\left[\boldsymbol{e}_{2} \quad \cdots \quad \boldsymbol{e}_{n} \quad \boldsymbol{a}\right] \equiv \boldsymbol{K}_{n} \boldsymbol{C}_{n}$

    而$\boldsymbol{a}=\boldsymbol{K}_{n}^{-1} \boldsymbol{x}_{n}$。注意到这里C是一个Hessenberg矩阵。这样，实际上就是$K_{n}^{-1} A K_{n}=C_{n}$。进一步，由于Kn不是一组好的基，我们应当做QR分解，$Q_{n} R_{n}=K_{n}$，这样仍然是$\boldsymbol{Q}_{n}^{H} \boldsymbol{A} \boldsymbol{Q}_{n}=\boldsymbol{R}_{n} \boldsymbol{C}_{n} \boldsymbol{R}_{n}^{-1} \equiv \boldsymbol{H}$，其中H是个hessenberg矩阵。

    我们使用Arnoldi迭代来计算Q的列向量q，并用q得到H，用H得到A的本征值和本征向量。H的本征值叫做Ritz values，它们近似等于A的本征值。用Qk乘以H的本征向量（Ritz向量），就得到近似的A的本征向量。计算q们所使用的办法就是一边乘以A，一边Gram-Schimidt正交化。这种办法每次需要动用此前算出的所有q，因此时间和内存上都比较贵。

    如果待分解的矩阵是对称或厄密的，那么计算起来会容易很多。我们使用Lanczos迭代（这就是个名字）来计算。为什么？因为这个时候H是tridiagonal的，我们并不需要每次正交化q的时候都让它和之前的所有q都垂直，而只需要让它和它之前的两个q垂直，剩下的就自动垂直。所以我们大量减少了计算量和内存。

    • 雅克比方法。这是一种计算对阵矩阵本征值的办法。它通过平面旋转来一次消去待求解矩阵的对称的两个元素。但是这种方法比较慢，而且也很难推广到非对称的矩阵。

第五章

这一章的主题是非线性问题。这样的问题可以归结为求一个非线性函数或函数组的零点。

  • 解的存在性和唯一性。一般来说，这是比较复杂的。一般我们在一个区间或区域内求解，这叫做bracket，而在bracket上使用中值定理、反函数定理可以得到一些关于解的数学性质（比如存在性）。但是总的来说，似乎没有什么好的办法。
  • 敏感性。absolute condition number，可以想见，是雅可比矩阵决定的。在某根$x\*$处的此数值为$\left\|\boldsymbol{J}_{f}^{-1}\left(\boldsymbol{x}^\*\right)\right\|$。
  • 迭代方法。
    • 人们用指数记录迭代方法收敛的速度，例如，线性收敛的convergence rate是1。
    • 二分法。它当然是线性收敛的。还有不动点法，一般来说，它的收敛速度是和泰勒展开有关的，如果在不动点处$g^{\prime}\left(x^{*}\right)=0$，那么就肯定至少是二次以上的快速收敛。
    • 牛顿法。这就是说，$x_{k+1}=x_{k}-\frac{f\left(x_{k}\right)}{f^{\prime}\left(x_{k}\right)}$，或$x_{k+1}=x_{k}-J\left(x_{k}\right)^{-1} f\left(x_{k}\right)$。这相当于解$g(x)=x-f(x) / f^{\prime}(x)$的不动点，可以知道，一般情况下，牛顿法至少是二次的收敛（但要离根够近，否则不收敛）。计算雅克比矩阵的成本是$n^{2}$，n是维数，而解线性系统的成本（回忆两章前）是$\mathcal{O}\left(n^{3}\right)$。
    • 割线法。牛顿法有个问题，就是每次都要在下一个点重新算一遍导数的值。割线法就是用$f^{\prime}\left(x_{k}\right) \approx \frac{f\left(x_{k}\right)-f\left(x_{k-1}\right)}{x_{k}-x_{k-1}}$来算这个导数。在二维以上，如果最初的Jacobian是$B_{0}$，那么每一轮的jacobian是$\boldsymbol{B}_{k+1}=\boldsymbol{B}_{k}+\left(\left(\boldsymbol{y}_{k}-\boldsymbol{B}_{k} \boldsymbol{s}_{k}\right) \boldsymbol{s}_{k}^{T}\right) /\left(\boldsymbol{s}_{k}^{T} \boldsymbol{s}_{k}\right)$（这个时候叫Broyden方法）。它的conv rate是1.618。
    • 反向内插。实际上，可以认为割线法里面用了一次函数去内插下一轮的解，而我们也可以用更高次的函数。例如，用二次函数，用前三轮去（lagrange法）逼近一个二次函数，找它的零点。它的conv rate是1.839。
    • 一些安全方法。安全是说，其他的一些方法虽然更快，但如果起点离实际值太远，则可能不收敛，而安全的方法不会有这个情况。特别是，如果某一轮使用快速方法的结果跑到bracket外面去了，那么这时应该做一轮安全方法。这些办法的核心基本上就是不断搜索适当的步长和方向，使得这一步迭代真的能提高精度。

第六章

优化问题。我愿意把优化问题想象成在一个（不是无穷大的）势能面上滚动着寻找最低点的小球。这样的问题的基本形式是$\min _{x} f(x)$，并保持 $\vec{g(x)}=\vec{o}$和$\vec{h(x)} \leq\vec{o}$。一般对于优化问题的分类取决于函数f的形式。另外，大多数优化算法（比如说，神经网络）都是主要找局部最小值的，因为全局最小值太难找了。

  • 一些概念。比如，level set和sublevel set。其意义如同地理上的等高线。再比如，coerciveness。这是说一个函数在任何无穷远点都趋于正无穷。

  • 敏感性。例如，一维情况下求解$f^{\prime}\left(x^{\*}\right)=0$，conditioning number可以到$\sqrt{2 \epsilon /\left|f^{\prime \prime}\left(x^{\*}\right)\right|}$。

  • 基本思路：解带约束的优化问题。例如，有一个限制，是$g(x)=0$。用拉格朗日乘子把约束写到待求解的问题中，是$-\nabla f\left(x^{\*}\right)=J_{g}^{T}\left(x^{\*}\right) \lambda$。为了方便，可以弄一个拉格朗日函数$\mathcal{L}: \mathbb{R}^{n+m} \rightarrow \mathbb{R}$，其中$\mathcal{L}(\boldsymbol{x}, \boldsymbol{\lambda})=f(\boldsymbol{x})+\boldsymbol{\lambda}^{T} \boldsymbol{g}(\boldsymbol{x})$，这样只需求此函数的极值。如果又来了一个约束，是$h(x) \leq 0$形式的，那么仍然可用乘子来纳入它，只是此时有额外的限制（这些额外的限制与其他基本的条件一起称作KKT条件）

    $\lambda_{i}^{*} \geq 0, \quad i=m+1, \ldots, m+p$

    $h_{i}\left(x^{*}\right) \lambda_{i}^{*}=0, \quad i=m+1, \ldots, m+p$.

    这里下标的范围表示这些限制只是针对对应约束h的乘子的。

  • 具体算法。——首先是不带约束的优化问题。

    • 一维优化。例如在一个unimodal的一维区域中，使用golden section search(这里和二分法差不多，只是用0.618取代1/2)可以达到线性的收敛速度。但这种办法似乎永远只利用函数的近似一阶导。一种更好的办法是successive parabolic interpolation，二分法的实质是只拟合直线，这一种办法则只拟合一些抛物线。可以想象，它进一步利用到了二阶导。它的收敛是比线性好一些的。

      还有牛顿法。这是利用截断到二阶的泰勒展开，得出迭代$x_{k+1}=x_{k}-f^{\prime}\left(x_{k}\right) / f^{\prime \prime}\left(x_{k}\right)$。它并不在全局收敛，但是一旦收敛就是quadratic的。

      推广到多维，这一类只比较函数值的方法称为direct search method（但多维之下，并不能保证安全性）。对于n维的函数，一般用Nelder-Mead法，这办法是从n+1个随机撒的点开始找起。

    • 另外一类我们熟悉的办法是梯度下降。例如，$\boldsymbol{x}_{k+1}=\boldsymbol{x}_{k}-\alpha_{k} \nabla f\left(\boldsymbol{x}_{k}\right)$。这里有许多的变种，可以养活一批机器学习的专家。其中一个变种是这样的：

      $x_{k+1}=x_{k}-H_{f}^{-1}\left(x_{k}\right) \nabla f\left(x_{k}\right)$

      其中H是Hessian阵。这事实上是牛顿法在高维的推广。好处是，H是对称的，这样算它的逆会省一半资源。坏处是，在远离极值的时候，这H不能确保是正定的。这种牛顿法的复杂度由于求逆的缘故仍然要到每步$\mathcal{O}\left(n^{3}\right)$。还有准牛顿法（从另一角度说，是割线法），$\boldsymbol{x}_{k+1}=\boldsymbol{x}_{k}-\alpha_{k} \boldsymbol{B}_{k}^{-1} \nabla f\left(\boldsymbol{x}_{k}\right)$，其中B是对H的近似，而系数是line search决定的。一种很有效的这类方法是BFGS法。

      考虑到$H_{i+1}^{-1}\left[\nabla f\left(X_{i+1}\right)-\nabla f\left(X_{i}\right)\right] \approx X_{i+1}-X_{i}$， 我们想要近似求解$B_{i+1}\left[\nabla f\left(X_{i+1}\right)-\nabla f\left(X_{i}\right)\right] \approx X_{i+1}-X_{i}$。如果要通过迭代来实现，需要每一步都校正B，即$B_{k+1}=B_{k}+E_{k}$。现在先定义$y_{k}=\nabla f\left(x_{k+1}\right)-\nabla f\left(x_{k}\right), \quad s_{k}=x_{k+1}-x_{k}$，然后可以更方便地推导E。让$rank(E)=2$，即$E_{k}=\alpha u_{k} u_{k}^{T}+\beta v_{k} v_{k}^{T}$。于是得到$\alpha\left(u_{k}^{T} s_{k}\right) u_{k}+\beta\left(v_{k}^{T} s_{k}\right) v_{k}=y_{k}-B_{k} s_{k}$。令$u_{k}=r B_{k} s_{k,} \quad v_{k}=\theta y_{k}$（让两个主要影响方向的因素平衡一下？这也是取E秩2的原因），则解出来所有的东西，在$\mathcal{O}\left(n^{2}\right)$之内（因为只涉及矩阵和向量的乘法）前进一步。

      $B_{k+1}=B_{k}-\frac{B_{k} s_{k} s_{k}^{T} B_{k}}{s_{k}^{T} B_{k} s_{k}}+\frac{y_{k} y_{k}^{T}}{y_{k}^{T} s_{k}}$。

    • 非线性最小二乘。仍然是从残差$r_{i}(x)=y_{i}-f\left(t_{i}, x\right), \quad i=1, \ldots, m$开始——目的是$\phi(x)=\frac{1}{2} r^{T}(x) r(x)$达到最小。不过，这次我们要拟合的是一个非线性函数。高斯-牛顿法就是从最小化$\phi$的角度进行迭代。求两次导，得到

      $\boldsymbol{H}_{\phi}(\boldsymbol{x})=\boldsymbol{J}^{T}(\boldsymbol{x}) \boldsymbol{J}(\boldsymbol{x})+\sum_{i=1}^{m} r_{i}(\boldsymbol{x}) \boldsymbol{H}_{i}(\boldsymbol{x})$

      其中J、H分别是r的jacobian和hessian。迭代时，取

      $\left(J^{T}\left(x_{k}\right) J\left(x_{k}\right)+\sum_{i=1}^{m} r_{i}\left(x_{k}\right) H_{r_{i}}\left(x_{k}\right)\right) s_{k}=-J^{T}\left(x_{k}\right) r\left(x_{k}\right)$

      忽略左边第二项，得到迭代的步长满足$J\left(x_{k}\right) s_{k} \cong-r\left(x_{k}\right)$。于是$x_{k+1}=x_{k}+s_{k}$。这一种方法是否收敛取决于起点。在计算软件中常用的变体是Levenberg-Marquardt方法，核心是$\left(\boldsymbol{J}^{T}\left(\boldsymbol{x}_{k}\right) \boldsymbol{J}\left(\boldsymbol{x}_{k}\right)+\mu_{k} \boldsymbol{I}\right) \boldsymbol{s}_{k}=-\boldsymbol{J}^{T}\left(\boldsymbol{x}_{k}\right) \boldsymbol{r}\left(\boldsymbol{x}_{k}\right)$， 而参数$\mu$的选取则似乎要tricky一些。

第七章

插值问题。插值问题是很有代表性的优化问题，它涉及精确性与普适性的平衡。忽略了前者是欠拟合，而忽略了后者可称为过拟合。同时，它还涉及一些对问题具体情况的认识，因为问题的物理背景常常决定插值函数的性质。

  • 最基本的思路是，给定一系列数据点$\left(t_{i}, y_{i}\right), i=1, \ldots, m$，拿一组基$\phi_i$去插值：

    $f\left(t_{i}\right)=\sum_{j=1}^{n} x_{j} \phi_{j}\left(t_{i}\right)=y_{i}, \quad i=1, \ldots, m$

    在每个数据点处对比，得到一个线性方程组。

    比如，最简单的一种情况是，用单项式作为基底去插值，得到线性方程组的系数矩阵就是范德蒙矩阵。范德蒙矩阵有两个特点，一个是求逆的时间复杂度$\mathcal{O}(n^3)$省不动，一个是它conditioning性质不好（这是因为幂函数，特别是在高次的时候，互相长得太像了）。

  • 还有一种基底，是拉格朗日的“fundamental polynomials” ，其定义为

    $\ell_{j}(t)=\frac{\prod_{k=1, k \neq j}^{n}\left(t-t_{k}\right)}{\prod_{k=1, k \neq j}^{n}\left(t_{j}-t_{k}\right)}, \quad j=1, \ldots, n$

    在这种情况下，$\left(t_{i}, y_{i}\right)$插出来就是$p_{n-1}(t)=y_{1} \ell_{1}(t)+y_{2} \ell_{2}(t)+\cdots+y_{n} \ell_{n}(t)$，而系数矩阵是$I$。

  • 单项式插值难于得到插值函数，而易于计算此函数在特定点的函数值。拉格朗日插值则恰好相反。牛顿插值可以看做这两种办法的一种平衡。它的基底是

    $\pi_{j}(t)=\prod_{k=1}^{j-1}\left(t-t_{k}\right), \quad j=1, \ldots, n$，

    对应的系数矩阵是下三角矩阵，因为对于$i

    $p_{n-1}(t)=x_{1}+\left(t-t_{1}\right)\left(x_{2}+\left(t-t_{2}\right)\left(x_{3}+\left(t-t_{3}\right)\left(\cdots\left(x_{n-1}+x_{n}\left(t-t_{n-1}\right)\right) \cdots\right)\right)\right)$

    这样可以减少做运算的次数。

  • 还有其他一些插值，使用正交的一组基，如勒让德多项式，切比雪夫多项式等等。

  • 插值的收敛性——有些时候，插出来的结果不一定准确，可能与真实的函数相差甚远。不过我们可以判定插值的误差。例如，若f是真实的函数，那么在t点的误差是：

    $f(t)-p_{n-1}(t)=\frac{f^{(n)}(\theta)}{n !}\left(t-t_{1}\right)\left(t-t_{2}\right) \cdots\left(t-t_{n}\right)$

    其中，$\theta$是插值区间里面任意的点。我们要想知道误差，还得知道它的上界。一般来讲，如果n阶导有上界$\left|f^{(n)}(t)\right| \leq M$， 而数据点的最大间隔为$h$，那么$\max _{t \in\left[t_{1}, t_{n}\right]}\left|f(t)-p_{n-1}(t)\right| \leq \frac{M h^{n}}{4 n}$。

    均匀分布的数据点往往在区间两端拟合得不好，而取切比雪夫多项式的零点（在区间两端更密一些）则使得误差更为平均。有时，这会导致定性的区别：Runge函数的多项式拟合，在前者不收敛，在后者收敛。当然数据点的分布常常不是计算家所能决定的。

  • 在有大量数据点的时候，可以做分段插值。最简单的分段插值就是用线段把数据点连起来。

  • 在数据更丰富的时候，可以做埃尔米特插值。这时候你还额外约束函数在每个数据点的导数，使他们连续。埃尔米特立方插值比较常见，这时候你分段用三次多项式去插值——如果分n段，那么有4(n-1)个参数，可以由数据点构造2(n-1)个方程，而导数的条件给出n-2个方程。剩下的n个自由度留着以备问题中其他可能的限制。

第八章

这一章的内容是数值积分与微分。

积分

  • 敏感性上，绝大多数数值积分天生就是well-conditioned的，因为积分本身就是一个带有取平均的过程。
  • 粗略地说，数值积分可以统一地用式子$Q_{n}(f)=\sum_{i=1}^{n} w_{i} f\left(x_{i}\right)$来表示，这叫quadrature rule。用来取函数值的点叫做nodes，而系数则叫做weights。可以用算法来选取这些参数或其中的一部分，使得积分的误差尽可能小。并且，不是为每一个积分选取一组参数，而是首先选定了参数，再用它们去估计所有待计算的积分。例如，给定一组nodes，可以用拉格朗日插值去插某一个函数，再用插出来的结果决定weights：$w_{i}=\int_{a}^{b} \ell_{i}(x) d x, \quad i=1, \ldots, n$。还可以用待定系数法，即考虑给定nodes的情况下认为我们的weights对前n个单项式是准的，以此算出这些weights。例如，$n=3$时，得到的结果是Simpson’s rule: $w_{1}=\frac{b-a}{6}, \quad w_{2}=\frac{2(b-a)}{3}, \quad w_{3}=\frac{b-a}{6}$。这样做会有一个粗糙的上界$\left|I(f)-Q_{n}(f)\right| \leq \frac{1}{4} h^{n+1}\left\|f^{(n)}\right\|_{\infty}$，并有condition number $\sum_{i=1}^{n}\left|w_{i}\right|$。
  • 有些办法来放置这些nodes。例如，Newton-Cotes quadrature。这就是最简单的等分区间。例如，取区间的两端点和中点时，对积分的近似是Simpson’s rule: $S(f)=\frac{b-a}{6}\left(f(a)+4 f\left(\frac{a+b}{2}\right)+f(b)\right)$。当然，就像我们之前在runge函数那里看到的，等分区间并不一定好，反而有可能取点越多误差越大。一种改进是用Chebyshev点，这样的办法叫Clenshaw-Curtis quadrature。之前我们说切比晓夫点是切比晓夫多项式的零点，但这里为方便我们也可以用它的极值点，因为这时的quadrature有progressive的性质，这是说旧nodes全数复用。
  • 更进一步，有一类叫做Gaussian quadrature的，就是说这时候同时把nodes和weights做待定系数，用前2n个单项式来确定。运用这方法时必须注意到确定参数时的区间和实际问题的区间的联系，注意变换。这方法不是progressive的。
  • 于是就有一种尝试想把它变progressive，这就是Kronrod quadrature。对于一个n点的Gaussian rule，可以算出一个2n+1点的Kronrod rule，其中复用了全部的n个nodes。这是现在比较成熟、广泛应用的办法。
  • 另外一个不同的思路是composite quadrature，这是指把区间分割成一些段落，然后用简单的办法去算这些段落上的积分。
  • 有一个参数叫degree of precision，指一个quadrature可以无误差计算其积分的最高阶多项式的阶数（必须对于任意这一阶多项式都可以）。

微分

  • 与积分相反，微分天然是敏感的。
  • 一些尝试：$\frac{f(x)-f(x-h)}{h}$是$\mathcal{O}(h)$准确的，而$\frac{f(x+h)-f(x-h)}{2 h}$则是$\mathcal{O}(h^2)$准确的。但是h本身是一个限制：把h弄太小的话，数据类型的精度上限又会介入误差的计算，造成误差反而变大。有一个办法，叫Richardson外插，可以推断$h=0$时的导数。具体是这样做的：取微分为$F(h)$，取$F(h)=a_{0}+a_{1} h^{p}+\mathcal{O}\left(h^{r}\right)$，并且从问题本身的性质中确定安全的$p, r$。然后，用两种步长$h$和$h/q$来拟合两个线性参数，得到$a_{0}=F(h)+\frac{F(h)-F(h / q)}{q^{-p}-1}+\mathcal{O}\left(h^{r}\right)$就是结果。选取更多种步长，可以提高精度。

第九章

ODE的初值问题。

  • 一些基本概念。我们可以把所有的ODE都先整理成explicit、一阶的方程组。ODE的Explicit指的是把最高阶单独放在等式的一侧，这很容易。我们可以做代换，$u_{1}(t)=y(t), u_{2}(t)=y^{\prime}(t), \ldots, u_{k}(t)=y^{(k-1)}(t)$，这样产生一个一阶的方程组。初值问题（IVP）这时只需要对每一个分量给出一个关于初值的约束。

    论及稳定性，最好的是渐近稳定性，就是说微扰之后的解在远处收敛到正确的解。只说稳定性，则是说微扰之后的解产生的误差不发散。最简单的例子是$y^{\prime}=\lambda y$的三种解：$\lambda$为正则不稳定，为负则渐近稳定，为零则只是稳定。现在把此方程升维，那么可以想见，全部本征值皆负则渐近稳定——其他以此类推。那么就可以推断，判定任何ODE的局部稳定性都可使用Jacobian——考虑一阶近似$z^{\prime}=J_{f}(t, y(t)) z$即可。

  • Euler法。解ODE一般的办法是把区间划分成小段来步进。考虑系统$\boldsymbol{y}^{\prime}=\boldsymbol{f}(t, \boldsymbol{y})$，Euler法就是直接看一阶近似$\boldsymbol{y}_{k+1}=\boldsymbol{y}_{k}+h_{k} \boldsymbol{f}\left(t_{k}, \boldsymbol{y}_{k}\right)$。对于不稳定的解，它的误差是发散的。

  • 在这里所说的误差不包括机器的rounding error，而只是截断到一阶的误差。它又分为全局（$\boldsymbol{e}_{k}=\boldsymbol{y}_{k}-\boldsymbol{y}\left(t_{k}\right)$）和局部（$\ell_{k}=\boldsymbol{y}_{k}-\boldsymbol{u}_{k-1}\left(t_{k}\right)$，$\boldsymbol{u}_{k-1}$指上一个点的函数形式）误差。注意到，对于不稳定的解，全局误差通常大于局部误差之和，因为在开头处弄错函数形式造成的局部误差到了后面会发散掉，称之为propagation。我们只能控制局部误差，但目标却是减小全局的误差。对于局部误差，一种量度是称（在应用于$y^{\prime}=\lambda y$时的）$\ell_{k}=\mathcal{O}\left(h_{k}^{p+1}\right)$的方法具有$p$阶的准确性。例如，Euler法是一阶准确的。其误差之propagation具有一growth factor，即$\boldsymbol{e}_{k+1}=\left(\boldsymbol{I}+h_{k} \boldsymbol{J}_{f}\right) \boldsymbol{e}_{k}+\boldsymbol{\ell}_{k+1}$中右第一项的系数。此系数与$e^{\lambda h}$在$\lambda$的一阶项吻合。

  • Euler法是explicit的，因为其每一步都使用当下确知的信息。另有一类是implicit方法，如backward Euler： $\boldsymbol{y}_{k+1}=\boldsymbol{y}_{k}+h_{k} \boldsymbol{f}\left(t_{k+1}, \boldsymbol{y}_{k+1}\right)$。这导致一个待解的代数方程，可用前章所述诸法来解它。有些implicit方法（不是全部）牺牲了速度，但提高了稳定性，因此人们有时会使用。例如，对于$y^{\prime}=\lambda y$，此法给出的$y_{k}=\left(\frac{1}{1-h \lambda}\right)^{k} y_{0}$即比Euler法给出的$y_{k}=\left({1+h \lambda}\right)^{k} y_{0}$在更大的$\lambda$范围稳定。不过，它的准确性也是一阶的。一阶就说明不太行。

  • 使用梯形法——$\boldsymbol{y}_{k+1}=\boldsymbol{y}_{k}+h_{k}\left(\boldsymbol{f}\left(t_{k}, \boldsymbol{y}_{k}\right)+\boldsymbol{f}\left(t_{k+1}, \boldsymbol{y}_{k+1}\right)\right) / 2$，可知Growth factor为$\frac{1+h \lambda / 2}{1-h \lambda / 2}$，与$e^{\lambda h}$在$\lambda$的二阶项吻合。因此是二阶准确。它的稳定性也不错。

  • 另有一概念为Stiffness。有些物理过程中存在差距较大的多种特征时间，这类过程可能会造成一种情况，即在一个缓变的解周围有许多快速收束的其他暂态解，这个收束的过程较为激烈。这时说方程stiff。例如，$\boldsymbol{J}_{f}$的本征值若差距过大，则可能产生stiffness。Euler法不适合它们，因为受到了步长限制，但backward Euler却很适合。

  • Euler法是Taylor级数法的一种。其外还有Runge-Kutta等别的单步的方法，也有多步的方法。R-K法通过计算小区间中几个点的f来试图模仿出高阶导的效果。其中最简单的是Heun法，即

    $\begin{aligned} \boldsymbol{y}_{k+1} &=\boldsymbol{y}_{k}+\frac{h_{k}}{2}\left(\boldsymbol{k}_{1}+\boldsymbol{k}_{2}\right) \\ \boldsymbol{k}_{1} &=\boldsymbol{f}\left(t_{k}, \boldsymbol{y}_{k}\right) \\ \boldsymbol{k}_{2} &=\boldsymbol{f}\left(t_{k}+h_{k}, \boldsymbol{y}_{k}+h_{k} \boldsymbol{k}_{1}\right) \end{aligned}$

    但较为常用的是模仿到了四阶导的R-K法：

    $\begin{aligned} \boldsymbol{y}_{k+1} &=\boldsymbol{y}_{k}+\frac{h_{k}}{6}\left(\boldsymbol{k}_{1}+2 \boldsymbol{k}_{2}+2 \boldsymbol{k}_{3}+\boldsymbol{k}_{4}\right) \\ \boldsymbol{k}_{1} &=\boldsymbol{f}\left(t_{k}, \boldsymbol{y}_{k}\right) \\ \boldsymbol{k}_{2} &=\boldsymbol{f}\left(t_{k}+h_{k} / 2, \boldsymbol{y}_{k}+\left(h_{k} / 2\right) \boldsymbol{k}_{1}\right) \\ \boldsymbol{k}_{3} &=\boldsymbol{f}\left(t_{k}+h_{k} / 2, \boldsymbol{y}_{k}+\left(h_{k} / 2\right) \boldsymbol{k}_{2}\right) \\ \boldsymbol{k}_{4} &=\boldsymbol{f}\left(t_{k}+h_{k}, \boldsymbol{y}_{k}+h_{k} \boldsymbol{k}_{3}\right) \end{aligned}$

    这个思路（优化的办法）和上一章的quadrature是一致的。它的一个大好处是可以容忍在运行中变化步长，又不损失稳定性准确性。它的坏处是没法对付stiff的方程，以及没法用误差分析来找一个经济的步长。但有些变种，例如implicit R-K可克服前一种困难，embedded R-K可以克服后一种。

  • 还有一些多步法，例如线性的：$\boldsymbol{y}_{k+1}=\sum_{i=1}^{m} \alpha_{i} \boldsymbol{y}_{k+1-i}+h \sum_{i=0}^{m} \beta_{i} \boldsymbol{f}\left(t_{k+1-i}, \boldsymbol{y}_{k+1-i}\right)$。可以看到，这也是implicit的。在启动这种方法的时候需要用explicit的方法先算出来前几个点，称之为Predictor-Corrector配合。多步法可以提升准确性的阶数，可以做stiff ODE，但是难以在运行中变化步长，因为那些希腊字母的系数常常是用步长优化出来的。

第十章

带有边界条件的ODE。

  • 一般形式的一阶BVP表现为$\boldsymbol{y}^{\prime}=\boldsymbol{f}(t, \boldsymbol{y}), \quad a

  • 对于线性的一阶BVP来说，$\boldsymbol{y}^{\prime}=\boldsymbol{A}(t) \boldsymbol{y}$的解可以根据a处的值$\boldsymbol{y}(a)=\boldsymbol{e}_{i}$产生好几个modes。这些modes$\boldsymbol{y}_{i}(t)$组成一个阵$Y(t)$，易知BVP有唯一解等价于$\boldsymbol{Q} \equiv \boldsymbol{B}_{a} \boldsymbol{Y}(a)+\boldsymbol{B}_{b} \boldsymbol{Y}(b)$非奇异。

    此时，可定义$\Phi(t)=Y(t) Q^{-1}$与格林函数

    $\boldsymbol{G}(t, s)=\left\{\begin{aligned} \boldsymbol{\Phi}(t) \boldsymbol{B}_{a} \boldsymbol{\Phi}(a) \boldsymbol{\Phi}^{-1}(s), & a \leq s \leq t \\-\boldsymbol{\Phi}(t) \boldsymbol{B}_{b} \boldsymbol{\Phi}(b) \boldsymbol{\Phi}^{-1}(s), & t

    解出$\boldsymbol{y}(t)=\boldsymbol{\Phi}(t) \boldsymbol{c}+\int_{a}^{b} \boldsymbol{G}(t, s) \boldsymbol{b}(s) d s$。此时，可知BVP的condition number是$\kappa=\max \left\{\|\Phi\|_{\infty},\|G\|_{\infty}\right\}$。

  • 对于BVP来说，不能从一个初值一步一步地进行积分。有一些其他的办法

    • Shooting。固定一端，把它当做IVP来解。导数是猜的。通过这个参数来调节解使得另一端能够接近真实值。这个方法在稳定性上不会有任何的帮助。实际上经常是把区间分成小段分别来做，叫multiple shooting。

    • 有限元。最简单的情况，把区间分成小段，得到近似的导数，例如$u^{\prime \prime}=f\left(t, u, u^{\prime}\right)$在一小段区间上转化为$\frac{y_{i+1}-2 y_{i}+y_{i-1}}{h^{2}}=f\left(t_{i}, y_{i}, \frac{y_{i+1}-y_{i-1}}{2 h}\right)$。如果f线性，就可转化为解一个tridiagonal矩阵的过程。

    • Collocation。就是用$u(t) \approx v(t, \boldsymbol{x})=\sum_{i=1}^{n} x_{i} \phi_{i}(t)$在一定数量的collocation points去拟合。如果这个基底是分段的，例如B-spline，那么这个就化归到我之前了解的有限元。如果是global的，那么就得到一类所谓谱方法。

    • Galerkin法。与其在一定量的点上拟合（使误差为0），不如考虑优化全局的误差。例如，对泊松方程

      $u^{\prime \prime}=f(t), \quad a $u(a)=0, \quad u(b)=0$

      误差为$r(t, \boldsymbol{x})=v^{\prime \prime}(t, \boldsymbol{x})-f(t)=\sum_{i=1}^{n} x_{i} \phi_{i}^{\prime \prime}(t)-f(t)$。

      可以取带权的误差$\int_{a}^{b} r(t, \boldsymbol{x}) w_{i}(t) d t=0, \quad i=1, \ldots, n$，但这样可能会造成需要解一个非对称的矩阵。Galerkin法中取$\int_{a}^{b} r(t, x) \phi_{i}(t) d t=0, \quad i=1, \ldots, n$，即将权重与基底取成一样的。即得到$\int_{a}^{b} v^{\prime \prime}(t, \boldsymbol{x}) \phi_{i}(t) d t=\int_{a}^{b} f(t) \phi_{i}(t) d t, \quad i=1, \ldots, n$ 。把基底$\phi$在区间两端钉到0，分部积分，变成$-\int_{a}^{b} v^{\prime}(t) \phi_{i}^{\prime}(t) d t=\int_{a}^{b} f(t) \phi_{i}(t) d t, \quad i=1, \ldots, n$。这样，就得到一个线性方程组$A x=b$，其中$a_{i j}=-\int_{a}^{b} \phi_{j}^{\prime}(t) \phi_{i}^{\prime}(t) d t, \quad b_{i}=\int_{a}^{b} f(t) \phi_{i}(t) d t$。