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

關(guān)于數(shù)學(xué)建模插值法與曲線擬合講課第一頁，共六十五頁，2022年，8月28日

一、問題的提出

在生產(chǎn)和實驗中，關(guān)于函數(shù)f(x)，經(jīng)常存在兩種情況：（1）其表達(dá)式不便于計算；（2）無表達(dá)式.

而只有函數(shù)在給定點的函數(shù)值，怎樣預(yù)測其它點的函數(shù)值？xx0x1x2…xnyy0y1y2…yn第二頁，共六十五頁，2022年，8月28日飛機(jī)機(jī)翼制造

下表給出的x、y數(shù)據(jù)位于機(jī)翼端面的輪廓線上，Y1和Y2分別對應(yīng)輪廓的上下線。假設(shè)需要得到x坐標(biāo)每改變0.1時的y坐標(biāo)，試完成加工所需數(shù)據(jù)，畫出曲線.x035791112131415Y101.82.22.73.03.12.92.52.01.6Y201.21.72.02.02.01.81.21.01.6第三頁，共六十五頁，2022年，8月28日山體地貌要在某山區(qū)方圓大約27平方公里范圍內(nèi)修建一條公路，從山腳出發(fā)經(jīng)過一個居民區(qū)，再到達(dá)一個礦區(qū)。橫向縱向分別每隔400米測量一次，得到一些地點的高程：試做出該山區(qū)的地貌圖.第四頁，共六十五頁，2022年，8月28日船在該海域會擱淺嗎？---作業(yè)

在某海域測得一些點(x,y)處的水深z由下表給出，船的吃水深度為5英尺，在矩形區(qū)域（75，200）*（-50，150）里的哪些地方船要避免進(jìn)入.第五頁，共六十五頁，2022年，8月28日

水深和流速的問題

在水文數(shù)據(jù)測量中，不同水深的流速是不同的.水文數(shù)據(jù)的測量時天天進(jìn)行的，為了減少測量的工作，希望得到確定的水深和水流之間的關(guān)系.為此測量了一系列不同水深和流速值.下表給出了對某河流的測量數(shù)據(jù)，其中水深和流速根據(jù)適當(dāng)?shù)膯挝贿M(jìn)行了規(guī)范化，共10個值.第六頁，共六十五頁，2022年，8月28日美國人口問題據(jù)美國人口普查局?jǐn)?shù)據(jù)：從1790每隔10年至2000年的總?cè)丝冢▎挝唬喊偃f）如下示

t=1790:10:2000;p=[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];

預(yù)測2001，2002年的美國人口數(shù)？并與調(diào)查數(shù)據(jù)285.318，288.369比較，選擇擬合較好的模型。第七頁，共六十五頁，2022年，8月28日農(nóng)作物施肥效果分析1992年A題

在農(nóng)業(yè)生產(chǎn)試驗研究中，對某地區(qū)土豆的產(chǎn)量與化肥的關(guān)系做了一實驗，得到了氮肥、磷肥的施肥量與土豆產(chǎn)量的對應(yīng)關(guān)系如下表：

1.根據(jù)上表數(shù)據(jù)分別給出土豆產(chǎn)量與氮、磷肥的關(guān)系式。

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第八頁，共六十五頁，2022年，8月28日

配藥方案---作業(yè)

一種新藥用于臨床之前,必須設(shè)計給藥方案.在快速靜脈注射的給藥方式下,所謂給藥方案是指,每次注射劑量多大,間隔時間多長.

藥物進(jìn)入機(jī)體后隨血液輸送到全身,在這個過程中不斷地被吸收,分布,代謝,最終排出體外.藥物在血液中的濃度,即單位體積血液中的藥物含量,稱血藥濃度.在最簡單的一室模型中,將整個機(jī)體看作一個房室,稱中心室,室內(nèi)的血藥濃度是均勻的.快速靜脈注射后,濃度立即上升;然后逐漸下降.當(dāng)濃度太低時,達(dá)不到預(yù)期的治療效果;血藥濃度太高,又可能導(dǎo)致藥物中毒或副作用太強(qiáng).臨床上,每種藥物有一個最小有效濃度c1和一個最大治療濃度c2.設(shè)計給藥方案時,要使血藥濃度保持在c1-c2之間.設(shè)本題所研究藥物的最小有效濃度c1=10,最大治療濃度c2=25().第九頁，共六十五頁，2022年，8月28日

