【数学建模算法】(26)插值和拟合：埃尔米特(Hermite)插值和样条插值

1.埃尔米特(Hermite)插值

1.1.Hermite插值多项式

如果对插值函数，不仅要求它在节点处与函数同值，而且要求它与函数有相同的一阶、二阶甚至更高阶的导数值，这就是 Hermite 插值问题。本节主要讨论在节点处插值函数与函数的值及一阶导数值均相等的 Hermite 插值。

设已知函数 $y=f(x)$ 在 $n+1$ 个互异节点 $x_{0}, x_{1}, \cdots, x_{n}$ 上的函数值 $y_{i}=f\left(x_{i}\right)(i=0,1, \cdots, n)$ 和导数值 $y_{i}^{\prime}=f^{\prime}\left(x_{i}\right) \quad(i=0,1, \cdots, n)$ ，要求一个至多 $2 n+1$ 次的多项式 $H(x)$ ，使得：
$H\left(x_{i}\right)=y_{i} \quad H^{\prime}\left(x_{i}\right)=y_{i}^{\prime} \quad(i=0,1, \cdots, n)$
满足上述条件的多项式 $H(x)$ 称为Hermite多项式。

Hermite 插值多项式为：
$H(x)=\sum_{i=0}^{n} h_{i}\left[\left(x_{i}-x\right)\left(2 a_{i} y_{i}-y_{i}^{\prime}\right)+y_{i}\right]$
其中：
$h_{i}=\prod_{j=0 \atop j \neq i}^{n}\left(\frac{x-x_{j}}{x_{i}-x_{j}}\right)^{2}, \quad a_{i}=\sum_{j=0 \atop j \neq i}^{n} \frac{1}{x_{i}-x_{j}}$

1.2.用Matlab实现Hermite插值

Matlab 中没有现成的 Hermite 插值函数，必须编写一个 M 文件实现插值。
设 $n$ 个节点的数据以数组 $x0$ （已知点的横坐标）， $y0$ （函数值）， $y1$ （导数值）输入（注意 Matlat 的数组下标从 1 开始）， $m$ 个插值点以数组 $x$ 输入，输出数组 $y$ 为 $m$ 个插值。编写一个名为herite.m的M文件：

function y=hermite(x0,y0,y1,x);
n=length(x0);m=length(x);
for k=1:m
    yy=0.0;
    for i=1:n
        h=1.0;
        a=0.0;
        for j=1:n
            if j~=i
                h=h*((x(k)-x0(j))/(x0(i)-x0(j)))^2;
                a=1/(x0(i)-x0(j))+a;
            end
        end
        yy=yy+h*((x0(i)-x(k))*(2*a*y0(i)-y1(i))+y0(i));
    end
    y(k)=yy;
end

2.样条插值

许多工程技术中提出的计算问题对插值函数的光滑性有较高要求，如飞机的机翼外形，内燃机的进、排气门的凸轮曲线，都要求曲线具有较高的光滑程度，不仅要连续，而且要有连续的曲率，这就导致了样条插值的产生。

2.1.样条函数的概念

所谓样条（Spline）本来是工程设计中使用的一种绘图工具，它是富有弹性的细木条或细金属条。绘图员利用它把一些已知点连接成一条光滑曲线（称为样条曲线），并使连接点处有连续的曲率。

