计算方法(b)上机实习报告

1、硕士课程报告 课程名称 计算方法B上机报告 任课老师 马军 老师 学 院 电气工程学院 班 级 硕4023 学生姓名 张鹏飞 学 号 3114164006 日 期 2021年12月24日 目 录 TOC o 1-3 h z u HYPERLINK l _Toc407162412 1.实验一 PAGEREF _Toc407162412 h 1 HYPERLINK l _Toc407162413 1.1实验题目 PAGEREF _Toc407162413 h 1 HYPERLINK l _Toc407162414 1.2算法组织 PAGEREF _Toc407162414 h 1 HYPERLIN

2、K l _Toc407162415 1.3算法的Matlab实现 PAGEREF _Toc407162415 h 1 HYPERLINK l _Toc407162416 1.4实验结果及分析 PAGEREF _Toc407162416 h 3 HYPERLINK l _Toc407162417 2.实验二 PAGEREF _Toc407162417 h 4 HYPERLINK l _Toc407162418 2.1实验题目 PAGEREF _Toc407162418 h 4 HYPERLINK l _Toc407162419 2.2三次样条插值拟合数据点 PAGEREF _Toc4071624

3、19 h 4 HYPERLINK l _Toc407162420 2.2.1算法组织 PAGEREF _Toc407162420 h 4 HYPERLINK l _Toc407162421 2.2.2算法的Matlab实现 PAGEREF _Toc407162421 h 6 HYPERLINK l _Toc407162422 2.2.3实验结果及分析 PAGEREF _Toc407162422 h 11 HYPERLINK l _Toc407162423 2.3复化Simpson公式近似计算光缆长度 PAGEREF _Toc407162423 h 14 HYPERLINK l _Toc4071

4、62424 2.3.1算法组织 PAGEREF _Toc407162424 h 14 HYPERLINK l _Toc407162425 2.3.2算法的Matlab实现 PAGEREF _Toc407162425 h 16 HYPERLINK l _Toc407162426 2.3.3实验结果及分析 PAGEREF _Toc407162426 h 16 HYPERLINK l _Toc407162427 3.实验三 PAGEREF _Toc407162427 h 17 HYPERLINK l _Toc407162428 3.1实验题目 PAGEREF _Toc407162428 h 17 H

5、YPERLINK l _Toc407162429 3.2算法组织 PAGEREF _Toc407162429 h 17 HYPERLINK l _Toc407162430 3.3算法的Matlab实现 PAGEREF _Toc407162430 h 18 HYPERLINK l _Toc407162431 3.3.1法方程最小二乘法 PAGEREF _Toc407162431 h 18 HYPERLINK l _Toc407162432 3.3.2QR分解 PAGEREF _Toc407162432 h 20 HYPERLINK l _Toc407162433 3.4实验结果及分析 PAGER

6、EF _Toc407162433 h 25 HYPERLINK l _Toc407162434 3.4.1法方程最小二乘法输出结果 PAGEREF _Toc407162434 h 25 HYPERLINK l _Toc407162435 3.4.2QR分解输出结果 PAGEREF _Toc407162435 h 28 HYPERLINK l _Toc407162436 3.4.3结果分析 PAGEREF _Toc407162436 h 29 HYPERLINK l _Toc407162437 4.实验四 PAGEREF _Toc407162437 h 30 HYPERLINK l _Toc40

7、7162438 4.1实验题目 PAGEREF _Toc407162438 h 30 HYPERLINK l _Toc407162439 4.2二分法算法组织 PAGEREF _Toc407162439 h 30 HYPERLINK l _Toc407162440 4.2.1算法思想及依据 PAGEREF _Toc407162440 h 30 HYPERLINK l _Toc407162441 4.2.2算法结构 PAGEREF _Toc407162441 h 30 HYPERLINK l _Toc407162442 4.3二分法算法的Matlab实现 PAGEREF _Toc40716244

8、2 h 31 HYPERLINK l _Toc407162443 4.3.1利用Matlab校验根的区间 PAGEREF _Toc407162443 h 31 HYPERLINK l _Toc407162444 4.3.2使用二分法Matlab程序 PAGEREF _Toc407162444 h 32 HYPERLINK l _Toc407162445 4.4二分法实验结果及分析 PAGEREF _Toc407162445 h 33 HYPERLINK l _Toc407162446 4.5割线法算法组织 PAGEREF _Toc407162446 h 34 HYPERLINK l _Toc4

