计算方法B上机报告

1、计算方法B上机实习报告 计算方法B上机实习报告1.对以下和式计算：，要求：(1)若只需保留11个有效数字，该如何进行计算；(2)若要保留30个有效数字，则又将如何进行计算；问题分析： 在该题中S的每一项存在两个相近的数相减的问题，因此为了避免有效数字损失，最好是改变运算顺序，分别将正数和负数相加，然后再将其和相加。另外，sn中有多个负数相加，可以按照绝对值递增的顺序求和，以减少舍入误差的影响。同时，为了避免大数吃小数的问题，本题先计算出保留目标有效数字所需要的迭代次数，然后采用倒序相加的方法实现。程序实现：clear;clc;m=input('请输入要保留的有效数字位数：');

2、s1=0;s2=0;k=0;s=1;%判断多需要的迭代次数while s>=0.5*10-(m-1) s=4/(16k*(8*k+1)-(2/(16k*(8*k+4)+1/(16k*(8*k+5)+1/(16k*(8*k+6); k=k+1;end%正负数分别按照绝对值递增的顺序倒序相加for n=(k-1):-1:0 a1=4/(16n*(8*n+1); a2=2/(16n*(8*n+4); a3=1/(16n*(8*n+5); a4=1/(16n*(8*n+6); s1=a1+s1; s2=a4+a3+a2+s2;end S=s1-s2; S=vpa(S,m)运算结果：总结心得：在计

3、算求和问题中，应特别注意相近数相减的问题，这样会造成有效数字灾难性的损失。另外在两个数量级相差较大的数字相加减时，较小数的有效数字会被丧失，因此要按照从小到大的顺序相加。在上题计算中分别对正负相采用倒序相加，这样就有效的避免了“大数吃小数”的问题。2.某通信公司在一次施工中，需要在水面宽度为20米的河沟底部沿直线走向铺设一条沟底光缆。在铺设光缆之前需要对沟底的地形进行初步探测，从而估计所需光缆的长度，为工程预算提供依据。已探测到一组等分点位置的深度数据(单位:米)如下表所示：分点0123456深度9.018.967.967.978.029.0510.13分点78910111213深度11.18

4、12.2613.2813.3212.6111.2910.22分点14151617181920深度9.157.907.958.869.8110.8010.93 (1)请用合适的曲线拟合所测数据点；(2)预测所需光缆长度的近似值，并作出铺设河底光缆的曲线图；问题分析：本题的主要目的是对测量数据进行拟合，同时对拟合曲线进行线积分即可得到河底光缆长度的近似值。由于数值点较多时，使用拉格朗日差值多项式会出现龙格现象。为了将所有的数据点都用上，采用分段差插法，本题使用三次样条插值。算法思想：样条函数在每个子区间上是三次多项式，它的二阶导数必是一次多项式。若用 记在 处的二阶导数。则在区间 上 式中 （1）

5、对上式进行两次积分得 （2）它的一阶导数为 （3）满足连续性条件，即 上式和（3）式得 （4）用差商记号，并记 （4）式可以写成 方程组可以写成如下形式 自然样条插值条件为在估计河底光缆长度时使用第一类线积分程序实现：clear;clc;x=0: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.90 7.95 8.86 9.81 10.80 10.93;d=y;plot(x,y,'k.','markersize',15)hold on

6、%计算牛顿二阶差商for k=1:2 for i=21:-1:(k+1) d(i)=(d(i)-d(i-1)/(x(i)-x(i-k); endend%假定d的边界条件，采用自然三次样条for i=2:20 d(i)=6*d(i+1);endd(1)=0;d(21)=0;%追赶法求解带状矩阵的m值a=0.5*ones(1,21);b=2*ones(1,21);c=0.5*ones(1,21);a(1)=0;c(21)=0;u=ones(1,21);u(1)=b(1);r=c;yy(1)=d(1);%追的过程for k=2:21 l(k)=a(k)/u(k-1); u(k)=b(k)-l(k)*

7、r(k-1); yy(k)=d(k)-l(k)*yy(k-1);end%赶的过程m(21)=yy(21)/u(21);for k=20:-1:1 m(k)=(yy(k)-r(k)*m(k+1)/u(k);end%利用插值点画出拟合曲线k=1;nn=100;xx=linspace(0,20,nn);l=0;for j=1:nn for i=2:20 if xx(j)<=x(i) k=i; break; else k=i+1; end endh=1;xbar=x(k)-xx(j);xmao=xx(j)-x(k-1);s(j)=(m(k-1)*xbar3/6+m(k)*xmao3/6+(y(k

8、-1)-m(k-1)*h2/6)*xbar+(y(k)-m(k)*h2/6)*xmao)/h;sp(j)=-m(k-1)*(x(k)-xx(j)2/(2*h)+m(k)*(xx(j)-x(k-1)2/(2*h)+(y(k)-y(k-1)/h-(m(k)-m(k-1)*h/6;l(j+1)=(1+sp(j)2)0.5*(20/nn)+l(j);%利用第一类线积分求河底光缆的长度 end%绘图plot(xx,s,'r-','linewidth',1.5)griddisp('所需光缆长度为',num2str(l(nn+1),'米')运行

9、结果：总结心得：采用三次样条插值对数据进行拟合时，可以有效避免龙格现象。在本题的计算中采用自然三次样条函数的边界条件。在解线性方程组时使用了追赶法求解带状矩阵，在求解三对角矩阵时追赶法计算速度快，是一种求解线性方程组的有效手段。在估计河底光缆长度时使用第一类线积分。本题计算中间变量非常多，在调试的过程中遇到了一些麻烦，这更加使我认识到在编程的过程中由不得一点儿马虎，在下面问题的处理中我变得更加小心。3.假定某天的气温变化记录如下表所示，试用数据拟合的方法找出这一天的气温变化的规律;试计算这一天的平均气温,并试估计误差。时刻0123456789101112平均气温1514141414151618

10、2020232528时刻131415161718192021222324平均气温313431292725242220181716问题分析：对一天的气温进行数据拟合，可以考虑使用最小二乘的二次函数、三次函数、四次函数以及指数函数。问题的难度是求解各种拟合函数的系数。利用法方程求解最小二乘系数时，方程的解不够稳定，本题利用正交化方法。程序实现：clear;clc;x=0:24;y=15 14 14 14 14 15 16 18 20 20 23 25 28 31 34 31 29 27 25 24 22 20 18 17 16;m=length(x);n=input('请输入函数的次数 &

11、#39;);plot(x,y,'k.',x,y,'-')grid;hold on;n=n+1;G=zeros(m,n+1);G(:,n+1)=y'c=zeros(1,n);%建立c来存放q=0;f=0;b=zeros(1,m);%建立b用来存放%形成矩阵Gfor j=1:n for i=1:m G(i,j)=x(1,i)(j-1); endend%建立矩阵Qkfor k=1:n for i=k:m c(k)=G(i,k)2+c(k); end c(k)=-sign(G(k,k)*(c(k)0.5); w(k)=G(k,k)-c(k);%建立w来存放 fo

12、r j=k+1:m w(j)=G(j,k); end b(k)=c(k)*w(k); %变换矩阵Gk-1到Gk G(k,k)=c(k); for j=k+1:n+1 q=0; for i=k:m q=w(i)*G(i,j)+q; end s=q/b(k); for i=k:m G(i,j)=s*w(i)+G(i,j); end endend%求解三角方程 Rx=h1a(n)=G(n,n+1)/G(n,n); for i=n-1:(-1):1 for j=i+1:n f=G(i,j)*a(j)+f; end a(i)=(G(i,n+1)-f)/G(i,i); %a（i）存放各级系数 f=0; e

13、nd a%拟合后的回代过程 p=zeros(1,m); for j=1:m for i=1:n p(j)=p(j)+a(i)*x(j)(i-1); end end plot(x,p,'r*',x,p,'-'); E2=0;%用E2来存放误差%误差求解for i=n+1:m E2=G(i,n+1)2+E2;endE2=E20.5;disp('误差为');disp(E2);t=0for i=1:m t=t+p(i);endt=t/m; %平均温度disp('平均温度为',num2str(t),'') 运行结果：二次函数

14、时，系数a1=8.3063,a2=2.6064,a3=-0.0938。误差为16.7433，平均温度为21.2。函数形式为 三次函数时，系数a1=13.3880,a2=-0.2273,a3=0.2075,a4=-0.0084,误差为11.4482，平均温度为21.2。函数形式 四次函数时，系数a1=16.7939,a2=-3.7050,a3=0.8909,a4=-0.0532,a5=0.0009,误差为7.6838，平均温度为21.2。函数形式为 当为指数函数时，假定指数函数的形式为 。只需将y值求对数，其主体部分不变，求得的系数为a1=2.3835,a2=0.1251,a3=-0.0044,

15、误差为14.6867，平均温度为20.9399。函数形式为 四种函数图像的对比：红色表示指数函数，绿色表示二次函数，青色表示三次函数，紫色表示四次函数，蓝色的折线表示数据点。总结心得：同上面的结果可以看出，随着次数的增加误差在逐渐减小，拟合曲线也更加接近数据点。因此，在使用最小二乘进行数据拟合时应尽量使用较高次数的拟合函数。通过对比图像，可以发现二次函数和指数型二次函数图像比较接近，那是因为指数型二次函数是对数据点求对数进行拟合的，因此二者比较接近。给定的数据点在中间的位置时波动较大，造成了拟合曲线的误差也较大，但并没有影响到整体的规律。在进行最小二次编程中，使用正交化方法中间变量较多容易造成

16、混乱，计算量也较大，但是比使用法方程方法更稳定。本题在编程过程中有一定的适用性，只需要对数据点进行更改就可以拟合出给定次数的多项式，为以后的数据拟合提供了方便。4.设计算法，求出非线性方程的所有根，并使误差不超过。问题分析：本题采用牛顿迭代法计算方程根。可令 ，其迭代格式为要使牛顿迭代法收敛，则必须寻找根的合理区间a,b,使得在该区间内,即有根。在选定的区间内函数的一二阶导数、均不变号。选定一个初值，使 则牛顿迭代法收敛。程序实现：clear;clc;%观察根的大致位置a=-10:0.1:10;for i=1:201y(i)=6*a(i)5-45*a(i)2+20;endplot(a,y);g

17、rid;k=1;y(202)=0;%使用二分法来确定迭代初值for i=1:201 if y(i)*y(i+1)<0 b(k)=a(i); k=k+1; continue endend x=b; n=length(x);%牛顿迭代法计算方程根dlt=1e-4;for i=1:n for j=1:100 X(i,j)=x(i)-(6*x(i)5-45*x(i)2+20)/(30*x(i)4-90*x(i); c=abs(X(i,j)-x(i); if c<dlt%误差判断 break end x(i)=X(i,j); endenddisp('计算结果为:');x运行结

18、果：函数图像如下，可以观察根的大致位置。计算结果：总结心得： 在该题的计算中，首先画出了函数的图像可以初步观察出根的大致位置，有效缩减了求根区间。为防止迭代进入死循环，还设置了最高迭代次数为100次，结果显示方程的三个根均在4次以内就可以迭代出，可以看出牛顿迭代法的求解速度相当快。5.编写程序实现大规模方程组的列主元高斯消去法程序，并对所附的方程组进行求解。针对本专业中所碰到的实际问题，提炼一个使用方程组进行求解的例子，并对求解过程进行分析、求解。问题分析：高斯消去法在整个消去过程中， 均不为零，但是若在某一步中，对角线上元的绝对值比其下方的某些元素的绝对值小得多，即那么，相应的数 将会是绝对

19、值很大的数，这可能引起上溢而中断计算，即使不发生上溢，也会带来很大的舍入误差。而列主元消去法的原则使中，绝对值最大的一个换到（k,k）位置，这样就有效的避免了大数除以小数而造成的上溢现象。列主元消去法的主要核心是查找最大的并交换到。程序实现：（1）非压缩格式dat141.dat;dat142.datclc;clear;%读取数据fid=fopen('dat141.dat','r');a=fread(fid,5,'long',0);n=a(3); a=fread(fid,n,n+1,'float'); for i=1:n for j

20、=1:n A(i,j)=a(j,i); end b(i)=a(i,n+1); end%查找主元并交换位置for k=1:n for m=k:n if abs(A(m,k)=max(abs(A(:,k) continue else r=A(k,:);A(k,:)=A(m,:);A(m,:)=r; r=b(k);b(k)=b(m);b(m)=r; continue end end for i=k+1:n m=A(i,k)/A(k,k); for j=1:n A(i,j)=A(i,j)-m*A(k,j); end b(i)=b(i)-m*b(k); endend%回代过程x(n)=b(n)/A(n,

21、n);b(n)=x(n); for k=n-1:(-1):1 for j=k+1:n b(k)=b(k)-A(k,j)*x(j); end x(k)=b(k)/A(k,k); b(k)=x(k); end x=x'（2）压缩格式dat143.datclear;clc;format short%读取数据fid=fopen('dat143.dat','r');d=fread(fid,5,'long',0);n=d(3);%矩阵维数p=d(5);%上带宽q=d(4);%下带宽a=fread(fid,p+q+1,n,'float'

22、);b=fread(fid,n,1,'float');a=a'd=b'%消去过程for k=1:n-1 for i=k+1:min(k+p,n) a(i,k-i+p+1)=a(i,k-i+p+1)/a(k,p+1); for j=k+1:min(k+q,n) a(i,j-i+p+1)=a(i,j-i+p+1)-a(i,k-i+p+1)*a(k,j-k+p+1); end b(i) = b(i) - b(k) * a(i,k-i+p+1); endend%回代过程x(n) = b(n)/a(n,p+1);for k=n-1:-1:1 sum=0; for j=k+1:min(k+q,n) sum=sum+a(k,j-k+p+1)*x(j); end x(k)=(b(k)-sum)/a(k,p+1);endx运算结果：1）非压缩格式dat141.dat;dat142.dat的运算结果如下（2）压缩格式dat143.dat总结心得：所给数据的方程的根都相等，其中dat141.dat，为