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

1、樣條插值函數(shù)及應(yīng)用摘要樣條函數(shù)具有廣泛的應(yīng)用，是現(xiàn)代函數(shù)論的一個(gè)十分活躍的分支，是計(jì)算方法的主要基礎(chǔ)和工具之一，由于生產(chǎn)和科學(xué)技術(shù)向前發(fā)展的推動(dòng)以及電子計(jì)算機(jī)廣泛應(yīng)用的需要，人們便更多地應(yīng)用這個(gè)工具，也更深刻的認(rèn)識(shí)了它的本質(zhì)。在實(shí)際問(wèn)題中所遇到許多函數(shù)往往很復(fù)雜，有些甚至是很難找到解析表達(dá)式的。根據(jù)函數(shù)已有的數(shù)據(jù)來(lái)計(jì)算函數(shù)在一些新的點(diǎn)處的函數(shù)值，就是插值法所需要解決的問(wèn)題。插值法是數(shù)值逼近的重要方法之一，它是根據(jù)給定的自變量值和函數(shù)值,求取未知函數(shù)的近似值。早在一千多年前,我國(guó)科學(xué)家就在研究歷法時(shí)就用到了線(xiàn)性插值和二次插值。而在實(shí)際問(wèn)題中,有許多插值函數(shù)的曲線(xiàn)要求具有較高的光滑性,在整個(gè)曲線(xiàn)中

2、,曲線(xiàn)不但不能有拐點(diǎn),而且曲率也不能有突變。因此，對(duì)于插值函數(shù)必須二次連續(xù)可微且不變號(hào) ,這就需要用到三次樣條插值。關(guān)鍵詞 三次樣條函數(shù)；插值法目 錄引 言1第一章三次樣條插值21.1 樣條插值函數(shù)簡(jiǎn)介21.2 三次樣條函數(shù)應(yīng)用3第二章AMCM91A 估計(jì)水塔水流量52.1 理論分析及計(jì)算62.2運(yùn)用MATLAB軟件計(jì)算9參考文獻(xiàn)142引 言樣條函數(shù)具有廣泛的應(yīng)用，是現(xiàn)代函數(shù)論的一個(gè)十分活躍的分支，是計(jì)算方法的主要基礎(chǔ)和工具之一，由于生產(chǎn)和科學(xué)技術(shù)向前發(fā)展的推動(dòng)以及電子計(jì)算機(jī)廣泛應(yīng)用的需要，人們便更多地應(yīng)用這個(gè)工具，也更深刻的認(rèn)識(shí)了它的本質(zhì)。上世紀(jì)四十年代，在研究數(shù)據(jù)處理的問(wèn)題中引出了樣條函數(shù)

3、，例如，在1946年Schoenberg將樣條引入數(shù)學(xué),即所謂的樣條函數(shù)，直到五十年代，還多應(yīng)用于統(tǒng)計(jì)數(shù)據(jù)的處理方面，從六十年代起，在航空、造船、汽車(chē)等行業(yè)中，開(kāi)始大量采用樣條函數(shù)。在我國(guó)，從六十年代末開(kāi)始，從船體數(shù)學(xué)放樣到飛機(jī)外形設(shè)計(jì)，逐漸出現(xiàn)了一個(gè)使用樣，逐漸出現(xiàn)了一個(gè)使用樣條函數(shù)的熱潮，并推廣到數(shù)據(jù)處理的許多問(wèn)題中。在實(shí)際生活中有許多計(jì)算問(wèn)題對(duì)插值函數(shù)的光滑性有較高的要求，例如飛機(jī)機(jī)翼外形、發(fā)動(dòng)機(jī)進(jìn)、排氣口都要求有連續(xù)的二階導(dǎo)數(shù)，用三次樣條繪制的曲線(xiàn)不僅有很好的光滑度，而且當(dāng)節(jié)點(diǎn)逐漸加密時(shí)其函數(shù)值整體上能很好地逼近被插函數(shù)，相應(yīng)的導(dǎo)數(shù)值也收斂于被插函數(shù)的導(dǎo)數(shù)值，不會(huì)發(fā)生“龍格現(xiàn)象”?，F(xiàn)在

