2017数值计算方法上机实习报告学院：专业：班级：姓名：学号：数值计算方法上机实习题1． 设⎰+=105dx x x I nn ， （1） 由递推公式nI I n n 151+-=-，从0=0.1823I ,1824.00=I 出发，计算20I ； （2） 20=0I ，20=10000I , 用n I I n n 515111+-=--，计算0I ； （3） 分析结果的可靠性及产生此现象的原因(重点分析原因)。

解：（1）分别令I 0的近似值为0.1823、0.1824，MATLAB 程序如下：I=0.1823; %题中的已知数据for n=1:20;I=(-5)*I+1/n; %由递推公式所得endfprintf('I20=%f\n',I)M=0.1824; %与I 的计算结果形成对比for i=1:20;M=(-5)*M+1/i; %由递推公式所得endfprintf('M20=%f\n',M)%% 输出结果I20=-2055816073.851284M20=7480927090.212283（2）分别令I 20的近似值为0、10000，MATLAB 程序如下：I=0; %赋予I20的初始值for n=0:19;I=(-1/5)*I+1/(5*(20-n)); %由递推公式所得end fprintf('I0=%f\n',I)M=10000;for i=0:19;M=(-1/5)*M+1/(5*(20-i));%由递推公式所得endfprintf('M0=%f\n',M)%% 输出结果I0=0.182322M0=0.182322（3）分析：由输出结果可看出第一种算法为不稳定算法，第二种算法为稳定算法。

