數(shù)學(xué)建模——擬合與插值

1、1 3/5/2022 數(shù)學(xué)建模教程數(shù)學(xué)建模教程擬擬 合與合與 插插 值值 2 3/5/2022 q在大量的應(yīng)用領(lǐng)域中，人們經(jīng)常面臨這樣的問題：在大量的應(yīng)用領(lǐng)域中，人們經(jīng)常面臨這樣的問題：給定一批數(shù)據(jù)點(diǎn)，需確定滿足特定要求的曲線或給定一批數(shù)據(jù)點(diǎn)，需確定滿足特定要求的曲線或曲面。對(duì)這個(gè)問題有兩種方法。曲面。對(duì)這個(gè)問題有兩種方法。q 一種是一種是插值法插值法，數(shù)據(jù)假定是正確的，要求以某種方法描述數(shù)，數(shù)據(jù)假定是正確的，要求以某種方法描述數(shù)據(jù)點(diǎn)之間所發(fā)生的情況。據(jù)點(diǎn)之間所發(fā)生的情況。q 另一種方法是另一種方法是曲線擬合或回歸曲線擬合或回歸。人們設(shè)法找出某條光滑曲線，。人們設(shè)法找出某條光滑曲線，它最佳地?cái)M

2、合數(shù)據(jù)，但不必要經(jīng)過任何數(shù)據(jù)點(diǎn)。它最佳地?cái)M合數(shù)據(jù)，但不必要經(jīng)過任何數(shù)據(jù)點(diǎn)。 q 本專題的主要目的是：了解插值和擬合的基本內(nèi)容；本專題的主要目的是：了解插值和擬合的基本內(nèi)容；掌握用掌握用MatlabMatlab求解插值與擬合問題的基本命令。求解插值與擬合問題的基本命令。 函數(shù)插值與曲線擬合都是要根據(jù)一組數(shù)據(jù)構(gòu)造一個(gè)函數(shù)作為近似，函數(shù)插值與曲線擬合都是要根據(jù)一組數(shù)據(jù)構(gòu)造一個(gè)函數(shù)作為近似，由于近似的要求不同，二者的數(shù)學(xué)方法上是完全不同的。由于近似的要求不同，二者的數(shù)學(xué)方法上是完全不同的。3 3/5/2022 內(nèi)容提綱內(nèi)容提綱1.擬合問題引例及基本理論2.Matlab求解擬合問題3.應(yīng)用實(shí)例4.插值問

3、題引例及基本理論5.Maltab求解插值問題6.應(yīng)用實(shí)例4 3/5/2022 擬擬 合合 問問 題題 引引 例例 1 1溫度溫度t(0C) 20.5 32.7 51.0 73.0 95.7電阻電阻R( ) 765 826 873 942 1032已知熱敏電阻數(shù)據(jù)：已知熱敏電阻數(shù)據(jù)：求求60600C時(shí)的電阻時(shí)的電阻R。2040608010070080090010001100 設(shè)設(shè) R=at+ba,b為待定系數(shù)為待定系數(shù)一、擬合問題一、擬合問題5 3/5/2022 擬擬 合合 問問 題題 引引 例例 2 2 t (h) 0.25 0.5 1 1.5 2 3 4 6 8c ( g/ml) 19.21

4、 18.15 15.36 14.10 12.89 9.32 7.45 5.24 3.01已知一室模型快速靜脈注射下的血藥濃度數(shù)據(jù)已知一室模型快速靜脈注射下的血藥濃度數(shù)據(jù)(t=0注射注射300mg)求血藥濃度隨時(shí)間的變化規(guī)律求血藥濃度隨時(shí)間的變化規(guī)律c(t).作半對(duì)數(shù)坐標(biāo)系作半對(duì)數(shù)坐標(biāo)系(semilogy)下的圖形下的圖形00( ),ktc tc eck為待定系數(shù)024681001011026 3/5/2022 曲曲 線線 擬擬 合合 問問 題題 的的 提提 法法已知一組（二維）數(shù)據(jù)，即平面上已知一組（二維）數(shù)據(jù)，即平面上 n個(gè)點(diǎn)個(gè)點(diǎn)（xi,yi) i=1,n, 尋求一個(gè)函數(shù)（曲線）尋求一個(gè)函數(shù)

