數(shù)學(xué)建模中的數(shù)值方法

1、微分方程模型微分方程模型建立與求解建立與求解如物質(zhì)流動(dòng)與擴(kuò)散問(wèn)題膜和弦的振動(dòng)問(wèn)題力學(xué)和運(yùn)動(dòng)問(wèn)題溫度分布病毒傳染種群變化1 1、用微分方程建模的大致步驟、思、用微分方程建模的大致步驟、思 路及注意點(diǎn)路及注意點(diǎn)2 2、微分方程求解方法概述、微分方程求解方法概述3 3、簡(jiǎn)單微分方程的、簡(jiǎn)單微分方程的matlabmatlab解法解法4 4、復(fù)雜方程的理論求解方法、復(fù)雜方程的理論求解方法5 5、模型參數(shù)辨識(shí)及函數(shù)擬合問(wèn)題、模型參數(shù)辨識(shí)及函數(shù)擬合問(wèn)題微分方程模型是建模中常見(jiàn)的建模方法大致建模思路是：大致建模思路是：根據(jù)對(duì)現(xiàn)實(shí)對(duì)象特性的認(rèn)識(shí)，分析其因果關(guān)系，找出反映內(nèi)部機(jī)理的規(guī)律，建立模型的基本結(jié)構(gòu)。然后根

2、據(jù)該結(jié)構(gòu)，通過(guò)給定數(shù)據(jù)用參數(shù)辨識(shí)的方法來(lái)確定模型的參數(shù)，從而最終確定模型。步驟：（1 1）建立基本結(jié)構(gòu)）建立基本結(jié)構(gòu)，如是偏微分、還是常微分是一階還是二階；（2 2）參數(shù)辨識(shí)）參數(shù)辨識(shí)，如最小二乘法建立模型結(jié)構(gòu)的大致思路：建立模型結(jié)構(gòu)的大致思路：結(jié)合問(wèn)題用微積分思想的實(shí)現(xiàn)。這個(gè)過(guò)程在建結(jié)合問(wèn)題用微積分思想的實(shí)現(xiàn)。這個(gè)過(guò)程在建模中很多時(shí)候不需要自己完全的創(chuàng)新，更多是模中很多時(shí)候不需要自己完全的創(chuàng)新，更多是需要找到相關(guān)的模型做一些改進(jìn)。（因此，需要找到相關(guān)的模型做一些改進(jìn)。（因此，快快查查到相關(guān)研究是重要的，平時(shí)多積累查找文查查到相關(guān)研究是重要的，平時(shí)多積累查找文獻(xiàn)資料的方法獻(xiàn)資料的方法）注意：注

3、意：需要針對(duì)問(wèn)題修改模型。需要針對(duì)問(wèn)題修改模型。這是建模中最這是建模中最關(guān)鍵的一步，注意要有針對(duì)問(wèn)題的結(jié)合分析，關(guān)鍵的一步，注意要有針對(duì)問(wèn)題的結(jié)合分析，為什么用這個(gè)模型及為什么要做修改，為什么用這個(gè)模型及為什么要做修改，對(duì)模型對(duì)模型的借用方式應(yīng)該是隨內(nèi)容的分析展開(kāi)而自然的的借用方式應(yīng)該是隨內(nèi)容的分析展開(kāi)而自然的引入，不應(yīng)該是沒(méi)有分析強(qiáng)行拷貝。引入，不應(yīng)該是沒(méi)有分析強(qiáng)行拷貝。計(jì)算參數(shù)的方法：計(jì)算參數(shù)的方法：最小二乘法最小二乘法情況情況1 1：如果模型中待求解函數(shù)可由所提的數(shù)如果模型中待求解函數(shù)可由所提的數(shù)學(xué)模型通過(guò)學(xué)模型通過(guò)精確精確理論理論推導(dǎo)推導(dǎo)求解（含有未知參求解（含有未知參數(shù)），數(shù)），則可

