第4講插值(簡化)

1、后勤工程學(xué)院數(shù)學(xué)教研室后勤工程學(xué)院數(shù)學(xué)教研室MATLAB與數(shù)學(xué)實(shí)驗(yàn)與數(shù)學(xué)實(shí)驗(yàn)數(shù)學(xué)教研室數(shù)學(xué)教研室 劉淑環(huán)劉淑環(huán) 第四講第四講 插值插值實(shí)例：海底探測問題實(shí)例：海底探測問題某公司用聲納對海底進(jìn)行測試，在某公司用聲納對海底進(jìn)行測試，在5 55 5海里的海里的坐標(biāo)點(diǎn)上測得海底深度的值，希望通過這些有坐標(biāo)點(diǎn)上測得海底深度的值，希望通過這些有限的數(shù)據(jù)了解更多處的海底情況。并繪出較細(xì)限的數(shù)據(jù)了解更多處的海底情況。并繪出較細(xì)致的海底曲面圖。致的海底曲面圖。在應(yīng)用領(lǐng)域中，由有限個已知數(shù)據(jù)點(diǎn)，構(gòu)造一個在應(yīng)用領(lǐng)域中，由有限個已知數(shù)據(jù)點(diǎn)，構(gòu)造一個解析表達(dá)式，由此計(jì)算數(shù)據(jù)點(diǎn)之間的函數(shù)值，稱解析表達(dá)式，由此計(jì)算數(shù)據(jù)點(diǎn)之