5、（曲線）y=f(x), 使使 f(x) 在某種準(zhǔn)則下與所有在某種準(zhǔn)則下與所有數(shù)據(jù)點(diǎn)最為接近，即曲線擬合得最好。數(shù)據(jù)點(diǎn)最為接近，即曲線擬合得最好。 +xyy=f(x)(xi,yi)i i 為點(diǎn)為點(diǎn)（xi,yi) 與與曲線曲線 y=f(x) 的距離的距離7 3/5/2022 線性最小二乘擬合線性最小二乘擬合 f(x)=a1r1(x)+ +amrm(x)中中函數(shù)函數(shù)rr1 1(x), r(x), rm m(x)(x)的選取的選取 1. 1. 通過機(jī)理分析建立數(shù)學(xué)模型來確定通過機(jī)理分析建立數(shù)學(xué)模型來確定 f(x)f(x)；+f=a1+a2xf=a1+a2x+a3x2f=a1+a2x+a3x2f=a1+

6、a2/xf=aebxf=ae-bx2. 2. 將數(shù)據(jù)將數(shù)據(jù) (xi,yi) i=1, n 作圖，通過直觀判斷確定作圖，通過直觀判斷確定 f(x)：8 3/5/2022 曲線擬合問題最常用的解法曲線擬合問題最常用的解法線性最小二乘法的基本思路線性最小二乘法的基本思路第一步: :先選定一組函數(shù)先選定一組函數(shù) r1(x), r2(x), rm(x), mn, 令令 f(x)=a1r1(x)+a2r2(x)+ +amrm(x) （1）其中其中 a1,a2, am 為待定系數(shù)。為待定系數(shù)。 第二步: 確定確定a1,a2, am 的準(zhǔn)則（最小二乘準(zhǔn)則）：的準(zhǔn)則（最小二乘準(zhǔn)則）：使使n個(gè)點(diǎn)個(gè)點(diǎn)（xi,yi)

