温度分布的曲线拟合

温度分布的曲线拟合1. 实验描述曲线拟合是指用连续曲线近似地刻画或比拟平面上离散点组所表示的坐标之间的函数关系。

更广泛地说，空间或高维空间中的相应问题亦属此范畴。

在数值分析中，曲线拟合就是用解析表达式逼近离散数据，即离散数据的公式化。

实践中，离散点组或数据往往是各种物理问题和统计问题有关量的多次观测值或实验值，它们是零散的，不仅不便于处理，而且通常不能确切和充分地体现出其固有的规律。

这种缺陷正可由适当的解析表达式来弥补。

2. 实验内容温度分布的曲线拟合，度数据采用下表中的数据：要求：1.2.曲线的最小二乘抛物线拟合； 3.三次样条插值拟合； 4.T7的三角多项式拟合。

5.有4个控制点的贝塞尔曲线拟合。

3. 实验结果及分析线性的最小二乘拟合:设 x k ，y k k =1N有N 个点，其中横坐标 x k k =1N 是确定的。

最小二乘拟合曲线y =Ax +B的系数是下列线性方程组的解，这些方程成为正规方程：x k 2Nk =1A + x k Nk =1B = x k y k Nk =1x k Nk =1A +NB = y k N k =1将x ，y 的值代入俩个方程解出系数A ，B 在画图可得下图如图所示红*号为原来实际数据，直线为所求式，可见其与实际情况差别甚大，不能反映实际的情况。

所以模拟效果不好。

曲线的最小二乘抛物线拟合:设 x k ，y k k =1N有N 个点，其中横坐标是确定的。

最小二乘的抛物线细数标示为y =f x =Ax 2+Bx +C求解A ，B 和C 的线性方程组为：x k 4Nk =1A + x k 3Nk =1B + x k 2Nk =1C = y k x k 2Nk =1510152025305959.56060.56161.56262.56363.564X(357 X)/2300 + 5445/92x k 3N k =1A + x k 2N k =1B + x k N k =1C = y k x k Nk =1x k 2N k =1 A + x k N k =1 B +N C = y k N k =1。

为求系数A 、B 、C 同理使(,,)E A B C 最小2(,,)()Nk k kk E A B C A B C y x x ==++-∑分别对(,,)E A B C 求偏导并使其等于零。

可得43221111()()()NNNNk k k kkk k k k A B yx x x x====++=∑∑∑∑321111()()()N N N Nk k k kkk k k k A B C yx x x x====++=∑∑∑∑2111()()NNNk k kk k k A B NC y x x ===++=∑∑∑带入数据得Y=-6713909909675507/144115188075855872*x^2+2909186376197643/2251799813685248*x+1923218375727969/35184372088832三次样条插值拟合:设 x k ，y k k =1N有N+1个点，其中a=x 0<x 1<…<x N =b 。

如果存在存在N 个三次多项式S k x ，系数为s k,0，s k,1，s k,2和s k,3，满足如下性质： ⅠS x =S k x =s k ,0+s k ,1 x −x k +s k ,2 x −x k 2+s k ,3 x −x k 3x ∈x k ,x k +1 ，k=0,1,…,N-1ⅡS x k =y k k=0,1,…,N ⅢS k x k +1 =S k +1 x k +1 k=0,1,…,N-2ⅣS k ′ x k +1 =S k +1′ x k +1 k=0,1,…,N-2 ⅤS k ′′ x k +1 =S k +1′′ x k +1 k=0,1,…,N-2则称函数S(X)为三次样条函数。

又可将S k x 写成如下形式： S k x =−m k k x k +1−x 3+−m k +1k x −x k 3+ y k k −m k ℎkx k +1−x + y k +1ℎk −m k +1ℎk 6x −x k且其中m k =S k ′′ x k ，m k +1=S k ′′ x k +1 ，ℎk =x k +1−x k ，y k =S k x k ，y k +1=S k x k +1 。

051015202552545658606264X- (135 X 2)/2392 + (93657 X)/59800 + 9765/184化简可得S k ′ x k =−m k 3ℎk −m k +16ℎk +d k （1）S k−1′x k =−m k 3ℎk−1−m k−16ℎk−1+d k−1（2）其中d k =y k +1−y kℎk。

利用上俩式及性质Ⅳ可得如下关系式：ℎk−1m k−1+2 ℎk−1+ℎk m k +ℎk m k +1=u k其中u k =6 d k −d k−1 ，k=1,2,…,N-1。

有图可得已知点都在所求曲线上，很好的反应了其的变化，但对以后的估计上，会存在问题，很难得到需要的结果。

三角多项式拟合:设2M+1<=N ，根据N 个等距离的横坐标值x k =−π+2πk N，k=1,2,…,N ，构造M 阶三角多项式，其形式为P x =a 02+ a j cos jx +b j sin jx Mj =1令2,0,1....j j j N Nx ππ=-+=其中 也同样用最小二乘法使下式最小得51015202512()cos(),0,1.....,Nj k k k f j j M N a x x ===∑其中 12()sin(),1,2.....,Nj k k k f j j M N b x x ===∑其中 带入所知数据，并将其横坐标该为得图原始数据图0510152025拟合图形图可见变化基本很模糊，不能反映真实情况，没有实际的意义。

