數(shù)學建模插值法與曲線擬合講課

1、 插值與曲線擬合插值與曲線擬合-5-4-3-2-1012345-0.500.511.52一、問題的提出 在生產(chǎn)和實驗中，關于函數(shù)f(x)，經(jīng)常存在兩種情況：（1）其表達式不便于計算；（2）無表達式. 而只有函數(shù)在給定點的函數(shù)值，怎樣預測其它點的函數(shù)值？xx0 x1x2xnyy0y1y2yn飛機機翼制造 下表給出的下表給出的x x、y y數(shù)據(jù)位于機翼端面的輪廓線上，數(shù)據(jù)位于機翼端面的輪廓線上，Y1Y1和和Y2Y2分分別對應輪廓的上下線。假設需要得到別對應輪廓的上下線。假設需要得到x x坐標每改變坐標每改變0.10.1時的時的y y坐標，試完成加工所需數(shù)據(jù)，畫出曲線坐標，試完成加工所需數(shù)據(jù)，畫出曲

2、線. .x035791112131415Y101.82.22.73.03.12.92.52.01.6Y201.21.72.02.02.01.81.21.01.6山體地貌n要在某山區(qū)方圓大約27平方公里范圍內(nèi)修建一條公路，從山腳出發(fā)經(jīng)過一個居民區(qū)，再到達一個礦區(qū)。橫向縱向分別每隔400米測量一次，得到一些地點的高程： 山區(qū)地貌：山區(qū)地貌： 在某山區(qū)測得一些地點的高程如下表。平面區(qū)域為在某山區(qū)測得一些地點的高程如下表。平面區(qū)域為 1200=x=4000,1200=y=3600) 試作出該山區(qū)的地貌圖和等高線圖，并對幾種插值方法進行比較。試作出該山區(qū)的地貌圖和等高線圖，并對幾種插值方法進行比較。 X

3、 Y 1200 1600 2000 2400 2800 3200 3600 4000 1200 1130 1250 1280 1230 1040 900 500 700 1600 1320 1450 1420 1400 1300 700 900 850 2000 1390 1500 1500 1400 900 1100 1060 950 2400 1500 1200 1100 1350 1450 1200 1150 1010 2800 1500 1200 1100 1550 1600 1550 1380 1070 3200 1500 1550 1600 1550 1600 1600 1600

4、1550 3600 1480 1500 1550 1510 1430 1300 1200 980 試做出該山區(qū)的地貌圖.船在該海域會擱淺嗎？船在該海域會擱淺嗎？-作業(yè)作業(yè) 在某海域測得一些點(x,y)處的水深z由下表給出，船的吃水深度為5英尺，在矩形區(qū)域（75，200）*（-50，150）里的哪些地方船要避免進入.xyz129 140 103.5 88 185.5 195 1057.5 141.5 23 147 22.5 137.5 85.54 8 6 8 6 8 8xyz157.5 107.5 77 81 162 162 117.5-6.5 -81 3 56.5 -66.5 84 -33.5

5、9 9 8 8 9 4 9水深和流速的問題 在水文數(shù)據(jù)測量中，不同水深的流速是不同的. 水文數(shù)據(jù)的測量時天天進行的，為了減少測量的工作，希望得到確定的水深和水流之間的關系. 為此測量了一系列不同水深和流速值. 下表給出了對某河流的測量數(shù)據(jù)，其中水深和流速根據(jù)適當?shù)膯挝贿M行了規(guī)范化，共10個值.2.983.063.133.203.193.233.253.263.223.19流速0.90.80.70.60.50.40.30.20.10水深美國人口問題n據(jù)美國人口普查局數(shù)據(jù)： 從1790每隔10年至2000年的總?cè)丝冢▎挝唬喊偃f）如下示n t = 1790:10:2000; np = 3.9, 5.

6、3 , 7.2 , 9.6 , 12.9 , 17.1 , 23.1 , 31.4 , 38.6 , 50.2 , 62.9 , 76 , 92 , 105.7 , 122.8 , 131.7 , 150.7 , 179 , 205 , 226.5 , 251.4 , 281.422; n預測2001，2002年的美國人口數(shù)？并與調(diào)查數(shù)據(jù)285.318，288.369比較，選擇擬合較好的模型。農(nóng)作物施肥效果分析1992年A題 在農(nóng)業(yè)生產(chǎn)試驗研究中，對某地區(qū)土豆的產(chǎn)量與化肥的關系做了一實驗，得到了氮肥、磷肥的施肥量與土豆產(chǎn)量的對應關系如下表： 1.根據(jù)上表數(shù)據(jù)分別給出土豆產(chǎn)量與氮、磷肥的關系式。