4、國(guó)內(nèi)外學(xué)者對(duì)這方面的研究也越來(lái)越重視，根據(jù)我們的需要來(lái)解決不同的問(wèn)題，而且函數(shù)的形式也在不斷地改進(jìn)，長(zhǎng)期以來(lái)很多學(xué)者致力于樣條插值的研究，對(duì)三次樣條的研究已相當(dāng)成熟。14第一章三次樣條插值1.1 樣條插值函數(shù)簡(jiǎn)介 在實(shí)際問(wèn)題中所遇到許多函數(shù)往往很復(fù)雜，有些甚至是很難找到解析表達(dá)式的。有時(shí)通過(guò)實(shí)驗(yàn)或者數(shù)值計(jì)算所得到的也只是一些離散的點(diǎn)上的函數(shù)值，即。根據(jù)函數(shù)已有的數(shù)據(jù)來(lái)計(jì)算函數(shù)在一些新的點(diǎn)處的函數(shù)值，就是插值法所需要解決的問(wèn)題。插值法的基本思想就是，首先根據(jù)已有的函數(shù)值來(lái)構(gòu)造一個(gè)簡(jiǎn)單的函數(shù)作為的近似表達(dá)式，然后用來(lái)計(jì)算新的點(diǎn)上的函數(shù)值作為的近似值。通常可以選多項(xiàng)式函數(shù)作為近似函數(shù)，因?yàn)槎囗?xiàng)式具有

5、各階導(dǎo)數(shù)，求值也比見(jiàn)方便。常用的有Lagrange插值、Newton插值、Hermite插值和樣條插值。線(xiàn)性插值在分段點(diǎn)上僅連續(xù)而不可導(dǎo)，三次埃爾米特插值有連續(xù)的一階導(dǎo)數(shù)，這樣的光滑程度常不能滿(mǎn)足物理問(wèn)題的需要，樣條函數(shù)可以同時(shí)解決這兩個(gè)問(wèn)題,使插值函數(shù)既是低階分段函數(shù),又是光滑的函數(shù)，并且只需在區(qū)間端點(diǎn)提供某些導(dǎo)數(shù)信息。三次樣條函數(shù)定義：設(shè)在區(qū)間上取個(gè)節(jié)點(diǎn)，函數(shù)在各個(gè)節(jié)點(diǎn)處的函數(shù)值為若滿(mǎn)足： (1) (2) 在區(qū)間上，具有連續(xù)的二階導(dǎo)數(shù)； (3) 在區(qū)間上,是三次的多項(xiàng)式；則稱(chēng)是函數(shù)在區(qū)間上的三次樣條插值函數(shù)。 由以上定義可以看出，雖然每個(gè)子區(qū)間上的多項(xiàng)式可以各不相同，但在相鄰子區(qū)間的連接處

6、卻是光滑的。因此，樣條插值也稱(chēng)為分段光滑插值。從定義知要求出，在每一個(gè)小區(qū)間上確定4個(gè)待定系數(shù)，共有n個(gè)小區(qū)間，故應(yīng)有4n個(gè)參數(shù)。根據(jù)在上二階導(dǎo)數(shù)連續(xù)，在節(jié)點(diǎn)出滿(mǎn)足連續(xù)性條件共有3n-3個(gè)條件，再加上滿(mǎn)足插值條件共有4n-2個(gè)條件，因此還需要2個(gè)條件才能確定。通?？稍趨^(qū)間端點(diǎn)上各加一個(gè)條件（稱(chēng)為邊界條件），可根據(jù)實(shí)際問(wèn)題的要求給定，通常有以下三種：（1）已知端點(diǎn)的一階導(dǎo)數(shù)值，即（2）倆端點(diǎn)的二階導(dǎo)數(shù)已知，即其特殊情況稱(chēng)為自然邊界條件。（3）當(dāng)是以為周期的函數(shù)時(shí)，則要求也是周期函數(shù)。這時(shí)邊界條件應(yīng)滿(mǎn)足而此時(shí)。這樣確定的樣條函數(shù)，稱(chēng)為周期函數(shù)。1.2 三次樣條函數(shù)應(yīng)用作函數(shù)在取間隔為0.1的點(diǎn)圖，

7、用插值進(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ì)的用水率以及

8、每天總的用水量,但許多社區(qū)并沒(méi)有測(cè)量水流入或流出當(dāng)?shù)厮乃康脑O(shè)備,他們只能代之以每小時(shí)測(cè)量水塔中的水位,精度在0.5%以?xún)?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í)間32284359

9、3539332394354331846636499535393657254水位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è)量開(kāi)始的以秒為單位的時(shí)刻，以及該時(shí)刻的高度單位為百分之一英尺的水塔中

10、水位的測(cè)量值，例如，3316秒后，水塔中的水位達(dá)到31.10英尺。水塔是一個(gè)垂直圓形柱體，高為40英尺，直徑57英尺，通常當(dāng)水塔的水位降至27.00英尺時(shí)水泵開(kāi)始向水塔充水，而當(dāng)水塔的水位升至35.50英尺時(shí)水泵停止工作。2.1 理論分析及計(jì)算1. 水塔充水時(shí)間的確定(1) 第一次充水時(shí)間的確定當(dāng)時(shí)間t=32284秒時(shí)，水位26.97英尺,約低于最低水位27英尺，因此可作為第一次開(kāi)始充水時(shí)刻。當(dāng)t=39435秒時(shí)，水塔水位35.5英尺，恰為最高水位，因此可作為第一次充水的結(jié)束時(shí)刻。充水時(shí)間為=(39435-32284)/3600=1.9864小時(shí)，也接近充水時(shí)間2小時(shí)。(2) 第二次充水時(shí)間的