4、利用此理論解和給定觀測(cè)數(shù)值，再則可利用此理論解和給定觀測(cè)數(shù)值，再借助借助最小二乘法進(jìn)行參數(shù)計(jì)算最小二乘法進(jìn)行參數(shù)計(jì)算；情況情況2 2：如果模型中待求解函數(shù)無(wú)法通過(guò)所提如果模型中待求解函數(shù)無(wú)法通過(guò)所提的數(shù)學(xué)模型精確求解，那就需通過(guò)數(shù)值方法獲的數(shù)學(xué)模型精確求解，那就需通過(guò)數(shù)值方法獲得模型的近似解（含有未知參數(shù)），然后再用得模型的近似解（含有未知參數(shù)），然后再用最小二乘法進(jìn)行參數(shù)計(jì)算。最小二乘法進(jìn)行參數(shù)計(jì)算。 重金屬在土壤中的傳播：重金屬在土壤中的傳播：（1 1）由于是在土壤中擴(kuò)散，由土壤傳播的特性）由于是在土壤中擴(kuò)散，由土壤傳播的特性（慢，相對(duì)于空氣或液體中），因此，這個(gè)題更（慢，相對(duì)于空氣或液體

5、中），因此，這個(gè)題更多的要求我們分析物質(zhì)的空間分布，而不側(cè)重各多的要求我們分析物質(zhì)的空間分布，而不側(cè)重各區(qū)域內(nèi)重金屬物質(zhì)隨機(jī)時(shí)間的變化規(guī)律。同時(shí)，區(qū)域內(nèi)重金屬物質(zhì)隨機(jī)時(shí)間的變化規(guī)律。同時(shí)，主要是數(shù)據(jù)中也沒(méi)有給出我們關(guān)于時(shí)間的數(shù)據(jù)；主要是數(shù)據(jù)中也沒(méi)有給出我們關(guān)于時(shí)間的數(shù)據(jù)；（2 2）物質(zhì)污染擴(kuò)散是源點(diǎn)濃度最大，然后向四）物質(zhì)污染擴(kuò)散是源點(diǎn)濃度最大，然后向四周空間區(qū)域擴(kuò)散，梯次減小。周空間區(qū)域擴(kuò)散，梯次減小。（3 3）假設(shè)物質(zhì)是由高濃度區(qū)向低濃度區(qū)擴(kuò)散）假設(shè)物質(zhì)是由高濃度區(qū)向低濃度區(qū)擴(kuò)散 最小二乘法：最小二乘法：111( ,), ( ,), (,)iiinnnu x y zu x y zu xyz

6、已知數(shù)據(jù)已知數(shù)據(jù)模型求解模型求解111(,), ( ,), (,)iiinnnu x y zu x y zu xyz如果模型正確的話，即模型能夠準(zhǔn)確刻畫(huà)現(xiàn)象如果模型正確的話，即模型能夠準(zhǔn)確刻畫(huà)現(xiàn)象，那也意味著，模型解和觀測(cè)值之間的誤差差，那也意味著，模型解和觀測(cè)值之間的誤差差異應(yīng)該很小，因此，我們選擇參數(shù)的方法為：異應(yīng)該很小，因此，我們選擇參數(shù)的方法為：21min ( ,)( ,)niiiiiiiu x y zu x y z待估計(jì)參數(shù)那現(xiàn)在的問(wèn)題是：那現(xiàn)在的問(wèn)題是：這樣模型的精確解好求嗎？這樣模型的精確解好求嗎？微分方程的解析解微分方程的解析解 求微分方程（組）解析解的命令（matlab）:d

7、solve(方程方程1,方程方程2,方程方程n,初始條件初始條件,自變量自變量) 結(jié) 果：u = tg(t+c1) 解解 輸入命令: y=dsolve(D2y+4*Dy+29*y=0,y(0)=0,Dy(0)=15,x)結(jié) 果 為 : y =3*exp(-2*x)*sin(5*x)解解 輸入命令 ： x,y,z=dsolve(Dx=2*x-3*y+3*z,Dy=4*x-5*y+3*z,Dz=4*x-4*y+2*z, t) 結(jié) 果 為：x = (c1-c2+c3+c2e -3t-c3e-3t)e2t y = -c1e-4t+c2e-4t+c2e-3t-c3e-3t+c1-c2+c3)e2t z

8、= (-c1e-4t+c2e-4t+c1-c2+c3)e2t 對(duì)于偏微分方程而言，一般借助對(duì)于偏微分方程而言，一般借助matlabmatlab偏微分偏微分方程工具箱方程工具箱PDEtoolPDEtool，有可視化操作窗口，也可，有可視化操作窗口，也可以自己編制命令。但是缺點(diǎn)是，只能處理二階以自己編制命令。但是缺點(diǎn)是，只能處理二階問(wèn)題。問(wèn)題。 我們重點(diǎn)要說(shuō)的是，這些命令不是萬(wàn)能我們重點(diǎn)要說(shuō)的是，這些命令不是萬(wàn)能的，如果可以求解當(dāng)然好，否則，可能的，如果可以求解當(dāng)然好，否則，可能需要自己根據(jù)計(jì)算方法自己編寫(xiě)程序。需要自己根據(jù)計(jì)算方法自己編寫(xiě)程序。因此，明白一點(diǎn)原理是必要的。因此，明白一點(diǎn)原理是必要

