數(shù)學建模與實驗嚴可頌主講教學項目四數(shù)值分析法模型4.1插值法建模

拉格朗日插值分段線性插值三次樣條插值一、插值的定義二、插值的方法三、用Matlab解插值問題已知n+1個節(jié)點其中互不相同，不妨設(shè)求任一插值點處的插值節(jié)點可視為由產(chǎn)生，表達式復(fù)雜，或無封閉形式，或未知。4.1.1插值法的定義

構(gòu)造一個(相對簡單的)函數(shù)通過全部節(jié)點,即再用計算插值，即4.1.2插值的方法

已知函數(shù)f(x)在n+1個點x0,x1,…,xn處的函數(shù)值為y0,y1,…,yn

。求一n次多項式函數(shù)Pn(x)，使其滿足：拉格朗日(Lagrange)插值稱為拉格朗日插值基函數(shù)。解決此問題的拉格朗日插值多項式公式如下其中Li(x)為n次多項式：特殊情況兩點一次(線性)插值多項式:三點二次(拋物)插值多項式:采用拉格朗日多項式插值：選取不同插值節(jié)點個數(shù)n+1，其中n為插值多項式的次數(shù)，當n取10時，繪出插值結(jié)果圖形.例分段線性插值計算量與n無關(guān);n越大，誤差越小.xjxj-1xj+1x0xnxoy比分段線線性插值值更光滑滑。xyxi-1xiab在數(shù)學上上，光滑滑程度的的定量描描述是：：函數(shù)(曲線)的k階導(dǎo)數(shù)存存在且連連續(xù)，則則稱該曲線具具有k階光滑性性。光滑性的的階次越越高，則則越光滑滑。是否否存在較較低次的的分段多多項式達達到較高高階光滑滑性的方方法？三三次樣條條插值就就是一個個很好的的例子。。三次樣條條插值三次樣條條插值g(x)為被插值函函數(shù)。用MATLAB作插值計計算一維插值值函數(shù)：：yi=interp1(x，y，xi，'method')插值方法被插值點插值節(jié)點xi處的插值結(jié)果‘nearest’：最鄰近近插值‘linear’’：線性插值值；‘spline’’：三次次樣條插插值；‘cubic’：立方方插值。。缺省時：：分分段線性性插值。。注意：所所有的插插值方法法都要求求x是單調(diào)的的，并且且xi不能夠超超過x的范圍。。例：在1-12的11小時內(nèi)，，每隔1小時測量量一次溫溫度，測測得的溫溫度依次次為：5，8，9，15，25，29，31，30，22，25，27，24。試估計計每隔1/10小時的溫溫度值。。hours=1:12;h=1:0.1:12;t=interp1(hours,temps,h,'spline');(直接輸出出數(shù)據(jù)將將是很多多的)plot(hours,temps,'+',h,t,hours,temps,'r:')%作圖xlabel('Hour'),ylabel('DegreesCelsius’)xy實驗習題題已知飛機機下輪廓廓線上數(shù)數(shù)據(jù)如下下，求x每改變0.1時的y值。機翼下輪廓線4.2擬合法簡介

2.擬合的基基本原理理1.擬合問題題引例擬合問問題題引例例1溫度t(0C)20.532.751.073.095.7電阻R()7658268739421032已知熱敏電阻數(shù)據(jù)：求600C時的電阻阻R。設(shè)R=at+ba,b為待定系系數(shù)擬合問問題題引例例2

t(h)0.250.511.523468c(g/ml)19.2118.1515.3614.1012.899.327.455.243.01已知一室模型快速靜脈注射下的血藥濃度數(shù)據(jù)(t=0注射300mg)求血藥濃濃度隨時時間的變變化規(guī)律律c(t).作半對數(shù)數(shù)坐標系系(semilogy)下的圖形形曲線擬擬合合問題題的的提法法已知一組組（二維維）數(shù)據(jù)據(jù)，即平平面上n個點（xi,yi)i=1,…n,尋求一個個函數(shù)（（曲線））y=f(x),使f(x)在某種準準則下與與所有數(shù)數(shù)據(jù)點最最為接近近，即曲曲線擬合合得最好好。+++++++++xyy=f(x)(xi,yi)ii為點（xi,yi)與曲線y=f(x)的距離擬合與插值的的關(guān)系函數(shù)插值與曲曲線擬合都是是要根據(jù)一組組數(shù)據(jù)構(gòu)造一一個函數(shù)作為為近似，由于于近似的要求求不同，二者者的數(shù)學方法法上是完全不不同的。實例：下面數(shù)據(jù)是某某次實驗所得得，希望得到到X和f之間的關(guān)系？？問題：給定一批數(shù)據(jù)據(jù)點，需確定定滿足特定要要求的曲線或或曲面解決方案：若不要求曲線線（面）通過過所有數(shù)據(jù)點點，而是要求求它反映對象象整體的變化化趨勢，這就就是數(shù)據(jù)擬合，又稱曲線擬擬合或曲面擬擬合。若要求所求曲曲線（面）通通過所給所有有數(shù)據(jù)點，就就是插值問題；最臨近插值、、線性插值、、樣條插值與與曲線擬合結(jié)結(jié)果：曲線擬合問題題最常用的解解法——線性最小二乘乘法的基本思思路第一步:先選定一組函函數(shù)r1(x),r2(x),……rm(x),m<n,令f(x)=a1r1(x)+a2r2(x)+……+amrm(x)（1）其中a1,a2,…am為待定系數(shù)。。第二步:確定a1,a2,…am的準則（最小小二乘準則））：使n個點（xi,yi)與曲線y=f(x)的距離i的平方和最小小。記