7、 2.施肥問題優(yōu)化策略氮肥量（公斤/公頃）03467101135202259336404471土豆產(chǎn)量（公斤）15.1821.3625.7232.293439.4543.1543.4640.8330.75磷肥量（公斤/公頃）024497398147196245294342土豆產(chǎn)量（公斤）33.4632.4736.0637.964140.141。342.240.442.7配藥方案-作業(yè) 一種新藥用于臨床之前, 必須設計給藥方案. 在快速靜脈注射的給藥方式下, 所謂給藥方案是指, 每次注射劑量多大, 間隔時間多長. 藥物進入機體后隨血液輸送到全身, 在這個過程中不斷地被吸收, 分布, 代謝, 最終

8、排出體外. 藥物在血液中的濃度, 即單位體積血液中的藥物含量, 稱血藥濃度. 在最簡單的一室模型中, 將整個機體看作一個房室, 稱中心室, 室內(nèi)的血藥濃度是均勻的. 快速靜脈注射后, 濃度立即上升; 然后逐漸下降. 當濃度太低時, 達不到預期的治療效果; 血藥濃度太高, 又可能導致藥物中毒或副作用太強. 臨床上, 每種藥物有一個最小有效濃度 c1 和一個最大治療濃度 c2. 設計給藥方案時, 要使血藥濃度保持在 c1-c2 之間. 設本題所研究藥物的最小有效濃度c1=10, 最大治療濃度 c2=25( ). /g ml 顯然, 要設計給藥方案, 必須知道給藥后血藥濃度隨時間變化的規(guī)律. 為此,

9、 從實驗和理論兩方面著手. 在實驗方面, 對某人用快速靜脈注射方式一次注入該藥物300mg后, 在一定時刻 t (小時)采集血樣, 測得血藥濃度c. 如表: 血藥濃度c(t) 的測試數(shù)據(jù) t0.250.511.523468c19.2118.1515.3614.1012.899.327.455.243.01問題問題 ：1. 1. 在快速靜脈注射的給藥方式下，研究血藥濃度（單位體積血液中的藥物含在快速靜脈注射的給藥方式下，研究血藥濃度（單位體積血液中的藥物含量）的變化規(guī)律；量）的變化規(guī)律；2. 2. 給定藥物的最小有效濃度和最大治療濃度，設計給藥方案：每次注射劑量給定藥物的最小有效濃度和最大治療濃

10、度，設計給藥方案：每次注射劑量多大；間隔時間多長？多大；間隔時間多長？配藥方案二、問題的解決（1）插值法；（2）曲線擬合法 1、問題的抽象、問題的抽象xx1x2xmyy1y2ym構(gòu)造一個構(gòu)造一個簡單易于計算簡單易于計算的近似函數(shù)的近似函數(shù) p(x) f (x) （精確函數(shù)）。（精確函數(shù)）。2、構(gòu)造近似函數(shù)，、構(gòu)造近似函數(shù)， p(x) 的方法有兩種：的方法有兩種：在實驗中經(jīng)常給出一組離散點，在實驗中經(jīng)常給出一組離散點，插值法定義：定義：當精確函數(shù)當精確函數(shù) y = f(x) 非常復雜非常復雜或或未知時未知時，在一系列節(jié)點，在一系列節(jié)點 x0 xn 處測得函數(shù)值處測得函數(shù)值 y0 = f(x0),