9、的。()( ), ( )dduduprqufaxbdxdxdxu au b()()( , ),( , )( , )( , ),( , )uuppquf x yx yxxyyu x yx yx y微分方程的差分法，即用微分方程的差分法，即用差分式代替微分式差分式代替微分式 11211,112221,21,12,11,1,;,;,ThNNNNNNuuuuuuuuuu21hHugh也可得到也可得到通過(guò)上面過(guò)程，我們可以求解問(wèn)題的真是解的通過(guò)上面過(guò)程，我們可以求解問(wèn)題的真是解的近似解（含有未知參數(shù)）近似解（含有未知參數(shù)）通過(guò)調(diào)整該函數(shù)中待定系數(shù)，使得該函數(shù)與已通過(guò)調(diào)整該函數(shù)中待定系數(shù)，使得該函數(shù)與已知

10、點(diǎn)集的差別知點(diǎn)集的差別( (最小二乘意義最小二乘意義) )最小最小函數(shù)的插值與擬合函數(shù)的插值與擬合曲曲 線線 擬擬 合合 問(wèn)問(wèn) 題題 的的 提提 法法已知一組（二維）數(shù)據(jù)，即平面上已知一組（二維）數(shù)據(jù)，即平面上 n n個(gè)點(diǎn)個(gè)點(diǎn)（x xi i, ,y yi i) ) i i=1,=1,n n, , 尋求一個(gè)函數(shù)（曲線）尋求一個(gè)函數(shù)（曲線）y y= =f f( (x x),), 使使 f f( (x x) ) 在某在某種準(zhǔn)則下與所有數(shù)據(jù)點(diǎn)最為接近，即曲線擬合得最好種準(zhǔn)則下與所有數(shù)據(jù)點(diǎn)最為接近，即曲線擬合得最好 +x xy yy y= =f f( (x x) )(x xi,y yi)i i i 為點(diǎn)

11、為點(diǎn)（x xi i, ,y yi i) ) 與與曲線曲線 y y= =f f( (x x) ) 的的距離距離第一步: :先選定一組函數(shù)先選定一組函數(shù) r r1 1( (x x), ), r r2 2( (x x), ,), ,r rm m( (x x), ), m m n n, , 令令 f f( (x x)=)=a a1 1r r1 1( (x x)+)+a a2 2r r2 2( (x x)+ +)+ +a am mr rm m( (x x) ) （1 1）其中其中 a a1 1, ,a a2 2, , ,a am m 為待定系數(shù)為待定系數(shù)第二步: : 確定確定a a1 1, ,a a2

