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

计算方法B上机报告姓名：学号：班级：学院：任课教师：2017年12月29日题目一：1.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)估算所需光缆长度的近似值，并作出铺设河底光缆的曲线图；1.2实现题目的思想及算法依据首先在题目（1）中要实现的是数据的拟合，显然用到的是我们在第三章中数据近似的知识内容。多项式插值时，这里有21个数据点，则是一个20次的多项式，但是多项式插值随着数据点的增多，会导致误差也会随之增大，插值结果会出现龙格现象，所以不适用于该题目中点数较多的情况。为了避免结果出现大的误差，同时又希望尽可能多地使用所提供的数据点，提高数据点的有效使用率，这里选择分段插值方法进行数据拟合。分段插值又可分为分段线性插值、分段二次插值和三次样条插值。由于题目中所求光缆的现实意义，而前两者在节点处的光滑性较差，因此在这里选择使用三次样条插值。根据课本SPLINEM算法和TSS算法，采用第三种真正的自然边界条件，在选定边界条件和选定插值点等距分布后，可以先将数据点的二阶差商求出并赋值给右端向量d，再根据TSS解法求解三对角线线性方程组从而解得M值。求出M后，对区间进行加密，计算200个点以便于绘图以及光缆长度计算。对于问题（2），使用以下的公式:200=()Lfxds20'20()1()fxfxdx191'20()1()kkkfxfxdx1.3算法结构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;)/(114.3iidM165.0000;;cMdMdnnnnnbab2;;206.1111,db7.Formk,,3,2!获取M的矩阵元素个数，存入m7.1kkkla1/7.2kkkkclb1-7.3kkkkld1-8.mmmM/9.For1,,2,1mmk9.1kkkkkMMc/)(110.k1!获取x的元素个数存入s11.For1,,2,1si11.1ifixx~thenki；breakelseki112.xxxxxxhxxkkkkˆ~;~;11yhxhMyxhMyxMxMkkkkkk~/]ˆ)6()6(6ˆ6[22113311.4matlab源程序n=20;x=0:n;y=[9.018.967.967.978.029.0510.1311.1812.2613.2813.3212.6111.2910.229.157.907.958.869.8110.8010.93];M=y;%用于存放差商，此时为零阶差商h=zeros(1,n+1);c=zeros(1,n+1);d=zeros(1,n+1);a=zeros(1,n+1);b=2*ones(1,n+1);h(2)=x(2)-x(1);fori=2:n%书本110页算法SPLINEMh(i+1)=x(i+1)-x(i);c(i)=h(i+1)/(h(i)+h(i+1));a(i)=1-c(i);enda(n+1)=-2;%计算边界条件c(0),a(n+1),采用的是第三类边界条件c(1)=-2;fork=1:3%计算k阶差商fori=n+1:-1:k+1M(i)=(M(i)-M(i-1))/(x(i)-x(i-k));endif(k==2)%计算2阶差商d(2:n)=6*M(3:n+1);%给d赋值endif(k==3)d(1)=(-12)*h(2)*M(4);%计算边界条件d(0),d(n),采用的是第三类边界条件d(n+1)=12*h(n+1)*M(n+1);endendl=zeros(1,n+1);r=zeros(1,n+1);u=zeros(1,n+1);q=zeros(1,n+1);u(1)=b(1);r(1)=c(1);q(1)=d(1);fork=2:n+1%利用书本49页算法TSS求解三对角线性方程组r(k)=c(k);l(k)=a(k)/u(k-1);u(k)=b(k)-l(k)*r(k-1);q(k)=d(k)-l(k)*q(k-1);endp(n+1)=q(n+1)/u(n+1);fork=n:-1:1p(k)=(q(k)-r(k)*p(k+1))/u(k);endfprintf('三对角线性方程组的解为：');disp(p);%求拟合曲线x1=0:0.1:20;%首先对区间进行加密,增加插值点n1=10*n;x2=zeros(1,n1+1);x3=zeros(1,n1+1);s=zeros(1,n1+1);fori=1:n1+1forj=1:nifx1(i)=x(j)&&x1(i)=x(j+1)%利用书本111页算法EVASPLINE求解拟合曲线s(x)h(j+1)=x(j+1)-x(j);x2(i)=x(j+1)-x1(i);x3(i)=x1(i)-x(j);s(i)=(p(j).*(x2(i)).^3/6+p(j+1).*(x3(i)).^3/6+(y(j)-p(j).*((h(j+1)).^2/6)).*x2(i)+...(y(j+1)-p(j+1).*(h(j+1)).^2/6).*x3(i))/h(j+1);endendendplot(x,-y,'x')%画出插值点holdonplot(x1,-s)%画出三次样条插值拟合曲线holdontitle('三次样条插值法拟合电缆曲线');xlabel('河流宽度/m');ylabel('河流深度/m');Length=0;fori=1:n1L=sqrt((x1(i+1)-x1(i))^2+(s(i+1)-s(i))^2);%计算电缆长度Length=Length+L;endfprintf('电缆长度(m)=');disp(Length);1.5结果与说明铺设海底光缆的曲线如图1.1所示由上图可以看出，所得到的曲线光滑，能够较好得反映实际的河沟底部地势形貌。电缆长度计算结果为26.6656m（图1.2）。图1.1三次样条插值法拟合海底光缆曲线图1.2海底光缆长度结果题目二2.1题目内容：已知非线性方程0)sincos(10dx在[2,3]中有根。请设计算法，求出该根，并使求出的根的误差不超过410。2.2实现题目的思想及算法依据对于该题的非线性方程，可以将其分解成两个部分：（1）求解数值积分；（2）求解非线性方程。首先求解数值积分，令()cos(sin)gx，则利用最简单的梯形公式可以得到11()[(()()]2niiihfxgg，其中,ihihn。于是就有了()0fx形式的非线性方程，这里选择二分法进行求解。算法参考课本BISECTION算法。2.3算法结构1.00fxf11fxf2.if010ffthenstop3.if02fthen输出0x作为根；stop4.if12fthen输出1x作为根；stop5.0112x+xx6.If111xxxthen输出x作为根；stop7.fxf8.if2fthen输出x作为根；9.if10ffthen9.10xx;0ffElse9.21xx;1ff10.goto52.4Matlab源程序************************functiony=theta(i)y=i*h;end*************************functiony=g(x,theta)%为了计算关于theta的数值积分，先令g=cos(x*sin(theta))y=cos(x*sin(theta));end**************************functionf=hsz(x)%计算数值积分n=10000;h=pi/n;%将区间分成n份f=0;fori=1:nf=f+h/(2*pi)*(g(x,i+1)+g(x,i));endend**************************error=1e-4;%误差允许值a=2;b=3;%初始区间f0=hsz(a);f1=hsz(b);iff0*f10%判断方程是否有解disp('该方程在[a,b]上无解');elseiff0==0x=a;elseiff1==0x=b;%判断方程解是否在区间两边界上else%二分法求解方程得解a0=a;b0=b;whileabs((b0-a0)/2)=errorhalf=(a0+b0)/2;fa=hsz(a0);fb=hsz(b0);fhalf=hsz(half);%计算中点处的函数值，用以判断解的位置iffhalf==0x=half;break;elseiffa*fhalf0b0=half;%定义新区间，为原区间的一半elsea0=half;%定义新区间endendx=(b0+a0)/2;%方程组的解endfprintf('方程组的解为：')disp(x)2.5结果与说明由于是利用梯形公式来求解，需要将区间划分为n个区间，而n的取值是很关键的，需要取得适当的n值才能满足误差精度。这里将区间分成10000份，可以得到结果图2.1，x=2.4048。图2.1n=10000方程的解由于变量的嵌套很多及自身水平的不足，所以我将很多变量通过调用子程序的方式实现，分别为theta.m:g.m;hsz.m。题目三3.1题目内容：编写程序实现大规模方程组的高斯消去法程序，并对所附的方程组进行求解。所附方程组的类型为对角占优的带状方程组。数据说明：(1)dat20171.dat为非压缩带状方程组，阶数为10阶，该方程组供调试程序使用,该方程组的根都为1;(2)dat20172.dat为压缩带状方程组，阶数为20阶，该方程组供调试程序使用,该方程组的根都为1;(2)dat20173.dat为非压缩带状方程组，阶数在3000阶左右;(3)dat20174.dat为压缩带状方程组，阶数在50000阶左右;3.2实现题目的思想及算法依据高斯消去法主要分为消去和回代过程，消去过程形成上三角矩阵，回代过程就是自下而上实现方程组求解。题目中系数矩阵是严格对角占优矩阵，所以无需进行列主元，极大的提高了程序的效率。对于n阶、上带宽q、下带宽p的非压缩格式带状矩阵，通过比较k+p、k+q及n的大小，对于行消去，选择k+p和n中的最小值作为循环终点。在用第k行进行消去时，只需要计算第k列的k+1ka，到k+pka，，每一行只需计算kk+1a，到kk+qa，；从第n-c+1到n行，第k列计算k+1ka，到nka，，每一行只需计算kk+1a，到kna，，于是减小运算量。压缩格式忽略了零元素，存储方式与非压缩格式有所不同，所以需要建立两种格式中元素之间的对应关系，即压缩kp+1A，=非压缩kkA，，根据此关系进行程序编写。3.3算法结构1.A的阶数n2．For12ik,k,,n…2.1ikkkik2.2For12jk,k,,n…2.2.1ijkjikij2.3ikiki-nnnnx/3.For121kn,n,,…