7、 與與曲線曲線 y=f(x) 的距離的距離 i 的平方和最小的平方和最小 。記記 )2()()(),(211211221iiknimkkininiiimyxrayxfaaaJ 問題歸結(jié)為，求問題歸結(jié)為，求 a1,a2, am 使使 J(a1,a2, am) 最小。最小。9 3/5/2022 線性最小二乘法的求解：預(yù)備知識(shí)線性最小二乘法的求解：預(yù)備知識(shí)超定方程組超定方程組：方程個(gè)數(shù)大于未知量個(gè)數(shù)的方程組：方程個(gè)數(shù)大于未知量個(gè)數(shù)的方程組)( 221111212111mnyarararyarararnmnmnnmm即即 Ra=ynmnmnnmyyyaaarrrrrrR112111211,其中其中超定

8、方程一般是不存在解的矛盾方程組。超定方程一般是不存在解的矛盾方程組。 如果有向量如果有向量a使得使得 達(dá)到最小，達(dá)到最小，則稱則稱a為上述為上述超定方程的最小二乘解超定方程的最小二乘解。 212211)(imniimiiyararar10 3/5/2022 線性最小二乘法的求解線性最小二乘法的求解定理：定理：當(dāng)當(dāng)R RT TR R可逆時(shí)，超定方程組（可逆時(shí)，超定方程組（3 3）存在最小二乘解，）存在最小二乘解， 且即為方程組且即為方程組 R RT TRa=RRa=RT Ty -y -正則（正規(guī)）方程組正則（正規(guī)）方程組的解：的解：a=(Ra=(RT TR)R)-1-1R RT Ty y所以，曲

9、線擬合的最小二乘法要解決的問題，實(shí)際上就是求以所以，曲線擬合的最小二乘法要解決的問題，實(shí)際上就是求以下超定方程組的最小二乘解的問題。下超定方程組的最小二乘解的問題。nmnmnmyyyaaaxrxrxrxrR111111,)()()()(其中其中Ra=y （3）11 3/5/2022 用用MATLAB解擬合問題解擬合問題1 1、線性最小二乘擬合、線性最小二乘擬合2 2、非線性最小二乘擬合、非線性最小二乘擬合12 3/5/2022 用用MATLAB作線性最小二乘擬合作線性最小二乘擬合1. 1. 作多項(xiàng)式作多項(xiàng)式f(x)=a1xm+ +amx+am+1擬合擬合, ,可利用已有命令可利用已有命令:a=

10、polyfit(x,y,m)3.3.對(duì)超定方程組對(duì)超定方程組)(11nmyaRnmmn可得最小二乘意義下的解?？傻米钚《艘饬x下的解。，用，用yRa2.2.多項(xiàng)式在多項(xiàng)式在x x處的值處的值y y的計(jì)算命令：的計(jì)算命令：y=polyvaly=polyval（a a，x x）輸出擬合多項(xiàng)式系數(shù)輸出擬合多項(xiàng)式系數(shù)a=a1,am,am+1 （數(shù)組）數(shù)組）輸入同長度輸入同長度數(shù)組數(shù)組X，Y擬合多項(xiàng)式擬合多項(xiàng)式次數(shù)次數(shù)13 3/5/2022 即要求即要求 出二次多項(xiàng)式出二次多項(xiàng)式:3221)(axaxaxf中中 的的123( ,)Aa a a 使得使得:1121() iiif xy最小例例 對(duì)下面一組數(shù)

11、據(jù)作二次多項(xiàng)式擬合對(duì)下面一組數(shù)據(jù)作二次多項(xiàng)式擬合xi 0.1 0.2 0.4 0.5 0.6 0.7 0.8 0.9 1 yi 1.978 3.28 6.16 7.34 7.66 9.58 9.48 9.30 11.2 211112222232111111,1,1,1xxyayxxaayxx14 3/5/2022 1）輸入命令）輸入命令：x=0:0.1:1; y=-0.447,1.978,3.28,6.16,7.08,7.34,7.66,9.56,9.48,9.30,11.2; R=(x.2), x, ones(11,1)； A=Ry11 11211121xxxxR此時(shí)MATLAB(zxec1

12、)解法解法1 1解超定方程的方法解超定方程的方法2）計(jì)算結(jié)果）計(jì)算結(jié)果： = -9.8108, 20.1293, -0.03170317.01293.208108.9)(2xxxfRAy15 3/5/2022 16 3/5/2022 2）計(jì)算結(jié)果：）計(jì)算結(jié)果： = -9.8108, 20.1293, -0.0317解法解法2用多項(xiàng)式擬合的命令用多項(xiàng)式擬合的命令MATLAB(zxec2)00.20.40.60.81-20246810120317.01293.208108.9)(2xxxf1）輸入命令：）輸入命令：x=0:0.1:1; y=-0.447,1.978,3.28,6.16,7.08,7

13、.34,7.66,9.56,9.48,9.30,11.2;A=polyfit(x,y,2)z=polyval(A,x);plot(x,y,k+,x,z,r) %作出數(shù)據(jù)點(diǎn)和擬合曲線的圖形作出數(shù)據(jù)點(diǎn)和擬合曲線的圖形17 3/5/2022 1. 1. lsqcurvefitlsqcurvefit已知數(shù)據(jù)點(diǎn)數(shù)據(jù)點(diǎn)：xdata=xdata=（xdata1，xdata2，xdataxdatan n） ydata=ydata=（ydataydata1 1，ydataydata2 2，ydataydatan n）用用MATLAB作非線性最小二乘擬合作非線性最小二乘擬合兩個(gè)求非線性最小二乘擬合的函數(shù)：兩個(gè)求非

14、線性最小二乘擬合的函數(shù)：lsqcurvefitlsqcurvefit、lsqnonlinlsqnonlin。相同點(diǎn)和不同點(diǎn)：兩個(gè)命令都要先建立相同點(diǎn)和不同點(diǎn)：兩個(gè)命令都要先建立M-M-文件文件fun.mfun.m，定義函，定義函數(shù)數(shù)f(x)f(x)，但定義，但定義f(x)f(x)的方式不同的方式不同。211( ,) 2niiiF x xdataydata最小 lsqcurvefitlsqcurvefit用以求含參量用以求含參量x x（向量）的向量值函數(shù)（向量）的向量值函數(shù)F(x,xdata)=F(x,xdata)=（F F（x x，xdataxdata1 1），），F(xiàn) F（x x，xdatax

15、datan n）T T使得使得 18 3/5/2022 輸入格式輸入格式: : (1) x = lsqcurvefit (fun,x0,xdata,ydata); (2) x =lsqcurvefit (fun,x0,xdata,ydata,lb, ub); (3) x =lsqcurvefit (fun,x0,xdata,ydata, lb, ub, options); (4) x, resnorm = lsqcurvefit (fun,x0,xdata,ydata,); (5) x, resnorm, residual = lsqcurvefit (fun,x0,xdata,ydata,)

16、; fun是一個(gè)事先建立的是一個(gè)事先建立的定義函數(shù)定義函數(shù)F(x,xdata) 的的M-文件文件, 自變量為自變量為x和和xdata說明：說明：x = lsqcurvefit (fun,x0,xdata,ydata,options);迭代初值迭代初值已知數(shù)據(jù)點(diǎn)已知數(shù)據(jù)點(diǎn)選項(xiàng)見無選項(xiàng)見無約束優(yōu)化約束優(yōu)化19 3/5/2022 lsqnonlin用以求含參量用以求含參量x x（向量）的向量值函數(shù)（向量）的向量值函數(shù) f(x)f(x)=(f=(f1 1(x),f(x),f2 2(x),f(x),fn n(x)(x)T T ，使得，使得 最小。最小。 其中其中 fi（x）=f（x，xdatai，yda

17、tai） =F(x,xdatai)-ydatai22212( ) ( )( )( )( )Tnfx f xf xfxfx2. lsqnonlin已知數(shù)據(jù)點(diǎn)：已知數(shù)據(jù)點(diǎn)： xdata=xdata=（xdata1，xdata2，xdataxdatan n） ydata=ydata=（ydataydata1 1，ydataydata2 2，ydataydatan n）20 3/5/2022 輸入格式：輸入格式： 1） x=lsqnonlin（fun，x0）； 2） x= lsqnonlin （fun，x0，lb，ub）； 3） x= lsqnonlin （fun，x0， ，lb，ub，options

18、）； 4） x， resnorm= lsqnonlin （fun，x0，）； 5） x， resnorm ， residual= lsqnonlin （fun，x0，）；說明：說明：x= lsqnonlinlsqnonlin （fun，x0，options）；）；fun是一個(gè)事先建立的是一個(gè)事先建立的定義函數(shù)定義函數(shù)f(x)的的M-文件，文件，自變量為自變量為x迭代初值迭代初值選項(xiàng)見無選項(xiàng)見無約束優(yōu)化約束優(yōu)化21 3/5/2022 100200 30040050060070080090010004.54 4.99 5.35 5.65 5.90 6.10 6.26 6.39 6.50 6.59j

19、t310jc100.022111min( , , )22jktjjF a b kabec例例2 用下面一組數(shù)據(jù)擬合用下面一組數(shù)據(jù)擬合 中的參數(shù)中的參數(shù)a，b，k0.0.2( )ktc tabe 該問題即解的最優(yōu)化問題：該問題即解的最優(yōu)化問題：22 3/5/2022 1 1）編寫）編寫M-M-文件文件 curvefun1.mcurvefun1.m function f=curvefun1(x,tdata) f=x(1)+x(2)*exp(-0.02*x(3)*tdata) %其中其中 x(1)=a; x(2)=b；x(3)=k;2）輸入命令）輸入命令tdata=100:100:1000tdata

20、=100:100:1000cdata=cdata=1e-03* *4.54,4.99,5.35,5.65,5.90,6.10,6.26,6.39,4.54,4.99,5.35,5.65,5.90,6.10,6.26,6.39,6.50,6.59;6.50,6.59; x0=0.2,0.05,0.05; x0=0.2,0.05,0.05; x=lsqcurvefit (curvefun1,x0,tdata,cdata)x=lsqcurvefit (curvefun1,x0,tdata,cdata) f= f= curvefun1(x,tdata) F(x，tdata)= ，x=(a，b，k)Tk

21、tktbeabea),(10102. 002. 0解法解法1 1. 用命令用命令lsqcurvefitlsqcurvefit100200 30040050060070080090010004.54 4.99 5.35 5.65 5.90 6.10 6.26 6.39 6.50 6.59jt310jc23 3/5/2022 3 3）運(yùn)算結(jié)果）運(yùn)算結(jié)果：f =0.0043 0.0051 0.0056 0.0059 0.0061 f =0.0043 0.0051 0.0056 0.0059 0.0061 0.0062 0.0062 0.0063 0.0063 0.0063 0.0062 0.0062

22、 0.0063 0.0063 0.0063 x =0.0063 -0.0034 0.2542 x =0.0063 -0.0034 0.25424）結(jié)論）結(jié)論：a=0.0063, b=-0.0034, k=0.2542 0.02 0.25420.0050840.00630.00340.00630.0034ttc tee24 3/5/2022 1）編寫編寫M-M-文件文件 curvefun2.mcurvefun2.m function f=curvefun2(x) tdata=100:100:1000; cdata=1e-03*4.54,4.99,5.35,5.65,5.90, 6.10,6.26

23、,6.39,6.50,6.59; f=x(1)+x(2)*exp(-0.02*x(3)*tdata)- cdata2）輸入命令）輸入命令: x0=0.2,0.05,0.05;x=lsqnonlin(curvefun2,x0)f= curvefun2(x)函數(shù)函數(shù)curvefun2的自變量是的自變量是x，cdata和和tdata是已知參數(shù)，故應(yīng)是已知參數(shù)，故應(yīng)將將cdata tdata的值寫在的值寫在curvefun2.m中中解法解法 2 2 用命令用命令lsqnonlinlsqnonlin 1010.020.02110,(,)ktktTf xF x tdata ctadaabecabec x=

24、 x=（a a，b b，k k）25 3/5/2022 3 3）運(yùn)算結(jié)果為）運(yùn)算結(jié)果為 f =1.0e-003 f =1.0e-003 * *（0.2322 -0.1243 -0.2495 -0.2413 0.2322 -0.1243 -0.2495 -0.2413 -0.1668 -0.0724 0.0241 0.1159 0.2030 0.2792-0.1668 -0.0724 0.0241 0.1159 0.2030 0.2792） x =0.0063 -0.0034 0.2542x =0.0063 -0.0034 0.2542可以看出，兩個(gè)命令的計(jì)算結(jié)果是相同的可以看出，兩個(gè)命令的計(jì)算

25、結(jié)果是相同的。4）結(jié)論）結(jié)論：即擬合得即擬合得a=0.0063 b=-0.0034 k=0.25420.0063 b=-0.0034 k=0.254226 3/5/2022 插值問題插值問題27 3/5/2022 擬合與插值的關(guān)系擬合與插值的關(guān)系說明：說明：函數(shù)插值與曲線擬合都是要根據(jù)一組數(shù)據(jù)構(gòu)造一個(gè)函數(shù)作為近函數(shù)插值與曲線擬合都是要根據(jù)一組數(shù)據(jù)構(gòu)造一個(gè)函數(shù)作為近似，由于近似的要求不同，二者的數(shù)學(xué)方法上完全不同。似，由于近似的要求不同，二者的數(shù)學(xué)方法上完全不同。 實(shí)例：實(shí)例：下面數(shù)據(jù)是某次實(shí)驗(yàn)所得，希望得到x和 f之間的關(guān)系？x124791 21 31 51 7f1 .53 .96 .611

26、.71 5 .61 8 .81 9 .62 0 .62 1 .1MATLAB(cn)問題：問題：給定一批數(shù)據(jù)點(diǎn)，需確定滿足特定要求的曲線或曲面解決方案：解決方案：若不要求曲線（面）通過所有數(shù)據(jù)點(diǎn)，而是要求它反映對(duì)象整體的變化趨勢，就是數(shù)據(jù)擬合數(shù)據(jù)擬合，又稱曲線擬合或曲面擬合。若要求所求曲線（面）通過所給所有數(shù)據(jù)點(diǎn)，就是插值問題插值問題；28 3/5/2022 最臨近插值、線性插值、樣條插值與曲線擬合結(jié)果：最臨近插值、線性插值、樣條插值與曲線擬合結(jié)果：0246810121416180510152025已知數(shù)據(jù)點(diǎn)spline三次多項(xiàng)式插值0246810121416180510152025已知數(shù)據(jù)點(diǎn)

27、linest三次多項(xiàng)式插值0246810121416180510152025已知數(shù)據(jù)點(diǎn)nearest三次多項(xiàng)式插值29 3/5/2022 一一 維維 插插 值值一、插值的定義一、插值的定義二、插值的方法二、插值的方法三、用三、用Matlab解插值問題解插值問題30 3/5/2022 返回返回二維插值二維插值一、二維插值定義一、二維插值定義二、網(wǎng)格節(jié)點(diǎn)插值法二、網(wǎng)格節(jié)點(diǎn)插值法三、用三、用MatlabMatlab解插值問題解插值問題最鄰近插值最鄰近插值分片線性插值分片線性插值雙線性插值雙線性插值網(wǎng)格節(jié)點(diǎn)數(shù)據(jù)的插值網(wǎng)格節(jié)點(diǎn)數(shù)據(jù)的插值散點(diǎn)數(shù)據(jù)的插值散點(diǎn)數(shù)據(jù)的插值31 3/5/2022 一維插值的定義一

28、維插值的定義已知已知 n+1個(gè)節(jié)點(diǎn)個(gè)節(jié)點(diǎn), 1 , 0(),(njyxjj其中其中jx互不相同，不妨設(shè)互不相同，不妨設(shè)),10bxxxan求任一插值點(diǎn)求任一插值點(diǎn))(*jxx 處的插值處的插值.*y0 x1xnx0y1y*x*y32 3/5/2022 構(gòu)造一個(gè)構(gòu)造一個(gè)(相對(duì)簡單的相對(duì)簡單的)函數(shù)函數(shù)),(xfy 通過全部節(jié)點(diǎn)通過全部節(jié)點(diǎn), 即即), 1 ,0()(njyxfjj再用再用)(xf計(jì)算插值，即計(jì)算插值，即).(*xfy 0 x1xnx0y1y*x*y返回返回33 3/5/2022 的函數(shù)值的函數(shù)值 nyyy,10已知已知 y=f(x) 在在n+1 個(gè)點(diǎn)個(gè)點(diǎn)01na xxxbnkyx

29、fxpkkkn, 1 , 0,)()( 構(gòu)造構(gòu)造n次多項(xiàng)式次多項(xiàng)式 pn(x) ,使得使得從而得到從而得到 f(x) 的近似計(jì)算式的近似計(jì)算式 ,),()(baxxpxfn 34 3/5/2022 求解求解 L1(x)=a1 x+a0ixiy0 x1x0y1y已知已知使得使得 L1(xi) = yi . (i=0,1)(1xLy )(xfy 0 x1xOxy如果令如果令1010)(xxxxxl 0101)(xxxxxl 0)(, 1)(1000 xlxl則稱則稱 l0(x) , l1(x)為為x0, x1上的上的線性插值線性插值基函數(shù)基函數(shù)。這時(shí)。這時(shí)10100101yxxxxyxxxx 根據(jù)

30、點(diǎn)斜式得到根據(jù)點(diǎn)斜式得到)()(0010101xxxxyyyxL 并稱其為一次并稱其為一次Lagrange插值多項(xiàng)式。插值多項(xiàng)式。1011()0, ( )1l xl xf(x) L1(x)=y0l0(x)+y1 l1(x)35 3/5/2022 求解求解 L2(x)=a2x2+a1 x+a0使得使得 L2(xi)=yi , i=0,1,2.關(guān)于關(guān)于二次多項(xiàng)式二次多項(xiàng)式的構(gòu)造采用如下方法：令的構(gòu)造采用如下方法：令已知已知ixiy0 x1x0y1y2x2y并由插值條件并由插值條件得到得到)(20100 xxxxyA )(21011xxxxyB )(12022xxxxyC L2(x)=A(x-x1)

31、(x-x2)+B(x-x0)(x-x2)+C(x-x0)(x-x1)L2(x0)=y0 , L2(x1)=y1 , L2(x)=y236 3/5/2022 于是得到于是得到2120210121012002010212)()()()()()()(yxxxxxxxxyxxxxxxxxyxxxxxxxxxL 這時(shí)：這時(shí)：f(x) L2(x)=y0l0(x)+y1 l1(x)+ y2 l2(x)如果令如果令則有則有 jijixlijij, 0, 1)( )()()(2010210 xxxxxxxxxl )()()(2101201xxxxxxxxxl )()()(1202102xxxxxxxxxl 并稱

32、其為并稱其為二次二次Lagrange插值多項(xiàng)式。插值多項(xiàng)式。 基函數(shù)表示基函數(shù)表示 x0, x1, x2上的二次插值基函數(shù)37 3/5/2022 稱為拉格朗日插值基函數(shù)拉格朗日插值基函數(shù)。n0iiiny)x(L)x(P 已知函數(shù)f(x)在n+1個(gè)點(diǎn)x0,x1,xn處的函數(shù)值為 y0,y1,yn 。求一n次多項(xiàng)式函數(shù)Pn(x)，使其滿足： Pn(xi)=yi,i=0,1,n. 解決此問題的拉格朗日插值多項(xiàng)式公式如下其中Li(x) 為n次多項(xiàng)式：)xx()xx)(xx()xx)(xx()xx()xx)(xx()xx)(xx()x(Lni1ii1ii1i0in1i1i10i拉格朗日拉格朗日(Lagr

33、ange)插值插值38 3/5/2022 拉格朗日拉格朗日(Lagrange)插值插值特別地特別地:兩點(diǎn)一次兩點(diǎn)一次(線性線性)插值多項(xiàng)式插值多項(xiàng)式: 101001011yxxxxyxxxxxL三點(diǎn)二次三點(diǎn)二次(拋物拋物)插值多項(xiàng)式插值多項(xiàng)式: 2120210121012002010212yxxxxxxxxyxxxxxxxxyxxxxxxxxxL .,滿足插值條件直接驗(yàn)證可知xLn39 3/5/2022 拉格朗日多項(xiàng)式插值的這種振蕩現(xiàn)象叫 Runge現(xiàn)象現(xiàn)象55,11)(2xxxg 采用拉格朗日多項(xiàng)式插值：選取不同插值節(jié)點(diǎn)個(gè)數(shù)n+1，其中n為插值多項(xiàng)式的次數(shù)，當(dāng)n分別取2,4,6,8,10時(shí)，

34、繪出插值結(jié)果圖形.例例40 3/5/2022 41 3/5/2022 分段線性插值分段線性插值其它,0,)()()(1111110jjjjjjjjjjjnjjjnxxxxxxxxxxxxxxxlxlyxL計(jì)算量與n無關(guān);n越大，誤差越小.nnnxxxxgxL0),()(limxjxj-1xj+1x0 xnxoy42 3/5/2022 66,11)(2xxxg例例用分段線性插值法求插值用分段線性插值法求插值,并觀察插值誤差并觀察插值誤差.在在-6,6中平均選取中平均選取41個(gè)點(diǎn)作插值個(gè)點(diǎn)作插值,結(jié)果如圖示結(jié)果如圖示43 3/5/2022 比分段線性插值更光滑。比分段線性插值更光滑。xyxi-1

35、xiab 在數(shù)學(xué)上，光滑程度的定量描述是：函數(shù)(曲線)的k階導(dǎo)數(shù)存在且連續(xù)，則稱該曲線具有k階光滑性。 光滑性的階次越高，則越光滑。是否存在較低次的分段多項(xiàng)式達(dá)到較高階光滑性的方法？三次樣條插值就是一個(gè)很好的例子。三次樣條插值三次樣條插值44 3/5/2022 , 1,),()(1nixxxxsxSiii,)()3), 1 ,0()()2), 1()()10223niiiiiiixxCxSniyxSnidxcxbxaxs) 1, 1()()(),()(),()(111 nixsxsxsxsxsxsiiiiiiiiiiii自然邊界條件）(0)()()40 nxSxS)(,)4)3)2xSdcba