12、2, , ,a am m 的準(zhǔn)則（最小二乘準(zhǔn)則）：的準(zhǔn)則（最小二乘準(zhǔn)則）：使使n n個(gè)點(diǎn)個(gè)點(diǎn)（x xi i, ,y yi i) ) 與與曲線曲線 y y= =f f( (x x) ) 的距離的距離 i i 的平方和最小的平方和最小 記記 )2()()(),(211211221iiknimkkininiiimyxrayxfaaaJ 問(wèn)題：求問(wèn)題：求 a a1 1, ,a a2 2, , ,a am m 使使 J J ( (a a1 1, ,a a2 2, , ,a am m) ) 最小最小超定方程組超定方程組：方程個(gè)數(shù)大于未知量個(gè)數(shù)的方程組：方程個(gè)數(shù)大于未知量個(gè)數(shù)的方程組11 1122111 1

13、22 ()mmnnnmmnr ar ar aynmr ar ar ay即即 Ra=y111112112,mnnnmmnayrrrRayrrray其中其中超定方程組一般不存在解的矛盾方程組超定方程組一般不存在解的矛盾方程組超定方程組一般不存在解的矛盾方程組超定方程組一般不存在解的矛盾方程組 如果有向量如果有向量a使得使得 達(dá)到最小，達(dá)到最小，則稱(chēng)則稱(chēng)a為上述為上述超定方程組的最小二乘解超定方程組的最小二乘解 212211)(imniimiiyararar 定理：定理：當(dāng)當(dāng)RTR可逆時(shí)，超定方程組（可逆時(shí)，超定方程組（3 3）存在最小二乘解，）存在最小二乘解，且即為方程組且即為方程組 RTRa=

14、=RTy的解：的解：a=(RTR)-1RTy 所以，曲線擬合的最小二乘法要解決的問(wèn)題，實(shí)際上就是所以，曲線擬合的最小二乘法要解決的問(wèn)題，實(shí)際上就是求以下超定方程組的最小二乘解的問(wèn)題求以下超定方程組的最小二乘解的問(wèn)題111111()(),()()mnmnmnayrxrxRayrxrxay其中其中Ra=y1. 1. 通過(guò)機(jī)理分析建立數(shù)學(xué)模型來(lái)確定通過(guò)機(jī)理分析建立數(shù)學(xué)模型來(lái)確定 f f( (x x) )；+f f= =a a1 1+ +a a2 2x xf f= =a a1 1+ +a a2 2x x+ +a a3 3x x2 2f f= =a a1 1+ +a a2 2x x+ +a a3 3x

15、x2 2f f= =a a1 1+ +a a2 2/ /x xf f= =a ae ebxbxf f= =a ae e- -bxbx 2. 2. 將數(shù)據(jù)將數(shù)據(jù) ( (x xi i, ,y yi i) )作圖，通過(guò)直觀判斷確定作圖，通過(guò)直觀判斷確定 f f( (x x) )1. 1. 作多項(xiàng)式作多項(xiàng)式f f( (x x)=)=a a1 1x xm m+ + +a am mx x+ +a am+1m+1擬合擬合, ,可利用已有程序可利用已有程序: :a=polyfit(x,y,m)2. 2. 對(duì)超定方程組對(duì)超定方程組)(11nmyaRnmmn可得最小二乘意義下的解可得最小二乘意義下的解，用，用yR

16、a輸出擬合多項(xiàng)式系數(shù)輸出擬合多項(xiàng)式系數(shù)a=a1, ,am , am+1 ( (數(shù)組數(shù)組) )）輸入同長(zhǎng)度輸入同長(zhǎng)度的數(shù)組的數(shù)組x x，y y擬合多項(xiàng)擬合多項(xiàng)式次數(shù)式次數(shù)即要求即要求 出二次多項(xiàng)式出二次多項(xiàng)式: :3221)(axaxaxf中中 的的),(321aaaA 使得使得: :最小 )(1112iiiyxf例題例題 對(duì)下面一組數(shù)據(jù)作二次多項(xiàng)式擬合對(duì)下面一組數(shù)據(jù)作二次多項(xiàng)式擬合1 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)

