函数逼近

基本概念

使用较简单的函数表示已知的复杂函数或未知函数的过程，称为函数逼近

函数逼近要求近似函数在整体上近似于所需函数，因此需要利用整体误差进行判断

函数逼近的目标函数 $f(x)$ ：有表达式的复杂函数，或者无表达式的表格函数（即已知部分数据点的值）

函数逼近的近似函数 $p(x)$ ：属于某种简单的函数类 $\Phi$

逼近的判据：某种度量意义下误差函数 $p(x)-f(x)$ 最小

即求解以下问题：

\min_{p(x)\in\Phi}||p(x)-f(x)||

为此，我们需要引入函数空间的概念，事实上，某定义域上的所有函数构成线性空间，而某个区间上的一定阶连续可微函数 $C^k[a,b]$ 也构成线性空间

可以使用函数范数来衡量函数的大小，函数范数是对向量范数的连续延拓

下面是 $C[a,b]$ 上的几种常用范数：

∞-范数： $||f(x)||_\infty=\max_{x\in[a,b]}|f(x)|$
1-范数： $||f(x)||_1=\int_a^b|f(x)|\mathrm dx$
2-范数： $||f(x)||_2=\sqrt{\int_a^b(f(x))^2\mathrm dx}$

同样地，我们定义 $C[a,b]$ 上的内积：

\left<u,v\right>=\int_a^bu(x)v(x)\mathrm dx

定理（函数内积的Cauchy-Schwarz不等式）

|\left<u,v\right>|^2\leq\left<u,u\right>\cdot\left<v,v\right>|

定理设 $S$ 为实内积空间， $u_1,u_2,\cdots,u_n\in S$ ，则Gram矩阵

\mathbf G=\begin{bmatrix}\left<u_1,u_1\right>&\cdots&\left<u_1,u_n\right>\\ \vdots&\ddots&\vdots\\ \left<u_n,u_1\right>&\cdots&\left<u_n,u_n\right> \end{bmatrix}

非奇异的充要条件是 $u_1,\cdots,u_n$ 线性无关

定义（权函数） 针对实函数空间 $C[a,b]$ ，如果满足

非负： $\rho(x)\geq 0,\ \forall x\in[a,b]$
多项式可积： $\int_a^b x^k\rho(x)\mathrm dx$ 对任意 $k=0,1,\cdots$ 存在
无局部恒为零：若 $\int_a^b g(x)\rho(x)\mathrm dx=0$ 则 $g(x)\equiv 0$

则此函数可以作为一个权函数

定义（加权内积） 对满足上述定义的权函数 $\rho(x)$ ，定义以下加权内积：

\left<u,v\right>=\int_a^b\rho(x)u(x)v(x)\mathrm dx

加权内积可以诱导出内积范数：

||u||=\sqrt{\left<u,u\right>}

用不同的误差度量方式，可以对函数逼近产生不同的意义：

使用∞-范数度量误差，则为最佳一致逼近，此误差具有最高的全局一致性
使用1-范数和2-范数则类似于一种平均误差，其中2-范数的情况称为最佳平方（最小二乘）逼近

事实上，求解最佳一致逼近很复杂，仅考虑最佳平方逼近

连续函数的最佳平方逼近

现对 $f(t)\in C[a,b]$ 进行函数逼近，选取近似空间为线性空间 $\Phi={\rm span}\{\phi_1(t),\cdots,\phi_n(t)\}$

则，问题描述为，求解 $S(t)=\sum\limits_{j=1}^nx_j\phi_j(t)\in\Phi$ ，使 $||S(t)-f(t)||_2$ 最小

法方程法

令

F\equiv ||S(t)-f(t)||_2^2=\left<\sum\limits_{j=1}^nx_j\phi_j-f,\sum\limits_{j=1}^nx_j\phi_j-f\right>

=\sum\limits_{j=1}^n\sum\limits_{l=1}^nx_jx_l\left<\phi_j,\phi_l\right>-2\sum\limits_{j=1}^nx_j\left<f,\phi_j\right>+\left<f,f\right>

这是一个含有自变量 $x_1,\cdots,x_n$ 的 $n$ 元二次函数（在函数空间的基给定的情况下，上式中所有内积均为常数）

我们要寻求其全局最小值，只需要令梯度为零

\frac{\partial F}{\partial x_k}=2\sum\limits_{j=1}^nx_j\left<\phi_j,\phi_k\right>-2\left<f,\phi_k\right>

上面的式子对 $k=1,2,\cdots,n$ 分别取 $0$ 构成了一个线性方程组，其等价于

\mathbf {Gx}=\mathbf b

其中 $\mathbf G$ 为之前提到的Gram矩阵：

