计算方法上机报告(姚威)

1、计算方法B上机作业报告 姓名：姚威 班级：硕5081班 学号：3115312022 学科专业：电气工程任课教师：苏剑老师2015年12月目录1. 上机题一11.1 解题源程序11.2 计算结果31.3 实验结果分析及总结32. 上机题二32.1 三次样条插值拟合数据点42.1.1算法组织42.1.2 算法52.1.3算法的Matlab实现62.1.4 实验结果及分析92.2 复化Simpson公式近似计算光缆长度122.2.1算法组织122.2.2算法的Matlab实现132.2.3实验结果及分析143. 上机题三143.1 算法思想143.2 算法组织143.3 算法的Matlab实现153

2、.4 实验结果及分析173.4.1法方程最小二乘法输出结果173.4.2 基于最小二乘法ployfit函数法193.5 QR分解243.5.1 结果分析294. 上机题四304.1 二分法算法组织304.1.1算法思想及依据304.1.2算法结构304.2 二分法算法的Matlab实现304.2.1利用Matlab校验根的区间304.2.2使用二分法Matlab程序314.3 二分法实验结果及分析324.4 割线法算法组织334.4.1算法思想及依据334.4.2算法结构334.5 割线法算法的Matlab实现334.5.1割线法实验结果及分析345. 上机题五345.1 算法思想355.2

3、算法组织355.2.1算法思想及依据355.2.2算法结构355.3 算法的Matlab实现365.3.1非压缩带状对角方程组365.3.2压缩带状对角方程组385.4 实验结果及分析405.4.1 Matlab运行结果405.4.2 Matlab运行结果分析435.5 实际问题445.6 分析总结476. 总结及致谢481. 上机题一 对以下和式计算：，要求： (1)若只需保留11个有效数字，该如何进行计算；(2)若要保留30个有效数字，则又将如何进行计算；计算方法任何一个进制、具有t位有效数字的实数总可以表示成，其中，是满足的整数。有效数字是指一个近似数的“有意义”的数字的数位，通常在十进

4、制中讨论，设,其中，近似数，若，则称是的具有t位有效数字的近似数，或称具有t位有效数字。根据有效数字的要求可以得到计算公式的精度要求。再利用后验误差估计式可以根据精度要求得到和式累加的项数。首先满足题目要求的情况下，利用后验误差估计计算出算法中的最大运行次数。即以一下条件作为迭代终止的准则：，其中为误差数值。该实验第一问利用后验误差估计截取的项数，使误差小于1e-10；第二问利用后验误差估计截取的项数，使误差小于1e-30。得出运行次数，作为累加的终止条件。特别地，为了减小舍入误差在计算S时所采用的方法是逆序相加，其依据是：S中各项的绝对值为递减的正数，而两个数量级相差较大的数字相加减时，较小

5、数的有效数字会丢失，从而导致最后的运算结果失真。为避免“大数吃小数”现象的发生，采用逆序相加；同时为了避免相近数相减和减小计算步骤，我们将化简为进行计算。另外，要限定计算结果的有效数字的个数，在MATLAB软件中可直接调用函数vpa(变量名,位数)来控制计算结果的位数，最终完成精度要求较高的计算。1.1 解题源程序function solve1clc;clear;sum1=0;sum2=0;n=0;N=0;%第一小问%根据精度要求利用后验误差式，估计要使误差小于，需要累加的项数n，所以要保证计算精度第n+1项要小于，才能保证计算的精度while abs(960*n2+1208*n+376)/(

6、8*n+1)*(8*n+4)*(8*n+5)*(8*n+6)/(16.n)>10(-11)%为了减小误差，我们对计算公式进行变形，一是防止相近数相减，二是为了减少计算步骤； n=n+1;endn;%累加进行的项数disp('累加的次数:')n%为减少舍入误差的影响，按绝对值递增的顺序求和。因此按n按从小到大的顺序相加,步距取为-1for i=n:-1:0 sum1=sum1+(960*i.2+1208*i+376)/(8*i+1)*(8*i+4)*(8*i+5)*(8*i+6)/(16.i);enddisp('保留11位有效数字时的近似值:')Sum1=v

