数值计算方法上机实习题.doc

数值计算方法上机实习题1 设，（1） 由递推公式，从I0=0.1824, 出发，计算；（2） ，, 用，计算；(1)由I0计算I20递推子程序：function f=fib(n,i)if n=1 f=fib(n-1,i)*(-5)+(1/(n);elseif n=0 f=i;end计算和显示程序：I=0.1824;I1=0.1823;fib1=fib(20,I);fib2=fib(20,I1);fprintf(I_0=0.1824时，I_20=%dn,fib1);fprintf(I_0=0.1823时，I_20=%dn,fib2);计算结果显示：I_0=0.1824时，I_20=7.480927e+09I_0=0.1823时，I_20=-2.055816e+09(2)由I20计算I0程序：n=21;i1=0;i2=10000;f1=i1;f2=i2;while n=0; f1=f1*(-1/5)+(1/(5*n); f2=f2*(-1/5)+(1/(5*n); n=n-1;endfprintf(I_20=0 时，I_0=%4.8fn,f1);fprintf(I_20=10000时，I_0=%4.8fn,f2);计算结果显示：I_20=0 时，I_0=0.18232156I_20=10000时，I_0=0.18232156（3） 分析结果的可靠性及产生此现象的原因(重点分析原因)。答：第一个算法可得出e0=I0-I0*en=In-In*=5ne0易知第一个算法每一步计算都把误差放大了5倍，n次计算后更是放大了5n倍，可靠性低。 第二个算法可得出en=In-In*e0=15nen可以看出第二个算法每一步计算就把误差缩小5倍，n次后缩小了5n倍，可靠性高。2 求方程的近似根，要求，并比较计算量。（1） 在0，1上用二分法；(1) 0，1上的二分法二分法子程序：function root,n=EFF3(f,x1,x2)%第二题(1)二分法f1=subs(f,symvar(f),x1);%函数在x=x1的值f2=subs(f,symvar(f),x2);%x=x2n=0;%步数er=5*10-4;%误差if(f1=0) root=x1; return;elseif(f2=0) root=x2; return;elseif(f1*f20) disp(两端点函数值乘积大于0！); return;else while(abs(x1-x2)er)%循环 x3=(x1+x2)/2; f3=subs(f,symvar(f),x3); n=n+1; if(f3=0) root=x3; break; elseif(f1*f30) x1=x3; else x2=x3; end end root=(x1+x2)/2;%while循环少一步需加上end计算根与步数程序：fplot(x) exp(x)+10*x-2,0,1);grid on;syms x;f=exp(x)+10*x-2;root,n=EFF3(f,0,1);fprintf(root=%6.8f ,n=%d n,root,n);计算结果显示：root=0.09057617 ,n=11 （2） 取初值，并用迭代；(2) 初值x0=0迭代迭代法子程序：function root,n=DDF(g,x0,err,max) （接下页）%root根，n+1步数，g函数，x0初值，err误差，max最大迭代次数X(1)=x0;for n=2:max X(n)=subs(g,symvar(g),X(n-1); c=abs(X(n)-X(n-1); root=X(n); if(cerr) break; end if n=max disp(超出预设迭代次数); endendn=n-1;%步数减1计算根与步数程序：syms x;f=(2-exp(x)/10; （接下页）x0=0; err=5*10(-4);max=100;root,n=DDF(f,x0,err,max);fprintf(root=%6.8f ,n=%d n,root,n);计算结果显示：root=0.09051262 ,n=4（3） 加速迭代的结果;(3) 加速迭代加速迭代计算程序：x0=0;err=5*10(-4);max=100;syms x;g=x-(f(x)-x)2/(f(f(x)-2*f(x)-x);root,n=DDF(g,x0,err,max);fprintf(root=%6.8f, n=%d,root,n);程序函数设置：function g=f(x)g=(2-exp(x)/10;计算结果显示：root=0.09048375, n=2（4） 取初值，并用牛顿迭代法；(4) 牛顿迭代法牛顿迭代法子程序：function root,n=NDDDFx(g,x0,err,max)%root根，n+1步数，g函数，x0初值，err误差，max最大迭代次数X(1)=x0;for n=2:max X(n)=subs(g,symvar(g),X(n-1); c=abs(X(n)-X(n-1); root=X(n); if(cerr)|(root=0) break; end if n=max disp(超出预设迭代次数); （接下页） endendn=n-1;%步数减1牛顿迭代法计算程序：x0=0;err=5*10(-4);max=100;syms x;g=exp(x)+10*x-2;g=x-(g/diff(g);root,n=NDDDFx(g,x0,err,max);fprintf(root=%6.8f, n=%d n,root,n)；计算结果显示：root=0.09052511, n=2（5） 分析绝对误差。答：可以看到，在同一精度下，二分法运算了11次，题设迭代算式下运算了4次，加速迭代下运算了2次，牛顿迭代下运算了2次。因不动点迭代法和二分法都是线性收敛的，但二分法压缩系数比题设迭代方法大，收敛速度较慢。加速迭代速度是超线性收敛，牛顿法是二阶，收敛速度快。3钢水包使用次数多以后，钢包的容积增大，数据如下：x23456789y6.428.29.589.59.7109.939.991011121314151610.4910.5910.6010.810.610.910.76试从中找出使用次数和容积之间的关系，计算均方差。（用拟合）拟合曲线程序：x=2: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;f,fval=fminsearch(delta1,0,0,0,optimset,x,y);%fminsearch为求解多元函数的最小值函数%f为多元函数初值x0附近的极小值点%fval为极小值k=2:100;y1=(f(1).*k+f(2)./(f(3)+k);figure1=figure(Color,1 1 1);gxt=plot(x,y,xr,k,y1,k-);legend(原数据,拟合曲线,Location,northwest); %y为数据点连接曲线，y1为拟合曲线title(函数y=(ax+b)/(x+c)的拟合,FontSize,14,FontWeight,Bold);xlabel(次数,FontSize,14,FontWeight,Bold);ylabel(容积,FontSize,14,FontWeight,Bold);set(gxt,LineWidth,1.5);grid on;%计算均方差for i=1:15y2(i)=(f(1).*x(i)+f(2)./(f(3)+x(i); （接下页）endj=0;for i=1:15j=j+(y(i)-y2(i)2;endjfc=sqrt(j/15);fprintf(拟合出的方程为:(x+%4.4f)y=%4.4fx+%4.4f n均方差为:%4.8f n,f(3),f(1),f(2),jfc);构造函数子程序：function delta=delta1(f,x,y)a=f(1);b=f(2);c=f(3);delta=0;for k=1:15 delta=delta+(x(k)+c)*y(k)-(a*x(k)+b)2;end计算结果显示：拟合出的方程为:(x+-0.7110)y=11.2657x+-15.5024 均方差为:0.33165089总结：指标选择，因题设方程为非线性的，要转化为线性方程故需提指标为：其驻点方程为:计算结果显示：4设，分析下列迭代法的收敛性，并求的近似解及相应的迭代次数。（1） JACOBI迭代；(1) Jacobi迭代Jacobi迭代子程序：function x,k=Jacobifl2(A,b,x0,eps,max1)%jacobi按矩阵的分量算法%x0初值，eps误差，max1最大迭代次数n=length(x0);for k=1:max1 x(1)=(b(1)-A(1,2:n)*x0(2:n)/A(1,1); for i=2:n x(i)=(b(i)-A(i,1:i-1)*x0(1:i-1)-A(i,i+1:n)*x0(i+1:n)/A(i,i); end if norm(x-x0,2)eps return else x0=x; endendJacobi迭代计算程序：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;x0=0 0 0 0 0 0;eps=10(-4);max1=200;X,k=Jacobifl2(A,b,x0,eps,max1);fprintf(近似解x=n%4.8fn %4.8f n %4.8f n %4.8fn %4.8f n %4.8fn迭代次数n=%dn,X(1),X(2),X(3),X(4),X(5),X(6),k);计算结果显示：近似解x=0.99998177 1.99995018 0.99997509 1.99995018 0.99997509 1.99996353迭代次数n=28（2） GAUSS-SEIDEL迭代；(2) GAUSS-SEIDEL迭代GAUSS-SEIDEL迭代子程序：function x,k=GSDD(A,b,x0,eps,max1)n=length(x0);for k=1:max1 x(1)=(b(1)-A(1,2:n)*x0(2:n)/A(1,1); for i=2:n x(i)=(b(i)-A(i,1:i-1)*x(1:i-1)-A(i,i+1:n)*x0(i+1:n)/A(i,i); end if norm(x-x0,2)=eps x0=x; x(1)=x0(1)+w(s)*(b(1)-A(1,1:n)*x0(1:n)/A(1,1); for i=2:n x(i)=x0(i)+w(s)*(b(i)-A(i,1:i-1)*x(1:i-1)-A(i,i:n)*x0(i:n)/A(i,i); end k(s)=k(s)+1; if (k(s)=max1) disp(超出迭代次数上限！n); return; end end Datas=x;endmin1,index=min(k);x=Dataindex;k=min1;m=w(index);fprintf(n解x=n%4.8fn%4.8fn%4.8fn%4.8fn%4.8fn%4.8f,x);fprintf(n步数最小omega=%2.2f, 步数n=%d 。n,m,k);SOR迭代计算程序：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;x0=0 0 0 0 0 0;eps=10(-4);max1=1000;X,k,m=SOR1(A,b,x0,eps,max1);fprintf(n近似解x=n%4.8fn%4.8fn%4.8fn%4.8fn%4.8fn%4.8f,x);fprintf(n步数最小omega=%2.1f, 步数n=%dn,m,k);计算结果显示：近似解x=1.000009771.999997820.999995672.000010710.999996621.99999923步数最小omega=1.20, 步数n=9 。总结：Jacobi迭代次数为28次；GAUSS-SEIDEL迭代次数为15次；SOR迭代次数在w=1.2时达到最小次数9次。系数矩阵严格对角占优，则Jacobi迭代法与Gauss-Seidel迭代法均收敛；又因系数矩阵对称正定，且0w2，则SOR迭代法也收敛。三种迭代法的结果分析及比较：（1）Jacobi法，其设计思想是将线性方程组的求解归结为对角方程组求解过程的重复，计算规则简单而易于编写计算程序。但收敛速度慢，谱半径相对偏大（2）Gauss-Seidel充分利用新信息进行计算，其迭代的效果比Jacobi迭代好，谱半径比Jacobi法小。（3）SOR迭代法，计算公式简单、编制程序容易，是求解大型稀疏方程组的一种有效方法，且若松弛因子选择得当超松弛法收敛速度比Gauss-Seidel迭代和Jacobi迭代都要快。 5用逆幂迭代法求最接近于11的特征值和特征向量，准确到。逆幂迭代法子程序：function lamda,V,k=NMSF2(A,v,eps,p)n,n=size(A);A=A-p*eye(n);v0=v;tmax,tindex=max(abs(v0);lamd0=v0(tindex);u0=v0/lamd0;flag=0;k=0;while(flag=0) V=Au0; tmax,tindex=max(abs(V); lamd1=V(tindex); u0=V/lamd1; if (abs(lamd0)(-1)-(lamd1)(-1)=eps flag=1; end lamd0=lamd1; k=k+1;end （接下页）lamda=(lamd1)(-1)+p;V=u0;逆幂迭代计算程序：A=6 3 1;3 2 1;1 1 1;v=1 1 1;eps=0.001;p=11;lamda,U,k=NMSF2(A,v,eps,p);fprintf(最接近11的特征值为:%4.8fn特征向量：n%4.8fn %4.8fn%4.8fn,lamda,U(1),U(2),U(3);计算结果显示：最接近于11的特征值为：7.87313712特征向量： 1.00000000 0.54922698 0.225456186用经典R-K方法求解初值问题（1）， ；（2）， 。和精确解比较，进行误差分析得到结论，图形显示精确解和数值解。经典四阶R-K方法R-K程序：f=(x,y) -2*y(1)+y(2)+2*sin(x);y(1)-2*y(2)+2*cos(x)-2*sin(x);g=(x,y) -2*y(1)+y(2)+2*sin(x);998*y(1)-999*y(2)+999*cos(x)-999*sin(x); （接下页）x=0:0.1:10;Y1=2*exp(-x)+sin(x);Y2=2*exp(-x)+cos(x);x1=0:10/10000:10;Y3=2*exp(-x1)+sin(x1);Y4=2*exp(-x1)+cos(x1);N=100;N1=10000;a1,b1=ode23(f,0:0.1:10,2 3);a2,b2=ode23(g,0:1/1000:10,2 3);a3,b3=ode45(f,0:0.1:10,2 3);a4,b4=ode45(g,0:1/1000:10,2 3);n=length(b1);figure1=figure(Color,1 1 1);dbt=plot(x,b1(:,1),r-,x,b1(:,2),k-,x,Y1,g-,x,Y2,y-);legend(y1,y2,精确解Y1,精确解Y2); title(四阶RK算法下的方程组的数值解&方程的精确解,FontSize,14,FontWeight,Bold);ylabel(Y,FontSize,14,FontWeight,Bold);xlabel(X,FontSize,14,FontWeight,Bold);set(dbt,LineWidth,1.5);grid on;Q=b1(:,1)-Y1;W=b1(:,2)-Y2;Q1=b2(:,1)-Y3;W1=b2(:,2)-Y4;Q2=b3(:,1)-Y1;W2=b3(:,2)-Y2;Q3=b4(:,1)-Y3;W3=b4(:,2)-Y4; ER1=max(abs(Q);ER2=max(abs(W);ER3=max(abs(Q1);ER4=max(abs(W1);ER5=max(abs(Q2);ER6=max(abs(W2);ER7=max(abs(Q3);ER8=max(abs(W3);fprintf(n（1）中y1在ode23下求出的最大误差为：%4.8fn（1）中y2在ode23下求出的最大误差为：%4.8fn,ER1,ER2);fprintf(n（2）中y1在ode23下求出的最大误差为：%4.8fn（2）中y2在ode23下求出的最大误差为：%4.8fn,ER3,ER4);fprintf(n（1）中y1在ode45下求出的最大误差为：%4.8fn（1）中y2在ode45下求出的最大误差为：%4.8fn,ER5,ER6);fprintf(n（2）中y1在ode45下求出的最大误差为：%4.8fn（2）中y2在ode45下求出的最大误差为：%4.8fn,ER7,ER8);计算与显示程序：（1）中y1在ode23下求出的最大误差为：0.00169850（1）中y2在ode23下求出的最大误差为：0.00207422（2）中y1在ode23下求出的最大误差为：0.00000583（2）中y2在ode23下求出的最大误差为：0.00581329（1）中y1在ode45下求出的最大误差为：0.00052155（1）中y2在ode45下求出的最大误差为：0.00045793（2）中y1在ode45下求出的最大误差为：0.00000276（2）中y2在ode45下求出的最大误差为：0.00275331总结：对比2种方程在ode23，ode45命令下，与精确值的误差可以发现，在此题中，ode45的误差相对于ode23的误差小，