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

2、值、三次樣條插值、端點約束。關 鍵 詞 曲線擬合、最小二乘法擬合、冪函數擬合、對數函數擬合、線性插值、三次樣條插值、端點約束一、實驗目的1. 掌握曲線擬合方式及其常用函數指數函數、冪函數、對數函數的擬合。2. 掌握最小二乘法、線性插值、三次樣條插值、端點約束等。3. 掌握實現曲線擬合的編程技巧。二、實驗原理1. 曲線擬合曲線擬合是平面上離散點組所表示的坐標之間的函數關系的一種數據處理方法。用解析表達式逼近離散數據的一種方法。在科學實驗或社會活動中，通過實驗或觀測得到量x與y的一組數據對（Xi，Yi)（i=1，2，.m），其中各Xi是彼此不同的 。人們希望用一類與數據的背景材料規律相適應的解析表

3、達式，y=f(x，c）來反映量x與y之間的依賴關系，即在一定意義下“最佳”地逼近或擬合已知數據。f(x，c)常稱作擬合模型 ，式中c=(c1，c2，cn)是一些待定參數。當c在f中線性出現時，稱為線性模型，否則稱為非線性模型。有許多衡量擬合優度的標準，最常用的一種做法是選擇參數c使得擬合模型與實際觀測值在各點的殘差(或離差)的加權平方和達到最小，此時所求曲線稱作在加權最小二乘意義下對數據的擬合曲線。有許多求解擬合曲線的成功方法，對于線性模型一般通過建立和求解方程組來確定參數，從而求得擬合曲線。至于非線性模型，則要借助求解非線性方程組或用最優化方法求得所需參數才能得到擬合曲線，有時稱之為非線性最

4、小二乘擬合。曲線擬合：貝塞爾曲線與路徑轉化時的誤差。值越大，誤差越大；值越小，越精確。2. 最小二乘法擬合：最小二乘法（又稱最小平方法）是一種數學優化技術。它通過最小化誤差的平方和尋找數據的最佳函數匹配。利用最小二乘法可以簡便地求得未知的數據，并使得這些求得的數據與實際數據之間誤差的平方和為最小。最小二乘法還可用于曲線擬合。其他一些優化問題也可通過最小化能量或最大化熵用最小二乘法來表達。 函數曲線為： y=Ax+B 其中系數滿足下列的正規方程： 3. 冪函數擬合：函數曲線為： 設有N個點，其中橫坐標是確定的。最小二乘冪函數擬合曲線的系數A為： 、 4. 對數函數擬合：對數函數（lograrit

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

6、f(x) = 1 / (1 + x2)如果在區間-5，5上取7個等距節點:xk=5*(k/3-1) (k=0,1,2,.,6),由lagrange插值公式可得到f(x)的次L7(x)。如圖所示：L7(x)僅在區間的中部能較好的逼近函數f(x), 在其它部位差異較大，而且越接近端點，逼近效果越差。可以證明，當節點無限加密時，Ln(x)也只能在很小的范圍內收斂，這一現象稱為Runge現象。它表明通過增加節點來提高逼近程度是不適宜的，因而不采用高次多項式插值。如果我們把以上的點用直線連接起來，顯然比L7(x)要更逼近f(x)。這就是分段線性插值。而在實際應用中通常采用分段低次插值來提高近似程度，比如

7、可用分段線性插值或分段三次埃爾米特插值來逼近已知函數，但它們的總體光滑性較差。為了克服這一缺點，一種全局化的分段插值方法三次樣條插值成為比較理想的工具。6. 三次樣條插值：設有N+1個點，其中。如果存在N個三次多項式，系數為滿足如下性質： 則成函數S(x)為三次樣條函數。7. 端點約束：緊壓樣條：存在唯一的三次樣條曲線，其一階導數的邊界條件是： natural樣條：存在唯一的三次樣條曲線，它的自由邊界條件是： 外推樣條：存在唯一的三次樣條曲線，其中通過對點x1和x2進行外推得到 ，同時通過對點X(n-1)和X(N-2)進行外推得到。端點曲率調整：存在唯一的三次樣條曲線，其中二階導數的邊界條件和

8、是確定的。拋物線終結樣條：存在唯一的三次樣條曲線，其中二階在區間X0，X1內，而在Xn-1,Xn內。3、 實驗內容1. P202 1胡克定律指出F=kx，其中F是拉伸彈簧的拉力（單位為盎司），x為拉伸的長度（單位為英寸）。根據下列試驗數據，求解拉伸常量k的近似值。xkFk0.23.60.47.30.610.90.814.51.018.2XkFk0.25.30.410.60.615.90.821.21.026.4 （a） （b） 2. P215 1洛杉磯（美國城市）郊區11月8日的溫度記錄入下表所示，其中共有24個數據點。（a）根據例5.5中的處理過程（使用fmins命令），對給定的數據求解最小

9、二乘曲線。（b）求 （c）在同一坐標系下畫出這些點集和(a)得出的最小二乘曲線。時間，p.m.溫度時間，a.m.溫度1661582662583653584644585635576636577627578618589609601060106411591167午夜58正午683. P229 1一個轎車在時間Tk時經過的距離dk，如下表所示。使用程序5.3，并根據一階導數邊界條件,求這些數據的三次緊壓樣條插值。 時間，tk02468距離，dk0401603004804. P238 5美國洛杉磯郊區11月9日的溫度（華氏溫度）如表5.10所示。采用24小時制。 （a）求三角多項式 （b）在同一坐標系下

10、，畫出圖和24個數據點。 （c）使用本地的溫度情況重新求解問題(a)和問題(b)。時間,p.m溫度時間,a.m溫度1661582662583653584644585635576636577627578618589609601060106411591167午夜58正午68 5. P246 1編寫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)的貝塞爾曲線。4、 實驗結果及分析1.P202 1實驗描述：由題意

