實驗六、擬合與插值問題

1、2021/6/161數學模型數學模型數學實驗數學實驗 六六一、實驗項目名稱一、實驗項目名稱 傳染病模型傳染病模型 插值與擬合插值與擬合2021/6/162進一步鞏固、加強進一步鞏固、加強Matlab的應用能力、的應用能力、學會用學會用MATLAB軟件進行數據擬合軟件進行數據擬合了解在最小二乘意義下數據擬合的理了解在最小二乘意義下數據擬合的理論和方法論和方法.通過對實際問題的分析和研究，初步通過對實際問題的分析和研究，初步掌握建立數據擬合數學模型的方法掌握建立數據擬合數學模型的方法插值的基本原理插值的基本原理二、實驗目的二、實驗目的2021/6/163三、實驗內容與步三、實驗內容與步驟驟1 1、

2、建模實例：傳染、建模實例：傳染病模型病模型2 2、MATLABMATLAB求函求函數的擬合與插值數的擬合與插值2021/6/164 目目 的的1、傳染病模型、傳染病模型 描述傳染病的傳播過程描述傳染病的傳播過程 分析受感染人數的變化規(guī)律分析受感染人數的變化規(guī)律 預報傳染病高潮到來的時刻預報傳染病高潮到來的時刻 預防傳染病蔓延的手段預防傳染病蔓延的手段 按照傳播過程的一般規(guī)律，按照傳播過程的一般規(guī)律，用機理分析方法建立模型用機理分析方法建立模型2021/6/1651.1問題重述問題重述n問題: 有一種傳染病（如SARS、甲型H1N1）正在流行。現在希望建立適當的數學模型，利用已經掌握的一些數據資

3、料對該傳染病進行有效地研究，以期對其傳播蔓延進行必要的控制，減少人民生命財產的損失。考慮如下的幾個問題，建立適當的數學模型，并進行一定的比較分析和評價展望。n1、不考慮環(huán)境的限制，設單位時間內感染人數的增長率是常數，建立模型求t時刻的感染人數。n2、假設環(huán)境條件下所允許的最大可感染人數為 。單位時間內感染人數的增長率是感染人數的線性函數，最大感染時的增長率為零。建立模型求t時刻的感染人數。mx2021/6/1663、現有衛(wèi)生防疫部門采集到的某地區(qū)一定時間內一定間隔區(qū)間的感染人數數據（見下表），利用該數據確定上述兩個模型中的相關參數，并將它們的預測值與實際數據進行比較分析（計算仿真偏差）并對兩個

4、模型進行適當的評價。（注：該問題中，設最大可感染人數為2000人）2021/6/167實驗問題實驗問題2021/6/168模型（一）模型（一）模型（二）模型（二）Matlab-Matlab-微分方程的求解微分方程的求解已解決已解決今天我們研究問今天我們研究問題題4 4的方法的方法擬合擬合插值插值2021/6/169據人口統(tǒng)計年鑒，知我國從據人口統(tǒng)計年鑒，知我國從19491949年至年至19941994年人年人口數據資料如下：口數據資料如下： ( (人口數單位為：百萬人口數單位為：百萬) )（1 1）在直角坐標系上作出人口數的圖象。）在直角坐標系上作出人口數的圖象。（2 2）建立人口數與年份的函

5、數關系，并估算）建立人口數與年份的函數關系，并估算19991999年年的人口數。的人口數。擬擬 合合 問問 題題 實實 例例 1年份年份19491954 1959 1964 1969人口數人口數 541.67602.66 672.09 704.99 806.71 年份年份 1974 1979 1984 1989 1994人口數人口數 908.59 975.42 1034.751106.761176.74 2021/6/1610如何確定如何確定21axay線性模型線性模型21,aa2021/6/16111 曲線擬合問題的提法曲線擬合問題的提法: 已知一組（二維）數據，即平面上的已知一組（二維）數

