数值计算day6-曲线拟合与插值

上节课主要介绍了特征值与特征向量的概念，低阶矩阵的特征值可以通过列出特征方程求解，高阶矩阵则可以通过幂法与反幂法迭代求解出最大特征值与最小特征值（模），要求出矩阵的全部特征值则需要借助矩阵的 QR分解来将矩阵相似化为一个上三角矩阵，相似化过程不改变矩阵的特征值，因此转化后的上三角矩阵的对角线元素即为原矩阵的特征值。本节课主要介绍曲线拟合（机器学习中的线性回归算法）与插值算法。

1. 拟合与插值（Fitting, Interpolation）

• 曲线拟合是给定一个数据集，构建一个数学方程能够最好的匹配这些数据点
• 插值则是构建一个多项式方程来完美地契合给定的数据集，当数据点比较少的时候，一个简单的多项式方程就可以实现插值，而数据点比较多的时候，更多地采用的是样条插值，在各个区间构建不同的多项式方程

2. 线性曲线拟合

线性曲线拟合是使用如下线性方程（一元一阶多项式，对应机器学习线性回归算法中只有一个特征的情况）来拟合给定的数据集：$$f(x)=a_1x+a_0$$ 若数据集中只有两个点，则该线性方程可以准确求解出来，通过这两个点；当数据集超过两个点时，该线性方程无法通过所有的数据点，此时需要找出拟合数据点最好的那个方程。

如何衡量拟合的准确度？对数据点\((x_i,y_i)\)，定义\(y_i\)和\(f(x_i)\)的差值为残差\(r_i\)： $$r_i=y_i-f(x_i)$$
一个可选的衡量测度是计算残差和：$$E=\sum^{n_{i=1}r_i=\sum}n_{i=1}[y_i-(a_1x_i+a_0)]$$ 然而这种测度有时候会不好用，因为当有些点残差为很大的正数，而有些点残差为很小的负数时，将会给出一个接近于\(0\)的残差和，

另外一个可选的测度是残差的绝对值和：$$E=\sum^{n_{i=1}|r_i|=\sum}n_{i=1}|y_i-(a_1x_i+a_0)|$$ 该测度总是正的，可以用来衡量拟合的精度，但无法用来确定哪个是最好的拟合曲线，因为不同的拟合曲线有可能给出相同的残差绝对值和。

因此，一个误差函数应该既能够衡量拟合曲线的好坏，又能够用来确定出唯一地拟合曲线，目前广泛采用的是残差的平方和：$$E=\sum^n_{i=1}r2_i=\sum^{n_{i=1}[y_i-(a_1x_i+a_0)]}2$$ 采用线性最小二乘回归算法就可以找出使得上述误差最小的那条拟合曲线。给定数据集，\(E\)是关于\(a_0,a_1\)的非线性方程，结合微积分相关概念，可知\(E\)取最小值的点\((a_1,a_0)\)满足：$$\frac{\partial E}{\partial a_0}=-2\sum^n_{i=1}(y_i-a_1x_i-a_0) = 0$$ $$\frac{\partial E}{\partial a_1}=-2\sum^n_{i=1}(y_i-a_1x_i-a_0) x_i= 0$$ 化简有：$$na_0+(\sum^n_{i=1}x_i) a_1 = \sum^n_{i=1}y_i$$ $$(\sum^{n_{i=1}x_i)a_0+(\sum}n_{i=1}x^2_i) a_1 = \sum^n_{i=1}x_iy_i$$ 其解为： $$a_1=\frac{n\sum^{n_{i=1}x_iy_i-(\sum}n_{i=1}x_i)(\sum^{n_{i=1}y_i)}{n\sum}n_{i=1}x^2_i -(\sum^n_{i=1}x_i)2}$$ $$a_0=\frac{(\sum^n_{i=1}x_i)2(\sum^{n_{i=1}y_i)-(\sum}n_{i=1}x_iy_i)(\sum^{n_{i=1}x_i)}{n\sum}n_{i=1}x^2_i -(\sum^n_{i=1}x_i)2}$$ 如果定义 $$S_x=\sum^{n_{i=1}x_i,S_y=\sum}n_{i=1}y_i,S_{xx}=\sum^n_{i=1}x2_i,S_{xy}=\sum^n_{i=1}x_iy_i$$ 则有 $$a_1=\frac{nS_{xy}-S_xS_y}{nS_{xx}-(S_x)^{2},a_0=\frac{S_{xx}S_{y}-S_{xy}S_x}{nS_{xx}-(S_x)}2}$$
MATLAB实现：