7、pa(sum1,11)%第二小问%根据精度要求利用后验误差式，估计要使误差小于0.5e-30需要累加的项数NWhileabs(960*N2+1208*N+376)/(8*N+1)*(8*N+4)*(8*N+5)*(8*N+6)/(16.N)>10(-30) N=N+1;endN;%累加进行的项数disp('累加的次数:')N%为减少舍入误差的影响，按绝对值递增的顺序求和。因此按N按从小到大的顺序相加,步距取为-1for j=N:-1:0 sum2=sum2+(960*j.2+1208*j+376)/(8*j+1)*(8*j+4)*(8*j+5)*(8*j+6)/(16.j

8、);enddisp('保留30位有效数字时的近似值:')Sum2=vpa(sum2,30)End1.2 计算结果 累加的次数:n = 8 保留11位有效数字时的近似值: sum1 =3.1415926536 累加的次数:N = 23保留30位有效数字时的近似值:sum2 =3.141592653589793238462643383281.3 实验结果分析及总结通过运行上述程序可知：本题通过计算可以发现，计算值和的值非常接近，由，可知这里的,当t=10时，我们为了得到有效数位为11，所以计算的最后一项t应该等于t=-11。同理，为了得到有效数位为30，所以计算的最后一项t应该等于

9、t=-30。第一问中，在满足误差小于1e-10时，截取的项数为8，保留11位有效数字估算值:3.1415926536 。第二问中，在满足误差小于1e-30时，截取的项数为23，保留30位有效数字估算值:3.14159265358979323846264338328。从上述的算法思想中可以看出，浮点数运算中不仅要满足误差要求，还要尽可能地减少计算量，而其计算量的大小很大程度上是由算式要求的精度来决定的；此外还要考虑舍入误差的影响，这时就要对所运算数据的性质进行分析，设置合适的算法，从而提高运算的精度。2. 上机题二某通信公司在一次施工中，需要在水面宽度为20米的河沟底部沿直线走向铺设一条沟底光缆

10、。在铺设光缆之前需要对沟底的地形进行初步探测，从而估计所需光缆的长度，为工程预算提供依据。已探测到一组等分点位置的深度数据(单位:米)如下表所示：分点0123456深度9.018.967.967.978.029.0510.13分点78910111213深度11.1812.2613.2813.3212.6111.2910.22分点141516171819 20深度9.157.907.958.869.8110.8010.93(1)请用合适的曲线拟合所测数据点；(2)预测所需光缆长度的近似值，并作出铺设河底光缆的曲线图；2.1 三次样条插值拟合数据点2.1.1算法组织使用分段插值多项式拟合数据的曲线

11、没有足够的光滑性，所以此题采用三次样条函数拟合所测数据点，即在每一个子区间上是一个三次多项式，这些多项式的曲线在节点处“光滑”地连接起来。这样的分段三次多项式称为三次样条函数。样条函数为：记则上式可写成：这是一个有n-1个方程，n+1个未知数的线性方程组，还需要两个边界条件才能确定唯一解。（1）自然样条（2）指定端点的一阶导数。此时得（3）若两端点的导数值均不知道时，可以通过推导得到最后可以得到这些方程组的统一形式本题采用边界条件三求解，其中2.1.2 算法SPLINEM()本算法计算三次样条的参数Mi，存放于数组M中，参数根据不同边界条件给定。 1.1 2.1 2.1.1 3. 4. 4.1

12、 4.2 4.3 用算法TSS(a,b,c,M,n,M)算法TTS（a，b，c，d，n，x）本算法解对角线性方程组，系数存于数组a、b和c中，d是有段向量。2.1.3算法的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=zeros(2,N);for i=1

