西安交通大学计算方法B上机报告

计算方法上机报告姓名：学号：班级：能动上课班级：题目及求解：一、对以下和式计算：∑∞⎪⎭⎫⎝⎛+-+-+-+=0681581482184161n n n n S n ，要求：① 若只需保留11个有效数字，该如何进行计算； ② 若要保留30个有效数字，则又将如何进行计算；1 算法思想（1）根据精度要求估计所加的项数，可以使用后验误差估计，通项为： 142111416818485861681n n na n n n n n ε⎛⎫=---<< ⎪+++++⎝⎭； （2）为了保证计算结果的准确性，写程序时，从后向前计算； （3）使用Matlab 时，可以使用以下函数控制位数： digits(位数)或vpa(变量，精度为数)2 算法结构;0=s⎪⎭⎫⎝⎛+-+-+-+=681581482184161n n n n t n ； for 0,1,2,,n i =⋅⋅⋅ if 10m t -≤end;for ,1,2,,0n i i i =--⋅⋅⋅;s s t =+3 Matlab 源程序clear; %清除工作空间变量 clc; %清除命令窗口命令 m=input('请输入有效数字的位数m='); %输入有效数字的位数 s=0;for n=0:50t=(1/16^n)*(4/(8*n+1)-2/(8*n+4)-1/(8*n+5)-1/(8*n+6));if t<=10^(-m) %判断通项与精度的关系break;endend;fprintf('需要将n值加到n=%d\n',n-1); %需要将n值加到的数值for i=n-1:-1:0t=(1/16^i)*(4/(8*i+1)-2/(8*i+4)-1/(8*i+5)-1/(8*i+6));s=s+t; %求和运算ends=vpa(s,m) %控制s的精度4 结果与分析若保留11位有效数字，则n=7，此时求解得：s =3.1415926536；若保留30位有效数字时，则n=22, 此时求解得：s =3.8。

通过上面的实验结果可以看出，通过从后往前计算，这种算法很好的保证了计算结果要求保留的准确数字位数的要求。

二、某通信公司在一次施工中，需要在水面宽度为20米的河沟底部沿直线走向铺设一条沟底光缆。

在铺设光缆之前需要对沟底的地形进行初步探测，从而估计所需光缆的长度，为工程预算提供依据。

已探测到一组等分点位置的深度数据(单位:米)如下表所示：分点0 1 2 3 4 5 6深度9.01 8.96 7.96 7.97 8.02 9.05 10.13分点7 8 9 10 11 12 13深度11.18 12.26 13.28 13.32 12.61 11.29 10.22分点14 15 16 17 18 19 20深度9.15 7.90 7.95 8.86 9.81 10.80 10.93①请用合适的曲线拟合所测数据点；②预测所需光缆长度的近似值，作出铺设河底光缆的曲线图；1 算法思想如果使用多项式差值，则由于龙格现象，误差较大，因此，用相对较少的插值数据点作插值，可以避免大的误差，但是如果又希望将所得数据点都用上，且所用数据点越多越好，可以采用分段插值方式，即用分段多项式代替单个多项式作插值。

分段多项式是由一些在相互连接的区间上的不同多项式连接而成的一条连续曲线，其中三次样条插值方法是一种具有较好“光滑性”的分段插值方法。

在本题中，假设所铺设的光缆足够柔软，在铺设过程中光缆触地走势光滑，紧贴地面，并且忽略水流对光缆的冲击。

海底光缆线的长度预测模型如图2-1所示，光缆从Ah。

