计算方法上机作业

计算方法上机报告姓名：学号：班级：上课班级：计算方法上机报告1说明：本次上机实验使用的编程语言是Matlab语言，编译环境为MATLAB7.11.0，运行平台为Windows7。1.对以下和式计算：0681581482184161nnnnSn，要求：①若只需保留11个有效数字，该如何进行计算；②若要保留30个有效数字，则又将如何进行计算；（1）算法思想1、根据精度要求估计所加的项数，可以使用后验误差估计，通项为：142111416818485861681nnnannnnn；2、为了保证计算结果的准确性，写程序时，从后向前计算；3、使用Matlab时，可以使用以下函数控制位数：digits(位数)或vpa(变量，精度为数)（2）算法结构1.;0s681581482184161nnnntn；2.for0,1,2,,niif10mtend;3.for,1,2,,0niii;sst计算方法上机报告2（3）Matlab源程序clear;%清除工作空间变量clc;%清除命令窗口命令m=input('请输入有效数字的位数m=');%输入有效数字的位数s=0;forn=0:50t=(1/16^n)*(4/(8*n+1)-2/(8*n+4)-1/(8*n+5)-1/(8*n+6));ift=10^(-m)%判断通项与精度的关系break;endend;fprintf('需要将n值加到n=%d\n',n-1);%需要将n值加到的数值fori=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值加到n=7，s=3.1415926536；当保留30位有效数字时，需要将n值加到n=22,s=3.14159265358979323846264338328。通过上面的实验结果可以看出，通过从后往前计算，这种算法很好的保证了计算结果要求保留的准确数字位数的要求。计算方法上机报告32.某通信公司在一次施工中，需要在水面宽度为20米的河沟底部沿直线走向铺设一条沟底光缆。在铺设光缆之前需要对沟底的地形进行初步探测，从而估计所需光缆的长度，为工程预算提供依据。已探测到一组等分点位置的深度数据(单位:米)如下表所示：分点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）算法思想如果使用多项式差值，则由于龙格现象，误差较大，因此，用相对较少的插值数据点作插值，可以避免大的误差，但是如果又希望将所得数据点都用上，且所用数据点越多越好，可以采用分段插值方式，即用分段多项式代替单个多项式作插值。分段多项式是由一些在相互连接的区间上的不同多项式连接而成的一条连续曲线，其中三次样条插值方法是一种具有较好“光滑性”的分段插值方法。在本题中，假设所铺设的光缆足够柔软，在铺设过程中光缆触地走势光滑，紧贴地面，并且忽略水流对光缆的冲击。海底光缆线的长度预测模型如下所示，光缆从A点铺至B点，在某点处的深度为ih。计算方法上机报告4海底光缆线的长度预测模型计算光缆长度时，用如下公式：200()Lfxds20'20()1()fxfxdx191'20()1()kkkfxfxdx22()()xy（2）算法结构1.Forni,,2,1,01.1iiMy2.For2,1k2.1Forknni,,1,2.1.1ikiiiiMxxMM)/()(13.101hxx4.For1-,,2,1ni4.111iiihxx4.2bacchhhiiiiii2;1;)/(11计算方法上机报告54.3iidM165.0000;;cMdMdnnnnnbab2;;206.1111,db7.获取M的矩阵元素个数，存入m8.Formk,,3,28.1kkkla1/8.2kkkkclb1-8.3kkkkld1-9.mmmM/10.For1,,2,1mmk10.1kkkkkMMc/)(111.获取x的元素个数存入s12.k113.For1,,2,1si13.1ifixx~thenki；breakelseki114.xxxxxxhxxkkkkˆ~;~;11yhxhMyxhMyxMxMkkkkkk~/]ˆ)6()6(6ˆ6[2211331（3）Matlab源程序clear;clc;x=0:1:20;%产生从0到20含21个等分点的数组计算方法上机报告6X=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;fork=2:3;%计算二阶差商并存放在M中fori=n:-1:k;M(i)=(M(i)-M(i-1))/(x(i)-x(i-k+1));endendh(1)=x(2)-x(1);%计算三对角阵系数a,b,c及右端向量dfori=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;计算方法上机报告7b(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);fork=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)fork=n-1:-1:1;M(k)=(y1(k)-c(k)*M(k+1))/u(k);ends=zeros(1,N);form=1:N;k=1;fori=2:n-1ifX(m)=x(i);计算方法上机报告8k=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;%计算所需光缆长度fori=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);%标注坐标轴含义计算方法上机报告9ylabel('深度/m','fontsize',16);title('铺设河底光缆的曲线图','fontsize',16);grid;（4）结果与分析铺设海底光缆的曲线图如下图所示：仿真结果表明，运用分段三次样条插值所得的拟合曲线能较准确地反映铺设光缆的走势图，计算出所需光缆的长度为L=26.4844m。3.假定某天的气温变化记录如下表所示，试用数据拟合的方法找出这一天的气温变化的规律；试计算这一天的平均气温，并试估计误差。时刻0123456789101112平均气温15141414141516182020232528时刻131415161718192021222324平均气温313431292725242220181716计算方法上机报告10（1）算法思想在本题中，数据点的数目较多。当数据点的数目很多时，用“多项式插值”方法做数据近似要用较高次的多项式，这不仅给计算带来困难，更主要的缺点是误差很大。用“插值样条函数”做数据近似，虽然有很好的数值性质，且计算量也不大，但存放参数Mi的量很大，且没有一个统一的数学公式来表示，也带来了一些不便。另一方面，在有的实际问题中，用插值方法并不合适。当数据点的数目很大时，要求)(xp通过所有数据点，可能会失去原数据所表示的规律。如果数据点是由测量而来的，必然带有误差，插值法要求准确通过这些不准确的数据点是不合适的。在这种情况下，不用插值标准而用其他近似标准更加合理。通常情况下，是选取i使2E最小，这就是最小二乘近似问题。在本题中，采用“最小二乘法”找出这一天的气温变化的规律，使用二次函数、三次函数、四次函数以及指数型函数2)(ctbaeC，计算相应的系数，估算误差，并作图比较各种函数之间的区别。（2）算法结构本算法用正交化方法求数据的最小二乘近似。假定数据以用来生成了G，并将y作为其最后一列（第1n列）存放。结果在a中，是误差22E。I、使用二次函数、三次函数、四次函数拟合时1.将“时刻值”存入ix，数据点的个数存入m2.输入拟合多项式函数)(xp的最高项次数1nh，则拟合多项式函数为)()()()(2211xgxgxgxpnn,根据给定数据点确定GFor1,,2,1,0njFormi,,2,12.11,jijigx2.21,niigy计算方法上机报告113.Fornk,,2,13.1[形成矩阵kQ]3.1.1mkiikkkgg2/12))(sgn(3.1.2kkkg3.1.3Formkkj,,2,13.1.3.1jjkg3.1.4k3.2[变换1kG到kG]3.2.1kkgFor1,,,2,1nnkkj3.2.2tgmkiiji/)(3.2.3Formkki,,1,3.2.3.1ijiijgtg4.[解三角方程1hRa]4.1nnnnnagg/1.4.2For1,,2,1nni4.2.1iiinijjijniagxgg/][11.5.[计算误差22E]mninig121,II、使用指数函数拟合时现将指数函数进行变形：计算方法上机报告12将yC,xt代入2)(ctbaeC得:2)(cxbaey对上式左右取对数得:222lnlnbxbcxbcay令bbcbcayz21202lnln，，，则可得多项式:2210xxz（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);%将数据点的个数存入mT=sum(y(1:m))/m;fprintf('一天的平均气温为T=%f\n',T);%求一天的平均气温%%二次、三次、四次函数的最小二乘近似h=input('请输入拟合多项式的最高项次数=');%根据给定数据点生成矩阵