插值:求过已知有限个数据点的近似函数
一维插值
插值多项式
已知函数
f
(
x
)
f(x)
f(x)在区间
[
a
,
b
]
[a,b]
[a,b]上n+1个不同点
x
0
,
x
1
,
.
.
.
,
x
n
x_0,x_1,...,x_n
x0,x1,...,xn处的函数值
y
i
=
f
(
x
i
)
(
i
=
0
,
1
,
.
.
.
,
n
)
y_i=f(x_i)(i=0,1,...,n)
yi=f(xi)(i=0,1,...,n),
求一个至多n次多项式
ϕ
n
(
x
i
)
=
a
0
+
a
1
x
+
.
.
.
+
a
n
x
n
\color{red}\phi_n(x_i)=a_0+a_1x+...+a_nx^n
ϕn(xi)=a0+a1x+...+anxn,使其满足
ϕ
n
(
x
i
)
=
f
(
x
i
)
(
i
=
0
,
1
,
.
.
.
,
n
)
\phi_n(x_i)=f(x_i)(i=0,1,...,n)
ϕn(xi)=f(xi)(i=0,1,...,n),
ϕ
n
(
x
)
\phi_n(x)
ϕn(x)称为插值多项式
插值多项式与被插函数之间的差
R
n
(
x
)
=
f
(
x
)
−
ϕ
n
(
x
)
\color{red}R_n(x)=f(x)-\phi_n(x)
Rn(x)=f(x)−ϕn(x),称为截断误差(插值余项)
f(x)充分光滑时有,
R
n
(
x
)
=
f
(
x
)
−
L
n
(
x
)
=
f
n
+
1
(
ξ
)
(
n
+
1
)
!
w
n
+
1
(
x
)
,
ξ
∈
(
a
,
b
)
\color{red} {R_n(x)=f(x)-L_n(x)=\frac{f^{n+1}(\xi)}{(n+1)!}w_{n+1}(x),\xi\in(a,b)}
Rn(x)=f(x)−Ln(x)=(n+1)!fn+1(ξ)wn+1(x),ξ∈(a,b),其中
w
n
+
1
(
x
)
=
∏
j
=
0
n
(
x
−
x
j
)
\color{red}w_{n+1}(x)=\prod_{j=0}^{n}(x-x_j)
wn+1(x)=∏j=0n(x−xj)
解下列方程即可求出插值多项式
拉格朗日插值
- 原理
拉格朗日插值实际并未求出插值多项式系数,而是构造了一组基函数
- matlab代码
%设n个节点数据以数组x0, y0输入
%m个插值点以数组x输入,输出数组y为m个插值
function y=lagrange(x0,y0,x)
n=length(x0);m=length(x);
for i=1:m
z=x(i);
s=0.0;
for k=1:n
p=1.0;
for j=1:n
if j~=k
p=p*(z-x0(j))/(x0(k)-x0(j));
end
end
s=p*y0(k)+s;
end
y(i)=s;
end
- 示例
牛顿插值
原理
特点:每增加一个节点,插值多项式只增加一项
- 差商
(2) 一阶差商: f [ x i , x j ] = f ( x i ) − f ( x j ) x i − x j \color{red}f[x_i,x_j]=\frac{f(x_i)-f(x_j)}{x_i-x_j} f[xi,xj]=xi−xjf(xi)−f(xj)
(3) 二阶差商: f [ x i , x j , x k ] = f [ x i , x j ] − f [ x j , x k ] x i − x k \color{red}f[x_i,x_j,x_k]=\frac{f[x_i,x_j]-f[x_j,x_k]}{x_i-x_k} f[xi,xj,xk]=xi−xkf[xi,xj]−f[xj,xk]
(4) k阶差商: f [ x 0 , x 1 , . . . , x k ] = f [ x 0 , x 1 , . . . , x k − 1 ] − f [ x 1 , x 2 , . . . , x k ] x 0 − x k \color{red}f[x_0,x_1,...,x_k]=\frac{f[x_0,x_1,...,x_k-1]-f[x_1,x_2,...,x_k]}{x_0-x_k} f[x0,x1,...,xk]=x0−xkf[x0,x1,...,xk−1]−f[x1,x2,...,xk] - Newton插值公式
(1)线性公式 : ϕ 1 ( x ) = f ( x 0 ) + ( x − x 0 ) f [ x 0 , x 1 ] \color{red}\phi_1(x)=f(x_0)+(x-x_0)f[x_0,x_1] ϕ1(x)=f(x0)+(x−x0)f[x0,x1]
(2) N n ( x ) = f ( x 0 ) + ( x − x 0 ) [ x 0 , x 1 ] + . . . + ( x − x 0 ) ( x − x 1 ) . . . ( x − x n − 1 ) f [ x 0 , x 1 , . . . , x n ] \color{red}N_n(x)=f(x_0)+(x-x_0)[x_0,x_1]+...+(x-x_0)(x-x_1)...(x-x_{n-1})f[x_0,x_1,...,x_n] Nn(x)=f(x0)+(x−x0)[x0,x1]+...+(x−x0)(x−x1)...(x−xn−1)f[x0,x1,...,xn]
插值余项: R n ( x ) = ( x − x 0 ) ( x − x 1 ) . . . ( x − x n ) f [ x , x 0 , x 1 , . . . , x n ] = w n + 1 ( x ) f [ x , x 0 , x 1 , . . . , x n ] \color{red}R_n(x)=(x-x_0)(x-x_1)...(x-x_n)f[x,x_0,x_1,...,x_n]=w_{n+1}(x)f[x,x_0,x_1,...,x_n] Rn(x)=(x−x0)(x−x1)...(x−xn)f[x,x0,x1,...,xn]=wn+1(x)f[x,x0,x1,...,xn] - 差分(前向差分)
结点等距(h)
(1) 一阶差分: Δ f k = f k + 1 − f k \color{red}\Delta{f_k}=f_{k+1}-f_k Δfk=fk+1−fk
(2) 二阶差分: Δ 2 f k = Δ f k + 1 − Δ f k \color{red}\Delta^2f_k=\Delta{f_{k+1}}-\Delta{f_k} Δ2fk=Δfk+1−Δfk
(3) m阶差分: Δ m f k = Δ m − 1 f k + 1 − Δ m − 1 f k \color{red}\Delta^mf_k=\Delta^{m-1}{f_{k+1}}-\Delta^{m-1}{f_k} Δmfk=Δm−1fk+1−Δm−1fk - 等距结点插值公式
牛顿前向插值公式: N n ( x 0 + t h ) = f 0 + t Δ f 0 + . . . + t ( t − 1 ) . . . ( t − n + 1 ) n Δ n f 0 \color{red}N_n(x_0+th)=f_0+t\Delta{f_0}+...+\frac{t(t-1)...(t-n+1)}{n}\Delta^nf_0 Nn(x0+th)=f0+tΔf0+...+nt(t−1)...(t−n+1)Δnf0, ( x = x 0 + t h ) \color{red}(x=x_0+th) (x=x0+th)
matlab代码
%newton.m
%求牛顿插值多项式、差商、插值及其误差估计的MATLAB主程序
%输入的量:X是n+1个节点(x_i,y_i)(i = 1,2, ... , n+1)横坐标向量,Y是纵坐标向量,
%x是以向量形式输入的m个插值点,M在[a,b]上满足|f~(n+1)(x)|≤M
%注:f~(n+1)(x)表示f(x)的n+1阶导数
%输出的量:向量y是向量x处的插值,误差限R,n次牛顿插值多项式L及其系数向量C,
%差商的矩阵A
function[y,R,A,C,L] = newton(X,Y,x,M)
n = length(X);
m = length(x);
for t = 1 : m
z = x(t);
A = zeros(n,n);
A(:,1) = Y';
s = 0.0; p = 1.0; q1 = 1.0; c1 = 1.0;
for j = 2 : n
for i = j : n
A(i,j) = (A(i,j-1) - A(i-1,j-1))/(X(i)-X(i-j+1));
end
q1 = abs(q1*(z-X(j-1)));
c1 = c1 * j;
end
C = A(n, n); q1 = abs(q1*(z-X(n)));
for k = (n-1): -1:1
C = conv(C, poly(X(k)));
d = length(C);
C(d) = C(d) + A(k,k);%在最后一维,也就是常数项加上新的差商
end
y(t) = polyval(C,z);
R(t) = M * q1 / c1;
end
L = poly2sym(C);
分段插值
- 原理
将每两个相邻的节点用直线连起来,形成的折线就是分段线性插值函数.
2. matlab实现
y
=
i
n
t
e
r
p
1
(
x
0
,
y
0
,
x
,
′
l
i
n
e
a
r
′
)
\color{red}y=interp1(x0,y0,x,'linear')
y=interp1(x0,y0,x,′linear′)
Hermite插值
- 原理
不仅要求节点处函数同值,而且要求与函数有相同的一阶、二阶、多阶导
已知: y = f ( x ) y=f(x) y=f(x)在n+1个互异结点 x 0 , x 1 , . . . , x n x_0,x_1,...,x_n x0,x1,...,xn上的函数值 y i = f ( x i ) , ( i = 0 , 1 , . . . , n ) y_i=f(x_i),(i=0,1,...,n) yi=f(xi),(i=0,1,...,n)和导数值 y i ′ = f ′ ( x i ) , ( i = 0 , 1 , . . . , n ) y'_i=f'(x_i),(i=0,1,...,n) yi′=f′(xi),(i=0,1,...,n)
求:至多2n+1次的多项式 H ( x ) H(x) H(x),使得 H ( x i ) = y i , H ′ ( x i ) = y i ′ , ( i = 0 , 1 , . . . , n ) H(x_i)=y_i,H'(x_i)=y'_i,(i=0,1,...,n) H(xi)=yi,H′(xi)=yi′,(i=0,1,...,n)
Hermite插值多项式: H ( x ) = ∑ i = 0 n h i [ ( x i − x ) ( 2 a i y i − y i ′ ) + y i ] \color{red}H(x)=\sum_{i=0}^nh_i[(x_i-x)(2a_iy_i-y'_i)+y_i] H(x)=∑i=0nhi[(xi−x)(2aiyi−yi′)+yi]
其中 h i = ∏ j = 0 , j ≠ i n ( x − x j x i − x j ) 2 , a i = ∑ j = 0 , j ≠ i n 1 x i − x j \color{red}h_i=\prod_{j=0,j\neq i}^{n}(\frac{x-x_j}{x_i-x_j})^2,a_i=\sum_{j=0,j\neq i}^n{\frac{1}{x_i-x_j}} hi=∏j=0,j=in(xi−xjx−xj)2,ai=∑j=0,j=inxi−xj1 - matlab代码
n个节点的数据以数组x0(已知点的横坐标),y0(函数值),y1(导数值)输入
% m个插值点以数组x输入,输出数组y为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
- 示例
样条插值
对插值函数的光滑性有较高的要求
原理
给定区间[a,b]的一个划分
Δ
:
a
=
x
0
<
x
1
<
.
.
.
<
x
n
−
1
<
x
n
=
b
\color{red}\Delta: a=x_0<x_1<...<x_{n-1}<x_n=b
Δ:a=x0<x1<...<xn−1<xn=b
如果
s
(
x
)
s(x)
s(x)满足:
(i) 在每个小区间
[
x
i
,
x
i
−
1
]
(
i
=
0
,
1
,
.
.
.
,
n
−
1
)
[x_i,x_{i-1}](i=0,1,...,n-1)
[xi,xi−1](i=0,1,...,n−1)上*s(x)*是
k
k
k次多项式;
(ii)
s
(
x
)
s(x)
s(x)在
[
a
,
b
]
[a,b]
[a,b]上具有
k
−
1
k-1
k−1阶连续导数
则称
s
(
x
)
s(x)
s(x)为关于分划
Δ
\Delta
Δ
的
k
k
k次样条函数
- 折线是一次样条函数
- 二次样条函数
- 三次样条函数
matlab中的样条插值
1. 三次样条曲线有4种边界条件:
(i)自然边界条件,二阶导数在边界处为0,可视为简支梁,是最常用的边界条件。
(ii)第一边界条件,二阶导在边界处已知的边界条件。自然边界条件可视为特例。
(iii)第二边界条件,一阶导在边界处已知的边界条件。
(iv)循环边界条件,一阶导与二阶导在边界处相等的边界条件,适用于封闭或循环的图形。
(1) y = i n t e r p 1 ( x 0 , y 0 , x , ′ s p l i n e ′ ) \color{red}y=interp1(x0,y0,x,'spline') y=interp1(x0,y0,x,′spline′)
(2) y = s p l i n e ( . . . ) \color{red}y=spline(...) y=spline(...)
spline函数只能实现自然边界条件和第二边界条件,可以实现一维或者高维的曲线插值
(3) p p = c s a p e ( x 0 , y 0 , c o n d s ) , y = p p v a l ( p p , x ) \color{red}pp=csape(x0,y0,conds),y=ppval(pp,x) pp=csape(x0,y0,conds),y=ppval(pp,x)
csape函数可以实现三次样条曲线的各种条件,包括第一边界、第二边界、循环边界、混合边界等。可以实现一维至多维的各种情况。
matlab中的一维插值函数interp1
B 样条函数插值方法
二维插值
示 例:
clc,clear
[x,y]=meshgrid(-4:0.8:4); %原始数据
z=peaks(x,y);
[xi,yi]=meshgrid(-4:0.2:4) %插值数据
zi_nearest=interp2(x,y,z,xi,yi,‘nearest’);
zi_linear=interp2(x,y,z,xi,yi);
zi_spline=interp2(x,y,z,xi,yi,‘spline’);
zi_cubic=interp2(x,y,z,xi,yi,‘cubic’);
figure;
hold on;
subplot(321);
surf(x,y,z);
title(‘原始数据’);
subplot(322);
surf(xi,yi,zi_nearest);
title(‘临近点插值’);
subplot(323);
surf(xi,yi,zi_linear);
title(‘线性插值’);
subplot(324);
surf(xi,yi,zi_spline);
title(‘三次样条插值’);
subplot(325);
surf(xi,yi,zi_cubic);
title(‘三次多项式插值’);