1.计算以下和式:0142118184858616nn S n n n n ∞=⎛⎫=--- ⎪++++⎝⎭∑,要求: (1)若保留11个有效数字,给出计算结果,并评价计算的算法;(2)若要保留30个有效数字,则又将如何进行计算。
(1)题目分析该题是对无穷级数求和,因此在使用matlab 进行累加时需要一个累加的终止条件。
这里令⎪⎭⎫ ⎝⎛+-+-+-+=681581482184161n n n n a nn ,则()()1.01616855844864816114851384128698161 681581482184161148113811282984161111<<⎪⎪⎭⎫ ⎝⎛⎪⎭⎫ ⎝⎛++++++⎪⎪⎭⎫ ⎝⎛⎪⎭⎫ ⎝⎛++++++=⎪⎭⎫⎝⎛⎪⎭⎫ ⎝⎛+-+-+-+⎪⎭⎫ ⎝⎛⎪⎭⎫ ⎝⎛+-+-+-+=+++n n n n n n n n n n n n n n n n a a n n n n n n 故近似取其误差为1+≈k a ε,并且有m-1m -11102121⨯=⨯=≈+βεk a ,(2)算法依据使用matlab 编程时用digits 函数和vpa 函数来控制位数。
(3)Matlab 运行程序%%保留11位有效数字 k1=11;s1=0;%用于存储这一步计算值 for n=0:50a=(1/16^n)*(4/(8*n+1)-2/(8*n+4)-1/(8*n+5)-1/(8*n+6)); n1=n-1;if a<=0.5*10^(1-k1) break end end;for i=0:1:n1t=(1/16^i)*(4/(8*i+1)-2/(8*i+4)-1/(8*i+5)-1/(8*i+6)); s1=s1+t; ends11=vpa(s1,k1);disp('保留11位有效数字的结果为:');disp(s11); disp('此时n 值为:');disp(n1);%%保留30位有效数字 clear all; k2=30;digits(k2+2);s2=vpa(0);%用于存储这一步计算值for n=0:50a=vpa((1/16^n)*(4/(8*n+1)-2/(8*n+4)-1/(8*n+5)-1/(8*n+6)));n2=n-1;if a<=0.5*10^(1-k2)breakendend;for i=0:1:n2t=vpa((1/16^i)*(4/(8*i+1)-2/(8*i+4)-1/(8*i+5)-1/(8*i+6)));s2=vpa(s2+t);ends30=vpa(s2,k2);disp('保留30位有效数字的结果为:');disp(s30);disp('此时n值为:');disp(n2);2.某通信公司在一次施工中,需要在水面宽度为20米的河沟底部沿直线走向铺设一条沟底光缆。
在铺设光缆之前需要对沟底的地形进行初步探测,从而估计所需光缆的长度,为工程预算提供依据。
已探测(1)(2)估算所需光缆长度的近似值,并作出铺设河底光缆的曲线图;(2)算法选择由于多项式插值只适用于数据点较少的插值案例,此时可以避免产生较大的误差;当数据点过多时,多项式插值中的龙格现象会使误差变大。
因此当所用数据点过多时,我们采用分段插值方式,即用分段多项式代替单个多项式作插值。
分段多项式是由一些在相互连接的区间上的不同多项式连接而成的一条连续曲线,其中三次样条插值方法是一种具有较好“光滑性”的分段插值方法。
(3)matlab算法clear all;x=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.90,7.95,8.86,9.81,10.80,10.93];%对应水面宽度处的水深n=length(x);%取样点数目M=y;for k=2:3;%计算二阶差商for i=n:-1:k;M(i)=(M(i)-M(i-1))/(x(i)-x(i-k+1));endendh(1)=x(2)-x(1);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);endx1=1:0.1:20; %确定取值点n1=length(x1); %确定取值点个数s=zeros(1,n1);for m=1:n1; k=1;for i=2:n-1if x1(m)<=x(i); k=i-1;break; else k=i; end endh1=x(k+1)-x(k);%在各区间用三次样条插值函数计算X 点处的值 x11=x(k+1)-x1(m);x12=x1(m)-x(k);sum(m)=(M(k)*(x11^3)/6+M(k+1)*(x12^3)/6+(y(k)-(M(k)*(h1^2)/6))*x11+(y (k+1)-(M(k+1)*(h1^2)/6))*x12)/h1; endl=0; %计算所需光缆长度 for i=2:n1l=l+sqrt((x1(i)-x1(i-1))^2+(sum(i)-sum(i-1))^2); enddisp('所需光缆长度为'); disp(l); figure(1)plot(x,y,'*',x1,sum,'-') xlabel('位置'); ylabel('深度/m');title('铺设河底光缆的曲线图'); grid on;(4)结果分析铺设海底光缆的曲线图如下图:位置深度/m铺设河底光缆的曲线图使用matlab 仿真计算出所需光缆的长度为 25.4525m 。
3.假定某天的气温变化记录如下表所示,试用数据拟合的方法找出这(1) 选择算法在本题中,采用“最小二乘法”来近似这一天的气温变化规律,分别使用二次函数、三次函数、五次函数以及指数型函数来对数据点进行拟合,并计算相应的误差。
(2)算法框架a )使用多项式进行最小二乘近似时,首先生成矩阵G ,然后利用算法LSS 求解多项式的系数和误差22E 。
b )使用指数函数进行最小二乘近似时,设其指数形式为2321d x d x d ey +=,再对函数形式进行变换332211a ,a ,ln a ,ln d d d y z ====,即2321a a a x x z ++=。
然后生成矩阵G ,再利用算法LSS 求解多项式的系数和误差22E 。
(3)matlab 程序%%多项式的最小二乘近似 clear all; 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); %将数据点存入ms=0; %求一天的平均气温for i=1:1:24s=s+y(i);endT=s/m;disp('一天的平均气温为');disp(T);n=h+1; %h次函数的最小二乘近似,h分别取2、3、5for j=1:n; %生成矩阵Gfor i=1:25G(i,j)=x(i).^(j-1);endendfor i=1:m; %将y作为G的最后一列存放G(i,j+1)=y(i);endfor k=1:nif G(k,k)>0 %形成矩阵Q(k)sgn=1;elseif G(k,k)==0;sgn=0;else sgn=-1;enda1=-sgn*sqrt(sum(G(k:m,k).^2));w=zeros(1,n);w(k)=G(k,k)-a1;for j=k+1:mw(j)=G(j,k);endb1=a1*w(k);G(k,k)=a1; %变换Gk-1到Gkfor j=k+1:n+1t=sum(w(k:m)*G(k:m,j))/b1;for i=k:mG(i,j)=G(i,j)+t*w(i);endendenda(n)=G(n,n+1)/G(n,n); %求解方程for 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); %计算误差t=0:0.01:24;a=fliplr(a); %将系数数组左右翻转y1=poly2sym(a); %将系数数组转化为多项式subs(y1,'x',t);y1=double(ans);figure(1)plot(x,y,'*',t,y1,'r-');xlabel('时刻');plot(x,y,'*',t,y1,'r-');xlabel('时刻');ylabel('温度');title('5次函数的最小二乘曲线');grid on;disp('误差为');disp(E);%%指数形式的最小二乘近似clear all;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=3;z=zeros(1,m);for i=1:mz(i)=log(y(i));endfor j=1:n; %生成矩阵Gfor i=1:m;G(i,j)=x(i).^(j-1);endendfor i=1:m; %将y作为G的最后一列存放G(i,j+1)=z(i);endfor k=1:nif G(k,k)>0; %形成矩阵Q(k)sgn=1;elseif G(k,k)==0;sgn=0;else sgn=-1;enda1=-sgn*sqrt(sum(G(k:m,k).^2));w=zeros(1,n);w(k)=G(k,k)-a1;for j=k+1:mw(j)=G(j,k);endb1=a1*w(k);G(k,k)=a1; %变换Gk-1到Gkfor j=k+1:n+1t=sum(w(k:m)*G(k:m,j))/b1;for i=k:m;G(i,j)=G(i,j)+t*w(i);endendenda(n)=G(n,n+1)/G(n,n); %解方程组for i=n-1:-1:1a(i)=(G(i,n+1)-sum(G(i,i+1:n).*a(i+1:n)))/G(i,i);endd3=a(3);d2=a(2);d1=exp(a(1));E=0;for i=1:25y1(i)=d1*exp(d2*x(i)+d3*x(i)*x(i));E=(y1(i)-y(i)).^2+E; %计算误差endt=0:0.01:24;y1=d1.*exp(d2.*t+d3.*t.*t);fprintf('指数函数拟合的系数是:d1=%f,d2=%f,d3=%f',d1,d2,d3); figure(1)plot(x,y,'*',t,y1,'r-');xlabel('时刻');ylabel('温度');title('指数函数的最小二乘曲线');grid on;fprintf('\n指数函数拟合的系数是:%f',E);(4)结果分析a)二次多项式的最小二乘近似一天的平均气温为:20.5600误差为:280.3395拟合曲线如下图所示:时刻温度2次函数的最小二乘曲线b ) 三次多项式的最小二乘近似 误差为:131.0618 近似函数曲线如下图所示:时刻温度3次函数的最小二乘曲线c )五次多项式的最小二乘近似 误差为:33.1446近似函数曲线如下图所示:0510152025时刻温度5次函数的最小二乘曲线d )指数形式的最小二乘近似指数函数拟合的系数是:d1=10.842462,d2=0.125093,d3=-0.004442 指数函数拟合的误差E 2是:215.700187 近似函数曲线如下图所示:0510152025时刻温度指数函数的最小二乘曲线从上述几种拟合可以得出,多项式的次数越高,拟合的效果越好,误差越小;指数函数的二次多项式拟合相对于同阶的多项式最小二乘近似其误差相对较小,但不如高阶多项式最小二乘近似精确。