36、iiiig g( (x x) )為被插值函數(shù)為被插值函數(shù)。lim( )( )nS xg x 三次樣條插值三次樣條插值45 3/5/2022 例例66,11)(2xxxg用三次樣條插值選取用三次樣條插值選取11個(gè)基點(diǎn)計(jì)算插值個(gè)基點(diǎn)計(jì)算插值46 3/5/2022 用用MATLABMATLAB作插值計(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

37、 ： 立方插值。立方插值。缺省時(shí)：缺省時(shí)： 分段線性插值。分段線性插值。 注意：所有的插值方法都要求注意：所有的插值方法都要求x x是單調(diào)的，并且是單調(diào)的，并且xi不不能夠超過能夠超過x的范圍。的范圍。47 3/5/2022 例：在例：在1-121-12的的1111小時(shí)內(nèi)，每隔小時(shí)內(nèi)，每隔1 1小時(shí)測量一次溫小時(shí)測量一次溫度，測得的溫度依次為：度，測得的溫度依次為：5 5，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

38、25 29 31 30 22 25 27 24;h=1:0.1:12;t=interp1(hours,temps,h,spline); (直接輸出數(shù)據(jù)將是很多的)plot(hours,temps,+,h,t,hours,temps,r:) %作圖xlabel(Hour),ylabel(Degrees Celsius)48 3/5/2022 49 3/5/2022 xy機(jī)翼下輪廓線X035791 11 21 31 41 5Y01 . 21 . 72 . 02 . 12 . 01 . 81 . 21 . 01 . 6例例 已知飛機(jī)下輪廓線上數(shù)據(jù)如下，求已知飛機(jī)下輪廓線上數(shù)據(jù)如下，求x每改變每改變0

39、.1時(shí)的時(shí)的y值。值。50 3/5/2022 51 3/5/2022 二維插值的定義二維插值的定義 xyO O第一種（網(wǎng)格節(jié)點(diǎn)）：第一種（網(wǎng)格節(jié)點(diǎn)）：52 3/5/2022 已知已知 m n個(gè)節(jié)點(diǎn)個(gè)節(jié)點(diǎn) ),2 , 1;,.,2 , 1(),(njmizyxijji 其中其中jiyx ,互不相同，不妨設(shè)互不相同，不妨設(shè)bxxxam 21dyyycn 21 構(gòu)造一個(gè)二元函數(shù)構(gòu)造一個(gè)二元函數(shù)),(yxfz 通過全部已知節(jié)點(diǎn)通過全部已知節(jié)點(diǎn),即即再用再用),(yxf計(jì)算插值，即計(jì)算插值，即).,(*yxfz ),1 ,0;,1 ,0(),(njmizyxfijji 53 3/5/2022 第二種（散

40、亂節(jié)點(diǎn)）：第二種（散亂節(jié)點(diǎn)）： yx0 054 3/5/2022 已知已知n個(gè)節(jié)點(diǎn)個(gè)節(jié)點(diǎn)),.,2 , 1(),(nizyxiii 其中其中),(iiyx互不相同，互不相同， 構(gòu)造一個(gè)二元函數(shù)構(gòu)造一個(gè)二元函數(shù)),(yxfz 通過全部已知節(jié)點(diǎn)通過全部已知節(jié)點(diǎn),即即),1 ,0(),(nizyxfiii 再用再用),(yxf計(jì)算插值，即計(jì)算插值，即).,(*yxfz 55 3/5/2022 注意：注意：最鄰近插值一般不連續(xù)。具有連續(xù)性的最簡單的插值是分片線性插值。最鄰近插值最鄰近插值x y(x1, y1)(x1, y2)(x2, y1)(x2, y2)O O 二維或高維情形的最鄰近插值，與被插值點(diǎn)

41、最鄰近的節(jié)點(diǎn)的函數(shù)值即為所求。56 3/5/2022 將四個(gè)插值點(diǎn)（矩形的四個(gè)頂點(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)=f457 3/5/2022 插值函數(shù)為：jii1ij1jy)xx(xxyyy)yy)(ff ()xx)(ff (f)y, x(fj23i121第二片(上三角形區(qū)域)：(x, y)滿足iii1ij1jy)xx(xxyyy插值函數(shù)為：)xx)(ff ()yy

42、)(ff (f)y, x(fi43j141注意注意：(x, y)當(dāng)然應(yīng)該是在插值節(jié)點(diǎn)所形成的矩形區(qū)域內(nèi)。顯然，分片線性插值函數(shù)是連續(xù)的；分兩片的函數(shù)表達(dá)式如下：第一片(下三角形區(qū)域)： (x, y)滿足58 3/5/2022 雙線性插值是一片一片的空間二次曲面構(gòu)成。雙線性插值函數(shù)的形式如下：)dcy)(bax()y, x(f其中有四個(gè)待定系數(shù)，利用該函數(shù)在矩形的四個(gè)頂點(diǎn)（插值節(jié)點(diǎn)）的函數(shù)值，得到四個(gè)代數(shù)方程，正好確定四個(gè)系數(shù)。雙線性插值雙線性插值x y(x1, y1)(x1, y2)(x2, y1)(x2, y2)O O59 3/5/2022 要求要求x0,y0 x0,y0單調(diào)；單調(diào)；x x，

43、y y可取可取為矩陣，或?yàn)榫仃?，或x x取取行向量，行向量，y y取為列向量，取為列向量，x,yx,y的值分別不能超出的值分別不能超出x0,y0 x0,y0的范圍。的范圍。z=interp2(x0,y0,z0,x,y,method)被插值點(diǎn)插值方法用用MATLAB作網(wǎng)格節(jié)點(diǎn)數(shù)據(jù)的插值作網(wǎng)格節(jié)點(diǎn)數(shù)據(jù)的插值插值節(jié)點(diǎn)被插值點(diǎn)的函數(shù)值nearest nearest 最鄰近插值最鄰近插值linear linear 雙線性插值雙線性插值cubic cubic 雙三次插值雙三次插值缺省時(shí)缺省時(shí), , 雙線性插值雙線性插值60 3/5/2022 例：測得平板表面例：測得平板表面3 3* *5 5網(wǎng)格點(diǎn)處的溫度

44、分別為：網(wǎng)格點(diǎn)處的溫度分別為： 82 81 80 82 84 82 81 80 82 84 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ù)，畫出粗糙的溫度分布曲圖.61 3/5/2022 62 3/5/2022 再輸入以下命

45、令:xi=1:0.2:5;yi=1:0.2:3;zi=interp2(x,y,temps,xi,yi,cubic);mesh(xi,yi,zi)畫出插值后的溫度分布曲面圖. 2以平滑數(shù)據(jù),在x、y方向上每隔0.2個(gè)單位的地方進(jìn)行插值.63 3/5/2022 64 3/5/2022 例例 山區(qū)地貌：山區(qū)地貌： 在某山區(qū)測得一些地點(diǎn)的高程如下表。平面區(qū)域?yàn)樵谀成絽^(qū)測得一些地點(diǎn)的高程如下表。平面區(qū)域?yàn)?1200=x=4000,1200=y0)k(0)模型假設(shè)模型假設(shè)1. 1. 機(jī)體看作一個(gè)房室，室內(nèi)血藥濃度均勻機(jī)體看作一個(gè)房室，室內(nèi)血藥濃度均勻一室模型一室模型模型建立模型建立d/c(0) 3得：由假設(shè)-kcdtdc 2得：由假設(shè)ktevdtc)( 在此，在此，d=300mg，t