17、；A=Ry211211111 1xxRxx此時(shí)解法解法1 1用解超定方程的方法用解超定方程的方法2 2）計(jì)算結(jié)果）計(jì)算結(jié)果： = -9.8108 20.1293 -0.03170317.01293.208108.9)(2xxxf1 1）輸入以下命令：）輸入以下命令： x=0:0.1: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; y=-0.447 1.978 3.28 6.16 7.08 7.34 7.66 9.56 9.48 9.30 11.2; A=polyfit(x,y,2) A=poly

18、fit(x,y,2) z= z=polyval(A,xpolyval(A,x);%);%計(jì)算點(diǎn)計(jì)算點(diǎn)x x的函數(shù)值的函數(shù)值 plot(x,y,k+,x,z,rplot(x,y,k+,x,z,r) %) %作出數(shù)據(jù)點(diǎn)和擬合曲線的圖形作出數(shù)據(jù)點(diǎn)和擬合曲線的圖形2 2）計(jì)算結(jié)果：）計(jì)算結(jié)果： = -9.8108 20.1293 -0.0317= -9.8108 20.1293 -0.0317解法解法2 2用多項(xiàng)式擬合的命令用多項(xiàng)式擬合的命令00.20.40.60.81-20246810120317.01293.208108.9)(2xxxf1. 1. lsqcurvefitlsqcurvefit已知

19、數(shù)據(jù)點(diǎn)數(shù)據(jù)點(diǎn)： xdataxdata= =（xdataxdata1 1,xdata,xdata2 2,xdataxdatan n）, , ydataydata= =（ydataydata1 1,ydata,ydata2 2,ydataydatan n） 用用MATLABMATLAB作非線性最小二乘擬合作非線性最小二乘擬合 MATLAB MATLAB提供了兩個(gè)求非線性最小二乘擬合的函數(shù)：提供了兩個(gè)求非線性最小二乘擬合的函數(shù)：lsqcurvefitlsqcurvefit和lsqnonlinlsqnonlin兩個(gè)命令都要先建立兩個(gè)命令都要先建立M M文件文件fun.mfun.m，在其中定義函數(shù)，在其

20、中定義函數(shù)f f( (x x) )，但，但兩者定義兩者定義f f( (x x) )的方式是不同的的方式是不同的, ,可參考例題可參考例題.21( ,) niiiF x xdataydata最小lsqcurvefitlsqcurvefit用以求含參量用以求含參量x x（向量）的向量值函數(shù)（向量）的向量值函數(shù) F(x,xdataF(x,xdata)=)=（F F（x,xdatax,xdata1 1）,F,F（x x，xdataxdatan n）T T中的參變量中的參變量x(x(向量向量),),使得使得 輸入格式為輸入格式為: : (1) x = x = lsqcurvefitlsqcurvefit

21、 (fun,x0,xdata,ydata); (fun,x0,xdata,ydata); (2) x =x =lsqcurvefitlsqcurvefit(fun,x0,xdata,ydata,options);(fun,x0,xdata,ydata,options); (3) x=x=lsqcurvefitlsqcurvefit(fun,x0,xdata,ydata,options,grad);(fun,x0,xdata,ydata,options,grad); (4) x,optionsx,options=lsqcurvefitlsqcurvefit(fun,x0,xdata,ydata,

22、);(fun,x0,xdata,ydata,); (5) x,options,funvalx,options,funval=lsqcurvefitlsqcurvefit(fun,x0,xdata,ydata,);(fun,x0,xdata,ydata,);(6) x,options,funval,Jacobx,options,funval,Jacob=lsqcurvefitlsqcurvefit(fun,x0,xdata, (fun,x0,xdata, ydataydata,);,);funfun是一個(gè)事先建立的是一個(gè)事先建立的定義函數(shù)定義函數(shù)F(x,xdata) 的的M M文件文件, , 自

23、變量為自變量為x和和xdataxdata說(shuō)明:x=lsqcurvefit(fun,x0,xdata,ydata,options);x=lsqcurvefit(fun,x0,xdata,ydata,options);迭代初值迭代初值已知數(shù)據(jù)點(diǎn)已知數(shù)據(jù)點(diǎn)選項(xiàng)見(jiàn)無(wú)選項(xiàng)見(jiàn)無(wú)約束優(yōu)化約束優(yōu)化2 2、lsqnonlinlsqnonlin 1) x=x=lsqnonlinlsqnonlin（funfun，x0 x0）； 2)x=x=lsqnonlinlsqnonlin（funfun，x0 x0，optionsoptions）； 3)x= x= lsqnonlinlsqnonlin（funfun，x0 x0，

24、optionsgradoptionsgrad）； 4)xx，options=options=lsqnonlinlsqnonlin （funfun，x0 x0，）； 5)xx，optionsoptions，funvalfunval=lsqnonlinlsqnonlin（funx0funx0，）；說(shuō)明：x= x= lsqnonlinlsqnonlin （funfun，x0 x0，optionsoptions）；）；funfun是一個(gè)事先建立的定義是一個(gè)事先建立的定義函數(shù)函數(shù)f(x)的的M M文件，文件，自變量為自變量為x x迭代初值迭代初值選項(xiàng)見(jiàn)無(wú)選項(xiàng)見(jiàn)無(wú)約束優(yōu)化約束優(yōu)化100200 300400

25、50060070080090010004.54 4.99 5.35 5.65 5.90 6.10 6.26 6.39 6.50 6.59jt3( 10 )jc 100.0221min( , , )ejktjjF a b kabc例題例題 用下面一組數(shù)據(jù)擬合用下面一組數(shù)據(jù)擬合 中的參數(shù)中的參數(shù)a a，b b，k k0.0.2( )ektc tab該問(wèn)題即解最優(yōu)化問(wèn)題：該問(wèn)題即解最優(yōu)化問(wèn)題： 1 1）編寫(xiě)）編寫(xiě)M M文件文件 curvefun1.m function f=curvefun1(x,tdata) f=x(1)+x(2)*exp(-0.02*x(3)*tdata) %其中其中 x(1)=

26、a; x(2)=b；x(3)=k;2）輸入命令）輸入命令tdata=100:100:1000cdata=1e-03*4.54,4.99,5.35,5.65,5.90,6.10,6.26,6.39,6.50,6.59;x0=0.2,0.05,0.05;x=lsqcurvefit (curvefun1,x0,tdata,cdata)f= curvefun1(x,tdata) F(x，tdata)= ，x=(a，b，k)1010.020.02T(e,e)ktktabab解法解法1 1. 用命令用命令lsqcurvefit1）編寫(xiě)編寫(xiě)M M文件文件 curvefun2.m function f=cur