11、確定當(dāng)時(shí)間t=75021秒時(shí)，水位26.97英尺,約低于最低水位27英尺，因此可作為第二次開(kāi)始充水時(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.037

12、55425540.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)表示

13、第一段開(kāi)始，(2) 表示第二段開(kāi)始，(3) 表示第三段開(kāi)始；單位：時(shí)間(小時(shí))，水流量：加侖/小時(shí)；3.計(jì)算各時(shí)刻點(diǎn)的水流量(加侖/小時(shí))水流量公式為： 以上25個(gè)時(shí)刻處的水流量采用差分的方法得到，共分三段分別處理。差分公式為：(1) 對(duì)每段前兩點(diǎn)采用向前差分公式(2) 對(duì)每段最后兩點(diǎn)采用向后差分公式(3) 對(duì)每段中間點(diǎn)采用中心差分公式得到各點(diǎn)水流量表表3 不同時(shí)刻水流量表時(shí)間水流量時(shí)間水流量時(shí)間水流量0 (1)144048.96781102319.0375166520.92111118210.954(2)1946919.9594144951.84311006312.03282019520.8

14、392146482.94971101212.95441894122.958(3)152203.8714879813.87581590323.8800152634.9781999114.98221805524.9869137115.9000812415.90391564625.908396347.00641016116.8261137427.9286848717.931714962其中(1)表示第一段開(kāi)始，(2) 表示第二段開(kāi)始，(3) 表示第三段開(kāi)始：?jiǎn)挝唬簳r(shí)間(小時(shí))，水流量：加侖/小時(shí)；4用三次樣條擬合流量數(shù)據(jù)對(duì)表3中25個(gè)時(shí)刻點(diǎn)的流量數(shù)據(jù)采用三次樣條插值得到一條光滑曲線(xiàn)，作為任意時(shí)刻的流

15、量曲線(xiàn),見(jiàn)圖2.1。圖2.1 水塔流量圖其中*表示數(shù)據(jù)點(diǎn)，實(shí)線(xiàn)為樣條曲線(xiàn)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ì)算1. 依據(jù)理論分析，編寫(xiě)程序代碼。MATLAB程序代

16、碼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

17、;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ì)算第一段for i=1:n1 if i<=2 %前兩點(diǎn)采用向前差分 f(i)=abs(-3*v(i)+4*v(i+1)-v(i+2)/(2*(t(i+1)-t(i); elseif i<=n1-2 %采

18、用三點(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);elseif i>=n1-1 f(i)=abs(3*v(i)-4*v(i-1)+v(i-2)/(2*(t(i)-t(i-1); end endn2=21; %計(jì)算第二段for i=n1+1:n2 if i<=n1+2 %前兩點(diǎn)采用向前差分 f(i)=abs(-3*v(i)+4*v(i+1)-v(i+2)/(2*(t(i+1)-t(i); elseif i<=n2-2 f(i)=abs(-v(i+2)+8*v(i+1)-8*v(i-1)+v(

19、i-2)/(12*(t(i+1)-t(i);elseif i>=n2-1 f(i)=abs(3*v(i)-4*v(i-1)+v(i-2)/(2*(t(i)-t(i-1); end endn3=25; %計(jì)算第三段for i=n2+1:n3 if i<=n2+2 %前兩點(diǎn)采用向前差分 f(i)=abs(-3*v(i)+4*v(i+1)-v(i+2)/(2*(t(i+1)-t(i); elseif i<=n3-2 f(i)=abs(-v(i+2)+8*v(i+1)-8*v(i-1)+v(i-2)/(12*(t(i+1)-t(i);elseif i>=n3-1 f(i)=ab

20、s(3*v(i)-4*v(i-1)+v(i-2)/(2*(t(i)-t(i-1); end endplot(t,f,'r*'); %畫(huà)原始點(diǎn)圖tmin=min(t); tmax=max(t);tt=tmin:0.1:tmax; %獲得離散的時(shí)間點(diǎn),用于作樣條曲線(xiàn)ff=spline(t,f,tt); %計(jì)算三次樣條插值hold onplot(tt,ff,'b'); %畫(huà)樣條曲線(xiàn)xlabel('時(shí)間(小時(shí))');ylabel('流量(加侖/小時(shí))');title('水塔流量圖');hold offdt=0.05;t2=

21、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.2fn',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)

22、+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.2fn'