插值:求過已知有限個數據點的近似函數
擬合:已知有限個數據點,求近似函數,不要求過已知數據點,只要求在某種意義下在這些點的誤差最小
(一)插值方法
一、拉格朗日多項式插值
1、插值多項式
就是做出一個多項式函數,經過給出的n個節點,並盡可能的接近原函數,將點帶入多項式函數得到一個線性方程組
當系數矩陣滿秩時,有唯一解。而,系數矩陣的行列式為
這是一個范德蒙德行列式,只要各個節點不同時,行列式就不為0,因此可得,一定能夠解出系數方程
還有些指標
2、拉格朗日插值多項式
3、MATLAB實現
function y=lagrange(x0,y0,x)
%n個數據以數組x0,y0輸入,m個插值點以數組x輸入,輸出數組y為m個插值。 n=length(x0);m=length(x); for i=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
二、牛頓插值
1、差商
2、牛頓插值公式
有點就是,多一個數據點,只多一項
PS:
3、差分
4、等距節點插值公式
三、分段線性插值
1、插值多項式的振盪
即如果插值多項式的次數越高,越容易發生振盪,不能很好的擬合。
2、分段線性插值
3、MATLAB實現
四、Hermite插值
1、Hermite插值多項式
2、MATLAB實現
function y-hermite(x0,y0,y1,x); %x0,y0為樣本點數據,y1為導數指,m個插值點以數組x輸入,輸出數組y為m個插值 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=i: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)-x0(k))*(2*a*y0(i)-y(i))+y0(i)); end y(k)=yy; end
五、樣條插值
1、概念
實際中最常用的是k=2和k=3的情況
二、二次樣條函數插值
二次樣條函數有n+2個待定常數,所以要有n+2個條件,才能有唯一解
一定要有一個條件為一階導數
三、三次樣本函數插值
二次樣條函數有n+3個待定常數,所以要有n+3個條件,才能有唯一解
四、三次插值在MATLAB中的實現
部分轉載
1、y=interp1(x0,y0,x,`spline`); % (spline改成linear,則變成線性插值)
2、y=spline(x0,y0,xi);%這個是根據己知的x,y數據,用樣條函數插值出xi處的值。即由x,y的值計算出xi對應的函數值。
3、pp=spline(x0,y0);%是由根據己知的x,y數據,求出它的樣條函數表達式,不過該表達式不是用矩陣直接表示,要求點x`的值,要用函數y`=ppval(pp,x`);
4、pp=csape(x,y,'變界類型','邊界值conds');生成各種邊界條件的三次樣條插值. 其中,(x,y)為數據向量,邊界類型可為:
'complete':給定邊界一階導數,即默認的邊界條件,Lagrange邊界條件
'not-a-knot':非扭結條件,不用給邊界值.
'periodic':周期性邊界條件,不用給邊界值.
'second':給定邊界二階導數.
'variational':自然樣條(邊界二階導數為[0,0]。
五、demo
%轉載= =
clear,clc 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]; t=0:0.05:15; %拉格朗日插值函數 y1=lagrange(x0,y0,t);%調用編寫的lagrange函數 dy1=(lagrange(x0,y0,0.0001)-lagrange(x0,y0,0))/0.0001%x=0處斜率 min1=min(lagrange(x0,y0,13:0.001:15))%13到15最小值 subplot(2,2,1); plot(x0,y0,'ro',t,y1);%畫出曲線 title('拉格朗日插值函數'); %分段線性插值 y2=interp1(x0,y0,t,'spline');%注意區分spline與linear Y2=interp1(x0,y0,t);%默認linear dy2=(interp1(x0,y0,0.0001,'spline')-interp1(x0,y0,0,'spline'))/0.0001%x=0處斜率 min2=min(interp1(x0,y0,13:0.001:15,'spline'))%13到15最小值 subplot(2,2,2); plot(t,y2,'b',t,Y2,'r',x0,y0,'ro');%畫出曲線 title('分段線性插值'); legend('邊條','線性');%顯示圖形圖例 %三次線條插值A y3=spline(x0,y0,t); dy3=(spline(x0,y0,0.0001)-spline(x0,y0,0))/0.0001%x=0處斜率 min3=min(spline(x0,y0,13:0.001:15))%13到15最小值 subplot(2,2,3); plot(x0,y0,'ro',t,y3);%畫出曲線 title('三次線條插值A'); %三次線條插值B pp1=csape(x0,y0);%默認的邊界條件,即給定邊界一階導數 pp2=csape(x0,y0,'second');%給定邊界二階導數 y4=ppval(pp1,t); Y4=ppval(pp2,t); dy4=(ppval(pp1,0.0001)-ppval(pp1,0))/0.0001%x=0處斜率 min4=min(ppval(pp1,13:0.001:15))%13到15最小值 subplot(2,2,4); plot(t,y4,'b',t,Y4,'r',x0,y0,'ro');%畫出曲線 title('三次線條插值B'); legend('一階','二階');
七、二維插值
如果節點是二維的,插值函數是二元函數的話(曲面),我們可以畫出三維的效果圖
1、插值節點為網絡節點
MatLab封裝程序
z=interp2(x0,y0,z0,x,y,'method')
a、x0,y0為節點坐標,z0為n*m維矩陣,表示節點的值
b、x0,y0要求一個為行向量一個為列向量
c、z為矩陣,n=length(y),m=length(x) 因為MATLAB是列優先
e、x,y為插值點坐標,z為函數值
然后如果是三次樣條插值,可以使用命令
pp=csape({x0,y0},z0,conds,valconds),z=fnval(pp,{x,y})
a、x0,y0為節點坐標,z0為n*m維矩陣,表示節點的值
b、x0,y0要求一個為行向量一個為列向量
c、“conds”與一維相同,一般默認
d、x,y為插值點坐標 ,z為函數值
2、demo
clear,clc
%樣本點信息
x=100:100:500;
y=100:100:400;
z=[636 697 624 478 450
698 712 630 478 420
680 674 598 412 400
662 626 552 334 310];
%錄入樣本點信息
pp=csape({x,y},z'); %注意z矩陣的行列所對應的向量
xi=100:10:500;
yi=100:10:400;
cz1=fnval(pp,{xi,yi});
cz2=interp2(x,y,z,xi,yi','spline');
[i,j]=find(cz1==max(max(cz1)))
subplot(1,2,1);
surf(xi,yi,cz1');
shading interp; %插入顏色插值
axis equal;
title('cz1');
subplot(1,2,2);
surf(xi,yi,cz2);
shading interp;
axis equal;
title('cz2');
二、插值節點為散亂節點
1、定義
MATLAB提供了一個函數
zi=griddata(x,y,z,xi,yi)
a、x,y,z為n維向量,就是數據點
b、xi,yi是給定的網格點橫縱坐標(插值點),返回zi的值
2、demo
clear,clc
%樣本點信息
x=[129,140,103.5,88,185.5 195 105 157.5 107.5 77 81 162 162 117.5];
y=[7.5 141.5 23 147 22.5 137.5 85.5 -6.5 -81 3 56.5 -66.5 84 -33.5];
z=-[4 8 6 8 6 8 8 9 9 8 8 9 4 9];
xi=75:200;
yi=-50:150;
zi=griddata(x,y,z,xi,yi','cubic');
subplot(1,2,1);
plot(x,y,'*');
title('xy');
subplot(1,2,2);
mesh(xi,yi,zi);
shading interp;
axis equal;
title('xyz');
