曲線擬合的數(shù)值計算方法實驗

1、曲線擬合的數(shù)值計算方法實驗（總18頁）-本頁僅作為文檔封面，使用時請直接刪除即可-內(nèi)頁可以根據(jù)需求調(diào)整合適字體及大小-11曲線擬合的數(shù)值計算方法實驗【摘要】實際工作中，變量間未必都有線性關(guān)系，如服藥后血藥濃度與時間的關(guān)系；疾病療效與療程長短的關(guān)系；毒物劑量與致死率的關(guān)系等常呈曲線關(guān)系。曲線擬合(curvefitting)是指選擇適當(dāng)?shù)那€類型來擬合觀測數(shù)據(jù)，并用擬合的曲線方程分析兩變量間的關(guān)系。曲線直線化是曲線擬合的重要手段之一。對于某些非線性的資料可以通過簡單的變量變換使之直線化，這樣就可以按最小二乘法原理求出變換后變量的直線方程，在實際工作中常利用此直線方程繪制資料的標準工作曲線，同時根據(jù)

2、需要可將此直線方程還原為曲線方程,實現(xiàn)對資料的曲線擬合。常用的曲線擬合有最小二乘法擬合、幕函數(shù)擬合、對數(shù)函數(shù)擬合、線性插值、三次樣條插值、端點約束。關(guān)鍵詞曲線擬合、最小二乘法擬合、幕函數(shù)擬合、對數(shù)函數(shù)擬合、線性插值、三次樣條插值、端點約束一、實驗?zāi)康? .掌握曲線擬合方式及其常用函數(shù)指數(shù)函數(shù)、幕函數(shù)、對數(shù)函數(shù)的擬合2 .掌握最小二乘法、線性插值、三次樣條插值、端點約束等。3 .掌握實現(xiàn)曲線擬合的編程技巧。實驗原理1 .曲線擬合曲線擬合是平面上離散點組所表示的坐標之間的函數(shù)關(guān)系的一種數(shù)據(jù)處理方法。用解析表達式逼近離散數(shù)據(jù)的一種方法。在科學(xué)實驗或社會活動中，通過實驗或觀測得到量x與y的一組數(shù)據(jù)對(

3、X,Y)(i=1,2,m),其中各X是彼此不同的。人們希望用一類與數(shù)據(jù)的背景材料規(guī)律相適應(yīng)的解析表達式，y=f(x,c)來反映量x與y之間的依賴關(guān)系，即在一定意義下最佳”地逼近或擬合已知數(shù)據(jù)。f(x,c)常稱作擬合模型，式中C=(Cl,C2,n)是一些待定參數(shù)。當(dāng)C在f中線性出現(xiàn)時，稱為線性模型，否則稱為非線性模型。有許多衡量擬合優(yōu)度的標準，最常用的一種做法是選擇參數(shù)C使得擬合模型與實際觀測值在各點的殘差(或離差)qy-(fk,c)的加權(quán)平方和達到最小，此時所求曲線稱作在加權(quán)最小二乘意義下對數(shù)據(jù)的擬合曲線。有許多求解擬合曲線的成功方法，對于線性模型一般通過建立和求解方程組來確定參數(shù)，從而求得擬

4、合曲線。至于非線性模型，則要借助求解非線性方程組或用最優(yōu)化方法求得所需參數(shù)才能得到擬合曲線，有時稱之為非線性最小二乘擬合。曲線擬合：貝塞爾曲線與路徑轉(zhuǎn)化時的誤差。值越大，誤差越大；值越小，越精確。2 .最小二乘法擬合：最小二乘法(又稱最小平方法)是一種數(shù)學(xué)優(yōu)化技術(shù)。它通過最小化誤差的平方和尋找數(shù)據(jù)的最佳函數(shù)匹配。利用最小二乘法可以簡便地求得未知的數(shù)據(jù)，并使得這些求得的數(shù)據(jù)與實際數(shù)據(jù)之間誤差的平方和為最小。最小二乘法還可用于曲線擬合。其他一些優(yōu)化問題也可通過最小化能量或最大化嫡用最小二乘法來表達。函數(shù)曲線為：y=Ax+B其中系數(shù)滿足下列的正規(guī)方程：NNN2XkAXkBXkykk1k1k1NNXk