\mathbf G=\begin{bmatrix}\left<\phi_1,\phi_1\right>&\cdots&\left<\phi_1,\phi_n\right>\\ \vdots&\ddots&\vdots\\ \left<\phi_n,\phi_1\right>&\cdots&\left<\phi_n,\phi_n\right> \end{bmatrix}

同时

\mathbf b=\begin{bmatrix}\left<f,\phi_1\right>\\ \vdots\\ \left<f,\phi_n\right> \end{bmatrix}

此方法被称为法方程法，法方程存在唯一正确解的前提是Gram矩阵对称正定

其对称性是显然的，正定性需要选取合适的基函数

法方程法的讨论

假设法方程的解为 $\mathbf x^*$ ，得到的逼近函数为 $S^*$

根据法方程

\sum\limits_{j=1}^n x_j^*\left<\phi_j,\phi_k\right>=\left<f,\phi_k\right>

根据内积的双线性性，这等价于

\left<S^*,\phi_k\right>=\left<f,\phi_k\right>

即

\left<S^*-f,\phi_k\right>=0

这说明 $S^*$ 正交于任何一个基向量，进而说明

S^*-f\in\Phi

即，误差函数与求解空间正交

法方程法的逼近误差

||\delta||_2^2=||S^*-f||_2^2

=\left<S^*-f,S^*-f\right>=\left<S^*-f,S^*\right>-\left<S^*-f,f\right>

由于 $S^*\in\Phi$ ， $S^*-f\perp S^*$ ，于是

||\delta||_2^2=\left<f-S^*,f\right>=||f||_2^2-\sum\limits_{j=1}^mx_j^*\left<\phi_j,f\right>

最佳平方逼近多项式

如果我们在 $C[0,1]$ 上求解，选取 $\Phi=\mathbb P_{n-1}$ ，即次数不超过 $n-1$ 的所有多项式函数的集合，则基函数集为

\{1,t,\cdots,t^{n-1}\}

于是

\left<\phi_j,\phi_k\right>=\int_0^1t^{j+k-2}\mathrm dt=\frac 1{j+k-1}

Gram矩阵为

\mathbf G_n=\begin{bmatrix} 1&\frac12&\cdots&\frac1n\\ \frac12&\frac13&\cdots&\frac1{n+1}\\ \vdots&\vdots&\ddots&\vdots\\ \frac1n&\frac1{n+1}&\cdots&\frac1{2n-1} \end{bmatrix}

这是一个高度病态的Hilbert矩阵，因此采用单项式作为基函数的逼近方法是不太好的

正交函数族

一组两两正交的函数组成一个正交函数族

函数的Fourier级数实际上就是在正交函数族上的逼近

然而，在这里我们不希望使用三角形状的正交函数族，因为它们在数值计算上并不好用，于是我们寻求一组正交多项式

一种简单的构造正交多项式的方式为：由基 $\{1,t,t^2,\cdots,t^{n-1}\}$ 出发，使用Gram-Schmidt正交化

\phi_1(t)=1

\phi_k(t)=t^{k-1}-\sum\limits_{j=1}^{k-1}\frac{\left<t^{k-1},\phi_j\right>}{\left<\phi_j,\phi_j\right>}\phi_j(t)

此时，由于选取的基正交，Gram矩阵为对角矩阵，于是法方程法的结果就非常简单

S^*(t)=\sum\limits_{k=1}^n\frac{\left<f,\phi_k\right>}{\left<\phi_k,\phi_k\right>}\phi_k(t)

上面这种操作得到的正交多项式族由以下几个初始变量决定：

定义域 $[a,b]$ ，即进行内积时的积分限
权函数 $\rho(t)$ ，可以选择加权内积作为Gram-Schmidt正交化时使用的内积

当上述初始变量取值

$[a,b]=[-1,1]$
$\rho(t)=1$

时，得到的正交多项式族为勒让德多项式，其递推公式如下

P_0(t)=1

P_1(t)=t

(k+1)P_{k+1}(t)=(2k+1)tP_k(t)-kP_{k-1}(t)

对于一般的定义域 $[a,b]$ ，只需做变量代换

s=\frac{b-a}{2}t+\frac{b+a}{2}

即可利用标准勒让德多项式求出此时的多项式

曲线拟合的最小二乘法

与上面的已知函数逼近不同，曲线拟合问题中需要逼近的目标函数是一系列散点数据构成的表格函数

设已知的数据点为

(t_i,f_i),\ i=1,\cdots,m

则拟合目标为最小化下述表达式

\sum\limits_{i=1}^m[S(t_i)-f_i]^2

一般求解空间为线性空间

S(t)\in\Phi={\rm span}\{\phi_1,\cdots,\phi_n\}

此时的问题称为线性最小二乘

