西安交大 计算方法B上机作业.docxVIP

一、三次样条拟合 计算方法（B）上机作业 某通信公司在一次施工中，需要在水面宽度为 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、算法实现的思想及依据 题目（1）为曲线拟合问题多项式插值、分段插值和最小二乘法。多项式插值，随着插值数据点的数目增多，误差也会随之增大，因此不选用。最小二乘法适于数据点较多的场合，在此也不适用。故选用分段插值。 分段插值又分为分段线性插值、分段二次插值、三次样条插值及更高阶的多项式插值。由本题的物理背景知，光缆正常工作时各点应该是平滑过渡，因此至少选用三次样条插值法。对于更高阶的多项式插值，由于“龙格现象”而不选用。 题目（2）求光缆长度，即求拟合曲线在0 到 20 的长度，对弧长进行积分即 可。光缆长度的第一型线积分表达式为l ? ?19 k ?0 ?k ?1 k 1? f (x)2 dx 。 2、算法实现的结构 参考教材给出的SPLINEM算法和TTS算法，在选定边界条件和选定插值点等距分布后，可以先将数据点的二阶差商求出来并赋值给右端向量d，再根据TSS解法 k求解M。光缆长度的第一型线积分表达式为l ? ?19 ?k ?1 1? f (x)2 dx 。 k 3、程序运行结果及分析 k ?0 图 1.1 三种边界条件下三次样条插值 图 1.2 光缆长度 4、MATLAB 代码： 1）自己编程实现代码 clear; clc; I=input(你想使用第几种边界条件？请输入1、2、3 之一: ); 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.8 10.93]; plot(x,-y,k.,markersize,15)%y为深度，取负号hold on %% 计算一阶差商 y1=ones(1,21); for i=2:1:21 y1(i)=(y(i)-y(i-1))/(x(i)-x(i-1)); end %% 计算二阶差商 y2=ones(1,21); for i=3:1:21 y2(i)=(y1(i)-y1(i-1))/(x(i)-x(i-2)); end %% 计算三阶差商 y3=ones(1,21); for i=4:1:21 y3(i)=(y2(i)-y2(i-1))/(x(i)-x(i-3)); end %% 选择边界条件（I） if I==1 d(1)=0;d(21)=0;a(21)=0;c(1)=0;% 第一个点和最后一个点的二阶差商为0 end if I==2 d(1)=6*y1(1); d(21)=-6*y1(21); a(1)=1;c(1)=1; end if I==3 d(1)=-12*y3(1); d(21)=-12*y3(21); a(21)=-2;c(1)=-2;% end for i=2:20 d(i)=6*y2(i+1); end %% 构造带状矩阵求解（追赶法） b=2*ones(1,21); a=0.5*ones(1,21);%a(21)=-2; c=0.5*ones(1,21);%c(1)=-2; u(1)=b(1);r(1)=c(1); %% 追 yz(1)=d(1); for i=2:21 l(i)=a(i)/u(i-1); u(i)=b(i)-l(i)*r(i-1); r(i)=c(i); yz(i)=d(i)-l(i)*yz(i-1); end %% 赶xg(21)=yz(21)/u(21); for i=20:-1:1 xg(i)=(yz(i)-r(i)*xg(i+1))/u(i); end M=xg;%%所有点的二阶导数值 %% 求函数表达式并积分 t=1; h=1; N=1000 x1=0:20/(N-1):20; length=0; for i=1:N for j=2:20 if x1(i)=x(j) t=j; break; end els