5、ANBykk1k13 .幕函數(shù)擬合：函數(shù)曲線為：設(shè)xkyk。有N個點，其中橫坐標是確定的。最小二乘幕函數(shù)擬合曲線的系數(shù)A為：NNM2M、A(Xky)(Xk)k1k14 .對數(shù)函數(shù)擬合：對數(shù)函數(shù)(lograrithmicfunction)的標準式形式為YablnX(X0)b>0時，Y隨X增大而增大，先快后慢；b<0時，Y隨X增大而減少，先快后慢，見圖(c)、(d)。當(dāng)以Y和lnX繪制的散點圖呈直線趨勢時，可考慮采用對數(shù)函數(shù)描述Y與X之間的非線性關(guān)系，式中的b和a分別為斜率和截距。更一般的對數(shù)函數(shù)Y=a+bln(X+k)式中k為一常量，往往未知。(a)lnY=lna+bX(b)lnY=

6、lna-bX(c)Y=a+blnX(d)Y=a-blnX5 .線性插值:在代數(shù)才S值中，為了提高插值多項式對函數(shù)的逼近程度一般是增加節(jié)點的個數(shù)，即提高多項式的次數(shù)，但這樣做往往不能達到預(yù)想的效果。如下圖所示：f(x)=1/(1+x2訓(xùn)果在區(qū)間-5,5上取7個等距節(jié)點:xk=5*(k/3-1)(k=0,1,2,6)，由lagrange插值公式可得到f(x)的次L7(x)。如圖所示：L7(x；僅在區(qū)問的中部能較好的逼近函數(shù)f(x),在其它部位差異較大，而且越接近端點，逼近效果越差?？梢宰C明，當(dāng)節(jié)點無限加密時，Ln(x也只能在很小的范圍內(nèi)收斂，這一現(xiàn)象稱為Runge現(xiàn)象。它表明通過增加節(jié)點來提高逼近

7、程度是不適宜的，因而不采用高次多項式插值。如果我們把以上的點用直線連接起來，顯然比L7(x度更逼近f(x)0這就是分段線性插值。而在實際應(yīng)用中通常采用分段低次插值來提高近似程度，比如可用分段線性插值或分段三次埃爾米特插值來逼近已知函數(shù)，但它們的總體光滑性較差。為了克服這一缺點，一種全局化的分段插值方法一一三次樣條插值成為比較理想的工具。6 .三次樣條插值：設(shè)xk,ykN0有N+1個點，其中ax0x1x2.xNbo如果存在N個三次多項式Sk(x),系數(shù)為Sko,Ski,Sk2,Sk3滿足如下性質(zhì)：、c/、/、/、2/、3S(x)Sk(x)sk,0Sk,l(xxk)Sk,2(xxk)sk,3(xx

8、k)xxk,xki,k0,1,.,N1S(xk)ykk0,1,.NS(xki)Ski(xki)k0,1,.N2S'(xki)S'k1(xk1)k0,1,.N2S''(xk1)S''k1(xk1)k0,1,.N2則成函數(shù)S(x)為三次樣條函數(shù)。7 .端點約束：緊壓樣條：存在唯一的三次樣條曲線，其一階導(dǎo)數(shù)的邊界條件是：S'(a)d0,S'(b)dNnatural樣條：存在唯一的三次樣條曲線，它的自由邊界條件是：S''(a)0,S"(b)0外推樣條：存在唯一的三次樣條曲線，其中通過對點xi和x2進行外推得到Sa

9、,同時通過對點X(n-1)ftX(N-2)進行外才t得到Sb。端點曲率調(diào)整：存在唯一的三次樣條曲線，其中二階導(dǎo)數(shù)的邊界條件Sa和Sb是確定的。拋物線終結(jié)樣條：存在唯一的三次樣條曲線，其中二階在區(qū)間X0,Xi內(nèi)S(x)0,而在Xn-1,Xn內(nèi)S(x)0。三、實驗內(nèi)容1.P2021胡克定律指出F=kx,其中F是拉伸彈簧的拉力(單位為盎司)，x為拉伸的長度(單位為英寸)。根據(jù)下列試驗數(shù)據(jù)，求解拉伸常量k的近似值。(a)xkFk(b)XkFk2.P2151洛杉磯(美國城市)郊區(qū)11月8日的溫度記錄入下表所示，其中共有24個數(shù)據(jù)點。(a)根據(jù)例中的處理過程(使用fmins命令)，對給定的數(shù)據(jù)求解最小二乘

10、曲線f(x)Acos(Bx)Csin(Dx)E。(b)求E2(f)o(c)在同一坐標系下畫出這些點集和(a)得出的最小二乘曲線。時間，.溫度時間，.溫度16615826625836535846445856355766365776275786185819609601060106411591167午夜58正午683.P2291一個轎車在時間Tk時經(jīng)過的距離dk,如下表所示。使用程序，并根據(jù)一階導(dǎo)數(shù)邊界條件S'(0)0,S'(8)98,求這些數(shù)據(jù)的三次緊壓樣條插值。時間，tk02468距離，dk0401603004804.P2385美國洛杉磯郊區(qū)11月9日的溫度(華氏溫度)如表所示。

11、采用24小時(a)求三角多項式T7(x)(b)在同一坐標系下，畫出圖T7(x)和24個數(shù)據(jù)點。(c)使用本地的溫度情況重新求解問題(a)和問題(b)。時間，溫度時間，溫度1661582662583653584644585635576636577627578618589609601060106411591167午夜58正午685.P2461編寫Matlab程序，生成并繪制組合貝塞爾曲線。利用該程序生成和繪制過3個控制點集(0,0),(1,2),(1,1),(3,0),(3,0),(4,-1),(5,-2),(6,1),(7,0),(7,0),(4,-3),(2,-1),(0,0)的貝塞爾曲線。四

12、、實驗結(jié)果及分析1實驗描述：由題意可知，此題需要用最小二乘法進行計算，因為已知函數(shù)的5個插點并且知道它們的x,y的值。且函數(shù)的表達式為F=kx,所以只需用方程中NN(xk)Ayk便可計算出k的數(shù)值。k1k1定義一個動態(tài)數(shù)組a，用來依次存取x和y的插值。其中x,y的插值通過鍵盤手動輸入并賦予給a中的元素。然后通過相應(yīng)的求和和基本運算便可以得到相應(yīng)插值下的最小二乘法的函數(shù)表達式。(其中k保留小數(shù)點后4位)具體函數(shù)運行效果如下：(a)請在此輸入x的各插值請在此輸入y的各插值此函數(shù)的最小二乘法曲線表達式為Y=*x請按任意鍵繼續(xù)。(b)請在此輸入x的各插值請在此輸入y的各插值此函數(shù)的最小二乘法曲線表達式

13、為Y=*x請按任意鍵繼續(xù)。結(jié)果分析：易得，最小二乘法多項式計算可以很好的做出較為精確的最小二乘法擬合曲線，并且程序通用性高，只要輸入相應(yīng)的插值便可以計算出結(jié)果，結(jié)果系數(shù)的小數(shù)點有效位同時也可以自己確定。實驗描述：給出的最小二乘曲線表達式為：f(x)Acos(Bx)Csin(Dx)E其中變量有5個，而給出的數(shù)據(jù)點有24個，在C語言中可以用牛頓-拉夫森算法迭代計算分別得出5個變量的值，但是方法繁瑣，且迭代計算量龐大，因此這里采用Matlab進行相關(guān)的計算，調(diào)用fminsearch函數(shù)，求得當(dāng)5個參量都為1附近時候多項式的最小值，此時便可求出此5個參變量的值.然后繼續(xù)通過Matlab,將得到的公式以

14、及各點，用plot函數(shù)，便可以求得。實驗結(jié)果：運行matlab結(jié)果如下：>>fminsearch('fiveOne',11111)ans=此時的所求值便為函數(shù)的待定系數(shù)。所以可得最小二乘曲線的表達式為：f(x)15.7225cos(1.3717x)15.5359sin(1.2768)60.3579然后進行相應(yīng)的繪制圖形便可以求出所要求出的結(jié)果。結(jié)果分析：通過最小二乘法多項式同樣可以繪制出三角函數(shù)的曲線。并且程序通用性高，只要輸入相應(yīng)的插值便可以計算出結(jié)果，結(jié)果系數(shù)的小數(shù)點有效位同時也可以自己確定。1實驗描述：由題意可知，由于這里涉及到了樣條線的運算，計算較為復(fù)雜。且

15、要涉及到畫圖的部分，所以此處采用matlab計算較為方便快捷。而書本上給出了一個用來計算三次緊壓樣條曲線的可調(diào)用函數(shù)，現(xiàn)在將其引用，并根據(jù)已知點得出相應(yīng)的曲線。實驗結(jié)果：代入后得出的結(jié)果如下所示:> >X=02468;> >Y=040160300480;> >S=csfit(X,Y,0,98)S=00由結(jié)果可知此插值為在區(qū)間0,8中有三個分別劃分了0,2,2,4,4,6,6,8四個區(qū)間的插點。且多項式分別為_3_2S0(x)0.8125x38.375x2.32S(x)2.4375(x2)13.25(x2)43.25(x2)40S2(x)1.4375(x4)3

16、1.375(x4)267(x4)160S3(x)0.8125(x6)37.25(x6)278.75(x6)300在matlab中通過polyval作出相應(yīng)的曲線，再用plot函數(shù)便可繪制圖線。繪制后圖線如下：平滑的樣條線此時便可以直觀的看到一個結(jié)果分析:通過給定的一些數(shù)據(jù)點，便可以繪制出過這些點的相應(yīng)的樣條線，通過觀察能發(fā)現(xiàn)樣條線的平滑度與你選擇的樣條線類型以及數(shù)據(jù)點的分布有一定關(guān)聯(lián)。不僅僅是緊壓樣條線，相關(guān)其它的線也可以用同樣的方法一一給出。5實驗描述：由題意可知，題目是想求出所繪制數(shù)據(jù)點的三角多項式的逼近，三角多項式逼近的公式為：a0MTm(x)(ajcos(jx)bjSin(jx)2ji

17、止匕外，X為離散傅里葉級數(shù)，滿足條件：Xjj,j0,1,.,NjN實驗結(jié)果：而所給的點X為1,24的離散數(shù)值點，所以無法直接對其作出逼近公式，需要進行尺度變換，將X點事專換為：Xjj,j0,1,.,24j24(1)通過matlab繪制出相關(guān)的三角多項式曲線，然后同樣通過matlab的繪制點，將點繪制到這個曲線之中，具體的matlab代碼如下：> >holdon;聯(lián)持圖形不動，繪制新的圖形入曲線中> >X=123456789101112131415161718192021222324;數(shù)據(jù)點的X的取值> >Y=5858585857575758606467686

18、66665646363626160605958;%數(shù)據(jù)點的y的取值> >plot(X,Y：xk)噬制出數(shù)據(jù)點(2)然后便可以畫出如圖所示的插值數(shù)據(jù)點。結(jié)果分析：三角多項式的曲線擬合度非常高，能很好的繪制出圖像的具體形式而且曲線平滑，但是它需要滿足x屬于-pi到pi的區(qū)間內(nèi)。1實驗描述：由題意可知，首先以33階貝塞爾函數(shù)為帶求函數(shù)。其求解格式如下：P(t)Bi,3(t)其中b為伯恩斯坦多項式，其定義如下：i03i3iBi,3(t)ti(1t)3ii用matlab先設(shè)置一個參數(shù)t,然后再根據(jù)公式，通過所給的數(shù)據(jù)點將伯恩斯坦多項式，以及x,y的關(guān)于參數(shù)t的多項式求出來。然后再把t設(shè)置為精

19、度足夠大的單位序列，繪制圖線即可得出貝塞爾的效果。實驗結(jié)果：其matlab運行后結(jié)果如下：【以第一組控制點為例】>>X=0113X=0113>>Y=0210Y=0210>>fiveTwo(X,Y)x=3*t*(t-1)A2-3*tA2*(t-1)+3*tA3y=6*t*(t-1)A2-3*tA2*(t-1)21且繪制出第一組控制點所在位置以及三階貝塞爾曲線如下所示:同理，第2,3組控制點所作的圖形如下所示:結(jié)果分析：可以通過Matlab函數(shù)繪制出三階的貝塞爾函數(shù)。只要給出控制點，便可自動繪制出所控制的三階貝塞爾函數(shù)以及控制點的位標。由此可以觀察到貝塞爾與控制

20、點的約束作用，以及所求得貝塞爾函數(shù)是個相對平滑的曲線。五、實驗結(jié)論曲線擬合對于實際的工程以及理論推導(dǎo)問題都有著重大的作用。在具體的實驗，數(shù)據(jù)分析上，往往我們只有巨量的已知離散點，想從離散點中得出函數(shù)表達式則就需要曲線擬合進行，根據(jù)不同的要求，我們可以選擇最小二乘法的幕函數(shù)擬合或者是一次函數(shù)，二次函數(shù)擬合，抑或是精度非常高的傅里葉變換的三角函數(shù)擬合。同時在建模方面，貝塞爾函數(shù)的引用也從數(shù)學(xué)層面解決了如何用計算機繪制出光滑圓潤的曲線，在一些設(shè)計軟件中，例如3dmax和maya的三維建模，AutoCAD的工程建模，貝塞爾運算對于曲線曲面的設(shè)計有著舉足輕重的作用。附件（代碼）1#include<

21、>#include<>#include<>voidmain()externfloatafunction(floatb);externfloatbfunction(floatz,floatv);float*a,y,q,c;inti;a=newfloat5;cos(i*B)+C.*sin(i*D)+E-y(i)A2+z;%函數(shù)的線性最小二乘法的多項式end之后在Matlab的命令窗口輸入如下語句：>>fminsearch(fiveOne,11111)(繪制圖形方程組)>>x1=0123456789101112131415161718192021

22、222324;%數(shù)值點的x取值%數(shù)值點的x 取值%繪制出24 個數(shù)%繪制出此函數(shù)的二乘>>y1=585858585757575860646768666665646363626160605958;>>plot(x1,y1,'x')值點的圖形>>x=linspace(0,25,100);>>holdon>>plot(x,y,'k')法后的函數(shù)曲線1(三次緊壓樣條線計算函數(shù))functionS=csfit(x,y,dx0,dxn)N=length(x)-1;H=diff(x);D=diff(y)./H;%d(

23、k)的值A(chǔ)=H(2:N-1);B=2*(H(1:N-1)+H(2:N);C=H(2:N);U=6*diff(D);%U(k)B值B(1)=B(1)-H(1)/2;U(1)=U(1)-3*(D(1)-dx0);B(N-1)=B(N-1)-H(N)/2;U(N-1)=U(N-1)-3*(dxn-D(N);fork=2:N-1temp=A(k-1)/B(k-1);B(k)=B(k)-temp*C(k-1);U(k)=U(k)-temp*U(k-1);endM(N)=U(N-1)/B(N-1);fork=N-2:-1:1M(k+1)=(U(k)-C(k)*M(k+2)/B(k);endM(1)=3*(

24、D(1)-dx0)/H(1)-M(2)%B次緊壓約束中S"(0)的值M(N+1)=3*(dxn-D(N)/H(N)-M(N)/2;%三次緊壓約束中S"(k)的值fork=0:N-1S(k+1,1)=(M(k+2)-M(k+1)/(6*H(k+1);%S(k,3)勺值S(k+1,2)=M(k+1)/2;%S(k,2)勺值S(k+1,3)=D(k+1)-H(k+1)*(2*M(k+1)+M(k+2)/6;%S(k,1)勺值S(k+1,4)=y(k+1);%S(k,0的值End(曲線繪制函數(shù)程序)> >x1=0:.01:2;y1=polyval(S(1,:),x1-X

25、(1);%第一段樣條線> >x2=2:.01:4;y2=polyval(S(2,:),x2-X(2);%第二段樣條線> >x3=4:.01:6;y3=polyval(S(3,:),x3-X(3);%第三段樣條線> >x4=6:.01:8;y4=polyval(S(4,:),x4-X(4);%第四段樣條線> >plot(x1,y1,x2,y2,x3,y3,x4,y4,X,Y,'x')%繪制連接成完整曲線，并且標注數(shù)據(jù)點的位置。#include<>#include<>#include<>voidma

26、in()externfloatafunction(floatX,floatY,inttemp);externfloatbfunction(floatX,floatY,inttemp);inti;floata24=1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24;floatb24=58,58,58,58,57,57,57,58,60,64,67,68,66,66,65,64,63,63,62,61,60,60,59,58;float*c;float*d;c=newfloat7,d=newfloat7;for(i=0;i&

27、lt;=24;i+)ci=afunction(a,b,i);di=bfunction(a,b,i);printf("由所給的數(shù)據(jù)點可以得出n三角多項式T7(x)中a的系數(shù)分別為n");printf("%f%f%f%f%f%f%fn",c0,c1,c2,c3,c4,c5,c6);printf("由所給的數(shù)據(jù)點可以得出n三角多項式T7(x)中b的系數(shù)分別為n");printf("%f%f%f%f%f%f%fn",d0,d1,d2,d3,d4,d5,d6);system("pause");return;floatafunction(floatX,floatY,inttemp)floatA=0;inti;for(i=0;i<=24;i+)