数值分析方法：C++的工程实现
《数值分析方法》想来是不少工科研究生的必修课程，最近发现不少工程问题的初始推导都来自数值分析的一些方法，全面梳理其中的算法和结论推导是想做很久的事情之一，虽然很耗时间，但是并不是什么很困难的事情（对比此前各种的Labs），毕竟教会一个engineer如何做researching，比教会一个reseacher如何做好engineering要来得简单（然而很多reseacher乐意去做engineering，但engineer有机会且沉下心做reseach的少之又少）。这本书更是有很多适合工程化的算法，因此本文会基于原教材的编排和定理结论，尤其注重结论的证明和相关算法demo的C++实现，提供：
数值分析的基本思想方法和手段。
遵循数理的结论推导与证明。
开箱即用的各种C++算法。
适量的数值优化/最优化理论拓展；
第一章：数值分析与科学计算引论
误差定性与分析
误差分类
模型误差：建立的数学模型和实际问题之间出现的误差，称为模型误差；只有简化和抽象问题时合理才能得到较好的结果，这种误差不在数值分析的讨论中，同时要衡量这种误差有时候需要根据观测物理量，如温度、电压等，也会引入观测误差；
截断误差/方法误差：在数学模型不能求解精确解时，需要使用数值方法求解它的近似解，其和精确解之间的误差称为截断误差或者方法误差，例如可微的f(x)与泰勒多项式误差表示为：
舍入误差：计算机字长不能完美表示浮点数，因此在数值计算时会产生新的误差，这就是舍入误差，例如带入Π的计算误差：
数值分析主要衡量截断误差和舍入误差；
舍入误差定性
目前对舍入误差的定量估计还没有有效的分析方法，为了确保数值计算正确性通常只进行定性分析；
算法的数值稳定性
计算
解：根据分部积分法：
我们把计算过程舍入误差不增长的算法称为数值稳定算法，否则为不稳定算法；对于上例，更好的估计方法是使用
数值计算中的算法设计艺术
最快速的多项式计算算法——秦九韶算法(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
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代表长除法产生的余数。所以使用
令
对上式子两边同时求导得：
值得注意的是，比起函数值系数，这里的末项是
看理论是很抽象的，以下是应用秦九韶进行函数值、导数值计算的例子：
计算
解：由递推可知：
故
由：
故
迭代法与开方求值
计算机能够加减乘除四则运算的前提是其能进行二进制的加减和移位，但是如果计算一个数的开方，没有相关的二进制计算方式，因此使用的是迭代法，并且将开方转化成四则运算；
假定a>0, 求解
此处忽略了
注：证明过程：
证明若
, 单调且有界，即 收敛： 基本不等式得：
有下界，即 , 即 单调递减，因此 收敛； 收敛至 : 两侧取极限： 得: 的收敛速度： 故：
因此这是一种误差平方迭代，收敛速度很快。
数值松弛法
数值松弛法能够通过更少的迭代次数，逼近真值
魏晋时期的刘徽在估计圆周率时就使用了类似的思想，估计出
大大降低了计算量。
关于松弛因子的确定，可以使用泰勒公式估计，首先给出一个半径为r，内接n边形的面积为：
根据泰勒展开：
使用右式表示
故仅考虑二阶因子情况下，能够计算96边形和192边形的松弛因子：
第二章：插值法
引言
很多实际问题使用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+1个，均可以表示为：
即对于分子， 用
故：
这就是拉格朗日插值多项式,值得注意的是，
插值余项和误差估计
在区间[a,b]使用
若
在[a,b]上连续， 在(a,b)存在，对任何 ，插值余项满足： 其中:
插值余项证明
余项证明如下：令K(x)是关于x的待定函数，有：
假设在[a,b]存在固定点x，使得：
注意，上式x、K(x)看成常数。根据假设，
因此有：
因为
故：
注意，当且仅当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
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常数)
勒贝格常数常被用于分析插值算法的数值稳定性，定义为：
对拉格朗日插值，由前述所知有：
勒贝格常数有这样的不等式：
可见是一种指数增长，其基函数在端点附近可能产生剧烈的震荡。
因为插值算法得到的多项式是唯一的，不同插值算法只是基函数不同、计算复杂度不同，使用拉格朗日来定义勒贝格常数只是因为其基函数相对简单，要优化勒贝格常数，需要从插值节点分布的方向进行优化。
切比雪夫节点(Chebyshev Points)
当区间上有一条光滑的曲线，要逼近这条线，很直觉的方法是通过在曲线上均匀地取点，然后构成一个多项式来逼近它，但从勒贝格常数来看，等距取点往往导致区间的两端疯狂摆动，远远偏离原始曲线的龙格现象，甚至随着均匀取点数目的增加，龙格现象会更加剧烈，可见均匀取点是一种错误的直觉；为了驯服这头在区间边缘摆动的野兽，我们只得在区间两端取更多的插值点将其钉住，实现两侧密集、中段稀疏的点分布，这是通往正确的第一缕曙光。
几何运动学中为这种点集的产生提供了思路。在
第一类切比雪夫多项式
第一类和第二类切比雪夫多项式是切比雪夫多项式族最常用的两种，它们在逼近理论、数值分析和谱方法中起核心作用；
第一类切比雪夫多项式通过三角恒等式定义：
等价的，令
由和差化积公式可知：
而且可知初始条件：
第一类切比雪夫多项式其有相关数学性质：
有界性：
； 奇偶性：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两个端点，比较适用于谱方法、边值问题；
C++实现：等距与切比雪夫插值的拉格朗日插值对比
1 |
|
输出： 1
2fixed lebesgue = 10772.4
chebyshev lebesgue = 2.90082
上述代码没有虽然使用标准的CG点的标准形式，实际上是等效的，也可以进一步对比CG、CGL点的插值效果，但对于该函数插值而言，效果是相近的，可经过改代码进行对比：
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
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
using namespace std;
double test_func(double x){
return 1.0 / (1 + 25*x*x);
}
double lagrange_interpolation(const std::vector<cv::Point2d>& input_points, double target_x){
double res = 0.0;
for(int i=0; i<input_points.size(); i++){
double li = 1.0;
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;
}
std::vector<double> getChebyshevsPoint(int num, double low, double high, int mode){
std::vector<double> res;
double theta = 0.0;
switch (mode) {
case 0: { // 一般切比雪夫点（Chebyshev points of the first kind）
for(int i = 0; i < num; i++){
theta = (2*i+1)*CV_PI/(2*num);
double xi = std::cos(theta);
double mapped_x = 0.5*((high-low)*xi + (low + high));
res.push_back(mapped_x);
}
break;
}
case 1: { // Chebyshev-Gauss points (CG)
for(int i = 0; i < num; i++){
theta = (2*i+1)*CV_PI/(2*num);
double xi = std::cos(theta);
double mapped_x = 0.5*((high-low)*xi + (low + high));
res.push_back(mapped_x);
}
break;
}
case 2: { // Chebyshev-Gauss-Lobatto points (CGL)
for(int i = 0; i <= num; i++){
theta = i*CV_PI/num;
double xi = std::cos(theta);
double mapped_x = 0.5*((high-low)*xi + (low + high));
res.push_back(mapped_x);
}
break;
}
}
return res;
}
int main(){
const int ratio = 600;
//200个参考点，表征原函数趋势
std::vector<cv::Point2d> vp_orig;
for(int i=0; i<=200; i++){
double xi = -1 + i/100.0;
vp_orig.push_back({xi, test_func(xi)});
}
//获取切比雪夫插值点
std::vector<double> vec_chebyshev_x = getChebyshevsPoint(20, -1, 1,0);
std::vector<cv::Point2d> vec_chebyshev_point;
for(const auto& pos : vec_chebyshev_x){
vec_chebyshev_point.push_back({pos, test_func(pos)});
}
//获取CG插值点
std::vector<double> vec_chebyshev_cg = getChebyshevsPoint(20, -1, 1,1);
std::vector<cv::Point2d> vec_chebyshev_point_cg;
for(const auto& pos : vec_chebyshev_cg){
vec_chebyshev_point_cg.push_back({pos, test_func(pos)});
}
//获取CGL插值点
std::vector<double> vec_chebyshev_cgl = getChebyshevsPoint(20, -1, 1,2);
std::vector<cv::Point2d> vec_chebyshev_point_cgl;
for(const auto& pos : vec_chebyshev_cgl){
vec_chebyshev_point_cgl.push_back({pos, test_func(pos)});
}
//切比雪夫点插值效果
std::vector<cv::Point2d> vp_chebyshev;
for(int i=0; i<=200; i++){
double xi = -1 + i/100.0;
vp_chebyshev.push_back({xi, lagrange_interpolation(vec_chebyshev_point, xi)});
}
//CG插值效果
std::vector<cv::Point2d> vp_chebyshev_cg;
for(int i=0; i<=200; i++){
double xi = -1 + i/100.0;
vp_chebyshev_cg.push_back({xi, lagrange_interpolation(vec_chebyshev_point_cg, xi)});
}
//CGL点插值效果
std::vector<cv::Point2d> vp_chebyshev_cgl;
for(int i=0; i<=200; i++){
double xi = -1 + i/100.0;
vp_chebyshev_cgl.push_back({xi, lagrange_interpolation(vec_chebyshev_point_cgl, xi)});
}
//画图
cv::Mat drawPic(512,1024,CV_8UC3); //底图
drawPic.setTo(255);
//绿色线--表征原函数趋势
for(const auto&pos : vp_orig){
circle(drawPic, cv::Point(ratio*pos.x, ratio*pos.y), 2, cv::Scalar(0,255,0), -1);
}
//红色线--切比雪夫点的拉格朗日插值
for(const auto&pos : vp_chebyshev){
circle(drawPic, cv::Point(ratio*pos.x, ratio*pos.y), 2, cv::Scalar(0,0,255), -1);
}
//蓝色线--CG点的拉格朗日插值
for(const auto&pos : vp_chebyshev_cg){
circle(drawPic, cv::Point(ratio*pos.x, ratio*pos.y), 2, cv::Scalar(255,0,0), -1);
}
//黑色线--CGL点的拉格朗日插值
for(const auto&pos : vp_chebyshev_cgl){
circle(drawPic, cv::Point(ratio*pos.x, ratio*pos.y), 2, cv::Scalar(0,0,0), -1);
}
// 添加标签
cv::putText(drawPic, "Green: Original Function", cv::Point(20, 30),
cv::FONT_HERSHEY_SIMPLEX, 0.7, cv::Scalar(0,255,0), 2);
cv::putText(drawPic, "Red: Normal Chebyshev Lagrange", cv::Point(20, 60),
cv::FONT_HERSHEY_SIMPLEX, 0.7, cv::Scalar(0,0,255), 2);
cv::putText(drawPic, "Blue: CG Chebyshev Lagrange", cv::Point(20, 90),
cv::FONT_HERSHEY_SIMPLEX, 0.7, cv::Scalar(255,0,0), 2);
cv::putText(drawPic, "Black: CGL Chebyshev Lagrange", cv::Point(20, 120),
cv::FONT_HERSHEY_SIMPLEX, 0.7, cv::Scalar(0,0,0), 2);
cv::imshow("drawPic", drawPic);
/**************************************************************勒贝格常数估算**************************************************************/
//计算基函数lk在target_x处的值, k = 0,1...,n
auto calculate_lk = [](const std::vector<cv::Point2d>& input_points, int k, double target_x){
double lk = 1.0;
int n = input_points.size()-1; //插值的阶数，三点二阶、四点三阶......
for(int j=0; j<=n; j++){
if(j != k){
lk *= (target_x - input_points[j].x) / (input_points[k].x - input_points[j].x);
}
}
return lk;
};
//估计勒贝格常数，按0.01的精度取最大值作勒贝格常数值
auto getLebesgue = [calculate_lk](const std::vector<cv::Point2d>& input_points, double low, double high){
double sum, maxsum;
maxsum = sum = 0.0;
int n = input_points.size()-1; //插值的阶数，三点二阶、四点三阶......
for(double x=low; x<=high; x += 0.01){
sum = 0.0;
for(int k=0; k<=n; k++){
sum += std::abs(calculate_lk(input_points, k, x));
}
maxsum = sum > maxsum ? sum : maxsum;
}
return maxsum;
};
//切比雪夫点勒贝格常数：
auto chebyshev_leb = getLebesgue(vec_chebyshev_point, -1, 1);
cout << "normal chebyshev lebesgue = " << chebyshev_leb << endl;
//cg点勒贝格常数：
auto chebyshev_leb_cg = getLebesgue(vec_chebyshev_point_cg, -1, 1);
cout << "cg chebyshev lebesgue = " << chebyshev_leb_cg << endl;
//cgl点勒贝格常数：
auto chebyshev_leb_cgl = getLebesgue(vec_chebyshev_point_cgl, -1, 1);
cout << "cgl chebyshev lebesgue = " << chebyshev_leb_cgl << endl;
cout <<"done" << endl;
cv::waitKey(0);
cv::destroyAllWindows();
return 0;
}
均差与牛顿插值多项式
从上节可知，拉格朗日插值是根据基函数推导的，推导方式简单，所以也用其基函数定义了勒贝格常数。勒贝格常数与插值算法实现无关，可以根据点分布进行优化，切比雪夫分布点是优化插值龙格现象的有效方法。但从其实现方法也可知，拉格朗日的计算复杂度比较高，计算一个插值点的插值函数值就需要遍历所有的点，牛顿多项式旨在设计一种可逐次生成插值多项式的方法，以下。
线性插值的点斜式可以表示为：
其可以看成
进一步推广到抛物线插值，有:
故：
均差(差商)定义
我们称f(x)关于点
如上式，令
结合k阶差商，插值多项式就能够表示为：
均差的相关重要性质如下：
- k阶均差可以表示为函数值
, ,..., 的线性组合，即：
关于线性组合这个关系可由数学归纳法证明，此处我们反向结合拉格朗日表达式可以粗略说明这个等式的推导（注意，非证明，因为我们不经证明承认了牛顿插值多项式），根据拉格朗日插值多项式表示：
由于我们在上一节已经证明，不同方法得到的插值多项式是唯一的，因此 ，因此其系数也是对应的，故对 系数是相同的,有： 表述为：k阶均差(差商)等于通过k+1个节点的插值多项式的最高次项系数。 - k阶均差可以表示为函数值
- 由性质1可知均差的计算具有对称性，k阶差商的表示可以推广为：
即分母减去什么，只要和分子的均差项对应，结果是相等的。
- 由性质1可知均差的计算具有对称性，k阶差商的表示可以推广为：
- 若
f(x)在[a,b]上存在n阶导数，且节点，则n阶均差和导数满足关系： 此结论可有罗尔定理证明；更简单的记忆是对应到泰勒公式（实际上泰勒公式是基于牛顿多项式推导，但我们更熟悉泰勒公式），泰勒公式展开的系数是使用导数关系描述的，而牛顿多项式的系数是使用均差关系描述的，根据插值多项式的唯一性，二者自然有上述关系。
- 若
综上所述，可以通过函数值的线性组合、低阶递推、导数三种方法计算均差；
牛顿插值余项
在拉格朗日插值余项中，我们的表示是：
其中
根据上述差商和导数关系，自然得到牛顿插值余项可以表示为：
牛顿插值余项证明
除了对应到拉格朗日，我们还可以正向证明一下，使用到性质1和性质3的定义式：
根据线性组合定义式子，我们可以得到关于f(x)的均差表示式：
故：
为了得到x看成区间[a,b]一点，构造有：
故有：
得证，以此类推，我们会有：
故余项：
得证,该余项等价于拉格朗日定义的余项，而且更具有一般性，当f(x)是离散给出或者导数不存在时，拉格朗日余项失效，但该牛顿插值余项仍然适用。
C++实现：牛顿插值
牛顿插值通过递推得到均差表，插值算法只使用均差表每列的首行，理论上存储优化可到n是参考插值点，m是需要插值的点数，对应的拉格朗日插值的复杂度为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
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
using namespace std;
//输入：插值点
//输出：均差表
std::vector<std::vector<double>> caldivided_difference(const std::vector<cv::Point2d>& points)
{
int n = points.size();
std::vector<std::vector<double>> res(n, std::vector<double>(n, 0.0));
for(int j=0; j<n; j++){
res[j][0] = points.at(j).y;
}
for(int i=1; i<n; i++){
for(int j=i; j<n; j++){
res[j][i] = (res[j][i-1] - res[j-1][i-1])/(points.at(j).x - points.at(j-i).x);
}
}
return res;
}
//输入：插值点、均差表和目标插值点，输出插值函数值
double newton_interpolation(const std::vector<cv::Point2d>& points, const std::vector<std::vector<double>> ÷d_diff, double target_x){
double res = divided_diff[0][0];
double term = 1.0;
for(int i=1; i<divided_diff.size(); i++) //牛顿插值只用到均差表每列的首行
{
term *= (target_x - points.at(i-1).x);
res += divided_diff[i][i]*term;
}
return res;
}
int main(){
//模拟数据
std::vector<cv::Point2d> 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::Point2d> res_points;
for(int x=0; x < 1024; x++){
res_points.push_back({(double)x,newton_interpolation(calibPoint,caldivided_difference(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("newton_interpolation",drawPic);
cout <<"done" << endl;
cv::waitKey(0);
cv::destroyAllWindows();
return 0;
}
差分形式牛顿插值
前述牛顿插值多项式适用于任意节点分布，但通常为了方便，实际应用经常使用等距节点（虽然可能导致数值不稳定），其牛顿公式还能进一步优化。
先引入一些概念：
令等距节点h为步长,因此有：
对应的二阶差分和n阶差分：
函数值与差分
为了描述方便，引入两个算子，有： I称为不变算子，E称为位移算子，由此有：
反过来，f(x)的函数值可以使用各阶差分表示：
算子关系：
导数和前向差分转换
由导数定义就可知：
均差和前向差分转换
设等距节点步长为h，有：
故：
牛顿前插公式
从牛顿多项式看，为：
设取了n个等距节点，步长为h，自变量步数为t,令
对应的其余项形式变为：
埃尔米特（Hermite）插值
插值算法均要求插值函数在插值点的值等于原函数值，有的实际问题不仅要求函数值相等，还要求其导数值、高阶导数值相等，满足这种要求的多项式就称为埃尔米特插值多项式。
重节点均差与泰勒插值
给出一个关于均差的结论（原书也没有证明）：
- 设
为[a,b]上相异节点，则均差 是其变量的连续函数，即：
由此结论，定义重节点均差，即重复节点的均差为：
故：
故得到泰勒插值多项式：
泰勒多项式的导数满足： m+1个插值条件，可造出次数不超过m次的埃尔米特插值多项式。
两个典型埃尔米特插值函数
三次埃尔米特插值
三次埃尔米特插值给定条件为：
给定条件数是
由
可计算得待定系数为：
此即为三次埃尔米特插值多项式。
其余项为：
该结论能从罗尔定理开始证明，但是我们在证明拉格朗日和均差时都用过了，直接推导也可以，如果是
此处因为一阶导数条件施加在
两点三次埃尔米特插值
两点三次埃尔米特插值是很常用的埃尔米特插值公式，其约束条件只有两个点，常用于某段区间插值，满足：
采用拉格朗日相同的基函数推导，两点三次埃尔米特插值的待定函数为：
基函数分别满足：
观察到
因为满足
故：
故有：
同理可计算其它三个基函数待定系数，此处略，两点三次埃尔米特插值最后解析式记为：
其余项形式和上述三次埃尔米特插值还是一样的，事实上n代表给出条件数，重根项取决于条件数中函数值、导数、高阶导数施加的条件多少，例如如果某个点函数值、导数、二阶导数都是0，那么这个项就是一个三阶项。
两点三次埃米尔特插值的余项表示为：
三次样条插值
分段高次插值会造成严重的龙格现象，低次插值的光滑程度不满足应用关系，三次的埃尔米特插值虽然比低次插值光滑程度有所改善，但需要提高导数信息，结果也仅仅是一阶连续的光滑程度，因此应用场景局限，此几节不再讨论。
基于这些缺点，三次样条插值被提出，由于该方法太著名，在很早的插值算法文章就提过并且给出相关实现，见数值分析：三次样条插值算法(Cubic Spline Interpolation)，本节略。
第三章：函数逼近与快速傅里叶变换
第二章插值法是函数逼近的一种，在本章，我们给出的函数逼近定义是：对函数类A中给定的函数
魏尔斯特拉斯(Weierstrass)定理
- 设
, 对任意 , 总是存在一个代数多项式p(x),使得： 在[a,b]区间上一致成立。
魏尔斯特拉斯定理主要告诉我们，任何连续函数总能找到一个代数多项式逼近，且在某个闭区间中误差小于任意给定的正数。
很多定理反过来证明了这点，如泰勒展开的余项总是一个与高阶导数相关的无限小，这里要介绍的是最经典的证明方法————伯恩斯坦多项式证明，它遵循我们很熟悉的分布，即二项分布。
伯恩斯坦多项式证明
伯恩斯坦多项式证明用到二次项分布的三个引理：
对引理1，是典型二项分布关系式，其概率和就是1。
由引理1可以证明引理2： 实际上引理2和二项分布的期望有关，即：
同理，使用引理1和2硬算来证明，有： 引理三显然和二项分布的方差计算相关，利用引理三计算二项分布的方差：
扯远了，总而言之，我们引入并且证明了上述三个引理，构造伯恩斯坦多项式为：
因为f(x)是连续的，严谨来说，任意连续区间映射到[0,1]区间，f(x)一致连续，即对
估计误差:
利用引理3对第二项误差进行控制，设
第二步的放大我们引入了引理三部分，此处必然有
故最终有：
证毕。故显然的，伯恩斯坦多项式能够逼近闭区间任意连续的f(x),这也是魏尔斯特拉斯逼近定理的经典证明，但二次项多项式收敛速度太慢，工程上一般不使用伯恩斯坦多项式作为函数逼近多项式。
函数逼近与函数空间数学理论
在后续的章节涉及逼近理论的一些术语，这里先按原书概念和额外拓展记录一些定义。
范数与赋范线性空间
为了定义线性空间中向量的长度，需要引入范数的概念，因为是描述向量长度的，它应该满足非负、线性延申、两点直线最短等直观的几何性质，从代数语言的角度，范数的标准定义是：
- 若S为线性空间，
，若存在唯一实数 （这个符号只是代表占位符， 既可以是一个变量或者函数）满足条件：
正定性(非负的、不能无中生有，有又变无)：
，当且仅当 时， ； 齐次性(延申是线性的)：
， 三角不等式(直线距离最短)：
其中
简单的，对n维的(欧几里得)向量空间
1-范数：
，即绝对值和； 2-范数：
，即平方和开根号； 范数（最大范数）: ,即绝对值最大者；
对连续函数空间的C[a,b],
-** 2-范数**：
范数（最大范数）: ,即绝对值最大者；
内积与内积空间
范数与赋范线性空间是用于定义向量长度的，本身没有方向信息，因此此处还要引入内积和内积空间来定义角度信息。
过去在线性代数中，我们定义向量空间
推广到线性空间
, ；当 为实数域 对应的线性空间时，有 ； , ；拓展一下，课本没有提到， ,这个在下面的柯西施瓦茨不等式复数域证明就用到了。 , ； , 当且仅当u=0，(u,u) = 0.
柯西-施瓦茨(Cauchy-Schwarz)不等式
要证明也很简单，但是书本的证明不那么严谨，以下。
当v=0时，等号成立； 当
令
即：
问题出在了展开(u,v)不等于(v,u)，因此直接写成2
故：
在复数域中此处的竖线代表模长，即令
柯西施瓦茨不等式在任意数域的线性空间成立，证明也比较抽象，如果放在二、三维空间很容易理解，例如在
