1、数值计算方法上机实习题NEW数值计算方法上机实习题1 设，（1） 由递推公式 ，从, 出发，计算；程序如下：function I=myhs(I0,n)if n=1 I=myhs(I0,n-1)*(-5)+1/n;elseif n=0 I=I0;end命令行窗口输入：I0=0.1822;I1=myhs(I0,20);I0=0.1823;I2=myhs(I0,20);运行结果：当I0=0.1822时，I20=-1.1593e+10。当I0=0.1823时，I20= -2.0558e+9。（2） ，, 用，计算；程序如下：function I=myhs2(I20,n)if(n5*1e-4) c=(a
2、+b)/2; y=exp(c)+10*c-2; if(y0) b=c; else a=c; end n=n+1;endc得到结果：c= 0.090332，n=11（2） 取初值，并用迭代；输入程序如下：xk=0;xki=(2-exp(xk)/10;n=1;while(abs(xki-xk)5*1e-4) xk=xki; xki=(2-exp(xk)/10; n=n+1;endxki程序运行结果为：xki= 0.090513，n=4（3） 加速迭代的结果；输入程序如下xk=0;yk=(2-exp(xk)/10;zk=(2-exp(yk)/10;xki=xk-(yk-xk)2/(zk-2*yk+x
3、k);n=1;while(abs(xki-xk)5*1e-4) xk=xki; yk=(2-exp(xk)/10; zk=(2-exp(yk)/10; xki=xk-(yk-xk)2/(zk-2*yk+xk); n=n+1;endxki程序运行的结果为xki= 0.090525，n=2（4） 取初值，并用牛顿迭代法；输入程序如下：xk=0;xki=xk-(exp(xk)+10*xk-2)/(exp(xk)+10);n=1;while(abs(xki-xk)5*1e-4) xk=xki; xki=xk-(exp(xk)+10*xk-2)/(exp(xk)+10); n=n+1;endxki程序运
4、行的结果为xki= 0.090525，n=2（5） 分析绝对误差。通过程序x=solve(exp(x)+10*x-2=0);解得x= 0.090525用二分法求得x的值为0.090332，迭代次数为11，绝对误差为0.000193。用不动点迭代法求得x的值为0.090513，迭代次数为4，绝对误差为0.000012。用加速迭代法和牛顿迭代法求得x的值均为0.090525，迭代次数均为2，绝对误差为0。3钢水包使用次数多以后，钢包的容积增大，数据如下：x23456789y6.428.29.589.59.7109.939.991011121314151610.4910.5910.6010.810.
5、610.910.76试从中找出使用次数和容积之间的关系，计算均方差。（用拟合）解：拟合曲线的模型为y=(ax+b)/(c+x)，将原模型变为(c+x)y=(ax+b)，采用非线性最小二乘法。应选取参数a、b、c使得达到极小值。具体方法：对Z分别求a、b、c的偏导，并使偏导为0，相应的方程组如下对上述方程组进行整理，写成AX=B的形式，利用MATLAB求解X，对应的就是a、b、c三个参数的值。具体程序如下：%方程组求解和画拟合曲线C=xlsread(test3.xlsx);x=C(1,:);y=C(2,:);A=sum(x.2),sum(x),-sum(x.*y);sum(x),15,-sum(
6、y);sum(x.*y),sum(y),-sum(y.2);B=sum(x.2.*y);sum(x.*y);sum(x.*y.2);Z=AB;a=Z(1);b=Z(2);c=Z(3);X=linspace(1,50,50);Y=zeros(1,50);for i=1:50 Y(i)=(a*X(i)+b)/(c+X(i);endplot(x,y,*,X,Y,-);%计算方差z1=0;for i=2:16 z1=z1+Y(i);endz1=z1/15;sum=0;for i=2:16 sum=sum+(Y(i)-z1)2;endfc=sqrt(sum/15);运行结果为：a = 11.3400，b
7、 = -12.4953，c = -0.3403拟合方程为y=(11.34x-12.4953)/(-0.3403+x),均方差为1.22854设，分析下列迭代法的收敛性，并求的近似解及相应的迭代次数。（1） JACOBI迭代；解：输入程序如下A=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;D=diag(diag(A);L=D-tril(A);U=D-triu(A);Bj=D(L+U);V,landa=eig(Bj);Maxlanda=
8、abs(landa(1,1);for i=2:size(landa) if(abs(landa(i,i)Maxlanda) Maxlanda=abs(landa(i,i); endend运行结果为：Maxlanda=0.6830即雅克比迭代矩阵半径=0.6830.0001) fprintf(%d %f %f %f %f %f %f %fn,iter,x(1),x(2),x(3),x(4),x(5),x(6),error); if(errorMaxlanda) Maxlanda=abs(landa(i,i); endend运行结果为Maxlanda = 0.4806即GAUSS-SEIDEL迭代
9、矩阵半径=0.48060.0001)fprintf(%d %f %f %f %f %f %f %fn,iter,x(1),x(2),x(3),x(4),x(5),x(6),error); if(errorMaxlanda) Maxlanda=abs(landa(j,j); end end if(Maxlandatol) x0=x; x=(D-w*L)(D*x0+w*(b-D*x0+U*x0); error=norm(x-x0); iter=iter+1;end在命令行输入程序：format short;A=4,-1,0,-1,0,0;-1,4,-1,0,-1,0;0,-1,4,-1,0,-1;
10、-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;tol=0.0001;x0=zeros(size(B);w=0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1;X,Iter,Error=sor(A,B,x0,0.1,tol);for i=2:length(w) x0=zeros(size(B); x,iter,error=sor(A,B,x0,w(i),tol); if(Iteriter) X=x; Iter=iter; W=w(i); Error=error; endend运行得到：I
11、ter =12 X =1.0000 2.0000 1.0000 2.0000 1.0000 2.0000Error =6.8688e-05 W =1.1000即SOR迭代收敛的条件下，当w*=1.1时，迭代次数最少为12，对应的近似解为x =1.0000 2.0000 1.0000 2.0000 1.0000 2.0000，误差为6.8688e-056用经典R-K方法求解初值问题（1）， ；（2）， 。和精确解比较，进行误差分析得到结论，图形显示精确解和数值解。解：1) 利用欧拉公式求解式（1）的数值解，然后将得到的数值解与精确解通过图形显示，步长设置为0.1，输入如下程序：clear;a=0;b=10;N=100; y1(1,:)=2,3;h=(b-a)/N;x=a:h:b;for n=1:N y1(n+1,1)=y1(n,1)+h*(-2*y1(n,1)+y1(n,2
