计算方法(B)上机作业

《计算方法B》上机实习报告学院：电信学院姓名：许梦芸班级：硕5026学号：3115032012代课老师：苏剑1.对以下和式计算：0142118184858616nnSnnnn，要求：(1)若只需保留11个有效数字，该如何进行计算；(2)若要保留30个有效数字，则又将如何进行计算；解答：A:算法实现的思想及依据①首先对问题分析可知：当一般项116𝑛(48𝑛+1−28𝑛+4−18𝑛+5−18𝑛+6)时，截断误差不会对保留数字位数和精度产生影响；②当n=0时，估计出10S01,这表明S要保留m个有效数字，那么小数点后位数应该取到m-1位，因此误差上限=0.1×10−（m−1）。由题可得S(1/16^n)*(4/(8*n+1))保留11个有效数字，即误差ε0.5*10^（1-11）保留30个有效数字，即误差ε0.5*10^（1-30）B:算法实现的结构分为两个部分：顺序结构和循环结构。由于不知道循环次数，因此采用if循环，利用通项是否大于等于误差上限来判断循环条件的真伪。Matalab中可以利用vpa函数控制最终结果的有效数位。C:源程序及相关的注释说明有效数字为11的情况digita=11;%有效数字位数fori=1:1:1000S0=(1/(16^(i-1)))*(4/(8*i-7));S1=(1/16^i)*(4/(8*i+1));if((S0-S1)0.5*10^(1-digita))n=i;break;endendS=0;digits(digita);fori=0:1:nS=S+vpa(1/16^i)*(4/(8*i+1)-2/(8*i+4)-1/(8*i+5)-1/(8*i+6));End（1）运行结果：算法同上所述%有效数字为30的情况digita=30;%有效数字位数fori=1:1:1000S0=(1/(16^(i-1)))*(4/(8*i-7));S1=(1/16^i)*(4/(8*i+1));if((S0-S1)0.5*10^(1-digita))n=i;break;endendS=0;digits(digita);fori=0:1:nS=S+vpa(1/16^i)*(4/(8*i+1)-2/(8*i+4)-1/(8*i+5)-1/(8*i+6));end（2）运行结果：D:结果分析及总结为保持计算精度，计算中对误差上取得较为严格：=0.1×10−（m−1）；由于事先不知道循环次数，采用if循环。编程中出现过错误：对vpa函数理解不当，vpa函数是计算精度控制函数，一般可与digits位控制函数配合使用，如果没有声明，系统默认digits控制精度20位；使用单独使用vpa函数，只控制显示结果的精度，函数内部变量的运算精度不受影响，matlab浮点运算按8字节双精度进行，本题的精度完全可以得到满足。最后，由于本例matlab运算中的精度很高，不会出现大数逐次与小数相加减而导致误差，也没有出现相近数相减而导致精度降低，因此没有进行特别的处理，直接运算，结果精度满足要求。2.某通信公司在一次施工中，需要在水面宽度为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)请用合适的曲线拟合所测数据点；(2)预测所需光缆长度的近似值，并作出铺设河底光缆的曲线图；解答：A:算法实现的思想及依据多项式插值，随着插值数据点的数目增多，误差也会随之增大，因而要用较少的数据点才可避免大的误差，但我们总希望将所得的数据点都用上，分段插值可以解决这一矛盾。分段插值低阶次多项式不具有“光滑性”，因此选用分段三次多项式样条插值。拟合数据点使用多项式插值时由于龙格现象误差较大，故考虑使用三次样条插值。B:算法实现的结构参考教材给出的SPLINEM算法和TTS算法，在选定边界条件和选定插值点等距分布后，可以先将数据点的二阶差商求出来并赋值给右端向量d，再根据TSS解法求解M。求出M后利用曲线积分公式带入复化梯形公式求解积分，求解各插值点处的值和导数值，采用内外两重积分的方法。C:源程序及相关的注释说明%------三次样条插值多项式------clc;clear;%-----插值函数和插值点---------n=20;%插值点个数11i=0:n;%子元素从0到nx=10.*i;f=[9.018.967.967.978.029.0510.1311.1812.2613.2813.3212.6111.2910.229.157.907.958.869.8110.8010.93];%--输出插值点---fprintf('当n=%d时\n',n)fprintf('插值点x坐标为：\n')disp(x);fprintf('插值点y坐标f（x）为：\n')disp(f);%----计算差商表，并将有用的差商进行选择和保存-----y=f;%存放差商d=zeros(1,n+1);%存放三弯矩方程组的右端项fori=1:n%计算第i阶差商forj=n+1:-1:i+1y(j)=(y(j)-y(j-1))/(x(j)-x(j-i));endifi==2%将第2阶差商赋值给dd(2:n)=6.*y(3:n+1);endifi==3%三次样条插值的第三类边界条件d(1)=-12*(x(2)-x(1))*y(i+1);d(n+1)=12*(x(n+1)-x(n))*y(n+1);endend%-------用追赶法求三弯矩方程组,Tss------a=ones(1,n+1);%存储下三角的非零元素b=2*ones(1,n+1);%存储对角元素c=ones(1,n+1);%存储上三角的非零元素fori=2:na(i)=(x(i)-x(i-1))/(x(i+1)-x(i-1));c(i)=1-a(i);enda(n+1)=-2;c(1)=-2;%-------计算矩阵L和U-------------%初始化L和Ul=zeros(1,n+1);u=zeros(1,n+1);r=zeros(1,n+1);u(1)=b(1);r(1)=c(1);fork=2:n+1l(k)=a(k)/u(k-1);u(k)=b(k)-l(k)*r(k-1);r(k)=c(k);end%--------使用追赶法进行求解------------M=zeros(1,n+1);%方程组的解M(1)=d(1);fork=2:n+1M(k)=d(k)-l(k)*M(k-1);endM(n+1)=M(n+1)/u(n+1);fork=n:-1:1M(k)=(M(k)-r(k)*M(k+1))/u(k);endfprintf('三弯矩方程组的解M为：\n')disp(M);%------计算f（x）、N(x)和S(x)---------Sy1=zeros(1,200);x1=zeros(1,200);fori=1:200x1(i)=i-1;endfori=1:200forj=1:nifx1(i)=x(j)&&x1(i)=x(j+1)temp=x(j+1)-x(j);Sy1(i)=M(j).*(x(j+1)-x1(i)).^3/6/temp+M(j+1).*(-x(j)+x1(i)).^3/6/temp+(f(j)-M(j)*temp^2/6).*(x(j+1)-x1(i))/temp+(f(j+1)-M(j+1)*temp^2/6).*(-x(j)+x1(i))/temp;endendendfprintf('由三次样条插值方法计算的各年的产量为：\n')disp(Sy1)plot(x,f,'*')%画出由插值点构成的曲线holdonplot(x1,Sy1,'r')%用红线画出三次样条插值拟合的曲线holdonlegend('三次样条插值');%%%%%%%%%%%求解电缆长度&&&&&&&&&&&&&&&&&&&&&&&&Cordinate=[0:20/(length(Sy1)-1):20;Sy1];diff_Cordinate=[diff(Cordinate(1,:));diff(Cordinate(2,:))];disp('电缆长度为：');anwser=sum(sqrt(diff_Cordinate(1,:).*diff_Cordinate(1,:)+diff(Cordinate(2,:)).*diff(Cordinate(2,:))))D:运行结果及分析：可以看出，三次样条插值拟合数据点比较好，没有龙格现象。3.假定某天的气温变化记录如下表所示，试用数据拟合的方法找出这一天的气温变化的规律;试计算这一天的平均气温,并试估计误差。时刻0123456789101112平均气温15141414141516182020232528时刻131415161718192021222324平均气温313431292725242220181716解答：A：算法实现的思想及依据本题共有25个数据点，数据量较大，因此选用多项式插值会导致较大的误差，插值样条函数虽然有较好的数值性质，但是Mi的计算量不小，而且没有统一的公式来表示。另外，本题要求找出温度变化规律，插值方法要求p(x)必须通已知数据点，只会导致拟合曲线失去本有的数据所表示的规律；从表中的数据点可以看出，温度并不准确（温度只测量到整数位），因此没有必要要求拟合曲线通过数据点。因此，选取最小二乘近似。B:算法实现的结构算法参考课本LSS算法。利用已有数据来生成了G，将y作为其最后一列（第n+1列）存放。先形成矩阵Qk，再变换Gk-1到Gk；求解三角方程Rx=h1；最后根据定义计算误差。在计算平均温度时，先将积分求出来，由于多项式为简单的4阶，因此直接根据牛顿莱布尼兹公式求解，避免数值积分的繁琐；再将积分求均值，获得平均温度。C:源程序及相关的注释说明clear%清除Workspace中的变量clc%清除CommandWindow的内容x=0:1:24;%时间序列矩阵Y=[151414141415161820202325283134....31292725242220181716];%原始温度值m=25;%一共25个数据点n=5;%定义多项式拟合的阶数为4G=zeros(m,n+1);forj=1:n%定义最小二乘法的G矩阵foru=1:mG(u,j)=x(u)^(j-1);endendfori=1:mG(i,n+1)=Y(i);%将y存放在矩阵G的最后一列endsigma=0;w=ones(m,1);fork=1:n%形成正交变换矩阵Qkfori=k:msigma=sigma+G(i,k)^2;endsigma=-sigma^0.5;w(k)=G(k,k)-sigma;forj=(k+1):mw(j)=G(j,k);endbeta=sigma*w(k);G(k,k)=sigma;%变换Gk-1到Gkforj=(k+1):mG(j,k)=0;endforj=(k+1):n+1t=0;fori=k:mt=t+w(i)*G(i,j)/beta;endfori=k:m%更新Gk-1矩阵右下角(m-k+1)(n-k+1)列的值G(i,j)=G(i,j)+t*w(i);endendsigma=0;%每经历一次循环，sigma强制归零，避免上一次循环累加到下一次循环endalpha=ones(n,1);alpha(n)=G(n,n+1)/G(n,n);fori=n-1:-1:1V=0;forj=i+1:nV=V+G(i,j)*alpha(j);endalpha(i)=(G(i,n+1)-