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使得擬合模型與實際觀測值在各點的殘差（或離差）ek yH f, c）的加權平方和達到最小，此時所求曲線稱作在加權最小二乘意義下對數據的擬合曲線。 有許多求解擬合曲線的成功方法，對于線 性模型一般通過建立和求解方程組來確定參數，從而求得擬合曲線。至于非線性 模型，則要借助求解非線性方程組或用最優化方法求得

4、所需參數才能得到擬合曲 線，有時稱之為非線性最小二乘擬合。曲線擬合：貝塞爾曲線與路徑轉化時的誤差。值越大，誤差越大；值越小， 越精確。2. 最小二乘法擬合:最小二乘法（又稱最小平方法）是一種數學優化技術。它通過最小化誤差的 平方和尋找數據的最佳函數匹配。利用最小二乘法可以簡便地求得未知的數據， 并使得這些求得的數據與實際數據之間誤差的平方和為最小。最小二乘法還可用于曲線擬合。其他一些優化問題也可通過最小化能量或最大化熵用最小二乘法來 表達。函數曲線為：y=Ax+B其中系數滿足下列的正規方程：NX：k 1ANXkk 1NBXkYkk 1NNXk aNByk3.幕函數擬合：k 1k 1函數曲線為:

5、設xk,yk N1有N個點，其中橫坐標是確定的。最小二乘幕函數擬合曲線的系數A為：NNM2M、A (xk yk)/(xk )k 1k 1、4. 對數函數擬合:對數函數(lograrithmic function )的標準式形式為Y a bl nX(X 0)b>0時，丫隨X增大而增大，先快后慢；b<0時，丫隨X增大而減少，先快 后慢，見圖12.4(c)、(d)。當以丫和InX繪制的散點圖呈直線趨勢時，可考慮采 用對數函數描述丫與X之間的非線性關系，式中的b和a分別為斜率和截距。更一般的對數函數Y=a+bln (X+k)式中k為一常量，往往未知。(a) lnY=l na+bX(b) l

6、nY=l na-bX(c) Y=a+blnX(d) Y=a-blnX5. 線性插值:在代數插值中，為了提高插值多項式對函數的逼近程度一般是增加節點的個 數，即提高多項式的次數，但這樣做往往不能達到預想的效果。如下圖所示：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),在其它部位差異較大，而且越接近端點，逼 近效果越差?？梢宰C明，當節點無限加密時，Ln(x)也只能在很小的范圍內收斂， 這一現象稱為R

7、unge現象。它表明通過增加節點來提高逼近程度是不適宜的， 因而不采用高次多項式插值。如果我們把以上的點用直線連接起來，顯然比L7(x)要更逼近f(x)。這就是分段線性插值。而在實際應用中通常采用分段低次插值來提高近似程度，比如可用分段線性插值或分段三次埃爾米特插值來逼近已知函數，但它們的總體光滑性較差。為了克服這一缺點，一種全局化的分段插值方法一一三次樣條插值成為比較理想的工 具。6. 三次樣條插值:設 xk, yk N0有N+1個點，其中ax0x x2. xNb。如果存在 N個三次多項式Sk(x),系數為SkJ0，Sk，Skj2，Skj3滿足如下性質:S(x)Sk ( x)Sk ,0Sk

