温度分布的曲线拟合.docx
- 文档编号:22992266
- 上传时间:2023-04-29
- 格式:DOCX
- 页数:20
- 大小:294.73KB
温度分布的曲线拟合.docx
《温度分布的曲线拟合.docx》由会员分享,可在线阅读,更多相关《温度分布的曲线拟合.docx(20页珍藏版)》请在冰豆网上搜索。
温度分布的曲线拟合
温度分布的曲线拟合
学号:
XX姓名:
XXX
1.实验描述
美国洛杉矶郊区11月8日的温度(华氏温度)如表1所示。
采用24小时制表1温度数据
时间,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
要求:
1.线性的最小二乘拟合
2.曲线的最小二乘抛物线拟合;
3.三次样条插值拟合
4.T7的三角多项式拟合
5.有4个控制点的贝塞尔曲线拟合
2.实验内容
、线性最小二乘拟合定理5.1(最小二乘拟合曲线)设{(x^yk)}:
#有N个点,其中横坐标{兀}仁是确定的
最小—乘拟合曲线
(1)
y=AxB
的系数是下列线性方程组的解,这些方程称为正规方程:
fN2)fN)N
lZXkIA+庄XkB=Ex NN 'XkANB八yk2k/ 核心代码为: %求方程组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) gridon legendf已知点','最小二乘拟合') 主要算法为: (1).输入x,y; NNNN (2).求正规方程的系数7x2? 'Xk^yk? 'Xkyk k二k二k斗心 (3).解正规方程组am=b(4).绘制拟合曲线 图1线性的最小二乘拟合流程图 、曲线的最小二乘抛物线拟合 定理5.3(最小二乘抛物线拟合)设{(Xk,yQ}Nm有N个点,横坐标是确定的 二乘抛物线的系数表示为 y=f(x)二Ax2BxC 求解A,B和C的线性方程组为 『N\「N、(“)N IZx: IA+! 送X;B+! EX: C=送ykx: 1心丿I心丿1心丿k¥ fN3)fN2)")N .区xkIA+gxkB+: 瓦xkC=£ykxk1心丿1心丿1心丿kT NNN |Hx: A亠―xk|BNC八yk kd? kdkm 最小 ⑶ ⑷ 根据式(4),核心代码为: a(1,1)=sum(x.A4);a(2,3)=sum(x);b (1)=(x.A2)*y'; b (2)=x*y'; %求方程组am=b的根 m=a\b; 算法流程图为: 图2抛物线的最小二乘拟合流程图 三、三次样条插值拟合 定义5.1设{(Xk,yj}打有N-1个点,其中a=X。 : : : 为: : : |1(: : : Xn: : : b。 如果存在N个三次多项式Sk(X),系数为Sk,o,Sk,i,Sk,2和Sk,3,满足如下性质: IS(x)=Sk,0+Sk,i(x-xk)+Sk,2(X-Xk)+Sk,3(X-Xk) II怒)二Yk III.Sk(xk1^-Sk1(xk1)iv.Sk(xk.J=q1(xk1) V.Sk'(xk1)=Sk1(xk1) x•[xk,xk1],k=0,11(,N-1k=O,1,|l|,N (5) k=O,1,|l|,N-2 k=0,1川,N—2 k=0,1,|l|,N-2 则称函数为三次样条函数。 令mk=s"(Xk),mk1=s"(Xk1),hk二x—-Xk和dk=山吐,可得包含口㈠皿和hk mki的重要关系式: hkjmkj2(hkj•hk)mkhEk1二山(6) 其中Uk=6(dk-dk」),k=1,2」|(,N-1 方程组(6)中的未知数是要求的值{叫},而且其他的项是可以通过数据点集{(Xk.Yk)}进行简单数学计算得到的常量。 因此方程组(6)是包含N1个未知数,具有N-1个线性方程组的不定方程组。 所以需要另外两个方程组才能求解。 可通过它们消去方程组(6)中的第一个方程的m0和第个方程的mN。 如果给定m。 ,则可以计算出m°ho,而且方程组⑹的第一个方程(当k=1时)为: 2(hoh)m1gm2二5-hom。 ⑺ 如果给定mN,则可以计算出hN^mN,而且方程组⑹的最后一个方程(当k=N-1时)为: hN_2mN-2'2(hNJ2hN4)mN4~UN4_hN-jmN(8) 考虑方程组⑹以及方程组⑺和方程组(8),其中k二2,3川I,N-2,可形成N-1阶线性方程组,包含系数m1,m2,,mN4。 重写方程组(6)中的方程1到方程N-1方程组HM=V,表示为: b1 a1 当得到系数{叫}后,可以用如下公式计算Sk(x)的样条系数{Sk,j} sk,o=yk, Sk’^dk-h^g叫1) sk,2- mik (10) 6hk 为了更有效地计算,每个三次多项式 Sk(x)可表示成嵌套形式: (11) Sk(x)=((Sk,3WSk,2)wSk,i)wyk,其中w^x- 其中Sjx)在给区间Xk^x乞Xki内使用。 核心代码为: fork=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); end %求m(0)和m(N) M (1)=2*(D (1)-dxO)/H (1)-M (2)/2; M(N+1)=3*(dxn-D(N))/H(N)-M(N)/2; %求样条系数s(k,j) fork=0: N-1 S(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 算法流程图为: 图3三次样条拟合流程图 四、T7的三角多项式拟合 定义5.4具有如下形式的级数; 称为M阶的三角多项式。 定理5.8(离散傅里叶级数)设有N1个点{(&』"}: £,其中yj二f(x),而且横 坐标之间等距,即: Xj=7M'j,其中j=0,1」川,N(13) Nn 如果f(x)的周期为2二,而且2M: : : N,则存在式(12)所示的三角多项式Tm(x), 使得下式的值最小 N 2 、(f(Xk)-TM(x)) k4 多项式的系数 j和bj可通过如下公式计算: 2N ajf(Xk)cos(jXk),其中j=0,1,|)|,M Nk4 2N bjf(Xk)sin(jXk),其中j=0,1,|l(,M (14) (15) (16) 图4三角多项式拟合流程图 核心代码为: %计算A和B forj=1: M Nk4 A(j+1)=cos(j*X)*Y'; B(j+1)=sin(j*X)*Y'; end %求三角多项式T T=A (1); forj=1: M T=T+A(j+1)*cos(j*x)+B(j+1)*sin(j*x); end 五、有4个控制点的贝塞尔曲线拟合 定义5.5N阶伯恩斯坦多项式定义为 N.B,N(t)二匚t,(1-t)亠(17) 廿+心N! i=0,1,2,N,其中VUi! (N-i)! 定义5.6给定一个控制点集,{R}寫,其中R,定义 N P(t)八RBi,N(t)(18) i=0 为N阶贝塞尔曲线,其中Bi,N(t),i=0,1,,N,是N阶伯恩斯坦多项式,"[0,1]。 公式(18)中的控制点是表示平面中x和y坐标的有序对。 可将控制点作为向量,儿对 应的伯恩斯坦多项式作为标量处理, 这样公式(18)可参数化表示为 P(t)=(x(t),yt其中 N X(t)二"xiBi,N(t) i=0 核心代码为: %B%3阶伯恩斯坦多项式系数矩阵 B=[(1-t)A3,3*t*(1-t)A2,3*tA2*(1-t),tA3]; %xx为3阶贝塞尔曲线横坐标,yy为纵坐标xx=0;yy=0;fori=n+1: n+4 xx=xx+x(i)*B(i-n); yy=yy+y(i)*B(i-n); end 图5贝塞尔曲线拟合流程图 3.实验结果及分析 从图6中看出由于温度变化快,数据比较分散,不宜用最小二乘直线进行拟合 图7最小二乘抛物线拟合曲线 比较图6和图7,显然图7中的拟合曲线更贴近于温度点坐标。 而图7能够直观的反映一天中的温度变化趋势。 但所拟合的曲线明显与数据点偏离较大。 从图8中可以很容易看出拟合曲线将经过样点,读出该地区一天的温度变化走势: 上午时段温度快速升高,到12点左右达到最高。 下午时段,温度慢慢降低,直到24时温度降至最低并保持一段时间。 从图9中可以很容易读出该地区一天的温度变化走势,与图8所反映的温度变化 规律大致相同。 与图三次样条曲线相比,三角多项式拟合出的曲线更光滑。 图10贝塞尔曲线拟合曲线 图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; x1=1: 0.1: 24; y1=m (1)*x1+m (2); %绘图,其中(x,y)为已知点,用红色的星号表示,y1为拟合曲线 plot(x,y,'*r',x1,y1) gridon legendf已知点','最小二乘拟合') 二、曲线的最小二乘抛物线拟合; clear %时间 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,b a=zeros(3,3); b=zeros(3,1); a(1,1)=sum(x.A4); a(1,2)=sum(x.A3); a(1,3)=sum(x.A2); a(2,1)=a(1,2); a(2,2)=a(1,3); a(2,3)=sum(x); a(3,1)=a(2,2); a(3,2)=a(2,3); a(3,3)=24; b (1)=(x.A2)*y'; b (2)=x*y'; b(3)=sum(y); %求方程组am=b的根 m=a\b; x1=1: 0.1: 24; y1=m (1)*x142+m (2).*x1+m(3); %绘图,其中(x,y)为已知点,用红色的星号表示,y1最小二乘抛物线 plot(x,y,'*r',x1,y1) gridon legendf已知点','最小二乘抛物线') 三、三次样条插值拟合 clear N=23; %时间 X=1: N+1; %温度 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]; %dxO=S'(xO),dxn=S'(xn),为自由边界条件 dx0=0;dxn=0; H=diff(X); D=diff(Y)./H; A=H(2: N-1); B=2*(H(1: N-1)+H(2: N)); C=H(2: N); U=6*diff(D); %clampedsplineendpointconstraints 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*(dx0-D(N)); fork=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); end M(N)=U(N-1)/B(N-1); fork=N-2: -1: 1 M(k+1)=(U(k)-C(k)*M(k+2))/B(k); end %求m(0)和m(N) M⑴=2*(D⑴-dxO)/H (1)-M (2)/2; M(N+1)=3*(dxn-D(N))/H(N)-M(N)/2; %求样条系数s(k,j) fork=0: N-1 S(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 %绘图,其中(X,Y)为已知点,用红色的星号表示,y为三次样条曲线 forj=1: N x=j: 0.01: j+1; y=polyval(S(j,: ),x-X(j)); plot(X,Y,'*r',x,y),holdon end gridon legend('已知点','三次样条曲线') 四、T7的三角多项式拟合 symsx N=23; M=7; %时间 X=-pi: 2*pi/N: pi; %温度 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为包含cos(jx)系数的矩阵 %B为包含sin(jx)系数的矩阵 A=zeros(1,M+1); B=zeros(1,M+1); %Yends为非连续点处的值 Yends=(丫⑴+Y(N))/2; Y (1)=Yends; Y(N)=Yends; %计算A和B A (1)=sum(Y); forj=1: M A(j+1)=cos(j*X)*Y'; B(j+1)=sin(j*X)*Y'; end A=2*A/N; B=2*B/N; A (1)=A (1)/2; %求三角多项式T T=A (1); forj=1: M T=T+A(j+1)*cos(j*x)+B(j+1)*sin(j*x); end %绘图,其中(X,Y)为已知点,用红色的星号表示,y为 x1=-pi: .01: pi; y1=subs(T,x,x1); plot(X,Y,'*r',x1,y1) gridon legendf已知点','三角多项式曲线') 五、有4个控制点的贝塞尔曲线拟合 clear %时间,第25个值为构造值 x=1: 24+1; %温度,第25个值为构造值 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,58]; symst %B%3阶伯恩斯坦多项式系数矩阵 B=[(1-t)A3,3*t*(1-t)A2,3*tA2*(1-t),tA3];n=0; fork=1: 8%xx为3阶贝塞尔曲线横坐标,yy为纵坐标 xx=0;yy=0; fori=n+1: n+4 xx=xx+x(i)*B(i-n); yy=yy+y(i)*B(i-n); end n=n+3; tt=0: 0.1: 1; bx=subs(xx,t,tt); by=subs(yy,t,tt); %绘图,其中(x,y)为已知点,用红色的星号表示,by为贝塞尔曲线 plot(x,y,'*r',bx,by),holdon legendf已知点','贝塞尔曲线') end
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 温度 分布 曲线拟合