《数值计算方法上机实习题(14页).docx》由会员分享，可在线阅读，更多相关《数值计算方法上机实习题(14页).docx（13页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。
1、-数值计算方法上机实习题-第 13 页数值计算方法上机实习题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
2、.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） 分析结果的可靠性
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
4、),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
5、计算根与步数程序：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)
6、,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);m
7、ax=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)=
8、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）
9、分析绝对误差。答：可以看到，在同一精度下，二分法运算了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.5
10、8 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为拟合曲
11、线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
12、=%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设，分析下列迭代法的收敛性，并求的近似解及相应的迭
13、代次数。（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; endendJac
14、obi迭代计算程序：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.9
15、9998177 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 no
16、rm(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.8f
17、n%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);f
18、printf(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法，其设计思想是将线性方程组
19、的求解归结为对角方程组求解过程的重复，计算规则简单而易于编写计算程序。但收敛速度慢，谱半径相对偏大（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;tma
20、x,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=NMS
21、F2(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)
22、+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
23、,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)
24、-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,ER
25、2);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
26、（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的误差小，较适合用在此类钢性问题中。计算结果显示：(1) Ode 23的题1(2) Ode 3的题2
27、(3) Ode 45的题1(4) Ode 45的题27用有限差分法求解边值问题（h=0.1），并图形显示。有限差分法子程序：function y=YXCFF_D(h)n=(1-(-1)/h;x(1)=-1;for i=1:n-1 % n-1 x(i)=x(1)+(i-1)*h; q(i)=(1+x(i)2); B(i)=2/(h2)+q(i);end for i=1:n-2 % n-2 C(i)=-1/(h2);endb(1)=0+1/(h2); b(n-1)=0+1/(h2); for i=2:n-2 b(i)=0;endb=b;% n-1 u(1)=b(1)/B(1);v(1)=C(1)/
28、B(1);for k=2:n-2; u(k)=(b(k)-u(k-1)*C(k-1)/(B(k)-v(k-1)*C(k-1); v(k)=C(k)/(B(k)-v(k-1)*C(k-1);endu(n-1)=(b(n-1)-u(n-2)*C(n-2)/(B(n-1)-v(n-2)*C(n-2);V(n-1)=u(n-1);for k=n-2:-1:1 V(k)=u(k)-v(k)*V(k+1);endy=1,V,1;计算与显示程序：h=0.0001;h1=0.1;L=YXCFF_D(h);L1=YXCFF_D(h1);X=-1:h:1;X1=-1:h1:1;figure1=figure(Color,1 1 1);gxt=plot(X,L,k-,X1,L1,r-);title(有限差分法,FontSize,14,FontWeight,Bold);xlabel(输入值x,FontSize,14,FontWeight,Bold);ylabel(值y,FontSize,14,FontWeight,Bold);set(gxt,LineWidth,1.5);grid on;总结：对比两个步长下，求边值问题，步长小，越精确。