《数值分析方法》想来是不少工科研究生的必修课程，最近发现不少工程问题的初始推导都来自数值分析的一些方法，全面梳理其中的算法和结论推导是想做很久的事情之一，虽然很耗时间，但是并不是什么很困难的事情（对比此前各种的Labs），毕竟教会一个engineer如何做researching，比教会一个reseacher如何做好engineering要来得简单（然而很多reseacher乐意去做engineering，但engineer有机会且沉下心做reseach的少之又少）。这本书更是有很多适合工程化的算法，因此本文会基于原教材的编排和定理结论，尤其注重结论的证明和相关算法demo的C++实现，提供：

  • 数值分析的基本思想方法和手段。

  • 遵循数理的结论推导与证明。

  • 开箱即用的各种C++算法。

  • 适量的数值优化/最优化理论拓展；

第一章：数值分析与科学计算引论

误差定性与分析

误差分类

  • 模型误差：建立的数学模型和实际问题之间出现的误差，称为模型误差；只有简化和抽象问题时合理才能得到较好的结果，这种误差不在数值分析的讨论中，同时要衡量这种误差有时候需要根据观测物理量，如温度、电压等，也会引入观测误差；

  • 截断误差/方法误差：在数学模型不能求解精确解时，需要使用数值方法求解它的近似解，其和精确解之间的误差称为截断误差或者方法误差，例如可微的f(x)与泰勒多项式误差表示为：

  • 舍入误差：计算机字长不能完美表示浮点数，因此在数值计算时会产生新的误差，这就是舍入误差，例如带入Π的计算误差：

数值分析主要衡量截断误差和舍入误差；

舍入误差定性

目前对舍入误差的定量估计还没有有效的分析方法，为了确保数值计算正确性通常只进行定性分析；

算法的数值稳定性

计算 并且估计误差。

解：根据分部积分法： 移项得： 故有： 假设使用该递推公式，并且使用： 推到第八项，会发现，而： 可见递推公式是正确的，但是却得到错误的估计结果，这就是初始误差和每步累计的误差造成的错误，其关系是： 故： 由递推可知： 由于阶乘的存在，随着n增大，误差会逐渐放大，因此这种数值计算方法是不稳定的。

我们把计算过程舍入误差不增长的算法称为数值稳定算法，否则为不稳定算法；对于上例，更好的估计方法是使用和的均值估计，然后往前递推，这样误差反而越来越小（往后递推有n的倍增，往前递推则有n分1的递减）。

数值计算中的算法设计艺术

最快速的多项式计算算法——秦九韶算法(Hernor算法)

多项式计算中，计算这样形式的式子：

反复地计算，需要次乘法计算和n次的加法计算，总复杂度为。

容易观察到，低次项的是能够被高次项复用的，这最早于十三世纪被南宋数学家秦九韶提出，在1819年西方提出了相同的Hernor算法(霍纳算法)，表示为：

从该式子，我们能算出每次迭代的系数表示，为：

C++实现

从实现来看，实际上就是复用每次的计算结果，如下： 

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
#include <iostream>
#include <vector>
#include <cmath>

using namespace std;

int Hernor(std::vector<int>& coeff, int value){  //计算coeff[0]x^{n-1} + coeff[1]x^{n-1} + ... + coeff[n-1] 
    int n = coeff.size() - 1;   //最高阶数
    int res = coeff[0];
    for(int i = 1; i < coeff.size(); i++){
        res = res*value + coeff[i];
    }
    cout << "res = " <<  res << endl;
    return res;
}

int main(){

    std::vector<int> vp{1,2,3,4,5,6};
    int value = 5;

    //普通算法
    int normal = 0;
    for(int i=0; i<vp.size(); i++){
        normal += vp.at(i)*pow(value, vp.size()-i-1);
    }
    cout << " normal = " << normal << endl;

    //秦九韶算法
    Hernor(vp, value);

    cout << "done" << endl;
    return 0;
}

拓展

由于秦九韶算法使用了递推的形式，其与乘数求导的链式法则不谋而合，该算法还能进一步推导用于求解多项式在某个点上的一阶导数，表示为： 该式子的推导也比较巧妙，我们知道，任意多项式P(X)总能表示为，Q(X)代表长除法的商，R代表长除法产生的余数。所以使用长除以，商的系数恰好为到，余数即为