注线性最小二乘要求拟合函数的空间是线性空间，不要求拟合函数是线性函数

最小二乘问题的矩阵描述

令

\mathbf A=\begin{bmatrix} \phi_1(t_1)&\cdots&\phi_n(t_1)\\ \vdots&\ddots&\vdots\\ \phi_1(t_m)&\cdots&\phi_n(t_m) \end{bmatrix}

而我们有观测点

\mathbf f=\begin{bmatrix}f_1\\\vdots\\f_m\end{bmatrix}

则此问题转化为对（一般而言的）超定方程组

\mathbf{Ax}=\mathbf f

求一个最佳近似解使

||\mathbf f-\mathbf{Ax}||_2

最小。此最佳解称为最小二乘解

线性最小二乘的法方程法

回忆函数逼近中的法方程法，我们求解方程 $\mathbf {Gx}=\mathbf b$ ，其中

\mathbf G=\begin{bmatrix}\left<\phi_1,\phi_1\right>&\cdots&\left<\phi_1,\phi_n\right>\\ \vdots&\ddots&\vdots\\ \left<\phi_n,\phi_1\right>&\cdots&\left<\phi_n,\phi_n\right> \end{bmatrix}

\mathbf b=\begin{bmatrix}\left<f,\phi_1\right>\\ \vdots\\ \left<f,\phi_n\right> \end{bmatrix}

然而，对于建立在一系列观测点 $\{t_1,\cdots,t_m\}$ 上的表格函数，我们无法用极限的形式定义内积，但可以定义如下形式的向量内积（表格函数是一种离散向量）：

\left<u,v\right>=\sum\limits_{i=1}^mu(t_i)v(t_i)

则我们根据上面定义的 $\mathbf A$ 和 $\mathbf f$ 可知，在此问题中

\mathbf G=\mathbf A^T\mathbf A

\mathbf b=\mathbf A^T\mathbf f

于是，法方程法即根据矩阵 $\mathbf A$ 解如下方程：

\mathbf A^T\mathbf{Ax}=\mathbf A^T\mathbf f

此时方法是否有唯一解取决于 $\mathbf A^T\mathbf A$ 是否可逆

$\mathbf A^T\mathbf A$ 可逆等价于 $\mathbf A$ 列满秩，亦即表格函数 $\phi_1,\cdots,\phi_n$ 线性无关

注这里的 表格函数 线性无关不等于说对应的 连续函数 线性无关，例如 $\sin t$ 和 $\sin 2t$ 线性无关，但在点集 $\{0,\pi,2\pi\}$ 上，它们是同一表格函数

因此，这里说的表格函数线性无关本质上仍然是向量线性无关

然而，对多项式拟合的情形，只要不同的观测点个数大于多项式的次数，就一定可以找到唯一的最小二乘解

线性最小二乘的正交变换法

法方程法解线性最小二乘的主要缺点有：

$n$ 较大时问题病态
计算 $\mathbf A^T\mathbf A$ 会产生较大误差

如果能够构造表格函数的正交基，则法方程会变成对角方程，可以很精确地解出

然而，表格函数是否正交依赖观测点取值，因此我们无法使用之前讨论过的正交多项式组（例如勒让德多项式）

问题可以如下解决：对 $\mathbf A$ 作QR分解

\mathbf A=\mathbf{QR}

于是

\mathbf f-\mathbf{Ax}=\mathbf Q(\mathbf Q^T\mathbf f-\mathbf {Rx})

由于正交变换不改变向量2-范数

||\mathbf f-\mathbf{Ax}||_2=||\mathbf Q^T\mathbf f-\mathbf {Rx}||_2

于是可以改为最小化右端的2-范数

然而，矩阵 $\mathbf A$ 一般行数（观测点数）远超过列数（基函数数），此时进行完整的QR分解冗余较多，可以改为精简的QR分解：

\mathbf A=\mathbf Q_1\mathbf R_1

其中 $\mathbf Q_1\in\mathbb R^{m\times n}$ 列正交， $\mathbf R_1$ 是上三角方阵，将 $\mathbf A$ 的QR分解中的剩余部分分别标记为 $\mathbf Q_2$ 和 $\mathbf R_2=\mathbf O$ ，则

\mathbf Q^T\mathbf f-\mathbf {Rx}=\begin{bmatrix}\mathbf Q_1^T\\\mathbf Q_2^T\end{bmatrix}\mathbf f-\begin{bmatrix}\mathbf R_1\\\mathbf O\end{bmatrix}\mathbf x=\begin{bmatrix}\mathbf Q_1^T\mathbf f-\mathbf R_1\mathbf x\\\mathbf Q_2^T\mathbf f\end{bmatrix}

于是