9、07162447 4.5.1算法思想及依据 PAGEREF _Toc407162447 h 34 HYPERLINK l _Toc407162448 4.5.2算法结构 PAGEREF _Toc407162448 h 34 HYPERLINK l _Toc407162449 4.6割线法算法的Matlab实现 PAGEREF _Toc407162449 h 34 HYPERLINK l _Toc407162450 4.7割线法实验结果及分析 PAGEREF _Toc407162450 h 35 HYPERLINK l _Toc407162451 5.实验五 PAGEREF _Toc407162

10、451 h 36 HYPERLINK l _Toc407162452 5.1实验题目 PAGEREF _Toc407162452 h 36 HYPERLINK l _Toc407162453 5.2算法组织 PAGEREF _Toc407162453 h 36 HYPERLINK l _Toc407162454 5.2.1算法思想及依据 PAGEREF _Toc407162454 h 36 HYPERLINK l _Toc407162455 5.2.2算法结构 PAGEREF _Toc407162455 h 36 HYPERLINK l _Toc407162456 5.3算法的Matlab实现

11、 PAGEREF _Toc407162456 h 37 HYPERLINK l _Toc407162457 5.3.1非压缩带状对角方程组 PAGEREF _Toc407162457 h 37 HYPERLINK l _Toc407162458 5.3.2压缩带状对角方程组 PAGEREF _Toc407162458 h 39 HYPERLINK l _Toc407162459 5.4实验结果及分析 PAGEREF _Toc407162459 h 41 HYPERLINK l _Toc407162460 5.4.1Matlab运行结果 PAGEREF _Toc407162460 h 41 HY

12、PERLINK l _Toc407162461 5.4.2Matlab运行结果分析 PAGEREF _Toc407162461 h 42 HYPERLINK l _Toc407162462 5.5实际问题 PAGEREF _Toc407162462 h 431.实验一1.1实验题目对以下和式计算：，要求：(1)假设只需保存11个有效数字，该如何进行计算；(2)假设要保存30个有效数字，那么又将如何进行计算；1.2算法组织1、在首先满足题目要求的情况下，利用后验误差估计计算出算法中的最大运行次数。即以一下条件作为迭代终止的准那么：其中为误差数值。该实验第一问利用后验误差估计截取的项数，使误差小于