顯然,要設(shè)計給藥方案,必須知道給藥后血藥濃度隨時間變化的規(guī)律.為此,從實驗和理論兩方面著手.在實驗方面,對某人用快速靜脈注射方式一次注入該藥物300mg后,在一定時刻t(小時)采集血樣,測得血藥濃度c.如表:血藥濃度c(t)的測試數(shù)據(jù)

t0.250.511.523468c19.2118.1515.3614.1012.899.327.455.243.01問題：1.在快速靜脈注射的給藥方式下，研究血藥濃度（單位體積血液中的藥物含量）的變化規(guī)律；2.給定藥物的最小有效濃度和最大治療濃度，設(shè)計給藥方案：每次注射劑量多大；間隔時間多長？

配藥方案第十頁，共六十五頁，2022年，8月28日二、問題的解決（1）插值法；（2）曲線擬合法．

1、問題的抽象xx1x2…xmyy1y2…ym構(gòu)造一個簡單易于計算的近似函數(shù)

p(x)f(x)

（精確函數(shù)）。2、構(gòu)造近似函數(shù)，p(x)的方法有兩種：在實驗中經(jīng)常給出一組離散點，第十一頁，共六十五頁，2022年，8月28日插值法定義：當(dāng)精確函數(shù)y=f(x)非常復(fù)雜或未知時，在一系列節(jié)點

x0…xn

處測得函數(shù)值y0

=f(x0),…，yn

=f(xn)，由此構(gòu)造一個簡單易算的近似函數(shù)p(x)

f(x)，滿足條件p(xi)=f(xi)(i=0,…n)，（插值條件）這里的p(x)

稱為f(x)的插值函數(shù)；構(gòu)造插值函數(shù)的方法為插值法。第十二頁，共六十五頁，2022年，8月28日曲線擬合

但是不要求使p(xi)=yi,而只要p(xi)yi

總體上盡可能小。這種構(gòu)造近似函數(shù)p(x)

的方法稱為曲線擬合法，p(x)稱為擬合函數(shù)。定義：當(dāng)精確函數(shù)y=f(x)非常復(fù)雜或未知時，在一系列節(jié)點x0…xn

處,測得函數(shù)值y0

,…，yn

，由此構(gòu)造一個簡單易算的近似函數(shù)p(x)

f(x)，第十三頁，共六十五頁，2022年，8月28日插值與擬合的相同點都需要根據(jù)已知數(shù)據(jù)構(gòu)造函數(shù)?？墒褂玫玫胶瘮?shù)計算未知點的函數(shù)值。xx1x2…xmyy1y2…ym求一個簡單易算的近似函數(shù)

p(x)f(x)

。第十四頁，共六十五頁，2022年，8月28日插值與擬合的不同點插值:過節(jié)點;；擬合:不過點,整體近似；第十五頁，共六十五頁，2022年，8月28日插值法拉格朗日插值牛頓插值三次埃爾米特插值法分段線性插值分段三次埃爾米特插值法三次樣條插值第十六頁，共六十五頁，2022年，8月28日1、拉格朗日插值公式（１）定義對給定的n+1個節(jié)點x0,x1,x2,…,xn及對應(yīng)的函數(shù)值y0,

y1,y2,…,yn，構(gòu)造一個n次插值多項式：即為拉格朗日插值公式，其中插值基函數(shù)第十七頁，共六十五頁，2022年，8月28日拉格朗日插值的matlab實現(xiàn)functiony=lagrange(x0,y0,x)%x0插值節(jié)點，y0插值節(jié)點處的函數(shù)值，x要計算函數(shù)值的點；n=length(x0);%計算x0的長度m=length(x);%計算x的長度fori=1:ms=0;z=x(i);

fork=1:np=1.0;forj=1:n ifj~=kp=p*(z-x0(j))/(x0(k)-x0(j));%計算插值基函數(shù)

endends=p*y0(k)+s;endy(i)=s;%計算在x(i)處的函數(shù)值（拉格朗日）end第十八頁，共六十五頁，2022年，8月28日2、牛頓插值法牛頓插值公式：Nn(x)=f(x0)+f[x0,x1](x-x0)+…+f[x0,x1,x2,…,xn](x-x0)(x-x1)…(x-xn)其中：

f[x0,x1]一階差商