||\mathbf Q^T\mathbf f-\mathbf {Rx}||_2^2=||\mathbf Q_1^T\mathbf f-\mathbf {R}_1\mathbf x||_2^2+||\mathbf Q_2^T\mathbf f||_2^2

若 $\mathbf A$ 列满秩（在实际情况中大概率如此！），则 $\mathbf R_1$ 是可逆的，于是上述等式右侧的第一项可以取到 $0$ ，则最小二乘解为 $\mathbf R_1^{-1}\mathbf Q_1^T\mathbf f$

注用这种方法计算的数值比较稳定，更实用

注在实际计算时，计算 $\mathbf Q_1^T\mathbf f$ 可以和计算 $\mathbf A$ 的精简QR分解同步进行（相当于对增广矩阵乘 $\mathbf Q_1^T$ ），而后一步也无需计算 $\mathbf R_1^{-1}$ ，直接利用其上三角特征进行回代法即可

函数插值

插值的基本概念：构造一个函数通过各个离散点

这与表格函数逼近的区别在于：表格函数逼近只要求误差尽量小，不要求函数严格通过每个表格点

因此，函数插值要求插值点的横坐标必须两两不同，否则形成竖直线段无法插值。这在最小二乘问题中得到允许

多项式插值

设插值点为 $(x_0,y_0),\cdots,(x_n,y_n)$ ，则构造多项式

P(x)=a_0+a_1x+\cdots+a_{n-1}x^{n-1}+a_nx^n

只需要求解以下方程

\begin{cases}a_0+a_1x_0+\cdots+a_{n-1}x_0^{n-1}+a_nx_0^n=y_0\\ a_0+a_1x_1+\cdots+a_{n-1}x_1^{n-1}+a_nx_1^n=y_1\\ \cdots\\ a_0+a_1x_n+\cdots+a_{n-1}x_n^{n-1}+a_nx_n^n=y_n\end{cases}

即

\begin{bmatrix}1&x_0&\cdots&x_0^n\\ 1&x_1&\cdots&x_1^n\\ \vdots&\vdots&\ddots&\vdots\\ 1&x_n&\cdots&x_n^n\end{bmatrix} \begin{bmatrix}a_0\\a_1\\\vdots\\a_n\end{bmatrix}= \begin{bmatrix}y_0\\y_1\\\vdots\\y_n\end{bmatrix}

此时的系数矩阵是Vendermonde矩阵，当 $x_0,\cdots,x_n$ 两两不同时可逆

求插值多项式可以采用待定系数法，解上述线性方程组，但这样不便于计算和理论分析

Lagrange插值法

首先考虑只有两个插值点，即 $n=1$ 的情况，此时插值的结果为一个两点式直线方程：

y=\frac{x-x_1}{x_0-x_1}y_0+\frac{x-x_0}{x_1-x_0}y_1

这满足如下的形式：

y=y_0l_0(x)+y_1l_1(x)

其中

l_k(x_i)=\begin{cases}1,\ i=k\\0,\ i\neq k\end{cases}

即：对插值点集，构造一组基函数，使得每个基函数仅在其对应的插值点处取值为 $1$ ，其他插值点处取值为 $0$ ，之后只需要将所有基函数以 $y_i$ 为系数线性组合即可

如何构造出这样的 $l_k$ ？首先，我们需要让它在 $x_k$ 以外的插值点都取 $0$ ：

l_k(x)=C(x-x_0)\cdots(x-x_{k-1})(x-x_{k+1})\cdots(x-x_n)

此时只需要适应 $C$ 的值使 $l_k(x_k)=1$ ，于是

l_k(x)=\frac{(x-x_0)\cdots(x-x_{k-1})(x-x_{k+1})\cdots(x-x_n)}{(x_k-x_0)\cdots(x_k-x_{k-1})(x_k-x_{k+1})\cdots(x_k-x_n)}

如果令

\omega_{n+1}(x)=\prod_{j=0}^n(x-x_j)

则上述表达式简化为