function [a1,a0] = LinearRegression(x, y)
% LinearRegration calculates the coefficients a1 and a0 of the linear
% equation y = a1*x + a0 that best fit n data points.
% Input variables:
% x    A vector with the coordinates x of the data points.
% y    A vector with the coordinates y of the data points.
% Output variable:
% a1   The coefficient a1.
% a0   The coefficient a0.

nx = length(x);
ny = length(y);
if nx ~= ny
    disp('ERROR: The number of elements in x must be the same as in y.')
    a1 = 'Error';
    a0 = 'Error';
else
Sx = sum(x);
Sy = sum(y);
Sxy = sum(x.*y);
Sxx = sum(x.^2);
a1 = (nx*Sxy - Sx*Sy)/(nx*Sxx - Sx^2);
a0 = (Sxx*Sy - Sxy*Sx)/(nx*Sxx - Sx^2);
end

3. 非线性曲线拟合

很多实际应用中，使用线性方程的拟合效果不够好，采用非线性拟合则能够达到比较好的效果。

3.1 常见的几种非线性方程

（这里同线性拟合类似，只考虑单变量的情况，对应线性回归算法中的一个特征）：$$y=bx^m$$ $$y=be^{mx}$$ $$y=\frac{1}{mx+b}$$ 对非线性拟合，很容易将其改写为线性形式：$$ln(y)=mln(x)+ln(b)$$ $$ln(y) = mx+ln(b)$$ $$\frac{1}{y}=mx+b$$ 将数据集进行相应的特征变换后，采用线性最小二乘法求解相应的参数。然而，对于给定的数据集，如何判断使用哪一种非线性方程进行拟合呢？一个可行的方案是首先通过绘制数据集中的点（当特征很多时，很难可视化），确定是选择线性还是非线性进行拟合，如果选择非线性拟合，具体选择哪一种非线性方程则要取决于其背后的物理现象和数学原理。同样的，可以通过以特定的方式绘制数据点并检查这些点是否符合直线来预测所提出的非线性函数是否具有很好的拟合能力。
MATLAB实现：

texp = 2:2:30;
vexp = [9.7 8.1 6.6 5.1 4.4 3.7 2.8 2.4 2.0 1.6 1.4 1.1 0.85 0.69 0.6];
vexpLOG = log(vexp);
R = 5E6;
[a1,a0] = LinearRegression(texp, vexpLOG)
b = exp(a0)
C = -1/(R*a1)
t = 0:0.5:30;
v = b*exp(a1*t);
plot(t,v,texp,vexp,'or')

3.2 二次及高阶多项式拟合

多项式是指具有如下形式的函数（单变量）：$$f(x)=a_nx^n+a_{n-1}x+...+a_1x+a_0$$ \(n\)称作多项式的阶数。如果给定数据集中有\(n\)个点，则可以确定\(n-1\)个从1阶到\(n-1\)阶的拟合多项式，其中\(n-1\)阶多项式可以实现通过每一个数据点，但并不是多项式的阶数越高越好，这取决于实际问题，当数据集点本身就具有噪声污染时，提高拟合的阶数并没有意义。（同机器学习中的过拟合，overfitting类似）

对多项式拟合，采用多项式回归算法来确定各个参数。给定\(n\)个数据点，想要使用\(m\)阶多项式对其进行拟合：$$f(x)=a_mx^m+a_{m-1}x+...+a_1x+a_0$$ 以\(m=2\)为例，首先定义总误差：$$E=\sum^{n_{i=1}[y_i-(a_2x}2_i+a_1x_i+a_0)]^2$$ 要让总误差最小，取对应的导数为\(0\)即可： $$\frac{\partial E}{\partial a_0}=-2\sum^{n_{i=1}(y_i-a_2x}2_i-a_1x_i-a_0)=0$$ $$\frac{\partial E}{\partial a_1}=-2\sum^{n_{i=1}(y_i-a_2x}2_i-a_1x_i-a_0)x_i=0$$ $$\frac{\partial E}{\partial a_2}=-2\sum^{n_{i=1}(y_i-a_2x}2_i-a_1x_i-a_0)x^2_i=0$$ 将其更改为线性系统形式：