6、據，即平面上的n 個點個點),(iiyx， ixni, 2 , 1L互不相同，尋求一個函數（曲線）互不相同，尋求一個函數（曲線）)(xfy ，使使)(xf在某種準則下與所有數據點最為接近，即曲線擬合在某種準則下與所有數據點最為接近，即曲線擬合得最好，如圖得最好，如圖: xy0+i ),(iiyx)(xfy 一、曲線擬合一、曲線擬合 niiiniiyxf1212)( 確定確定f(x)使得使得 達到最小達到最小 最小二乘準則最小二乘準則 2021/6/1612. 用什么樣的曲線擬合已知數據用什么樣的曲線擬合已知數據?常用的曲線函數系類型：常用的曲線函數系類型：畫圖觀察；畫圖觀察；理論分析理論分析x

7、aeay21 指數曲線：指數曲線： 21axay 雙曲線（一支雙曲線（一支):): 11 mmmaxaxayL多項式：多項式： 21axay直線：直線： 2021/6/1613線性最小二乘擬合線性最小二乘擬合 f(x)=a1r1(x)+ +amrm(x)中中函數函數 r1 1( (x),), ,rm( (x)的選取的選取 1. 1. 通過機理分析建立數學模型來確定通過機理分析建立數學模型來確定 f( (x) )；+f=a1+a2xf=a1+a2x+a3x2f=a1+a2x+a3x2f=a1+a2/xf=aebxf=ae-bx 2. 2. 將數據將數據 (xi,yi) i=1, ,n 作圖，通過

8、直觀判斷確定作圖，通過直觀判斷確定 f(x)：2021/6/1614 擬合函數組中系數的確定擬合函數組中系數的確定達達到到最最小小。 niiiniiyaxaaaJ12211221),( 21,aa,)(21為例為例以以axaxf 使使得得，即即確確定定21aa 0) 20) 2121121niiiiniiiyaxaxyaxa2021/6/1615用用MATLAB作線性最小二乘擬合作線性最小二乘擬合1. 1. 作多項式作多項式f(x)=a1xm+ +amx+am+1擬合擬合, ,可利用已有程序可利用已有程序:a=polyfit(x,y,m)2.2.多項式在多項式在x處的值處的值y可用以下命令計算

9、：可用以下命令計算： y=polyval（a，x）輸出擬合多項式系數輸出擬合多項式系數a=a1, ,am , am+1 (數組）數組）輸入同長度輸入同長度的數組的數組x，y擬合多項擬合多項式次數式次數2021/6/1616人口預測線性模型人口預測線性模型對于開始提出的實驗問題對于開始提出的實驗問題, 代如數據，計算得代如數據，計算得27754,1521 aa從而得到人口數與年份的函數關系為從而得到人口數與年份的函數關系為2775415 xy把把x=1999代如，估算出代如，估算出1999年的人口數為年的人口數為 y=1252.1（百萬）（百萬）12.52億億1999年實際人口數量為年實際人口數

10、量為.億。億。線性預測模型線性預測模型2021/6/1617syms xx=1949:5:1994;y=541.67 602.66 672.09 704.99 806.71 908.59 975.42 1034.75 1106.76 1176.74;plot(x,y,r*)z=polyfit(x,y,1)y1=polyval(z,x)hold onplot(x,y1,b+-)19451950195519601965197019751980198519901995500600700800900100011001200人口模型的解人口模型的解2021/6/1618nsyms xnx=0:1:14n

11、y=39 53 72 96 129 171 232 314 386 502 629 760 920 1065 1232;nplot(x,y,r*)nz=polyfit(x,y,2)ny1=polyval(z,x)nhold onnplot(x,y1,b+-)傳染病模型傳染病模型3 3的解的解0246810121402004006008001000120014002021/6/1619二、插二、插 值值 2. 面面 對對 一一 個個 實實 際際 問問 題，應題，應 該該 用用 插插 值，還值，還 是是 擬擬 合。合。1. 插值的基本原理；三種插值方法：拉格朗日插值的基本原理；三種插值方法：拉格朗

12、日插值，分段線性插值，三次樣條插值。插值，分段線性插值，三次樣條插值。2021/6/1620插插 值值 問問 題題 實實 例例 1 1xy機翼下輪廓線2021/6/1621插插 值值 問問 題題 的的 提提 法法已知已知 n+1個節(jié)點個節(jié)點, 1 , 0(),(njyxjjL其中其中jx互不相同，不妨設互不相同，不妨設),10bxxxanL求任一插值點求任一插值點)(*jxx 處的插值處的插值.*y0 x1xnx0y1y節(jié)點可視為由節(jié)點可視為由)(xgy 產生產生,g表達式復雜表達式復雜或無封閉形式或無封閉形式,，或未知或未知.。*x*y2021/6/16220 x1xnx0y1y求求 解解

13、插插 值值 問問 題題 的的 基基 本本 思思 路路 構造一個構造一個(相對簡單的相對簡單的)函數函數),(xfy 通過全部節(jié)點通過全部節(jié)點, 即即), 1 ,0()(njyxfjjL再用再用)(xf計算插值，即計算插值，即).(*xfy *x*y2021/6/1623插值的基本原理插值的基本原理選講選講PPT（24-30）2021/6/16241.1.拉格朗日(Lagrange)多項式插值1.1 1.1 插值多項式插值多項式) 1 ()(0111axaxaxaxLnnnnnLnnnnnnnnyyYaaAxxxxXLLL001100,11)2(,0)det(X有唯一解有唯一解), 1 , 0(

14、)(njyxLjjnL)2(YXA 2021/6/16251.1.拉格朗日(Lagrange)多項式插值1.2 拉格朗插值多項式nixxxxxxxxxxxxxxxxxlniiiiiiniiiLLLLL1 , 0,)()()()()()()(110110)3()()(0 xlyxLiniinjjnjiyxLjijixl)(,0, 1)(又（2）有唯一解，故（3）與（1）相同。2021/6/16261.1.拉格朗日(Lagrange)多項式插值1.3 1.3 誤差估計誤差估計),(),()!1()()()()(0)1(baxxngxLxgxRnjjnnn1)1()(nnMgnjjnnxxnMxR0

15、1)!1()()(xRnn)(xRgn光滑)(xRxxnj接近2021/6/16271.1.拉格朗日(Lagrange)多項式插值1.4 1.4 例例 將0,/2 n等分，用g(x)=cos(x)產生n+1個節(jié)點，作Ln(x)（取n=1,2) ,計算cos(/6), 估計誤差。解解: n=1, (x0,y0)=(0,1), (x1,y1)=(/2,0), L1(x)=y0l0+y1l1=1-2x/, cos(cos( /6)=0.6667/6)=0.6667n=2, (x0,y0)=(0,1), (x1,y1)=(/4,0.7071),(x2,y2)=(/2,0),L2(x)=y0l0+y1l