問題歸結(jié)為，，求a1,a2,…am使J(a1,a2,…am)最小。線性最小二乘乘法的求解：：預(yù)備知識超定方程組：方程個數(shù)大大于未知量個個數(shù)的方程組組即Ra=y其中超定方程一般般是不存在解解的矛盾方程程組。如果有向量a使得達到最小，則稱a為上述超定方程的最小二乘解。線性最小二乘乘法的求解定理：當RTR可逆時，超定定方程組（3）存在最小二二乘解，且即即為方程組RTRa=RTy的解：a=(RTR)-1RTy所以，曲線擬擬合的最小二二乘法要解決決的問題，實實際上就是求求以下超定方方程組的最小小二乘解的問問題。其中Ra=y（3）線性最小二乘乘擬合f(x)=a1r1(x)+……+amrm(x)中函數(shù){r1(x),……rm(x)}的選取1.通過機理分析析建立數(shù)學模模型來確定f(x)；++++++++++++++++++++++++++++++f=a1+a2xf=a1+a2x+a3x2f=a1+a2x+a3x2f=a1+a2/xf=aebxf=ae-bx2.將數(shù)據(jù)(xi,yi)i=1,…n作圖，通過直直觀判斷確定定f(x)：用MATLAB解擬合問題1、線性最小二二乘擬合2、非線性最小小二乘擬合用MATLAB作線性最小二二乘擬合1.作多項式f(x)=a1xm+…+amx+am+1擬合,可利用已有程程序:a=polyfit(x,y,m)2.對超定方程組可得最小二乘意義下的解。，用3.多項式在x處的值y可用以下命令令計算：y=polyval（a，x）輸入同長度的數(shù)組X，Y擬合多項式次數(shù)即要求出二次多項式:中的使得:例對下面一組數(shù)據(jù)作二次多項式擬合1）輸入以下命命令：x=0:0.1:1;y=[-0.4471.9783.286.167.087.347.669.569.489.3011.2];R=[(x.^2)'x'ones(11,1)]；A=R\y'解法1．用解超定方程程的方法2）計算結(jié)果：Ａ1）輸入以下命命令：x=0:0.1:1;y=[-0.4471.9783.286.167.087.347.669.569.489.3011.2];A=polyfit(x,y,2)z=polyval(A,x);plot(x,y,'k+',x,z,'r')%作出數(shù)據(jù)點和和擬合曲線的的圖形2）計算結(jié)果：：Ａ解法2．用多項式擬合合的命令1.lsqcurvefit已知數(shù)據(jù)點：xdata=（xdata1，xdata2，…，xdatan），ydata=（ydata1，ydata2，…，ydatan）用MATLAB作非線性最小小二乘擬合Matlab的提供了兩個個求非線性最最小二乘擬合合的函數(shù)：lsqcurvefit和lsqnonlin。兩個命令都都要先建立M-文件fun.m，在其中定義義函數(shù)f(x)，但兩者定義義f(x)的方式是不同同的,可參考例題.

