建模算法(八)——插值


插值:求過已知有限個數據點的近似函數

擬合:已知有限個數據點,求近似函數,不要求過已知數據點,只要求在某種意義下在這些點的誤差最小

(一)插值方法

一、拉格朗日多項式插值

1、插值多項式

        就是做出一個多項式函數,經過給出的n個節點,並盡可能的接近原函數,將點帶入多項式函數得到一個線性方程組

image

        當系數矩陣滿秩時,有唯一解。而,系數矩陣的行列式為

image

        這是一個范德蒙德行列式,只要各個節點不同時,行列式就不為0,因此可得,一定能夠解出系數方程

        還有些指標

image

image

2、拉格朗日插值多項式

image

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、差商

image

2、牛頓插值公式

image

      有點就是,多一個數據點,只多一項

PS:

image

3、差分

image

image

4、等距節點插值公式

image

三、分段線性插值

1、插值多項式的振盪

       即如果插值多項式的次數越高,越容易發生振盪,不能很好的擬合。

2、分段線性插值

image

image

3、MATLAB實現

image

四、Hermite插值

1、Hermite插值多項式

image

image

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、概念

image

image

        實際中最常用的是k=2和k=3的情況

二、二次樣條函數插值

二次樣條函數有n+2個待定常數,所以要有n+2個條件,才能有唯一解

image

         一定要有一個條件為一階導數

三、三次樣本函數插值

二次樣條函數有n+3個待定常數,所以要有n+3個條件,才能有唯一解

image

image

四、三次插值在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]。

image

image

五、demo

image

%轉載= =
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('一階','二階');

image

七、二維插值

      如果節點是二維的,插值函數是二元函數的話(曲面),我們可以畫出三維的效果圖

1、插值節點為網絡節點

image

       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是列優先

     d、image

     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');

image

二、插值節點為散亂節點

1、定義

image

    MATLAB提供了一個函數

zi=griddata(x,y,z,xi,yi)

    a、x,y,z為n維向量,就是數據點

    b、xi,yi是給定的網格點橫縱坐標(插值點),返回zi的值

2、demo

image

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');
image


免責聲明!

本站轉載的文章為個人學習借鑒使用,本站對版權不負任何法律責任。如果侵犯了您的隱私權益,請聯系本站郵箱yoyou2525@163.com刪除。



 
粵ICP備18138465號   © 2018-2025 CODEPRJ.COM