\[na_0+(\sum^n_{i=1}x_i)a_1+(\sum^n_{i=1}x^2_i)a_2=\sum^n_{i=1}y_i$$ $$(\sum^n_{i=1}x_i)a_0+(\sum^n_{i=1}x^2_i)a_1+(\sum^n_{i=1}x^3_i)a_2=\sum^n_{i=1}x_iy_i$$ $$(\sum^n_{i=1}x^2_i)a_0+(\sum^n_{i=1}x^3_i)a_1+(\sum^n_{i=1}x^4_i)a_2=\sum^n_{i=1}x^2_iy_i$$ 该线性系统的解即为多项式的各个参数。 MATLAB实现： ``` clear, clc x = 0:0.4:6; %x=1:4 y = [0 3 4.5 5.8 5.9 5.8 6.2 7.4 9.6 15.6 20.7 26.7 31.1 35.6 39.3 41.5]; n = length(x); m = 4; for i = 1:2*m xsum(i) = sum(x.^(i)); end % Beginning of Step 3 a(1,1) = n; b(1,1)= sum(y); for j= 2:m+1 a(1,j)=xsum(j-1); end for i = 2:m+1 for j=1:m+1 a(i,j)= xsum(j+i-2); end b(i,1) = sum(x.^(i-1).*y); end % Step 4 p=(a\b)' for i = 1:m+1 Pcoef(i) = p(m+2-i); end epsilon=0:0.1:6; stressfit=polyval(Pcoef,epsilon); plot(x,y,'ro',epsilon,stressfit,'k','linewidth',2) xlabel('Strain','fontsize',20) ylabel('Stress (MPa)','fontsize',20) %%%%%%%%%%%%%%%%%%%% >> help polyval polyval - 多项式计算 此 MATLAB 函数 返回在 x 处计算的 n 次多项式的值。输入参数 p 是长为 n+1 的向量， 其元素是按要计算的多项式降幂排序的系数。 y = polyval(p,x) [y,delta] = polyval(p,x,S) y = polyval(p,x,[],mu) [y,delta] = polyval(p,x,S,mu) See also polyder, polyfit, polyint, polyvalm Reference page for polyval Other functions named polyval ```  上述多项式拟合可以扩展到更一般的情况：$$f(x)=c_1f_1(x)+c_2f_2(x)+...+c_mf_m(x)$$ 定义总误差为 $$E=\sum^n_{i=1} [y_i-\sum^m_{j=1}c_jf_j(x_i)]^2$$ 令偏导为$0$ $$\frac{\partial E}{\partial c_k} = 2\sum^n_{i=1}(y_i-\sum^m_{j=1}C_jf_j(x_i))(-f_k(x_i))=0,k=1,...,m$$ 改为线性系统形式，有$$\sum^n_{i=1}\sum^m_{j=1}c_{j}f_j(x_i)f_k(x_i)=\sum^n_{i=1}y_if_k(x_i),k=1,...,m$$ 矩阵形式为 $$\begin{bmatrix}\sum^n_{i=1}f^2_1(x_i) & \sum^n_{i=1}f_1(x_i)f_2(x_i) & \cdots & \sum^n_{i=1}f_m(x_i)f_1(x_i) \\ \sum^n_{i=1}f_1(x_i)f_2(x_i) & \sum^n_{i=1}f^2_2(x_i) & \cdots & \sum^n_{i=1}f_m(x_i)f_2(x_i)\\ \vdots & \vdots & \ddots & \vdots\\ \sum^n_{i=1}f_1(x_i)f_m(x_i) & \sum^n_{i=1}f_2(x_i)f_m(x_i) & \cdots & \sum^n_{i=1}f^2_m(x_i) \end{bmatrix}\begin{bmatrix}c_1\\c_2\\ \vdots \\ c_m \end{bmatrix}=\begin{bmatrix}\sum^n_{i=1}y_if_1(x_i)\\\sum^n_{i=1}y_if_2(x_i) \\ \vdots \\ \sum^n_{i=1}y_if_m(x_i) \end{bmatrix}\]