8、,1(xXk)Sk,2(XXk)2Sk,:3(Xx Xk,xk1,k0,1,.,N1S(xJykk0,1,.NS(xk 1)Sk 1 (xk 1)k0,1,.N2S'(Xk1 ) S'k 1 (xk 1)k0,1,.N2S''(xk1 )S''k 1 (xk 1)k0,1,.N2則成函數S(x)為三次樣條函數。Xk)'7. 端點約束：緊壓樣條：存在唯一的三次樣條曲線，其一階導數的邊界條件是：S'(a) d°,S(b) dNnatural樣條：存在唯一的三次樣條曲線，它的自由邊界條件是：S''(a)0,S&

9、quot;(b)0外推樣條：存在唯一的三次樣條曲線，其中通過對點x1和x2進行外推得到S a，同時通過對點X(n-1)和X(N-2)進行外推得到S b。端點曲率調整：存在唯一的三次樣條曲線，其中二階導數的邊界條件S a和S b是確定的。拋物線終結樣條：存在唯一的三次樣條曲線，其中二階在區間X0，X1內S (x)0,而在Xn-1,Xn內S (x)0。三、實驗內容1. P202 1胡克定律指出F=kx，其中F是拉伸彈簧的拉力（單位為盎司），x為拉伸 的長度（單位為英寸）。根據下列試驗數據，求解拉伸常量k的近似值。(a)XkFk0.23.6 10.47.30.610.90.814.5 T1.018.

10、2(b)XkFk0.2 ”5.30.410.60.615.90.8 n21.21.026.42. P215 1洛杉磯（美國城市）郊區11月8日的溫度記錄入下表所示，其中共有 24個 數據點。（a）根據例5.5中的處理過程（使用fmins命令），對給定的數據求解最小 二乘曲線f（x） Acos（Bx） Csin（Dx） E。（b）求 E2（f）。（c）在同一坐標系下畫出這些點集和（a）得出的最小二乘曲線。時間，p.m.溫度時間，a.m.溫度166158266258365358464458563557663657762757:86185896096010601064:11591167午夜58正午6

11、83. P229 1一個轎車在時間Tk時經過的距離dk,如下表所示。使用程序5.3,并根據 階導數邊界條件S' （0）0, S'（8）98，求這些數據的三次緊壓樣條插值。時間，tk02468距離，dk0401603004804. P238 5美國洛杉磯郊區11月9日的溫度（華氏溫度）如表 5.10所示。采用24小 時制。（a）求三角多項式T7 （x）（b）在同一坐標系下，畫出圖T7（x）和24個數據點。（c）使用本地的溫度情況重新求解問題（a）和問題（b）。時間,p.m溫度時間,a.m溫度166158266258365358464458563557663657762757861

12、8589609601060106411591167午夜58正午685. 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）的貝塞爾曲線。四、實驗結果及分析1. P202 1實驗描述：由題意可知，此題需要用最小二乘法進行計算，因為已知函數的5個插點并 且知道它們的x,y的值。且函數的表達式為F=kx，所以只需用方程中NN(Xk)Ayk便可計算出k的數值。k 1k 1定

13、義一個動態數組a ，用來依次存取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.2 0.4 0.6 0.8 1.0請在此輸入y的各插值5.3 10.6 15.9 21.2 26.4 此函數的最小二乘法曲線表達式為Y=26.466

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

15、續通過Matlab, 將得到的公式以及各點，用plot函數，便可以求得。實驗結果：運行 matlab結果如下：>> fmin search('five On e',1 1 1 1 1)ans =15.72251.371715.53591.276860.3579此時的所求值便為函數的待定系數。所以可得最小二乘曲線的表達式為：f (x)15.7225cos(1.3717x) 15.5359sin(1.2768) 60.3579然后進行相應的繪制圖形便可以求出所要求出的結果。結果分析：通過最小二乘法多項式同樣可以繪制出三角函數的曲線。并且程序通用性 高，只要輸入相應的插值

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

17、001.4375-1.375067.0000160.0000-0.81257.250078.7500300.0000由結果可知此插值為在區間0,8中有三個分別劃分了 0,2,2,4,4,6,6,8四 個區間的插點。且多項式分別為S°(x)0.8125x38.375x2S(x)32.4375(x 2)13.25(x2)243.25(x2) 40S2(x)1.4375(x 4)3 1.375(x4)267(x 4)160S3(x)0.8125(x 6)37.25(x6)278.75(x6) 300在matlab中通過polyval作出相應的曲線，再用plot函數便可繪制圖線制后圖線如下:

18、此時便可以直觀的看到一個結果分析:通過給定的一些數據點，便可以繪制出過這些點的相應的樣條線，通過觀察 能發現樣條線的平滑度與你選擇的樣條線類型以及數據點的分布有一定關聯。 不 僅僅是緊壓樣條線，相關其它的線也可以用同樣的方法一一給出。4. P238 5實驗描述：由題意可知，題目是想求出所繪制數據點的三角多項式的逼近，三角多項式逼近的公式為：a0MTm(x) (ajcos(jx) bjSin(jx)2 j i此外，x為離散傅里葉級數，滿足條件：Xjj,j 0,1,NjN實驗結果：而所給的點X為1,24的離散數值點，所以無法直接對其作出逼近公式，需 要進行尺度變換，將X點轉換為：Xjj,j 0,1

19、,.,24j24(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的取值>>Y=58 58 58 58 57 57 57 58 60 64 67 68 66 66 65 64 63 63 62 61 60 60 5958;%數據點的y的取值>>plot(X,Y

20、,'xk')%繪制出數據點(2) 然后便可以畫出如圖所示的插值數據點。結果分析:三角多項式的曲線擬合度非常高，能很好的繪制出圖像的具體形式而且曲線 平滑，但是它需要滿足x屬于-pi到pi的區間內。5. P246 1實驗描述：由題意可知，首先以3階貝塞爾函數為帶求函數。其求解格式如下：P(t)Bi,3(t)其中b為伯恩斯坦多項式，其定義如下：3 i 3 iBi,3(t). ti(1 t)3ii用matlab先設置一個參數t,然后再根據公式，通過所給的數據點將伯恩斯坦 多項式，以及x,y的關于參數t的多項式求出來。然后再把t設置為精度足夠大 的單位序列，繪制圖線即可得出貝塞爾的效果

21、。實驗結果：其matlab運行后結果如下：【以第一組控制點為例】>>X=0 11 3X =01 13>>Y=0 21 0Y =02 10>> fiveTwo(X,Y)x = 3*t*(t - 1)A2 - 3*tA2*(t - 1) + 3*tA36*t*(t - 1)A2 - 3*tA2*(t - 1)且繪制出第一組控制點所在位置以及三階貝塞爾曲線如下所示:同理，第2,3組控制點所作的圖形如下所示:結果分析:可以通過Matlab函數繪制出三階的貝塞爾函數。只要給出控制點，便可自 動繪制出所控制的三階貝塞爾函數以及控制點的位標。由此可以觀察到貝塞爾與控制點的

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

23、 1#i nclude<stdio.h>#in clude<math.h>#i nclude<stdlib.h>void mai n()extern float afunction(float b);exter n float bfun cti on( float z,float v);float *a,y,q,c;int i;a=new float5;/定義動態數組 aprintf("請在此輸入x的各插值n");scan f("%f %f %f %f %f",&a0,&a1,&a2,&

24、a3,&a4);/ 手動輸入x各值q=afunction(a);對 x 值進行求和a=new float5;/重新為a進行數據分配printf("請在此輸入y的各插值n");scan f("%f %f %f %f %f",&a0,&a1,&a2,&a3,&a4);/ 手動輸入y各值c=afunction(a);對y值進行求和y=bfunction(q,c);計算出 k的系數的大小printf("此數據的最小二乘法曲線表達式為n Y=%f*xn",y);return; float afun

25、ction( 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 fun cti on z=five On e(u,y)A=u(1);B=u( 2);C=u(3);D=u ；E=u (5);Z=0;for i=1:1:24z=(A.*cos(i*B)+C.*s in (i*D)+E-y(i)A2+z;% 函數的線性最小二乘法的多項式end之后在Matlab的命令窗口輸入如下語句：>

26、>fmi nsearch( five On e'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取值>> 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; %數值點 的x取值>> plot(x1,y1,'x')>> x=li nspace(0,25,100);>>

27、hold on>> plot(x,y,'k')%繪制出24個數值點的圖形%繪制出此函數的二乘法后的函數曲線3. P229 1(三次緊壓樣條線計算函數) fun cti on S=csfit(x,y,dxO,dx n)N=le ngth(x)-1;H=diff(x);D=diff(y)./H;A=H(2:N-1);%d(k)的值B=2*(H(1:N-1)+H(2:N);C=H(2:N);U=6*diff(D);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*(

28、dx n-D(N);for k=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);for k=N-2:-1:1M(k+1)=(U(k)-C(k)*M(k+2)/B(k); endM(1)=3*(D(1)-dxO)/H(1)-M(2) /2; M(N+1)=3*(dx n-D(N)/H(N)-M(N)/2; 值%U(k)的值%三次緊壓約束中S"(0)的值%三次緊壓約束中S"(k)的for k=0:N-1S(k+1,1)=(M(k+2)-M(k+1

29、)/(6*H(k+1); S(k+1,2)=M(k+1)/2;S(k+1,3)=D(k+1)-H(k+1)*(2*M(k+1)+M(k+2)/6; S(k+1,4)=y(k+1);End%S( k,3)的值%S( k,2)的值%S( k,1)的值 %S( k,0)的值(曲線繪制函數程序)>> x1=0:.01:2; y1=polyval(S(1,:),x1-X(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:.

30、01:8; y4=polyval(S(4,:),x4-X(4);>> plot(x1,y1,x2,y2,x3,y3,x4,y4,X,Y ,'x') 注數據點的位置。%第一段樣條線%第二段樣條線%第三段樣條線%第四段樣條線%繪制連接成完整曲線，并且標#i nclude<stdio.h>#i nclude<stdlib.h>#in clude<math.h>void mai n()exter n float afun cti on( float X,float Y,i nt temp);extern float bfunction(f

31、loat 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;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=new float7,d=new float7;for(i=0;i<=24;i+)ci=afu nctio n(a,b,i);di=bfu nctio n(a,b,i);printf("由所給的數據點可以得出n三角多項式T7(x)中a的系數分別為n");prin tf("%f %f %f %f %f %f %fn",c0,c1,c2,c3,c4,c5,c6);printf("由所給的數據點可以得出n三角多項式T7(x)中b的系數分別為n");prin tf("%f %f %f %f %f %f %fn",d0,d1,d2,d3,d4,d5,d6);system("pause");return;float afu