《计算方法上机作业》由会员分享，可在线阅读，更多相关《计算方法上机作业（34页珍藏版）》请在金锄头文库上搜索。
1、计算方法上机报告姓 名：学 号：班 级：上课班级：计算方法上机报告1说明：本次上机实验使用的编程语言是 Matlab 语言，编译环境为 MATLAB 7.11.0，运行平台为 Windows 7。1. 对以下和式计算： ，要求： 0 681548216nnSn 若只需保留 11 个有效数字，该如何进行计算； 若要保留 30 个有效数字，则又将如何进行计算；（1） 算法思想1、根据精度要求估计所加的项数，可以使用后验误差估计，通项为：；421146818568n nan 2、为了保证计算结果的准确性，写程序时，从后向前计算；3、使用 Matlab 时，可以使用以下函数控制位数 ：digits(位
2、数)或 vpa(变量，精度为数)（2）算法结构1. ;0s；681548216nntn2. for 0,iif mtend;3. for ,12,0ni;st计算方法上机报告2（3）Matlab 源程序clear; %清除工作空间变量clc; %清除命令窗口命令m=input(请输入有效数字的位数 m=); %输入有效数字的位数s=0; for n=0:50t=(1/16n)*(4/(8*n+1)-2/(8*n+4)-1/(8*n+5)-1/(8*n+6);if t0;sgn=1;elseif G(k,k)=0;sgn=0;else sgn=-1;endsgm=-sgn*sqrt(sum(G(
3、k:m,k).2);w=zeros(1,n); w(k)=G(k,k)-sgm;for j=k+1:mw(j)=G(j,k);endbt=sgm*w(k);G(k,k)=sgm; %变换Gk-1到Gkfor j=k+1:n+1t=sum(w(k:m)*G(k:m,j)/bt;for i=k:m;G(i,j)=G(i,j)+t*w(i);end计算方法上机报告14endend A (n)=G(n,n+1)/G(n,n); %解三角方程求系数Afor i=n-1:-1:1A (i)=(G(i,n+1)-sum(G(i,i+1:n).*A (i+1:n)/G(i,i);end e=sum(G(n+1
4、:m,n+1).2); %计算误差efprintf(%d次函数的系数是： ,h); %输出系数a及误差edisp(A);fprintf(使用%d次函数拟合的误差是%f：,h,e);t=0:0.05:24;A=fliplr(A); %将系数数组左右翻转 Y=poly2sym(A); %将系数数组转化为多项式subs(Y,x,t);Y=double(ans);figure(1)plot(x,y,k*,t,Y,r-); %绘制拟合多项式函数图形xlabel(时刻); %标注坐标轴含义ylabel(平均气温);title(num2str(n-1),次函数的最小二乘曲线);grid;% 指数函数的最小二
5、乘近似计算方法上机报告15yy=log(y);n=3;G=;GG=;for j=0:(n-1)g=x.j; %g(x)按列排列G=vertcat(G,g); %g垂直连接Ggg=t.j; %g(x)按列排列GG=vertcat(GG,gg); %g垂直连接GendG=G; %转置得到矩阵Gfor i=1:m %将数据y作为G 的最后一列（n+1列）G(i,n+1)=yy(i);endG; for k=1:n %形成矩阵Q(k)if G(k,k)0;sgn=1;elseif G(k,k)=0;sgn=0;else sgn=-1;计算方法上机报告16endsgm=-sgn*sqrt(sum(G(k
6、:m,k).2);w=zeros(1,n); w(k)=G(k,k)-sgm;for j=k+1:mw(j)=G(j,k);endbt=sgm*w(k);G(k,k)=sgm; %变换Gk-1到Gkfor j=k+1:n+1t=sum(w(k:m)*G(k:m,j)/bt;for i=k:m;G(i,j)=G(i,j)+t*w(i);endendend A(n)=G(n,n+1)/G(n,n); %解三角方程求系数Afor i=n-1:-1:1A (i)=(G(i,n+1)-sum(G(i,i+1:n).*A (i+1:n)/G(i,i);end b=-A(3);c=A(2)/(2*b);计算
7、方法上机报告17a=exp(A(1)+b*(c2);G(n+1:m,n+1)=exp(sum(G(n+1:m,n+1).2);e=sum(G(n+1:m,n+1).2); %计算误差efprintf(n指数函数的系数是：a=%f,b=%f,c=%f,a,b,c); %输出系数及误差efprintf(n使用指数函数拟合的误差是：%f,e);t=0:0.05:24;YY=a.*exp(-b.*(t-c).2);figure(2)plot(x,y,k*,t,YY,r-); %绘制拟合指数函数图形xlabel(时刻); %标注坐标轴含义ylabel(平均气温);title(指数函数的最小二乘曲线);g
8、rid;（4）结果与分析a、二次函数：一天的平均气温为： 21.20002 次函数的系数: 8.3063 2.6064 -0.0938使用 2 次函数拟合的误差是：280.339547二次函数的最小二乘曲线如下图所示：计算方法上机报告18b、三次函数：一天的平均气温为： 21.20003 次函数的系数: 13.3880 -0.2273 0.2075 -0.0084使用 3 次函数拟合的误差是： 131.061822 三次函数的最小二乘曲线如下图所示：计算方法上机报告19c、四次函数：一天的平均气温为： 21.20004 次函数的系数: 16.7939 -3.7050 0.8909 -0.053
9、2 0.0009使用 4 次函数拟合的误差是：59.04118四次函数的最小二乘曲线如下图所示：d、指数函数：一天的平均气温为： 21.2000指数函数的系数是： a=26.160286，b=0.004442，c=14.081900使用指数函数拟合的误差是： 57.034644指数函数的最小二乘曲线如下图所示：计算方法上机报告20通过上述几种拟合可以发现，多项式的次数越高，计算拟合的效果越好，误差越小，说明结果越准确；同时，指数多项式拟合的次数虽然不高，但误差最小，说明结果最准确。4.设计算法，求出非线性方程 52640x的所有根，并使误差不超过410。（1）算法思想首先，研究函数的形态，确定
10、根的范围；通过剖分区间的方法确定根的位置，然后利用二分法的基本原理进行求解，找到满足精度要求的解。二分法是产生一串区间，使新区间 是旧区间 的一个子区间，其长度是)1(kI)(kI的一半，且有一个端点是 的一个端点。由区间 确定区间 的)(kI )(kI ,)1()(kx)1(kI方法是计算区间 的中点)(kI )(211()( kkxx计算方法上机报告21若 ，则取 ，否则取 ，重复这一0)(2)kxf ,)2()1(kkxI ,)1()2()1(kkxI过程即可。显然，每次迭代使区间长度减小一半，故二分法总是收敛的。（2）算法结构1. ；0)(fxf1)(fx2. If then stop
11、103. If then 输出 作为根； stop2|f)0(x4. If then 输出 作为根； stop1| )1(5. xx)(2()06. If then 输出 作为根； stop| )1()1(x7. fxf8. If then 输出 作为根； 2|x9. If then 01f9.1 ；)(x0felse9.2 ；)1(1f10. go to 5 （3）Matlab 源程序x=-100:100;y=6*(x.5)-45*(x.2)+20; %非线性方程组的表达式g=; for i=-100:1:100 %确定根所在的区间k=i+1;if (y(x=i).*y(x=k)=10(-4) if subs(f,x,x0)*subs(f,x,(x0+x1)/2)epsx0=(x0+x1)/2;elsex1=(x0+x1)/2;endendroot=x0 %输出方程的根
