西安交大计算方法b大作业

西安交大计算方法b大作业 《计算方法B》上机 实验 学院: 机械学院 班级: 姓名: 学号: 2015年12月22日 1 ,14211,,1.计算以下和式:，要求: S,,,,,,,nnnnn81848586，，，，160,,n, (1)若保留11个有效数字，给出计算结果，并评价计算的算法; (2)若要保留30个有效数字，则又将如何进行计算。 实现思想: 以上问题出现了近似数相减的问题，为了减小误差，可分别求得减数之和以及被减数之和，最后将两者相减。另外，减数与被减数求和均为同号计算，按照绝对值递增顺序相加可减小舍入误差。此题中对有效数字有要求，因而计算时首先需要根据有效数字位数计算得出迭代次数，以保证计算值的精度。 源程序: m=input('输入有效数字个数m='); s0=1;s1=0;s2=0;n=0; %判断迭代次数 while s0>=0.5*10^-(m-1) s0=4/(16^n*(8*n+1))-2/(16^n*(8*n+4))-1/(16^n*(8*n+5))-1/(16^n*(8*n+6) ); n=n+1; end %分别求解各项并求和 for k=n-1:-1:0 a1=4/(16^k*(8*k+1)); a2=2/(16^k*(8*k+4)); a3=1/(16^k*(8*k+5)); a4=1/(16^k*(8*k+6)); s1=a1+s1; s2=a4+a3+a2+s2; end S=vpa(s1-s2,m) 2 实验结果:11位有效数字计算结果如图1所示;30为有效数字计算结果如图2所示。 图1.11位有效数字计算结果 图2.30为有效数字计算结果 3 1. 某通信公司在一次施工中，需要在水面宽度为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)估算所需光缆长度的近似值，并作出铺设河底光缆的曲线图; 算法思想:由于题中所给点数为20，若采用高次多项式插值将产生很大的误差，所以拉格朗日或牛顿并不适用。题中光缆为柔性，可光滑铺设于水底，鉴于此特性，采用三次样条插值插值法较为合适。 算法结构:三次样条算法结构见《计算方法教程》P110; 光缆长度计算公式: 202019k，1'2'2l=1+(f())1+(f())xdxxdx, ,,,,k,k000 源程序: clear; clc; 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)); end end %%%设定d的边界条件 for i=2:20 4 d(i)=6*d(i+1); end d(1)=0; d(21)=0; %%%带状矩阵求解(追赶法) 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=linspace(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 end h=1; xbar=x(k)-xx(j); xmao=xx(j)-x(k-1); s(j)=(m(k-1)*xbar^3/6+m(k)*xmao^3/6+(y(k-1)-m(k-1)*h^2/6)*xbar+(y(k)- m(k)*h^2/6)*xmao)/h; 5 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(j+1)=(1+sp(j)^2)^0.5*(20/nn)+l(j); %求解光缆长度 end %%%绘图 plot(xx,s,'r-','linewidth',1.5) disp(['?光缆长度为ª',num2str(l(nn+1)),'Ã×']) 曲线图如图2-1所示，计算光缆长度如图2-2所示。 图2-1光缆插值曲线图 图2-1光缆计算长度显示 6 3.假定某天的气温变化记录如下表所示，试用数据拟合的方法找出这一天的气温变化的规律;试计算这一天的平均气温,并试估计误差。 时刻 0 1 2 3 4 5 6 7 8 9 10 11 12 平均气温 15 14 14 14 14 15 16 18 20 20 23 25 28 时刻 13 14 15 16 17 18 19 20 21 22 23 24 平均气温 31 34 31 29 27 25 24 22 20 18 17 16 实现思想:此题中所给数据点数目较多，采用拉格朗日插值法或者牛顿插值法需要很高次的多项式，计算困难，误差大;采用样条插值计算量虽然不大，但是存放参数Mi的量很大，且没有一个统一的公式来表示，也不是很方便。所以可考虑用最小二乘法进行拟合。计算过程中，分别使用二次函数、三次函数以及四次函数，计算其相应的系数，估算误差并作图比较各个函数之间的区别。 算法结构:(参考课本P123) 1.1[形成矩阵Qk] 1.2[变换Gk-1到Gk] 2.[求解三角方程] 3.[计算误差] 源代码: 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=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用来存放β %%%形成矩阵G 7 for j=1:n for i=1:m G(i,j)=x(1,i)^(j-1); end end %%%建立矩阵Qk for k=1:n for 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 end end %%%求解三角方程 Rx=h1 a(n)=G(n,n+1)/G(n,n); for i=n-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); %a(i)存放各级系数 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 8 plot(x,p,'r*',x,p,'-'); E2=0;%用E2来存放误差 %%%误差求解 for i=n+1:m E2=G(i,n+1)^2+E2; end E2=E2^0.5; disp('误差为'); disp(E2); t=0; for i=1:m t=t+p(i); end t=t/m; %%%平均温度 disp(['平均温度为',num2str(t),'?']) 实验结果: 二次函数拟合，结果如下图所示 图3-1 二次函数拟合结果 9 三次函数拟合，结果如下图所示 图3-2 三次函数拟合结果 四次函数拟合，结果如下图所示 图3-3 四次函数拟合结果 10 结果对比:将二次函数、三次函数和四次函数拟合结果绘制在同一个坐标内，如图3-4所示。其计算误差结果见表3-1所示。 图3-4 拟合结果对比分析 11 524.设计算法，求出非线性方程的所有实根，并使645200xx,，, ,4误差不超过。 10 52算法思想:本题可采用牛顿法迭代求解,令 ,得带格式为 f(x)=6x-45x+20 f(x)k x=x-k+1k'f(x)k 根据函数图像可以找出根的大致分布区间,带入不同的初值即可解出不同的根. 源代码: function y=f2(x) y=6*x.^5-45*x.^2+20;%定义原函数 function y=f3(x) y=30*x^4-90*x;%定义原函数倒数 i=-5:0.1:5; y=f2(i); plot(i,y) hold on plot(i,0,'-')%画出原函数图像 %%Newton法求根 x1=input('输入初值'); e=10^(-4);%误差设定 Nmax=1000;%迭代最大次数限定 for n=1:Nmax f0=f2(x1); if abs(f2(x1))