关键字：Matalb；拉格朗日插值；线性插值；曲线拟合。
系列文章目录
数值分析——牛顿插值多项式
数值分析——埃尔米特（Hermit）插值
数值分析——分段低次插值
数值分析——三次样条插值
前言
函数、曲线、图形是我们在进行数学分析时的重要工具。但是，在现实生活中，数据中给出的往往是分散的数个数据点，有时我们需要基于数据集拟合曲线，用于分析数据变化的趋势，或者求解未知点的近似值。这种数据分析的方式就是今天要介绍的插值法。
本文对一些复杂的流程与概念进行了简化与省略，如果想详细了解的小伙伴可以参考李庆扬老师的《数值分析 -第5版》的第2章（需要电子版请私戳）
一、插值问题的提出
在许多实际的问题中工程师或者科研人员使用函数
y
=
f
(
x
)
y=f \left( x \right)
y=f(x)来表示某种数量关系，其中相当一部分函数是通过实验或者观测得到的，虽然
f
(
x
)
f \left( x \right)
f(x)在某个区间
[
a
,
b
]
\left[a,b\right]
[a,b]上是存在的，而且某些情况下还是连续的。但往往我们在实际工作中获得的是一张函数表
y
k
=
f
(
x
k
)
(
k
=
0
,
1
,
2
,
3
⋯
,
n
)
y_k=f\left( x_k \right)(k=0,1,2,3\cdots,n)
yk=f(xk)(k=0,1,2,3⋯,n)。根据给定函数表得到一个既能反应函数
f
(
x
)
f\left( x \right)
f(x)特性，又便于计算的插值函数
P
(
x
)
P\left( x \right)
P(x)，就是插值所解决的问题。
设函数
y
=
(
x
)
y=\left( x \right)
y=(x)在区间
[
a
,
b
]
\left[a,b\right]
[a,b]上有定义，且已知在点
a
≤
x
0
<
x
1
<
⋯
<
x
n
−
1
<
x
n
≤
b
a\le x_0<x_1<\cdots<x_{n-1}<x_{n}\le b
a≤x0<x1<⋯<xn−1<xn≤b上的函数值
y
0
,
y
1
,
⋯
,
y
n
−
1
,
y
n
y_0,y_1,\cdots,y_{n-1},y_{n}
y0,y1,⋯,yn−1,yn，存在一个插值函数：
P
(
x
)
=
a
0
+
a
1
x
+
a
2
x
2
+
⋯
+
a
n
−
1
x
n
−
1
+
a
n
x
n
(1)
P\left( x \right)=a_0+a_1x+a_2x^2+\cdots+a_{n-1}x^{n-1}+a_nx^n \tag{1}
P(x)=a0+a1x+a2x2+⋯+an−1xn−1+anxn(1)其中
a
k
,
k
=
0
,
1
,
⋯
,
n
−
1
,
n
a_k,k=0,1,\cdots,n-1,n
ak,k=0,1,⋯,n−1,n是待定常数，且该函数满足：
P
(
x
k
)
=
y
k
,
k
=
0
,
1
,
⋯
,
n
−
1
,
n
(2)
P\left( x_k \right)=y_k,k=0,1,\cdots,n-1,n \tag{2}
P(xk)=yk,k=0,1,⋯,n−1,n(2)上述插值函数就是多项式插值中的插值函数，当然插值函数也可以是三角函数组成的多项式，该插值方法称为三角插值。
二、理论推导
1.线性插值
设有一个待定的线性插值函数
L
(
x
)
L\left( x \right)
L(x)，该插值函数可以表示为两点式：
L
1
(
x
)
=
x
k
+
1
−
x
x
k
+
1
−
x
k
y
k
+
x
−
x
k
x
k
+
1
−
x
k
y
k
+
1
(3)
L_1\left( x \right)=\frac{x_{k+1}-x}{x_{k+1}-x_k}y_k+\frac{x-x_k}{x_{k+1}-x_k}y_{k+1} \tag{3}
L1(x)=xk+1−xkxk+1−xyk+xk+1−xkx−xkyk+1(3)其中
x
k
x_k
xk、
x
k
+
1
x_{k+1}
xk+1分别是原函数上的两个插值点的横坐标，
y
k
y_k
yk、
y
k
+
1
y_{k+1}
yk+1是两个插值点分别对应的函数值。插值函数与原函数的关系如图2.1：
由函数形式可知：
L
1
(
x
k
)
=
y
k
,
L
1
(
x
k
+
1
)
=
y
k
+
1
(4)
L_1\left( x_k \right)=y_k,L_1\left( x_{k+1} \right)=y_{k+1} \tag{4}
L1(xk)=yk,L1(xk+1)=yk+1(4)根据公式(3)设两个基函数：
l
1
(
x
)
=
x
k
+
1
−
x
x
k
+
1
−
x
k
l
2
(
x
)
=
x
−
x
k
x
k
+
1
−
x
k
(5)
\begin{align} l_1\left( x \right) &=\frac{x_{k+1}-x}{x_{k+1}-x_k} \\ l_2\left( x \right) &=\frac{x-x_k}{x_{k+1}-x_k} \end{align} \tag{5}
l1(x)l2(x)=xk+1−xkxk+1−x=xk+1−xkx−xk(5)其中多项式满足
l
1
(
x
k
)
=
1
,
l
1
(
x
k
+
1
)
=
0
l
2
(
x
k
)
=
0
,
l
2
(
x
k
+
1
)
=
1
(6)
\begin{align} l_1\left( x_{k} \right)=1 &, l_1\left( x_{k+1} \right)=0 \\ l_2\left( x_{k} \right)=0 &, l_2\left( x_{k+1} \right)=1 \end{align} \tag{6}
l1(xk)=1l2(xk)=0,l1(xk+1)=0,l2(xk+1)=1(6)公式(3)可以写为：
L
1
(
x
)
=
y
k
l
1
(
x
)
+
y
k
+
1
l
2
(
x
)
(7)
L_1\left( x \right)=y_k l_1\left( x \right)+y_{k+1} l_2\left( x \right) \tag{7}
L1(x)=ykl1(x)+yk+1l2(x)(7)
2.抛物线插值
基于线性插值基函数的设计方式，我们也可以推测抛物线基函数的性质：
l
k
−
1
(
x
k
−
1
)
=
1
l
k
−
1
(
x
k
)
=
0
l
k
−
1
(
x
k
+
1
)
=
0
l
k
(
x
k
−
1
)
=
0
l
k
(
x
k
)
=
1
l
k
(
x
k
+
1
)
=
0
l
k
+
1
(
x
k
−
1
)
=
0
l
k
+
1
(
x
k
)
=
0
l
k
+
1
(
x
k
+
1
)
=
1
(8)
\begin{matrix} l_{k-1}\left( x_{k-1} \right)=1 & l_{k-1}\left( x_{k} \right)=0 & l_{k-1}\left( x_{k+1} \right)=0 \\ l_{k}\left( x_{k-1} \right)=0 & l_{k}\left( x_{k} \right)=1 & l_{k}\left( x_{k+1} \right)=0 \\ l_{k+1}\left( x_{k-1} \right)=0 & l_{k+1}\left( x_{k} \right)=0 & l_{k+1}\left( x_{k+1} \right)=1 \end{matrix} \tag{8}
lk−1(xk−1)=1lk(xk−1)=0lk+1(xk−1)=0lk−1(xk)=0lk(xk)=1lk+1(xk)=0lk−1(xk+1)=0lk(xk+1)=0lk+1(xk+1)=1(8)其中
x
k
−
1
x_{k-1}
xk−1、
x
k
x_k
xk、
x
k
+
1
x_{k+1}
xk+1分别是原函数上的三个插值点的横坐标，
y
k
−
1
y_{k-1}
yk−1、
y
k
y_k
yk、
y
k
+
1
y_{k+1}
yk+1是三个插值点分别对应的函数值。由基函数有两个零点可以设其形式为：
l
k
−
1
(
x
)
=
A
k
−
1
(
x
−
x
k
)
(
x
−
x
k
+
1
)
l
k
(
x
)
=
A
k
(
x
−
x
k
−
1
)
(
x
−
x
k
+
1
)
l
k
+
1
(
x
)
=
A
k
+
1
(
x
−
x
k
−
1
)
(
x
−
x
k
)
(9)
\begin{align} l_{k-1}\left( x \right)&=A_{k-1}\left( x-x_{k} \right)\left( x-x_{k+1} \right) \\ l_{k}\left( x \right)&=A_{k}\left( x-x_{k-1} \right)\left( x-x_{k+1} \right) \\ l_{k+1}\left( x \right)&=A_{k+1}\left( x-x_{k-1} \right)\left( x-x_{k} \right) \end{align} \tag{9}
lk−1(x)lk(x)lk+1(x)=Ak−1(x−xk)(x−xk+1)=Ak(x−xk−1)(x−xk+1)=Ak+1(x−xk−1)(x−xk)(9)其中
A
k
−
1
、
A
k
、
A
k
+
1
A_{k-1}、A_{k}、A_{k+1}
Ak−1、Ak、Ak+1分别是三个待定系数，可以根据
l
k
−
1
(
x
k
−
1
)
=
1
、
l
k
(
x
k
)
=
1
、
l
k
+
1
(
x
k
+
1
)
=
1
l_{k-1}\left( x_{k-1} \right)=1、l_{k}\left( x_{k} \right)=1、l_{k+1}\left( x_{k+1} \right)=1
lk−1(xk−1)=1、lk(xk)=1、lk+1(xk+1)=1确定：
A
k
−
1
(
x
)
=
1
(
x
k
−
1
−
x
k
)
(
x
k
−
1
−
x
k
+
1
)
A
k
(
x
)
=
1
(
x
k
−
x
k
−
1
)
(
x
k
−
x
k
+
1
)
A
k
+
1
(
x
)
=
1
(
x
k
+
1
−
x
k
−
1
)
(
x
k
+
1
−
x
k
)
(10)
\begin{align} A_{k-1}\left( x \right)&=\frac{1}{\left(x_{k-1}-x_{k}\right)\left(x_{k-1}-x_{k+1}\right)} \\ A_{k}\left( x \right)&=\frac{1}{\left(x_{k}-x_{k-1}\right)\left(x_{k}-x_{k+1}\right)} \\ A_{k+1}\left( x \right)&=\frac{1}{\left(x_{k+1}-x_{k-1}\right)\left(x_{k+1}-x_{k}\right)} \end{align} \tag{10}
Ak−1(x)Ak(x)Ak+1(x)=(xk−1−xk)(xk−1−xk+1)1=(xk−xk−1)(xk−xk+1)1=(xk+1−xk−1)(xk+1−xk)1(10)可得抛物线插值的基函数：
l
k
−
1
(
x
)
=
(
x
−
x
k
)
(
x
−
x
k
+
1
)
(
x
k
−
1
−
x
k
)
(
x
k
−
1
−
x
k
+
1
)
l
k
(
x
)
=
(
x
−
x
k
−
1
)
(
x
−
x
k
+
1
)
(
x
k
−
x
k
−
1
)
(
x
k
−
x
k
+
1
)
l
k
+
1
(
x
)
=
(
x
−
x
k
−
1
)
(
x
−
x
k
)
(
x
k
+
1
−
x
k
−
1
)
(
x
k
+
1
−
x
k
)
(11)
\begin{align} l_{k-1}\left( x \right)&=\frac{\left( x-x_{k} \right)\left( x-x_{k+1} \right)}{\left(x_{k-1}-x_{k}\right)\left(x_{k-1}-x_{k+1}\right)} \\ l_{k}\left( x \right)&=\frac{\left( x-x_{k-1} \right)\left( x-x_{k+1} \right)}{\left(x_{k}-x_{k-1}\right)\left(x_{k}-x_{k+1}\right)} \\ l_{k+1}\left( x \right)&=\frac{\left( x-x_{k-1} \right)\left( x-x_{k} \right)}{\left(x_{k+1}-x_{k-1}\right)\left(x_{k+1}-x_{k}\right)} \end{align} \tag{11}
lk−1(x)lk(x)lk+1(x)=(xk−1−xk)(xk−1−xk+1)(x−xk)(x−xk+1)=(xk−xk−1)(xk−xk+1)(x−xk−1)(x−xk+1)=(xk+1−xk−1)(xk+1−xk)(x−xk−1)(x−xk)(11)设抛物线的插值函数：
L
2
(
x
)
=
l
k
−
1
(
x
)
y
k
−
1
+
l
k
(
x
)
y
k
+
l
k
+
1
(
x
)
y
k
+
1
(12)
L_2\left( x \right)=l_{k-1}\left( x \right)y_{k-1}+l_{k}\left( x \right)y_{k}+l_{k+1}\left( x \right)y_{k+1} \tag{12}
L2(x)=lk−1(x)yk−1+lk(x)yk+lk+1(x)yk+1(12)由插值基函数的性质公式(8)得到抛物线插值函数满足：
L
2
(
x
k
−
1
)
=
y
k
−
1
L
2
(
x
k
)
=
y
k
L
2
(
x
k
+
1
)
=
y
k
+
1
(13)
\begin{align} L_2\left( x_{k-1} \right) &=y_{k-1} \\ L_2\left( x_{k} \right) &=y_{k} \\ L_2\left( x_{k+1} \right) &=y_{k+1} \end{align} \tag{13}
L2(xk−1)L2(xk)L2(xk+1)=yk−1=yk=yk+1(13)将式(11)代入式(12)：
L
2
(
x
)
=
(
x
−
x
k
)
(
x
−
x
k
+
1
)
(
x
k
−
1
−
x
k
)
(
x
k
−
1
−
x
k
+
1
)
y
k
−
1
+
(
x
−
x
k
−
1
)
(
x
−
x
k
+
1
)
(
x
k
−
x
k
−
1
)
(
x
k
−
x
k
+
1
)
y
k
+
(
x
−
x
k
−
1
)
(
x
−
x
k
)
(
x
k
+
1
−
x
k
−
1
)
(
x
k
+
1
−
x
k
)
y
k
+
1
(14)
\begin{align} L_2\left( x \right)= &\frac{\left( x-x_{k} \right)\left( x-x_{k+1} \right)}{\left(x_{k-1}-x_{k}\right)\left(x_{k-1}-x_{k+1}\right)}y_{k-1}+\frac{\left( x-x_{k-1} \right)\left( x-x_{k+1} \right)}{\left(x_{k}-x_{k-1}\right)\left(x_{k}-x_{k+1}\right)}y_{k} \\ &+\frac{\left( x-x_{k-1} \right)\left( x-x_{k} \right)}{\left(x_{k+1}-x_{k-1}\right)\left(x_{k+1}-x_{k}\right)}y_{k+1} \end{align} \tag{14}
L2(x)=(xk−1−xk)(xk−1−xk+1)(x−xk)(x−xk+1)yk−1+(xk−xk−1)(xk−xk+1)(x−xk−1)(x−xk+1)yk+(xk+1−xk−1)(xk+1−xk)(x−xk−1)(x−xk)yk+1(14)
3.拉格朗日插值
基于抛物线插值基函数的性质，可以推测对于
n
+
1
n+1
n+1个插值点的插值基函数：
l
j
(
x
k
)
=
{
1
,
k
=
j
,
0
,
k
≠
j
,
j
,
k
=
0
,
1
,
⋯
,
n
−
1
,
n
(15)
l_j\left( x_k \right)= \begin{cases} 1,\quad k=j, \\ 0,\quad k\neq j, \end{cases} \quad j,k=0,1,\cdots,n-1,n \tag{15}
lj(xk)={1,k=j,0,k=j,j,k=0,1,⋯,n−1,n(15)这
n
+
1
n+1
n+1个多项式
l
0
(
x
)
、
l
1
(
x
)
、
⋯
、
l
n
−
1
(
x
)
、
l
n
(
x
)
l_0\left( x \right)、l_1\left( x \right)、\cdots、l_{n-1}\left( x \right)、l_{n}\left( x \right)
l0(x)、l1(x)、⋯、ln−1(x)、ln(x)就是插值点
x
0
、
x
1
、
⋯
、
x
n
−
1
、
x
n
x_{0}、x_{1}、\cdots、x_{n-1}、x_{n}
x0、x1、⋯、xn−1、xn上的n次插值基函数。该基函数的形式确定思路与抛物线插值基函数相同，这里直接给出：
l
k
(
x
)
=
(
x
−
x
0
)
(
x
−
x
1
)
⋯
(
x
−
x
k
−
1
)
(
x
−
x
k
+
1
)
⋯
(
x
−
x
n
−
1
)
(
x
−
x
n
)
(
x
k
−
x
0
)
(
x
k
−
x
1
)
⋯
(
x
k
−
x
k
−
1
)
(
x
k
−
x
k
+
1
)
⋯
(
x
k
−
x
n
−
1
)
(
x
k
−
x
n
)
(16)
l_{k}\left( x \right)= \frac {\left( x-x_0 \right) \left( x-x_1 \right) \cdots \left( x-x_{k-1} \right) \left( x-x_{k+1} \right) \cdots \left( x-x_{n-1} \right) \left( x-x_{n} \right)} {\left( x_k-x_0 \right) \left( x_k-x_1 \right) \cdots \left( x_k-x_{k-1} \right) \left( x_k-x_{k+1} \right) \cdots \left( x_k-x_{n-1} \right) \left( x_k-x_{n} \right)} \tag{16}
lk(x)=(xk−x0)(xk−x1)⋯(xk−xk−1)(xk−xk+1)⋯(xk−xn−1)(xk−xn)(x−x0)(x−x1)⋯(x−xk−1)(x−xk+1)⋯(x−xn−1)(x−xn)(16)设拉格朗日插值函数：
L
n
(
x
)
=
∑
k
=
0
n
y
k
l
k
(
x
)
(17)
L_{n}\left( x \right)=\sum_{k=0}^n y_kl_{k}\left( x \right) \tag{17}
Ln(x)=k=0∑nyklk(x)(17)由插值基函数的性质可知：
L
n
(
x
j
)
=
∑
k
=
0
n
y
k
l
k
(
x
j
)
=
y
j
,
j
=
0
,
1
,
⋯
,
n
−
1
,
n
(18)
L_{n}\left( x_j \right)=\sum_{k=0}^n y_{k} l_{k}\left( x_j \right)=y_{j},\quad j=0,1,\cdots,n-1,n \tag{18}
Ln(xj)=k=0∑nyklk(xj)=yj,j=0,1,⋯,n−1,n(18)为了方便公式的表达，引入记号：
w
n
+
1
(
x
)
=
(
x
−
x
0
)
(
x
−
x
1
)
⋯
(
x
−
x
n
−
1
)
(
x
−
x
n
)
(19)
w_{n+1}\left( x \right)= \left( x-x_0 \right) \left( x-x_1 \right) \cdots \left( x-x_{n-1} \right) \left( x-x_{n} \right) \tag{19}
wn+1(x)=(x−x0)(x−x1)⋯(x−xn−1)(x−xn)(19)对该函数求导得：
w
n
+
1
′
(
x
)
=
(
x
−
x
0
)
(
x
−
x
1
)
⋯
(
x
−
x
k
−
1
)
(
x
−
x
k
+
1
)
⋯
(
x
−
x
n
−
1
)
(
x
−
x
n
)
(20)
w_{n+1}^{'}\left( x \right)= \left( x-x_0 \right) \left( x-x_1 \right) \cdots \left( x-x_{k-1} \right) \left( x-x_{k+1} \right) \cdots \left( x-x_{n-1} \right) \left( x-x_{n} \right) \tag{20}
wn+1′(x)=(x−x0)(x−x1)⋯(x−xk−1)(x−xk+1)⋯(x−xn−1)(x−xn)(20)将公式(19)、(20)代入公式(17)可得：
L
n
(
x
)
=
∑
k
=
0
n
y
k
l
k
(
x
)
=
∑
k
=
0
n
y
k
w
n
+
1
(
x
)
(
x
−
x
k
)
w
n
+
1
′
(
x
k
)
(21)
L_n\left( x \right) =\sum_{k=0}^n y_kl_k\left( x \right) =\sum_{k=0}^n y_k\frac{w_{n+1}\left( x \right)} {\left( x-x_k \right)w_{n+1}^{'}\left( x_k \right)} \tag{21}
Ln(x)=k=0∑nyklk(x)=k=0∑nyk(x−xk)wn+1′(xk)wn+1(x)(21)
三、MATLAB实现
基于线性插值的推导原理与公式编写函数func_LinearInterpolation：
% 利用两个已知点进行一次线性插值
function [aArr]=func_LinearInterpolation(x0,y0,x1,y1)
aArr(2)=-x1*y0/(x0-x1)-x0*y1/(x1-x0);
aArr(1)=y0/(x0-x1)+y1/(x1-x0);
end
基于抛物线插值的推导原理与公式编写函数func_QuadraticInterpolation：
% 利用三个已知点进行二次线性插值
function [aArr]=func_QuadraticInterpolation(x0,y0,x1,y1,x2,y2)
syms x;
y=((x-x1)*(x-x2))*y0/((x0-x1)*(x0-x2))+...
((x-x0)*(x-x2))*y1/((x1-x0)*(x1-x2))+...
((x-x0)*(x-x1))*y2/((x2-x0)*(x2-x1));
aArr=sym2poly(y);
end
基于拉格朗日插值的推导原理与公式编写函数func_QuadraticInterpolation：
% 拉格朗日插值多项式
function [aArr]=func_LagrangeInterpolation(xInArr,yInArr,n)
xInLength=length(xInArr);
yInLength=length(yInArr);
% ************************************************************
% 识别异常
% ************************************************************
% 插值数组不匹配
if xInLength~=yInLength
error('Error:The length of xArr=%d and yArr=%d are not equal in length.',xInLength,yInLength);
elseif xInLength~=n
error('Error:The length of xArr=%d and n=%d are not equal.',xInLength,n);
elseif yInLength~=n
error('Error:The length of yArr=%d and n=%d are not equal.',yInLength,n);
elseif n==1
error('Error:The n should gainer than %d.',n);
end
% 插值数组不为递增数组
for i=1:1:xInLength-1
if xInArr(i)>xInArr(i+1)
error('Error:The xInArr should be a incremental array.');
end
end
% ************************************************************
% 插值计算
% ************************************************************
x=sym('x');
% 计算w
for k1=1:1:n
w(k1)=sym(1);
for k2=1:1:n
if k2==k1
continue;
else
w(k1)=w(k1)*(x-sym(xInArr(k2)));
end
end
end
% 计算w的导数
for k1=1:1:n
dw(k1)=1;
for k2=1:1:n
if k2==k1
continue;
else
dw(k1)=dw(k1)*(xInArr(k1)-xInArr(k2));
end
end
end
L=0;
% 计算拉格朗日插值函数
for k=1:1:n
L=L+yInArr(k)*w(k)/dw(k);
end
L=vpa(L,5);
% 将syms格式的拉格朗日插值函数转换为poly格式
aArr=sym2poly(L);
end
以上的函数都返回的是MATLAB代码中poly格式的多项式系数，比如syms格式的多项式
x
2
+
2
x
+
3
x^2+2x+3
x2+2x+3用poly形式表示为[1,2,3]，在MATLAB中两者可以通过poly2sym或sym2poly进行相互转换：
针对函数 l n ( x ) ln\left( x \right) ln(x)的测试例程如下：
% 测试插值函数——func_LinearInterpolation、func_QuadraticInterpolation、func_LagrangeInterpolation
syms x;
f_x=[-0.916291,-0.693147,-0.510826,-0.356675,-0.223144];
x_s=[0.4 0.5 0.6 0.7 0.8];
L1=func_LinearInterpolation(x_s(2),f_x(2),x_s(3),f_x(3));
L2=func_QuadraticInterpolation(x_s(1),f_x(1),x_s(2),f_x(2),x_s(3),f_x(3));
polyval(L1,0.54)
polyval(L2,0.54)
L3=func_LagrangeInterpolation([x_s(2),x_s(3)],[f_x(2),f_x(3)],2);
L4=func_LagrangeInterpolation([x_s(1),x_s(2),x_s(3)],[f_x(1),f_x(2),f_x(3)],3);
polyval(L3,0.54)
polyval(L4,0.54)
log(0.54)
可得测试结果：
可以看到线性插值与量点的拉格朗日插值具有相同的结果、二次插值与三点的拉格朗日插值具有相同的结果，由公式可知拉格朗日插值分别在两个点与三个点的情况下就是线性插值与抛物线插值。
总结
以上在本文中对拉格朗日插值函数的推导进行了简单的介绍，并提供了MATLAB的实现代码与测试用例。
参考文献
[1]李庆扬,王能超,易大义. 数值分析.第5版[M]. 清华大学出版社, 2008.
扫一扫
5060
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