令

对上式子两边同时求导得： 即： 这种形式还是符合多项式，一样可以展开为：

值得注意的是，比起函数值系数，这里的末项是而不是，因此c也只递推到,即: 可见我们求出函数值b，通过相同的秦九韶递推，能够得到某一点的一阶导数值，该结论能继续递推，继续求解某点的二阶、任意高阶等导数值。

看理论是很抽象的，以下是应用秦九韶进行函数值、导数值计算的例子：

计算在x = 2函数值和导数值：

解：由递推可知：

故 ；

由：

故

迭代法与开方求值

计算机能够加减乘除四则运算的前提是其能进行二进制的加减和移位，但是如果计算一个数的开方，没有相关的二进制计算方式，因此使用的是迭代法，并且将开方转化成四则运算；

假定a>0, 求解等价于求解方程： 方程求根的相关迭代会在第七章具体介绍，此处构造简单迭代法：令，即：

此处忽略了，这是一个小量的高阶项，则： 故有： 此处，是一个基于初始值的校正量，不断重复这个过程： 且能够证明： 而且收敛速度很快，该方法迭代4-5次就能满足精度。

注：证明过程：

  1. 证明若, 单调且有界，即收敛：

    基本不等式得： 有下界，即, 即单调递减，因此收敛；

  2. 收敛至:

    两侧取极限： 得:

  3. 的收敛速度： 故：

因此这是一种误差平方迭代，收敛速度很快。

数值松弛法

数值松弛法能够通过更少的迭代次数，逼近真值，其公式就是一个加权平均： 细网格细网格粗网格 其中细网格指迭代次数更高、更接近真值的估计值，粗网格则是迭代次数较少的粗值，被称为松弛因子，是一个比细网格更逼近真值的值。

魏晋时期的刘徽在估计圆周率时就使用了类似的思想，估计出，等效于使用半径为10的割圆法到3072边形，而使用松弛法，只需要使用96和192边形，松弛系数为:

大大降低了计算量。

关于松弛因子的确定，可以使用泰勒公式估计，首先给出一个半径为r，内接n边形的面积为：

根据泰勒展开：

使用右式表示，令,有：

故仅考虑二阶因子情况下，能够计算96边形和192边形的松弛因子：

对估计： 得： 这个值是来自二阶的理论计算值，与刘徽根据递增关系得到的比较接近，都能加快割圆法的收敛速度。一个思想类似，但现代的理论是理查森外推法(Richardson Extrapolation)，将在后续的章节介绍。

第二章：插值法

引言

很多实际问题使用y=f(x)表示某种内在规律的数量关系，其中相当一部分是通过实验或者观测得到的，虽然f(x)在区间[a, b]上存在，但可能不连续，或者只给出了函数表，或者是连续，但解析式计算复杂，因此希望使用一个简单函数P(x),使得对成立，这个P(x)就是f(x)对应的插值函数，求解这个函数的过程就是插值法，对应的就是插值节点。

本章仅讨论多项式插值和分段插值，即P(x)为多项式或者分段多项式。

多项式插值

根据原函数的n+1个点，能够建立一个P(x)多项式方程，求解这个多项式方程系数的过程就是多项式插值，y = f(x)取n+1个点：

求解n+1个系数的矩阵表示为：

这样的矩阵被称为范德蒙德(Vandermonde)矩阵,其行列式计算为： 例如矩阵：

范德蒙德行列式计算证明

范德蒙德行列式计算证明也比较简单，转置矩阵的行列式等于原矩阵：

从倒数第二行开始，乘()加到下一行得到：

按第一个数余子式展开得：

继续以此类推，从倒数第二行开始，乘()加到下一行,有：

得证；

定理：满足范德蒙德系数矩阵的插值多项式P(x)是存在且唯一的。

由上述分析可知，选取n+1个点进行多项式插值，系数矩阵是一个范德蒙德矩阵，其行列式的值非零：

即系数矩阵必然满秩，求解、、... 、均存在且唯一，故多项式存在且唯一。

但直接求解n+1个点组成的多项式解是求解插值多项式最繁杂的方式，一般不会使用，后续两节将提出构造多项式更简单的方法；

拉格朗日插值

线性插值

给定区间和端点值和，通过线性拟合函数易得： 点斜式两点式