贝塞尔曲线拟合:给定一个控点集， P i i =0N ，其中P i = x i ,y i ，定义P t = P i B i ,N t Ni =0为N 阶贝塞耳曲线，其中B i ,N t ,i=0,1…,N 是N 阶伯恩斯坦多项式，t ∈ 0,1 。

另外P(t)也可表示为P(t)=(x(t),y(t))，其中x t = x i B i ,N (t )N i =0，y t = y i B i ,N (t )N i =0运用3阶伯恩斯坦多项式为30,3(1)t B=-，21,33(1)t t B =-，22,33(1)t t B =-，33,3t B = 5.2可令,()0()Ni N t i i x t x B ==∑，,()()Ni N t ii y t yB==∑ 5.3将24个已知数据分成8个4个控制点的贝塞尔曲线拟合。

前一个点与后一个重合。

-6-4-20246-500-400-300-200-1000100200300400500x-2816/23 cos(x)-...+1467/230510152025可见逼近曲线与实际重合的很好，很好的反映了实际数据的情况。

效果很好。

4. 结论五种方法都描述了数据的变化，其中第三、第五种方法很好的反映了数据，且在不知道点能很好的估计出未知的数据。

第一第二种方法同属多项式逼近。

在数据成某种变化时，因为简单，常被采用。

但在本实验中效果不好。

