西安交通大学计算方法上机课后复习

^算方法上机作业.对以下和式计算：S—Xf上一上一，一，］，要求:n=016n18n+18n+48n+58〃+6^⑴若只需保留11个有效数字，该如何进行计算；⑵若要保留30个有效数字，则又将如何进行计算；(1)解题思想和算法实现:根据保留有效位数的要求，可以由公式：三12E一得出计算精度要求。只需要很少内存，时间复杂度和d呈线性，不需要高浮点支持。先根据while语句求出符合精度要求的n值的大小，然后利用for语句对这n项进行求和，输出计算结果及n值大小即可。2)matlab源程序:保留11位有效数字时;clearclcformatlongn=0;sum=1/(16An)*(4/(8*n+1)-2/(8*n+4)-1/(8*n+5)-1/(8*n+6));whilesum>=5*10A(-11)；n=n+1;sum=1/(16An)*(4/(8*n+1)-2/(8*n+4)-1/(8*n+5)-1/(8*n+6));endfori=0:n-1;sum=sum+1/(16Ai)*(4/(8*i+1)-2/(8*i+4)-1/(8*i+5)-1/(8*i+6));endvpa(sum,11)n保留30位有效数字时；clearclcformatlongn=0;sum=1/(16An)*(4/(8*n+1)-2/(8*n+4)-1/(8*n+5)-1/(8*n+6));whilesum>=5*10A(-30);n=n+1;sum=1/(16An)*(4/(8*n+1)-2/(8*n+4)-1/(8*n+5)-1/(8*n+6));endfori=0:n-1;sum=sum+1/(16Ai)*(4/(8*i+1)-2/(8*i+4)-1/(8*i+5)-1/(8*i+6));endvpa(sum,30)n(3)实验结果分析anw=she= 3.L4159265358379323846204r33832S3A»I士»I图1.1保留11位有效数字的n值及计算结果图 图1.2保留30位有效数字的n值及计算结果图由计算结果可知，通过合理的误差控制，分别通过7次和22次循环，可以实现题目所要求的精确度。2.某通信公司在一次施工中，需要在水面宽度为20米的河沟底部沿直线走向铺设一条沟底光缆。在铺设光缆之前需要对沟底的地形进行初步探测，从而估计所需光缆的长度，为工程预算提供依据。已探测到一组等分点位置的深度数据（单位:米）如下表所示：^点0123456深度9.018.967.967.978.029.0510.13^点78910111213深度11.1812.2613.2813.3212.6111.2910.22^点14151617181920深度9.157.907.958.869.8110.810.9⑵预测所需光缆长度的近似值，并作出铺设河底光缆的曲线图;(1)解题思想和算法原理给定区间［a,b］一个分划/:a=x0<x1<_<xN=b若函数S(x)满足下列条件:S(X)在每个区间［XyXj］上是不高于3次的多项式。5仅汲其2阶导数在［a,b］上连续。则称S(x)使关于分划/的三次样条函数。(2)matlab源程序:clc,clearx=0:1:20;y=[9.018.967.967.978.029.0510.1311.1812.2613.2813.3212.6111.2910.229.157.97.958.869.8110.8010.93];n=length(x);l(1)=0;m(n)=0;h=diff(x);df=diff(y)./diff(x);d(1)=0;d(n)=0;forj=2:n-1l(j)=h(j)/(h(j-1)+h(j));m(j)=h(j-1)/(h(j-1)+h(j));d(j)=6*(df(j)-df(j-1))/(h(j-1)+h(j));endm=m(2:end);u=diag(m,-1);r=diag(l,1);a=diag(2*ones(1,n));A=u+r+a;M=inv(A)*d';symsgforj=1:n-1s(j)=M(j)*(x(j+1)-g)八3/(6*h(j))+M(j+1)*((g-x(j))八3/(6*h(j)))+(y(j)-M(j)*h(j)八2/6)*(x(j+1)-g)/h(j)+(y(j+1)-M(j+1)*h(j)八2/6)*(g-x(j))/h(j);endsr=0;forj=1:n-1df=diff(s(j),g);warningoffall;q=int(sqrt(1+df.A2),g,j-1,j);r=r+q;endL=vpa(r,8);disp('thelengthofthelabelisL=');disp(L);forj=1:n-1S(j,:)=sym2poly(s(j));endforj=1:n-1x1=x(j):0.1:x(j+1);y1=polyval(S(j,:),x1);ifj==1y2=y1;elsefori=1:11k=(j-1)*10+i;y2(k)=y1(i);endendendx2=x(1):0.1:x(n);plot(x,y,'o')gridholdonplot(x2,y2,'r')(3)实验结果分析

[(3958520385365547*e)/14073748835532800-(186483313085687^"3)/562919953421312■+901/100.(5067619118220011*(g-thelengthafthelabelisL=26.498:514图2.1铺设河底电缆长度图2.2铺设河底光缆的曲线图由三次样条插值得出的函数曲线的长度和即铺设河底电缆的长度为26.498514。为了提高插值精度，用三次样条插值可以增加插值节点的办法来满足要求，而且在给定节点数的条件下，三次样条插值的精度要优于多项式插值以及线性分段插值，虽然舍弃了降低误差这个优点，但是其曲线的光滑性要好一些。3.假定某天的气温变化记录如下表所示，试用数据拟合的方法找出这一天的气温变化的规律;试计算这一天的平均气温，并试估计误差。时刻012345678911111111222平均气温5144144156180203258111122时刻314516718920122324332221平均气温134129725422018716(1)解题思想和数学原理：对于具体实验时，通常不是先给出函数的解析式，再进行实验，而是通过实验的观察和测量给出离散的一些点，再来求出具体的函数解析式。又因为测量误差的存在，实际真实的解析式曲线并不一定通过测量给出的所有点。最小二乘法，形成法方程是求解这一问题的很好的方法，本实验运用这一方法实现对给定数据的拟合。对于给定的测量数据(xfi)(i=12…，n)，设函数分布为y(x)=Ea①(x)j=o特别的，取甲j(x)为多项式①C)=xj(j=0,1,...，m)则根据最小二乘法原理，可以构造泛函i=1ji=1(k=0,1,..,m)(f，P(f，P0)一30,50)

(P4)

0:1■(P4)0m求该解方程组，则可以得到解a0,a1，…,am，因此可得到数据的最小二乘解f(x)xEa①(x)j=0(2)matlab源程序：x=[0123456789101112131415161718192021222324];%给出题目数据(时间)y=[15141414141516182020232528313231292725242220181716]；%给出题目数据(温度)plot(x,y,'m*')%画出各个离散数据点holdonforn=2:4;%2、3、4代表拟合函数最高项次数alltemp=25;%alltemp代表数据点总共有25个A=zeros(n+1,n+1);%定义初始正规方程组的系数矩阵AC=ones(n+1,1);%定义初始正规方程组的系数向量CD=zeros(n+1,1);%定义初始正规方程组的向量Dfori=1:n+1forj=1:n+1fork1=1:alltempA(i,j)=A(i,j)+(x(k1).人(i-1+j-1));endendfork2=1:alltempD(i,1)=D(i,1)+(x(k2).人(i-1)).*(y(k2));endend%以上为计算出正规方程组矩阵A、D的所有元素的程序tol=1.0e-12;maxit=1000;C=bicg(A,D,tol,maxit);%使用bicg迭代算出正规方程组的系数向量Cp=0;%误差分量E=0;%误差总量ifn==2b=0:24;f=C(1)+C(2).*b+C(3).*(b.人2);p=y(b+1)-f;forv=1:25E=E+(p(v)).人2;endplot(b,f,'r-')end%以上是对2阶拟合函数的图形处理与误差计算ifn==3b=0:24;f=C(1)+C(2).*b+C(3).*(b.人2)+C(4).*(b.人3);p=y(b+1)-f;forv=1:25E=E+(p(v)).人2;endplot(b,f,'g-')end%以上是对3阶拟合函数的图形处理与误差计算ifn==4b=0:24;f=C(1)+C(2).*b+C(3).*(b.A2)+C(4).*(b.A3)+C(5).*(b.A4);p=y(b+1)-f;forv=1:25E=E+(p(v)).A2;endplot(b,f,'b-')end%以上是对4阶拟合函数的图形处理与误差计算C,Eendn=2;%重新对n赋值，进行指数函数拟合A=zeros(n+1,n+1);%重新对A矩阵赋初值C=zeros(n+1,1);%重新对C向量赋初值D=zeros(n+1,1);%重新对D向量赋初值fori=1:n+1forj=1:n+1fork=1:alltempA(i,j)=A(i,j)+(x(k).人(i-1+j-1));endendforl=1:alltempD(i,1)=D(i,1)+((x(l).人(i-1)).*(log(y(l))));endend%计算出A矩阵、D向量各元素数值C=bicg(A,D,tol,maxit);%利用bicg迭代求解系数b=0:24;p=0;E=0;f=exp(C(1)+C(2).*b+C(3).*(b.人2));p=y(b+1)-f;forv=1:25E=E+(p(v)).人2;endplot(b,f,'c-')%对指数拟合函数进行图形处理和误差计算b=-C⑶；c-C(2)/(2*b);a=exp(b*(c^2)+C(1));%算出题设要求的指数拟合函数的各个系数a,b,c,Egridonlegend('测量数据',‘二次函数‘,‘三次函数‘,‘四次函数’,’指数拟合,'Location','NorthWest')holdoff%holdon与holdoff联合使用，表示将各个曲线画在同一个图中8.42742.5606-0.0920E=253.6476图3.1二次曲线拟合系数与2范数误差C=13.4072-32161CL2032-0.0082E=110.2954图3.2三次曲线拟合系数与2范数误差-3.55470.S592-0.0E130.0009F=43.9299图3.3四次曲线拟合系数与2范数误差25.024Tb=0.004414.09139E=191a0771图3.4指数曲线拟合系数与2范数误差图3.5数据原始点与拟合曲线对比图(3)实验结果分析：根据国家有关规定，用的是02、08、14、20时四个观测时次的数据做平均，最有代表性。从图中可以看出并不是多项式次数越高越好，随着次数的增高，曲线所呈现出的给定点处和实际的吻合度越好，但对于其他地方的吻合度降低了。4.设计算法，求出非线性方程6x5_45x2+20-0的所有根，并使误差不超过10-4。解：(1)解题思想和算法实现：对于一个非线性方程的数值解法很多。在此介绍两种最常见的方法：二分法和Newton法。首先要研究函数的形态，确定根的数量和大致区间的位置。对于二分法，其数学实质就是说对于给定的待求解的方程f(x)，其在［4切上连续，f(a)f(b)<0,且f(x)在［a,b］内仅有一个实根x*，取区间中点c，若，则c恰为其根，否则根据f(a)f(c)<0是否成立判断根在区间［a,c］和［c,b］中的哪一个，

从而得出新区间，仍称为［a,b］。重复运行计算，直至满足精度为止。这就是二分法的计算思想。xk+1Newton法通常预先要给出一个猜测初值xk+1f1(x)

k产生逼近解X*的迭代数列｛xk｝，这就是Newton法的思想。当x0接近x*时收敛很快，但是当x0选择不好时，可能会发散，因此初值的选取很重要。另外，若将该迭代公式改进为__f(x)*k+1 Xkrf-(x)k其中r为要求的方程的根的重数，这就是改进的Newton法，当求解已知重数的方程的根时，在同种条件下其收敛速度要比Newton法快的多。本题中用matlab画曲线y本题中用matlab画曲线y=6x5_45x2+20如下:图4.1y=6*x.A5-45*x.A2+20的曲线由曲线可以看出，该方程有三个实根，根所在区间依次为：[-1,-0.5] [0.5,1] [1.5,2]采用Newton迭代法，依据迭代式：(8-1)(8-2)X(k+1)-X(k)_〃x(k))，k=0,12…(8-1)(8-2)X(K+)—X(KJf'(xck))该方程的迭代式为：6X(k)5一45X(k)2+20X(k+1)=X(k)一 30x(k)4一90x(k)分别选用迭代初值X(0)二一1、x(0)=1、x(0)=2进行迭代，分别求得满足精度的三个实根。(2)matlab源程序：clc;clearx=-2:0.1:2;y=6*x.人5-45*x.人2+20;plot(x,y,'g')%画相应的曲线gridtitle('y=6*x.人5-45*x.人2+20曲线')formatlongroots([600-45020])%根的真实解clearx0=input('请输入迭代初值x0=');formatlong

i=1；e=10人-4;%精度x=x0-(6*x0.人5-45*x0.人2+20)/(30*x0.人4-90*x0);whileabs(x-x0)>ei=i+1;x0=x;x=x0-(6*x0.人5-45*x0.人2+20)/(30*x0.人4-90**0);%迭代enddisplay('方程的根’)display(迭代的次数’)(3)实验结果分析»Q4_2清输;<迭代初值»Q4_2清输;<迭代初值xQ=-L方程的根»Q4_2清输；(迭代初值如=L方程的根>>Q4_2请输；C迭代初值知=2方程的根-0.6545423853973500.-0.6545423853973500.6811741073263151.870799017267090迭代的次数迭代的次数迭代的次数迭代的次数迭代的次数迭代的次数图4.2运行结果对于Newton迭代法，三个初值x0都使得迭代收敛，这是非常重要的。考虑Newton法迭代的收敛性条件：(1)存在一个区间，满足。由曲线和所选的三个区间可知这一条件满足。（2）f（x）是［a,b］上的单调函数，即对一切不变号。经验证所选的三个区间满足这一条件。（3）f（x）的凹向在［a,b］上不变，即在［a,b］上不改变符号。经验证所选的三个区间满足这一条件。（4）-111.5>0>0>0另外，可以看出Newton迭代法收敛速度也很快，且很快达到很高的精度，源于它一般是超线性收敛的。5.编写程序实现大规模方程组的列主元高斯消去法程序，并对所附的方程组进行求解。针对本专业中所碰到的实际问题，提炼一个使用方程组进行求解的例子，并对求解过程进行分析、求解。解：（1）算法原理由于一般线性方程在使用Gauss消去法求解时从求解的过程中可以看到，若ak厂=0厕必须进行行交换，才能使消去过程进行下去。有的时候即使a『）牛0，但是其绝对值非常小，由于机器舍入误差的影响，消去过程也会出现不稳定得现象，导致结果不正确。因此有必要进行列主元技术，以最大可能的消除这种现象。这一技术要寻找行「，使得Ia（I）|二maxa（k-1）ik并将第r行和第k行的元素进行交换以使得当前的akk厂的数值比0要大的多。这种列主元的消去法的主要步骤如下:

.消元过程对k=1,2,…,n-1,进行如下步骤。ik选主元，记Ia1=maxaikrk i>k'1a3很小，这说明方程的系数矩阵严重病态，给出警告，提示结果可能不对。2）交换增广阵A的r,k两行的元素。a-ajj o=k-n+1)3)计算消元jjakakj/,(i=k+1,...,n;j=k+1,……,n+1).回代过程对k=n,n-1,…,1,进行如下计算kjjkjj=k-1jkk至此，完成了整个方程组的求解。(2)matlab源程序:%非压缩,dat51.dat、dat53.datclear;clcfp=fopen('dat53.dat','rb');id=fread(fp,1,'int32');ver=fread(fp,1,'int32');N=fread(fp,1,'int32');q=fread(fp,1,'int32');p=fread(fp,1,'int32');fori=1:NA(i,:)=fread(fp,N,'float');endfori=1:Nd(i)=fread(fp,1,'float');end%正向消去过程fori=1:N-qfork=1:pll=A(i+k,i)/A(i,i);forj=i:i+qA(i+k,j)=A(i+k,j)-ll*A(i,j);endd(i+k)=d(i+k)-ll*d(i);endendfori=N-q+1:Nfork=1:N-ill=A(i+k,i)/A(i,i);forj=i:NA(i+k,j)=A(i+k,j)-ll*A(i,j);endd(i+k)=d(i+k)-ll*d(i);endend%回代过程x(N)=d(N)/A(N,N);fori=N-1:-1:1S=0;ifi+q>Ncv=N;%cv-criticalvalueelsecv=i+q;endforj=i+1:cvS=S+A(i,j)*x(j);endx(i)=(d(i)-S)/A(i,i);endx%压缩,dat52.dat、dat54.datclear;clcfp=fopen('dat54.dat','rb');id=fread(fp,1,'int32');ver=fread(fp,1,'int32');N=fread(fp,1,'int32');q=fread(fp,1,'int32');p=fread(fp,1,'int32');fori=1:NA(i,:)=fread(fp,p+q+1,'float');endfori=1:Nd(i)=fread(fp,1,'float');end%正向消去过程fori=1:Nifi+p<Ncv=p;%cv-criticalvalueelsecv=N-i;endfork=1:cvll=A(i+k,p+1-k)/A(i,p+1);forj=p+1:p+q+1A(i+k,j-k)=A(i+k,j-k)-ll*A(i,j);endd(i+k)=d(i+k)-ll*d(i);endend%回代过程x(N)=d(N)/A(N,p+l);fori=N-l:-l:l;S=0;ii=i+l;jj=P+2；while(ii<=N)&&(jj<=p+q+l)S=S+A(i,jj)*x(ii);11=11+1;♦♦ ♦♦Tendx(i)=(d(i)-S)/A(i,p+l);endx⑶实验结果：

,CcimmendAa.1400 1100 乱1印制 3,L-LOO2.1400 工180 3.110(1 2IQ凶 XL10Q3.1400 ，“忖 3,110(1 1 L*同Cabjina即而throueti20703.1300 1L4D02.WXi2.L400 2.140。 3.3300J.L4D0 皇1⑪0 3.L400W.1400 2.3100 2.liDOJ.1300 1L400Cal'Jina307]throuch2OS53.1400 11400 2.l*Xi 3.L400 3.14C-fl 2.3400 3.14吨 2.14u£i 3.L40Q 第 3.3400 2.140Q 3.1400 3. L4Q0ColiJinz2036thrDUfiti2L0D3.1400 3.1400 2.l*Xi 2.L400 2.1400 3.]JQO 2.l40fl 3L1400 3.L400 2.14<i0 3.3100 2.140Q 3. ]JOO 2. L4D0Cabjin32101through2id33.1440 3. L4DD S.]-d«X' ，LM。 M140。 3.3400 ML40口 2.1400 3.L4DQ 3.1400 ，“口口 3.140口 3. 34QO 3, L4DQ匚乳皿匕2116>±irw曲2L?D3.1400 3. L4DD H.1JCC 3.L-LDO M140。 3.]400 3.L4DH Z.1400 3.L4D0 3.140D 3.]4iD0 3.1400 3. ]4Q0 3. L4D0C■口Ljoms213]>±irD-u£h2L153m"QQ3,L4D03,1-SM 3,1400工3g3,1400 3,1400 3.L400 3,1400 3・“QQ 3.14003门却0 3,L4dOCabjins2146throuih2L603”1切0 3,L4Q03.1.iWILMQ3,140。 1Ng3,L400 3.14W 3,L400 3,1400 1HE3,1400 1：NQO 3,L400>1非压缩矩阵求解结果(部分)州nd日出'L5-21005.2LOO5»2100 5.2L0D 521005.2LOO52100 5.2100 5.2100 5.2M5.2100 5„2100 5.2100 5„2100Colimis$2141匚kruu曲力[前□-Z100 5.2LOO 5,2100 5.2100 52100 5.2100 521Q0 5.2100 5.21W 5,£100 5.2100 a.ZIQO 5.2100 5,2100Colimis4131S6匚hrnu狗4J]705-2100 5.2LOO 5.2100 5.2100 52100 5.2LOO 52100 5.2100 5.2100 5,£100 5.2LOO a-ElOO 5.2100 5.2100Ccdumai1217L匚krnugk门】的a-2100 5.2LOO 5,2100 5.2100 52100 5.2100 521Q0 5.2100 5.21W 5,£100 5.2100 a.ZIQO 5.2100 5,2100Celuixia121£I6匚hrnugh422005-2100 5.2LOO 5.2100 5.2100 52100 5.2LOO 52100 5.2100 5.2100 5,£100 5.2LOO a-ElOO 5.2100 5.2100Ccduma1220L匚krnugk42215a-2100 5.2LOO 5,2100 5.2100 52100 5.2100 521Q0 5.2100 5.21W 5,£100 5.2100 a.ZIQO 5.2100 5,2100roluma132L6uhrnugti4323D5.2100 5.210D 5u210O 5,2LO