11、可知，此題需要用最小二乘法進行計算，因為已知函數的5個插點并且知道它們的x,y的值。且函數的表達式為F=kx，所以只需用方程中便可計算出k的數值。定義一個動態數組,用來依次存取x和y的插值。其中x,y的插值通過鍵盤手動輸入并賦予給a中的元素。然后通過相應的求和和基本運算便可以得到相應插值下的最小二乘法的函數表達式。（其中k保留小數點后4位）具體函數運行效果如下：（a）請在此輸入x的各插值 0.2 0.4 0.6 0.8 1.0 請在此輸入y的各插值 3.6 7.3 10.9 14.5 8.2 此函數的最小二乘法曲線表達式為 Y=14.8333*x 請按任意鍵繼續。（b）請在此輸入x的各插值 0

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

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

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

15、75 -1.3750 67.0000 160.0000 -0.8125 7.2500 78.7500 300.0000由結果可知此插值為在區間0,8中有三個分別劃分了0,2,2,4,4,6,6,8四個區間的插點。 且多項式分別為 在matlab中通過polyval作出相應的曲線，再用plot函數便可繪制圖線。繪制后圖線如下：此時便可以直觀的看到一個平滑的樣條線。結果分析：通過給定的一些數據點，便可以繪制出過這些點的相應的樣條線，通過觀察能發現樣條線的平滑度與你選擇的樣條線類型以及數據點的分布有一定關聯。不僅僅是緊壓樣條線，相關其它的線也可以用同樣的方法一一給出。4.P238 5實驗描述：由題意

16、可知，題目是想求出所繪制數據點的三角多項式的逼近，三角多項式逼近的公式為： 此外，x為離散傅里葉級數，滿足條件： 實驗結果：而所給的點x為1,24的離散數值點，所以無法直接對其作出逼近公式，需要進行尺度變換，將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; %數據點的x的取值&

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

18、恩斯坦多項式，以及x,y的關于參數t的多項式求出來。然后再把t設置為精度足夠大的單位序列，繪制圖線即可得出貝塞爾的效果。實驗結果：其matlab運行后結果如下：【以第一組控制點為例】>> X=0 1 1 3X = 0 1 1 3>> Y=0 2 1 0Y = 0 2 1 0>> fiveTwo(X,Y)x = 3*t*(t - 1)2 - 3*t2*(t - 1) + 3*t3 y = 6*t*(t - 1)2 - 3*t2*(t - 1) 且繪制出第一組控制點所在位置以及三階貝塞爾曲線如下所示：同理，第2,3組控制點所作的圖形如下所示： 結果分析：可以通過

