樣條插值函數(shù)及應(yīng)用

樣條插值函數(shù)及應(yīng)用摘要樣條函數(shù)具有廣泛的應(yīng)用，是現(xiàn)代函數(shù)論的一個(gè)十分活潑的分支，是計(jì)算方法的主要根底和工具之一，由于生產(chǎn)和科學(xué)技術(shù)向前開展的推動(dòng)以及電子計(jì)算機(jī)廣泛應(yīng)用的需要，人們便更多地應(yīng)用這個(gè)工具，也更深刻的認(rèn)識(shí)了它的本質(zhì)。在實(shí)際問題中所遇到許多函數(shù)往往很復(fù)雜，有些甚至是很難找到解析表達(dá)式的。根據(jù)函數(shù)已有的數(shù)據(jù)來(lái)計(jì)算函數(shù)在一些新的點(diǎn)處的函數(shù)值，就是插值法所需要解決的問題。插值法是數(shù)值逼近的重要方法之一，它是根據(jù)給定的自變量值和函數(shù)值,求取未知函數(shù)的近似值。早在一千多年前,我國(guó)科學(xué)家就在研究歷法時(shí)就用到了線性插值和二次插值。而在實(shí)際問題中,有許多插值函數(shù)的曲線要求具有較高的光滑性,在整個(gè)曲線中,曲線不但不能有拐點(diǎn),而且曲率也不能有突變。因此，對(duì)于插值函數(shù)必須二次連續(xù)可微且不變號(hào),這就需要用到三次樣條插值。關(guān)鍵詞三次樣條函數(shù)；插值法目錄TOC\o"1-2"\h\u24884引言124033第一章三次樣條插值2188261.1樣條插值函數(shù)簡(jiǎn)介2105481.2三次樣條函數(shù)應(yīng)用32869第二章AMCM91A估計(jì)水塔水流量5264682.1理論分析及計(jì)算6319232.2運(yùn)用MATLAB軟件計(jì)算928131參考文獻(xiàn) 14引言樣條函數(shù)具有廣泛的應(yīng)用，是現(xiàn)代函數(shù)論的一個(gè)十分活潑的分支，是計(jì)算方法的主要根底和工具之一，由于生產(chǎn)和科學(xué)技術(shù)向前開展的推動(dòng)以及電子計(jì)算機(jī)廣泛應(yīng)用的需要，人們便更多地應(yīng)用這個(gè)工具，也更深刻的認(rèn)識(shí)了它的本質(zhì)。上世紀(jì)四十年代，在研究數(shù)據(jù)處理的問題中引出了樣條函數(shù)，例如，在1946年Schoenberg將樣條引入數(shù)學(xué),即所謂的樣條函數(shù)，直到五十年代，還多應(yīng)用于統(tǒng)計(jì)數(shù)據(jù)的處理方面，從六十年代起，在航空、造船、汽車等行業(yè)中，開始大量采用樣條函數(shù)。在我國(guó)，從六十年代末開始，從船體數(shù)學(xué)放樣到飛機(jī)外形設(shè)計(jì)，逐漸出現(xiàn)了一個(gè)使用樣，逐漸出現(xiàn)了一個(gè)使用樣條函數(shù)的熱潮，并推廣到數(shù)據(jù)處理的許多問題中。在實(shí)際生活中有許多計(jì)算問題對(duì)插值函數(shù)的光滑性有較高的要求，例如飛機(jī)機(jī)翼外形、發(fā)動(dòng)機(jī)進(jìn)、排氣口都要求有連續(xù)的二階導(dǎo)數(shù)，用三次樣條繪制的曲線不僅有很好的光滑度，而且當(dāng)節(jié)點(diǎn)逐漸加密時(shí)其函數(shù)值整體上能很好地逼近被插函數(shù)，相應(yīng)的導(dǎo)數(shù)值也收斂于被插函數(shù)的導(dǎo)數(shù)值，不會(huì)發(fā)生“龍格現(xiàn)象”?，F(xiàn)在國(guó)內(nèi)外學(xué)者對(duì)這方面的研究也越來(lái)越重視，根據(jù)我們的需要來(lái)解決不同的問題，而且函數(shù)的形式也在不斷地改良，長(zhǎng)期以來(lái)很多學(xué)者致力于樣條插值的研究，對(duì)三次樣條的研究已相當(dāng)成熟。第一章三次樣條插值1.1樣條插值函數(shù)簡(jiǎn)介在實(shí)際問題中所遇到許多函數(shù)往往很復(fù)雜，有些甚至是很難找到解析表達(dá)式的。有時(shí)通過(guò)實(shí)驗(yàn)或者數(shù)值計(jì)算所得到的也只是一些離散的點(diǎn)上的函數(shù)值，即。根據(jù)函數(shù)已有的數(shù)據(jù)來(lái)計(jì)算函數(shù)在一些新的點(diǎn)處的函數(shù)值，就是插值法所需要解決的問題。插值法的根本思想就是，首先根據(jù)已有的函數(shù)值來(lái)構(gòu)造一個(gè)簡(jiǎn)單的函數(shù)作為的近似表達(dá)式，然后用來(lái)計(jì)算新的點(diǎn)上的函數(shù)值作為，因?yàn)槎囗?xiàng)式具有各階導(dǎo)數(shù)，求值也比見方便。常用的有Lagrange插值、Newton插值、Hermite線性插值在分段點(diǎn)上僅連續(xù)而不可導(dǎo)，三次埃爾米特插值有連續(xù)的一階導(dǎo)數(shù)，這樣的光滑程度常不能滿足物理問題的需要，樣條函數(shù)可以同時(shí)解決這兩個(gè)問題,使插值函數(shù)既是低階分段函數(shù),又是光滑的函數(shù)，并且只需在區(qū)間端點(diǎn)提供某些導(dǎo)數(shù)信息。三次樣條函數(shù)定義：設(shè)在區(qū)間上取個(gè)節(jié)點(diǎn)，函數(shù)假設(shè)滿足：(1)(2)在區(qū)間上，(3)在區(qū)間上,是那么稱是函數(shù)在區(qū)間上由以上定義可以看出，雖然每個(gè)子區(qū)間上的多項(xiàng)式可以各不相同，但在相鄰子區(qū)間的連接處卻是光滑的。因此，樣條插值也稱為分段光滑插值。從定義知要求出，在每一個(gè)小區(qū)間上確定4個(gè)待定系數(shù)，共有n個(gè)小區(qū)間，故應(yīng)有4n個(gè)參數(shù)。根據(jù)在上二階導(dǎo)數(shù)連續(xù)，在節(jié)點(diǎn)出滿足連續(xù)性條件共有3n-3個(gè)條件，再加上滿足插值條件共有4n-2個(gè)條件，因此還需要2個(gè)條件才能確定。通常可在區(qū)間端點(diǎn)上各加一個(gè)條件〔稱為邊界條件〕，可根據(jù)實(shí)際問題的要求給定，通常有以下三種：〔1〕端點(diǎn)的一階導(dǎo)數(shù)值，即〔2〕倆端點(diǎn)的二階導(dǎo)數(shù)，即其特殊情況稱為自然邊界條件。〔3〕當(dāng)是以為周期的函數(shù)時(shí)，那么要求也是周期函數(shù)。這時(shí)邊界條件應(yīng)滿足而此時(shí)。這樣確定的樣條函數(shù)，稱為周期函數(shù)。1.2三次樣條函數(shù)應(yīng)用作函數(shù)在取間隔為0.1的點(diǎn)圖，用插值進(jìn)行實(shí)驗(yàn)。使用MATLAB軟件程序代碼如下：%產(chǎn)生原始數(shù)據(jù)x=0:0.1:1;y=(x.^2-3*x+7).*exp(-4*x).*sin(2*x);%作圖subplot(1,2,1);plot(x,y,x,y,'ro')%待求插值點(diǎn)xx=0:0.02:1;yy=interp1(x,y,xx,'spline');%作圖subplot(1,2,2)plot(x,y,'ro',xx,yy,'b')運(yùn)行截圖圖1.1運(yùn)行結(jié)果第二章AMCM91A估計(jì)水塔水流量美國(guó)某洲的各用水管理機(jī)構(gòu)要求各社區(qū)提供以每小時(shí)多少加侖計(jì)的用水率以及每天總的用水量,但許多社區(qū)并沒有測(cè)量水流入或流出當(dāng)?shù)厮乃康脑O(shè)備,他們只能代之以每小時(shí)測(cè)量水塔中的水位,精度在0.5%以內(nèi),更為重要的是,無(wú)論什么時(shí)候,只要水塔中的水位下降到某一最低水位L時(shí),水泵就啟動(dòng)向水塔重新充水至某一最高水位H,但也無(wú)法得到水泵的供水量的測(cè)量數(shù)據(jù)。因此，在水泵工作時(shí)，人們?nèi)菀捉⑺械乃慌c水泵工作時(shí)的用水量之間的關(guān)系。水泵每天向水塔充水一次或兩次，每次約兩小時(shí)。表1白某小鎮(zhèn)某天的水塔水位時(shí)間033166635106191393717921212402522328543水位317531103054299429472892285027972752時(shí)間322843593539332394354331846636499535393657254水位2697水泵工作水泵工作355034453350326031673087時(shí)間605746455468535718547502179154826498596889953水位30122927284227672697水泵工作水泵工作34753397時(shí)間93270水位3340單位：時(shí)間/秒；水位/0.01英尺試估計(jì)在任何時(shí)刻，甚至包括水泵正在工作期間內(nèi)，水從水塔流出的流量，并估計(jì)一天的總用水量，表1中給出了某個(gè)真實(shí)小鎮(zhèn)某一天的真實(shí)數(shù)據(jù)。表1中給出了從第一次測(cè)量開始的以秒為單位的時(shí)刻，以及該時(shí)刻的高度單位為百分之一英尺的水塔中水位的測(cè)量值，例如，3316秒后，水塔中的水位到達(dá)31.10英尺。水塔是一個(gè)垂直圓形柱體，高為40英尺，直徑57英尺，通常當(dāng)水塔的水位降至27.00英尺時(shí)水泵開始向水塔充水，而當(dāng)水塔的水位升至35.50英尺時(shí)水泵停止工作。2.1理論分析及計(jì)算水塔充水時(shí)間確實(shí)定(1)第一次充水時(shí)間確實(shí)定當(dāng)時(shí)間t=32284秒時(shí)，水位26.97英尺,約低于最低水位27英尺，因此可作為第一次開始充水時(shí)刻。當(dāng)t=39435秒時(shí)，水塔水位35.5英尺，恰為最高水位，因此可作為第一次充水的結(jié)束時(shí)刻。充水時(shí)間為=(39435-32284)/3600=1.9864小時(shí)，也接近充水時(shí)間2小時(shí)。(2)第二次充水時(shí)間確實(shí)定當(dāng)時(shí)間t=75021秒時(shí)，水位26.97英尺,約低于最低水位27英尺，因此可作為第二次開始充水時(shí)刻。當(dāng)t=82649秒時(shí)，水泵在工作，但充水時(shí)間到達(dá)=(8264-7502)/3600=2.1189小時(shí)；但下一時(shí)刻t=85968時(shí)，水塔水位34.75英尺，低于最高水位35.50英尺。因此可將t=82649秒作為第二次充水的結(jié)束時(shí)刻,且該時(shí)刻水位為最大充水高度35.50英尺。2.計(jì)算各時(shí)刻塔內(nèi)水的體積單位轉(zhuǎn)換為1英尺=0.3048米,1升=1/3.785411加侖體積計(jì)算公式為表2不同時(shí)刻水體積表時(shí)間水體積時(shí)間水體積時(shí)間水體積0(1)6061258.967851487219.03755425540.921159371610.9542(2)67771519.95945282361.843158302612.032865767020.83925148722.949757157112.954463953422.9581(3)6777153.871456259913.875862235223.88006633974.978155209914.982260459824.98696485065.900054408115.903958932525.90836376257.006453396316.82615750087.928652537217.9317558781其中(1)表示第一段開始，(2)表示第二段開始，(3)表示第三段開始；單位：時(shí)間(小時(shí))，水流量：加侖/小時(shí)；3.計(jì)算各時(shí)刻點(diǎn)的水流量(加侖/小時(shí))水流量公式為：以上25個(gè)時(shí)刻處的水流量采用差分的方法得到，共分三段分別處理。差分公式為：對(duì)每段前兩點(diǎn)采用向前差分公式對(duì)每段最后兩點(diǎn)采用向后差分公式對(duì)每段中間點(diǎn)采用中心差分公式得到各點(diǎn)水流量表表3不同時(shí)刻水流量表時(shí)間水流量時(shí)間水流量時(shí)間水流量0(1)144048.96781102319.0375166520.92111118210.954(2)1946919.9594144951.84311006312.03282019520.8392146482.94971101212.95441894122.958(3)152203.8714879813.87581590323.8800152634.9781999114.98221805524.9869137115.9000812415.90391564625.908396347.00641016116.8261137427.9286848717.931714962其中(1)表示第一段開始，(2)表示第二段開始，(3)表示第三段開始：?jiǎn)挝唬簳r(shí)間(小時(shí))，水流量：加侖/小時(shí)；4．用三次樣條擬合流量數(shù)據(jù)對(duì)表3中25個(gè)時(shí)刻點(diǎn)的流量數(shù)據(jù)采用三次樣條插值得到一條光滑曲線，作為任意時(shí)刻的流量曲線,見圖2.1。圖2.1水塔流量圖其中*表示數(shù)據(jù)點(diǎn)，實(shí)線為樣條曲線5.一天總用水量計(jì)算一天流水總量計(jì)算:方法1:直接積分法:方法2:分段計(jì)算法第一次充水前用水(加侖)第一次充水后第二次充水前用水(加侖)[22.9581，23.88]期間用水(加侖)第一次充水期間用水:第二次充水期間用水:[23.88，24]期間用水:總共用水兩種方法結(jié)果相差6.水泵水流量計(jì)算第一次充水期間水塔體積增加充水時(shí)間:第一次充水期間水泵平均流量第二次充水期間水塔體積增加充水時(shí)間:第二次充水期間水泵平均流量那么整個(gè)充水期間水泵平均流量2.2運(yùn)用MATLAB軟件計(jì)算依據(jù)理論分析，編寫程序代碼。MATLAB程序代碼c=0.3048;%1英尺等于0.3048米p=1.0/3.785;%1升=1/3.785411加侖d=57*c;h=31.75*c;v=pi*d*d*h/4*1000*p;data=[0,3175;3316,3110;6635,3054;10619,2994;13937,2947;17921,2892;21240,2850;25223,2797;28543,2752;32284,2697;39435,3550;43318,3445;46636,3350;49953,3260;53936,3167;57254,3087;60574,3012;64554,2927;68535,2842;71854,2767;75021,2697;82649,3550;85968,3475;89953,3397;93270,3340];%原始數(shù)據(jù)t=data(:,1)/3600;%計(jì)算時(shí)間(小時(shí)為單位)v=pi*d*d*data(:,2)/100*c/4*1000*p;%計(jì)算體積%計(jì)算差分n=length(v);f=zeros(n,1);%存儲(chǔ)差分值n1=10;%計(jì)算第一段fori=1:n1ifi<=2%前兩點(diǎn)采用向前差分f(i)=abs(-3*v(i)+4*v(i+1)-v(i+2))/(2*(t(i+1)-t(i)));elseifi<=n1-2%采用三點(diǎn)中心差分公式f(i)=abs(-v(i+2)+8*v(i+1)-8*v(i-1)+v(i-2))/(12*(t(i+1)-t(i)));elseifi>=n1-1f(i)=abs(3*v(i)-4*v(i-1)+v(i-2))/(2*(t(i)-t(i-1)));endendn2=21;%計(jì)算第二段fori=n1+1:n2ifi<=n1+2%前兩點(diǎn)采用向前差分f(i)=abs(-3*v(i)+4*v(i+1)-v(i+2))/(2*(t(i+1)-t(i)));elseifi<=n2-2f(i)=abs(-v(i+2)+8*v(i+1)-8*v(i-1)+v(i-2))/(12*(t(i+1)-t(i)));elseifi>=n2-1f(i)=abs(3*v(i)-4*v(i-1)+v(i-2))/(2*(t(i)-t(i-1)));endendn3=25;%計(jì)算第三段fori=n2+1:n3ifi<=n2+2%前兩點(diǎn)采用向前差分f(i)=abs(-3*v(i)+4*v(i+1)-v(i+2))/(2*(t(i+1)-t(i)));elseifi<=n3-2f(i)=abs(-v(i+2)+8*v(i+1)-8*v(i-1)+v(i-2))/(12*(t(i+1)-t(i)));elseifi>=n3-1f(i)=abs(3*v(i)-4*v(i-1)+v(i-2))/(2*(t(i)-t(i-1)));endendplot(t,f,'r*');%畫原始點(diǎn)圖tmin=min(t);tmax=max(t);tt=tmin:0.1:tmax;%獲得離散的時(shí)間點(diǎn),用于作樣條曲線ff=spline(t,f,tt);%計(jì)算三次樣條插值holdonplot(tt,ff,'b');%畫樣條曲線xlabel('時(shí)間(小時(shí))');ylabel('流量(加侖/小時(shí))');title('水塔流量圖');holdoffdt=0.05;t2=0.5:dt:24.5;%獲得離散的時(shí)間點(diǎn),用于積分nn=length(t2);f2=spline(t,f,t2);s=(f2(1)+f2(nn)+2*sum(f2(2:nn-1)))*dt/2;%計(jì)算24小時(shí)用水量,采用復(fù)化梯形公式fprintf('(全部積分法)1天總水流量s=%8.2f\n',s);v10=v(11)-v(10);%第一次水塔增加的水dt1=t(11)-t(10);tp=t(10):dt:t(11);%第一次充水其間流出的水nn=length(tp);yp=spline(t,f,tp);%計(jì)算三次樣條插值v11=(yp(1)+yp(nn)+2*sum(yp(2:nn-1)))*dt/2;v1=v10+v11;%第一次充水總量p1=v1/dt1;%第一次充水的平均水流量v20=v(22)-v(21);%第二次水塔增加的水dt2=t(22)-t(21);tp1=t(21):dt:t(22);%第二次充水其間流出的水nn=length(tp1);yp1=spline(t,f,tp1);%計(jì)算三次樣條插值v21=(yp1(1)+yp1(nn)+2*sum(yp1(2:nn-1)))*dt/2;v2=v20+v21;%第二次充水總量p2=v2/dt2;%第二次充水的平均水流量p=(p1+p2)/2;%兩次充水平均水流量fprintf('兩次充水平均水流量p=%8.2f\n',p);v