l_k(x)=\frac{\omega_{n+1}(x)}{(x-x_k)\omega_{n+1}'(x_k)}

插值函数为

L_n(x)=\sum\limits_{k=0}^ny_k\frac{\omega_{n+1}(x)}{(x-x_k)\omega_{n+1}'(x_k)}

多项式插值误差估计

假设我们的插值是对一个已知若干插值点的未知函数进行插值来估计未知函数，则此时的函数插值就是一种函数逼近的手段，可以分析其逼近误差：

R_n(x):=f(x)-L_n(x)

定理（Langrange插值余项估计） 设 $f(x)\in C^n[a,b]$ ，且其 $n+1$ 阶导数存在，则

\forall x\in[a,b],\ \exists\xi\in(a,b),\ R_n(x)=\frac{f^{(n+1)}(\xi)}{(n+1)!}\omega_{n+1}(x)

证明

显然，余项 $R_n(x)$ 在所有插值点上都取 $0$ ，因此设

R_n(x)=g(x)\omega_{n+1}(x)

此时对每个不是插值点的 $x$ 构造辅助函数

\phi(t)=R_n(t)-g(x)\omega_{n+1}(t)=f(t)-L_n(t)-g(x)\omega_{n+1}(t)

由于 $x$ 不是插值点，此函数至少有 $n+2$ 个零点：

t=x,x_0,x_1,\cdots,x_n

根据连续性，此时 $\phi'(t)$ 至少有 $n+1$ 个互不相同的零点，类似可知 $\phi^{(n+1)}(t)$ 至少有 $1$ 个零点，记为 $\xi$

于是

\phi^{(n+1)}(\xi)=f^{(n+1)}(\xi)-g(x)(n+1)!=0

注意上述表达式中的 $^{(n+1)}$ 为对 $t$ 求导，因此 $L_n(t)$ 作为 $n$ 次多项式求导后消失，而 $g(x)$ 对 $t$ 是常数， $\omega_{n+1}(t)$ 是最高次项系数为 $1$ 的 $n+1$ 次多项式，因此求导得到 $(n+1)!$

故

g(x)=\frac{f^{(n+1)}(\xi)}{(n+1)!}

这给出了一种估计余项大小的方式，因为 $f^{(n+1)}(\xi)$ 往往可以给出一个上界

注 Lagrange插值公式结构对称，计算相对简单，但当需要增加或减少一个节点时，公式变化较大，必须重新计算，不能实现增量插值

Newton插值

我们试图用增量的方式来构造插值公式，即增加一个插值点，公式应该发生什么变化？

首先考虑只有一个点时

P_0(x)=y_0

此时引入插值点 $(x_1,y_1)$ ，则

P_1(x)=y_0+\frac{y_1-y_0}{x_1-x_0}(x-x_0)

这是一个点斜式直线方程，注意到它和Lagrange插值的结果其实是完全一样的，只是不同的看问题的角度

于是，假设我们已经有 $P_{n-1}(x)$ ，引入 $(x_n,y_n)$ ，则新引入的部分不能改变前 $n$ 个插值点处的值，因此应该形如

P_n(x)=P_{n-1}(x)+c_n(x-x_0)\cdots(x-x_{n-1})

得到Newton插值的公式

N_n(x)=c_0+c_1(x-x_0)+\cdots+c_n(x-x_0)\cdots(x-x_{n-1})

下面我们需要解决的是系数 $c_n$ 的取值

定义（差商） 函数关于一系列互不相等的点可以定义一个差商 $f[x_0,\cdots,x_k]$ ，差商的阶数为差商点数目减一：

f[x_i]=f(x_i)

f[x_0,x_i]=\frac{f(x_i)-f(x_0)}{x_i-x_0}

递归地定义高阶差商：

f[x_0,x_1,\cdots,x_k]=\frac{f[x_0,\cdots,x_{k-2},x_k]-f[x_0,x_1,\cdots,x_{k-1}]}{x_k-x_{k-1}}

定理（差商的非递归公式）

f[x_0,\cdots,x_k]=\sum\limits_{j=0}^k\frac{f(x_j)}{\prod_{l=0,l\neq j}^k(x_j-x_l)}=\sum\limits_{j=0}^k\frac{f(x_j)}{\omega_{k+1}'(x_j)}

证明

对递归基础，显然成立

假设对 $k-1$ 阶差商成立，则

f[x_0,\cdots,x_{k-2},x_k]=\sum\limits_{j=0}^{k-2}\frac{f(x_j)}{(x_j-x_{k})\prod_{l=0,l\neq j}^{k-2}(x_j-x_l)}+\frac{f(x_k)}{(x_k-x_0)\cdots(x_k-x_{k-2})}

f[x_0,\cdots,x_{k-1}]=\sum\limits_{j=0}^{k-1}\frac{f(x_j)}{\prod_{l=0,l\neq j}^{k-1}(x_j-x_l)}

于是

f[x_0,\cdots,x_k]=\frac{1}{x_k-x_{k-1}}(f[x_0,\cdots,x_{k-2},x_k]-f[x_0,\cdots,x_{k-1}])

=\sum\limits_{j=0}^{k-2}\frac{f(x_j)}{(x_k-x_{k-1})(x_j-x_{k})\prod_{l=0,l\neq j}^{k-2}(x_j-x_l)}+\frac{f(x_k)}{(x_k-x_0)\cdots(x_k-x_{k-2})(x_k-x_{k-1})}

-\sum\limits_{j=0}^{k-1}\frac{f(x_j)}{(x_k-x_{k-1})\prod_{l=0,l\neq j}^{k-1}(x_j-x_l)}

=\sum_{j=0}^{k-2}\left(\frac{f(x_j)}{(x_k-x_{k-1})(x_j-x_{k})\prod_{l=0,l\neq j}^{k-2}(x_j-x_l)}-\frac{f(x_j)}{(x_k-x_{k-1})\prod_{l=0,l\neq j}^{k-1}(x_j-x_l)}\right)

-\frac{f(x_{k-1})}{(x_{k-1}-x_0)\cdots(x_{k-1}-x_{k-2})(x_k-x_{k-1})}+\frac{f(x_k)}{(x_k-x_0)\cdots(x_k-x_{k-2})(x_k-x_{k-1})}

=\sum\limits_{j=0}^{k-2}\frac{(x_j-x_{k-1})f(x_j)-(x_j-x_k)f(x_j)}{(x_k-x_{k-1})\prod_{l=0,l\neq j}^{k}(x_j-x_l)}

+\frac{f(x_{k-1})}{(x_{k-1}-x_0)\cdots(x_{k-1}-x_{k-2})(x_{k-1}-x_k)}+\frac{f(x_k)}{(x_k-x_0)\cdots(x_k-x_{k-2})(x_k-x_{k-1})}

=\sum\limits_{j=0}^k\frac{f(x_j)}{\prod_{l=0,l\neq j}^k(x_j-x_l)}

定理（差商的对称性） 根据上面公式可知，差商的自变量顺序可以任意交换

定理（Newton插值的系数） 在Newton插值中，系数 $c_k$ 即为差商 $f[x_0,\cdots,x_k]$

证明略，可以归纳证明

注实际计算Newton插值时，不应总是使用非递归的公式计算，而是应该保存下来已经计算好的低阶差商，然后利用递推公式计算高阶差商，使用类似动态规划的思想

可以使用如下形式的差商表辅助计算：

对任何一个非插值点的 $x$ ，可以得到Newton插值的余项：

f(x)=N_n(x)+f[x,x_0,\cdots,x_n](x-x_0)\cdots(x-x_n)

这与Lagrange余项有潜在的等价关系：

f[x,x_0,\cdots,x_n]=\frac{f^{(n+1)}(\xi)}{(n+1)!}

由于每次增加插值节点增加的多项式系数为 $f[x_0,\cdots,x_k]$ ，可以用此差商来判断是否需要继续增加插值节点

Hermite插值

Hermite插值是一种多项式插值，但其插值条件除了要求函数值，还要求导数值：

P(x_i)=f(x_i),\ P'(x_j)=f'(x_j)

其存在性和唯一性与普通的多项式插值类似，常常使用的是函数值插值点和导数值插值点相同的情形：

P(x_i)=f_i,\ P'(x_i)=f'_i,\ i=0,1,\cdots,n

共 $2n+2$ 个条件，可以唯一确定一个次数不超过 $2n+1$ 的多项式 $H_{2n+1}(x)$

插值多项式的求法与Lagrange插值类似，为若干个基函数的线性组合：

H_{2n+1}(x)=\sum\limits_{j=0}^n[f_j\alpha_j(x)+f'_j\beta_j(x)]

并要求

\alpha_j(x_i)=\begin{cases}1&i=j\\0&i\neq j\end{cases},\quad \alpha_j'(x_i)=0

\beta_j(x_i)=0,\quad\beta_j'(x_i)=\begin{cases}1&i=j\\0&i\neq j\end{cases}

即，将函数值与导数值拆开独立考虑

得到的基函数为

\alpha_j(x)=[1-2(x-x_j)l_j'(x_j)]l_j^2(x)

\beta_j(x)=(x-x_j)l_j^2(x)

多项式插值的稳定性问题

问题本身比较病态，高次插值的数值稳定性差
- 这体现在一个插值点函数值的误差可能影响整个区间
保凸性差：有多余的拐点
Runge现象：次数越高有可能逼近效果越差，在某个范围外可能不收敛

为了防止一个插值点影响非常遥远的区域的函数值，考虑采用**分段插值

分段线性插值

使用三角波函数进行插值，每个三角波函数满足它只在对应插值点附近非零，具体表达式为

l_j(x)=\begin{cases}\frac{x-x_{j-1}}{x_j-x_{j-1}}&x_{j-1}\leq x\leq x_j\\ \frac{x-x_{j+1}}{x_j-x_{j+1}}&x_{j}\leq x\leq x_{j+1}\\ 0&\rm else \end{cases}

得到的插值函数为

I_h(x)=\sum\limits_{j=0}^n y_jl_j(x)

在每个插值区间内相当于一次Lagrange插值，于是误差为

|f(x)-I_h(x)|\leq\frac{M_2}{2}\max_{x_j\leq x\leq x_{j+1}}|(x-x_j)(x-x_{j+1})|\leq \frac{M_2}{8}h^2

其中

M_2=\max_{a\leq x\leq b}|f''(x)|

h=\max_j(x_{j+1}-x_j)

当最大插值区间长度 $h$ 趋于零时，分段线性插值收敛于真实函数

分段Hermite插值

进行完整Hermite插值的复杂度很高且没有太大的意义，而分段线性插值得到的函数具有不可导点

因此可以考虑通过Hermite插值构造一种处处可导、各个插值区间内插值过程独立的分段插值：两点三次Hermite插值

即：对每个插值区间进行 $n=1$ 的整体Hermite插值，利用 $2$ 个函数值和 $2$ 个导数值，得到 $2n+1=3$ 次多项式

可以得到整体一阶连续可微的插值函数 $H_h(x)$

保形分段插值

分段Hermite插值得到的整体一阶连续可微的性质很好，但这需要提供实际上不易得到的真实导数值

因此，当真实导数值未知时，可以人工设定导数值，然后进行Hermite插值，即可得到光滑的曲线

我们采取以下的逻辑设置插值点处的导数：

先计算插值点两侧的割线斜率

d_{k-1}=\frac{f(x_k)-f(x_{k-1})}{x_k-x_{k-1}},\ d_{k}=\frac{f(x_{k+1})-f(x_k)}{x_{k+1}-x_k}

若 $d_{k-1}d_k\leq 0$ ，两侧割线斜率不同，则设定 $f_k'=0$

否则取加权调和平均：

\frac{w_{k-1}+w_k}{f_k'}=\frac{w_{k-1}}{d_{k-1}}+\frac{w_{k}}{d_k}

其中

w_{k-1}=h_{k-1}+2h_k,\ w_k=2h_{k-1}+h_k

对于两端的插值点而言（例如 $x_0$ ）， $f_0'$ 不能直接设为 $d_0$ ，而是与 $d_0,d_1$ 都有关

样条函数插值

三次样条插值函数：要求对每个插值区间是三次多项式且整体二阶连续可微

注传统意义上的样条曲线不要求曲线经过每个控制点，此时不构成插值，也对曲线的形状没有任何限制（允许出现并不是函数形态的曲线）

下面讨论如何确定三次样条插值函数，设其分段形式为

S(x)=\begin{cases}s_0(x)&x\in[x_0,x_1]\\ s_1(x)&x\in[x_1,x_2]\\\cdots\\ s_{n-1}(x)&x\in[x_{n-1},x_n] \end{cases}

于是可以得到以下方程：

s_j(x_j)=f_j,\ s_j(x_j)=f_{j+1},\ j=0,1,\cdots,n-1

s'_{j-1}(x_j)=s'_j(x_j),\ s''_{j-1}(x_j)=s''_j(x_j),\ j=1,2,\cdots,n-1

上述共 $4n-2$ 个方程，但每个 $s_j$ 需要 $4$ 个方程才能确定，因此还缺少两个方程，可以发现方程缺少是来自于边界点的导数缺乏约束，可以使用以下的额外条件来约束：

给定端点一阶导数： $S'(x_0)=f_0',\ S'(x_n)=f_n'$
给定端点二阶导数： $S''(x_0)=f''_0,\ S''(x_n)=f''_n$ ，一般情况要求 $f''_0=f''_n=0$ ，此时称为自然样条插值
设定周期性：即要求 $x_0$ 与 $x_n$ 可以回绕，此时额外条件为 $S'(x_0)=S'(x_n),\ S''(x_0)=S''(x_n)$ ，一般同时要求 $f_0=f_n$
要求前两个插值区间和后两个插值区间的多项式一致

样条插值函数的构造

我们首先设每个插值点处的二次导数为参数：

S''(x_j)=M_j,\quad j=0,1,\cdots,n

由于 $S(x)$ 在 $[x_j,x_{j+1}]$ 上为三次多项式， $S''(x)$ 为一次多项式，故

S''(x)=M_j\left(\frac{x-x_{j+1}}{x_j-x_{j+1}}\right)+M_{j+1}\left(\frac{x-x_{j}}{x_{j+1}-x_{j}}\right)=M_j\frac{x_{j+1}-x}{h_j}+M_{j+1}\frac{x-x_j}{h_j}

于是

S(x)=\frac{M_j}{6h_j}(x_{j+1}-x)^3+\frac{M_{j+1}}{6h_j}(x-x_{j})^3+a_jx+b_j

由于

S(x_j)=f_j,\quad S(x_{j+1})=f_{j+1}

得到方程

\frac{1}{6}M_jh_j^2+a_jx_j+b_j=f_j

\frac{1}{6}M_{j+1}h_j^2+a_jx_{j+1}+b_j=f_{j+1}

解得

a_j=\frac{1}{h_j}\left(f_{j+1}-f_j+\frac{h_j^2}{6}(M_{j+1}-M_j)\right)

b_j=f_j-\frac16M_jh_j^2-a_jx_j

代入原方程，整理得到更对称的形式：

S(x)=\frac{M_j}{6h_j}(x_{j+1}-x)^3+\frac{M_{j+1}}{6h_j}(x-x_{j})^3

+\frac{x}{h_j}\left(f_{j+1}-f_j+\frac{h_j^2}{6}(M_{j+1}-M_j)\right)-\frac{x_j}{h_j}\left(f_{j+1}-f_j+\frac{h_j^2}{6}(M_{j+1}-M_j)\right)

+f_j-\frac16 M_jh_j^2

将所有左端点值合并、右端点值合并，得到

S(x)=\frac{M_j}{6h_j}(x_{j+1}-x)^3+\frac{M_{j+1}}{6h_j}(x-x_{j})^3

+\frac{x_{j+1}-x}{h_j}\left(f_j-\frac16M_jh_j^2\right)+\frac{x-x_j}{h_j}\left(f_{j+1}-\frac16M_{j+1}h_j^2\right)

我们此时已经满足了函数值给定、二阶导数连续的要求，接下来还需要满足一阶导数连续和边界条件，下面以第一类边界条件为例推导

对 $S(x)$ 求导，有

S'(x)=-\frac{M_j}{2h_j}(x_{j+1}-x)^2+\frac{M_{j+1}}{2h_j}(x-x_j)^2+\frac{f_{j+1}-f_j}{h_j}-\frac{M_{j+1}-M_j}{6}h_j,\quad x\in[x_j,x_{j+1}]

则可以求出边界点的两侧导数

S'(x_j+)=-\frac12M_jh_j+\frac{f_{j+1}-f_j}{h_j}-\frac{M_{j+1}-M_j}{6}h_j

=-\frac13M_jh_j-\frac16M_{j+1}h_j+\frac{f_{j+1}-f_j}{h_j}

S'(x_j-)=\frac12M_jh_{j-1}+\frac{f_j-f_{j-1}}{h_{j-1}}-\frac{M_j-M_{j-1}}{6}h_{j-1}

=\frac13M_jh_{j-1}+\frac16M_{j-1}h_{j-1}+\frac{f_{j}-f_{j-1}}{h_{j-1}}

于是

-\frac13M_jh_j-\frac16M_{j+1}h_j+\frac{f_{j+1}-f_j}{h_j}=\frac13M_jh_{j-1}+\frac16M_{j-1}h_{j-1}+\frac{f_{j}-f_{j-1}}{h_{j-1}}

\frac16h_{j-1}M_{j-1}+\frac13(h_{j-1}+h_j)M_j+\frac16h_jM_{j+1}=\frac{f_{j+1}-f_j}{h_j}-\frac{f_{j}-f_{j-1}}{h_{j-1}}

写成更便于理解的形式：

\mu_jM_{j-1}+2M_j+\lambda_jM_{j+1}=d_j,\quad j=1,\cdots,n-1

其中

\mu_j=\frac{h_{j-1}}{h_{j-1}+h_j}\quad\lambda_j=\frac{h_{j}}{h_{j-1}+h_j}

d_j=\frac{6}{h_{j-1}+h_j}\left(\frac{f_{j+1}-f_j}{h_j}-\frac{f_{j}-f_{j-1}}{h_{j-1}}\right)

此时有 $M_0,\cdots,M_n$ 共 $n+1$ 个待定参数，当前有 $n-1$ 个参数，剩余两个参数由边界条件给出，第一种边界条件的情形如下：

S'(x_0+)=f_0'

-\frac13M_0h_0-\frac16M_{1}h_0+\frac{f_{1}-f_0}{h_0}=f_0'

2M_0+M_1=\frac{6}{h_0}\left(\frac{f_1-f_0}{h_0}-f_0'\right)

另一侧的边界条件类似，最终得到方程：

\begin{bmatrix} 2&1\\ \mu_1&2&\lambda_1\\ &\ddots&\ddots&\ddots\\ &&\mu_{n-1}&2&\lambda_{n-1}\\ &&&1&2 \end{bmatrix}\begin{bmatrix} M_0\\M_1\\\vdots\\M_{n-1}\\M_n \end{bmatrix}=\begin{bmatrix} d_0\\d_1\\\vdots\\d_{n-1}\\d_n \end{bmatrix}

这种方程被称为三弯矩方程，解此方程可以得到插值函数中的所有二阶导数（即“弯矩”）参数