lsqcurvefit用以求含參量x（向量）的向量值函數(shù)F(x,xdata)=（F（x，xdata1），…，F(xiàn)（x，xdatan））T中的參變量x(向量),使得輸入格式為:(1)x=lsqcurvefit(‘fun’’,x0,xdata,ydata);(2)x=lsqcurvefit(‘fun’’,x0,xdata,ydata,options);(3)x=lsqcurvefit(‘fun’,x0,xdata,ydata,options,’grad’);(4)[x,options]=lsqcurvefit(‘fun’,x0,xdata,ydata,…);(5)[x,options,funval]=lsqcurvefit(‘‘fun’,x0,xdata,ydata,…);(6)[x,options,funval,Jacob]=lsqcurvefit(‘fun’,x0,xdata,ydata,…);fun是一個事先建立的定義函數(shù)F(x,xdata)

的M-文件,自變量為x和xdata說明：x=lsqcurvefit(‘fun’,x0,xdata,ydata,options);迭代初值已知數(shù)據(jù)點選項見無約束優(yōu)化

lsqnonlin用以求含參量x（向量）的向量值函數(shù)

f(x)=(f1(x),f2(x),…,fn(x))T

中的參量x，使得

最小。其中fi（x）=f（x，xdatai，ydatai）

=F(x,xdatai)-ydatai

2.lsqnonlin已知數(shù)據(jù)點：：xdata=（xdata1，xdata2，…，xdatan）ydata=（ydata1，ydata2，…，ydatan）輸入格式為：：1）x=lsqnonlin（‘fun’，x0）；2）x=lsqnonlin（‘fun’，x0，options）；3）x=lsqnonlin（‘fun’，x0，options，‘grad’）；4）[x，options]=lsqnonlin（‘fun’，x0，…）；5）[x，options，funval]=lsqnonlin（‘fun’，x0，…）；說明：x=lsqnonlin（‘fun’，x0，options）；fun是一個事先建立的定義函數(shù)f(x)的M-文件，自變量為x迭代初值選項見無約束優(yōu)化

例2用下面一組數(shù)據(jù)擬合中的參數(shù)a，b，k該問題即解最最優(yōu)化問題：：1）編寫M-文件curvefun1.mfunctionf=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: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)解法1.用命令令lsqcurvefit3）運算算結(jié)果果為：4）結(jié)論論：a=0.0063,b=-0.0034,k=0.2542

解法2

用命令lsqnonlin

f(x)=F(x,tdata,ctada)=x=（a，b，k）1）編寫M-文件curvefun2.mfunctionf=curvefun2(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)3）運算算結(jié)果果為可以看看出,兩個命命令的的計算算結(jié)果果是相相同的的.4）結(jié)論論：即擬合合得a=0.0063b=-0.0034k=0.2542MATLAB解應(yīng)用用問題題實例例1、電阻阻問題題2、給藥藥方案案問題題*3、水塔塔流量量估計計問題題電阻問問題溫度t(0C)20.532.751.073.095.7電阻R()7658268739421032例.由數(shù)據(jù)擬合R=a1t+a2方法1.用命令polyfit(x,y,m)得到a1=3.3940,a2=702.4918方法2.直接用結(jié)果相相同。。一室模模型：將整整個機機體看看作一一個房房室，，稱中心室室，室內(nèi)內(nèi)血藥藥濃度度是均均勻的的。快快速靜靜脈注注射后后，濃濃度立立即上上升；；然后后迅速速下降降。當當濃度度太低低時，，達不不到預(yù)預(yù)期的的治療療效果果；當當濃度度太高高，又又可能能導(dǎo)致致藥物物中毒毒或副副作用用太強強。臨臨床上上，每每種藥藥物有有一個個最小小有效效濃度度c1和一個個最大大有效效濃度度c2。設(shè)計計給藥藥方案案時，，要使使血藥藥濃度度保保持在在c1~c2之間。。本題題設(shè)c1=10，c2=25(ug/ml).擬合問題實例2給藥方案——一種新新藥用用于臨臨床之之前，，必須須設(shè)計計給藥藥方案案.藥物進進入機機體后后血液液輸送送到全全身，，在這這個過過程中中不斷斷地被被吸收收、分分布、、代謝謝，最最終排排出體體外，，藥物物在血血液中中的濃濃度，，即單單位體體積血血液中中的藥藥物含含量，，稱為為血藥濃濃度。。

在實驗方面,對某人用快速靜脈注射方式一次注入該藥物300mg后,在一定時刻t(小時)采集血藥,測得血藥濃度c(ug/ml)如下表:

