温度分布的曲线拟合

《数值计算方法》实验报告1温度分布的曲线拟合学号:XX姓名：XXX1.实验描述美国洛杉矶郊区11月8日的温度（华氏温度）如表1所示。采用24小时制。表1温度数据时间，p.m.温度时间，a.m.温度1661582662583653584644585635576636577627578618589609601060106411591167午夜58正午68要求：1.线性的最小二乘拟合2.曲线的最小二乘抛物线拟合；3.三次样条插值拟合4.T7的三角多项式拟合5.有4个控制点的贝塞尔曲线拟合2.实验内容一、线性最小二乘拟合定理5.1（最小二乘拟合曲线）设1{(,)}Nkkkxy有N个点，其中横坐标1{}Nkkx是确定的。《数值计算方法》实验报告2最小二乘拟合曲线yAxB(1)的系数是下列线性方程组的解，这些方程称为正规方程：211111NNNkkkkkkkNNkkkkxAxBxyxANBy(2)核心代码为：%求方程组am=b的根m=a\b;x1=1:0.1:24;y1=m(1)*x1+m(2);%绘图，其中(x,y)为已知点，用红色的星号表示，y1为拟合曲线plot(x,y,'*r',x1,y1)gridonlegend('已知点','最小二乘拟合')主要算法为：(1).输入x,y；(2).求正规方程的系数21Nkkx,1Nkkx,1Nkky,1Nkkkxy(3).解正规方程组am=b(4).绘制拟合曲线《数值计算方法》实验报告3二、曲线的最小二乘抛物线拟合定理5.3（最小二乘抛物线拟合）设1{(,)}Nkkkxy有N个点，横坐标是确定的。最小二乘抛物线的系数表示为2()yfxAxBxC(3)求解,AB和C的线性方程组为432211113211112111NNNNkkkkkkkkkNNNNkkkkkkkkkNNNkkkkkkxAxBxCyxxAxBxCyxxAxBNCy(4)根据式(4)，核心代码为：a(1,1)=sum(x.^4);a(2,3)=sum(x);b(1)=(x.^2)*y';开始输入x，y21Nkkx,1Nkkx,1Nkky,1Nkkkxy解正规方程组am=b绘图结束图1线性的最小二乘拟合流程图《数值计算方法》实验报告4b(2)=x*y';%求方程组am=b的根m=a\b;算法流程图为：三、三次样条插值拟合定义5.1设1{(,)}Nkkkxy有1N个点，其中01Naxxxb。如果存在N个三次多项式()kSx，系数为,0,1,2,,kkksss和,3ks，满足如下性质：开始输入x，y21Nkkx,1Nkkx,1Nkky,1Nkkkxy解正规方程组am=b绘图结束2()fxAxBxC图2抛物线的最小二乘拟合流程图《数值计算方法》实验报告523kk,0k,1kk,2kk,3k1111''111''''111I.S(x)=s+s(x-x)+s(x-x)+s(x-x)[,],0,1,,1.()0,1,,.()()0,1,,2.()()0,1,,2.()()0,1,,2kkkkkkkkkkkkkkkkxxxkNIISxykNIIISxSxkNIVSxSxkNVSxSxkN(5)则称函数为三次样条函数。令''''11(),()kkkkmSxmSx，1kkkhxx和1kkkkyydh，可得包含1,kkmm和1km的重要关系式：11112()kkkkkkkkhmhhmhmu(6)其中16(),1,2,,1kkkuddkN方程组(6)中的未知数是要求的值{}km，而且其他的项是可以通过数据点集{(,)}kkxy进行简单数学计算得到的常量。因此方程组(6)是包含1N个未知数，具有1N个线性方程组的不定方程组。所以需要另外两个方程组才能求解。可通过它们消去方程组(6)中的第一个方程的0m和第个方程的Nm。如果给定0m，则可以计算出00mh，而且方程组(6)的第一个方程（当k=1时）为：011121002()hhmhmuhm(7)如果给定Nm，则可以计算出1NNhm，而且方程组(6)的最后一个方程（当k=N-1时）为：22211112()NNNNNNNNhmhhmuhm(8)考虑方程组(6)以及方程组(7)和方程组(8)，其中2,3,,2kN，可形成1N阶线性方程组，包含系数121,,,Nmmm。重写方程组(6)中的方程1到方程1N，得到一个包含121,,,Nmmm的三角线性方程组HMV，表示为：111222232222211111NNNNNNNNNbcmvabcmvabcmvabmv(9)当得到系数{}km后，可以用如下公式计算()kSx的样条系数,{}kjS。《数值计算方法》实验报告61,0,11,2,3(2),6,26kkkkkkkkkkkkkhkmmsysdmmmssdh(10)为了更有效地计算，每个三次多项式()kSx可表示成嵌套形式：,3,2,1()(()),其中kkkkkkSxswswswywxx(11)其中()kSx在给区间1kkxxx内使用。核心代码为：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);end%求m(0)和m(N）M(1)=2*(D(1)-dx0)/H(1)-M(2)/2;M(N+1)=3*(dxn-D(N))/H(N)-M(N)/2;%求样条系数s(k,j)fork=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算法流程图为：《数值计算方法》实验报告7四、T7的三角多项式拟合定义5.4具有如下形式的级数；0()(cos()sin())2MjjaTxajxbjx(12)称为M阶的三角多项式。定理5.8（离散傅里叶级数）设有1N个点0{(,)}Nkkjxy，其中()jyfx，而且横坐标之间等距，即：21,0,1,,其中jjxjNNn(13)如果()fx的周期为2，而且2MN，则存在式(12)所示的三角多项式()MTx，使得下式的值最小。开始结束绘图1kkkkyydh16()kkkudd构造HM=V样条系数,{}kjS1kkkhxx输入X，Y图3三次样条拟合流程图《数值计算方法》实验报告821(()())NkMkfxTx(14)多项式的系数ja和jb可通过如下公式计算：12()cos(),0,1,,NjkkkafxjxjMN其中(15)12()sin(),0,1,,NjkkkbfxjxjMN其中(16)核心代码为：%计算A和Bforj=1:M开始结束输入X和YjMA(j+1)=cos(j*X)*Y'B(j+1)=sin(j*X)*Y';B(j+1)=sin(j*X)*Y';三角多项式T绘图j=1j=j+1图4三角多项式拟合流程图YN《数值计算方法》实验报告9A(j+1)=cos(j*X)*Y';B(j+1)=sin(j*X)*Y';end%求三角多项式TT=A(1);forj=1:MT=T+A(j+1)*cos(j*x)+B(j+1)*sin(j*x);end五、有4个控制点的贝塞尔曲线拟合定义5.5N阶伯恩斯坦多项式定义为,()(1)NiNiiNiBttt(17)!0,1,2,,!()!NiNiNiNi其中定义5.6给定一个控制点集，0{}NiiP，其中(,)iiiPxy，定义,()0()NiiNtiPtPB(18)为N阶贝塞尔曲线，其中,(),0,1,,iNtBiN，是N阶伯恩斯坦多项式，[0,1]t。公式(18)中的控制点是表示平面中x和y坐标的有序对。可将控制点作为向量，儿对应的伯恩斯坦多项式作为标量处理，这样公式(18)可参数化表示为()((),())Ptxtyt，其中,(),()00()()NNiiNtiiNtiixtxBytyB(19)核心代码为：%B为3阶伯恩斯坦多项式系数矩阵B=[(1-t)^3,3*t*(1-t)^2,3*t^2*(1-t),t^3];%xx为3阶贝塞尔曲线横坐标，yy为纵坐标xx=0;yy=0;fori=n+1:n+4xx=xx+x(i)*B(i-n);《数值计算方法》实验报告10yy=yy+y(i)*B(i-n);end开始结束n+3=nk+1=k0=n,1=k绘图xx+x(i)*B(i-n)=xxyy+y(i)*B(i-n)=yyk8输入X和Yi+1=in+1=iin+4YNNY图5贝塞尔曲线拟合流程图《数值计算方法》实验报告113.实验结果及分析从图6中看出由于温度变化快，数据比较分散，不宜用最小二乘直线进行拟合。图7最小二乘抛物线拟合曲线图6线性最小二乘拟合曲线《数值计算方法》实验报告12比较图6和图7，显然图7中的拟合曲线更贴近于温度点坐标。而图7能够直观的反映一天中的温度变化趋势。但所拟合的曲线明显与数据点偏离较大。从图8中可以很容易看出拟合曲线将经过样点，读出该地区一天的温度变化走势：上午时段温度快速升高，到12点左右达到最高。下午时段，温度慢慢降低，直到24时温度降至最低并保持一段时间。图8三次样条拟合曲线《数值计算方法》实验报告13从图9中可以很容易读出该地区一天的温度变化走势，与图8所反映的温度变化规律大致相同。与图三次样条曲线相比，三角多项式拟合出的曲线更光滑。图9三角多项式拟合曲线图10贝塞尔曲线拟合曲线《数值计算方法》实验报告14图8与图10所反映的温度变化规律大致相同,能够很容易地看出该地区一天温度变化的走势。图10采用4个控制点的贝塞尔拟合曲线在每段曲线交点出，其斜率变化较大4.结论1.根据五种不同的拟合结果可知，由于线性最小二乘拟合和抛物线最小二乘拟合误差很大，三角多项式拟合次之，而三次样条拟合和贝塞尔曲线拟合最为精确，误差最小。2.从本次实验可知，在选择最佳的拟合方式前可以先绘制出已知点的离散分布图像，根据其变化趋势合理选择最佳拟合方式。3.在本实验中，三次样条拟合和贝塞尔曲线拟合最为精确，但两者相比，三次样条拟合程序要复杂和繁琐些，而贝塞尔曲线拟合较为简单，故本实验中的最佳拟合方式为贝塞尔曲线拟合。附件（代码）一、线性的最小二乘拟合%时间x=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];%初始化正规方程的系数矩阵a,ba=zeros(2,2);b=zeros(2,1);a(1,1)=x*x';a(1,2)=sum(x);a(2,1)=a(1,2);a(2,2)=24;b(1)=x*y';b(2)=sum(y);%求方程组am=b的根m=a\b;《数值计算方法》实验报告15x1=1:0.1:24;y1=m(1)*x1+m(2);%绘图，其中(x,y)为已知点，用红色的星号表示，y1为拟合曲线plot(x,y,'*r',x1,y1)gridonlegend('已知点','最小二乘拟合')二、曲线的最小二乘抛物线拟合；clear%时间x=1:24;%温度y=[58,5