计算方法b上机报告.docx

《计算方法b上机报告.docx》由会员分享，可在线阅读，更多相关《计算方法b上机报告.docx（19页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

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

2、=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深度分点78910111213深度分点14151617181920深度 (1)请用适宜的曲线拟合所测数据点；(2)预测所需光缆

4、长度的近似值，并作出铺设河底光缆的曲线图；问题分析：此题的主要目的是对测量数据进行拟合，同时对拟合曲线进行线积分即可得到河底光缆长度的近似值。由于数值点较多时，使用拉格朗日差值多项式会出现龙格现象。为了将所有的数据点都用上，采用分段差插法，此题使用三次样条插值。算法思想：样条函数在每个子区间上是三次多项式，它的二阶导数必是一次多项式。假设用 记在 处的二阶导数。那么在区间 上 式中 1对上式进行两次积分得 2它的一阶导数为 3满足连续性条件，即 上式和3式得 4用差商记号，并记 4式可以写成 方程组可以写成如下形式 自然样条插值条件为在估计河底光缆长度时使用第一类线积分程序实现：clear;c

5、lc;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%计算牛顿二阶差商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;

6、%追赶法求解带状矩阵的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)*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=linsp

7、ace(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-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(

8、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平均气温15141414141516182020232528时刻131415161718192021222324平均气温313431292725242220181716问题分析：对一天的气温进行数据拟合，可以考虑使用最小二乘的二次函数、三次函数、四次函数以及指数函数。问题的难度是求解各种拟合函数的系数。利用法方程求解最小二乘系数时，方程的解不够稳定，此题利用正交化方法。程序实现：clear;clc

10、;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(请输入函数的次数 );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

11、 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来存放 for 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-

12、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); %ai存放各级系数 f=0; end 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; %平均温

13、度disp(平均温度为,num2str(t),) 运行结果：二次函数时，系数a1=8.3063,a2=2.6064,a3=-。误差为，平均温度为。函数形式为 三次函数时，系数,误差为，平均温度为。函数形式 四次函数时，系数a1=16.7939,a2=-3.7050,a3=0.8909,a4=-0.0532,a5=0.0009,误差为，平均温度为。函数形式为 当为指数函数时，假定指数函数的形式为 。只需将y值求对数，其主体局部不变，求得的系数为a1=2.3835,a2=0.1251,a3=-0.0044,误差为，平均温度为。函数形式为 四种函数图像的比照：红色表示指数函数，绿色表示二次函数，青色

14、表示三次函数，紫色表示四次函数，蓝色的折线表示数据点。总结心得：同上面的结果可以看出，随着次数的增加误差在逐渐减小，拟合曲线也更加接近数据点。因此，在使用最小二乘进行数据拟合时应尽量使用较高次数的拟合函数。通过比照图像，可以发现二次函数和指数型二次函数图像比拟接近，那是因为指数型二次函数是对数据点求对数进行拟合的，因此二者比拟接近。给定的数据点在中间的位置时波动较大，造成了拟合曲线的误差也较大，但并没有影响到整体的规律。在进行最小二次编程中，使用正交化方法中间变量较多容易造成混乱，计算量也较大，但是比使用法方程方法更稳定。此题在编程过程中有一定的适用性，只需要对数据点进行更改就可以拟合出给定次

15、数的多项式，为以后的数据拟合提供了方便。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);grid;k=1;y(202)=0;%使用二分法来确定迭代初值for i=1:201 if y(i)*y(i+1)0

16、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 cabs(A(max,i) max=k; end end temp=A(i,i:m);A(i,i:m)=A(max,i:m);A(max,i:m)=temp; for j=(i+1):m A(j,:)=A(j,:)-A(i,:)*A(j,i)/A(i,i); endenddisp(回代求解) x(m)=A(m,n)/A(m,m); for k=(m-1):-1:1 x(k)=(A(k,n)-A(k,k+1:m)*x(k+1:m)/A(k,k); end x运算结果为，。运算结果和正确结果吻合。