27、vefun2(x) tdata=100:100:1000; cdata=1e-03*4.54,4.99,5.35,5.65,5.90, 6.10,6.26,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)解法解法 2 2 用命令用命令lsqnonlinlsqnonlina=0.0063, b=-0.0034, k=0.2542 根據(jù)給定的數(shù)據(jù)（有限數(shù)據(jù)點(diǎn)），去找符合該數(shù)據(jù)規(guī)律特征的原始函數(shù)的特

28、征（全局特征），除了通過(guò)機(jī)理分析以外，我們還可以通過(guò)插值的方法。一維插值的定義一維插值的定義已知已知 n+1個(gè)節(jié)點(diǎn)個(gè)節(jié)點(diǎn)(,) (0,1, ,jjxyjn 其中其中jx互不相同，不妨設(shè)互不相同，不妨設(shè)01,nxxx求任一插值點(diǎn)求任一插值點(diǎn))(*jxx 處的插值處的插值.*y0 x1xnx0y1y節(jié)點(diǎn)可視為由節(jié)點(diǎn)可視為由)(xgy 產(chǎn)生產(chǎn)生,g表達(dá)式復(fù)雜表達(dá)式復(fù)雜,或無(wú)封閉形式或無(wú)封閉形式,或未知或未知.*x*y 構(gòu)造一個(gè)構(gòu)造一個(gè)(相對(duì)簡(jiǎn)單的相對(duì)簡(jiǎn)單的)函數(shù)函數(shù)),(xfy 通過(guò)全部節(jié)點(diǎn)通過(guò)全部節(jié)點(diǎn), 即即()(0,1,)jjfxyjn再用再用)(xf計(jì)算插值，即計(jì)算插值，即).(*xfy 0

29、x1xnx0y1y*x*y 稱(chēng)為拉格朗日插值基函數(shù)拉格朗日插值基函數(shù)0( )( )nniiiP xL xy 已知函數(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. 解決此問(wèn)題的拉格朗日插值多項(xiàng)式公式如下其中Li(x) 為n次多項(xiàng)式：01110111()()()()()( )()()()()()iiniiiiiiiinxxxxxxxxxxL xxxxxxxxxxx拉格朗日拉格朗日(Lagrange)插值插值0111111()(),(),0,nnjjjjjjjjjjjjjjLxy lxxxxxxxx

30、xxlxxxxxx其 他計(jì)算量與n無(wú)關(guān);n越大，誤差越小.nnnxxxxgxL0),()(limxjxj-1xj+1x0 xnxy比分段線性插值更光滑比分段線性插值更光滑xyxi-1 xiab數(shù)學(xué)上，光滑程度的定量描述是：函數(shù)(曲線)的k階導(dǎo)數(shù)存在且連續(xù)，則稱(chēng)該曲線具有k階光滑性光滑性的階次越高，則越光滑是否存在較低次的分段多項(xiàng)式達(dá)到較高階光滑性的方法？三次樣條插值就是一個(gè)很好的例子三次樣條插值三次樣條插值三次樣條插值三次樣條插值最鄰近插值最鄰近插值立方插值立方插值用用MATLAB作插值計(jì)算作插值計(jì)算一維插值函數(shù)：一維插值函數(shù)：yi=interp1(x，y，xi，method)插值方法插值方法

31、被插值點(diǎn)被插值點(diǎn)插值節(jié)點(diǎn)插值節(jié)點(diǎn)xi處的處的插值結(jié)果插值結(jié)果nearest 最鄰近插值；最鄰近插值；linear 線性插值；線性插值；spline 三次樣條插值；三次樣條插值；cubic 立方插值；立方插值； 缺省時(shí)缺省時(shí) 分段線性插值分段線性插值注意：所有的插值方法都要注意：所有的插值方法都要求求x是單調(diào)的，并且是單調(diào)的，并且xi不能夠不能夠超過(guò)超過(guò)x的范圍的范圍例題：從例題：從1 1點(diǎn)點(diǎn)1212點(diǎn)點(diǎn)的的1111小時(shí)內(nèi)，每隔小時(shí)內(nèi)，每隔1 1小時(shí)測(cè)量一次小時(shí)測(cè)量一次溫度，測(cè)得的溫度的數(shù)值依次為：溫度，測(cè)得的溫度的數(shù)值依次為：5 5，8 8，9 9，1515，2525，2929，3131，30

32、30，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); plot(hours,temps,+,h,t,r:)xlabel(Hour),ylabel(Degrees Celsius)0246810125101520253035HourDegrees Celsiusxy機(jī)翼下輪廓線例題例題. . 已知飛機(jī)下輪廓線上數(shù)據(jù)如下，求已知飛機(jī)下輪廓線上數(shù)據(jù)如下，求x每改變每

33、改變0.1時(shí)的時(shí)的y值值051015-20020lagrange051015024piecewise linear051015024spline第一種（網(wǎng)格節(jié)點(diǎn)）第一種（網(wǎng)格節(jié)點(diǎn)）第二種（散亂節(jié)點(diǎn)）第二種（散亂節(jié)點(diǎn)） 已知已知 m n個(gè)網(wǎng)格節(jié)點(diǎn)個(gè)網(wǎng)格節(jié)點(diǎn) (,) (1,2,.,;1,2, )ijijxyzim jn其中其中(,)ijx y互不相同，不妨設(shè)互不相同，不妨設(shè)bxxxam 21dyyycn 21 構(gòu)造一個(gè)二元函數(shù)構(gòu)造一個(gè)二元函數(shù)),(yxfz 通過(guò)全部已知節(jié)點(diǎn)通過(guò)全部已知節(jié)點(diǎn),即即再用再用),(yxf計(jì)算插值，即計(jì)算插值，即).,(*yxfz (,), (0,1,;0,1,)ijij

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

35、值即為所求將四個(gè)插值點(diǎn)（矩形的四個(gè)頂點(diǎn)）處的函數(shù)值依次簡(jiǎn)記為： 分片線性插值分片線性插值xy (xi, yj)(xi, yj+1)(xi+1, yj)(xi+1, yj+1)Of (xi, yj)=f1，f (xi+1, yj)=f2，f (xi+1, yj+1)=f3，f (xi, yj+1)=f4插值函數(shù)為：11()jjijiiyyyxxyxx1213211( , )()()jiiijjyyxxf x yfffffxxyy第二片(上三角形區(qū)域)：(x, y)滿足11()jjiiiiyyyxxyxx插值函數(shù)為：1413411( , )()()jijjijyyxxf x yfffffyyxy注

36、意注意：(x, y)當(dāng)然應(yīng)該是在插值節(jié)點(diǎn)所形成的矩形區(qū)域內(nèi)顯然，分片線性插值函數(shù)是連續(xù)的；分兩片的函數(shù)表達(dá)式如下：第一片(下三角形區(qū)域)： (x, y)滿足雙線性插值是一片一片的空間二次曲面構(gòu)成雙線性插值函數(shù)的形式如下：( , )()()f x yaxb cyd其中有四個(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要求要求x0, ,y0單調(diào)；單調(diào)；x，y可取可取為矩陣，或?yàn)榫仃?，或x取行取行向量，向量，y取為列向量，取為列向量，x,y的值分別不能超出的值

37、分別不能超出x0, ,y0 0的范圍的范圍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 最鄰近插值；最鄰近插值；linear 雙線性插值；雙線性插值；cubic 雙三次插值；雙三次插值；缺省時(shí)缺省時(shí) 雙線性插值雙線性插值. .例題：測(cè)得平板表面例題：測(cè)得平板表面3 35 5網(wǎng)格點(diǎn)處的溫度分別為：網(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 8484 82 85

38、86 82 85 86 試作出平板表面的溫度分布曲面試作出平板表面的溫度分布曲面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)畫(huà)出原始數(shù)據(jù)，畫(huà)出粗糙的溫度分布曲線圖.2以平滑數(shù)據(jù),在 x、y方向上每隔0.2個(gè)單位的地方進(jìn)行插值.再輸入以下命令:xi=1:0.2:5;yi=1:0.2:3;zi=interp2(x,y,temps,xi,yi,cubic);mesh(xi,yi,zi)畫(huà)出插值后的溫度分布曲面圖. 1

39、234511.522.53606570758085901234511.522.5360657075808590通過(guò)此例對(duì)最近鄰點(diǎn)插值、雙線性插值方法和雙三次插值方法的插值效果進(jìn)行比較x=0:400:5600;x=0:400:5600;y=0:400:4800;y=0:400:4800;z=370 470 550 600 670 690 670 620 580 450 400 300 100 150 250;.z=370 470 550 600 670 690 670 620 580 450 400 300 100 150 250;. 510 620 730 800 850 870 850 78

40、0 720 650 500 200 300 350 320;. 510 620 730 800 850 870 850 780 720 650 500 200 300 350 320;. 650 760 880 970 1020 1050 1020 830 900 700 300 500 550 480 350;. 650 760 880 970 1020 1050 1020 830 900 700 300 500 550 480 350;. 740 880 1080 1130 1250 1280 1230 1040 900 500 700 780 750 650 550;. 740 880

41、1080 1130 1250 1280 1230 1040 900 500 700 780 750 650 550;. 830 980 1180 1320 1450 1420 1400 1300 700 900 850 840 380 780 750;. 830 980 1180 1320 1450 1420 1400 1300 700 900 850 840 380 780 750;. 880 1060 1230 1390 1500 880 1060 1230 1390 1500 15001500 1400 900 1100 1060 950 870 900 930 950;. 1400 9

42、00 1100 1060 950 870 900 930 950;. 910 1090 1270 1500 1200 1100 1350 1450 1200 1150 1010 880 1000 1050 1100;. 910 1090 1270 1500 1200 1100 1350 1450 1200 1150 1010 880 1000 1050 1100;. 950 1190 1370 1500 1200 1100 1550 1600 1550 1380 1070 900 1050 1150 1200;. 950 1190 1370 1500 1200 1100 1550 1600 1

43、550 1380 1070 900 1050 1150 1200;. 1430 1430 14301430 1460 1500 1550 1600 1550 1600 1460 1500 1550 1600 1550 1600 16001600 16001600 1550 1500 1550 1500 15001500 1550 1550 15501550;.;. 1420 1430 1450 1480 1500 1550 1510 1430 1300 1200 980 850 750 550 500;. 1420 1430 1450 1480 1500 1550 1510 1430 1300

44、 1200 980 850 750 550 500;. 1380 1410 1430 1450 1470 1320 1280 1200 1080 940 780 620 460 370 350;. 1380 1410 1430 1450 1470 1320 1280 1200 1080 940 780 620 460 370 350;. 1370 1390 1410 1430 1440 1140 1110 1050 950 820 690 540 380 300 210;. 1370 1390 1410 1430 1440 1140 1110 1050 950 820 690 540 380

45、300 210;. 1350 1370 1390 1400 1410 960 940 880 800 690 570 430 290 210 150; 1350 1370 1390 1400 1410 960 940 880 800 690 570 430 290 210 150; figure(1); figure(1); meshz(x,y,zmeshz(x,y,z) ) xlabel(X),ylabel(Y),zlabel(Zxlabel(X),ylabel(Y),zlabel(Z)xi=0:50:5600;xi=0:50:5600;yiyi=0:50:4800;=0:50:4800;f

46、igure(2)figure(2)z1i=interp2(x,y,z,xi,yi,nearest);z1i=interp2(x,y,z,xi,yi,nearest);surfc(xi,yi,z1i)surfc(xi,yi,z1i)xlabel(X),ylabel(Y),zlabel(Zxlabel(X),ylabel(Y),zlabel(Z)figure(3)figure(3)z2i=interp2(x,y,z,xi,yi);% z2i=interp2(x,y,z,xi,yi);% 缺省時(shí)為雙線性缺省時(shí)為雙線性surfc(xi,yi,z2i)surfc(xi,yi,z2i)xlabel(X),

47、ylabel(Y),zlabel(Zxlabel(X),ylabel(Y),zlabel(Z)figure(4)figure(4)z3i=interp2(x,y,z,xi,yi,cubic);z3i=interp2(x,y,z,xi,yi,cubic);surfc(xi,yi,z3i)surfc(xi,yi,z3i)xlabel(X),ylabel(Y),zlabel(Zxlabel(X),ylabel(Y),zlabel(Z)figure(5)figure(5)subplot(1,3,1),contour(xi,yi,z1i,10);% subplot(1,3,1),contour(xi,yi,z1i,10);% 畫(huà)等高線畫(huà)等高線subplot(1,3,2),contour(xi,yi,z2i,10);subplot(1,3,2),contour(xi,yi,z2i,10);subplot(1,3,3),contour(xi,yi,z3i,10);subplot(1,3,3),contour(xi,yi,z3i,10);figure(6)figure(6)subplot(1,3,1),contourf(xi,yi,z1i,10);% subplot(1,3,1),contourf(xi,yi,z1i,10);% 畫(huà)等高面畫(huà)等高面subplot(1,3,2),