MATLAB实现：
\(f(x)=\frac{A}{x}+\frac{Be^{-2x^2}}{x}\)

clear all
x = [0.6 0.8 0.85 0.95 1.0 1.1 1.2 1.3 1.45 1.6 1.8];
y = [0.08 0.06 0.07 0.07 0.07 0.06 0.06 0.06 0.05 0.05 0.04];
a(1,1) = sum(1./x.^2);
a(1,2) = sum(exp(-2*x.^2)./x.^2);
a(2,1) = a(1,2);
a(2,2) = sum(exp(-4*x.^2)./x.^2);
b(1,1) = sum(y./x);
b(2,1) = sum((y.*exp(-2*x.^2))./x);
AB = a\b
xfit = 0.6:0.02:1.8;
yfit = AB(1)./xfit + AB(2)*exp(-2*xfit.^2)./xfit;
plot(x,y,'o',xfit,yfit)

4. 多项式插值

如前所述，插值是让拟合的曲线通过所有给定的数据点（总误差为0，机器学习中的\(E_{in}\)为\(0\)），对\(n\)个点，总是能找到一个\(n-1\)阶的多项式函数通过所有的点：$$P(x) = a_0+a_1x+...+a_mx^m=y$$ 矩阵形式 $$\begin{bmatrix}x_1^m & x_1^{m-1} & \cdots & x_1 & 1\ x_2^m & x_2^{m-1} & \cdots & x_2 & 1\ \vdots & \vdots & \ddots & \vdots & \vdots\ x_n^m & x_n^{m-1} & \cdots & x_n & 1 \end{bmatrix}\begin{bmatrix}a_m \ a_{m-1}\ \vdots \ a_0\end{bmatrix}=\begin{bmatrix}y_1 \ y_2\ \vdots \ y_n\end{bmatrix}$$

4.1 多项式插值的拉格朗日形式：

考虑两个点\((x_1,y_1)\), \((x_2,y_2)\)，其一阶拉格朗日多项式具有如下形式：$$y=a_1(x-x_2)+a_2(x-x_1)$$ 令该函数通过这两个点，得 $$a_1=\frac{y_1}{x_1-x_2}, a_2=\frac{y_2}{x_2-x_1}$$ 因此多项式为：$$y=\frac{y_1(x-x_2)}{x_1-x_2}+\frac{y_2(x-x_1)}{x_2-x_1}$$ 类似地，对三个点\((x_1,y_1),(x_2,y_2),(x_3,y_3)\)，可得其二阶拉格朗日具有形式$$y=\frac{(x-x_2)(x-x_3)}{(x_2-x_1)(x_3-x_1)}y_1+\frac{(x-x_1)(x-x_3)}{(x_1-x_2)(x_3-x_2)}y_2+\frac{(x-x_1)(x-x_2)}{(x_1-x_3)(x_1-x_3)}y_3$$ 因此，拉格朗日多项式的通解形式为：$$f(x)=\sum^n_{i=1}y_i\prod_{j\neq i}\frac{(x-x_j)}{(x_i-x_j)}$$ MATLAB实现：

function Yint = LagrangeINT(x,y,Xint)
% LagrangeINT fits a Lagrange polynomial to a set of given points and
% uses the polynomial to determines the interpolated value of a point.
% Input variables:
% x  A vector with the x coordinates of the given points.
% y  A vector with the y coordinates of the given points.
% Xint  The x coordinate of the point to be interpolated.
% Output variable:
% Yint  The interpolated value of Xint.

n = length(x);
for i = 1:n
    L(i) = 1;
    for j = 1:n
        if j ~= i
            L(i)= L(i)*(Xint-x(j))/(x(i)-x(j));
        end
    end
end
y.*L
Yint = sum(y.*L);
%%%%%%%%
>> x = [1 2 4 5 7];
>> y = [52 5 -5 -40 10];
>> Yinterpolated = LagrangeINT(x,y,3)