13、: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(i,5)=(D(i,4)-D

14、(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);end%利用TSS方法解方程

15、组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.001:xx(i+1); sy=Sy(sx); str=repmat('(',21,1) num2str(x') repmat(',',21,1) num2str(y'); plot(x,y,'o') %原探测数据点及数值 text(x,-y,cellstr(str) stem(x,-y,'fille

16、d','LineWidth',3)%花火柴杆图 plot(sx,sy,'r','LineWidth',3) %按插值函数表达式作图 hold on end%添加水平线l=line(0,20,0,0);%设置直线的宽度和颜色set(l,'Linewidth',1.5)set(l,'color','k')%设置坐标轴刻度set(gca,'Xtick',0:1:20)%设置坐标轴刻度set(gca,'Ytick',-15:2:15)%添加图形标题和坐标轴名称titl

17、e('作出的铺设河底光缆的曲线图')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; v=i; sl=sl+simpson(g,u,v,20);endSy;g;disp('所需光缆长度的近似值为: ');slTSS解方程组的

18、子程序%追赶法解方程组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)-c(i)*x(i+1)/u(i);endx;end得到子区间函数表达式的子程序function S

19、y=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)2*m(i-1)/6)*(xx(i)-x)/h(i-1)+(y(i)-h(i-1)2*m(i)/

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

21、41055744 - (7514075829091499*(x - 2)3)/54043195528445952 - (277217198887389*x)/562949953421312在时的函数表达式为S =(650311338860789*(x - 2)3)/1688849860263936 - (418602078066937*x)/562949953421312 + (2489533575770641*(x - 3)3)/6755399441055744 + 2762873458340869/281474976710656在时的函数表达式为S =(2624680819031*x)/4

22、398046511104 - (1092494410535677*(x - 3)3)/6755399441055744 - (650311338860789*(x - 4)3)/1688849860263936 + 3262063247973023/562949953421312在时的函数表达式为S =(79721223946041*x)/140737488355328 + (8155793057367127*(x - 4)3)/27021597764222976 + (1092494410535677*(x - 5)3)/6755399441055744 + 3330360244180239

23、/562949953421312在时的函数表达式为S =(814815571273233*x)/562949953421312 - (3544057556774495*(x - 5)3)/54043195528445952 - (8155793057367127*(x - 6)3)/27021597764222976 + 425353433367447/281474976710656在时的函数表达式为S =(274137988850561*x)/281474976710656 + (566803888785965*(x - 6)3)/54043195528445952 + (354405755

24、6774495*(x - 7)3)/54043195528445952 + 306243053520945/70368744177664在时的函数表达式为S =(617478217955525*x)/562949953421312 - (1377815456891219*(x - 7)3)/216172782113783808 - (566803888785965*(x - 8)3)/54043195528445952 + 1965528746386739/562949953421312在时的函数表达式为(545284355633681*x)/562949953421312 + (405384

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

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

27、x)/562949953421312 + 14925062697634997/562949953421312在时的函数表达式为S =(5376237321301877*(x - 12)3)/54043195528445952 - (745487614849087*x)/562949953421312 + (1394059081462219*(x - 13)3)/9007199254740992 + 15388685044907045/562949953421312在时的函数表达式为S =(493605447569857*(x - 13)3)/72057594037927936 - (55021

28、0270623049*x)/562949953421312 - (5376237321301877*(x - 14)3)/54043195528445952 + 12850079569968551/562949953421312在时的函数表达式为S =13944784877622509/562949953421312 - (857131708001443*(x - 14)3)/6755399441055744 - (493605447569857*(x - 15)3)/72057594037927936 - (314201753442023*x)/281474976710656在时的函数表达式

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

30、)/9007199254740992 - 10982567962964529/1125899906842624在时的函数表达式为S =(637035986695777*x)/562949953421312 - (1348101747182487*(x - 17)3)/36028797018963968 - (3896133174998587*(x - 18)3)/27021597764222976 - 5923044627661189/562949953421312在时的函数表达式为S =(255325723948709*x)/281474976710656 + (61451798855892

31、1*(x - 18)3)/13510798882111488 + (1348101747182487*(x - 19)3)/36028797018963968 - 3648122929290727/562949953421312在时的函数表达式为S =(157625986957967*x)/562949953421312 - (1412101843757803*(x - 19)3)/13510798882111488 - (614517988558921*(x - 20)3)/13510798882111488 + 1529680414279421/281474976710656所需光缆长度的

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

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

34、上的简便，复化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*simpson2);end2.2.3实验结果

35、及分析针对上述Matlab程序，运行结果如下：sl=25.4533，即为所求近似计算光缆长度。拟合曲线足够光滑，说明采用三样条插值多项式拟合效果很好，误差较小，拟合曲线比较接近于离散点的分布图像;复化Simpson求积公式，既简单又使用，使用时必须选取适当的步长，以保证计算的精度。本实验的步长直接选取，在实际应用中可由复化公式的余项来估计，但最好用计算过程中自动选取步长的方法。为了提高计算精度，可以采用Romberg积分法与自适应积分法。3. 上机题三 假定某天的气温变化记录如下表所示，试用数据拟合的方法找出这一天的气温变化的规律;试计算这一天的平均气温,并试估计误差。时刻0123456789