附件（代码）线性的最小二乘拟合:clearx=[1:1:24];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];syms XYN=length(x);a=0;b=0;c=0;d=0;for k=1:Na=x(:,k)^2+a;b=x(:,k)+b;c=x(:,k)*y(:,k)+c;d=y(:,k)+d;endA=[a b;b 24];P=A\B;Y=P(1)*X+P(2);ezplot(Y,[-1,30])hold onplot(x,y,'r*')hold off曲线的最小二乘抛物线拟合:clearx=[1:1:24];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];syms XYh=x.^4;j=x.^3;k=x.^2a=sum(h);b=sum(j);c=sum(k);d=y*(k)';e=sum(x);f=y*(x)';g=sum(y);%Çó½â·½³ÌϵÊýA=[a b c;b c e;c e 24]B=[d; f; g]P=A\B;%½â³öϵÊýY=P(1)*X^2+P(2)*X+P(3)ezplot(Y,[-1,25])hold onplot(x,y,'g',x,y,'r*')%»-ͼhold off三次样条插值拟合:clearX=[0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 1718 19 20 21 22 23];Y=[58 58 58 58 58 57 57 57 58 60 64 67 68 66 6665 64 63 63 62 61 60 60 59];%ÒÑÖªµãdx0=(Y(2)-Y(1))/(X(2)-X(1));dxn=(Y(24)-Y(23))/(X(24)-X(23));%±ß½çÌõ¼þN=length(X)-1;H=diff(X);D=diff(Y)./H;A=H(2:N-1);B=2*(H(1:N-1)+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*(dxn-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)-dx0)/H(1)-M(2)/2;M(N+1)=3*(dxn-D(N))/H(N)-M(N)/2;for k=0:N-1S(k+1,1)=(M(k+2)-M(k+1))/(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%Çó½âÑùÌõϵÊýSx1=0:0.01:1;y1=polyval(S(1,:),x1-X(1));x2=1:0.01:2;y2=polyval(S(2,:),x2-X(2));x3=2:0.01:3;y3=polyval(S(3,:),x3-X(3));x4=3:0.01:4;y4=polyval(S(4,:),x4-X(4));x5=4:0.01:5;y5=polyval(S(5,:),x5-X(5));x6=5:0.01:6;y6=polyval(S(6,:),x6-X(6));x7=6:0.01:7;y7=polyval(S(7,:),x7-X(7));x8=7:0.01:8;y8=polyval(S(8,:),x8-X(8));x9=8:0.01:9;y9=polyval(S(9,:),x9-X(9));x10=9:0.01:10;y10=polyval(S(10,:),x10-X(10)) ;x11=10:0.01:11;y11=polyval(S(11,:),x11-X(11) );x12=11:0.01:12;y12=polyval(S(12,:),x12-X(12) );x13=12:0.01:13;y13=polyval(S(13,:),x13-X(13) );x14=13:0.01:14;y14=polyval(S(14,:),x14-X(14) );x15=14:0.01:15;y15=polyval(S(15,:),x15-X(15));x16=15:0.01:16;y16=polyval(S(16,:),x16-X(16));x17=16:0.01:17;y17=polyval(S(17,:),x17-X(17));x18=17:0.01:18;y18=polyval(S(18,:),x18-X(18));x19=18:0.01:19;y19=polyval(S(19,:),x19-X(19));x20=19:0.01:20;y20=polyval(S(20,:),x20-X(20));x21=20:0.01:21;y21=polyval(S(21,:),x21-X(21));x22=21:0.01:22;y22=polyval(S(22,:),x22-X(22));x23=22:0.01:23;y23=polyval(S(23,:),x23-X(23));%Çó½â¶ÔÓ¦µãplot(x1,y1,x2,y2,x3,y3,x4,y4,x5,y5,x6,y6,x7,y7,x8,y8,x9,y9,x10,y10,x11,y11,x12,y12,x13,y13,x14,y14,x15,y15,x16,y16,x17,y17,x18,y18,x19,y19,x20,y20,x21,y21,x22,y22,x23,y23)%»-ͼhold onplot(X,Y,'r*')hold off三角多项式拟合:clearsyms xyy=0;M=input('input M=')Y=[58 58 58 58 58 57 57 57 58 60 64 67 68 66 6665 64 63 63 62 61 60 60 59];X=[0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 1718 19 20 21 22 23];%ÒÑÖªµãN=length(X)-1;max1=fix((N-1)/2);if M>max1M=max1;end%×î¸ß½×A=zeros(1,M+1);B=zeros(1,M+1);for i=1:M+1for k=1:NA(i)=A(i)+Y(k)*cos(i*X(i));B(i)=B(i)+Y(k)*sin(i*X(i));endendA=2*A/N;B=2*B/N;T=2*sum(Y)/N;%Çó½â¸µÀïÒ¶Õ¹¿ªÏµÊýfor i=1:My=y+A(i)*cos(i*x)+B(i)*sin(i*x);end%µÃÈý½Ç¶àÏîʽy=y+T/2;figure(1)plot(X,Y,'r*')figure(2)%»-ͼezplot(y)贝塞尔曲线拟合:clearY=[58 58 58 58 58 57 57 57 58 60 64 67 68 66 6665 64 63 63 62 61 60 60 59 58];X=[0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 1718 19 20 21 22 23 24];%ÒÑÖªµãt=0:0.01:1;x1=X(1)*(1-t).^3+X(2)*3*t.*(1-t).^2+X(3)*3*t.^2.*(1-t)+X(4)*t.^3;y1=Y(1)*(1-t).^3+Y(2)*3*t.*(1-t).^2+Y(3)*3*t.^2.*(1-t)+Y(4)*t.^3;x2=X(4)*(1-t).^3+X(5)*3*t.*(1-t).^2+X(6)*3*t.^2.*(1-t)+X(7)*t.^3;y2=Y(4)*(1-t).^3+Y(5)*3*t.*(1-t).^2+Y(6)*3*t.^2.*(1-t)+Y(7)*t.^3;x3=X(7)*(1-t).^3+X(8)*3*t.*(1-t).^2+X(9)*3*t.^2.*(1-t)+X(10)*t.^3;y3=Y(7)*(1-t).^3+Y(8)*3*t.*(1-t).^2+Y(9)*3*t.^2.*(1-t)+Y(10)*t.^3;x4=X(10)*(1-t).^3+X(11)*3*t.*(1-t).^2+X(12)*3*t.^2.*(1-t)+X(13)*t.^3;y4=Y(10)*(1-t).^3+Y(11)*3*t.*(1-t).^2+Y(12)*3*t.^2.*(1-t)+Y(13)*t.^3;x5=X(13)*(1-t).^3+X(14)*3*t.*(1-t).^2+X(15)*3*t.^2.*(1-t)+X(16)*t.^3;y5=Y(13)*(1-t).^3+Y(14)*3*t.*(1-t).^2+Y(15)*3*t.^2.*(1-t)+Y(16)*t.^3;x6=X(16)*(1-t).^3+X(17)*3*t.*(1-t).^2+X(18)*3*t.^2.*(1-t)+X(19)*t.^3;y6=Y(16)*(1-t).^3+Y(17)*3*t.*(1-t).^2+Y(18)*3*t.^2.*(1-t)+Y(19)*t.^3;x7=X(19)*(1-t).^3+X(20)*3*t.*(1-t).^2+X(21)* 3*t.^2.*(1-t)+X(22)*t.^3;y7=Y(19)*(1-t).^3+Y(20)*3*t.*(1-t).^2+Y(21)* 3*t.^2.*(1-t)+Y(22)*t.^3;x8=X(22)*(1-t).^3+X(23)*3*t.*(1-t).^2+X(24)* 3*t.^2.*(1-t)+X(24)*t.^3;y8=Y(22)*(1-t).^3+Y(23)*3*t.*(1-t).^2+Y(24)* 3*t.^2.*(1-t)+Y(24)*t.^3;%µÃ±´Èû¶ûº¯Êýplot(x1,y1)hold onplot(x2,y2)hold onplot(x3,y3)hold onplot(x4,y4)hold onplot(x5,y5)hold onplot(x6,y6)hold onplot(x7,y7)hold onplot(x8,y8)%»-ͼhold onplot(X,Y,'r*')。