ans =

   -5.7778    2.6667   -4.4444   13.3333    0.2222

Yinterpolated =

    6.0000

>> i=1;
>> for x1=0:0.1:10
y1(i)=LagrangeINT(x,y,x1);
i = i+1;
end
>> x1 = 0:0.1:10;
>> plot(x,y,'*',x1,y1)

4.2 牛顿多项式插值

上述拉格朗日多项式插值有一个缺陷，当数据集增加一个插值点时，每个系数都需要重新计算，而牛顿多项式插值则不需要重新计算：$$f(x)=a_1+a_2(x-x_1)+a_3(x-x_1)(x-x_2)+...+a_n(x-x_1)(x-x_2)...(x-x_{n-1})$$
以两个数据点为例 \((x_1,y_1)\) \((x_2,y_2)\), $$f(x)=a_1+a_2(x-x_1)$$ 容易计算出$$a_1=y_1,a_2=\frac{y_2-y_1}{x_2-x_1}$$ 当增加一个数据点\((x_3,y_3)\)时，$$f(x)=a_1+a_2(x-x_1)+a_3(x-x_1)(x-x_2)$$ 可以看到 \(f(x_1)=a_1=y_1\), \(f(x_2)=a_1+a_2(x_2-x_1)=y_2 \rightarrow a_2= \frac{y_2-y_1}{x_2-x_1}\)，这两个参数都没有发生变化，因此只需计算\(a_3\)，通过\(f(x_3)=y_3\)得，$$ a_1+a_2(x_3-x_1)+a_3(x_3-x_1)(x_3-x_2)=y_3 $$ $$a_3=\frac{y_3-y_1-\frac{y_2-y_1}{x_2-x_1}(x_3-x_1)}{(x_3-x_1)(x_3-x_2)}=\frac{\frac{y_3-y_1}{x_3-x_2}-\frac{y_2-y_1}{x_2-x_1}\frac{x_3-x_1}{x_3-x_2}}{x_3-x_1}$$ $$=\frac{\frac{y_3-y_2}{x_3-x_2}+\frac{y_2-y_1}{x_3-x_2}-\frac{y_2-y_1}{x_3-x_2}\frac{x_3-x_1}{x_2-x_1}}{x_3-x_1}=\frac{\frac{y_3-y_2}{x_3-x_2}-\frac{y_2-y_1}{x_2-x_1}}{x_3-x_1}$$ 一直进行下去，就可以得到牛顿多项式插值的通解形式：
第一层：\(y_1,y_2,y_3,...\) 其中\(a_1=y_1\)
第二层：\(f[x_2,x_1],f[x_3,x_2],f[x_4,x_3],...\) 其中\(f[x_2,x_1]=\frac{y_2-y_1}{x_2-x_1}=a_2,...\)
第三层：\(f[x_3,x_2,x_1],f[x_4,x_3,x_2],f[x_5,x_4,x_3],...\) 其中\(f[x_3,x_2,x_1]=\frac{f[x_3,x_2]-f[x_2,x_1]}{x_3-x_1}=a_3,...\)

MATLAB实现：

function Yint = NewtonsINT(x,y,Xint)
% NewtonsINT fits a Newtons polynomial to a set of given points and
% uses the polynomial to determines the interpolated value of a point.
% Input variables:
% x  A vector with the x coordinates of the given points.
% y  A vector with the y coordinates of the given points.
% Xint  The x coordinate of the point to be interpolated.
% Output variable:
% Yint  The interpolated value of Xint.

n = length(x);
a(1) = y(1);
for i = 1:n-1
    divDIF(i,1) = (y(i+1)-y(i))/(x(i+1)-x(i));
end
for j = 2:n-1
    for i = 1:n-j
        divDIF(i,j) = (divDIF(i+1,j-1) - divDIF(i,j-1))/(x(j+i) - x(i));
    end
end
for j=2:n
    a(j)=divDIF(1,j-1);
end
Yint=a(1);
xn=1;
for k=2:n
    xn=xn*(Xint-x(k-1));
    Yint=Yint+a(k)*xn;
end

5. 样条插值

5.1 线性样条插值

