matlab上机实验报告

1上海电力学院数值计算方法上机实习报告院系：数理学院专业年级：可再生能源科学与工程2013级学生姓名：司晓东学号：ys1310701006指导教师：黄建雄2012年12月26日2数值计算方法上机实习题1．设105dxxxInn，（1）由递推公式nIInn151，从0I的几个近似值出发，计算20I；解：I0=1015xdx=0.1823计算I20编辑matlab命令如下：I=0.1823forn=1:1:20,I=-5*I+1/n;fprintf('%.1d%.4f\n',n,I);end结果：（2）粗糙估计20I，用nIInn51511，计算0I；解：I20=20105xxdx使用复合中点公式进行积分，相应的matlab程序如下：3I=0;forh=0:0.001:1,m=h+0.0005;I=I+0.001*m^20/(5+m);fprintf('%.1d%.4f\n',m,I);enddisp(I);fork=1:20,n=21-k;I=0.2*(1/n-I);fprintf('%.1d%.4f\n',n,I);enddisp(I)结果：程序结束时输出两个I值，第一个表示I20，第二个表示I0；分别为I20=0.0082I0=0.1823（3）分析结果的可靠性及产生此现象的原因(重点分析原因)。从上述计算中分析得到如果先得到I0，再从I0由递推公式得到I20，I20结果跟精确值相比误差很大；如果先估算I20，在从I20有递推公式得到I0，I0的结果跟精确值相比近似相等。4原因分析：如果从I0推I20的近似值，需要用到递推公式In=-5In-1+1/n，I0本身结果是有误差的；经过递推公式计算20次，就等于误差被认为的放大5的20次方倍，所以得到的I20与其精确值相差甚远。如果从I20推I0的近似值，需要用到In-1=0.2（1/n-In），尽管I20本身有误差，但是经过20次运算，其误差缩小到原来的0.2的20次方倍，所以得到的I0与其精确值比较相近。2．求方程0210xex的近似根，要求41105kkxx，并比较计算量。（1）在[0，1]上用二分法；Matlab程序如下：a=0;b=1;c=b-a;n=0whilec0.0005,x=(a+b)/2;f=exp(x)+10*x-2;iff0,b=x;c=b-a;elseiff0,a=x;c=b-a;elsex=x;c=0;endn=n+1;fprintf('%.1d%.4f%.4f\n',n,x,c);end结果如下：解得到;x=0.09035（2）取初值00x，并用迭代1021xkex；采用matlab进行迭代的程序如下：x=0;c=1;n=0;whilec0.0005,m=x;m=(2-exp(m))/10;c=abs(m-x);x=m;n=n+1;fprintf('%.1d%.4f%.4f\n',n,x,c);end结果：解得x=0.0905（3）加速迭代的结果；采用matlab进行迭代的程序如下：x=0;n=0;a=0;b=1;whileabs(a-b)0.0005,n=n+1;a=x;y=(2-exp(x))/10;z=(2-exp(y))/10;x=x-(y-x)^2/(z-2*y+x);b=x;fprintf('%.1d%.4f%.4f\n',n,x,abs(a-b));end结果如下：（4）取初值00x，并用牛顿迭代法；Matlab程序如下：6x=0;a=1;n=0;whileabs(a)0.0005,n=n+1;a=(exp(x)+10*x-2)/(exp(x)+10);x=x-a;fprintf('%.1d%.4f%.4f\n',n,x,abs(a));end运行结果：（5）分析绝对误差。迭代次数二分法代数式迭代加速迭代牛顿迭代X（k）ErroeX（k）ErroeX（k）ErroeX（k）Erroe10.50000.50000.10000.10000.09050.09050.09090.090920.25000.25000.08950.01050.09050．00000.09050.000430.12500.12500.09060.001240.06250.06250.09050.000150.09380.031360.07810.015670.08590.007880.08980.003990.09180.0020100.09080.0010110.09030.0005我们可以看到，在运算要求到同一精度的情况下，采用（1）的二分法运算了11次，采用（2）的方法运算了4次，采用（3）的加速迭代法运算了2次，采用（4）的牛顿迭代法也需运算2次。也就是说牛顿的迭代的收敛速度与加速迭代速度都是超线性收敛的，而简单迭代法是线性收敛的。而二分法收敛速度较慢,所以在工程上不经常使用。3．钢水包使用次数多以后，钢包的容积增大，数据如下：x23456789y6.428.29.589.59.7109.939.991011121314151610.4910.5910.6010.810.610.910.76试从中找出使用次数和容积之间的关系，计算均方差。（注：增速减少，用何种模型）解：将使用次数x与体积y的关系用matlab采用如下程序绘制在二维坐标系:x=[2345678910111213141516];y=[6.428.29.589.59.7109.939.9910.4910.5910.6010.810.610.910.76];plot(x,y,'b*-');7结果如下：由数据点分布图可知，拟合曲线y=f(x)随着x的增加而上升，但上升速度由快到慢，当x趋于无穷大时，y趋于某个常数，故曲线有一水平渐进线。根据上述特征很容易想到用Logistic模型来拟合该曲线。设y=f(x)的形式为y=aeb/x（a0,b0），两边取对数，得lny=lna+b/x。记A=lna，B=b，并引入新变量m=lny，n=1/x。引入新变量后的数据表如下：x2345678910111213141516n=1/x0.50000.33330.25000.20000.16670.14290.12500.11110.10000.09090.08330.07690.07140.06670.0625m=lny1.85942.10412.25972.25132.27212.30262.29562.30162.35042.35992.36092.37952.36092.38882.3758Matlab拟合程序如下：n=[0.50000.33330.25000.20000.16670.14290.12500.11110.10000.09090.08330.07690.07140.06670.0625];m=[1.85942.10412.25972.25132.27212.30262.29562.30162.35042.35992.36092.37952.36092.38882.3758];polyfit(n,m,1)运行的结果：由此可得A=2.4578B=-1.1107a=eA=11.6791b=B=-1.1107由此得到使用次数与容积的函数为1.110711.6791xye将统计表和函数用matlab绘制在同一坐标系内程序如下：x1=[2345678910111213141516];y1=[6.428.29.589.59.7109.939.9910.4910.5910.6010.810.610.910.76];x=2:1:16;y=11.679*exp(-1.1107*x.^(-1));holdon;plot(x,y,'ro-');plot(x1,y1,'b*-');8结果如图：计算均方差s，matlab程序如下：y=[6.428.29.589.59.7109.939.9910.4910.5910.6010.810.610.910.76];s=0;forn=2:1:16,a=abs(11.679*exp(-1.1107*n.^(-1))-y(n-1));s=s+a.^2;ends=(s/15).^(1/2);disp(s);运算结果均方差S=0.2438小结：根据已给的条件计算函数是十分困难的，但通过对离散点的分析及变化规律找出其中的规律，并通过计算来得到实际的函数是十分有用的方法。本题就是这样做的一个典型，在n=1/x和m=lny的基础上找到了它们之间的关系并通过这种关系来拟合原函数，并最终验证计算结果。94．设410100141010014101101410010141001014A，625250b，bxA分析下列迭代法的收敛性，并求42110kkxx的近似解及相应的迭代次数。（1）JACOBI迭代；解matlab计算程序如下：A=[4-10-100;-14-10-10;0-14-10-1;-10-14-10;0-10-14-1;00-10-14];b=[0;5;-2;5;-2;6];error=1;D=diag(diag(A));L=D-tril(A);U=D-triu(A);X=zeros(size(b));whileerror0.0001,X=D\(b+L*X+U*X);error=norm(b-A*X)/norm(b);enddisp(x);disp(error);解得X=[0.9999;1.9999;0.9998;1.9999;0.9998;1.9999]error=7.0206e-05（2）GAUSS-SEIDEL迭代；A=[4-10-100;-14-10-10;0-14-10-1;-10-14-10;0-10-14-1;00-10-14];b=[0;5;-2;5;-2;6];error=1;D=diag(diag(A));L=D-tril(A);U=D-triu(A);X=zeros(size(b));whileerror0.0001,X=(D-L)\(b+U*X);error=norm(b-A*X)/norm(b);enddisp(x);disp(error);解得X=[0.9998;1.9998;0.9998;1.9999;0.9999;1.9999]error=5.5892e-0510（3）SOR迭代（95.0,95.1,334.1）。N=1.334使用matlab求解程序如下：A=[4-10-100;-14-10-10;0-14-10-1;-10-14-10;0-10-14-1;00-10-14];b=[0;5;-2;5;-2;6];error=1;D=diag(diag(A));L=D-tril(A);U=D-triu(A);X=zeros(size(b));whileerror0.0001,n=1.334;X=(D-n*L)\[(1-n)*D+n*U]*X+n*[(D-n*L)\b];error=norm(b-A*X)/norm(b);disp(X);enddisp(error);此循环得到的X=[0.9999;2.0000;1.0000;1.9999;1.0000;2.0000]error=6.3632e-05N=1.95使用matlab求解程序如下：A=[4-10-100;-14-10-10;0-14-10-1;-10-14-10;0-10-14-1;00-10-14];b=[0;5;-2;5;-2;6];error=1;D=diag(diag(A));L=D-tril(A);U=D-triu(A);X=zeros(size(b));whileerror0.0001,n=1.95;X=(D-n*L)\[(1-n)*D+n*U]*X+n*[(D-n*L)\b];error=norm(b-A*X)/norm(b);disp(X);enddisp