11、 ，yn = f(xn)，由此構(gòu)造一個由此構(gòu)造一個簡單易算簡單易算的近似函數(shù)的近似函數(shù) p(x) f(x)，滿足條件，滿足條件p(xi) = f(xi) (i = 0, n)，（插值條件），（插值條件）這里的這里的 p(x) 稱為稱為f(x) 的的插值函數(shù)；插值函數(shù)； 構(gòu)造插值函數(shù)的方法為構(gòu)造插值函數(shù)的方法為插值法插值法。曲線擬合 但是不要求使但是不要求使 p(xi) = yi ,而只要而只要 p(xi) yi 總體上總體上盡可能小。這種構(gòu)盡可能小。這種構(gòu)造近似函數(shù)造近似函數(shù)p(x) 的方法稱為的方法稱為曲線擬合法曲線擬合法， p(x) 稱為稱為擬合函數(shù)。擬合函數(shù)。定義：定義： 當精確函數(shù)當精

12、確函數(shù) y = f(x) 非常復雜非常復雜或或未知時未知時，在一系列節(jié)點，在一系列節(jié)點x0 xn 處處,測得函數(shù)值測得函數(shù)值 y0 , ，yn ，由此構(gòu)造一個，由此構(gòu)造一個簡單易簡單易 算算的近似函的近似函數(shù)數(shù) p(x) f(x)，插值與擬合的相同點n都需要根據(jù)已知數(shù)據(jù)構(gòu)造函數(shù)。n可使用得到函數(shù)計算未知點的函數(shù)值。xx1x2xmyy1y2ym求一個求一個簡單易算簡單易算的近似函數(shù)的近似函數(shù) p(x) f (x) 。插值與擬合的不同點n插值: 過節(jié)點; ；n擬合: 不過點, 整體近似；插值法n拉格朗日插值 n牛頓插值 n三次埃爾米特插值法n分段線性插值n分段三次埃爾米特插值法n三次樣條插值1、