当数据点比较少时，使用低阶多项式就可以实现插值，然而当数据点比较多时，插值的多项式需要很高的阶数，这会造成插值点不准确（同回归中预测不准确）。因此，这个时候可以采用区间插值的方法，即每两个或三个相邻点之间，构建具有不同参数的插值多项式（阶数一致）。

线性插值就是每两个相邻的点之间构建各自的线性函数进行插值。此时，第一个点\((x_1,y_1)\)与第二个点\((x_2,y_2)\)之间的插值函数(拉格朗日插值公式)为：$$f_1(x)=\frac{x-x_2}{x_1-x_2}y_1+\frac{x-x_1}{x_2-x_1}y_2$$，同样地，很容易得到第\(i\)个点\((x_i,y_i)\)与第\(i+1\)个点之间的插值函数为：$$f_i(x)=\frac{x-x_{i+1}}{x_i-x_{i+1}}y_i+\frac{x-x_i}{x_{i+1}-x_i}y_{i+1},i=1,...,n-1$$ 可以看到，线性插值是一种连续插值，因为相邻的两个多项式在同一节点的值相同，但在该节点处的斜率（导数）是不连续的。
MATLAB实现：

function Yint = LinearSpline(x, y, Xint)
% LinearSpline calculates interpolation using linear splines.
% Input variables:
% x    A vector with the coordinates x of the data points.
% y    A vector with the coordinates y of the data points.
% Xint The x coordinate of the interpolated point. 
% Output variable:
% Yint The y value of the interpolated point.

n = length(x);
for i = 2:n
    if Xint < x(i)
        break
    end
end
Yint = (Xint - x(i))*y(i-1)/(x(i-1)-x(i)) + (Xint - x(i-1))*y(i)/(x(i)-x(i-1));

5.2 二次样条插值

给定\(n\)个数据点，二次样条插值是在\(n-1\)个区间内各自构建一个二次曲线插值数据：$$f_i(x)=a_ix^2+b_ix+c_i,i=1,2,...,n-1$$ 可以看到每个多项式有\(3\)个待定系数，共有\(3n-3\)个待定系数，因此求解规则包括：

• 每个二次曲线在区间两端点\(x_i\)与\(x_{i+1}\)的值为\(y_i,y_{i+1}\)：$$f_i(x_i)=a_ix^2_i+b_ix_i+c_i=y_i,i=1,2,...,n-1$$ $$f_{i}(x_{i+1})=a_{i}x^2_{i+1}+b_ix_{i+1}+c_{i}=y_{i+1},i=1,2,...,n-1$$ 共构建出\(2n-2\)个方程
• 每两个相邻区间的节点处，对应多项式的斜率相同（保证导数连续）：$$f'i(x)=2a_ix_i+b_i=f'{i+1}(x)=2a_{i+1}x_{i+1}+b_{i+1},i=1,2,...,n-2$$ 共构建\(n-2\)个方程
• 上述条件共构建了\(3n-4\)个方程，想要求解\(3n-3\)个参数，还少一个条件，因此，通常假设第一个点或者最后一个点的二阶导数为零：$$f''_1(x_1)=2a_1=0$$ 即第一个或最后一个多项式为线性多项式
  

5.3 三次样条插值

对\(n\)个数据点，三次样条插值是在\(n-1\)个区间内各自构建一个三次多项式进行插值。三次样条插值包含标准形式与拉格朗日形式。考虑如下标准形式：$$f_i(x)=a_ix^3+b_ix2+c_ix+d_i,i=1,2,...,n-1$$ 共有\(4n-4\)个待定参数。求解规则：

