最新上海电力学院数值计算方法上机实习题

2017数值计算方法上机实习报告学院：专业：班级：姓名：学号：数值计算方法上机实习题1．设105dxxxInn，（1）由递推公式nIInn151，从0=0.1823I,1824.00I出发，计算20I；（2）20=0I，20=10000I,用nIInn515111，计算0I；（3）分析结果的可靠性及产生此现象的原因(重点分析原因)。解：（1）分别令I0的近似值为0.1823、0.1824，MATLAB程序如下：I=0.1823;%题中的已知数据forn=1:20;I=(-5)*I+1/n;%由递推公式所得endfprintf('I20=%f\n',I)M=0.1824;%与I的计算结果形成对比fori=1:20;M=(-5)*M+1/i;%由递推公式所得endfprintf('M20=%f\n',M)%%输出结果I20=-2055816073.851284M20=7480927090.212283（2）分别令I20的近似值为0、10000，MATLAB程序如下：I=0;%赋予I20的初始值forn=0:19;I=(-1/5)*I+1/(5*(20-n));%由递推公式所得endfprintf('I0=%f\n',I)M=10000;fori=0:19;M=(-1/5)*M+1/(5*(20-i));%由递推公式所得endfprintf('M0=%f\n',M)%%输出结果I0=0.182322M0=0.182322（3）分析：由输出结果可看出第一种算法为不稳定算法，第二种算法为稳定算法。由于误差*000IIe02211*1*11*555)(5)15(15eeeIInInIIIennnnnnnnnn第一种算法为正向迭代算法，每计算一步误差增长5倍，虽然所给的初始值很接近，随着n的增大，误差也越来越大。*nnnIIennnnnnnnneIInInIeIIe51)(51)5151(5151***111nnee)51(0第二种算法为倒向迭代算法，每计算一步误差缩小5倍，虽然所给的初始值之间差很多，随着n的增大，误差也越来越小。算法趋近稳定，收敛，可以选择这种算法。2．求方程0210xex的近似根，要求41105kkxx，并比较计算量。（1）在[0，1]上用二分法；（2）取初值00x，并用迭代1021xkex；（3）加速迭代的结果；（4）取初值00x，并用牛顿迭代法；（5）分析绝对误差。解：（1）利用二分法，MATLAB程序如下：%%二分法程序clearallclca=0;b=1;f=@(x)(exp(x)+10*x-2);%@是定义函数句柄的运算符c=(a+b)/2;%取区间中点i=0;%分割次数whileabs(f(c))5*10^(-4)%判断f(x)的精度是否满足要求iff(a)*f(c)0b=c;c=(a+b)/2;elseiff(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程序如下：%%不动点迭代clearallclcx0=0;x=x0;fork=1:10000%规定迭代次数上限y=(2-exp(x))/10;%迭代结果存到y中ifabs(x-y)5*10^(-4)fprintf('初始值为x0%i\n迭代次数为%i\n',x0,k);breakendx=y;ifk==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;whileabs(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，MATLAB程序如下：%%清空环境变量clearallclc%%牛顿迭代法x=0;y=x-(exp(x)+10*x-2)./(exp(x)+10);i=1;whileabs(x-y)0.0005x=y;y=x-(exp(x)+10*x-2)./(exp(x)+10);i=i+1;end%%输出结果y=0.090525i=2（5）分析绝对误差根据方程解求得对应的绝对误差如下表所示：迭代方式迭代次数迭代结果绝对误差二分法迭代110.09057625.1099e-05不动点迭代40.0905126-1.2501e-05加速迭代20.0905251-1.3073e-09牛顿迭代20.0905251-1.3073e-09通过上述表格可以看出，二分法运算了11次，不动点迭代方法运算了4次，Atiken加速迭代法运算了2次，牛顿迭代法运算2次.比较绝对误差可以发现Atiken加速迭代和牛顿迭代的计算结果的绝对误差较小。下面就其原因进行分析：我们知道1212()()gxgxLxx,其中L为压缩常数，并且01L。进行误差估计：10*1kkLxxxxL，当L较小时，收敛较快，反之，当L很靠近1时，收敛很慢。若1L时，则迭代不收敛。由()1gxL，()1()lim=(0,1)!pkkkegxCCpep常数，收敛阶为p。二分法和不动点迭代为线性收敛；Atiken迭代和牛顿迭代是平方收敛，是超线性收敛的。其中，二分法L=0.5；用迭代1021xkex，2()==0.27101010xxeeegx，比0.5小，因此收敛比二分法快。不动点迭代为线性收敛，Atiken迭代速度与牛顿迭代速度都是超线性收敛的，因而收敛速度较快。相比之下，只有二分法的收敛速度较慢,绝对误差最大。3．钢水包使用次数多以后，钢包的容积增大，数据如下：x23456789y6.428.29.589.59.7109.939.991011121314151610.4910.5910.6010.810.610.910.76试从中找出使用次数和容积之间的关系，计算均方差。（用axbycx拟合）解：拟合曲线的模型是axbycx，将原模型变为)()(baxxcy，采用非线性最小二程法。按照最小二乘原理，应选取参数a，b，c使得表达式151k2))()((Zbaxxcykkk达到极小值。具体的方法就是对Z关于a、b、c求偏导数，并置这些偏导数等于0，相对应的方程组如下所示：0))()((20))()((20))()((2151151151kkkkkkkkkkkkkkbaxyxcycZbaxyxcbZbaxyxcxaZ通过对上述的偏导数方程组进行整理，可以将之写成的bAX的形式，利用MATLA充当计算器的作用，可以进行X的求解，对应的就是a、b、c三个参数的值。下面的程序就是进行求解方程组、画拟合曲线以及求取均方差的程序，具体如下：%求拟合方程，画拟合曲线x=[2345678910111213141516];y=[6.428.29.589.59.7109.939.9910.4910.5910.6010.810.610.910.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));%求均方差fori=1:15y1(i)=(f(2)+f(1)*x(i))/(f(3)+x(i));endc=0;fori=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，bxA分析下列迭代法的收敛性，并求42110kkxx的近似解及相应的迭代次数。（1）JACOBI迭代；（2）GAUSS-SEIDEL迭代；（3）SOR迭代（取0.1:0.1:1.9，找到迭代步数最少的*）。解：（1）JACOBI迭代%%JACOBI函数functiontx=jacobi(A,b,imax,x0,tol)%初始值x0，次数imax，精度toldel=10^(-10);tx=[x0];n=length(x0);fori=1:ndg=A(i,i);ifabs(dg)deldisp('对角元太小');%防止出现溢出现象returnendendfork=1:imax%jacobi循环fori=1:nsm=b(i);forj=1:nifj~=ism=sm-A(i,j)*x0(j);endendx(i)=sm/A(i,i);%x(1)到x(n)即完成一次迭代endtx=[tx;x];%矩阵中又加一行ifnorm(x-x0)tolreturnelsex0=x;endend%%主函数程序clearallclcA=[4-10-100;-14-10-10;0-14-10-1;-10-14-10;0-10-14-1;00-10-14];b=[05-25-26];x0=[000000];imax=100;%迭代次数tol=10^(-4);%精度tx=jacobi(A,b,imax,x0,tol);forj=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.00000020.0000001.250000-0.5000001.250000-0.5000001.50000030.6250001.0000000.5000001.0000000.5000001.25000040.5000001.6562500.3125001.6562500.3125001.75000050.8281251.5312500.7656251.5312500.7656251.65625060.7656251.8398440.6796881.8398440.6796881.88281370.91992