f[x0,x1,x2,…,xn]n階差商注：牛頓插值法與拉格朗日插值法，同一個多項式，不同的表達(dá)方式，但是計算量不一樣，牛頓插值法的計算量小。第十九頁，共六十五頁，2022年，8月28日龍格現(xiàn)象Runge在上個世紀(jì)初發(fā)現(xiàn)：在[-5,5]上用n+1個等距節(jié)點作n次插值多項式Pn(x)，當(dāng)在n→∞時，插值多項式Pn(x)在區(qū)間中部趨于f(x)=1/(1+x2),但對于3.63≤∣x∣≤1的x，Pn(x)嚴(yán)重發(fā)散。用圖形分析問題。

第二十頁，共六十五頁，2022年，8月28日forn=10:2:20%從10等份到20等份x0=[-5:10/n:5];%插值節(jié)點y0=1./(1+x0.^2);%插值節(jié)點處的精確函數(shù)值x=[-5:0.1:5];%要進(jìn)行計算函數(shù)值的點y=lagrange(x0,y0,x);%調(diào)用函數(shù)計算x點的函數(shù)值plot(x0,y0,‘*’,x,1./(1+x.^2),‘r’,x,y)%繪制圖形pause%等待，按任意鍵end第二十一頁，共六十五頁，2022年，8月28日3、分段低次插值法（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)=第二十二頁，共六十五頁，2022年，8月28日（2）分段三次埃爾米特插值法定義：已知n+1個不同節(jié)點x0,x1,…,xn,構(gòu)造分段多項式I(x)，使之滿足：l

I(x)在[a，b]上二階連續(xù)導(dǎo)數(shù)；l

I(xk)=yk，I’(xk)=y’k，；l

I(x)在[xi,xi+1]上是三次次多項式。第二十三頁，共六十五頁，2022年，8月28日4、三次樣條插值法

對于給定n+1個不同節(jié)點x0,x1,…,xn及函數(shù)值y0,y1,…,yn，其中a=x0<x1<…<xn=b，構(gòu)造三次樣條插值函數(shù)S(x)。

S(x)稱為三次樣條函數(shù)時需滿足：l

S(x)在[a,b]上二階導(dǎo)數(shù)連續(xù)；l

S(xk)=yk(k=0,1,…,n)；l

每個子區(qū)間[xk,xk+1]上S(x)是三次多項式(k=0,1,…,n)。第二十四頁，共六十五頁，2022年，8月28日插值法的matlab實現(xiàn)—一維插值

命令：interp1(x0,y0,x，’method’)

其中：x0：插值節(jié)點；

y0：插值節(jié)點處的函數(shù)值；

x：要計算函數(shù)值的點；

method：

linear

：分段線性插值；

cubic：分段三次埃爾米特插值；

spline：三次樣條插值。第二十五頁，共六十五頁，2022年，8月28日插值法的應(yīng)用一水庫上游河段降暴雨，根據(jù)預(yù)報測算上游流入水庫的流量為Q(t)(102立方米/秒)

：

t(時)81216243044485660Q（t）3654789210135251613

通過這個預(yù)報值，分別用不同的數(shù)值方法插值法來估計14和20時上游流入水庫的流量。第二十六頁，共六十五頁，2022年，8月28日二維插值的MATLAB實現(xiàn)

在MATLAB中，二維插值命令常用的有兩個，

1、一個是網(wǎng)格節(jié)點插值：

z=interp2(x0,y0,z0,x,y,’method’)

其中，

z：被插值點處的函數(shù)值；

x0,y0,z0：插值節(jié)點，x0,y0為向量，z0是矩陣，其列數(shù)等于x0的長度，行數(shù)等于y0的長度；

x,y：要計算函數(shù)值的點；