• 每个三次曲线在区间两端点\(x_i\)与\(x_{i+1}\)的值为\(y_i,y_{i+1}\)：$$f_i(x_i)=a_ix^3_i+b_ix2_i+c_ix_i+d_i=y_i,i=1,2,...,n-1$$ $$f_{i}(x_{i+1})=a_{i}x^3_{i+1}+b_ix2_{i+1}+c_{i}x_{i+1}+d_i=y_{i+1},i=1,2,...,n-1$$ 共构建出\(2n-2\)个方程
• 每两个相邻区间的节点处，对应多项式的斜率相同（保证一阶导数连续）：$$f'i(x)=3a_ix^{2_i+2b_ix_i+c_i=f'_{i+1}(x_{i+1})=3a_{i+1}x}2_{i+1}+2b_{i+1}x_{i+1}+c_{i+1},i=1,2,...,n-2$$ 共构建\(n-2\)个方程
• 每两个相邻区间的节点处，对应多项式的二阶导数相同（当从一条曲线转换为下一条曲线时，斜率的变化是连续的）：$$f''i(x)=6a_ix_i+2b_i=f''{i+1}(x)=6a_{i+1}x_{i+1}+2b_{i+1},i=1,2,...,n-2$$ 共构建\(n-2\)个方程
• 上述条件共构建了\(4n-6\)个方程，还差两个条件，通常要求第一个点和最后一个点的二阶导数为零(自然三次样条插值，natural cubic splines)：$$6a_1x_1+2b_1=0$$ $$6a_{n-1}x_n+2b_{n-1}=0$$

三次样条插值的标准形式易于理解，使用简单，但需要求解\(4n-4\)个线性方程，计算复杂。考虑三次样条的拉格朗日形式：

• 从三次多项式的二阶导数出发，三次多项式的二阶导数\(f''_i(x)\)是一个线性函数，采用拉格朗日多项式插值形式（点\((x_i,f''_i(x_i)),(x_{i+1},f''_i(x_{i+1}))\)）：$$f''i(x)=\frac{x-x{i+1}}{x_i-x_{i+1}}f''i(x_i)+\frac{x-x{i}}{x_{i+1}-x_i}f''i(x)$$ 积分可得：$$f'i(x)=\frac{(x-x)^{2}{2(x_i-x_{i+1})}f''_i(x_i)+\frac{(x-x_i)}2}{2(x_{i+1}-x_{i})}f''i(x)+C_1$$ 再次积分，有 $$f_i(x)=\frac{(x-x_{i+1})^{3}{6(x_i-x_{i+1})}f''_i(x_i)+\frac{(x-x_{i})}3}{6(x_{i+1}-x_{i})}f''i(x)+C_1x+C_2$$ 令\(f_i(x_i)=y_i,f_i(x_{i+1})=y_{i+1}\) 解得

\[C_1=\frac{y_{i+1}-y_i}{x_{i+1}-x_i}+\frac{x_{i+1}-x_{i}}{6}f''_i(x_i)-\frac{x_{i+1}-x_{i}}{6}f''_i(x_{i+1}) \]