13、拉格朗日插值公式（）定義對給定的n+1個節(jié)點x0 , x1,x2,xn及對應的函數(shù)值y0 , y1,y2,yn，構(gòu)造一個n次插值多項式：即為拉格朗日插值公式拉格朗日插值公式，其中)(0 xlyyknkk)()()()()()()(110110nkkkkkknkkkxxxxxxxxxxxxxxxxxl插值基函插值基函數(shù)數(shù)拉格朗日插值的matlab實現(xiàn)function y=lagrange(x0,y0,x) % x0插值節(jié)點， y0插值節(jié)點處的函數(shù)值，x要計算函數(shù)值的點；n=length(x0); %計算x0的長度m=length(x); %計算x的長度for i=1:m s=0;z=x(i);

14、nfor k=1:nn p=1.0;n for j=1:nn if j=knp=p*(z-x0(j)/(x0(k)-x0(j); n%計算插值基函數(shù)n endn endn s=p*y0(k)+s;nendny(i)=s; %計算在x(i)處的函數(shù)值（拉格朗日）nend2、牛頓插值法牛頓插值公式：Nn(x)=f(x0)+fx0,x1(x-x0)+fx0,x1,x2,xn(x-x0)(x-x1)(x-xn)其中： fx0,x1 一階差商 fx0,x1,x2,xn n階差商n注：牛頓插值法與拉格朗日插值法，同一個多項式，不同的表達方式，但是計算量不一樣，牛頓插值法的計算量小。龍格現(xiàn)象 Runge在上

15、個世紀初發(fā)現(xiàn)： 在-5,5上用n+1個等距節(jié)點作n次插值多項式Pn(x)，當在n時，插值多項式Pn(x)在區(qū)間中部趨于f(x)=1/(1+x2) ,但對于3.63x1的x，Pn(x)嚴重發(fā)散。用圖形分析問題。 for n=10:2:20 %從10等份到20等份x0=-5:10/n:5; %插值節(jié)點y0=1./(1+x0.2); %插值節(jié)點處的精確函數(shù)值x=-5:0.1:5; %要進行計算函數(shù)值的點y=lagrange(x0,y0,x); %調(diào)用函數(shù)計算x點的函數(shù)值 plot(x0,y0,*,x,1./(1+x.2),r,x,y) %繪制圖形pause %等待，按任意鍵end3、分段低次插值法（

16、1）分段線性插值 定義： 已知n+1個不同節(jié)點x0,x1,xn ,構(gòu)造分段多項式I(x)，使之滿足l I(x)在a，b上連續(xù)；l I(xk)=yk；l I(x)在xi,xi+1上是一次多項式； I(x)=,11111kkkkkkkkkkxxxxxxxyxxxxy（2）分段三次埃爾米特插值法定義：已知n+1個不同節(jié)點x0,x1,xn ,構(gòu)造分段多項式I(x)，使之滿足：l I(x)在a，b上二階連續(xù)導數(shù)；l I(xk)=yk， I(xk)=yk， ；l I(x)在xi,xi+1上是三次次多項式。4、三次樣條插值法 對于給定n+1個不同節(jié)點x0,x1,xn及函數(shù)值y0,y1,yn， 其中a=x0

17、x1xn=b，構(gòu)造三次樣條插值函數(shù)S(x)。 S(x)稱為三次樣條函數(shù)時需滿足：l S(x)在a,b上二階導數(shù)連續(xù)；l S(xk)=yk (k=0,1,n)；l 每個子區(qū)間xk,xk+1上S(x)是三次多項式(k=0,1,n)。插值法的matlab實現(xiàn)一維插值 命令：interp1(x0,y0,x，method) 其中：x0：插值節(jié)點； y0：插值節(jié)點處的函數(shù)值； x：要計算函數(shù)值的點； method： l i n e a r ：分段線性插值；：分段線性插值； c u b i c ：分段三次埃爾米特插值； s p l i n e ：三次樣條插值。插值法的應用n一水庫上游河段降暴雨，根據(jù)預報測算

18、上游流入水庫的流量為一水庫上游河段降暴雨，根據(jù)預報測算上游流入水庫的流量為Q(tQ(t) (10) (102 2立方米立方米/ /秒秒) ) ： t (時時) 8 12 16 24 30 44 48 56 60 Q（t） 36 54 78 92 101 35 25 16 13 通過這個預報值，分別用不同的數(shù)值方法插值法來估計通過這個預報值，分別用不同的數(shù)值方法插值法來估計14 14 和和2020時上游流入水庫的流量。時上游流入水庫的流量。二維插值的MATLAB實現(xiàn) 在MATLAB中，二維插值命令常用的有兩個， 1、一個是網(wǎng)格節(jié)點插值：、一個是網(wǎng)格節(jié)點插值： z=interp2(x0,y0,z0

19、,x,y,method) 其中，n z：被插值點處的函數(shù)值；n x0,y0,z0：插值節(jié)點， x0,y0為向量，z0是矩陣，其列數(shù)等于x0的長度，行數(shù)等于y0的長度；n x,y： 要計算函數(shù)值的點； interp1(x0,y0,x，method)山體地貌n要在某山區(qū)方圓大約27平方公里范圍內(nèi)修建一條公路，從山腳出發(fā)經(jīng)過一個居民區(qū)，再到達一個礦區(qū)。橫向縱向分別每隔400米測量一次，得到一些地點的高程： 山區(qū)地貌：山區(qū)地貌： 在某山區(qū)測得一些地點的高程如下表。平面區(qū)域為在某山區(qū)測得一些地點的高程如下表。平面區(qū)域為 1200=x=4000,1200=y=3600) 試作出該山區(qū)的地貌圖和等高線圖，并

20、對幾種插值方法進行比較。試作出該山區(qū)的地貌圖和等高線圖，并對幾種插值方法進行比較。 X Y 1200 1600 2000 2400 2800 3200 3600 4000 1200 1130 1250 1280 1230 1040 900 500 700 1600 1320 1450 1420 1400 1300 700 900 850 2000 1390 1500 1500 1400 900 1100 1060 950 2400 1500 1200 1100 1350 1450 1200 1150 1010 2800 1500 1200 1100 1550 1600 1550 1380 10

21、70 3200 1500 1550 1600 1550 1600 1600 1600 1550 3600 1480 1500 1550 1510 1430 1300 1200 980 試做出該山區(qū)的地貌圖，并對幾種插值法進行比較.n程序設計：程序設計：nclearnx0=1200:400:4000; ny0=1200:400:3600; nz0=1130 1250 1280 1230 1040 900 500 700;n 1320 1450 1420 1400 1300 700 900 850;n 1390 1500 1500 1400 900 1100 1060 950;n 1500 120

22、0 1100 1350 1450 1200 1150 1010;n 1500 1200 1100 1550 1600 1550 1180 1070;n 1500 1550 1600 1550 1600 1600 1600 1550;n 1480 1500 1550 1510 1430 1300 1200 980;nxi=1200:10:4000; %加密數(shù)據(jù)點nyi=1200:10:3600; nzil=interp2(x0,y0,z0,xi,yi,linear); %線性插值nzic=interp2(x0,y0,z0,xi,yi,cubic); %三次插值nzis=interp2(x0,y0

23、,z0,xi,yi,spline); %樣條插值nsubplot(2,2,1) nmesh(x0,y0,z0)nsubplot(2,2,2) nmesh(xi,yi,zil)nsubplot(2,2,3) nmesh(xi,yi,zic)nsubplot(2,2,4) nmesh(xi,yi,zis)二維插值的MATLAB實現(xiàn)2、另一個是離散數(shù)據(jù)節(jié)點的插值命令： z=griddata(x0,y0,z0,x,y,method) 其中，n z：被插值點處的函數(shù)值；n x0,y0,z0：插值節(jié)點， x0,y0,z0均為向量；n x,y：被插值點； method：插值方法，包括： n linear線性

24、插值；n cubic 三次插值；船在該海域會擱淺嗎？船在該海域會擱淺嗎？ 在某海域測得一些點(x,y)處的水深z由下表給出，船的吃水深度為5英尺，在矩形區(qū)域（75，200）*（-50，150）里的哪些地方船要避免進入.xyz129 140 103.5 88 185.5 195 1057.5 141.5 23 147 22.5 137.5 85.54 8 6 8 6 8 8xyz157.5 107.5 77 81 162 162 117.5-6.5 -81 3 56.5 -66.5 84 -33.59 9 8 8 9 4 9解決問題的步驟：1.作出測量點的分布圖；2.求出矩形區(qū)域（75，200）

25、*（-50，150）的細分網(wǎng)格節(jié)點之橫、縱坐標向量；3.利用MATLAB中的散點插值函數(shù)求網(wǎng)格節(jié)點的水深；4.作出海底曲面圖形和等高線圖；5.作出水深小于5的海域范圍.程序nclearnx0=129 140 103.5 88 185.5 195 105 157.5 107.5 77 81 162 162 117.5; ny0=7.5 141.5 23 147 22.5 137.5 85.5 -6.5 -81 3 56.5 -66.5 84 -33.5; nz0=-4 -8 -6 -8 -6 -8 -8 -9 -9 -8 -8 -9 -4 -9; nsubplot(2,2,1)nplot(x0,

26、y0,+); %作出測量點的分布圖；作出測量點的分布圖； nx=75:1:200; %加密加密ny=-50:1:150;nx,y=meshgrid(x,y); nz=griddata(x0,y0,z0,x,y,cubic);nsubplot(2,2,2)nmesh(x,y,z), %用插值方法求出網(wǎng)格節(jié)點處的用插值方法求出網(wǎng)格節(jié)點處的z坐標矩陣坐標矩陣,繪制出三繪制出三維圖形維圖形nsubplot(2,2,3)nmeshc(x,y,z), %繪制等高線繪制等高線nsubplot(2,2,4)ncontour(x,y,z,-5 -5); %水深水深5英尺處海底曲面的等高線英尺處海底曲面的等高線n

27、grid on擬合的標準n（1）用各點誤差絕對值的和表示n（2）用各點誤差按絕對值的最大值表示n（3）用各點誤差的平方和表示 11()miiiRp xy1( )maxiii mRp xy 221()miiiRp xy最小二乘擬合n式中 R2 稱為均方誤差。由于計算均方誤差的最小值的原則容易實現(xiàn)而被廣泛采用。n按均方誤差達到極小構(gòu)造擬合曲線的方法稱為最小二乘法。221( ()miiiRp xy+p=a1+a2xp=a1+a2x+a3x2p=a1+a2x+a3x2p=a1+a2/xp=aebxp=ae-bx將數(shù)據(jù)將數(shù)據(jù) (xi,yi) i=1, ,n 作圖，通過直觀判斷確定作圖，通過直觀判斷確定

28、p(x)：MATLAB-曲線擬合工具箱nMatlab有一個功能強大的曲線擬合工具箱 （Curve Fitting Toolbox ） cftool， 使用方便，能實現(xiàn)多種類型的線性、非線性曲線擬合。n調(diào)用：cftooln界面如下所示“Data”按鈕數(shù)據(jù)的選取n點擊“Data”按鈕，彈出“Data”窗口；n利用X data和Y data的下拉菜單讀入數(shù)據(jù)x,y，可修改數(shù)據(jù)集名“Data set name”，然后點擊“Create data set”按鈕，退出“Data”窗口，返回工具箱界面，這時會自動畫出數(shù)據(jù)集的曲線圖；“Fitting”按鈕 曲線擬合n點擊“Fitting”按鈕，彈出“Fitt

29、ing”窗口；n點擊“New fit”按鈕，可修改擬合項目名稱“Fit name”，通過“Data set”下拉菜單選擇數(shù)據(jù)集，然后通過下拉菜單“Type of fit”選擇擬合曲線的類型。SSEnThe sum of squares due to error. This statistic measures the deviation of the responses from the fitted values of the responses. A value closer to 0 indicates a better fit.n偏差平方和，越接近0越好niiyxpSSE1i2)(R-

30、squaren The coefficient of multiple determination. This statistic measures how successful the fit is in explaining the variation of the data. A value closer to 1 indicates a better fit.n復相關系數(shù)平方(決定系數(shù))，越接近1越好Adjusted R-squaren The degree of freedom adjusted R-square. A value closer to 1 indicates a be

31、tter fit. It is generally the best indicator of the fit quality when you add additional coefficients to your model.n修正的復相關系數(shù)平方，越接近1越好Adjusted R-squaren下列公式中的m為擬合函數(shù)中待估參數(shù)個數(shù)，如：對一元一次多項式擬合，f(x) = a + bx，此時m=2，n為數(shù)據(jù)點個數(shù)。該修正類似修正的樣本方差使其為總體方差的無偏估計。RMSEnThe root mean squared error. A value closer to 0 indicates

32、 a better fit.n偏差平方的均值的算術平方根，越接近0越好曲線擬合好壞如何評價n首要指標是目標函數(shù)誤差最小（擬合度最大）；n其次是應考慮關鍵點的吻合，這些關鍵點包括：初始點（有時是原點）、拐點、峰值點、極值點、中間點、漸近點、終值點等，在這些關鍵點上，數(shù)據(jù)觀察值點與函數(shù)值點應盡可能一致；n再次是擬合的模型應盡可能簡單（模型的形式簡單，參數(shù)數(shù)少）。實踐中如何選擇模型？n在數(shù)據(jù)擬合實踐中，理性模型畢竟是少數(shù)，大多數(shù)的情形是根據(jù)數(shù)據(jù)的趨勢尋找合適的模型，有時好幾個模型對數(shù)據(jù)都有較好的擬合，但通過對關鍵點的比較總會找到一種最合適的模型。n在選擇不同的模型時，合理性和可解釋性是首要考慮的因素

33、。美國人口問題n據(jù)美國人口普查局數(shù)據(jù)： 從1790每隔10年至2000年的總?cè)丝冢▎挝唬喊偃f）如下示n t = 1790:10:2000; np = 3.9, 5.3 , 7.2 , 9.6 , 12.9 , 17.1 , 23.1 , 31.4 , 38.6 , 50.2 , 62.9 , 76 , 92 , 105.7 , 122.8 , 131.7 , 150.7 , 179 , 205 , 226.5 , 251.4 , 281.422; 美國人口問題nt = 1790:10:2000; np = 3.9, 5.3 , 7.2 , 9.6 , 12.9 , 17.1 , 23.1 ,

34、31.4 , 38.6 , 50.2 , 62.9 , 76 , 92 , 105.7 , 122.8 , 131.7 , 150.7 , 179 , 205 , 226.5 , 251.4 , 281.422; n預測2001，2002年的美國人口數(shù)？并與調(diào)查數(shù)據(jù)285.318，288.369比較，選擇擬合較好的模型。Matlab求解n在命令窗口輸入命令ncftooln回車，得擬合的圖形用戶界面結(jié)果分析nLinear model Poly2: f(x) = p1*x2 + p2*x + p3nCoefficients (with 95% confidence bounds):np1 = 0.006757 (0.006369, 0.007144)np2 = -24.32 (-25.78, -22.85)np3 = 2.188e+004 (2.049e+004, 2.327e+004)nGoodness of fit:n SSE: 184.5n R-square: