东南大学数值分析上机题答案

数值分析上机题第一章17.（上机题）舍入误差与有效数设NjNjS2211，其精确值为）111-23（21NN。（1）编制按从大到小的顺序1-1···1-311-21222NSN，计算NS的通用程序；（2）编制按从小到大的顺序121···1)1(111222NNSN，计算NS的通用程序；（3）按两种顺序分别计算210S,410S,610S,并指出有效位数（编制程序时用单精度）；（4）通过本上机题，你明白了什么？解：程序：（1）从大到小的顺序计算1-1···1-311-21222NSN：functionsn1=fromlarge(n)%从大到小计算sn1formatlong;sn1=single(0);form=2:1:nsn1=sn1+1/(m^2-1);endend（2）从小到大计算121···1)1(111222NNSNfunctionsn2=fromsmall(n)%从小到大计算sn2formatlong;sn2=single(0);form=n:-1:2sn2=sn2+1/(m^2-1);endend（3）总的编程程序为：functionp203()clearallformatlong;n=input('pleaseenteranumberasthen:')sn=1/2*(3/2-1/n-1/(n+1));%精确值为snfprintf('精确值为%f\n',sn);sn1=fromlarge(n);fprintf('从大到小计算的值为%f\n',sn1);sn2=fromsmall(n);fprintf('从小到大计算的值为%f\n',sn2);functionsn1=fromlarge(n)%从大到小计算sn1formatlong;sn1=single(0);form=2:1:nsn1=sn1+1/(m^2-1);endendfunctionsn2=fromsmall(n)%从小到大计算sn2formatlong;sn2=single(0);form=n:-1:2sn2=sn2+1/(m^2-1);endendend运行结果：从而可以得到N值真值顺序值有效位数2100.740050从大到小0.7400495从小到大0.74005064100.749900从大到小0.7498523从小到大0.74990066100.749999从大到小0.7498523从小到大0.7499996（4）感想：通过本上机题，我明白了，从小到大计算数值的精确位数比较高而且与真值较为接近，而从大到小计算数值的精确位数比较低。机器数在进行加法运算时，用从大到小的顺序容易出现大数吃小数的情况，容易产生较大的误差，是因为对于相加的两个数值，计算机首先提供与大数相一致的位数，此时将小数的尾数向右移位，并进行四舍五入，之后对尾数进行依次相加。从大到小时，越往后计算，相加的数越小，从而出现大数吃小数的情况。相比之下。从小到大计算时，每次小数与大数相加，都会增加位数，从而精确度比较高。第二章20.（上机题）Newton迭代法（1）给定初值0x及容许误差ε，编制Newton法解方程)(xf=0根的通用程序。（2）给定方程xxxf3/)(3，易知其有三个根1x=-3，2x=0,3x=3。①有Newton方法的局部收敛性可知存在δ0，当0x（-δ,δ）时Newton迭代序列收敛于根2x，试确定尽可能大的δ；②试取若干初始值，观察当0x（-∞，-1），（-1，-δ），（-δ，δ），（δ，1），（1，+∞）时Newton序列是否收敛以及收敛于哪一个根。（3）通过本上机题，你明白了什么？解：（1）程序先编写函数function文件：文件fx.m%定义函数f(x)functionFx=fx(x)Fx=x^3/3-x;文件dfx.m%定义导函数df(x)%functionfx=dfx(x)fx=x^2-1;接下来是具体步骤文件newton1.m求尽可能大的delta值%%课本56页计算最大delta值clearflag=1;k=1;x0=0;whileflag==1delta=k*10^-6;x0=delta;k=k+1;m=0;flag1=1;whileflag1==1&&m=10^3x1=x0-fx(x0)/dfx(x0);ifabs(x1-x0)10^-6flag1=0;endm=m+1;x0=x1;endifflag1==1||abs(x0)=10^-6flag=0;endendfprintf('%f\n',delta);文件newton2.m求方程的根%%课本56页newton法求方程的根，确定收敛于哪个根formatlong;ef=1e-6;k=0;x0=input('pleaseentertheinitialnumberasthex0:');whilek1000x1=x0-fx(x0)/dfx(x0);ifabs(x1-x0)efbreakendx0=x1;k=k+1;endfprintf('方程的根为%f\n',x0);（2）运行结果①求尽可能大的delta值②判断收敛于哪个根。已知有三个根，1x=-3，2x=0,3x=3。，（-1，-δ），（-δ，δ），（δ，1），（1，+∞）0x（-∞，-1）收敛于1x=-30x（-1，-δ）可见部分收敛于1x，部分收敛于3x。0x（-δ，δ）可见，收敛于2x=0。0x（δ，1）可见部分收敛于1x，部分收敛于3x。0x（1，∞）收敛于3x（3）感想：通过自行编写NEWTON的程序，加深了我对于NEWTON迭代法具体的运行方式、结果、利弊的了解。在不同的区间上，根据给定的x0的不同，x会以不同的速度收敛于某个根，也会产生意想不到的跳动，如在（1，∞）区间是。因此给定初值的重要性很大，即在某些区间上，收敛于某个根是由一定限制的，要取得合适的初值，才能够正确地求得收敛的根。总体上来说，通过上机，巩固学习了第二章关于迭代的运用。第三章39.列主元Gauss消去法对于某电路的分析，归结为求解线性方程组RIV。其中3113000100001335901100000931100000000107930000900030577050000074730000000030410000005002720009000229R15,27,23,0,20,12,7,7,10TTV（1）编制解n阶线性方程组Axb的列主元高斯消去法的通用程序；（2）用所编程序线性方程组RIV，并打印出解向量，保留5位有效数；（3）本章编程之中，你提高了哪些编程能力？解：（1）程序Gauss函数的function程序%用列主元Gauss消去法求解线性方程组Ax=bfunction[x]=gauss1(A,b)n=length(b);fork=1:na=max(A(k:n,k));[pm]=find(A(:,k)==a);ifpkA([p,k],:)=A([k,p],:);b([p,k],:)=b([k,p],:);endm=A(k+1:n,k)/A(k,k);A(k+1:n,k+1:n)=A(k+1:n,k+1:n)-m*A(k,k+1:n);b(k+1:n)=b(k+1:n)-m*b(k);A(k+1:n,k)=zeros(n-k,1);endx=zeros(n,1);x(n)=b(n)/A(n,n);fori=n-1:-1:1m=0;fork=i+1:1:nm=m+A(i,k)*x(k);endx(i)=(b(i)-m)/A(i,i);endend下面是主程序的m文件%%课本127页用列主元gauss消去法求解方程clearR=[31-13000-10000;-1335-90-110000;0-931-1000000;00-1079-30000-9;000-3057-70-50;0000-747-3000;00000-304100;0000-50027-2;000-9000-229];V=[-15;27;-23;0;-20;12;-7;7;10];I=gauss1(R,V);I（2）运行结果（3）感想：本次上机重点学习了列主元GUASS消去法的应用。列主元GUASS消去法在进行第K步消元前，先选出位于第K列中位于对角线及其以下元素绝对值中的最大者，然后将他们互相交换，在进行接下来的一般消元过程。较少了误差，而且一般保证舍入误差不增加，基本上是稳定的。通过上机程序设计，加深了对该方法的了解，对于矩阵、编程的了解也更深入。第四章37.三次样条插值函数（1）编制求第一型3次样条插值函数的通用程序（2）已知汽车门曲线型值点的数据如下i012345678910ix012345678910iy2.513.304.044.705.225.545.785.405.575.705.80端点条件为,0y=0.8，,10y=0.2，用所编程序求出门的三次样条差值函数S(x)，打印出S(i+0.5),i=1,2,```9。（1）程序（n是x的总个数）%%样条插值函数clearclc%%输入相关参数xi,yin=input('enterthen:');%%输入X的总个数xn=zeros(1,n);yn=zeros(1,n);xn(1,:)=input('enterthex:');yn(1,:)=input('enterthey:');%%求h,mu,lambda,d值d=zeros(n,1);h=zeros(1,n-1);mu=zeros(1,n-2);%%lambda=zeros(1,n-2);f1=zeros(1,n-1);%%一阶导数f2=zeros(1,n-2);%%二阶导数dy0=input('enterthevalueofdy0:');dyn=input('enterthevalueofdyn:');fori=1:n-1h(i)=xn(i+1)-xn(i);f1(i)=(yn(i+1)-yn(i))/h(i);endfori=1:n-2mu(i)=h(i)/(h(i)+h(i+1));lambda(i)=1-mu(i);endd(1)=6*(f1(1)-dy0)/h(1);d(n)=6*(dyn-f1(n-1))/h(n-1);fori=2:n-1f2(i)=(f1(i)-f1(i-1))/(xn(i+1)-xn(i-1));d(i)=6*f2(i);endA=zeros(n);A(1,2)=1;A(n,n-1)=1;fori=1:nA(i,i)=2;endfori=2;n-1A(i,i-1)=mu(i-1);A(i,i+1)=lambda(i-1);endM=A\d;%%回代求插值函数symsx;disp('sn(x)=')%%sn(i)的表达式sn(i)仅是关于x的函数值fori=1:n-1sn(i)=collect(yn(i)+(f1(i)-(M(i)/3+M(i+1)/6)*h(i))*(x-xn(i))+M(i)/2*(x-xn(i))^2+(M(i+1)-M(i))/(6*h(i))*(x-xn(i))^3);sn(i)=vpa(sn(i),4);fprintf('sn(%d)=%s,%d¡Üxn¡Ü%d\n',i,char(sn(i)),xn(i),xn(i+1));enddisp('s(i+0.5)');%%计算在i+0.5处的插值disp('ix(i+0.5)s(i+0.5)')fori=1:n-1t=xn(i)+0.5;s(i)=yn(i)+(f1(i)-(M(i)/3+M(i+1)/6)*h(i))*(t-xn(i))+M(i)/2*(t-xn(i))^2+(M(i+1)-M(i))/(6*h(i))*(t-xn(i))^3;fprintf('%d%.4f%.4f\n',i,xn(i)+0.5,s(i));