2、間的函數(shù)值，稱之為插值。之為插值。3實(shí)驗(yàn)?zāi)康膶?shí)驗(yàn)?zāi)康膶?shí)驗(yàn)內(nèi)容實(shí)驗(yàn)內(nèi)容2、掌握用、掌握用MATLAB求解插值問題。求解插值問題。1、了解插值的基本內(nèi)容。、了解插值的基本內(nèi)容。11一維插值一維插值22二維插值二維插值33實(shí)驗(yàn)作業(yè)實(shí)驗(yàn)作業(yè)4一一 維維 插插 值值一、插值的定義一、插值的定義二、插值的方法二、插值的方法三、用三、用Matlab解插值問題解插值問題5一維插值的定義一維插值的定義已知已知 n+1個節(jié)點(diǎn)個節(jié)點(diǎn), 1 , 0(),(njyxjj其中其中jx互不相同，不妨設(shè)互不相同，不妨設(shè)),10bxxxan求任一插值點(diǎn)求任一插值點(diǎn))(*jxx 處的插值處的插值.*y 0 x1xnx0y1y節(jié)點(diǎn)

3、可視為由節(jié)點(diǎn)可視為由)(xgy 產(chǎn)生，產(chǎn)生，g表達(dá)式復(fù)雜表達(dá)式復(fù)雜,或無封閉形式或無封閉形式,或未知或未知.。*x*y6 構(gòu)造一個構(gòu)造一個(相對簡單的相對簡單的)函數(shù)函數(shù)),(xfy 通過全部節(jié)點(diǎn)通過全部節(jié)點(diǎn), 即即), 1 ,0()(njyxfjj再用再用)(xf計(jì)算插值，即計(jì)算插值，即).(*xfy 0 x1xnx0y1y*x*y7 稱為稱為拉格朗日插值基函數(shù)拉格朗日插值基函數(shù)。n0iiiny)x(L)x(P 已知函數(shù)已知函數(shù)f(x)在在n+1個點(diǎn)個點(diǎn)x0,x1,xn處的函數(shù)值為處的函數(shù)值為 y0,y1,yn 。求一。求一n次多項(xiàng)式函數(shù)次多項(xiàng)式函數(shù)Pn(x)，使其滿足：，使其滿足： Pn(

4、xi)=yi,i=0,1,n.解決此問題的拉格朗日插值多項(xiàng)式公式如下解決此問題的拉格朗日插值多項(xiàng)式公式如下其中其中Li(x) 為為n次多項(xiàng)式：次多項(xiàng)式：)xx()xx)(xx()xx)(xx()xx()xx)(xx()xx)(xx()x(Lni1ii1ii1i0in1i1i10i拉格朗日拉格朗日(Lagrange)插值插值8拉格朗日拉格朗日(Lagrange)插值插值特別地特別地:兩點(diǎn)一次兩點(diǎn)一次(線性線性)插值多項(xiàng)式插值多項(xiàng)式: 101001011yxxxxyxxxxxL三點(diǎn)二次三點(diǎn)二次(拋物拋物)插值多項(xiàng)式插值多項(xiàng)式: 2120210121012002010212yxxxxxxxxyxxx

5、xxxxxyxxxxxxxxxL拉格朗日插值的拉格朗日插值的MatlabMatlab實(shí)現(xiàn)實(shí)現(xiàn)MatlabMatlab中沒有現(xiàn)成的拉格朗日插值函數(shù)，必須編寫中沒有現(xiàn)成的拉格朗日插值函數(shù)，必須編寫一個一個M M文件實(shí)現(xiàn)拉格朗日插值。文件實(shí)現(xiàn)拉格朗日插值。 設(shè)設(shè)n n個節(jié)點(diǎn)數(shù)據(jù)以數(shù)組個節(jié)點(diǎn)數(shù)據(jù)以數(shù)組x0, y0 x0, y0輸入（注意輸入（注意MatlabMatlab的的數(shù)組下標(biāo)從數(shù)組下標(biāo)從1 1開始），開始），m m個插值點(diǎn)以數(shù)組個插值點(diǎn)以數(shù)組x x輸入，輸輸入，輸出數(shù)組出數(shù)組y y為為m m個插值。個插值。% %編寫一個名為編寫一個名為lagrange.mlagrange.m的的M M文件：文件

6、：function y=lagrange(x0,y0,x);function y=lagrange(x0,y0,x); n=length(x0);m=length(x); n=length(x0);m=length(x); for i=1:m for i=1:m z=x(i); z=x(i); s=0.0; s=0.0; for k=1:n for k=1:n p=1.0; p=1.0; for j=1:n for j=1:n if j=k if j=k p=p p=p* *(z-x0(j)/(x0(k)-x0(j);(z-x0(j)/(x0(k)-x0(j); end end end end

7、 s=p s=p* *y0(k)+s;y0(k)+s; end end y(i)=s; y(i)=s; end end1155,11)(2xxxg 采用拉格朗日多項(xiàng)式插值：選取不同插值采用拉格朗日多項(xiàng)式插值：選取不同插值節(jié)點(diǎn)個數(shù)節(jié)點(diǎn)個數(shù)n+1，其中，其中n為插值多項(xiàng)式的次數(shù)，當(dāng)為插值多項(xiàng)式的次數(shù)，當(dāng)n分別取分別取2,4,6,8,10時(shí)，繪出插值結(jié)果圖形時(shí)，繪出插值結(jié)果圖形.例例12x0=-5:2:5;x0=-5:2:5;y0=1./(1+x0.2);y0=1./(1+x0.2);x=-5:0.1:5;x=-5:0.1:5;y=lagrange(x0,y0,x);y=lagrange(x0,y0

8、,x);y1=1./(1+x.2);y1=1./(1+x.2);plot(x,y,-r)plot(x,y,-r)hold onhold onplot(x,y1,-b)plot(x,y1,-b)-5-4-3-2-1012345-0.200.20.40.60.811.2插值的目的是函數(shù)逼近的一種手段。函數(shù)逼近的目的是得到一插值的目的是函數(shù)逼近的一種手段。函數(shù)逼近的目的是得到一個數(shù)學(xué)問題的精確解或足夠精確地解。那么，是否插值多項(xiàng)式個數(shù)學(xué)問題的精確解或足夠精確地解。那么，是否插值多項(xiàng)式的次數(shù)越高（即插值節(jié)點(diǎn)越多），越能夠達(dá)到這個目的？下面的次數(shù)越高（即插值節(jié)點(diǎn)越多），越能夠達(dá)到這個目的？下面看看n=1

9、0n=10的情況。的情況。55,11)(2xxxgx0=-5:2:5;x0=-5:2:5;y0=1./(1+x0.2);y0=1./(1+x0.2);x=-5:0.1:5;x=-5:0.1:5;y=lagrange(x0,y0,x);y=lagrange(x0,y0,x);y1=1./(1+x.2);y1=1./(1+x.2);plot(x,y,-r)plot(x,y,-r)hold onhold onplot(x,y1,-b)plot(x,y1,-b)x0=-5:1:5;x0=-5:1:5;y0=1./(1+x0.2);y0=1./(1+x0.2);x=-5:0.1:5;x=-5:0.1:5

10、;y=lagrange(x0,y0,x);y=lagrange(x0,y0,x);y1=1./(1+x.2);y1=1./(1+x.2);plot(x,y,-g)plot(x,y,-g)hold onhold onplot(x,y1,-b)plot(x,y1,-b)n=5n=5n=10n=10-5-4-3-2-1012345-0.500.511.52n=5和 n=10時(shí) 的 插 值 效 果 對 比 圖 n=10data2n=5data415由圖看出，在由圖看出，在0 0附近插值效果較好，即余項(xiàng)較小。附近插值效果較好，即余項(xiàng)較小。另一種現(xiàn)象是插值多項(xiàng)式隨節(jié)點(diǎn)增多而振動增大另一種現(xiàn)象是插值多項(xiàng)式隨

11、節(jié)點(diǎn)增多而振動增大; ;雖然雖然P P1010(x)(x)在在1111個插值節(jié)點(diǎn)上與所逼近函數(shù)個插值節(jié)點(diǎn)上與所逼近函數(shù)f(x)f(x)有相同值，有相同值，但整體逼近效果較差，越靠近端點(diǎn)，逼近的效果就越差。但整體逼近效果較差，越靠近端點(diǎn)，逼近的效果就越差。這種現(xiàn)象稱為龍格（這種現(xiàn)象稱為龍格（RungeRunge）現(xiàn)象。）現(xiàn)象。這也是高次插值多項(xiàng)式的病態(tài)現(xiàn)象。這也是高次插值多項(xiàng)式的病態(tài)現(xiàn)象。16x0=-1:0.5:1;x0=-1:0.5:1;y0=1./(1+25y0=1./(1+25* *x0.2);x0.2);x1=-1:0.2:1;x1=-1:0.2:1;x2=-1:0.1:1;x2=-1:

12、0.1:1;y1=lagrange(x0,y0,x1);y1=lagrange(x0,y0,x1);y2=lagrange(x0,y0,x2);y2=lagrange(x0,y0,x2);y3=1./(1+25y3=1./(1+25* *x2.2)x2.2)plot(x0,y0,-r)plot(x0,y0,-r)hold onhold onplot(x1,y1,-b)plot(x1,y1,-b)hold onhold onplot(x2,y2,-g)plot(x2,y2,-g)hold onhold onplot(x2,y3,-y)plot(x2,y3,-y)-1-0.8-0.6-0.4-0.

13、200.20.40.60.81-0.4-0.200.20.40.60.81 步 長 0.5步 長 0.2步 長 0.1原 函 數(shù)從圖中看出，在靠近從圖中看出，在靠近-1-1或或1 1時(shí)，時(shí)，余項(xiàng)會隨余項(xiàng)會隨n n的增大而增大的增大而增大18上述現(xiàn)象和定理告訴我們，并不是插值多項(xiàng)式的次數(shù)越高上述現(xiàn)象和定理告訴我們，并不是插值多項(xiàng)式的次數(shù)越高（即插值節(jié)點(diǎn)越多）精度越高。從數(shù)值計(jì)算上可解釋為高（即插值節(jié)點(diǎn)越多）精度越高。從數(shù)值計(jì)算上可解釋為高次插值多項(xiàng)式的計(jì)算會帶來舍入誤差的增大。從而引起計(jì)次插值多項(xiàng)式的計(jì)算會帶來舍入誤差的增大。從而引起計(jì)算失真。因此，實(shí)際應(yīng)用做插值時(shí)一般只用一次、二次最算失真。因

14、此，實(shí)際應(yīng)用做插值時(shí)一般只用一次、二次最多用三次插值多項(xiàng)式。那么如何提高插值精度呢？采用分多用三次插值多項(xiàng)式。那么如何提高插值精度呢？采用分段插值是一種辦法段插值是一種辦法19練習(xí)實(shí)驗(yàn)練習(xí)實(shí)驗(yàn)2021分段線性插值分段線性插值 xjxj-1xj+1x0 xnxoy可以看出，如果以每個可以看出，如果以每個小區(qū)間的端點(diǎn)為插值點(diǎn)小區(qū)間的端點(diǎn)為插值點(diǎn)的話，那么相鄰區(qū)間的的話，那么相鄰區(qū)間的兩個插值函數(shù)在節(jié)點(diǎn)處兩個插值函數(shù)在節(jié)點(diǎn)處將保持連續(xù)。即將保持連續(xù)。即P(xP(x）是以連續(xù)函數(shù)，但在節(jié)是以連續(xù)函數(shù)，但在節(jié)點(diǎn)處一階二階導(dǎo)數(shù)并不點(diǎn)處一階二階導(dǎo)數(shù)并不一定連續(xù)。一定連續(xù)。22其它,0,)()()(11111

15、10jjjjjjjjjjjnjjjnxxxxxxxxxxxxxxxlxlyxL計(jì)算量與計(jì)算量與n n無關(guān)無關(guān); ;n n越大，誤差越小越大，誤差越小. .nnnxxxxgxL0),()(lim xjxj-1xj+1x0 xnxoy23求下列數(shù)據(jù)的多項(xiàng)式插值函數(shù)求下列數(shù)據(jù)的多項(xiàng)式插值函數(shù) x=0 1 2 3 4 5 6 7 8 9 10; x=0 1 2 3 4 5 6 7 8 9 10; y=1 2.3 2.1 2 4.6 4.7 4.3 8.1 9.2 9.8 10.3; y=1 2.3 2.1 2 4.6 4.7 4.3 8.1 9.2 9.8 10.3; plot(x,y, plot(x

16、,y,* *)012345678910123456789101124求下列數(shù)據(jù)的多項(xiàng)式插值函數(shù)求下列數(shù)據(jù)的多項(xiàng)式插值函數(shù) x=0 1 2 3 4 5 6 7 8 9 10; x=0 1 2 3 4 5 6 7 8 9 10; y=1 2.3 2.1 2 4.6 4.7 4.3 8.1 9.2 9.8 10.3; y=1 2.3 2.1 2 4.6 4.7 4.3 8.1 9.2 9.8 10.3; plot(x,y, plot(x,y,* *)0123456789101234567891011 plot(x,y, plot(x,y,* *,x,y),x,y)思考：如何作分段的多項(xiàng)式插值函數(shù)？思

17、考：如何作分段的多項(xiàng)式插值函數(shù)？25 x0=-5:2.5:5;%n=4 x0=-5:2.5:5;%n=4 y0=5./(1+x0.2); y0=5./(1+x0.2); x=-5:0.1:5; x=-5:0.1:5; y1=5./(1+x.2); y1=5./(1+x.2); plot(x0,y0, plot(x0,y0,* *,x0,y0,x,y1,r),x0,y0,x,y1,r)-5-4-3-2-101234500.511.522.533.544.55 散 點(diǎn) 圖分 段 線 性 插 值f(x）26 x0=-5:1:5;%n=10 x0=-5:1:5;%n=10 y0=5./(1+x0.2)

18、; y0=5./(1+x0.2); x=-5:0.1:5; x=-5:0.1:5; y1=5./(1+x.2); y1=5./(1+x.2); plot(x0,y0, plot(x0,y0,* *,x0,y0,x,y1,r),x0,y0,x,y1,r)-5-4-3-2-101234500.511.522.533.544.55 散 點(diǎn) 圖分 段 線 性 插 值f(x）x1=0:3;y1=1 2.3 2.1 2;x0=0:0.1:3;x1=0:3;y1=1 2.3 2.1 2;x0=0:0.1:3;y0=lagrange(x1,y1,x0);y0=lagrange(x1,y1,x0);x2=3:6

19、;y2=2 4.6 4.7 4.3;x20=3:0.1:6;x2=3:6;y2=2 4.6 4.7 4.3;x20=3:0.1:6;y20=lagrange(x2,y2,x20);y20=lagrange(x2,y2,x20);x3=6:8;y3=4.3 8.1 9.2;x30=6:0.1:8;x3=6:8;y3=4.3 8.1 9.2;x30=6:0.1:8;y30= lagrange(x3,y3,x30);y30= lagrange(x3,y3,x30);x4=8:10;y4=9.2 9.8 10.3;x40=8:0.1:10;x4=8:10;y4=9.2 9.8 10.3;x40=8:0

20、.1:10;y40= lagrange(x4,y4,x40);y40= lagrange(x4,y4,x40);plot(x0,y0,-r,x1,y1,plot(x0,y0,-r,x1,y1,* *,x20,y20,-r,x2,y2,x20,y20,-r,x2,y2,* *)hold onhold onplot(x30,y30,-r,x3,y3,plot(x30,y30,-r,x3,y3,* *, x40,y40,-r,x4,y4, x40,y40,-r,x4,y4,* *)29比分段線性插值更光滑。比分段線性插值更光滑。xyxi-1 xiab 在數(shù)學(xué)上，光滑程度的定量描述是：函數(shù)在數(shù)學(xué)上，光

21、滑程度的定量描述是：函數(shù)(曲曲線線)的的k階導(dǎo)數(shù)存在且連續(xù)，則稱階導(dǎo)數(shù)存在且連續(xù)，則稱該曲線具有該曲線具有k階光階光滑性滑性。 光滑性的階次越高，則越光滑。是否存在較低光滑性的階次越高，則越光滑。是否存在較低次的分段多項(xiàng)式達(dá)到較高階光滑性的方法？三次次的分段多項(xiàng)式達(dá)到較高階光滑性的方法？三次樣條插值就是一個很好的例子。樣條插值就是一個很好的例子。三次樣條插值三次樣條插值30 三次樣條插值三次樣條插值, 1,),()(1nixxxxsxSiii,)()3), 1 ,0()()2), 1()()10223niiiiiiixxCxSniyxSnidxcxbxaxs) 1, 1()()(),()(),

22、()(111 nixsxsxsxsxsxsiiiiiiiiiiii自然邊界條件）(0)()()40 nxSxS)(,)4)3)2xSdcbaiiiig g( (x x) )為被插值函數(shù)為被插值函數(shù)。lim ( )( )nS xg x31例例66,11)(2xxxg用三次樣條插值選取用三次樣條插值選取11個基點(diǎn)計(jì)算插值個基點(diǎn)計(jì)算插值3266,11)(2xxxgx=-6:6;y=1./(1+x.2);%插值節(jié)點(diǎn)共插值節(jié)點(diǎn)共11個個h=-6:0.1:6;%設(shè)置新的被插值點(diǎn)設(shè)置新的被插值點(diǎn)y1=interp1(x,y,h,spline); %用三次樣條方法求得相應(yīng)插值用三次樣條方法求得相應(yīng)插值plot

23、(x,y,+,x,y) %作原函數(shù)圖形作原函數(shù)圖形xlabel(x),ylabel(y) %x軸軸y軸標(biāo)簽軸標(biāo)簽hold on %繼續(xù)繪圖繼續(xù)繪圖plot(h,y1,*,h,y1) %作插值點(diǎn)圖形作插值點(diǎn)圖形-6-4-2024600.10.20.30.40.50.60.70.80.91-6-4-2024600.10.20.30.40.50.60.70.80.91xyh=-6:0.1:6;%設(shè)置新的被插值點(diǎn)設(shè)置新的被插值點(diǎn)-6-4-2024600.10.20.30.40.50.60.70.80.91xyh=-6:0.5:6;%設(shè)置新的被插值點(diǎn)設(shè)置新的被插值點(diǎn)36用用MATLABMATLAB作插值

24、計(jì)算作插值計(jì)算一維插值函數(shù)：一維插值函數(shù)：yi=interp1(x，y，xi，method)插值方法插值方法被插值點(diǎn)被插值點(diǎn)插值節(jié)點(diǎn)插值節(jié)點(diǎn)xixi處的插處的插值結(jié)果值結(jié)果nearest ：最鄰近插值：最鄰近插值linear ： 線性插值；線性插值；spline ： 三次樣條插值；三次樣條插值；cubic ： 立方插值。立方插值。缺省時(shí)：缺省時(shí)： 分段線性插值。分段線性插值。 注意：所有的插值方法都要求注意：所有的插值方法都要求x x是單調(diào)的，并且是單調(diào)的，并且xi不不能夠超過能夠超過x的范圍。的范圍。3738程序：程序：x=0:.1:1;x=0:.1:1;y=.3 .5 1 1.4 1.6

25、1 .6 .4 .8 1.5 2;y=.3 .5 1 1.4 1.6 1 .6 .4 .8 1.5 2;yi0=interp1(x,y,0.025,linear)yi0=interp1(x,y,0.025,linear)xi=0:.02:1;xi=0:.02:1;yi=interp1(x,y,xi,linear);yi=interp1(x,y,xi,linear);zi=interp1(x,y,xi,spline);zi=interp1(x,y,xi,spline);wi=interp1(x,y,xi,cubic);wi=interp1(x,y,xi,cubic);plot(x,y,o,xi,

26、yi,r+,xi,zi,gplot(x,y,o,xi,yi,r+,xi,zi,g* *,xi,wi,k.-),xi,wi,k.-)legend(legend(原始點(diǎn)原始點(diǎn),線性點(diǎn)線性點(diǎn),三次樣條三次樣條,三次多項(xiàng)式三次多項(xiàng)式)結(jié)果：結(jié)果：yi0=0.3500yi0=0.3500例例1 1：已知數(shù)據(jù)如表所示，求當(dāng)：已知數(shù)據(jù)如表所示，求當(dāng)x xi i=0.25=0.25時(shí)的時(shí)的y yi i的值。的值。x0.1.2.3.4.5.6.7.8.91y.3.511.4 1.6 1.9 .6.4.81.5 23900.10.20.30.40.50.60.70.80.910.20.40.60.811.21.4

27、1.61.82 原 始 點(diǎn)線 性 點(diǎn)三 次 樣 條三 次 多 項(xiàng) 式40要得到給定的幾個點(diǎn)的對應(yīng)函數(shù)值，可用：要得到給定的幾個點(diǎn)的對應(yīng)函數(shù)值，可用：xi=0.2500 xi=0.2500 0.35000.3500 0.45000.4500yi=interp1(x,y,xi,spline)yi=interp1(x,y,xi,spline)結(jié)果：結(jié)果：yi =1.2088 1.5802 1.3454yi =1.2088 1.5802 1.345441 例：在例：在1-121-12的的1111小時(shí)內(nèi)，每隔小時(shí)內(nèi)，每隔1 1小時(shí)測量一次小時(shí)測量一次溫度，測得的溫度依次為：溫度，測得的溫度依次為：5 5

28、，8 8，9 9，1515，2525，2929，3131，3030，2222，2525，2727，2424。試估計(jì)每隔。試估計(jì)每隔1/101/10小時(shí)的小時(shí)的溫度值。溫度值。hours=1:12;temps=5 8 9 15 25 29 31 30 22 25 27 24;h=1:0.1:12;t=interp1(hours,temps,h,spline); % (直接輸出數(shù)據(jù)將是很多的直接輸出數(shù)據(jù)將是很多的)plot(hours,temps,+,h,t,hours,temps,r:) %作圖作圖xlabel(Hour),ylabel(Degrees Celsius)4202468101251

29、01520253035HourDegrees Celsius43hours=1:12;temps=5 8 9 15 25 29 31 30 22 25 27 24;h=1:0.1:12;t1=interp1(hours,temps,h,nearest); % 最臨近插值最臨近插值t2=interp1(hours,temps,h,linear); % 線性插值線性插值t3=interp1(hours,temps,h,spline); % 三次樣條插值三次樣條插值t4=interp1(hours,temps,h,cubic); % 三次立方插值三次立方插值subplot(2,2,1),plot(h

30、ours,temps,+,h,t1,hours,temps,r:),title(最鄰近差值法最鄰近差值法)subplot(2,2,2), plot(hours,temps,+,h,t2,hours,temps,r:),title(線性插值法線性插值法)subplot(2,2,3),plot(hours,temps,+,h,t3,hours,temps,r:)title(三次樣條插值法三次樣條插值法)subplot(2,2,4),plot(hours,temps,+,h,t4,hours,temps,r:),title(三次立方插值法三次立方插值法)051015010203040最 鄰 近 差

31、值 法051015010203040線 性 插 值 法051015010203040三 次 樣 條 插 值 法051015010203040三 次 立 方 插 值 法45xy機(jī)翼下機(jī)翼下輪廓線輪廓線X035791 11 21 31 41 5Y01 . 21 . 72 . 02 . 12 . 01 . 81 . 21 . 01 . 6例例 已知飛機(jī)下輪廓線上數(shù)據(jù)如下，求已知飛機(jī)下輪廓線上數(shù)據(jù)如下，求x每改變每改變0.1時(shí)的時(shí)的y值。值。46x=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;h=0:0.1:15;y

32、1=interp1(x,y,h,spline); % (直接輸出數(shù)據(jù)將是很多的直接輸出數(shù)據(jù)將是很多的)plot(x,y,+,h,y1,r:) %作圖作圖xlabel(x),ylabel(y)05101500.511.522.5xy47x=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;h=0:0.1:15;y1=interp1(x,y,h,spline); % (直接輸出數(shù)據(jù)將是很多的直接輸出數(shù)據(jù)將是很多的)plot(x,y,+,h,y1,x,y,r:) %作圖作圖xlabel(x),ylabel(y)title

33、(飛機(jī)下輪廓線飛機(jī)下輪廓線)05101500.511.522.5xy飛 機(jī) 下 輪 廓 線48x=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;x0=0:0.1:15;y1=interp1(x,y,x0,); y2=interp1(x,y,x0,cubic); y3=interp1(x,y,x0,spline); subplot(3,1,1),plot(x,y,k+,x0,y1,r)grid title()subplot(3,1,2),plot(x,y,k+,x0,y2,r)gridtitle(linear)

34、subplot(3,1,3),plot(x,y,k+,x0,y3,r)title(spline)49x=0 2 4 5 8 12 12.8 17.2 19.9 20;x=0 2 4 5 8 12 12.8 17.2 19.9 20;y=exp(x).y=exp(x).* *sin(x);sin(x);xx=0:0.25:20;xx=0:0.25:20;yy=spline(x,y,xx);yy=spline(x,y,xx);plot(x,y,o,xx,yy)plot(x,y,o,xx,yy)例：對離散分布在例：對離散分布在y=exp(x)sin(x)y=exp(x)sin(x)函數(shù)曲線上的函數(shù)曲

35、線上的數(shù)據(jù)點(diǎn)進(jìn)行樣條插值計(jì)算數(shù)據(jù)點(diǎn)進(jìn)行樣條插值計(jì)算三次樣條數(shù)據(jù)插值三次樣條數(shù)據(jù)插值MATLABMATLAB實(shí)現(xiàn)還可用函數(shù)實(shí)現(xiàn)還可用函數(shù)splinespline格式：格式：yy=spline(x,y,xx)yy=spline(x,y,xx)數(shù)據(jù)插值數(shù)據(jù)插值5002468101214161820-1012345x 1085166,11)(2xxxg例例用分段線性插值法求插值用分段線性插值法求插值,并觀察插值誤差并觀察插值誤差.1.在在-6,6中平均選取中平均選取5個點(diǎn)作插值；個點(diǎn)作插值；4.在在-6,6中平均選取中平均選取41個點(diǎn)作插值；個點(diǎn)作插值；2.在在-6,6中平均選取中平均選取11個點(diǎn)作插

36、值；個點(diǎn)作插值；3.在在-6,6中平均選取中平均選取21個點(diǎn)作插值；個點(diǎn)作插值；52n=input(輸入輸入n=);x1=-6:12/n:6;%在四個問題中在四個問題中n的取值分別為的取值分別為4,10,20,40y1=1./(1+x1.2);x=-6:1:6;y=1./(1+x.2);y2=interp1(x1,y1,x,linear);y3=interp1(x1,y1,x,spline);y4=interp1(x1,y1,x,nearest);y5=interp1(x1,y1,x,cublic);subplot(2,2,1),plot(x,y,k,x,y2,:.)title(分段線性插值分

37、段線性插值),grid on,legend(g(x),y2)subplot(2,2,2),plot(x,y,k,x,y3,:.)title(3次樣條插值次樣條插值),grid on,legend(g(x),y3)subplot(2,2,3),plot(x,y,k,x,y4,:.)title(最近點(diǎn)插值最近點(diǎn)插值),grid on,legend(g(x),y4)subplot(2,2,4),plot(x,y,k,x,y5,:.)title(3次多項(xiàng)式插值次多項(xiàng)式插值),grid on,legend(g(x),y5)n=4時(shí)-10-5051000.51分 段 線 性 插 值 g(x)y2-10-5

38、0510-0.500.513次 樣 條 插 值 g(x)y3-10-5051000.51最 近 點(diǎn) 插 值 g(x)y4-10-5051000.513次 多 項(xiàng) 式 插 值 g(x)y5n=10時(shí)-10-5051000.51分 段 線 性 插 值 -10-5051000.513次 樣 條 插 值 -10-5051000.51最 近 點(diǎn) 插 值 -10-5051000.513次 多 項(xiàng) 式 插 值 g(x)y2g(x)y3g(x)y4g(x)y5n=20時(shí)-10-5051000.51分 段 線 性 插 值 g(x)y2-10-5051000.513次 樣 條 插 值 g(x)y3-10-5051

39、000.51最 近 點(diǎn) 插 值 g(x)y4-10-5051000.513次 多 項(xiàng) 式 插 值 g(x)y5n=40時(shí)-10-5051000.51分 段 線 性 插 值 g(x)y2-10-5051000.513次 樣 條 插 值 g(x)y3-10-5051000.51最 近 點(diǎn) 插 值 g(x)y4-10-5051000.513次 多 項(xiàng) 式 插 值 g(x)y557二維插值二維插值一、二維插值定義一、二維插值定義二、網(wǎng)格節(jié)點(diǎn)插值法二、網(wǎng)格節(jié)點(diǎn)插值法三、用三、用MatlabMatlab解插值問題解插值問題最鄰近插值最鄰近插值分片線性插值分片線性插值雙線性插值雙線性插值網(wǎng)格節(jié)點(diǎn)數(shù)據(jù)的插值網(wǎng)

40、格節(jié)點(diǎn)數(shù)據(jù)的插值散點(diǎn)數(shù)據(jù)的插值散點(diǎn)數(shù)據(jù)的插值58二維插值的定義二維插值的定義 xyO O第一種（網(wǎng)格節(jié)點(diǎn)）：第一種（網(wǎng)格節(jié)點(diǎn)）：59( (二二) ) 二元插值二元插值二元插值與一元插值的基本思想一致，對原始數(shù)據(jù)點(diǎn)二元插值與一元插值的基本思想一致，對原始數(shù)據(jù)點(diǎn)(x,y,z)(x,y,z)構(gòu)造見世面函構(gòu)造見世面函數(shù)求出插值點(diǎn)數(shù)據(jù)數(shù)求出插值點(diǎn)數(shù)據(jù)(xi,yi,zi)(xi,yi,zi)。一、單調(diào)節(jié)點(diǎn)插值函數(shù)，即一、單調(diào)節(jié)點(diǎn)插值函數(shù)，即x,yx,y向量是單調(diào)的。向量是單調(diào)的。調(diào)用格式調(diào)用格式1 1：zi=interp2(x,y,z,xi,yi,linear)zi=interp2(x,y,z,xi,yi

41、,linear)liner liner 是雙線性插值是雙線性插值 ( (缺省缺省) )調(diào)用格式調(diào)用格式2 2：zi=interp2(x,y,z,xi,yi,nearest)zi=interp2(x,y,z,xi,yi,nearest)nearest nearest 是最近鄰域插值是最近鄰域插值調(diào)用格式調(diào)用格式3 3：zi=interp2(x,y,z,xi,yi,spline)zi=interp2(x,y,z,xi,yi,spline)splinespline是三次樣條插值是三次樣條插值說明：這里說明：這里x x和和y y是兩個獨(dú)立的向量，它們必須是單調(diào)的。是兩個獨(dú)立的向量，它們必須是單調(diào)的。z

42、 z是矩陣，是由是矩陣，是由x x和和y y確定的點(diǎn)上的值。確定的點(diǎn)上的值。z z和和x,yx,y之間的關(guān)系是之間的關(guān)系是z(i,:)=f(x,y(i) z(i,:)=f(x,y(i) z(:,j)=f(x(j),y) z(:,j)=f(x(j),y) 即：當(dāng)即：當(dāng)x x變化時(shí)，變化時(shí)，z z的第的第i i行與行與y y的第的第i i個元素相關(guān)，當(dāng)個元素相關(guān)，當(dāng)y y變變化時(shí)化時(shí)z z的第的第j j列與列與x x的第的第j j個元素相關(guān)。如果沒有對個元素相關(guān)。如果沒有對x,yx,y賦值，則默認(rèn)賦值，則默認(rèn)x=1:n, x=1:n, y=1:my=1:m。n n和和m m分別是矩陣分別是矩陣z

43、z的行數(shù)和列數(shù)。的行數(shù)和列數(shù)。60例例2 2：已知某處山區(qū)地形選點(diǎn)測量坐標(biāo)數(shù)據(jù)為：已知某處山區(qū)地形選點(diǎn)測量坐標(biāo)數(shù)據(jù)為：x=0 x=0 0.50.5 1 1 1.51.5 2 2 2.52.5 3 3 3.53.5 4 4 4.54.5 5 5y=0y=0 0.50.5 1 1 1.51.5 2 2 2.52.5 3 3 3.53.5 4 4 4.54.5 5 5 5.55.5 6 6海拔高度數(shù)據(jù)為：海拔高度數(shù)據(jù)為：z=89 90 87 85 92 91 96 93 90 87 82z=89 90 87 85 92 91 96 93 90 87 82 92 96 98 99 95 91 89 8

44、6 84 82 8492 96 98 99 95 91 89 86 84 82 84 96 98 95 92 90 88 85 84 83 81 8596 98 95 92 90 88 85 84 83 81 85 80 81 82 89 95 96 93 92 89 86 8680 81 82 89 95 96 93 92 89 86 86 82 85 87 98 99 96 97 88 85 82 8382 85 87 98 99 96 97 88 85 82 83 82 85 89 94 95 93 92 91 86 84 8882 85 89 94 95 93 92 91 86 84

45、88 88 92 93 94 95 89 87 86 83 81 9288 92 93 94 95 89 87 86 83 81 92 92 96 97 98 96 93 95 84 82 81 8492 96 97 98 96 93 95 84 82 81 84 85 85 81 82 80 80 81 85 90 93 9585 85 81 82 80 80 81 85 90 93 95 84 86 81 98 99 98 97 96 95 84 8784 86 81 98 99 98 97 96 95 84 87 80 81 85 82 83 84 87 90 95 86 8880 81

46、 85 82 83 84 87 90 95 86 88 80 82 81 84 85 86 83 82 81 80 8280 82 81 84 85 86 83 82 81 80 82 87 88 89 98 99 97 96 98 94 92 8787 88 89 98 99 97 96 98 94 92 87繪出其地貌圖。繪出其地貌圖。61% %程序程序x=0:.5:5;x=0:.5:5;%11%11個數(shù)據(jù)個數(shù)據(jù)y=0:.5:6;y=0:.5:6;%13%13個數(shù)據(jù)個數(shù)據(jù)z z1 1=89 90 87 85 92 91 96 93 90 87 82=89 90 87 85 92 91 96

47、 93 90 87 82;z2=z2=92 96 98 99 95 91 89 86 84 82 8492 96 98 99 95 91 89 86 84 82 84;z3=z3=96 98 95 92 90 88 85 84 83 81 8596 98 95 92 90 88 85 84 83 81 85;z4=z4=80 81 82 89 95 96 93 92 89 86 8680 81 82 89 95 96 93 92 89 86 86;z5=z5=82 85 87 98 99 96 97 88 85 82 8382 85 87 98 99 96 97 88 85 82 83;z6=

48、z6=82 85 89 94 95 93 92 91 86 84 8882 85 89 94 95 93 92 91 86 84 88;z7=z7=88 92 93 94 95 89 87 86 83 81 9288 92 93 94 95 89 87 86 83 81 92;z8=z8=92 96 97 98 96 93 95 84 82 81 8492 96 97 98 96 93 95 84 82 81 84;z9=z9=85 85 81 82 80 80 81 85 90 93 9585 85 81 82 80 80 81 85 90 93 95;z10=z10=84 86 81 98

49、 99 98 97 96 95 84 8784 86 81 98 99 98 97 96 95 84 87;z11=z11=80 81 85 82 83 84 87 90 95 86 8880 81 85 82 83 84 87 90 95 86 88;z12=z12=80 82 81 84 85 86 83 82 81 80 8280 82 81 84 85 86 83 82 81 80 82;z13=z13=87 88 89 98 99 97 96 98 94 92 87;87 88 89 98 99 97 96 98 94 92 87;z=z1;z2;z3;z4;z5;z6;z7;z8;

50、z9;z10;z11;z12;z13;z=z1;z2;z3;z4;z5;z6;z7;z8;z9;z10;z11;z12;z13;%143%143個數(shù)據(jù)個數(shù)據(jù), ,必須是矩陣，得用分號連接必須是矩陣，得用分號連接62x y=meshgrid(x,y);x y=meshgrid(x,y); % %生成網(wǎng)格數(shù)據(jù)生成網(wǎng)格數(shù)據(jù)mesh(x,y,z)mesh(x,y,z) % %加密前的地貌圖加密前的地貌圖01234502468085909510063xi=linspace(0,5,50);%xi=linspace(0,5,50);%加密橫坐標(biāo)數(shù)據(jù)到加密橫坐標(biāo)數(shù)據(jù)到5050個個yi=linspace(0,

51、6,80);%yi=linspace(0,6,80);%加密縱坐標(biāo)數(shù)據(jù)到加密縱坐標(biāo)數(shù)據(jù)到8080個個xii,yii=meshgrid(xi,yi);%xii,yii=meshgrid(xi,yi);%生成新網(wǎng)格數(shù)據(jù)生成新網(wǎng)格數(shù)據(jù)zii=interp2(x,y,z,xii,yii,cubic);%zii=interp2(x,y,z,xii,yii,cubic);%立方插值立方插值mesh(xii,yii,zii)%mesh(xii,yii,zii)%加密后的地貌加密后的地貌012345024675808590951006401234502467580859095100hold onhold on

52、 % % 保持圖形保持圖形plot3(x,y,z+0.1,ob)plot3(x,y,z+0.1,ob) % %原始數(shù)據(jù)用原始數(shù)據(jù)用OO繪出繪出65 已知已知 m n個節(jié)點(diǎn)個節(jié)點(diǎn) ),2 , 1;,.,2 , 1(),(njmizyxijji 其中其中jiyx ,互不相同，不妨互不相同，不妨設(shè)設(shè)bxxxam 21dyyycn 21 構(gòu)造一個二元函數(shù)構(gòu)造一個二元函數(shù)),(yxfz 通過全部已知節(jié)點(diǎn)通過全部已知節(jié)點(diǎn),即即再用再用),(yxf計(jì)算插值，即計(jì)算插值，即).,(*yxfz ),1 ,0;,1 ,0(),(njmizyxfijji 66例：用隨機(jī)數(shù)據(jù)生成地貌圖再進(jìn)行插值例：用隨機(jī)數(shù)據(jù)生成地貌

53、圖再進(jìn)行插值x=rand(100,1)x=rand(100,1)* *4-2;4-2;y=rand(100,1)y=rand(100,1)* *4-2;4-2;z=x.z=x.* *exp(-x.2-y.2);exp(-x.2-y.2);ti=-2:.25:2;ti=-2:.25:2;xi,yi=meshgrid(ti,ti);%xi,yi=meshgrid(ti,ti);%加密數(shù)據(jù)加密數(shù)據(jù)zi=griddata(x,y,z,xi,yi);%zi=griddata(x,y,z,xi,yi);%默認(rèn)線性插值默認(rèn)線性插值mesh(xi,yi,zi)mesh(xi,yi,zi)hold onhold

54、 onplot3(x,y,z,o) plot3(x,y,z,o) 67-2-1012-2-1012-0.500.568第二種（散亂節(jié)點(diǎn)）：第二種（散亂節(jié)點(diǎn)）： yx0 069已知已知n個節(jié)點(diǎn)個節(jié)點(diǎn)),.,2 , 1(),(nizyxiii 其中其中),(iiyx互不相同，互不相同， 構(gòu)造一個二元函數(shù)構(gòu)造一個二元函數(shù)),(yxfz 通過全部已知節(jié)點(diǎn)通過全部已知節(jié)點(diǎn),即即),1 ,0(),(nizyxfiii 再用再用),(yxf計(jì)算插值，即計(jì)算插值，即).,(*yxfz 70 注意：最鄰近插值一般不連續(xù)。具有連續(xù)性的最注意：最鄰近插值一般不連續(xù)。具有連續(xù)性的最簡單的插值是分片線性插值。簡單的插值

55、是分片線性插值。最鄰近插值最鄰近插值x y(x1, y1)(x1, y2)(x2, y1)(x2, y2)O O 二維或高維情形的最鄰近插值，與被插值點(diǎn)最二維或高維情形的最鄰近插值，與被插值點(diǎn)最鄰近的節(jié)點(diǎn)的函數(shù)值即為所求。鄰近的節(jié)點(diǎn)的函數(shù)值即為所求。71 將四個插值點(diǎn)（矩形的四個頂點(diǎn)）處的函數(shù)值依次將四個插值點(diǎn)（矩形的四個頂點(diǎn)）處的函數(shù)值依次簡記為：簡記為： 分片線性插值分片線性插值xy (xi, yj)(xi, yj+1)(xi+1, yj)(xi+1, yj+1)O Of (xi, yj)=f1，f (xi+1, yj)=f2，f (xi+1, yj+1)=f3，f (xi, yj+1)=

56、f472插值函數(shù)為：插值函數(shù)為：jii1ij1jy)xx(xxyyy)yy)(ff ()xx)(ff (f)y, x(fj23i121第二片第二片(上三角形區(qū)域上三角形區(qū)域)：(x, y)滿足滿足iii1ij1jy)xx(xxyyy插值函數(shù)為：插值函數(shù)為：)xx)(ff ()yy)(ff (f)y, x(fi43j141注意：注意：(x, y)當(dāng)然應(yīng)該是在插值節(jié)點(diǎn)所形成的矩形區(qū)當(dāng)然應(yīng)該是在插值節(jié)點(diǎn)所形成的矩形區(qū)域內(nèi)。顯然，分片線性插值函數(shù)是連續(xù)的；域內(nèi)。顯然，分片線性插值函數(shù)是連續(xù)的；分兩片的函數(shù)表達(dá)式如下：分兩片的函數(shù)表達(dá)式如下：第一片第一片(下三角形區(qū)域下三角形區(qū)域)： (x, y)滿足滿

57、足73 雙線性插值是一片一片的空間二次曲面構(gòu)成。雙線性插值是一片一片的空間二次曲面構(gòu)成。雙線性插值函數(shù)的形式如下：雙線性插值函數(shù)的形式如下：)dcy)(bax()y, x(f其中有四個待定系數(shù)，利用該函數(shù)在矩形的四個頂點(diǎn)（插值其中有四個待定系數(shù)，利用該函數(shù)在矩形的四個頂點(diǎn)（插值節(jié)點(diǎn)）的函數(shù)值，得到四個代數(shù)方程，正好確定四個系數(shù)。節(jié)點(diǎn)）的函數(shù)值，得到四個代數(shù)方程，正好確定四個系數(shù)。雙線性插值雙線性插值x y(x1, y1)(x1, y2)(x2, y1)(x2, y2)O O74 要求要求x0,y0 x0,y0單調(diào)；單調(diào)；x x，y y可取為矩陣，或可取為矩陣，或x x取取行向量，行向量，y y

58、取為列向量，取為列向量，x,yx,y的值分別不能超出的值分別不能超出x0,y0 x0,y0的范圍。的范圍。z=interp2(x0,y0,z0,x,y,method)被插值點(diǎn)被插值點(diǎn)插值方法插值方法用用MATLAB作網(wǎng)格節(jié)點(diǎn)數(shù)據(jù)的插值作網(wǎng)格節(jié)點(diǎn)數(shù)據(jù)的插值插值插值節(jié)點(diǎn)節(jié)點(diǎn)被插值點(diǎn)被插值點(diǎn)的函數(shù)值的函數(shù)值nearestnearest 最鄰近插值最鄰近插值linearlinear 雙線性插值雙線性插值cubiccubic 雙三次插值雙三次插值缺省時(shí)缺省時(shí), , 雙線性插值雙線性插值75二者均是常用的二維差值方法，兩者的區(qū)別是，二者均是常用的二維差值方法，兩者的區(qū)別是，interp2interp2的插

59、值數(shù)據(jù)必須是矩形域，即已知數(shù)據(jù)的插值數(shù)據(jù)必須是矩形域，即已知數(shù)據(jù)點(diǎn)（點(diǎn)（x,y)x,y)組成規(guī)則的矩陣，或稱之為柵格，可使組成規(guī)則的矩陣，或稱之為柵格，可使用用meshgidmeshgid生成。而生成。而griddatagriddata函數(shù)的已知數(shù)據(jù)點(diǎn)函數(shù)的已知數(shù)據(jù)點(diǎn)（X X，Y Y）不要求規(guī)則排列，特別是對試驗(yàn)中隨機(jī)）不要求規(guī)則排列，特別是對試驗(yàn)中隨機(jī)沒有規(guī)律采取的數(shù)據(jù)進(jìn)行插值具有很好的效果。沒有規(guī)律采取的數(shù)據(jù)進(jìn)行插值具有很好的效果。76例：測得平板表面例：測得平板表面3 3* *5 5網(wǎng)格點(diǎn)處的溫度分別為：網(wǎng)格點(diǎn)處的溫度分別為： 82 81 80 82 84 82 81 80 82 84

60、79 63 61 65 81 79 63 61 65 81 84 84 82 85 86 84 84 82 85 86 試作出平板表面的溫度分布曲面試作出平板表面的溫度分布曲面z=f(x,y)z=f(x,y)的圖形。的圖形。輸入以下命令：輸入以下命令：x=1:5;y=1:3;temps=82 81 80 82 84;79 63 61 65 81;84 84 82 85 86;mesh(x,y,temps)1.先在三維坐標(biāo)畫出原始數(shù)據(jù)，畫出粗糙的溫度分布曲圖先在三維坐標(biāo)畫出原始數(shù)據(jù)，畫出粗糙的溫度分布曲圖.1234511.522.536065707580859078輸入以下命令輸入以下命令:xi