《西安交大计算方法B2017大作业》由会员分享，可在线阅读，更多相关《西安交大计算方法B2017大作业（26页珍藏版）》请在金锄头文库上搜索。
1、word计算方法B上机报告某某：学号：班 级：学院：任课教师：2017年12月29日题目一：1.1 题目内容某通信公司在一次施工中，需要在水面宽度为20米的河沟底部沿直线走向铺设一条沟底光缆。在铺设光缆之前需要对沟底的地形进展初步探测，从而估计所需光缆的长度，为工程预算提供依据。已探测到一组等分点位置的深度数据(单位:米)如下表所示：分点0123456深度分点78910111213深度分点14151617181920深度(1)请用适宜的曲线拟合所测数据点；(2)估算所需光缆长度的近似值，并作出铺设河底光缆的曲线图；1.2 实现题目的思想与算法依据 首先在题目1中要实现的是数据的拟合，显然用到的
2、是我们在第三章中数据近似的知识内容。多项式插值时，这里有21个数据点，如此是一个20次的多项式，但是多项式插值随着数据点的增多，会导致误差也会随之增大，插值结果会出现龙格现象，所以不适用于该题目中点数较多的情况。为了防止结果出现大的误差，同时又希望尽可能多地使用所提供的数据点，提高数据点的有效使用率，这里选择分段插值方法进展数据拟合。分段插值又可分为分段线性插值、分段二次插值和三次样条插值。由于题目中所求光缆的现实意义，而前两者在节点处的光滑性较差，因此在这里选择使用三次样条插值。 根据课本SPLINEM算法和TSS算法，采用第三种真正的自然边界条件，在选定边界条件和选定插值点等距分布后，可以
3、先将数据点的二阶差商求出并赋值给右端向量d，再根据TSS解法求解三对角线线性方程组从而解得M值。求出M后，对区间进展加密，计算200个点以便于绘图以与光缆长度计算。 对于问题2，使用以下的公式:1.3 算法结构1.For 1.1 2. For 2.1 For 2.1.1 3. 4. For 4.1 4.2 4.3 5. 6. 7. For ! 获取M的矩阵元素个数，存入m7.1 7.2 7.3 8. 9. For 9.1 10. ! 获取x的元素个数存入s11. For 11.1 if then ；break else 12. 1.4 matlab源程序n=20;x=0:n;y=9.01 8.
4、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;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);for i=2:n %书本110页算法SPLINEM h(i+1)=x(i+1)-x(i); c(i)=h(i+1)/(h(i)+h(i+1); a(i)=1-c
5、(i);enda(n+1)=-2; %计算边界条件c(0),a(n+1),采用的是第三类边界条件c(1)=-2;for k=1:3 %计算k阶差商 for i=n+1:-1:k+1 M(i)=(M(i)-M(i-1)/(x(i)-x(i-k); end if(k=2) %计算2阶差商 d(2:n)=6*M(3:n+1); %给d赋值 end if(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
6、(1,n+1);q=zeros(1,n+1);u(1)=b(1);r(1)=c(1);q(1)=d(1);for k=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);for k=n:-1:1 p(k)=(q(k)-r(k)*p(k+1)/u(k);endfprintf(三对角线性方程组的解为：);disp(p);%求拟合曲线x1=0:0.1:20; %首先对区间进展加密,增加插值点n1=1
7、0*n;x2=zeros(1,n1+1);x3=zeros(1,n1+1);s=zeros(1,n1+1);for i=1:n1+1 for j=1:n if x1(i)=x(j)&x1(i)0 %判断方程是否有解 disp(该方程在a,b上无解); elseif f0=0 x=a;elseif f1=0 x=b; %判断方程解是否在区间两边界上else %二分法求解方程得解 a0=a; b0=b; while abs(b0-a0)/2)=error half=(a0+b0)/2; fa=hsz(a0); fb=hsz(b0); fhalf=hsz(half); %计算中点处的函数值，用以判断
8、解的位置 if fhalf=0 x=half; break; elseif fa*fhalf0 b0=half; %定义新区间，为原区间的一半 else a0=half; %定义新区间 end end x=(b0+a0)/2; %方程组的解endfprintf(方程组的解为：)disp(x)2.5 结果与说明由于是利用梯形公式来求解，需要将区间划分为n个区间，而n的取值是很关键的，需要取得适当的n值才能满足误差精度。这里将区间分成10000份，可以得到结果图2.1，x=2.4048。图2. 1 n=10000方程的解由于变量的嵌套很多与自身水平的不足，所以我将很多变量通过调用子程序的方式实现，分别为theta.m : g.m; hsz.m。
