西安交通大学计算方法上机作业.docx

1、西安交通大学计算方法上机作业计算方法上机作业 1.对以下和式计算：，要求：(1)若只需保留11个有效数字，该如何进行计算；(2)若要保留30个有效数字，则又将如何进行计算；（1）解题思想和算法实现：根据保留有效位数的要求，可以由公式得出计算精度要求。只需要很少内存，时间复杂度和d呈线性，不需要高浮点支持。先根据while语句求出符合精度要求的n值的大小，然后利用for语句对这n项进行求和，输出计算结果及n值大小即可。（2）matlab源程序：保留11位有效数字时；clearclcformat longn=0;sum=1/(16n)*(4/(8*n+1)-2/(8*n+4)-1/(8*n+5)-

2、1/(8*n+6);while sum=5*10(-11);n=n+1;sum=1/(16n)*(4/(8*n+1)-2/(8*n+4)-1/(8*n+5)-1/(8*n+6);endfor i=0:n-1; sum=sum+1/(16i)*(4/(8*i+1)-2/(8*i+4)-1/(8*i+5)-1/(8*i+6);endvpa(sum,11)n保留30位有效数字时；clearclcformat longn=0;sum=1/(16n)*(4/(8*n+1)-2/(8*n+4)-1/(8*n+5)-1/(8*n+6);while sum=5*10(-30);n=n+1;sum=1/(16n

3、)*(4/(8*n+1)-2/(8*n+4)-1/(8*n+5)-1/(8*n+6);endfor i=0:n-1; sum=sum+1/(16i)*(4/(8*i+1)-2/(8*i+4)-1/(8*i+5)-1/(8*i+6);endvpa(sum,30)n（3）实验结果分析 图1.1 保留11位有效数字的n值及计算结果图 图1.2 保留30位有效数字的n值及计算结果图由计算结果可知，通过合理的误差控制，分别通过7次和22次循环，可以实现题目所要求的精确度。2.某通信公司在一次施工中，需要在水面宽度为20米的河沟底部沿直线走向铺设一条沟底光缆。在铺设光缆之前需要对沟底的地形进行初步探测，从

4、而估计所需光缆的长度，为工程预算提供依据。已探测到一组等分点位置的深度数据(单位:米)如下表所示：分点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.8010.93 (1)请用合适的曲线拟合所测数据点；(2)预测所需光缆长度的近似值，并作出铺设河底光缆的曲线图；（1）解题思想和算法原理给定区间a, b一个分划 ：a=x0x1xN=b 若函数S(x)满足下列条件：1)S(x)在每个区间xi

5、, xj上是不高于3次的多项式。2)S(x)及其2阶导数在a, b上连续。则称S(x)使关于分划的三次样条函数。（2）matlab源程序：clc,clearx=0:1:20;y=9.01 8.96 7.96 7.97 8.02 9.05 10.13 11.18 12.26 13.28 13.32 12.61 11.29 10.22 9.15 7.9 7.95 8.86 9.81 10.80 10.93;n=length(x);l(1)=0;m(n)=0;h=diff(x);df=diff(y)./diff(x);d(1)=0;d(n)=0; for j=2:n-1 l(j)=h(j)/(h(j

6、-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; syms gfor j=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=

7、0;for j=1:n-1 df=diff(s(j),g); warning off all; q=int(sqrt(1+df.2),g,j-1,j); r=r+q;endL=vpa(r,8);disp(the length of the label is L=);disp(L);for j=1:n-1 S(j,:)=sym2poly(s(j);end for j=1:n-1 x1=x(j):0.1:x(j+1); y1=polyval(S(j,:),x1); if j=1 y2=y1; else for i=1:11 k=(j-1)*10+i; y2(k)=y1(i); end endend

8、x2=x(1):0.1:x(n); plot(x,y,o)gridhold onplot(x2,y2,r)（3）实验结果分析图2.1 铺设河底电缆长度图2.2 铺设河底光缆的曲线图由三次样条插值得出的函数曲线的长度和即铺设河底电缆的长度为26.498514。为了提高插值精度，用三次样条插值可以增加插值节点的办法来满足要求，而且在给定节点数的条件下，三次样条插值的精度要优于多项式插值以及线性分段插值，虽然舍弃了降低误差这个优点，但是其曲线的光滑性要好一些。3.假定某天的气温变化记录如下表所示，试用数据拟合的方法找出这一天的气温变化的规律;试计算这一天的平均气温,并试估计误差。时刻01234567

9、89101112平均气温15141414141516182020232528时刻131415161718192021222324平均气温313431292725242220181716（1）解题思想和数学原理：对于具体实验时，通常不是先给出函数的解析式，再进行实验，而是通过实验的观察和测量给出离散的一些点，再来求出具体的函数解析式。又因为测量误差的存在，实际真实的解析式曲线并不一定通过测量给出的所有点。最小二乘法，形成法方程是求解这一问题的很好的方法，本实验运用这一方法实现对给定数据的拟合。对于给定的测量数据(xi,fi)(i=1,2,，n)，设函数分布为特别的，取为多项式 (j=0, 1,，

10、m)则根据最小二乘法原理，可以构造泛函令 (k=0, 1,，m)则可以得到法方程求该解方程组，则可以得到解，因此可得到数据的最小二乘解（2）matlab源程序：x=0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24; %给出题目数据（时间）y=15 14 14 14 14 15 16 18 20 20 23 25 28 31 32 31 29 27 25 24 22 20 18 17 16; %给出题目数据（温度）plot(x,y, m*) %画出各个离散数据点 hold onfor n=2:4; %2、3、4代表拟

11、合函数最高项次数alltemp=25; % alltemp代表数据点总共有25个A=zeros(n+1,n+1); %定义初始正规方程组的系数矩阵AC=ones(n+1,1); %定义初始正规方程组的系数向量CD=zeros(n+1,1); %定义初始正规方程组的向量Dfor i=1:n+1 for j=1:n+1 for k1=1: alltemp A(i,j)=A(i,j)+(x(k1).(i-1+j-1); end end for k2=1: alltemp D(i,1)=D(i,1)+(x(k2).(i-1).*(y(k2); endend %以上为计算出正规方程组矩阵A、D的所有元素

12、的程序tol=1.0e-12;maxit=1000;C=bicg(A,D,tol,maxit); %使用bicg迭代算出正规方程组的系数向量Cp=0; %误差分量E=0; %误差总量if n=2 b=0:24; f=C(1)+C(2).*b+C(3).*(b.2); p=y(b+1)-f; for v=1:25 E=E+(p(v).2; end plot(b,f, r-)end %以上是对2阶拟合函数的图形处理与误差计算if n=3 b=0:24; f=C(1)+C(2).*b+C(3).*(b.2)+C(4).*(b.3); p=y(b+1)-f; for v=1:25 E=E+(p(v).

13、2; end plot(b,f, g-)end %以上是对3阶拟合函数的图形处理与误差计算if n=4 b=0:24; f=C(1)+C(2).*b+C(3).*(b.2)+C(4).*(b.3)+C(5).*(b.4); p=y(b+1)-f; for v=1:25 E=E+(p(v).2; end plot(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向量赋初值for

14、 i=1:n+1 for j=1:n+1 for k=1: alltemp A(i,j)=A(i,j)+(x(k).(i-1+j-1); end end for l=1: alltemp D(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; for v=1:25 E=E+(p(v).2; endplot(b,f, c-) %对指数拟合函数进

15、行图形处理和误差计算b=-C(3);c=C(2)/(2*b);a=exp(b*(c2)+C(1); %算出题设要求的指数拟合函数的各个系数a,b,c,Egrid onlegend(测量数据,二次函数,三次函数,四次函数,指数拟合,Location,NorthWest)hold off %hold on与hold off联合使用，表示将各个曲线画在同一个图中图3.1 二次曲线拟合系数与2范数误差图3.2 三次曲线拟合系数与2范数误差图3.3 四次曲线拟合系数与2范数误差图3.4 指数曲线拟合系数与2范数误差图3.5 数据原始点与拟合曲线对比图（3）实验结果分析：根据国家有关规定，用的是02、08

16、、14、20时四个观测时次的数据做平均，最有代表性。从图中可以看出并不是多项式次数越高越好，随着次数的增高，曲线所呈现出的给定点处和实际的吻合度越好，但对于其他地方的吻合度降低了。4.设计算法，求出非线性方程的所有根，并使误差不超过。解：（1）解题思想和算法实现：对于一个非线性方程的数值解法很多。在此介绍两种最常见的方法：二分法和Newton法。首先要研究函数的形态，确定根的数量和大致区间的位置。对于二分法，其数学实质就是说对于给定的待求解的方程f(x)，其在a,b上连续，f(a)f(b)0，且f(x)在a,b内仅有一个实根x*，取区间中点c，若，则c恰为其根，否则根据f(a)f(c)e i=

17、i+1; x0=x; x=x0-(6*x0.5-45*x0.2+20)/(30*x0.4-90*x0);%迭代enddisplay(方程的根)xdisplay(迭代的次数)i（3）实验结果分析：图4.2 运行结果对于Newton迭代法，三个初值x0都使得迭代收敛，这是非常重要的。考虑Newton法迭代的收敛性条件：（1）存在一个区间，满足。由曲线和所选的三个区间可知这一条件满足。（2）f(x)是a,b上的单调函数，即对一切不变号。经验证所选的三个区间满足这一条件。（3）f(x)的凹向在a,b上不变，即在a,b上不改变符号。经验证所选的三个区间满足这一条件。（4）-111.5000 另外，可以看

18、出Newton迭代法收敛速度也很快，且很快达到很高的精度，源于它一般是超线性收敛的。5.编写程序实现大规模方程组的列主元高斯消去法程序，并对所附的方程组进行求解。针对本专业中所碰到的实际问题，提炼一个使用方程组进行求解的例子，并对求解过程进行分析、求解。解:（1）算法原理由于一般线性方程在使用Gauss消去法求解时，从求解的过程中可以看到，若=0，则必须进行行交换，才能使消去过程进行下去。有的时候即使0，但是其绝对值非常小，由于机器舍入误差的影响，消去过程也会出现不稳定得现象，导致结果不正确。因此有必要进行列主元技术，以最大可能的消除这种现象。这一技术要寻找行r，使得并将第r行和第k行的元素进

19、行交换，以使得当前的的数值比0要大的多。这种列主元的消去法的主要步骤如下：1.消元过程对k=1,2,n-1,进行如下步骤。1)选主元，记很小，这说明方程的系数矩阵严重病态，给出警告，提示结果可能不对。2)交换增广阵A的r，k两行的元素。 (j=k,n+1) 3）计算消元(i=k+1,n; j=k+1,n+1)2.回代过程对k= n, n-1,1,进行如下计算至此，完成了整个方程组的求解。（2）matlab源程序：%非压缩,dat51.dat、dat53.datclear;clcfp=fopen(dat53.dat,rb);id=fread(fp,1,int32);ver=fread(fp,1,

20、int32);N=fread(fp,1,int32);q=fread(fp,1,int32);p=fread(fp,1,int32);for i=1:N A(i,:)=fread(fp,N,float);endfor i=1:N d(i)=fread(fp,1,float);end%正向消去过程for i=1:N-q for k=1:p ll=A(i+k,i)/A(i,i); for j=i:i+q A(i+k,j)=A(i+k,j)-ll*A(i,j); end d(i+k)=d(i+k)-ll*d(i); endendfor i=N-q+1:N for k=1:N-i ll=A(i+k,i

21、)/A(i,i); for j=i:N A(i+k,j)=A(i+k,j)-ll*A(i,j); end d(i+k)=d(i+k)-ll*d(i); endend%回代过程x(N)=d(N)/A(N,N);for i=N-1:-1:1 S=0; if i+qN cv=N;%cv-critical value else cv=i+q; end for j=i+1:cv S=S+A(i,j)*x(j); end x(i)=(d(i)-S)/A(i,i);endx%压缩,dat52.dat、dat54.datclear;clcfp=fopen(dat54.dat,rb);id=fread(fp,1

22、,int32);ver=fread(fp,1,int32);N=fread(fp,1,int32);q=fread(fp,1,int32);p=fread(fp,1,int32);for i=1:N A(i,:)=fread(fp,p+q+1,float);endfor i=1:N d(i)=fread(fp,1,float);end%正向消去过程for i=1:N if i+pN cv=p;%cv-critical value else cv=N-i; end for k=1:cv ll=A(i+k,p+1-k)/A(i,p+1); for j=p+1:p+q+1 A(i+k,j-k)=A(

23、i+k,j-k)-ll*A(i,j); end d(i+k)=d(i+k)-ll*d(i); endend%回代过程x(N)=d(N)/A(N,p+1);for i=N-1:-1:1; S=0; ii=i+1; jj=p+2; while (ii=N)&(jj=p+q+1) S=S+A(i,jj)*x(ii); ii=ii+1; jj=jj+1; end x(i)=(d(i)-S)/A(i,p+1);end x（3）实验结果：非压缩矩阵求解结果（部分）压缩矩阵求解结果（部分）（4）分析心得：采用Gauss消去法时，如果在消元时对角线上的元素始终较大，那么本方法不需要进行列主元计算，计算结果一般

24、就可以达到要求，否则必须进行列主元这一步，以减少机器误差带来的影响，使方法得出的结果正确。（5）实例化学反应方程式严格遵守质量守恒定律，书写化学反应方程式写出反应物和生成物后，往往左右两边各原子数目不相等，不满足质量守恒定律，这就需要通过计算配平来解决。对于化学反应方程式X1KMnO4+x2MnSO4+x3H2O=x4MnO2+x5K2SO4+x6H2SO4 ，可按以下方式配平：上述化学反应式中包含5种不同的原子（钾、锰、氧、硫、氢），按照原子守恒有:K: x1=2x5Mn: x1+x2=x4O: 4x1+4x2+x3=2x4+4x5+4x6S: x2=x5+x6H: 3x2=2x6上述方程组

25、中有5个方程，6个未知量，因此有无数解。可先令x6=1，于是形成方程组。使用Gauss消去法求解此方程组得:x =1.0000,1.5000,1.0000,2.5000,0.5000,1.0000；由于化学方程式通常取最简的正整数，因此将x乘2得配平的化学方程式：2KMnO4+3MnSO4+2H2O=5MnO2+K2SO4+2H2SO4 程序如下：clear;clc;a=1 0 0 0 -2 0; 1 1 0 -1 0 0; 4 4 1 -2 -4 -4; 0 1 0 0 -1 -1; 0 0 2 0 0 -2; 0 0 0 0 0 -1;b=0 0 0 0 0 -1;n=length(b);for k=1:n-1 if a(k,k)=0 disp(error); return; end for i=k+1:n a(i,k)=a(i,k)/a(k,k); for j=k+1:n a(i,j)=a(i,j)-a(i,k)*a(k,j); end b(i)=b(i)-a(i,k)*b(k); endendx(n)=b(n)/a(n,n);for k=n-1:-1:1 S=b(k); for j=k+1:n S=S-a(k,j)*x(j); end x(k)=S/a(k,k);endx