36、101112平均气温15141414141516182020232528时刻131415161718192021222324平均气温3134312927252422201817163.1 算法思想本题可以用多项式插值方法作数据近似拟合，也可以用插值样条函数拟合数据，当数据点数目较大时，要用很高的多项式，会带来计算上的困难，误差大,用插值样条函数，虽然数据性质好，计算量也不大，但是存放参数 M i 的量很大。本题数据量较少，因此在本题中采用最小二乘近似拟合的方法即可得到较好的拟合效果。3.2 算法组织给定数据点和一组函数。求数(假定)，使函数满足为计算方便，只需考虑最小二乘近似问题的解应使上式成

37、立；在中，它是ak的二次函数，因此其取极小值的必要条件为，即经过推导知，若记则最小二乘问题的法方程为解此法方程，就可以得到各项系数。继续计算得到：把称为最小二乘拟合多项式的平方误差，记作因此，可以有两种方法求解1、求解法方程法直接进行求解解方程，解法比较简单，但是一方面运算量较大，并且由于条件数的影响可能会有较大误差。本设计采用了高斯消去法，减少了误差。2、正交分解法对生成的矩阵进行QR分解，可写成矩阵形式，利用矩阵的QR分解，可以获得解最小二乘问题可靠的解。具体实现程序见Matlab程序实现。3.3 算法的Matlab实现法一：%法方程法最小二乘法clcclear% x=input('

38、;请输入横坐标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%使用一次函数、二次函数、三次函数至八次函数的最小二乘法拟合%利用法方程法，其中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)=

39、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,'r');title('最小二乘法近似');endsubplot(2,4,1);grid on

40、title('最小二乘法一阶拟合近似曲线');subplot(2,4,2);grid ontitle('最小二乘法二阶拟合近似曲线');subplot(2,4,3);grid ontitle('最小二乘法三阶拟合近似曲线');subplot(2,4,4);title('最小二乘法四阶拟合近似曲线');grid onsubplot(2,4,5);grid ontitle('最小二乘法五阶拟合近似曲线');subplot(2,4,6);grid ontitle('最小二乘法六阶拟合近似曲线');subp

41、lot(2,4,7);grid ontitle('最小二乘法七阶拟合近似曲线');subplot(2,4,8);grid ontitle('最小二乘法八阶拟合近似曲线');3.4 实验结果及分析3.4.1法方程最小二乘法输出结果这一天平均温度为Tav =21.2000拟合一次多项式如下及相应的系数p =0.3554 16.9354Px =(231*x)/650 + 5504/325估算误差e =753.8123拟合二次多项式如下及相应的系数p =-0.0938 2.6064 8.3063Px =(5869189440802493*x)/2251799813685

42、248 - (1262*x2)/13455 + 24296/2925估算误差e =280.3395拟合三次多项式如下及相应的系数p =-0.0084 0.2075 -0.2273 13.3880Px =- (4824449974797021*x3)/576460752303423488 + (6142*x2)/29601 - (4095195769401027*x)/18014398509481984 + 7832/585估算误差e =131.0618拟合四次多项式如下及相应的系数p =0.0009 -0.0532 0.8909 -3.7050 16.7939Px =(431122204607

43、2087*x4)/4611686018427387904 - (959118195350923*x3)/18014398509481984 + (8024210619242803*x2)/9007199254740992 - (8342830140798239*x)/2251799813685248 + 590881761455037/35184372088832估算误差e =59.0412拟合五次多项式如下及相应的系数p =0.0001 -0.0045 0.0623 -0.1202 -0.4997 14.9504Px =(3360357177283195*x5)/36893488147419

44、103232 - (5222864195889131*x4)/1152921504606846976 + (4487849523002499*x3)/72057594037927936 - (8657808537992275*x2)/72057594037927936 - (4500908134762693*x)/9007199254740992 + 8416302560899915/562949953421312估算误差e =33.1446拟合六次多项式如下及相应的系数p =-0.0000 0.0005 -0.0160 0.2069 -0.9506 1.2804 14.3314Px =- (

45、6967708126795581*x6)/1180591620717411303424 + (2379712557858713*x5)/4611686018427387904 - (2311817739804005*x4)/144115188075855872 + (3726966007444799*x3)/18014398509481984 - (4280910250199293*x2)/4503599627370496 + (5766485188992375*x)/4503599627370496 + 8067842731307123/562949953421312估算误差e =29.11