16、1+y2l2=8(x-/4)(x-/2)/2-16x(x-/2)0.7071/2cos(cos( /6)=L/6)=L2 2( ( /6)=0.8508 /6)=0.8508 精確值：精確值：cos (cos ( /6)=0.8660/6)=0.8660nhhhhxxxxxnhMxRnjjjjnnL324,2, 1:)(2011112)2)(1(4324)!1(1)(nnnnnnhhhhnxRLn1234)( xRn0 .30 .0 44 .7 1 0-34 .7 1 0-42021/6/16281.1.拉格朗日(Lagrange)多項式插值1.51.5 拉格朗日插值多項式的振蕩?)(?)(x

17、RxLnnn55,11)(2xxxg63. 363. 3),()(limxxgxLnnTo MATLAB(runge)Runge現象現象: :-505-1.5-1-0.500.511.52y=1/(1+x2)n=2n=4n=6n=8n=102021/6/16292.2.分段線性插值分段線性插值xjxj-1xj+1x0 xn其它,0,)()()(1111110jjjjjjjjjjjnjjjnxxxxxxxxxxxxxxxlxlyxI計算量與n無關;n越大，誤差越小.nnnxxxxgxI0),()(lim2021/6/16303. 3. 三次樣條插值三次樣條插值, 1,),()(1nixxxxsx

18、SiiiL,)()3), 1 ,0()()2), 1()()10223niiiiiiixxCxSniyxSnidxcxbxaxsLL) 1, 1()()(),()(),()(111 nixsxsxsxsxsxsiiiiiiiiiiii自然邊界條件）(0)()()40 nxSxS)(,)4)3)2xSdcbaiiii)()(limxgxsn2021/6/16311. 1. 拉格朗日插值拉格朗日插值: :自編程序自編程序, ,如名為如名為 lagr1.m 的的M文件，文件， 第一行為第一行為 function y=lagr1(x0,y0,x) 輸入輸入: :節(jié)點節(jié)點x0,y0, 插值點插值點x (

19、 (均為數組，長度自定義均為數組，長度自定義) )）；）； 輸出輸出: :插值插值y ( (與與x同長度數組同長度數組) )）。）。 應用時輸入應用時輸入x0,y0,x后后, ,運行運行 y=lagr1(x0,y0,x)2. 2. 分段線性插值分段線性插值: :已有程序已有程序 y=interp1(x0,y0,x)3. 3. 三次樣條插值三次樣條插值: :已有程序已有程序 y=interp1(x0,y0,x, spline) 或或 y=spline(x0,y0,x)用MATLAB作插值計算2021/6/1632對表格給出的函數，求出沒有給出的函數值 例例2表中是待加工零件下輪廓線的一組數據，現

20、需要得到x坐標每改變0.1時所對應的y的坐標. x03571112131415y01.21.72.12.01.81.21.01.6下面是關于插值的兩條命令(專門用來解決這類問題)：y=interp1(x0,y0,x) 分段線性插值y=spline(x0,y0,x) 三次樣條插值 2021/6/1633n其中x0,y0是已知的節(jié)點坐標，是同維向量。y對應于x處的插值。y與x是同維向量。n解決上述問題,我們可分兩步:n一 用原始數據繪圖作為選用插值方法的參考.n二 確定插值方法進行插值計算x0=0,3,5,7,9,11,12,13,14,15;y0=0,1.2,1.7,2.0,2.1,2.0,1.

21、8,1.2,1.0,1.6plot(x0,y0) %完成第一步工作x=0:0.1:15;y=interp1(x0,y0,x); %用分段線性插值完成第二步工作plot(x,y)y=spline(x0,y0,x); plot(x,y) %用三次樣條插值完成第二步工作2021/6/1634n對y=1/(1+x2),-5x5，用n（=11）個節(jié)點（等分）作上述兩種插值，用m（=21）個插值點（等分）作圖，比較結果。n解：鍵入并運行如下命令nn=11;m=21;x=-5:10/(m-1):5;y=1./(1+x.2);nxo=-5:10/(n-1):5;yo=1./(1+xo.2);ny1=inter

22、p1(xo,yo,x);ny2=spline(xo,yo,x);nplot(x,y,r,x,y1,b,x,y2,k)練習：2021/6/1635例 3 一水庫上游河段降暴雨.，根據預報測算上游流入水庫的流量為Q(t) (102立方米/秒) ： t (時) 8 12 16 24 30 44 48 56 60Q( t ) 36 54 78 92 101 35 25 16 13 利用這個預報值估計14：30 和 20：30 時上游流入水庫的流量。假設: 1 已知數據準確。 2 相鄰兩個時刻之間的流量沒有突然的變化。2021/6/16362021/6/1637nt=8,12,16,24,30,44,4

23、8,56,60;nq=36,54,78,92,101,35,25,16,13;nt1=8:0.5:60;nq1=interp1(t,q,t1,linear);nplot(t,q,b,t1,q1);nhold on;nq2=interp1(t,q,t1,spline);nplot(t,q,b,t1,q1,t1,q2,r)nq1nq22021/6/16382021/6/1639擬合與插值的關系擬合與插值的關系 函數插值與曲線擬合都是要根據一組數據構造一個函數作函數插值與曲線擬合都是要根據一組數據構造一個函數作為近似，由于近似的要求不同，二者在數學方法上是完全不同為近似，由于近似的要求不同，二者在數

24、學方法上是完全不同的的 實例：實例：下面數據是某次實驗所得，希望得到X和 f之間的關系？x124791 21 31 51 7f1 .53 .96 .611 .71 5 .61 8 .81 9 .62 0 .62 1 .1問題：問題：給定一批數據點，需確定滿足特定要求的曲線或曲面解決方案：解決方案：若不要求曲線（面）通過所有數據點，而是要求它反映對象整體的變化趨勢，這就是數據擬合數據擬合，又稱曲線擬合或曲面擬合若要求所求曲線（面）通過所給所有數據點，就是插值問題插值問題；2021/6/1640最臨近插值、線性插值、樣條插值與曲線擬合結果：最臨近插值、線性插值、樣條插值與曲線擬合結果：0246810121416180510152025已知數據點spline三次多項式插值0246810121416180510152025已知數據點linest三次多項式插值0246810121416180510152025已知數據點nearest三次多項式插值2021/6/1641即要求即要求 出二次多項式出二次多項式:3221)(axaxaxf中中 的的),(321aaaA 使得使得:最小 )(1112iiiyxf1、對