是通过两个端点的拟合的直线，从两点式可看出，该拟合直线是两个线性函数的加权和，权值系数是两个端点值，两个线性函数分别为： 故： 和称为线性插值基函数，基函数满足：

抛物线插值

对抛物线插值，给定端点、、，满足：

求解抛物线方程宜从基函数角度入手，其三个基函数均需要满足在一个端点值为1，在另外两个端点值为0，故每个基函数都有两个已知零点，故： 故 有：

故第一个基函数为：

同理可知第二、三个基函数, 为：

故已知三个端点值，其抛物线插值的方程表示为：

拉格朗日插值多项式

拉格朗日插值是线性插值、抛物线插值推广到一般情形，给定n+1个端点(、、、)，可以求出其n次的插值多项式，插值多项式满足：

其方法还是n次的插值基函数，对于n次项，这样的基函数有n+1个，均可以表示为：

即对于分子， 用减去不含的所有项，分母是减去不含的所有项。

故：

这就是拉格朗日插值多项式,值得注意的是，通常是n次的，但是也可能小于n次，例如在抛物线插值中，如果三点共线，那么实际上是一条直线，即一次多项式；

插值余项和误差估计

在区间[a,b]使用近似f(x)，则其截断误差为，截断误差也称为插值多项式的余项，有定理：

  • 若在[a,b]上连续，在(a,b)存在，对任何，插值余项满足：

    其中:

插值余项证明

余项证明如下：令K(x)是关于x的待定函数，有：

假设在[a,b]存在固定点x，使得：

注意，上式x、K(x)看成常数。根据假设，在[a,b]n阶连续，n+1阶存在，后两项均为多项式，而且是n次多项式，其n+1阶导数为0；是n+1次多项式，其n+1阶导数就是最高阶数求导数，为(n+1)!;

因此有：

因为可微，存在n+2个零点，由罗尔定理可知，存在n+1个零点，一次类推，存在至少一个零点，有：

故：

得证。

注意，当且仅当f(x)的高阶导数存在时，该插值余项才能被使用，假设找出：

则得到截断误差限为：

拉格朗日插值C++实现

根据拉格朗日多项式：

实现方法比较简单，拟合一个值都需要遍历所有的参考点，其复杂度是： 

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
#include <iostream>
#include <opencv2/core/matx.hpp>
#include <vector>
#include <opencv2/opencv.hpp>

using namespace std;

//输入已知点(2个以上)，需要插值的x坐标,输出插值x坐标对应的逼近值
float lagrange_interpolation(const std::vector<cv::Point2f>& input_points, float target_x){
    float res = 0.0f;
    for(int i=0; i<input_points.size(); i++){
        float li = 1.0f;
        for(int j=0; j<input_points.size(); j++){
            if(i != j){
                li *=  (target_x - input_points[j].x )/ (input_points[i].x - input_points[j].x);
            }
        }
        res += input_points[i].y*li;
    }
    return res;
}

int main(){

    //模拟数据
    std::vector<cv::Point2f> calibPoint{
        {100,100.5},
        {200,298.8},
        {300,399.6},
        {400,402.1},
        {500,596.6},
        {600,698.2},
        {700,799.5},
        {800,900.2},
    };

    //8个点插值至1024个点
    std::vector<cv::Point2f> res_points;

    for(int x=0; x < 1024; x++){
        res_points.push_back({(float)x,lagrange_interpolation(calibPoint, x)});
    }

    //画图示意
    cv::Mat drawPic(1024,1024,CV_8UC3); //底图
    drawPic.setTo(255);

    //结果点
    for(const auto&pos : res_points){
        circle(drawPic, cv::Point(pos.x, pos.y), 2, cv::Scalar(0,0,255), -1);
    }
    
    //参考点
    for(const auto& pos : calibPoint){
        circle(drawPic, cv::Point(pos.x, pos.y), 2, cv::Scalar(255,0,0), -1);
    }

    cv::imshow("drawPic",drawPic);

    cout <<"done" << endl;
    cv::waitKey(0);
    cv::destroyAllWindows();
    return 0;
}

效果： 拉格朗日插值

多点低阶插值余项为0

从上述证明可知：

这意味着：如果要拟合的原本就是一个多项式，而且给定的点数多于该多项式的阶数，其n+1阶导数为0，余项也为0，这也再次印证了拉格朗日插值多项式的唯一性；