13、1e-10；第二问利用后验误差估计截取的项数，使误差小于1e-30。得出运行次数，作为累加的终止条件。2、为减小舍入误差在计算S时所采用的方法是逆序相加，其依据是：S中各项的绝对值为递减的正数，而两个数量级相差较大的数字相加减时，较小数的有效数字会被丧失，从而导致最后的运算结果失真。为防止“大数吃小数现象的发生，采用逆序相加。3、对于实现30位有效数字，那么调用从工具箱中vpa变量，精度位数即可实现。1.3算法的Matlab实现clc;clear;sum1=0;sum2=0;n=0;N=0;%第一问利用后验误差估计截取的项数，使误差小于1e-10while abs(4/(8*(n+1)+1)-

14、2/(8*(n+1)+4)-1/(8*(n+1)+5)-1/(8*(n+1)+6)/(16.(n+1)-(4/(8*(n)+1)-2/(8*(n)+4)-1/(8*(n)+5)-1/(8*(n)+6)/(16.n)1e-10 n=n+1endn;%限定误差disp(最大误差:)e1=abs(4/(8*(n+1)+1)-2/(8*(n+1)+4)-1/(8*(n+1)+5)-1/(8*(n+1)+6)/(16.(n+1)-(4/(8*(n)+1)-2/(8*(n)+4)-1/(8*(n)+5)-1/(8*(n)+6)/(16.n)%截取的项数disp(截取的项数:)n%按从小到大的顺序相加for

15、 i=n:-1:0 sum1=sum1+(4/(8*(i)+1)-2/(8*(i)+4)-1/(8*(i)+5)-1/(8*(i)+6)/(16i);enddisp(保存11位有效数字估算值:)fprintf(%5.10fnn,sum1)%第二问利用后验误差估计截取的项数，使误差小于1e-30while abs(4/(8*(N+1)+1)-2/(8*(N+1)+4)-1/(8*(N+1)+5)-1/(8*(N+1)+6)/(16.(N+1)-(4/(8*(N)+1)-2/(8*(N)+4)-1/(8*(N)+5)-1/(8*(N)+6)/(16.N)1e-30 N=N+1endN;%限定误差d

16、isp(最大误差:)e2=abs(4/(8*(N+1)+1)-2/(8*(N+1)+4)-1/(8*(N+1)+5)-1/(8*(N+1)+6)/(16.(N+1)-(4/(8*(N)+1)-2/(8*(N)+4)-1/(8*(N)+5)-1/(8*(N)+6)/(16.N)%截取的项数disp(截取的项数:)N%按从小到大的顺序相加for j=N:-1:0 sum2=sum2+(4/(8*(j)+1)-2/(8*(j)+4)-1/(8*(j)+5)-1/(8*(j)+6)/(16j);enddisp(保存30位有效数字估算值:)vpa(sum2,30)1.4实验结果及分析通过运行上述程序可知

17、：第一问中，在满足误差小于1e-10时，截取的项数为7，最大误差:e1 =1.5312e-11，保存11位有效数字估算值:3.1415926536。第二问中，在满足误差小于1e-30时，截取的项数为23，最大误差: e2 =8.1687e-32，保存30位有效数字估算值:3.14159265358979323846264338328。从上述的算法思想中可以看出，浮点数运算中不仅要满足误差要求，还要尽可能地减少计算量，而其计算量的大小很大程度上是由算式要求的精度来决定的；此外还要考虑舍入误差的影响，这时就要对所运算数据的性质进行分析，设置适宜的算法，从而提高运算的精度。2.实验二2.1实验题目某

18、通信公司在一次施工中，需要在水面宽度为20米的河沟底部沿直线走向铺设一条沟底光缆。在铺设光缆之前需要对沟底的地形进行初步探测，从而估计所需光缆的长度，为工程预算提供依据。已探测到一组等分点位置的深度数据(单位:米)如下表所示：分点0123456深度9.018.967.967.978.029.0510.13分点78910111213深度11.1812.2613.2813.3212.6111.2910.22分点14151617181920深度9.157.907.958.869.8110.8010.93(1)请用适宜的曲线拟合所测数据点；(2)预测所需光缆长度的近似值，并作出铺设河底光缆的曲线图；2

19、.2三次样条插值拟合数据点2.2.1算法组织1、算法思想及依据在分段插值多项式中，常用的是三次多项式。即在每一个子区间上是一个三次多项式，这些多项式的曲线在节点处“光滑地连接起来。这样的分段三次多项式称为三次样条函数。样条函数为：记那么上式可写成：这是一个有n-1个方程，n+1个未知数的线性方程组，还需要边界条件才能确定唯一解。1自然样条2指定端点的一阶导数。此时得3假设两端点的导数值均不知道时，可以通过推导得到最后可以得到这些方程组的统一形式此题采用边界条件三求解，其中2、算法结构其算法：SPLINEM()本算法计算三次样条的参数Mi，存放于数组M中，参数根据不同边界条件给定。 1.1 2.

20、1 2.1.1 3. 4. 4.1 4.2 4.3 用算法TSS(a,b,c,M,n,M)算法TTSa，b，c，d，n，x本算法解对角线性方程组，系数存于数组a、b和c中，d是有段向量。2.2.2算法的Matlab实现三次样条插值及复化Simpson公式主程序clcclear%构造数据矩阵By=9.01;8.96;7.96;7.97;8.02;9.05;10.13;11.18;12.26;13.28;13.32;12.61;11.29;10.22;9.15;7.90;7.95;8.86;9.81;10.80;10.93;N=length(y);x=0:20/(N-1):20;Data=zero

21、s(2,N);for i=1:N Data(1,i)=x(i); Data(2,i)=y(i);endData;% 构造步长矩阵hh = zeros(1,N-1);for k =1:N-1h(k)=x(k+1)-x(k);end%利用差分表计算d1，dN及dix=Data(1,:);y=Data(2,:);for i=1:N D(i,1)=x(i); D(i,2)=y(i);endfor i=2:N D(i,3)=(y(i)-y(i-1)/(x(i)-x(i-1);endfor i=3:N D(i,4)=(D(i,3)-D(i-1,3)/(x(i)-x(i-2);endfor i=4:N D(

22、i,5)=(D(i,4)-D(i-1,4)/(x(i)-x(i-3);endD;% 构造M的系数矩阵AA = zeros(N);% 由第三类边界条件可得A(1,2) = -2; % 由第三类边界条件可得A(N,N-1)=-2; for i=1:NA(i,i)=2;endfor i=2:N-1A(i,i+1)=h(i)/(h(i-1)+h(i);A(i,i-1)=h(i-1)/(h(i-1)+h(i);end% 构造M的结果矩阵dd =zeros(N,1);d(1)=-12*h(1)*D(4,5);d(N) =12*h(N-1)*D(N,5);for j=2:N-1d(j) =6*D(j,4);

23、end%利用TSS方法解方程组M=TSS(A,d);%得到各段三次样条插值的函数表达式，并作出曲线图y=Data(2,:);P,N=size(Data);for i=1:N-1 xx=x; Sy=getS(i,xx,M,h,y); sx=xx(i):0.01:xx(i+1); sy=Sy(sx); plot(x,y,o) %原探测数据点及数值 plot(sx,sy,r,LineWidth,1.5) %按插值函数表达式作图 hold on sy=-sy; plot(sx,sy,k,LineWidth,1.5) %作出铺设河底光缆的曲线图end%添加水平线l=line(0,20,0,0);%设置直

24、线的宽度和颜色set(l,Linewidth,1.5)set(l,color,k)%设置坐标轴刻度set(gca,Xtick,0:1:20)%设置坐标轴刻度set(gca,Ytick,-15:1:15)%添加图形标题和坐标轴名称title(一组等分点位置的深度探测数据)gridxlabel(一组等分点的位置)ylabel(深度数据)%利用复化Simpson公式计算光缆的长度y=Data(2,:);P,N=size(Data);sl=0;for i=1:N-1 Sy=getS(i,xx,M,h,y); syms t dsy=diff(Sy(t),t); g=sqrt(1+dsy2); u=i-1

25、; v=i; sl=sl+simpson(g,u,v,20) endSy;g;slTSS解方程组的子程序%追赶法解方程组function x=TSS(A,d)M,N=size(A);b(N)=A(N,N);for i=1:N-1 a(i+1)=A(i+1,i); b(i)=A(i,i); c(i)=A(i,i+1); endu(1)=b(1);y(1)=d(1);for i=2:N l(i)=a(i)/u(i-1); u(i)=b(i)-l(i)*c(i-1); y(i)=d(i)-l(i)*y(i-1);endx(N)=y(N)/u(N);for i=N-1:-1:1 x(i)=(y(i)-

26、c(i)*x(i+1)/u(i);endx;end得到子区间函数表达式的子程序function Sy=getS(k,xx,m,h,y)syms S x Syi=k+1;S=(xx(i)-x)3*m(i-1)/(6*h(i-1)+(x-xx(i-1)3*m(i)/(6*h(i-1)+(y(i-1)-h(i-1)2*m(i-1)/6)*(xx(i)-x)/h(i-1)+(y(i)-h(i-1)2*m(i)/6)*(x-xx(i-1)/h(i-1)Sy=inline(xx(i)-x)3*m(i-1)/(6*h(i-1)+(x-xx(i-1)3*m(i)/(6*h(i-1)+(y(i-1)-h(i-1

27、)2*m(i-1)/6)*(xx(i)-x)/h(i-1)+(y(i)-h(i-1)2*m(i)/6)*(x-xx(i-1)/h(i-1);end2.2.3实验结果及分析由上面程序可以得到如下结果：在时的函数表达式为S =(845002892516739*(x - 1)3)/4503599627370496 - (3313195038365*x)/8796093022208 + (7514075829091499*x3)/54043195528445952 + 5177804441890613/562949953421312在时的函数表达式为S =2621488579161321/281474

28、976710656 - (2489533575770641*(x - 1)3)/6755399441055744 - (7514075829091499*(x - 2)3)/54043195528445952 - (277217198887389*x)/562949953421312在时的函数表达式为S =(650311338860789*(x - 2)3)/1688849860263936 - (418602078066937*x)/562949953421312 + (2489533575770641*(x - 3)3)/6755399441055744 + 276287345834086

29、9/281474976710656在时的函数表达式为S =(2624680819031*x)/4398046511104 - (1092494410535677*(x - 3)3)/6755399441055744 - (650311338860789*(x - 4)3)/1688849860263936 + 3262063247973023/562949953421312在时的函数表达式为S =(79721223946041*x)/140737488355328 + (8155793057367127*(x - 4)3)/27021597764222976 + (10924944105356

30、77*(x - 5)3)/6755399441055744 + 3330360244180239/562949953421312在时的函数表达式为S =(814815571273233*x)/562949953421312 - (3544057556774495*(x - 5)3)/54043195528445952 - (8155793057367127*(x - 6)3)/27021597764222976 + 425353433367447/281474976710656在时的函数表达式为S =(274137988850561*x)/281474976710656 + (56680388

31、8785965*(x - 6)3)/54043195528445952 + (3544057556774495*(x - 7)3)/54043195528445952 + 306243053520945/70368744177664在时的函数表达式为S =(617478217955525*x)/562949953421312 - (1377815456891219*(x - 7)3)/216172782113783808 - (566803888785965*(x - 8)3)/54043195528445952 + 1965528746386739/562949953421312在时的函数表

32、达式为(545284355633681*x)/562949953421312 + (405384572326449*(x - 8)3)/9007199254740992 + (1377815456891219*(x - 9)3)/216172782113783808 + 2543079644961491/562949953421312在时的函数表达式为S =(179389613108493*x)/562949953421312 - (3156841900829689*(x - 9)3)/13510798882111488 - (405384572326449*(x - 10)3)/900719

33、9254740992 + 5836132327688183/562949953421312在时的函数表达式为S =(3156841900829689*(x - 11)3)/13510798882111488 - (814194773093437*(x - 10)3)/9007199254740992 - (480342372812029*x)/562949953421312 + 12433452186893403/562949953421312在时的函数表达式为S =(814194773093437*(x - 12)3)/9007199254740992 - (1394059081462219

34、*(x - 11)3)/9007199254740992 - (706852419243083*x)/562949953421312 + 14925062697634997/562949953421312在时的函数表达式为S =(5376237321301877*(x - 12)3)/54043195528445952 - (745487614849087*x)/562949953421312 + (1394059081462219*(x - 13)3)/9007199254740992 + 15388685044907045/562949953421312在时的函数表达式为S =(49360

35、5447569857*(x - 13)3)/72057594037927936 - (550210270623049*x)/562949953421312 - (5376237321301877*(x - 14)3)/54043195528445952 + 12850079569968551/562949953421312在时的函数表达式为S =13944784877622509/562949953421312 - (857131708001443*(x - 14)3)/6755399441055744 - (493605447569857*(x - 15)3)/720575940379279

36、36 - (314202153442023*x)/281474976710656在时的函数表达式为S =(2888372562541421*(x - 15)3)/9007199254740992 - (447606859642455*x)/1125899906842624 + (857131708001443*(x - 16)3)/6755399441055744 + 15751567443360463/1125899906842624在时的函数表达式为S =(1223276603252857*x)/1125899906842624 + (3896133174998587*(x - 16)3)

37、/27021597764222976 - (2888372562541421*(x - 17)3)/9007199254740992 - 10982567962964529/1125899906842624在时的函数表达式为S =(637035986695777*x)/562949953421312 - (1348101747182487*(x - 17)3)/36028797018963968 - (3896133174998587*(x - 18)3)/27021597764222976 - 5923044627661189/562949953421312在时的函数表达式为S =(2553

38、25723948709*x)/281474976710656 + (614517988558921*(x - 18)3)/13510798882111488 + (1348101747182487*(x - 19)3)/36028797018963968 - 3648122929290727/562949953421312在时的函数表达式为S =(157625986957967*x)/562949953421312 - (1412101843757803*(x - 19)3)/13510798882111488 - (614517988558921*(x - 20)3)/p>

39、111488 + 1529680414279421/281474976710656最后得到用三次样条插值拟合所测数据点的曲线如下列图所示。由此可以看出三次样条插值在节点处的光滑性得到了改善，但是要使用方程组计算增大了计算量，同时还要附加边界条件，分段三次样条插值对图形的控制能力还不够灵活。图中的圆圈是取的基点值，从图中可以看到三次样条插值到达较好的效果。2.3复化Simpson公式近似计算光缆长度铺设河底光缆的曲线图如上图所示。本实验采用复化Simpson公式近似计算光缆长度。2.3.1算法组织1、算法依据及思想由以上三次样条插值结果可以得到在20个区间上的是三次插值函数，结果如上述运行结果所

40、示。此题要求近似计算光缆长度，即在每个区间上采用弧长积分公式，便可得到所在区间的近似光缆长度。各个区间的近似光缆长度之和即为所求近似计算光缆长度。其中为对应区间上的三次插值函数。利用复化Simpson公式近似计算光缆长度的根本思想如下。利用定积分对区间的区间可加性，将求积公式分别应用于每一个小区间上，设是区间上a,b上满足的n+1个点，由于。取等距地分布在区间a,b上，此时，=,.。得到复化梯形公式为：误差为：得到复化simpson公式为：误差为：此题采用在每个区间上20等分的方法，近似求解各个区间的定积分。2、算法结构首先本实验将所在的每个区间分别进行20等分，然后采用如下算法结构：2.3.

41、2算法的Matlab实现为了实现在每个区间积分上的简便，复化Simpson公式的主函数与三次样条插值函数一起调用。复化Simpson公式的子函数的Matlab程序为：function simpson=simpson(g,a,b,n)%Simpson积分公式子程序f=inline(g);simpson1=0;simpson2=0;for i=1:n-1; simpson1=simpson1+f(a+(2*i-1)*(b-a)/(2*n); simpson2=simpson2+f(a+i*(b-a)/n);endsimpson=(b-a)/(6*n)*(f(a)+f(b)+4*simpson1+2

42、*simpson2);end2.3.3实验结果及分析针对上述Matlab程序，运行结果如下：sl=25.4533，即为所求近似计算光缆长度。复化Simpson求积公式，既简单又使用，使用时必须选取适当的步长，以保证计算的精度。本实验的步长直接选取，在实际应用中可由复化公式的余项来估计，但最好用计算过程中自动选取步长的方法。为了提高计算精度，可以采用Romberg积分法与自适应积分法。3.实验三3.1实验题目假定某天的气温变化记录如下表所示，试用数据拟合的方法找出这一天的气温变化的规律;试计算这一天的平均气温,并试估计误差。时刻0123456789101112平均气温15141414141516

43、182020232528时刻131415161718192021222324平均气温3134312927252422201817163.2算法组织给定数据点和一组函数。求数(假定)，使函数满足为计算方便，只需考虑最小二乘近似问题的解应使上式成立；在中，它是ak的二次函数，因此其取极小值的必要条件为，即经过推导知，假设记那么最小二乘问题的法方程为解此法方程，就可以得到各项系数。继续计算得到：把称为最小二乘拟合多项式的平方误差，记作因此，可以有两种方法求解1、求解法方程法直接进行求解解方程，解法比拟简单，但是一方面运算量较大，并且由于条件数的影响可能会有较大误差。本设计采用了高斯消去法，减少了误差

44、。2、正交分解法对生成的矩阵进行QR分解，可写成矩阵形式，利用矩阵的QR分解，可以获得解最小二乘问题可靠的解。具体实现程序见Matlab程序实现。3.3算法的Matlab实现3.3.1法方程最小二乘法%法方程法最小二乘法clcclear% x=input(请输入横坐标x=:)% y=input(请输入纵坐标y=:)x=0:24;y=15 14 14 14 14 15 16 18 20 20 23 25 28 31 34 31 29 27 25 24 22 20 18 17 16;s=0;n=length(x);%求平均温度Tavfor i=1:n s=s+y(i);endTav=s/n%使用一

45、次函数、二次函数、三次函数至八次函数的最小二乘法拟合%利用法方程法，其中A=G*G，b=G*yfor m=1:8A=zeros(m+1,m+1);for i=0:mfor j=0:mA(i+1,j+1)=sum(x.(i+j);endbb(i+1)=sum(x.i.*y);endaa=Abb; %左右翻转disp(拟合多项式如下及相应的系数)p=fliplr(aa)Px=poly2sym(p)xi=0:0.1:24;yi=polyval(p,xi); y1=polyval(p,x);disp(估算误差)e=sum(y-y1).2)subplot(2,4,m);plot(x,y,.,xi,yi,

46、r);title(最小二乘法近似);end3.3.2QR分解为了有效的看出正交分解的规律，报告依次使用二次函数、三次函数、四次函数、五次函数的最小二乘法拟合。1、正交分解法主程序%正交方法最小二乘法clcclear%x=input(请输入横坐标x=:)%y=input(请输入纵坐标y=:)x=0:24;y=15 14 14 14 14 15 16 18 20 20 23 25 28 31 34 31 29 27 25 24 22 20 18 17 16;%依次使用二次函数、三次函数、四次函数、五次函数的最小二乘法拟合%二次函数%-disp(拟合二次函数如下及相应的系数)z=G_wlg(x,y,

47、2);zz=QG_wlg(z);%q,r=qr(z)aa=jief_wlg(zz);p2=fliplr(aa)Px2=poly2sym(p2)xi=0:0.1:24;yi=polyval(p2,xi); y1=polyval(p2,x);disp(估算误差)e2=wucha_wlg(zz)subplot(2,2,1);plot(x,y,.,xi,yi,r);title(二次函数);%-%三次函数%-disp(拟合三次函数如下及相应的系数)z=G_wlg(x,y,3);zz=QG_wlg(z);%q,r=qr(z)aa=jief_wlg(zz);p3=fliplr(aa)Px3=poly2sym

48、(p3)xi=0:0.1:24;yi=polyval(p3,xi); y1=polyval(p3,x);disp(估算误差)e3=wucha_wlg(zz)subplot(2,2,2);plot(x,y,.,xi,yi,r);title(三次函数);%-%四次函数%-disp(拟合四次函数如下及相应的系数)z=G_wlg(x,y,4);zz=QG_wlg(z);%q,r=qr(z)aa=jief_wlg(zz);p4=fliplr(aa)Px4=poly2sym(p3)xi=0:0.1:24;yi=polyval(p4,xi); y1=polyval(p4,x);disp(估算误差)e4=wu

49、cha_wlg(zz)subplot(2,2,3);plot(x,y,.,xi,yi,r);title(四次函数);%-%四次函数%-disp(拟合五次函数如下及相应的系数)z=G_wlg(x,y,5);zz=QG_wlg(z);%q,r=qr(z)aa=jief_wlg(zz);p5=fliplr(aa)Px5=poly2sym(p5)xi=0:0.1:24;yi=polyval(p5,xi); y1=polyval(p5,x);disp(估算误差)e5=wucha_wlg(zz)subplot(2,2,4);plot(x,y,.,xi,yi,r);title(五次函数);2、G矩阵求解函数

50、%G矩阵function G = G_wlg(x,y,n)%得到G矩阵，最后一列存放y值，n表示拟合次数，1,x,x2.m=length(x);G=zeros(m,n+2);for i=1:mfor j=1:n+1 G(i,j)=(x(i)(j-1);endendfor i=1:m G(i,n+2)=y(i);end3、RO与h1、h2矩阵函数%RO与h1、h2矩阵function GG = QG_wlg(G)%依次变换得到Gk，Qk，最终得到RO与h1、h2m,n=size(G);n=n-1;w=zeros(m,1);for k=1:n %依次变换得到Qk gik=0; for i=k:m

51、gik=gik+G(i,k)2; end dlt=-(sign(G(k,k)*sqrt(gik); w(k)=G(k,k)-dlt; for j=k+1:m w(j)=G(j,k); end blt=dlt*w(k); %依次变换得到Gk G(k,k)=dlt; %k取值 for j=k+1:n+1 wg=0; for ii=k:m wg=wg+w(ii)*G(ii,j); end t=wg/blt; for i=k:m G(i,j)=G(i,j)+t*w(i); end endendGG=G;4、求解系数函数%回代法求解方程组，求解系数function x = jief_wlg(GG) m,

52、n=size(GG); n=n-1; x=zeros(n,1); x(n)=GG(n,n+1)/GG(n,n); for i=n-1:-1:1 gx=0; for j=i+1:n gx=gx+GG(i,j)*x(j); end x(i)=(GG(i,n+1)-gx)/GG(i,i); end5、误差函数%误差function rou = wucha_wlg(GG) m,n=size(GG); n=n-1; rou=0; for i=n+1:m rou=rou+GG(i,n+1)2;end%误差3.4实验结果及分析3.4.1法方程最小二乘法输出结果这一天平均温度为Tav =21.2000拟合一次

53、多项式如下及相应的系数p =0.3554 16.9354Px =(231*x)/650 + 5504/325估算误差e =753.8123拟合二次多项式如下及相应的系数p =-0.0938 2.6064 8.3063Px =(5869189440802493*x)/2251799813685248 - (1262*x2)/13455 + 24296/2925估算误差e =280.3395拟合三次多项式如下及相应的系数p =-0.0084 0.2075 -0.2273 13.3880Px =- (4824449974797021*x3)/576460752303423488 + (6142*x2

54、)/29601 - (4095195769401027*x)/18014398509481984 + 7832/585估算误差e =131.0618拟合四次多项式如下及相应的系数p =0.0009 -0.0532 0.8909 -3.7050 16.7939Px =(4311222046072087*x4)/4611686018427387904 - (959118195350923*x3)/18014398509481984 + (8024210619242803*x2)/9007199254740992 - (8342830140798239*x)/2251799813685248 + 5

55、90881761455037/35184372088832估算误差e =59.0412拟合五次多项式如下及相应的系数p =0.0001 -0.0045 0.0623 -0.1202 -0.4997 14.9504Px =(3360357177283195*x5)/36893488147419103232 - (5222864195889131*x4)/1152921504606846976 + (4487849523002499*x3)/72057594037927936 - (8657808537992275*x2)/72057594037927936 - (4500908134762693

56、*x)/9007199254740992 + 8416302560899915/562949953421312估算误差e =33.1446拟合六次多项式如下及相应的系数p =-0.0000 0.0005 -0.0160 0.2069 -0.9506 1.2804 14.3314Px =- (6967708126795581*x6)/1180591620717411303424 + (2379712557858713*x5)/4611686018427387904 - (2311817739804005*x4)/144115188075855872 + (3726966007444799*x3)

57、/18014398509481984 - (4280910250199293*x2)/4503599627370496 + (5766485188992375*x)/4503599627370496 + 8067842731307123/562949953421312估算误差e =29.1139拟合七次多项式如下及相应的系数p =-0.0000 0.0001 -0.0035 0.0507 -0.3597 1.3563 -2.2558 15.0741Px =- (6901230370958207*x7)/4722366482869645213696 + (8622383106910967*x6)

58、/73786976294838206464 - (4084605813955375*x5)/1152921504606846976 + (7309441973132021*x4)/144115188075855872 - (6479824387411487*x3)/18014398509481984 + (6108410266124833*x2)/4503599627370496 - (2539858743572183*x)/1125899906842624 + 8485987341360957/562949953421312估算误差e =20.1716拟合八次多项式如下及相应的系数p =0.

59、0000 -0.0000 0.0003 -0.0083 0.1043 -0.6909 2.3648 -3.4106 15.2169Px = (2341429260780705*x8)/37778931862957161709568 - (4374797679237899*x7)/590295810358705651712 + (3220983078808941*x6)/9223372036854775808 - (4772103004746209*x5)/576460752303423488 + (1879800419770211*x4)/18014398509481984 - (311169

60、7315246825*x3)/4503599627370496 + (41601748837387*x2)/17592186044416 - (7679972061905371*x)/2251799813685248 + 8566378569142467/562949953421312 估算误差e =19.6054Matlab输出图形如下图3.4.2QR分解输出结果拟合二次函数如下及相应的系数p2 =-0.0938 2.6064 8.3063Px2 =(5869189440802495*x)/2251799813685248 - (1262*x2)/13455 + 24296/2925估算误差