\[C_2=\frac{y_ix_{i+1}-y_{i+1}x_i}{x_{i+1}-x_i}-\frac{(x_{i+1}-x_i)}{6}f''_i(x_i)x_{i+1}+\frac{(x_{i+1}-x_i)}{6}f''_i(x_{i+1})x_i$$ 因此 $$\begin{aligned}f_i(x)&=\frac{(x_{i+1}-x)^3}{6(x_{i+1}-x_i)}f''_i(x_i)+\frac{(x-x_{i})^3}{6(x_{i+1}-x_{i})}f''_i(x_{i+1})+C_1x+C_2\\ &=\frac{(x_{i+1}-x)^3}{6(x_{i+1}-x_i)}f''_i(x_i)+\frac{(x-x_{i})^3}{6(x_{i+1}-x_{i})}f''_i(x_{i+1})\\ &+[\frac{y_i}{x_{i+1}-x_i}-\frac{(x_{i+1}-x_i)f''_i(x_i)}{6}](x_{i+1}-x)\\ & +[\frac{y_{i+1}}{x_{i+1}-x_i}-\frac{(x_{i+1}-x_i)f''_i(x_{i+1})}{6}](x-x_{i})\end{aligned},i=1,2,...,n-1$$ 定义 $h_i=x_{i+1}-x_i,a_i=f''(x_i)$，则$$f_i(x)=\frac{a_i}{6h_i}(x_{i+1}-x)^3+\frac{a_{i+1}}{6h_i}(x-x_i)^3+[\frac{y_i}{h_i}-\frac{h_ia_i}{6}](x_{i+1}-x)+[\frac{y_{i+1}}{h_i}-\frac{h_ia_{i+1}}{6}](x-x_{i})$$ 进一步利用 $f'_i(x_{i+1})=f'_{i+1}(x_{i+1}),i=1,...,n-2$，可得：$$6[\frac{y_{i+2}-y_{i+1}}{h_{i+1}}-\frac{y_{i+1}-y_i}{h_i}]=h_ia_i+2(h_i+h_{i+1})a_{i+1}+h_{i+1}a_{i+2}$$ 注意$a_1,a_2,...,a_{n-1},a_n$是待定数，然而这里只有$n-2$个方程，还需要两个条件，利用标准形式中的最后一个条件，即第一个点和最后一个点的二阶导数为零，可得$a_1=0,a_n=0$。 MATLAB实现： ``` function Yint = CubicSplines(x,y,xint) n = length(x); a1=0; an=0; A = zeros(n-2,n-2); b = zeros(n-2,1); for i=1:n-2 h(i) = x(i+1)-x(i); h(i+1) = x(i+2)-x(i+1); if i == 1 A(i,i) = 2*(h(i)+h(i+1)); A(i,i+1) = h(i+1); b(i,1) = 6*((y(i+2)-y(i+1))/h(i+1)-(y(i+1)-y(i))/h(i)); elseif i == n-2 A(i,i-1) = h(i); A(i,i) = 2*(h(i)+h(i+1)); b(i,1) = 6*((y(i+2)-y(i+1))/h(i+1)-(y(i+1)-y(i))/h(i)); else A(i,i-1) = h(i); A(i,i) = 2*(h(i)+h(i+1)); A(i,i+1) = h(i+1); b(i,1) = 6*((y(i+2)-y(i+1))/h(i+1)-(y(i+1)-y(i))/h(i)); end end p = A\b; a(1) = 0; a(n) = 0; for i=1:n-2 a(1+i) = p(i); end for i = 2:n if xint <= x(i) break end end h(i) = x(i)-x(i-1); Yint = a(i-1)/(6*h(i))*(x(i)-xint)^3+a(i)/(6*h(i))*(xint-x(i-1))^3+(y(i-1)/h(i)-a(i-1)*h(i)/6)*(x(i)-xint)+(y(i)/h(i)-a(i)*h(i)/6)*(xint-x(i-1)); end %%%%%%%%%%%%%%%%%%%%%%%%%%%% x = [5 7.5 10 12.5 15 17.5 20 22.5 25 27.5 30]; h = [3.3 7.5 41.8 51.8 61 101.1 132.9 145.5 171.4 225.8 260.9]; xint = 5:0.5:30; l = length(xint); yint = zeros(l,1); for i = 1:l yint(i) = CubicSplines(x,h,xint(i)); end figure subplot(1,2,1) plot(x,h,'k*',xint,yint,'r-') title('cubic') p = interp1(x,h,xint,'spline'); subplot(1,2,2) plot(x,h,'k*',xint,p,'r-') title('spline') ```  #### 6. 总结 本节课主要介绍了曲线的拟合与插值，与机器学习中的回归算法相似。给定一组数据集，拟合是找到最好的曲线来匹配数据集中的点（线性回归，非线性回归），插值则是找到通过所有数据点的曲线。线性回归是最简单的曲线拟合算法，非线性回归可以通过特征变换转化为线性回归，而其中的高阶多项式拟合还可以进一步扩展到更一般的情形，求解方法均是基于最优点处的导数为零得到。对于插值，任意$n$个点都可以找到一个$n-1$阶的多项式函数通过所有数据点，求解方法有拉格朗日法和牛顿法， 拉格朗日法计算简单，但增加数据点时，所有参数均需重新计算，而牛顿法则不需要重新计算，但计算较为复杂。另一方面，当数据点比较少时，多项式插值较为简单，但对于数据点较多的情况，多项式的阶数需要很高，计算复杂，并且插值时并不可靠（过拟合）。样条插值是在每个区间内构建各自的多项式插值函数，一般有线性样条插值、二次样条插值以及三次样条插值。其中为求解简便，三次样条插值除了标准形式外，还有一个便于求解的拉格朗日形式。\]