数学上将具有一定光滑性的分段多项式称为样条函数。具体地说，给定区间 $[a, b]$ 的一个分划：
$\Delta : \quad a=x_{0}<x_{1}<\cdots<x_{n-1}<x_{n}=b$
如果函数 $s(x)$ 满足：
（1）在每个小区间 $\left[x_{i}, x_{i-1}\right](i=0,1, \cdots, n-1)$ 上 $s(x)$ 是 $k$ 次多项式。
（2） $s(x)$ 在 $[a, b]$ 上具有 $k-1$
阶连续导数。
则称 $s(x)$ 为关于分划 $\Delta$ 的 $k$ 次样条函数，其图形称为 $k$ 次样条函数。 $x_{0}, x_{1}, \cdots, x_{n}$ 称为样条节点， $x_{1}, x_{2}, \cdots, x_{n-1}$ 称为内节点， $x_{0}, x_{n}$ 称为边界点，这类样条函数的全体记做 $S_{P}(\Delta, k)$ ，则 $s(x)$ 是关于分划 $\Delta$ 的 $k$ 次多项式样条函数。 $k$ 次多项式样条函数的一般形式为：
$s_{k}(x)=\sum_{i=0}^{k} \frac{\alpha_{i} x^{i}}{i !}+\sum_{j=1}^{n-1} \frac{\beta_{j}}{k !}\left(x-x_{j}\right)_{+}^{k}$
其中 $\alpha_{i}(i=0,1, \cdots, k)$ 和 $\beta_{j}(j=1,2, \cdots, n-1)$ 均为任意常数，而：
$\left(x-x_{j}\right)_{+}^{k}=\left\{\begin{array}{l}{\left(x-x_{j}\right)^{k}, x \geq x_{j}} \\ {0, \quad x<x_{j}}\end{array}, \quad(j=1,2, \cdots, n-1)\right.$
在实际中最常用的是 $k=2$ 或3的情况，即为二次样条函数和三次样条函数。

二次样条函数：
对于 $[a, b]$ 上的分划 $\Delta : a=x_{0}<x_{1}<\cdots<x_{n}=b$ ，则：
$s_{2}(x)=\alpha_{0}+\alpha_{1} x+\frac{\alpha_{2}}{2 !} x^{2}+\sum_{j=1}^{n-1} \frac{\beta_{j}}{2 !}\left(x-x_{j}\right)_{+}^{2} \in S_{P}(\Delta, 2)$
其中：
$\left(x-x_{j}\right)_{+}^{2}=\left\{\begin{array}{l}{\left(x-x_{j}\right)^{2}, x \geq x_{j}} \\ {0, \quad x<x_{j}}\end{array}, \quad(j=1,2, \cdots, n-1)\right.$

三次样条函数：
对于 $[a, b]$ 上的分划 $\Delta : a=x_{0}<x_{1}<\cdots<x_{n}=b$ ，则：
$s_{2}(x)=\alpha_{0}+\alpha_{1} x+\frac{\alpha_{2}}{2 !} x^{2}+\sum_{j=1}^{n-1} \frac{\beta_{j}}{2 !}\left(x-x_{j}\right)_{+}^{2} \in S_{P}(\Delta, 2)$
其中：
$\left(x-x_{j}\right)_{+}^{3}=\left\{\begin{array}{l}{\left(x-x_{j}\right)^{3}, x \geq x_{j}} \\ {0, \quad x<x_{j}}\end{array}, \quad(j=1,2, \cdots, n-1)\right.$

利用样条函数进行插值，即取插值函数为样条函数，称为样条插值。例如分段线性插值是一次样条插值。下面我们介绍二次、三次样条插值。

2.2.二次样条函数插值

首先，我们注意到 $s_{2}(x) \in S_{P}(\Delta, 2)$ 中含有 $n+2$ 个特定常数，故应需要 $n+2$ + n 个插值条件，因此，二次样条插值问题可分为两类：

问题（1）：
已知插值节点 $x_{i}$ 和相应的函数值 $y_{i}(i=0,1, \cdots, n)$ 以及端点 $x_{0}$ （或 $x_{n}$ ）处的导数值 $y_{0}^{\prime}$ （或 $y_{n}^{\prime}$ ），求 $s_{2}(x) \in S_{p}(\Delta, 2)$ 使得：
$\left\{\begin{array}{l}{s_{2}\left(x_{i}\right)=y_{i}(i=0,1,2, \cdots, n)} \\ {s_{2}^{\prime}\left(x_{0}\right)=y_{0}\left(或s_{n}^{\prime}\left(x_{n}\right)=y_{n}\right)}\end{array}\right.$

问题（2）：
已知插值节点 $x_{i}$ 和相应的导数值 $y_{i}^{\prime}(i=0,1,2, \cdots, n)$ 以及端点 $x_{0}$ （或 $x_{n}$ ）处的函数值 $y_{0}$ （或 $y_{n}$ ），求 $s_{2}(x) \in S_{p}(\Delta, 2)$ 使得：
$\left\{\begin{array}{l}{s_{2}^{\prime}\left(x_{i}\right)=y_{i}^{\prime}(i=0,1,2, \cdots, n)} \\ {s_{2}\left(x_{0}\right)=y_{0}\left(或s_{n}\left(x_{n}\right)=y_{n}\right)}\end{array}\right.$

事实上，可以证明这两类插值问题都是唯一可解的。
对于问题（1）有：
$\left\{\begin{array}{l}{s_{2}\left(x_{0}\right)=\alpha_{0}+\alpha_{1} x_{0}+\frac{1}{2} \alpha_{2} x_{0}^{2}=y_{0}} \\ {s_{2}\left(x_{1}\right)=\alpha_{0}+\alpha_{1} x_{1}+\frac{1}{2} \alpha_{2} x_{1}^{2}=y_{1}} \\ {s_{2}\left(x_{j}\right)=\alpha_{0}+\alpha_{1} x_{j}+\frac{1}{2} \alpha_{2} x_{j}^{2}+\frac{1}{2} \sum_{i=1}^{k-1} \beta_{i}\left(x_{j}-x_{i}\right)^{2}=y_{j} \quad(j=2,3, \cdots, n)} \\ {s_{2}^{\prime}\left(x_{0}\right)=\alpha_{1}+\alpha_{2} x_{0}=y_{0}^{\prime}}\end{array}\right.$

引入记号 $X=\left(\alpha_{0}, \alpha_{1}, \alpha_{2}, \beta_{1}, \cdots, \beta_{n-1}\right)^{T}$ 为未知向量， $C=\left(y_{0}, y_{1}, \cdots, y_{n}, y_{0}^{\prime}\right)$ 为已知向量。

$A=\left[\begin{array}{cccccc}{1} & {x_{0}} & {\frac{1}{2} x_{0}^{2}} & {0} & {\cdots} & {0} \\ {1} & {x_{1}} & {\frac{1}{2} x_{1}^{2}} & {0} & {\cdots} & {0} \\ {1} & {x_{2}} & {\frac{1}{2} x_{2}^{2}} & {\frac{1}{2}\left(x_{2}-x_{1}\right)^{2}} & {\cdots} & {0} \\ {\vdots} & {\vdots} & {\vdots} & {0} & {\vdots} & {\vdots} \\ {1} & {x_{n}} & {\frac{1}{2} x_{n}^{2}} & {\frac{1}{2}\left(x_{n}-x_{1}\right)^{2}} & {\cdots} & {\frac{1}{2}\left(x_{n}-x_{n-1}\right)^{2}} \\ {0} & {1} & {x_{0}} & {0} & {\cdots} & {0}\end{array}\right]$

于是，问题转化为求方程组 $AX=C$ 的解 $X=\left(\alpha_{0}, \alpha_{1}, \alpha_{2}, \beta_{1}, \cdots, \beta_{n-1}\right)^{T}$ 的问题，即可得到二次样条函数 $s_{2}(x)$ 的表达式。
对于问题（2）的情况类似。

2.3.三次样条函数插值

由于 $s_{3}(x) \in S_{p}(\Delta, 3)$ 中含有 $n+3$ 个待定系数。故应需要 $n+3$ 个待定系数，已知插值节点 $x_{i}$ 和相应的函数值 $f\left(x_{i}\right)=y_{i}(i=0,1,2, \cdots, n)$ ，这里提供了 $n+1$ 个条件，还需要2个边界条件。

常用的三次样条函数的边界条件有 3 种类型：
（1） $s_{3}^{\prime}(a)=y_{0}^{\prime}, s_{3}^{\prime}(b)=y_{n}^{\prime}$ 。由这种边界条件建立的样条插值函数称为 $f(x)$ 的完备三次样条插值函数。
特别地， $y_{0}^{\prime}=y_{n}^{\prime}=0$ 时，样条曲线在端点处呈水平状态。
如果 $f^{\prime}(x)$ 不知道，我们可以要求 $s_{3}^{\prime}(x)$ 与 $f^{\prime}(x)$ 在端点处近似相等。这时以 $x_{0}, x_{1}, x_{2}, x_{3}$ 为节点作一个三次Newton插值多项式 $N_{a}(x)$ ，以 $x_{n}, x_{n-1}, x_{n-2}, x_{n-3}$ 作一个三次Newton插值多项式 $N_{b}(x)$ ，要求：
$s^{\prime}(a)=N_{a}^{\prime}(a), s^{\prime}(b)=N_{b}^{\prime}(b)$

（2） $s^{\prime \prime}(a)=y^{\prime \prime}_{0}, s_{3}^{\prime \prime}(b)=y^{\prime \prime}_{3}$ 。特别地 $y_{n}^{\prime \prime}=y_{n}^{\prime \prime}=0$ ，称为自然边界条件。

（3） $s_{3}^{\prime}(a+0)=s_{3}^{\prime}(b-0), s_{3}^{\prime \prime}(a+0)=s_{3}^{\prime \prime}(b-0)$ ，（这里要求） $s_{3}(a+0)=s_{3}(b-0)$ 此条件称为周期条件。

2.4.三次样条插值在Matlab中的实现

在 Matlab 中数据点称之为断点。如果三次样条插值没有边界条件，最常用的方法，就是采用非扭结（not-a-knot）条件。这个条件强迫第 1 个和第 2 个三次多项式的三阶导数相等。对最后一个和倒数第 2 个三次多项式也做同样地处理。

Matlab 中三次样条插值也有现成的函数：

y=interp1(x0,y0,x,'spline')；
y=spline(x0,y0,x)；
pp=csape(x0,y0,conds)，y=ppval(pp,x);

其中 x0,y0 是已知数据点，x 是插值点，y 是插值点的函数值。
对于三次样条插值，我们提倡使用函数 csape，csape 的返回值是 pp 形式，要求插值点的函数值，必须调用函数 ppval。

pp=csape(x0,y0)：使用默认的边界条件，即 Lagrange 边界条件。
pp=csape(x0,y0,conds)中的 conds 指定插值的边界条件，其值可为：

conds	作用
'complete'	边界为一阶导数，即默认的边界条件
'not-a-knot'	非扭结条件
'periodic'	周期条件
'second'	边界为二阶导数，二阶导数的值[0,0]

对于一些特殊的边界条件，可以通过 conds 的一个 $1 \times 2$ 矩阵来表示，conds 元素的取值为1，2。此时，使用命令：

pp=csape(x0,y0_ext,conds)

其中 y0_ext=[left, y0, right]，这里 left 表示左边界的取值，right 表示右边界的取值。
conds(i)=j 的含义是给定端点 $i$ 的 $j$ 阶导数，即 conds 的第一个元素表示左边界的条件，第二个元素表示右边界的条件，conds=[2,1]表示左边界是二阶导数，右边界是一阶导数，对应的值由 left 和 right 给出。

例：机床加工
待加工零件的外形根据工艺要求由一组数据 $(x, y)$ 给出（在平面情况下），用程控铣床加工时每一刀只能沿 $x$ 方向和 $y$ 方向走非常小的一步，这就需要从已知数据得到加工所要求的步长很小的 $(x, y)$ 坐标。
下表给出的数据 $x, y$ 数据位于机翼断面的下轮廓线上，假设需要得到 $x$ 坐标每改变0.1时的 $y$ 坐标。试完成加工所需数据，画出曲线，并求出 $x=0$ 处的曲线斜率和 $13 \leq x \leq 15$ 范围内 $y$ 的最小值。
数据表：

$x$	0	3	5	7	9	11	12	13	14	15
$y$	0	1.2	1.7	2.0	2.1	2.0	1.8	1.2	1.0	1.6

利用Matlab编程，使用Lagrange，分段线性和三次样条三种插值方法计算。

clc,clear
x0=[0 3 5 7 9 11 12 13 14 15];
y0=[0 1.2 1.7 2.0 2.1 2.0 1.8 1.2 1.0 1.6];
x=0:0.1:15;
y1=lagrange(x0,y0,x); %调用前面编写的Lagrange插值函数
y2=interp1(x0,y0,x);
y3=interp1(x0,y0,x,'spline');
pp1=csape(x0,y0); y4=ppval(pp1,x);
pp2=csape(x0,y0,'second'); y5=ppval(pp2,x);
fprintf('比较一下不同插值方法和边界条件的结果:\n')
fprintf('x y1 y2 y3 y4 y5\n')
xianshi=[x',y1',y2',y3',y4',y5'];
fprintf('%f\t%f\t%f\t%f\t%f\t%f\n',xianshi')
subplot(2,2,1), plot(x0,y0,'+',x,y1), title('Lagrange')
subplot(2,2,2), plot(x0,y0,'+',x,y2), title('Piecewise linear')
subplot(2,2,3), plot(x0,y0,'+',x,y3), title('Spline1')
subplot(2,2,4), plot(x0,y0,'+',x,y4), title('Spline2')
dyx0=ppval(fnder(pp1),x0(1)) %求x=0处的导数
ytemp=y3(131:151);
index=find(ytemp==min(ytemp));
xymin=[x(130+index),ytemp(index)]

（Lagrange函数请参见我的简书文章【数学建模算法】(23)插值和拟合：拉格朗日插值）

最后编辑于：2019.08.21 20:56:13

人面猴
序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
沈念sama阅读 158,847评论 4赞 362
死咒
序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
沈念sama阅读 67,208评论 1赞 292
救了他两次的神仙让他今天三更去死
文/潘晓璐我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
开封第一讲书人阅读 108,587评论 0赞 243
道士缉凶录：失踪的卖姜人
文/不坏的土叔我叫张陵，是天一观的道长。经常有香客问我，道长，这世上最难降的妖魔是什么？我笑而不...
开封第一讲书人阅读 43,942评论 0赞 205
港岛之恋（遗憾婚礼）
正文为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
茶点故事阅读 52,332评论 3赞 287
恶毒庶女顶嫁案：这布局不是一般人想出来的
文/花漫我一把揭开白布。她就那样静静地躺着，像睡着了一般。火红的嫁衣衬着肌肤如雪。梳的纹丝不乱的头发上，一...
开封第一讲书人阅读 40,587评论 1赞 218
城市分裂传说
那天，我揣着相机与录音，去河边找鬼。笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。我是一名探鬼主播，决...
沈念sama阅读 31,853评论 2赞 312
双鸳鸯连环套：你想象不到人心有多黑
文/苍兰香墨我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
开封第一讲书人阅读 30,568评论 0赞 198
万荣杀人案实录
序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
沈念sama阅读 34,273评论 1赞 242
护林员之死
正文独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋以下内容为张勋视角年9月15日...
茶点故事阅读 30,542评论 2赞 246
白月光启示录
正文我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
茶点故事阅读 32,033评论 1赞 260
活死人
序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
沈念sama阅读 28,373评论 2赞 253
日本核电站爆炸内幕
正文年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
茶点故事阅读 33,031评论 3赞 236
男人毒药：我在死后第九天来索命
文/蒙蒙一、第九天我趴在偏房一处隐蔽的房顶上张望。院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
开封第一讲书人阅读 26,073评论 0赞 8
一桩弑父案，背后竟有这般阴谋
文/苍兰香墨我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。一阵脚步声响...
开封第一讲书人阅读 26,830评论 0赞 195
情欲美人皮
我被黑心中介骗来泰国打工，没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。一个月前我还...
沈念sama阅读 35,628评论 2赞 274
代替公主和亲
正文我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
茶点故事阅读 35,537评论 2赞 269