t(h)0.250.511.523468c(g/ml)19.2118.1515.3614.1012.899.327.455.243.01要設(shè)計計給藥藥方案案,必須知知道給給藥后后血藥藥濃度度隨時時間變變化的的規(guī)律律。從從實驗驗和理理論兩兩方面面著手手：給藥方方案1.在快速速靜脈脈注射射的給給藥方方式下下，研研究血血藥濃濃度（（單位位體積積血液液中的的藥物物含量量）的的變化化規(guī)律律。tc2cc10問題題2.給定定藥藥物物的的最最小小有有效效濃濃度度和和最最大大治治療療濃濃度度，，設(shè)設(shè)計計給給藥藥方方案案：：每每次次注注射射劑劑量量多多大大；；間間隔隔時時間間多多長長。。分析析理論論：：用用一一室室模模型型研研究究血血藥藥濃濃度度變變化化規(guī)規(guī)律律實驗驗：：對對血血藥藥濃濃度度數(shù)數(shù)據(jù)據(jù)作作擬擬合合，，符符合合負負指指數(shù)數(shù)變變化化規(guī)規(guī)律律3.血液液容容積積v,t=0注射射劑劑量量d,血藥藥濃濃度度立立即即為為d/v.2.藥物物排排除除速速率率與與血血藥藥濃濃度度成成正正比比，，比比例例系系數(shù)數(shù)k(>0)模型型假假設(shè)設(shè)1.機體體看看作作一一個個房房室室，，室室內(nèi)內(nèi)血血藥藥濃濃度度均均勻勻———一室室模模型型模型型建建立立在此此，，d=300mg，t及c（t）在在某某些些點點處處的的值值見見前前表表，，需需經(jīng)經(jīng)擬擬合合求求出出參參數(shù)數(shù)k、v用線線性性最最小小二二乘乘擬擬合合c(t)計算結(jié)果：d=300;t=[0.250.511.523468];c=[19.2118.1515.3614.1012.899.327.455.243.01];y=log(c);a=polyfit(t,y,1)k=-a(1)v=d/exp(a(2))程序：給藥藥方方案案設(shè)設(shè)計計cc2c10t設(shè)每每次次注注射射劑劑量量D,間隔隔時時間間血藥藥濃濃度度c(t)應(yīng)c1c(t)c2初次次劑劑量量D0應(yīng)加加大大給藥方案記為：2、1、計算結(jié)果：給藥方案：c1=10,c2=25k=0.2347v=15.02故可可制制定定給給藥藥方方案案：：即:首次次注注射射375mg，其余余每每次次注注射射225mg，注射射的的間間隔隔時時間間為為4小時時。。估計計水水塔塔的的流流量量2、解解題題思思路路3、算算法法設(shè)設(shè)計計與與編編程程1、問問題題某居居民民區(qū)區(qū)有有一一供供居居民民用用水水的的園園柱柱形形水水塔塔，，一一般般可可以以通通過過測測量量其其水水位位來來估估計計水水的的流流量量，，但但面面臨臨的的困困難難是是，，當當水水塔塔水水位位下下降降到到設(shè)設(shè)定定的的最最低低水水位位時時，，水水泵泵自自動動啟啟動動向向水水塔塔供供水水，，到到設(shè)設(shè)定定的的最最高高水水位位時時停停止止供供水水，，這這段段時時間間無無法法測測量量水水塔塔的的水水位位和和水水泵泵的的供供水水量量．．通通常常水水泵泵每每天天供供水水一一兩兩次次，，每每次次約約兩兩小小時時.水塔塔是是一一個個高高12.2米，，直直徑徑17.4米的的正正園園柱柱．．按按照照設(shè)設(shè)計計，，水水塔塔水水位位降降至至約約8.2米時時，，水水泵泵自自動動啟啟動動，，水水位位升升到到約約10.8米時時水水泵泵停停止止工工作作．．表1是某某一一天天的的水水位位測測量量記記錄錄，，試試估估計計任任何何時時刻刻（（包包括括水水泵泵正正供供水水時時））從從水水塔塔流流出出的的水水流流量量，，及及一一天天的的總總用用水水量量．．流量量估估計計的的解解題題思思路路擬合合水水位位~時間間函函數(shù)數(shù)確定定流流量量~時間間函函數(shù)數(shù)估計計一一天天總總用用水水量量擬合合水水位位~時間間函函數(shù)數(shù)測量量記記錄錄看看，，一一天天有有兩兩個個供供水水時時段段（（以以下下稱稱第第1供水水時時段段和和第第2供水水時時段段）），，和和3個水水泵泵不不工工作作時時段段（（以以下下稱稱第第1時段段t=0到t=8.97，第第2次時時段段t=10.95到t=20.84和第第3時段段t=23以后后））．．對對第第1、2時段段的的測測量量數(shù)數(shù)據(jù)據(jù)直直接接分分別別作作多多項項式式擬擬合合，，得得到到水水位位函函數(shù)數(shù)．．為為使使擬擬合合曲曲線線比比較較光光滑滑，，多多項項式式次次數(shù)數(shù)不不要要太太高高，，一一般般在在3~6．由由于于第第3時段段只只有有3個測測量量記記錄錄，，無無法法對對這這一一時時段段的的水水位位作作出出較較好好的的擬擬合合．．2、確定定流流量量~時間間函函數(shù)數(shù)對于于第第1、2時段段只只需需將將水水位位函函數(shù)數(shù)求求導(dǎo)導(dǎo)數(shù)數(shù)即即可可，，對對于于兩兩個個供供水水時時段段的的流流量量，，則則用用供供水水時時段段前前后后（（水水泵泵不不工工作作時時段段））的的流流量量擬擬合合得得到到，，并并且且將將擬擬合合得得到到的的第第2供水時段段流量外外推，將將第3時段流量量包含在在第2供水時段段內(nèi)．3、一天總用用水量的的估計總用水量量等于兩兩個水泵泵不工作作時段和和兩個供供水時段段用水量量之和，，它們都都可以由由流量對對時間的的積分得得到。算法設(shè)計計與編程程1、擬合第1、2時段的水水位，并并導(dǎo)出流流量2、擬合供水水時段的的流量3、估計一天天總用水水量4、流量及及總用水水量的檢檢驗1、擬合第1時段的水水位，并并導(dǎo)出流流量設(shè)t，h為已輸入入的時刻刻和水位位測量記記錄（水水泵啟動動的4個時刻不不輸入）），第1時段各時刻的的流量可可如下得得：1）c1=polyfit（t（1：10），h（1：10），3）；%用3次多項式式擬合第第1時段水位位，c1輸出3次多項式式的系數(shù)數(shù)2）a1=polyder（c1）；%a1輸出多項項式（系系數(shù)為c1）導(dǎo)數(shù)的的系數(shù)3）tp1=0：0.1：9；x1=-polyval（a1，tp1）；%x1輸出多項項式（系系數(shù)為a1）在tp1點的函數(shù)數(shù)值（取取負后邊邊為正值值），即即tp1時刻的流流量4）流量函數(shù)數(shù)為：2、擬合第2時段的水水位，并并導(dǎo)出流流量設(shè)t，h為已輸入入的時刻刻和水位位測量記記錄（水水泵啟動動的4個時刻不不輸入）），第2時段各時刻的的流量可可如下得得：1）c2=polyfit(t(10.9:21),h(10.9:21),3)；%用3次多項式式擬合第第2時段水位位，c2輸出3次多項式式的系數(shù)數(shù)2）a2=polyder(c2)；%a2輸出多項項式（系系數(shù)為c2）導(dǎo)數(shù)的的系數(shù)3）tp2=10.9:0.1:21;x2=-polyval(a2,tp2);%x2輸出多項項式（系系數(shù)為a2）在tp2點的函數(shù)數(shù)值（取取負后邊邊為正值值），即即tp2時刻的流流量4）流量函數(shù)數(shù)為：3、擬合供水水時段的的流量在第1供水時段段（t=9~11）之前（（即第1時段）和和之后（（即第2時段）各各取幾點點，其流流量已經(jīng)經(jīng)得到，，用它們們擬合第第1供水時段段的流量量．為使使流量函函數(shù)在t=9和t=11連續(xù)，我我們簡單單地只取取4個點，擬擬合3次多項式式（即曲曲線必過過這4個點），，實現(xiàn)如如下：xx1=-polyval（a1，[89]）；%取第1時段在t=8,9的流量xx2=-polyval（a2，[1112]）；%取第2時段在t=11,12的流量xx12=[xx1xx2]；c12=polyfit（[891112]，xx12，3）；%擬合3次多項式式tp12=9：0.1：11；x12=polyval（c12，tp12）；%x12輸出第1供水時段段各時刻的的流量擬合的流流量函數(shù)數(shù)為：在第2供水時段段之前取取t=20，20.8兩點的流流水量，，在該時時刻之后后（第3時段）僅僅有3個水位記記錄，我我們用差差分得到