计算方法B上机报告

计算方法B上机报告第1题某通信公司在一次施工中，需要在水面宽度为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)估算所需光缆长度的近似值，并作出铺设河底光缆的曲线图；问题分析和算法思想：本题的主要目的是对 21个测量数据进行拟合，同时对拟合曲线进行线积分即可得到河底光缆长度的近似值，可以用的插值方法很多：多项式插值、Lagrange插值、Newton插值、三次样条插值等。由于数值点较多时，采用高次多项式插值将产生很大的误差，用拉格朗日插值多项式会出现龙格现象。故为了将所有的数据点都用上，且题中光缆为柔性，可光滑铺设于水底，鉴于此特性，采用三次样条插值的方法较为合适。计算光缆长度近似值，只需将每两点之间的距离算出，然后依次相加，所得的折线长度，即为光缆长度的近似值。光缆长度计算公式：20219k1202l1fxdx1fxdx0k0k0算法结构：三次样条算法结构见《计算方法教程》 P110。源程序：clear;clc;x=0:20;y=[9.018.967.967.978.029.0510.1311.1812.2613.2813.3212.6111.2910.229.157.907.958.869.8110.8010.93];d=y;plot(x,y,'k.','markersize',15)holdon精品文档交流%%%计算二阶差商fork=1:2fori=21:-1:(k+1)d(i)=(d(i)-d(i-1))/(x(i)-x(i-k));endend%%%假定d的边界条件，采用自然三次样条fori=2:20d(i)=6*d(i+1);endd(1)=0;d(21)=0;%%%追赶法求解带状矩阵的 m值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);%%%追的过程fork=2:21l(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);fork=20:-1:1m(k)=(yy(k)-r(k)*m(k+1))/u(k);end%%%利用插值点画出拟合曲线k=1;nn=100;xx=linspace(0,20,nn);l=0;forj=1:nnfori=2:20ifxx(j)<=x(i)k=i;break;elsek=i+1;endend精品文档交流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;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%%%绘图title('光缆插值曲线')xlabel('分点')ylabel('深度')plot(xx,s,'r-','linewidth',1.5)griddisp(['所需光缆长度为',num2str(l(nn+1)),'米'])运行结果：光缆插值曲线14131211度深1098724681012141618200分点第2题假定某天的气温变化记录如下表所示，试用数据拟合的方法找出这一天的气温变化的规律;试计算这一天的平均气温,并试估计误差。时刻0123456789101112平均气温15141414141516182020232528精品文档交流时刻131415161718192021222324平均气温313431292725242220181716问题分析及算法思路：由于题中所给数据点有25个，此时采用多项式插值的方法做数据近似需要用较高次的多项式，这不仅给计算带来困难，而且有余项带来的误差很大；若采用样条插值计算量虽然不大，但是参数的个数很多，而且没有一个统一的数学公式来表示，对计算带来了很大的麻烦。所以可考虑使用最小二乘法进行拟合。在本题中，采用最小二乘法找出这一天的气温变化的规律，依次使用二次函数、三次函数以及四次函数进行拟合，计算其相应的系数，估算误差，并作图比较三个函数之间的区别。算法结构：最小二乘法算法结构见《计算方法教程》 P123。源程序：x=0:24;y=[15141414141516182020232528313431292725242220181716];m=length(x);n=input('请输入函数的次数 ');plot(x,y,'k.',x,y,'-')grid;holdon;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用来存放β%%%形成矩阵Gforj=1:nfori=1:mG(i,j)=x(1,i)^(j-1);endend%%%建立矩阵Qkfork=1:nfori=k:mc(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来存放ωforj=k+1:mw(j)=G(j,k);endb(k)=c(k)*w(k);%%%变换矩阵Gk-1到GkG(k,k)=c(k);forj=k+1:n+1q=0;fori=k:mq=w(i)*G(i,j)+q;ends=q/b(k);fori=k:mG(i,j)=s*w(i)+G(i,j);endendend%%%求解三角方程 Rx=h1a(n)=G(n,n+1)/G(n,n);fori=n-1:(-1):1forj=i+1:nf=G(i,j)*a(j)+f;enda(i)=(G(i,n+1)-f)/G(i,i);%%%a(i)存放各级系数f=0;end%%%拟合后的回代过程p=zeros(1,m);forj=1:mfori=1:np(j)=p(j)+a(i)*x(j)^(i-1);endendplot(x,p,'r*',x,p,'-');精品文档交流E2=0;%用E2来存放误差%%%误差求解fori=n+1:mE2=G(i,n+1)^2+E2;endE2=E2^0.5;disp('误差为');disp(E2);t=0fori=1:mt=t+p(i);endt=t/m;%%%平均温度disp(['平均温度为',num2str(t),'℃'])运行结果：二次函数时，系数 a0=8.3063，a1=2.6064，a2=-0.0938，误差为16.7433，平均温度为21.2℃。函数形式为： p(x) 8.3063 2.6064x 0.0938x2三次函数时，系数a0=13.3880，a1=-0.2273，a2=0.2075，a3=-0.0084，误差为11.4482，平均温度为21.2℃。函数形式： p(x) 13.3880 0.2273x 0.2075x2 0.0084x3精品文档交流四次函数时，系数a0=16.7939，a1=-3.7050，a2=0.8909，a3=-0.0532，a4=0.0009，误差为7.6838，平均温度为21.2℃。函数形式为： p(x) 16.7939 3.7050x 0.8909x2 0.0532x3 0.0009x4第3题10cos(xsin)d已知非线性方程0在[2,3]中有根。请设计算法，求出该根，并使求出的根的误差不超过104。问题分析以及算法思路：精品文档交流f(x)1cosxsind0本题采用复化梯形求积公式计算方程根，可令。算法结构：牛顿迭代法算法结构见《计算方法教程》 P196。源程序：clear;clc;x0=2; %设置求解区间x1=3;n=20;%复化梯形求积公式的步数e=0.0004;f0=T3Sub(x0,n);f1=T3Sub(x1,n);if(f0*f1)>0 %判断在区间内是否有error('Noanswer,becausef0*f1>0');endfori=1:201x=(x0+x1)/2;ifabs(x1-x)<e %判断解是否符合误差disp(['方程的解',num2str(x)]);disp(['迭代步数',num2str(i)]);disp(['误差',num2str(abs(x1-x))]);break;endf=T3Sub(x,n);if(f*f0)<0x1=x;f1=f;elsex0=x;f0=f;endendfunction[f]=T3Sub(x,n)%复化梯形求积公式h=pi/n;f=0;fori=1:nu=i*h;f=f+cos(x*sin(u-h))+cos(x*sin(u));endf=(h*f)/(2*pi);end运行结果：精品文档交流第4题线性方程组求解1）编写程序实现大规模方程组的高斯消去法程序，并对所附的方程组进行求解。所附方程组的类型为对角占优的带状方程组。2）针对本专业中所碰到的实际问题，提炼一个使用方程组进行求解的例子，并对求解过程进行分析、求解。算法思想:高斯消去法是利用现行方程组初等变换中的一种变换 ,将一个不为零的数乘到一个方程后加到另一个方程,使方程组变成同解的上三角方程组 ,然后再自下而上对上三角方程组求解。算法结构：1.读取二进制文件，存入计算矩阵2.对矩阵进行初等变换，然后求解 (计算方法教程高斯消去法及列主元高斯消去法算法)源代码：clear;clc;读取系数矩阵[f,p]=uigetfile('*.dat','选择数据文件');%读取数据文件num=5;%输入系数矩阵文件头的个数name=strcat(p,f);file=fopen(name,'r');head=fread(file,num,'uint');%读取二进制头文件id=dec2hex(head(1));%读取标识符fprintf('文件标识符为：');idver=dec2hex(head(2));%读取版本号fprintf('文件版本号为：');vern=head(3);%读取阶数fprintf('矩阵A的阶数：');nq=head(4);%上带宽精品文档交流fprintf('矩阵A的上带宽：');qp=head(5);%下带宽fprintf('矩阵A的下带宽：');pdist=4*num;fseek(file,dist,'bof');%把句柄值转向第六个元素开头处[A,count]=fread(file,inf,'float');%读取二进制文件，获取系数矩阵fclose(file);%关闭二进制头文件对非压缩带状矩阵进行求解ifver=='102',a=zeros(n,n);fori=1:n,forj=1:n,a(i,j)=A((i-1)*n+j);%求系数矩阵a(i,j)endendb=zeros(n,1);fori=1:n,b(i)=A(n*n+i);endfork=1:n-1,%列主元高斯消去法m=k;fori=k+1:n,%寻找主元ifabs(a(m,k))<abs(a(i,k))m=i;endendifa(m,k)==0%遇到条件终止disp('错误！')returnendforj=1:n,%交换元素位置得主元t=a(k,j);a(k,j)=a(m,j);a(m,j)=t;t=b(k);b(k)=b(m);b(m)=t;精品文档交流endfori=k+1:n,%计算l(i,k)并将其放到a(i,k)中a(i,k)=a(i,k)/a(k,k);forj=k+1:na(i,j)=a(i,j)-a(i,k)*a(k,j);endb(i)=b(i)-a(i,k)*b(k);endendx=zeros(n,1);%回代过程x(n)=b(n)/a(n,n);fork=n-1:-1:1,x(k)=(b(k)-sum(a(k,k+1:n)*x(k+1:n)))/a(k,k);endend对压缩带状矩阵进行求解if ver=='202', %高斯消去法m=p+q+1;a=zeros(n,m);fori=1:1:nforj=1:1:ma(i,j)=A((i-1)*m+j);endendb=zeros(n,1);fori=1:1:nb(i)=A(n*m+i);%求b(i)endfork=1:1:(n-1)%开始消去ifa(k,(p+1))==0disp('错误！');break;endst1=n;if(k+p)<nst1=k+p;end精品文档交流fori=(k+1):1:st1a(i,(k+p-i+1))=a(i,(k+p-i+1))/a(k,(p+1));forj=(k+1):1:(k+q)a(i,j+p-i+1)=a(i,j+p-i+1)-a(i,k+p-i+1)*a(k,j+p-k+1);endb(i)=b(i)-a(i,k+p-i+1)*b(k);endendx=zeros(n,1);%回代x(n)=b(n)/a(n,p+1);sum=0;fork=(n-1):-1:1sum=b(k);st2=n;if(k+q)<nst2=k+q;endforj=(k+1):1:st2sum=sum-a(k,j+p-k+1)*x(j);endx(k)=sum/a(k,p+1);sum=0;endenddisp('方程组的的解为：')%输出解disp(x)运行结果：首先是对测试文件 dat20171/dat20172的求解结果，具体如下：精品文档交流经过测试矩阵的初步测试，利用程序对接下来的两个数据文件进行处理， 由于dat20173/dat20174解的个数较多，在这里只截取不分解作为示意，具体求解如下：精品文档交流本专业算例一底边绝热的正方形导热物体，内无热源，其余三边温度如下，求解该正方形内节点1、2、3、4的温度。40 3030T1T22015T3T410yx x y 绝热解：导热物体无内