interp1(x0,y0,x，’method’)第二十七頁，共六十五頁，2022年，8月28日山體地貌要在某山區(qū)方圓大約27平方公里范圍內(nèi)修建一條公路，從山腳出發(fā)經(jīng)過一個居民區(qū)，再到達(dá)一個礦區(qū)。橫向縱向分別每隔400米測量一次，得到一些地點的高程：試做出該山區(qū)的地貌圖，并對幾種插值法進(jìn)行比較.第二十八頁，共六十五頁，2022年，8月28日程序設(shè)計：clearx0=[1200:400:4000];y0=[1200:400:3600];z0=[11301250128012301040900500700;13201450142014001300700900850;139015001500140090011001060950;15001200110013501450120011501010;15001200110015501600155011801070;15001550160015501600160016001550;1480150015501510143013001200980];xi=1200:10:4000;%加密數(shù)據(jù)點yi=1200:10:3600;zil=interp2(x0,y0,z0,xi',yi,'linear');%線性插值zic=interp2(x0,y0,z0,xi',yi,'cubic');%三次插值zis=interp2(x0,y0,z0,xi',yi,'spline');%樣條插值subplot(2,2,1)mesh(x0,y0,z0)subplot(2,2,2)mesh(xi,yi,zil)subplot(2,2,3)mesh(xi,yi,zic)subplot(2,2,4)mesh(xi,yi,zis)第二十九頁，共六十五頁，2022年，8月28日第三十頁，共六十五頁，2022年，8月28日二維插值的MATLAB實現(xiàn)2、另一個是離散數(shù)據(jù)節(jié)點的插值命令：

z=griddata(x0,y0,z0,x,y,’method’)

其中，

z：被插值點處的函數(shù)值；

x0,y0,z0：插值節(jié)點，x0,y0,z0均為向量；

x,y：被插值點；

method：插值方法，包括：

'linear'——線性插值；

'cubic'——三次插值；第三十一頁，共六十五頁，2022年，8月28日船在該海域會擱淺嗎？

在某海域測得一些點(x,y)處的水深z由下表給出，船的吃水深度為5英尺，在矩形區(qū)域（75，200）*（-50，150）里的哪些地方船要避免進(jìn)入.第三十二頁，共六十五頁，2022年，8月28日解決問題的步驟：1.作出測量點的分布圖；2.求出矩形區(qū)域（75，200）*（-50，150）的細(xì)分網(wǎng)格節(jié)點之橫、縱坐標(biāo)向量；3.利用MATLAB中的散點插值函數(shù)求網(wǎng)格節(jié)點的水深；4.作出海底曲面圖形和等高線圖；5.作出水深小于5的海域范圍.第三十三頁，共六十五頁，2022年，8月28日程序clearx0=[129140103.588185.5195105157.5107.57781162162117.5];y0=[7.5141.52314722.5137.585.5-6.5-81356.5-66.584-33.5];z0=[-4-8-6-8-6-8-8-9-9-8-8-9-4-9];subplot(2,2,1)plot(x0,y0,'+');%作出測量點的分布圖；x=75:1:200;%加密y=-50:1:150;[x,y]=meshgrid(x,y);z=griddata(x0,y0,z0,x,y,'cubic');subplot(2,2,2)mesh(x,y,z),%用插值方法求出網(wǎng)格節(jié)點處的z坐標(biāo)矩陣,繪制出三維圖形subplot(2,2,3)meshc(x,y,z),%繪制等高線subplot(2,2,4)contour(x,y,z,[-5-5]);%水深5英尺處海底曲面的等高線gridon第三十四頁，共六十五頁，2022年，8月28日第三十五頁，共六十五頁，2022年，8月28日擬合的標(biāo)準(zhǔn)（1）用各點誤差絕對值的和表示（2）用各點誤差按絕對值的最大值表示（3）用各點誤差的平方和表示第三十六頁，共六十五頁，2022年，8月28日最小二乘擬合式中R2稱為均方誤差。由于計算均方誤差的最小值的原則容易實現(xiàn)而被廣泛采用。按均方誤差達(dá)到極小構(gòu)造擬合曲線的方法稱為最小二乘法。第三十七頁，共六十五頁，2022年，8月28日++++++++++++++++++++++++++++++p=a1+a2xp=a1+a2x+a3x2p=a1+a2x+a3x2p=a1+a2/xp=aebxp=ae-bx將數(shù)據(jù)(xi,yi)i=1,…,n

作圖，通過直觀判斷確定p(x)：第三十八頁，共六十五頁，2022年，8月28日MATLAB---曲線擬合工具箱Matlab有一個功能強(qiáng)大的曲線擬合工具箱（CurveFittingToolbox

）cftool，使用方便，能實現(xiàn)多種類型的線性、非線性曲線擬合。調(diào)用：cftool界面如下所示第三十九頁，共六十五頁，2022年，8月28日第四十頁，共六十五頁，2022年，8月28日“Data”按鈕數(shù)據(jù)的選取點擊“Data”按鈕，彈出“Data”窗口；利用Xdata和Ydata的下拉菜單讀入數(shù)據(jù)x,y，可修改數(shù)據(jù)集名“Datasetname”，然后點擊“Createdataset”按鈕，退出“Data”窗口，返回工具箱界面，這時會自動畫出數(shù)據(jù)集的曲線圖；

第四十一頁，共六十五頁，2022年，8月28日第四十二頁，共六十五頁，2022年，8月28日“Fitting”按鈕曲線擬合點擊“Fitting”按鈕，彈出“Fitting”窗口；點擊“Newfit”按鈕，可修改擬合項目名稱“Fitname”，通過“Dataset”下拉菜單選擇數(shù)據(jù)集，然后通過下拉菜單“Typeoffit”選擇擬合曲線的類型。第四十三頁，共六十五頁，2022年，8月28日SSEThesumofsquaresduetoerror.Thisstatisticmeasuresthedeviationoftheresponsesfromthefittedvaluesoftheresponses.Avaluecloserto0indicatesabetterfit.偏差平方和，越接近0越好第四十四頁，共六十五頁，2022年，8月28日R-square

Thecoefficientofmultipledetermination.Thisstatisticmeasureshowsuccessfulthefitisinexplainingthevariationofthedata.Avaluecloserto1indicatesabetterfit.復(fù)相關(guān)系數(shù)平方(決定系數(shù))，越接近1越好第四十五頁，共六十五頁，2022年，8月28日AdjustedR-squareThedegreeoffreedomadjustedR-square.Avaluecloserto1indicatesabetterfit.Itisgenerallythebestindicatorofthefitqualitywhenyouaddadditionalcoefficientstoyourmodel.修正的復(fù)相關(guān)系數(shù)平方，越接近1越好第四十六頁，共六十五頁，2022年，8月28日AdjustedR-square下列公式中的m為擬合函數(shù)中待估參數(shù)個數(shù)，如：對一元一次多項式擬合，f(x)=a+bx，此時m=2，n為數(shù)據(jù)點個數(shù)。該修正類似修正的樣本方差使其為總體方差的無偏估計。第四十七頁，共六十五頁，2022年，8月28日RMSETherootmeansquarederror.Avaluecloserto0indicatesabetterfit.偏差平方的均值的算術(shù)平方根，越接近0越好第四十八頁，共六十五頁，2022年，8月28日曲線擬合好壞如何評價首要指標(biāo)是目標(biāo)函數(shù)誤差最?。〝M合度最大）；其次是應(yīng)考慮關(guān)鍵點的吻合，這些關(guān)鍵點包括：初始點（有時是原點）、拐點、峰值點、極值點、中間點、漸近點、終值點等，在這些關(guān)鍵點上，數(shù)據(jù)觀察值點與函數(shù)值點應(yīng)盡可能一致；再次是擬合的模型應(yīng)盡可能簡單（模型的形式簡單，參數(shù)數(shù)少）。第四十九頁，共六十五頁，2022年，8月28日實踐中如何選擇模型？在數(shù)據(jù)擬合實踐中，理性模型畢竟是少數(shù)，大多數(shù)的情形是根據(jù)數(shù)據(jù)的趨勢尋找合適的模型，有時好幾個模型對數(shù)據(jù)都有較好的擬合，但通過對關(guān)鍵點的比較總會找到一種最合適的模型。在選擇不同的模型時，合理性和可解釋性是首要考慮的因素。第五十頁，共六十五頁，2022年，8月28日美國人口問題據(jù)美國人口普查局?jǐn)?shù)據(jù)：從1790每隔10年至2000年的總?cè)丝冢▎挝唬喊偃f）如下示

t=1790:10:2000;p=[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];第五十一頁，共六十五頁，2022年，8月28日美國人口問題t=1790:10:2000;p=[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];預(yù)測2001，2002年的美國人口數(shù)？并與調(diào)查數(shù)據(jù)285.318，288.369比較，選擇擬合較好的模型。第五十二頁，共六十五頁，2022年，8月28日Matlab求解在命令窗口輸入命令cftool回車，得擬合的圖形用戶界面第五十三頁，共六十五頁，2022年，8月28日第五十四頁，共六十五頁，2022年，8月28日第五十五頁，共六十五頁，2022年，8月28日結(jié)果分析LinearmodelPoly2:f(x)=p1*x^2+p2*x+p3Coefficients(with95%confidencebounds):p1=0.006757(0.006369,0.007144)p2=-24.32(-25.78,-22.85)p3=2.188e+004(2.049e+004,2.327e+004)Goodnessoffit:SSE: