《曲線擬合的數(shù)值計算方法實驗.doc》由會員分享，可在線閱讀，更多相關《曲線擬合的數(shù)值計算方法實驗.doc（19頁珍藏版）》請在裝配圖網(wǎng)上搜索。

1、19 《數(shù)值方法》實驗報告 曲線擬合的數(shù)值計算方法實驗 【摘要】實際工作中，變量間未必都有線性關系，如服藥后血藥濃度與時間的關系；疾病療效與療程長短的關系；毒物劑量與致死率的關系等常呈曲線關系。曲線擬合（curve fitting）是指選擇適當?shù)那€類型來擬合觀測數(shù)據(jù)，并用擬合的曲線方程分析兩變量間的關系。曲線直線化是曲線擬合的重要手段之一。對于某些非線性的資料可以通過簡單的變量變換使之直線化，這樣就可以按最小二乘法原理求出變換后變量的直線方程，在實際工作中常利用此直線方程繪制資料的標準工作曲線，同時根據(jù)需要可將此直線方程還原為曲線方程，實現(xiàn)對資料的曲線擬合。常用的曲

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

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

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

5、中系數(shù)滿足下列的正規(guī)方程： 3. 冪函數(shù)擬合： 函數(shù)曲線為： 設有N個點，其中橫坐標是確定的。最小二乘冪函數(shù)擬合曲線的系數(shù)A為： 、 4. 對數(shù)函數(shù)擬合： 對數(shù)函數(shù)（lograrithmic function）的標準式形式為 b>0時，Y隨X增大而增大，

6、先快后慢；b<0時，Y隨X增大而減少，先快后慢，見圖12.4(c)、(d)。當以Y和lnX繪制的散點圖呈直線趨勢時，可考慮采用對數(shù)函數(shù)描述Y與X之間的非線性關系，式中的b和a分別為斜率和截距。 更一般的對數(shù)函數(shù) Y=a+bln(X+k) 式中k為一常量，往往未知。 (a) lnY=lna+bX (b) lnY=lna-bX (c) Y=a+blnX (d) Y=a-blnX 5. 線性插值： 在代數(shù)插值中,為了提高插值多項式對函數(shù)的逼近程度一般是增加節(jié)點的個數(shù)，即提高多項式的次數(shù)，但這樣做往往不能達到預想的效果。 如下圖所示：

7、 f(x) = 1 / (1 + x2)如果在區(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明，當節(jié)點無限加密時，Ln(x)也只能在很小的范圍內收斂，這一現(xiàn)象稱為Runge現(xiàn)象。它表明通過增加節(jié)點來提高逼近程度是不適宜的，因而不采用高次多項式插值。 如果我們把以上的點用直線連接起來，顯然比L7(x)要更逼近f(x)。這就是分段線性插值。 而在實際應用中通常采用分

8、段低次插值來提高近似程度，比如可用分段線性插值或分段三次埃爾米特插值來逼近已知函數(shù)，但它們的總體光滑性較差。為了克服這一缺點，一種全局化的分段插值方法——三次樣條插值成為比較理想的工具。 6. 三次樣條插值： 設有N+1個點，其中。如果存在N個三次多項式，系數(shù)為滿足如下性質： 則成函數(shù)S(x)為三次樣條函數(shù)。 7. 端點約束： 緊壓樣條：存在唯一的三次樣條曲線，其一階導數(shù)的邊界條件是：

9、 natural樣條：存在唯一的三次樣條曲線，它的自由邊界條件是： 外推樣條：存在唯一的三次樣條曲線，其中通過對點x1和x2進行外推得到 ，同時通過對點X(n-1)和X(N-2)進行外推得到。 端點曲率調整：存在唯一的三次樣條曲線，其中二階導數(shù)的邊界條件和是確定的。 拋物線終結樣條：存在唯一的三次樣條曲線，其中二階

10、在區(qū)間[X0，X1]內，而在[Xn-1,Xn]內。 3、 實驗內容 1. P202 1 胡克定律指出F=kx，其中F是拉伸彈簧的拉力（單位為盎司），x為拉伸的長度（單位為英寸）。根據(jù)下列試驗數(shù)據(jù)，求解拉伸常量k的近似值。 xk Fk 0.2 3.6 0.4 7.3 0.6 10.9 0.8 14.5 1.0 18.2 Xk Fk 0.2 5.3 0.4 10.6 0.6 15.9 0.8 21.2 1.0 26.4 （a） （b）

11、 2. P215 1 洛杉磯（美國城市）郊區(qū)11月8日的溫度記錄入下表所示，其中共有24個數(shù)據(jù)點。 （a）根據(jù)例5.5中的處理過程（使用fmins命令），對給定的數(shù)據(jù)求解最小二乘曲線 。 （b）求 （c）在同一坐標系下畫出這些點集和(a)得出的最小二乘曲線。 時間，p.m. 溫度 時間，a.m. 溫度 1 66 1 58 2 66 2 58 3 65 3 58 4 64 4 58 5 63 5 57 6 63 6 57 7 62 7 57 8 61 8 58 9 60 9 6

12、0 10 60 10 64 11 59 11 67 午夜 58 正午 68 3. P229 1 一個轎車在時間Tk時經過的距離dk，如下表所示。使用程序5.3，并根據(jù)一階導數(shù)邊界條件,求這些數(shù)據(jù)的三次緊壓樣條插值。 時間，tk 0 2 4 6 8 距離，dk 0 40 160 300 480 4. P238 5 美國洛杉磯郊區(qū)11月9日的溫度（華氏溫度）如表5.10所示。采用24小時制。 （a）求三角多項式 （b）在同一坐標系下，畫出圖和24個數(shù)據(jù)點。 （c）

13、使用本地的溫度情況重新求解問題(a)和問題(b)。 時間,p.m 溫度 時間,a.m 溫度 1 66 1 58 2 66 2 58 3 65 3 58 4 64 4 58 5 63 5 57 6 63 6 57 7 62 7 57 8 61 8 58 9 60 9 60 10 60 10 64 11 59 11 67 午夜 58 正午 68 5. P246 1 編寫Matlab程序，生成并繪制組合貝塞爾曲線。利用該程序生成和繪制過3

14、個控制點集{（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)}的貝塞爾曲線。 4、 實驗結果及分析 1.P202 1 實驗描述： 由題意可知，此題需要用最小二乘法進行計算，因為已知函數(shù)的5個插點并且知道它們的x,y的值。且函數(shù)的表達式為F=kx，所以只需用方程中便可計算出k的數(shù)值。 定義一個動態(tài)數(shù)組,用來依次存取x和y的插值。其中x,y的插值通過鍵盤手動輸入并賦予給a中的元素。然后通過相應的求和和基本運算便可以得到相應插值下的最小二乘法的函數(shù)

15、表達式。（其中k保留小數(shù)點后4位） 具體函數(shù)運行效果如下： （a）請在此輸入x的各插值 0.2 0.4 0.6 0.8 1.0 請在此輸入y的各插值 3.6 7.3 10.9 14.5 8.2 此函數(shù)的最小二乘法曲線表達式為 Y=14.8333*x 請按任意鍵繼續(xù)。。。 （b）請在此輸入x的各插值 0.2 0.4 0.6 0.8 1.0 請在此輸入y的各插值 5.3 10.6 15.9 21.2 26.4 此函數(shù)的最小二乘法曲線表達式為 Y=

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

17、的最小值，此時便可求出此5個參變量的值.然后繼續(xù)通過Matlab,將得到的公式以及各點，用plot函數(shù)，便可以求得。 實驗結果： 運行matlab結果如下： >> fminsearch(fiveOne,[1 1 1 1 1]) ans = 15.7225 1.3717 15.5359 1.2768 60.3579 此時的所求值便為函數(shù)的待定系數(shù)。所以可得最小二乘曲線的表達式為： 然后進行相應的繪制圖形便可以求出所要求出的結果。 結果分析： 通過最小二乘法多項式同樣可以繪制出三角函數(shù)的曲線。并且程序通用性高，只要輸入相應的插值便

18、可以計算出結果，結果系數(shù)的小數(shù)點有效位同時也可以自己確定。 3.P229 1 實驗描述： 由題意可知，由于這里涉及到了樣條線的運算，計算較為復雜。且要涉及到畫圖的部分，所以此處采用matlab計算較為方便快捷。 而書本上給出了一個用來計算三次緊壓樣條曲線的可調用函數(shù)，現(xiàn)在將其引用，并根據(jù)已知點得出相應的曲線。 實驗結果： 代入后得出的結果如下所示： >> X=[0 2 4 6 8]; >> Y=[0 40 160 300 480]; >> S=csfit(X,Y,0,98) S = 0.8125 8.3750 0

19、 0 -2.4375 13.2500 43.2500 40.0000 1.4375 -1.3750 67.0000 160.0000 -0.8125 7.2500 78.7500 300.0000 由結果可知此插值為在區(qū)間[0,8]中有三個分別劃分了[0,2],[2,4],[4,6],[6,8]四個區(qū)間的插點。 且多項式分別為 在matlab中通過polyval作出相應的曲線，再用

20、plot函數(shù)便可繪制圖線。繪制后圖線如下： 此時便可以直觀的看到一個平滑的樣條線。 結果分析： 通過給定的一些數(shù)據(jù)點，便可以繪制出過這些點的相應的樣條線，通過觀察能發(fā)現(xiàn)樣條線的平滑度與你選擇的樣條線類型以及數(shù)據(jù)點的分布有一定關聯(lián)。不僅僅是緊壓樣條線，相關其它的線也可以用同樣的方法一一給出。 4.P238 5 實驗描述： 由題意可知，題目是想求出所繪制數(shù)據(jù)點的三角多項式的逼近，三角多項式逼近的公式為： 此外，x為離散傅里葉級數(shù)，滿足條件： 實驗結果： 而所給的點x為[1,24]的離散數(shù)值點，所以無法直接對其作

21、出逼近公式，需要進行尺度變換，將x點轉換為: （1）通過matlab繪制出相關的三角多項式曲線，然后同樣通過matlab的繪制點，將點繪制到這個曲線之中，具體的matlab代碼如下： >>hold on; %保持圖形不動，繪制新的圖形入曲線中 >>X=[1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24]; %數(shù)據(jù)點的x的取值 >>Y=[58 58 58 58 57 57 57 58 60 64 67 68 66 66 65 64 6

22、3 63 62 61 60 60 59 58]; %數(shù)據(jù)點的y的取值 >>plot(X,Y,’xk’) %繪制出數(shù)據(jù)點 （2）然后便可以畫出如圖所示的插值數(shù)據(jù)點。 結果分析： 三角多項式的曲線擬合度非常高，能很好的繪制出圖像的具體形式而且曲線平滑，但是它需要滿足x屬于-pi到pi的區(qū)間內。 5.P246 1 實驗描述： 由題意可知，首先以3階貝塞爾函數(shù)為帶求函數(shù)。其求解格式如下： 其中b為伯恩斯坦多項式，其定義如下： 用matlab先設置一個參數(shù)t,然后再根據(jù)公式，通

23、過所給的數(shù)據(jù)點將伯恩斯坦多項式，以及x,y的關于參數(shù)t的多項式求出來。然后再把t設置為精度足夠大的單位序列，繪制圖線即可得出貝塞爾的效果。 實驗結果： 其matlab運行后結果如下：【以第一組控制點為例】 >> X=[0 1 1 3] X = 0 1 1 3 >> Y=[0 2 1 0] Y = 0 2 1 0 >> fiveTwo(X,Y) x = 3*t*(t - 1)^2 - 3*t^2*(t - 1) + 3*t^3 y = 6*t*(t - 1)^2 - 3*t^2*

24、(t - 1) 且繪制出第一組控制點所在位置以及三階貝塞爾曲線如下所示： 同理，第2,3組控制點所作的圖形如下所示： 結果分析： 可以通過Matlab函數(shù)繪制出三階的貝塞爾函數(shù)。只要給出控制點，便可自動繪制出所控制的三階貝塞爾函數(shù)以及控制點的位標。由此可以觀察到貝塞爾與控制點的約束作用，以及所求得貝塞爾函數(shù)是個相對平滑的曲線。 5、 實驗結論 曲線擬合對于實際的工程以及理論推導問題都有著重大的作用。在具體的實驗，數(shù)據(jù)分析上，往往我們只有巨量的已知離散點，想從離散點中得出函數(shù)表達式則就需要曲線擬合進行，根據(jù)

25、不同的要求，我們可以選擇最小二乘法的冪函數(shù)擬合或者是一次函數(shù)，二次函數(shù)擬合，抑或是精度非常高的傅里葉變換的三角函數(shù)擬合。同時在建模方面，貝塞爾函數(shù)的引用也從數(shù)學層面解決了如何用計算機繪制出光滑圓潤的曲線，在一些設計軟件中，例如3dmax和maya的三維建模，Auto CAD的工程建模，貝塞爾運算對于曲線曲面的設計有著舉足輕重的作用。 附件（代碼）: 1.P202 1 #include #include #include void main() { extern float

26、 afunction(float b[]); extern float bfunction(float z,float v); float *a,y,q,c; int i; a=new float[5]; //定義動態(tài)數(shù)組a[] printf("請在此輸入x的各插值\n"); scanf("%f %f %f %f %f",&a[0],&a[1],&a[2],&a[3],&a[4]); //手動輸入x各值 q=afunction(a);

27、 //對x值進行求和 a=new float[5]; //重新為a[]進行數(shù)據(jù)分配 printf("請在此輸入y的各插值\n"); scanf("%f %f %f %f %f",&a[0],&a[1],&a[2],&a[3],&a[4]); //手動輸入y各值 c=afunction(a); //對y值進行求和 y=bfunction(q,c); //計算出k的系數(shù)的大小 printf("此數(shù)據(jù)的最小二乘法曲線表達式為\n

28、 Y=%f*x\n",y); return; } float afunction(float b[]) { float x=0; int i; for(i=0;i<5;i++) x=b[i]+x; return x; } float bfunction(float z,float v) { float x; x=v/z; return x; } 2.P215 1 function z=fiveOne(u,y) A=u(1); B=u(2); C=u(3); D=u(4

29、); E=u(5); Z=0; for i=1:1:24 z=(A.*cos(i*B)+C.*sin(i*D)+E-y(i)).^2+z; %函數(shù)的線性最小二乘法的多項式 end 之后在Matlab的命令窗口輸入如下語句： >>fminsearch(‘fiveOne’,[1 1 1 1 1]) （繪制圖形方程組） >> x1=[0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24]; %數(shù)值點的x取值 >> y1=[58 58 58 58 57 57 57 58

30、60 64 67 68 66 66 65 64 63 63 62 61 60 60 59 58]; %數(shù)值點的x取值 >> plot(x1,y1,x) %繪制出24個數(shù)值點的圖形 >> x=linspace(0,25,100); >> hold on >> plot(x,y,k) %繪制出此函數(shù)的二乘法后的函數(shù)曲線 3.P229 1 （三次緊壓樣條線計算函數(shù)） function S=csfit(x,y,dx0,dxn) N=length(x)-1; H=diff(x); D=d

31、iff(y)./H; %d(k)的值 A=H(2:N-1); B=2*(H(1:N-1)+H(2:N)); C=H(2:N); U=6*diff(D); %U(k)的值 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)); for k=2:N-1 temp=A(k-1)/B(k-1); B(k)=B

32、(k)-temp*C(k-1); U(k)=U(k)-temp*U(k-1); end M(N)=U(N-1)/B(N-1); for k=N-2:-1:1 M(k+1)=(U(k)-C(k)*M(k+2))/B(k); end M(1)=3*(D(1)-dx0)/H(1)-M(2)/2; %三次緊壓約束中S"(0)的值 M(N+1)=3*(dxn-D(N))/H(N)-M(N)/2; %三次緊壓約束中S"(k)的值 for k=0:N-1 S(k+1,1)=(M(k+2)-M(k+1))/

33、(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(1)); %第一段樣條線 >> x2=2:.01:4; y

34、2=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 void main() { exte

35、rn float afunction(float X[],float Y[],int temp); extern float bfunction(float X[],float Y[],int temp); int i; float a[24]={1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24}; float b[24]={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; fl

36、oat *d; c=new float[7],d=new float[7]; for(i=0;i<=24;i++) { c[i]=afunction(a,b,i); d[i]=bfunction(a,b,i); } printf("由所給的數(shù)據(jù)點可以得出\n三角多項式T7(x)中a的系數(shù)分別為\n"); printf("%f %f %f %f %f %f %f\n",c[0],c[1],c[2],c[3],c[4],c[5],c[6]); printf("由所給的數(shù)據(jù)點可以得出\n三角多項式T7(x)中b的系數(shù)分別為\n"); printf("%f

37、%f %f %f %f %f %f\n",d[0],d[1],d[2],d[3],d[4],d[5],d[6]); system("pause"); return; } float afunction(float X[],float Y[],int temp) { float A=0; int i; for(i=0;i<=24;i++) A=Y[i]*cos(temp*X[i])/6+A; return A; } float bfunction(float X[],float Y[],int temp) { float A=0;

38、int i; for(i=0;i<=24;i++) A=Y[i]*sin(temp*X[i])/6+A; return A; } 4.P238 5 function []= fiveTwo() aa=0; bb=0; y=0; syms x for i=1:24 X(i)=-3.1415926+2.*i*3.1415926/24; end Y=[58 58 58 58 57 57 57 58 60 64 67 68 66 66 65 64 63 63 62 61 60 60 59 58]; for i

39、=1:24 for j=1:24 aa=Y(j).*cos(j.*X(i))+aa; end a(i)=12*aa; %求出多項式aj; end for i=1:24 for j=1:24 bb=Y(j).*sin(j.*X(i))+aa; end b(i)=12*bb; %求出多項式bj; end for i=1:24 y=a(i)*cos(i.*x)+b(i)*

40、sin(i.*x)+y; %求出三角多項式逼近 end x=-3.1415926:0.01:3.1415926; plot(x,y,g); %繪制圖形 5.P246 1 function z=fiveTwo(X,Y) syms t %定義t為參數(shù)變量 a=0; b=0; for i=0:1:3 B(i+1)=nchoosek(3,i)*t^i*(1-t)^(3-i); %計算伯恩斯坦多項式 end for i=1:1:4 a=X(i)

41、.*B(i)+a; %計算三階貝塞爾函數(shù)下x關于t的參數(shù)方程 b=Y(i).*B(i)+b; %計算三階貝塞爾函數(shù)下x關于t的參數(shù)方程 end t=0:.01:1; x=a; y=b; plot(X,Y,x) %在圖形中繪制出已知的4個控制點 X %輸出三階貝塞爾函數(shù)下x關于t的參數(shù)方程 Y %輸出三階貝塞爾函數(shù)下x關于t的參數(shù)方程 為了繪制出最終的貝塞爾函數(shù)，需要在matlab命令面板中輸入如下指令: >>x=X(轉換成元素的點乘除形式) >>y=Y(轉換成元素的點乘除形式) >>t=0:0.01:1; >>hold on %保持在一張圖內繪制 >>plot(x,y) %在圖形中繪制出所求得貝塞爾函數(shù)曲線