19、Matlab函數繪制出三階的貝塞爾函數。只要給出控制點，便可自動繪制出所控制的三階貝塞爾函數以及控制點的位標。由此可以觀察到貝塞爾與控制點的約束作用，以及所求得貝塞爾函數是個相對平滑的曲線。5、 實驗結論曲線擬合對于實際的工程以及理論推導問題都有著重大的作用。在具體的實驗，數據分析上，往往我們只有巨量的已知離散點，想從離散點中得出函數表達式則就需要曲線擬合進行，根據不同的要求，我們可以選擇最小二乘法的冪函數擬合或者是一次函數，二次函數擬合，抑或是精度非常高的傅里葉變換的三角函數擬合。同時在建模方面，貝塞爾函數的引用也從數學層面解決了如何用計算機繪制出光滑圓潤的曲線，在一些設計軟件中，例如3dm

20、ax和maya的三維建模，Auto CAD的工程建模，貝塞爾運算對于曲線曲面的設計有著舉足輕重的作用。附件（代碼）:1.P202 1#include<stdio.h>#include<math.h>#include<stdlib.h>void main() extern float afunction(float b); extern float bfunction(float z,float v); float *a,y,q,c; int i; a=new float5; /定義動態數組a printf("請在此輸入x的各插值n");

21、scanf("%f %f %f %f %f",&a0,&a1,&a2,&a3,&a4); /手動輸入x各值 q=afunction(a); /對x值進行求和 a=new float5;/重新為a進行數據分配 printf("請在此輸入y的各插值n"); scanf("%f %f %f %f %f",&a0,&a1,&a2,&a3,&a4); /手動輸入y各值 c=afunction(a);/對y值進行求和 y=bfunction(q,c); /計算出k的系

22、數的大小 printf("此數據的最小二乘法曲線表達式為n Y=%f*xn",y); return;float afunction(float b)float x=0;int i;for(i=0;i<5;i+)x=bi+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);E=u(5);Z=0;for i=1:1:24z=(A.*cos(i*B)+C.*sin(i*D)+E

23、-y(i).2+z;%函數的線性最小二乘法的多項式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; %數值點的x取值>> y1=58 58 58 58 57 57 57 58 60 64 67 68 66 66 65 64 63 63 62 61 60 60 59 58; %數值點的x取值>> plot(x1,y1,'

24、;x')%繪制出24個數值點的圖形>> x=linspace(0,25,100);>> hold on>> plot(x,y,'k')%繪制出此函數的二乘法后的函數曲線3.P229 1 （三次緊壓樣條線計算函數）function S=csfit(x,y,dx0,dxn)N=length(x)-1;H=diff(x);D=diff(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*

25、(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(k)-temp*C(k-1); U(k)=U(k)-temp*U(k-1);endM(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);endM(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&q

26、uot;(k)的值for k=0:N-1 S(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（曲線繪制函數程序）>> x1=0:.01:2; y1=polyval(S(1,:),x1-X(1); %第一段樣條線>> x2=2:.01:4; y2=polyval(S(2,:),x2-X(2);%第二段樣條線>

27、> 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')%繪制連接成完整曲線，并且標注數據點的位置。 #include<stdio.h>#include<stdlib.h>#include<math.h>void main()extern float afunction(float X,float Y

28、,int temp);extern float bfunction(float X,float Y,int temp);int i;float a24=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 b24=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=new float7,d=new float7;for(i=0;i<=24;i+)ci=afunction(a,

29、b,i);di=bfunction(a,b,i);printf("由所給的數據點可以得出n三角多項式T7(x)中a的系數分別為n");printf("%f %f %f %f %f %f %fn",c0,c1,c2,c3,c4,c5,c6);printf("由所給的數據點可以得出n三角多項式T7(x)中b的系數分別為n");printf("%f %f %f %f %f %f %fn",d0,d1,d2,d3,d4,d5,d6);system("pause");return;float afunction(float X,float Y,int temp)float A=0;int i;for(i=0;i<=24;i+)A=Yi*cos(temp*Xi)/6+A; return A;float bfunction(float X,f