数值计算方法上机实习题

数值计算方法上机实习题1．设105dxxxInn，（1）由递推公式nIInn151，从I0=0.1824,0=0.1823I出发，计算20I；（2）20=0I，20=10000I,用nIInn515111，计算0I；(1)由I0计算I20递推子程序：functionf=fib(n,i)ifn=1f=fib(n-1,i)*(-5)+(1/(n));elseifn==0f=i;end计算和显示程序：I=0.1824;I1=0.1823;fib1=fib(20,I);fib2=fib(20,I1);fprintf('I_0=0.1824时，I_20=%d\n',fib1);fprintf('I_0=0.1823时，I_20=%d\n',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;whilen~=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.8f\n',f1);fprintf('I_20=10000时，I_0=%4.8f\n',f2);计算结果显示：I_20=0时，I_0=0.18232156I_20=10000时，I_0=0.18232156（3）分析结果的可靠性及产生此现象的原因(重点分析原因)。答：第一个算法可得出𝑒0=|𝐼0−𝐼0∗|𝑒𝑛=|𝐼𝑛−𝐼𝑛∗|=5𝑛|𝑒0|易知第一个算法每一步计算都把误差放大了5倍，n次计算后更是放大了5n倍，可靠性低。第二个算法可得出e𝑛=|𝐼𝑛−𝐼𝑛∗|𝑒0=(15)𝑛|𝑒𝑛|可以看出第二个算法每一步计算就把误差缩小5倍，n次后缩小了5n倍，可靠性高。2．求方程0210xex的近似根，要求41105kkxx，并比较计算量。（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;elsewhile(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;elsex2=x3;endendroot=(x1+x2)/2;%while循环少一步需加上end计算根与步数程序：fplot(@(x)exp(x)+10*x-2,[0,1]);gridon;symsx;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）取初值00x，并用迭代1021xkex；(2)初值x0=0迭代迭代法子程序：function[root,n]=DDF(g,x0,err,max)（接下页）计算根与步数程序：symsx;f=(2-exp(x))/10;（接下页）%root根，n+1步数，g函数，x0初值，err误差，max最大迭代次数X(1)=x0;forn=2:maxX(n)=subs(g,symvar(g),X(n-1));c=abs(X(n)-X(n-1));root=X(n);if(cerr)break;endifn==maxdisp('超出预设迭代次数');endendn=n-1;%步数减1x0=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;symsx;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);程序函数设置：functiong=f(x)g=(2-exp(x))/10;计算结果显示：root=0.09048375,n=2（4）取初值00x，并用牛顿迭代法；(4)牛顿迭代法牛顿迭代法子程序：function[root,n]=NDDDFx(g,x0,err,max)%root根，n+1步数，g函数，x0初值，err误差，max最大迭代次数X(1)=x0;forn=2:maxX(n)=subs(g,symvar(g),X(n-1));c=abs(X(n)-X(n-1));root=X(n);if(cerr)||(root==0)break;endifn==maxdisp('超出预设迭代次数');（接下页）牛顿迭代法计算程序：x0=0;err=5*10^(-4);max=100;symsx;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=2endendn=n-1;%步数减1（5）分析绝对误差。答：可以看到，在同一精度下，二分法运算了11次，题设迭代算式下运算了4次，加速迭代下运算了2次，牛顿迭代下运算了2次。因不动点迭代法和二分法都是线性收敛的，但二分法压缩系数比题设迭代方法大，收敛速度较慢。加速迭代速度是超线性收敛，牛顿法是二阶，收敛速度快。3．钢水包使用次数多以后，钢包的容积增大，数据如下：x23456789y6.428.29.589.59.7109.939.991011121314151610.4910.5910.6010.810.610.910.76试从中找出使用次数和容积之间的关系，计算均方差。（用axbycx拟合）拟合曲线程序：x=2:16;y=[6.428.29.589.59.7109.939.9910.4910.5910.6010.810.610.910.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',[111]);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);gridon;%计算均方差fori=1:15构造函数子程序：functiondelta=delta1(f,x,y)a=f(1);b=f(2);c=f(3);delta=0;fork=1:15delta=delta+((x(k)+c)*y(k)-(a*x(k)+b))^2;end计算结果显示：拟合出的方程为:(x+-0.7110)y=11.2657x+-15.5024均方差为:0.33165089总结：指标选择，因题设方程为非线性的，要转化为线性方程故需提指标为：y2(i)=(f(1).*x(i)+f(2))./(f(3)+x(i));（接下页）endj=0;fori=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);其驻点方程为:计算结果显示：4．设410100141010014101101410010141001014A，625250b，bxA分析下列迭代法的收敛性，并求42110kkxx的近似解及相应的迭代次数。（1）JACOBI迭代；(1)Jacobi迭代Jacobi迭代子程序：function[x,k]=Jacobifl2(A,b,x0,eps,max1)%jacobi按矩阵的分量算法%x0初值，eps误差，max1最大迭代次数n=length(x0);fork=1:max1x(1)=(b(1)-A(1,2:n)*x0(2:n))/A(1,1);fori=2:nx(i)=(b(i)-A(i,1:i-1)*x0(1:i-1)-A(i,i+1:n)*x0(i+1:n))/A(i,i);endifnorm(x'-x0,2)epsreturnelsex0=x';endendJacobi迭代计算程序：A=[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]';eps=10^(-4);max1=200;[X,k]=Jacobifl2(A,b,x0,eps,max1);fprintf('近似解x=\n[%4.8f\n%4.8f\n%4.8f\n%4.8f\n%4.8f\n%4.8f]\n迭代次数n=%d\n',X(1),X(2),X(3),X(4),X(5),X(6),k);计算结果显示：近似解x=[0.999981771.999950180.999975091.999950180.999975091.99996353]迭代次数n=28（2）GAUSS-SEIDEL迭代；(2)GAUSS-SEIDEL迭代GAUSS-SEIDEL迭代子程序：function[x,k]=GSDD(A,b,x0,eps,max1)n=length(x0);fork=1:max1x(1)=(b(1)-A(1,2:n)*x0(2:n)')/A(1,1);fori=2:nx(i)=(b(i)-A(i,1:i-1)*x(1:i-1)'-A(i,i+1:n)*x0(i+1:n)')/A(i,i);GAUSS-SEIDEL迭代计算程序：A=[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];eps=10^(-4);max1=200;[X,k]=GSDD(A,b,x0,eps,max1);fprintf('近似解x=\n[%4.8f\n%4.8f\n%4.8fendifnorm(x-x0,2)epsreturnelsex0=x;endend\n%4.8f\n%4.8f\n%4.8f]\n迭代次数n=%d\