46、39拟合七次多项式如下及相应的系数p =-0.0000 0.0001 -0.0035 0.0507 -0.3597 1.3563 -2.2558 15.0741Px =- (6901230370958207*x7)/4722366482869645213696 + (8622383106910967*x6)/73786976294838206464 - (4084605813955375*x5)/1152921504606846976 + (7309441973132019*x4)/144115188075855872 - (6479824387411487*x3)/1801439850948

47、1984 + (6108410266124833*x2)/4503599627370496 - (2539858743572183*x)/1125899906842624 + 8485987341360957/562949953421312估算误差e =20.1716拟合八次多项式如下及相应的系数p =0.0000 -0.0000 0.0003 -0.0083 0.1043 -0.6909 2.3648 -3.4106 15.2169Px = (2341429260780705*x8)/37778931862957161709568 - (4374797679237899*x7)/590295

48、810358705651712 + (3220983078808941*x6)/9223372036854775808 - (4772103004746209*x5)/576460752303423488 + (1879800419770211*x4)/18014398509481984 - (3111697315246825*x3)/4503599627370496 + (41601748837387*x2)/17592186044416 - (7679972061905371*x)/2251799813685248 + 8566378569142467/562949953421312 估算

49、误差e =19.6054Matlab输出图形如图所示法二：多项式拟合使用一个多项式逼近一组给定的数据，是数据分析上常用的方法，使用polyfit函数实现。拟合的准则是最小二乘法，即找出使最小的。3.4.2 基于最小二乘法ployfit函数法 %法方程法最小二乘法clearclcx=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;for m=1:4p=polyfit(x,y,m);p1=polyval(p,x);y1=polyval(p,x);xi=0:0.01:24;yi=pol

50、yval(p,xi);figure(1)subplot(2,2,m);plot(xi,yi);plot(x,y,'.',xi,yi,'r','LineWidth',2);grid ondisp('各阶估算误差')e(m)=sum(y-y1).2);end subplot(2,2,1);title('最小二乘法一阶拟合近似曲线');subplot(2,2,2);title('最小二乘法二阶拟合近似曲线');subplot(2,2,3);title('最小二乘法三阶拟合近似曲线');su

51、bplot(2,2,4);title('最小二乘法四阶拟合近似曲线'); for m=5:8p=polyfit(x,y,m);p1=polyval(p,x);y1=polyval(p,x);xi=0:0.01:24;yi=polyval(p,xi);figure(2)m=m-4;subplot(2,2,m);plot(xi,yi);plot(x,y,'.',xi,yi,'r','LineWidth',2);grid ondisp('各阶估算误差')e(m+4)=sum(y-y1).2);end subplot(2,

52、2,1);title('最小二乘法五阶拟合近似曲线');subplot(2,2,2);title('最小二乘法六阶拟合近似曲线');subplot(2,2,3);title('最小二乘法七阶拟合近似曲线');subplot(2,2,4);title('最小二乘法八阶拟合近似曲线'); for m=9:12p=polyfit(x,y,m)p1=polyval(p,x);y1=polyval(p,x);xi=0:0.01:24;yi=polyval(p,xi);figure(3)m=m-8;subplot(2,2,m);plot(xi

53、,yi);plot(x,y,'.',xi,yi,'r','LineWidth',2);grid ondisp('各阶估算误差')e(m+8)=sum(y-y1).2);end subplot(2,2,1);title('最小二乘法九阶拟合近似曲线');subplot(2,2,2);title('最小二乘法十阶拟合近似曲线');subplot(2,2,3);title('最小二乘法十一阶拟合近似曲线');subplot(2,2,4);title('最小二乘法十二阶拟合近似曲线&

54、#39;); for m=13:16p=polyfit(x,y,m);p1=polyval(p,x);y1=polyval(p,x);xi=0:0.01:24;yi=polyval(p,xi);figure(4)m=m-12;subplot(2,2,m);plot(xi,yi);plot(x,y,'.',xi,yi,'r','LineWidth',2);grid ondisp('各阶估算误差')e(m+12)=sum(y-y1).2);end subplot(2,2,1);title('最小二乘法十三阶拟合近似曲线');subplot(2,2,2);title('最小二乘法十四阶拟合近似曲线');subplot(2,2,3);title('最小二乘法十五阶拟合近似曲线');sub