由于误差*000I I e -=02211*1*11*555)(5)15(15e e e I I nI n I I I e n n n n n n n n n n ===-=+--+-=-=------第一种算法为正向迭代算法，每计算一步误差增长5倍，虽然所给的初始值很接近，随着n 的增大，误差也越来越大。

*nn n I I e -= n n n n n n n n n e I I n I n I e I I e 51)(51)5151(5151***111=-=+--+-==-=--- n n e e )51(0= 第二种算法为倒向迭代算法，每计算一步误差缩小5倍，虽然所给的初始值之间差很多，随着n 的增大，误差也越来越小。

算法趋近稳定，收敛，可以选择这种算法。

2． 求方程0210=-+x e x 的近似根，要求41105-+⨯<-k k x x ，并比较计算量。

（1） 在[0，1]上用二分法；（2） 取初值00=x ，并用迭代1021x k e x -=+； （3） 加速迭代的结果；（4） 取初值00=x ，并用牛顿迭代法；（5） 分析绝对误差。

解：（1）利用二分法，MATLAB 程序如下：%% 二分法程序clear allclca=0;b=1;f=@(x)(exp(x)+10*x-2);%@是定义函数句柄的运算符c=(a+b)/2;%取区间中点i=0;%分割次数while abs(f(c))>5*10^(-4)%判断f(x)的精度是否满足要求if f(a)*f(c)<0b=c;c=(a+b)/2;elseif f(b)*f(c)<0a=c;c=(b+a)/2;endi=i+1;endfprintf('二分法运算次数为%i\n',i)fprintf('二分法计算结果为%f\n',c)%% 输出结果 二分法运算次数为13二分法计算结果为0.090515（2）用题目中给出的迭代法，取初始值x(1)=0，并用迭代x(i)=(2-exp(x(i-1)))/10，MATLAB 程序如下：%% 不动点迭代clear allclcx0=0;x=x0;for k=1:10000 %规定迭代次数上限y=(2-exp(x))/10; %迭代结果存到y中if abs(x-y)<5*10^(-4)fprintf('初始值为x0%i\n迭代次数为%i\n',x0,k);breakendx=y;if k==10000;fprintf('迭代次数超出上限%i\n',k);endendfprintf('迭代法计算结果为%f\n',y);%% 输出结果初始值x0为0迭代次数为4迭代法计算结果为0.090513（3）利用加速迭代法，MATLAB程序如下：%% 加速迭代法x=0;f=@(x)(( 2-exp(x))./10);y=f(x); y0=f(y);y1=x-((y-x)^2)/(y0-2*y+x);i=1;while abs(y1-x)> 5*10^(-4)x=y1;y=f(x);y0=f(y);y1=x-((y-x)^2)/(y0-2*y+x);i=i+1;end%% 输出结果y1= 0.090525i=2（4）牛顿迭代法：取初始值x=0，MA TLAB程序如下：%% 清空环境变量clear allclc%% 牛顿迭代法x=0;y=x-(exp(x)+10*x-2)./(exp(x)+10);i=1;while abs(x-y)>0.0005x=y;y=x-(exp(x)+10*x-2)./(exp(x)+10);i=i+1;end%% 输出结果y= 0.090525i=2（5）分析绝对误差根据方程解求得对应的绝对误差如下表所示：通过上述表格可以看出，二分法运算了11次，不动点迭代方法运算了4次， Atiken 加速迭代法运算了2次，牛顿迭代法运算2次.比较绝对误差可以发现 Atiken 加速迭代和牛顿迭代的计算结果的绝对误差较小。

下面就其原因进行分析： 我们知道1212()()g x g x L x x -≤-,其中L 为压缩常数，并且01L ≤<。

进行误差估计：10*1kk L x x x x L-≤--，当L 较小时，收敛较快，反之，当L 很靠近1时，收敛很慢。

若1L ≥时，则迭代不收敛。

由()1g x L '≤< ，()1()lim =(0,1)!p k k ke g x C C p e p *+→∞=≠≥常数，收敛阶为p 。

二分法和不动点迭代为线性收敛；Atiken 迭代和牛顿迭代是平方收敛，是超线性收敛的。

其中，二分法L=0.5；用迭代1021xk e x -=+，2()==0.27101010x x e e e g x '⎛⎫-'≤≈ ⎪⎝⎭，比0.5小，因此收敛比二分法快。

不动点迭代为线性收敛，Atiken 迭代速度与牛顿迭代速度都是超线性收敛的，因而收敛速度较快。

相比之下，只有二分法的收敛速度较慢,绝对误差最大。

3．钢水包使用次数多以后，钢包的容积增大，数据如下：试从中找出使用次数和容积之间的关系，计算均方差。

（用ax b y c x+=+拟合） 解：拟合曲线的模型是ax b y c x+=+，将原模型变为)()(b ax x c y +=+，采用非线性最小二程 法。

按照最小二乘原理，应选取参数a ，b ，c 使得表达式∑=+-+=151k 2))()((Z b ax x c y k k k 达到极小值。

具体的方法就是对Z 关于a 、b 、c 求偏导数，并置这些偏导数等于0，相对应的方程组如下所示：⎪⎪⎪⎩⎪⎪⎪⎨⎧=+-+-=∂∂=+-+-=∂∂=+-+-=∂∂∑∑∑---0))()((20))()((20))()((2151151151k k k k k k k k k k k k k k b ax y x c y c Z b ax y x c b Z b ax y x c x a Z 通过对上述的偏导数方程组进行整理，可以将之写成的b AX =的形式，利用 MA TLA 充当计算器的作用，可以进行X 的求解，对应的就是a 、b 、c 三个参数的值。

下面的程序就是进行求解方程组、画拟合曲线以及求取均方差的程序，具体如下：%求拟合方程，画拟合曲线x=[2 3 4 5 6 7 8 9 10 11 12 13 14 15 16];y=[6.42 8.2 9.58 9.5 9.7 10 9.93 9.99 10.49 10.59 10.60 10.8 10.6 10.9 10.76];A=[sum(x.^2),sum(x),-sum(x.*y); sum(x),15,-sum(y); sum(x.*y),sum(y),-sum(y.^2)]; B=[sum(x.^2.*y);sum(x.*y);sum(y.^2.*x)];f=A\B;a=f(1);b=f(2);c=f(3);z=linspace(2,50,48);Y=(f(1)*z+f(2))./(f(3)+z);plot(x,y,'r*',z,Y ,'b-');fprintf('a=%f\nb=%f\nc=%f\n',a,b,c);fprintf('拟合出来的方程式为:\n(%7.4f+x)y=%7.4f+%7.4fx\n',f(3),f(2),f(1));%求均方差for i=1:15y1(i)=(f(2)+f(1)*x(i))/(f(3)+x(i));endc=0;for i=1:15c=c+(y(i)-y1(i))^2;endjfc=sqrt(c/15);fprintf('均方差为%f\n',jfc)%%运行结果：a=11.340048b=-12.495325c=-0.340291拟合出来的方程式为:(-0.3403+x)y=-12.4953+11.3400x均方差为0.220812对应的拟合曲线如下图：分析总结：此方法利用线性方法求解非线性问题，避免了非线性误差较大的问题。

4．设⎪⎪⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛----------------=410100141010014101101410010141001014A ，⎪⎪⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛--=625250b ，b x =A 分析下列迭代法的收敛性，并求42110-+≤-k k x x 的近似解及相应的迭代次数。

（1） JACOBI 迭代；（2） GAUSS-SEIDEL 迭代；（3） SOR 迭代（取0.1:0.1:1.9ω=，找到迭代步数最少的*ω）。

解：（1）JACOBI 迭代%% JACOBI函数function tx=jacobi(A,b,imax,x0,tol)%初始值x0，次数imax，精度toldel=10^(-10);tx=[x0];n=length(x0);for i=1:ndg=A(i,i);if abs(dg)<deldisp('对角元太小');%防止出现溢出现象returnendendfor k=1:imax %jacobi循环for i=1:nsm=b(i);for j=1:nif j~=ism=sm-A(i,j)*x0(j);endendx(i)=sm/A(i,i);%x(1)到x(n)即完成一次迭代endtx=[tx;x];%矩阵中又加一行if norm(x-x0)<tolreturnelsex0=x;endend%% 主函数程序clear allclcA=[4 -1 0 -1 0 0;-1 4 -1 0 -1 0;0 -1 4 -1 0 -1;-1 0 -1 4 -1 0;0 -1 0 -1 4 -1;0 0 -1 0 -1 4]; b=[0 5 -2 5 -2 6];x0=[0 0 0 0 0 0];imax=100;%迭代次数tol=10^(-4);%精度tx=jacobi(A,b,imax,x0,tol);for j=1:length(tx)fprintf('%d%f%f%f%f%f%f\n',j,tx(j,1),tx(j,2),tx(j,3),tx(j,4),tx(j,5),tx(j,6));end%% 运行结果10.0000000.0000000.0000000.0000000.0000000.000000 20.0000001.250000-0.5000001.250000-0.5000001.500000 30.6250001.0000000.5000001.0000000.5000001.250000 40.5000001.6562500.3125001.6562500.3125001.750000 50.8281251.5312500.7656251.5312500.7656251.656250 60.7656251.8398440.6796881.8398440.6796881.882813 70.9199221.7812500.8906251.7812500.8906251.839844 80.8906251.9252930.8505861.9252930.8505861.945313 90.9626461.8979490.9489751.8979490.9489751.925293 100.9489751.9651490.9302981.9651490.9302981.974487 110.9825741.9523930.9761961.9523930.9761961.965149 120.9761961.9837420.9674841.9837420.9674841.988098 130.9918711.9777910.9888951.9777910.9888951.983742 140.9888951.9924150.9848311.9924150.9848311.994448 150.9962081.9896390.9948201.9896390.9948201.992415 160.9948201.9964620.9929231.9964620.9929231.997410 170.9982311.9951670.9975831.9951670.9975831.996462 180.9975831.9983490.9966991.9983490.9966991.998792 190.9991751.9977450.9988731.9977450.9988731.998349 200.9988731.9992300.9984601.9992300.9984601.999436 210.9996151.9989480.9994741.9989480.9994741.999230 220.9994741.9996410.9992821.9996410.9992821.999737 230.9998201.9995090.9997551.9995090.9997551.999641 240.9997551.9998320.9996651.9998320.9996651.999877 250.9999161.9997710.9998861.9997710.9998861.999832 260.9998861.9999220.9998441.9999220.9998441.999943 270.9999611.9998930.9999471.9998930.9999471.999922 280.9999471.9999640.9999271.9999640.9999271.999973 290.9999821.9999500.9999751.9999500.9999751.999964 （2）GAUSS-SEIDEL迭代%% GAUSS-SEIDEL函数function tx=gseidel(A,b,imax,x0,tol)%初始值x0，次数imax，精度toldel=10^-10; tx=[x0];n=length(x0);%tx是个二维矩阵，存储着每一步迭代的结果for i=1:ndg=A(i,i);if abs(dg)<deldisp('对角元太小');returnendendfor k=1:imax %G-S循环x=x0;for i=1:nsm=b(i);for j=1:nif j~=ism=sm-A(i,j)*x(j);endendx(i)=sm/A(i,i);endtx=[tx;x];if norm(x-x0)<tolreturnelsex0=x;endend%% 主函数程序clear allclcA=[4 -1 0 -1 0 0;-1 4 -1 0 -1 0;0 -1 4 -1 0 -1;-1 0 -1 4 -1 0;0 -1 0 -1 4 -1;0 0 -1 0 -1 4];b=[0 5 -2 5 -2 6];x0=[0 0 0 0 0 0];imax=100;%迭代次数tol=10^(-4);%精度tx=gseidel(A,b,imax,x0,tol);for j=1:length(tx)fprintf('%d %f %f %f %f %f %f\n',j,tx(j,1),tx(j,2),tx(j,3),tx(j,4),tx(j,5),tx(j,6)) end%% 运行结果1 0.000000 0.000000 0.000000 0.000000 0.000000 0.0000002 0.000000 1.250000 -0.187500 1.203125 0.113281 1.4814453 0.613281 1.384766 0.517334 1.560974 0.606796 1.7810334 0.736435 1.715141 0.764287 1.776880 0.818263 1.8956385 0.873005 1.863889 0.884102 1.893843 0.913342 1.9493616 0.939433 1.934219 0.944356 1.949283 0.958216 1.9756437 0.970875 1.968362 0.973322 1.975603 0.979902 1.9883068 0.985991 1.984804 0.987178 1.988268 0.990344 1.9943819 0.993268 1.992698 0.993837 1.994362 0.995360 1.99729910 0.996765 1.996490 0.997038 1.997291 0.997770 1.99870211 0.998445 1.998313 0.998577 1.998698 0.998928 1.99937612 0.999253 1.999189 0.999316 1.999374 0.999485 1.99970013 0.999641 1.999610 0.999671 1.999699 0.999752 1.99985614 0.999827 1.999813 0.999842 1.999855 0.999881 1.99993115 0.999917 1.999910 0.999924 1.999931 0.999943 1.99996716 0.999960 1.999957 0.999964 1.999967 0.999973 1.999984（3）SOR迭代%% SOR函数function tx=sor(A,b,imax,x0,tol,w)%初始值x0，次数imax，精度toldel=10^-10;tx=[x0];n=length(x0);for i=1:ndg=A(i,i);if abs(dg)<deldisp('对角元太小');returnendendfor k=1:imax %SOR循环x=x0;for i=1:nsm=b(i);for j=1:nif j~=ism=sm-A(i,j)*x(j);%先高斯迭代endendx(i)=sm/A(i,i);x(i)=w*x(i)+(1-w)*x0(i);%比较高斯与超松弛迭代公式，补上省缺的项endtx=[tx;x];if norm(x-x0)<tolreturnelsex0=x;endend%% 主函数程序clear allclcA=[4 -1 0 -1 0 0;-1 4 -1 0 -1 0;0 -1 4 -1 0 -1;-1 0 -1 4 -1 0;0 -1 0 -1 4 -1;0 0 -1 0 -1 4]; b=[0 5 -2 5 -2 6];x0=[0 0 0 0 0 0];imax=100;迭代次数tol=10^(-4);%精度w=1.2;%给定不同的w,w=0.1:0.1:1.9，找出使SOR 迭代法收敛速度最快tx=sor(A,b,imax,x0,tol,w);for j=1:size(tx,1)fprintf('%d %f%f%f%f%f%f\n',j,tx(j,1),tx(j,2),tx(j,3),tx(j,4),tx(j,5),tx(j,6))end%% 运行结果1 0.0000000.0000000.0000000.0000000.0000000.0000002 0.0000001.500000-0.1500001.4550000.2865001.8409503 0.8865001.5069000.8708551.8221560.8937021.9611774 0.8214171.9744120.9531531.9360500.9827511.9885365 1.0088551.9885450.9833092.0052650.9981531.9967326 0.9963721.9956411.0026291.9980940.9975092.0006957 0.9988462.0005670.9992811.9990721.0005991.9998258 1.0001231.9998870.9997792.0003360.9998951.9999379 1.0000421.9999371.0001071.9999460.9999672.00003510 0.9999572.0000220.9999791.9999821.0000181.99999211 1.0000101.9999980.9999962.0000110.9999971.999999W 值的范围在0~2之间，SOR 迭代才会收敛，令w=0.1:0.1:1.9，找出w*，使得SOR 迭代的收敛速度最快。