勒贝格常数(Lebesgue常数)

勒贝格常数常被用于分析插值算法的数值稳定性，定义为：

对拉格朗日插值，由前述所知有：

勒贝格常数有这样的不等式： 这说明插值函数的误差不仅取决于自身的可逼近性，而且微小的扰动可能被勒贝格常数放大，造成龙格(Runge)现象，对于等距节点，可以计算其勒贝格常数为：

可见是一种指数增长，其基函数在端点附近可能产生剧烈的震荡。

因为插值算法得到的多项式是唯一的，不同插值算法只是基函数不同、计算复杂度不同，使用拉格朗日来定义勒贝格常数只是因为其基函数相对简单，要优化勒贝格常数，需要从插值节点分布的方向进行优化。

切比雪夫节点(Chebyshev Points)

当区间上有一条光滑的曲线，要逼近这条线，很直觉的方法是通过在曲线上均匀地取点，然后构成一个多项式来逼近它，但从勒贝格常数来看，等距取点往往导致区间的两端疯狂摆动，远远偏离原始曲线的龙格现象，甚至随着均匀取点数目的增加，龙格现象会更加剧烈，可见均匀取点是一种错误的直觉；为了驯服这头在区间边缘摆动的野兽，我们只得在区间两端取更多的插值点将其钉住，实现两侧密集、中段稀疏的点分布，这是通往正确的第一缕曙光。

几何运动学中为这种点集的产生提供了思路。在的区间上，一个点绕上半圆周以恒定角速度进行运动，其投影几产生的点集，恰好对应了这种两侧密集、中段稀疏的分布特性，其投影为，这些投影点集，称为切比雪夫（Chebyshev）节点。

第一类切比雪夫多项式

第一类和第二类切比雪夫多项式是切比雪夫多项式族最常用的两种，它们在逼近理论、数值分析和谱方法中起核心作用；

第一类切比雪夫多项式通过三角恒等式定义：

等价的，令，即：

由和差化积公式可知： 即：

而且可知初始条件：

第一类切比雪夫多项式其有相关数学性质：

  • 有界性：，；

  • 奇偶性：n为奇数时为奇函数，为偶数时为偶函数，即

  • 首项系数：：

    首项系数证明：

    令, 有：

    故第一类切比雪夫多项式的最高阶系数：

  • 零点：为切比雪夫-高斯(CG)点，即:

  • 极值点:为切比雪夫-高斯-洛巴托点(CGL)，即： 且

  • 极小极大性质：在所有的首一多项式中在[-1,1]上的最大模最小。

    注：这个性质很有趣，如果在所有多项式中，要找到最高阶系数是1的多项式(首一多项式)，使得其在[-1,1]区间上波动最小，那么最优的选择就是第一类切比雪夫多项式，因为其在[-1,1]区间等幅震荡，幅度就是，其他任意首一多项式震荡幅度不可能小于该多项式，这也说明了切比雪夫插值点插值能最小化误差的本质；

第一类切比雪夫多项式还有一些复杂性质和边值问题、微分方程等研究相关，不在本文的讨论范围内：

  • 正交性：关于权函数正交：

  • 微分方程（施图姆-刘维尔问题）的解：

第二类切比雪夫多项式

第二类切比雪夫多项式的三角恒等式定义为： 等价的：

其递推公式为： 可见第二类切比雪夫多项式和第一类具有相同的递推公式，但初始条件不同：

第二类切比雪夫多项式有以下性质：

  • 首项系数:为，证明略；

  • 零点：, k = 1, 2, ... , n

  • 正交性：关于权函数正交

  • 微分方程：

与第一类切比雪夫多项式存在关系：

第一个差分关系由三角关系容易证明，此处略；

第二项求导证明如下：令: 即： 证毕。

切比雪夫-高斯点(Chebyshev-Gauss Points, CG)和切比雪夫-高斯-洛巴托点(Chebyshev-Gauss-Lobatto Poitns, CGL)

如前文所述，两类切比雪夫点都来自三角函数，其中段稀疏、端点密集的插值点特性使得勒贝格常数更优，对CG点： 即的基数倍，所以其不含-1和1两个端点，取的值都是内部点，广用于插值和积分等不重视边值条件的情况；

而CGL点： 其同时包含了-1和1两个端点，比较适用于谱方法、边值问题；

未完待续