点铺至B点，在某点处的深度为i图2-1 海底光缆线的长度预测模型计算光缆长度时，用如下公式：20()L f x ds =⎰200(f x =⎰ 191(k kk f x +==∑⎰=2 算法结构1） For n i ,,2,1,0⋅⋅⋅=1.1 i i M y ⇒ 2） For 2,1=k2.1 For k n n i ,,1, -=2.1.1 i k i i i i M x x M M ⇒----)/()(13）101h x x ⇒-4）For 1-,,2,1n i =4.1 11++⇒-i i i h x x4.2 b a c c h h h i i i i i i ⇒⇒-⇒+++2;1;)/(11 4.3 i i d M ⇒+165）0000;;c M d M d n n ⇒⇒⇒λn n n b a b ⇒⇒⇒2;;20μ6）1111,γμ⇒⇒d b7）获取M 的矩阵元素个数，存入m 8）For m k ,,3,2 =8.1 k k k l a ⇒-1/μ8.2 k k k k c l b μ⇒⋅-1- 8.3 k k k k l d γγ⇒⋅-1- 9）m m m M ⇒μγ/10）For 1,,2,1 --=m m k10.1 k k k k k M M c ⇒⋅-+μγ/)(1 11）获取x 的元素个数存入s 12）k ⇒113） For 1,,2,1-=s i13.1 if i x x ≤~then k i ⇒；break else k i ⇒+114）xx x x x x h x x k k k k ˆ~;~;11⇒-⇒-⇒--- y h x h M y x h M y x M x M k k k k k k ~/]ˆ)6()6(6ˆ6[2211331⇒-+-++---3 Matlab 源程序clear; %清除工作空间变量 clc; %清除命令窗口命令x=0:1:20; %产生从0到20含21个等分点的数组 X=0:0.2: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]; %等分点位置的深度数据 n=length(x); %等分点的数目 N=length(X);%% 求三次样条插值函数s(x) M=y;for k=2:3; %计算二阶差商并存放在M 中 for i=n:-1:k;M(i)=(M(i)-M(i-1))/(x(i)-x(i-k+1)); end endh(1)=x(2)-x(1); %计算三对角阵系数a,b,c及右端向量d for i=2:n-1;h(i)=x(i+1)-x(i);c(i)=h(i)/(h(i)+h(i-1));a(i)=1-c(i);b(i)=2;d(i)=6*M(i+1);endM(1)=0; %选择自然边界条件M(n)=0;b(1)=2;b(n)=2;c(1)=0;a(n)=0;d(1)=0;d(n)=0;u(1)=b(1); %对三对角阵进行LU分解y1(1)=d(1);for k=2:n;l(k)=a(k)/u(k-1);u(k)=b(k)-l(k)*c(k-1);y1(k)=d(k)-l(k)*y1(k-1);endM(n)=y1(n)/u(n); %追赶法求解样条参数M(i)for k=n-1:-1:1;M(k)=(y1(k)-c(k)*M(k+1))/u(k);ends=zeros(1,N);for m=1:N;k=1;for i=2:n-1if X(m)<=x(i);k=i-1;break;elsek=i;endendH=x(k+1)-x(k); %在各区间用三次样条插值函数计算X点处的值x1=x(k+1)-X(m);x2=X(m)-x(k);s(m)=(M(k)*(x1^3)/6+M(k+1)*(x2^3)/6+(y(k)-(M(k)*(H^2)/6))*x1+(y(k+1)-(M(k+1)*(H^2)/6 ))*x2)/H;end%% 计算所需光缆长度L=0; %计算所需光缆长度for i=2:NL=L+sqrt((X(i)-X(i-1))^2+(s(i)-s(i-1))^2);enddisp('所需光缆长度为L=');disp(L);figureplot(x,y,'*',X,s,'-') %绘制铺设河底光缆的曲线图xlabel('位置','fontsize',16); %标注坐标轴含义ylabel('深度/m','fontsize',16);title('铺设河底光缆的曲线图','fontsize',16);grid;4 结果与分析铺设海底光缆的曲线图如图2-2所示：图2-2 铺设河底光缆曲线图仿真结果表明，运用分段三次样条插值所得的拟合曲线能较准确地反映铺设光缆的走势图，计算出所需光缆的长度为 L=26.4844m 。

三、 假定某天的气温变化记录如下表所示，试用数据拟合的方法找出这一天的气温变化的规律；试计算这一天的平均气温，并试估计误差。

1 算法思想在本题中，数据点的数目较多。

当数据点的数目很多时，用“多项式插值”方法做数据近似要用较高次的多项式，这不仅给计算带来困难，更主要的缺点是误差很大。

用“插值样条函数”做数据近似，虽然有很好的数值性质，且计算量也不大，但存放参数Mi 的量很大，且没有一个统一的数学公式来表示，也带来了一些不便。

另一方面，在有的实际问题中，用插值方法并不合适。

当数据点的数目很大时，要求)(x p 通过所有数据点，可能会失去原数据所表示的规律。

如果数据点是由测量而来的，必然带有误差，插值法要求准确通过这些不准确的数据点是不合适的。

在这种情况下，不用插值标准而用其他近似标准更加合理。

通常情况下，是选取i α使2E 最小，这就是最小二乘近似问题。

在本题中，采用“最小二乘法”找出这一天的气温变化的规律，使用二次函数、三次函数、四次函数以及指数型函数2)(c t b ae C --=，计算相应的系数，估算误差，并作图比较各种函数之间的区别。

2 算法结构本算法用正交化方法求数据的最小二乘近似。

假定数据以用来生成了G ，并将y 作为其最后一列（第1+n 列）存放。

结果在a 中，ρ是误差22E 。

I 、使用二次函数、三次函数、四次函数拟合时 1. 将“时刻值”存入i x ，数据点的个数存入m2. 输入拟合多项式函数)(x p 的最高项次数1-=n h ，则拟合多项式函数为)()()()(2211x g x g x g x p n n ααα+⋅⋅⋅++= , 根据给定数据点确定G For1,,2,1,0-⋅⋅⋅=n jFor m i ,,2,1⋅⋅⋅=2.1 1,+⇒j i ji g x 2.2 1,+⇒n i i g y 3. For n k ,,2,1⋅⋅⋅=3.1 [形成矩阵k Q ]3.1.1 σ⇒-∑=mk i ikkk g g 2/12))(sgn( 3.1.2 k kk g ωσ⇒-3.1.3 For m k k j ,,2,1⋅⋅⋅++=3.1.3.1 j jk g ω⇒ 3.1.4 βσω⇒k 3.2 [变换1-k G 到k G ] 3.2.1 kk g ⇒σFor 1,,,2,1+⋅⋅⋅++=n n k k j 3.2.2 t g mk i ij i ⇒∑=βω/)(3.2.3 For m k k i ,,1,⋅⋅⋅+=3.2.3.1 ij i ij g t g ⇒+ω4. [解三角方程1h Ra =] 4.1nnn n n a g g ⇒+/1.4.2 For 1,,2,1⋅⋅⋅--=n n i 4.2.1iii ni j j ijn i a g x gg ⇒-∑+=+/][11.5. [计算误差22E ]ρ⇒∑+=+mn i n i g121,II 、使用指数函数拟合时现将指数函数进行变形： 将y C =,x t=代入2)(c t b aeC --=得:2)(c x b aey --=对上式左右取对数得: 222ln ln bx bcx bc a y -+-=令b bc bc a y z -==-==21202ln ln ααα，，，则可得多项式: 2210x x z ααα++=（3）Matlab 源程序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]=size(x); %将数据点的个数存入m T=sum(y(1:m))/m;fprintf('一天的平均气温为T=%f\n',T); %求一天的平均气温 %% 二次、三次、四次函数的最小二乘近似h=input('请输入拟合多项式的最高项次数='); %根据给定数据点生成矩阵G n=h+1; G=[];for j=0:(n-1)g=x.^j; %g(x)按列排列G=vertcat(G,g); %g垂直连接GendG=G'; %转置得到矩阵Gfor i=1:m %将数据y作为G的最后一列（n+1列）G(i,n+1)=y(i);endG;for k=1:n %形成矩阵Q(k)if G(k,k)>0;sgn=1;elseif G(k,k)==0;sgn=0;else sgn=-1;endsgm=-sgn*sqrt(sum(G(k:m,k).^2));w=zeros(1,n);w(k)=G(k,k)-sgm;for j=k+1:mw(j)=G(j,k);endbt=sgm*w(k);G(k,k)=sgm; %变换Gk-1到Gkfor j=k+1:n+1t=sum(w(k:m)*G(k:m,j))/bt;for i=k:m;G(i,j)=G(i,j)+t*w(i);endendendA (n)=G(n,n+1)/G(n,n); %解三角方程求系数Afor i=n-1:-1:1A (i)=(G(i,n+1)-sum(G(i,i+1:n).*A (i+1:n)))/G(i,i);ende=sum(G(n+1:m,n+1).^2); %计算误差efprintf('%d次函数的系数是：',h); %输出系数a及误差edisp(A);fprintf('使用%d次函数拟合的误差是%f：',h,e);t=0:0.05:24;A=fliplr(A); %将系数数组左右翻转Y=poly2sym(A); %将系数数组转化为多项式subs(Y,'x',t);Y=double(ans);figure(1)plot(x,y,'k*',t,Y,'r-'); %绘制拟合多项式函数图形xlabel('时刻'); %标注坐标轴含义ylabel('平均气温');title([num2str(n-1),'次函数的最小二乘曲线']);grid;%% 指数函数的最小二乘近似yy=log(y);n=3;G=[];GG=[];for j=0:(n-1)g=x.^j; %g(x)按列排列G=vertcat(G,g); %g垂直连接Ggg=t.^j; %g(x)按列排列GG=vertcat(GG,gg); %g垂直连接GendG=G'; %转置得到矩阵Gfor i=1:m %将数据y作为G的最后一列（n+1列）G(i,n+1)=yy(i);endG;for k=1:n %形成矩阵Q(k)if G(k,k)>0;sgn=1;elseif G(k,k)==0;sgn=0;else sgn=-1;endsgm=-sgn*sqrt(sum(G(k:m,k).^2));w=zeros(1,n);w(k)=G(k,k)-sgm;for j=k+1:mw(j)=G(j,k);endbt=sgm*w(k);G(k,k)=sgm; %变换Gk-1到Gkfor j=k+1:n+1t=sum(w(k:m)*G(k:m,j))/bt;for i=k:m;G(i,j)=G(i,j)+t*w(i);endendendA(n)=G(n,n+1)/G(n,n); %解三角方程求系数A for i=n-1:-1:1A (i)=(G(i,n+1)-sum(G(i,i+1:n).*A (i+1:n)))/G(i,i);endb=-A(3);c=A(2)/(2*b);a=exp(A(1)+b*(c^2));G(n+1:m,n+1)=exp(sum(G(n+1:m,n+1).^2));e=sum(G(n+1:m,n+1).^2); %计算误差efprintf('\n指数函数的系数是：a=%f,b=%f,c=%f',a,b,c); %输出系数及误差e fprintf('\n使用指数函数拟合的误差是：%f',e);t=0:0.05:24;YY=a.*exp(-b.*(t-c).^2);figure(2)plot(x,y,'k*',t,YY,'r-'); %绘制拟合指数函数图形xlabel('时刻'); %标注坐标轴含义ylabel('平均气温');title(['指数函数的最小二乘曲线']);grid;4 结果与分析a 二次函数：一天的平均气温为：21.2000二次函数的系数: 8.3063 2.6064 -0.0938使用二次函数拟合的误差是：280.339547二次函数的最小二乘曲线如下图3-1所示：图3-1 二次函数的最小二乘曲线b 三次函数：一天的平均气温为：21.2000三次函数的系数: 13.3880 -0.2273 0.2075 -0.0084 使用三次函数拟合的误差是：131.061822三次函数的最小二乘曲线如图3-2所示：图3-2 三次函数的最小二乘曲线c四次函数：一天的平均气温为：21.2000四次函数的系数: 16.7939 -3.7050 0.8909 -0.0532 0.0009 使用四次函数拟合的误差是：59.04118四次函数的最小二乘曲线如图3-3所示：图3-3 四次函数的最小二乘曲线d 指数函数：一天的平均气温为：21.2000指数函数的系数是： a=26.160286，b=0.004442，c=14.081900 使用指数函数拟合的误差是： 57.034644 指数函数的最小二乘曲线如下图所示：图3-4 指数函数的最小二乘曲线通过上述几种拟合可以发现，多项式的次数越高，计算拟合的效果越好，误差越小，说明结果越准确；同时，指数多项式拟合的次数虽然不高，但误差最小，说明结果最准确。