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

计算方法B上机报告

姓名：

学号：

班级：

学院：

任课教师：

2017年12月29日

题目一：

1.1

题目内容

某通信公司在一次施工中，需要在水面宽度为20米的河沟底部沿直线走向铺设一条沟底光缆。在铺设光缆之前需要对沟底的地形进行初步探测，从而估计所需光缆的长度，为工程预算提供依据。已探测到一组等分点位置的深度数据(单位:

(1)(2)估算所需光缆长度的近似值，并作出铺设河底光缆的曲线图；

1.2 实现题目的思想及算法依据

首先在题目（1）中要实现的是数据的拟合，显然用到的是我们在第三章中数据近似的知识内容。多项式插值时，这里有21个数据点，则是一个20次的多项式，但是多项式插值随着数据点的增多，会导致误差也会随之增大，插值结果会出现龙格现象，所以不适用于该题目中点数较多的情况。为了避免结果出现大的误差，同时又希望尽可能多地使用所提供的数据点，提高数据点的有效使用率，这里选择分段插值方法进行数据拟合。分段插值又可分为分段线性插值、分段二次插值和三次样条插值。由于题目中所求光缆的现实意义，而前两者在节点处的光滑性较差，因此在这里选择使用三次样条插值。

根据课本SPLINEM 算法和TSS 算法，采用第三种真正的自然边界条件，在选定边界条件和选定插值点等距分布后，可以先将数据点的二阶差商求出并赋值给右端向量d ，再根据TSS 解法求解三对角线线性方程组从而解得M 值。求出M 后，对区间进行加密，计算200个点以便于绘图以及光缆长度计算。 对于问题（2），使用以下的公式:

20

=()L f x ds ⎰

20

(f x =⎰ 19

1

0(k k

k f x +==∑⎰

1.3 算法结构

1. For n i ,,2,1,0⋅⋅⋅=

1.1 i i M y ⇒

2. For 2,1=k

2.1 For k n n i ,,1, -=

2.1.1 i k i i i i M x x M M ⇒----)/()(1

3. 101h x x ⇒-

4. For 1-,,2,1n i =

4.1 11++⇒-i i i h x x

4.2 b a c c h h h i i i i i i ⇒⇒-⇒+++2;1;)/(11 4.3 i i d M ⇒+16

5. 0000;;c M d M d n n ⇒⇒⇒λ

n n n b a b ⇒⇒⇒2;;20μ

6. 1111,γμ⇒⇒d b

7. For m k ,,3,2 = ! 获取M 的矩阵元素个数，存入m

7.1 k k k l a ⇒-1/μ 7.2 k k k k c l b μ⇒⋅-1- 7.3 k k k k l d γγ⇒⋅-1- 8. m m m M ⇒μγ/

9. For 1,,2,1 --=m m k

9.1 k k k k k M M c ⇒⋅-+μγ/)(1

10. k ⇒1 ! 获取x 的元素个数存入s 11. For 1,,2,1-=s i

11.1 if i x x ≤~

then k i ⇒；break else k i ⇒+1

12. x

x x x x x h x x k k k k ˆ~;~;

11⇒-⇒-⇒--- y h x h M y x h M y x M x M k k k k k k ~/]ˆ)6

()6(6ˆ6[2

211331

⇒-+-++---

1.4 matlab 源程序

n=20; x=0:n;

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];

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(i); end

a(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); end end

l=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);

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);

end

p(n+1)=q(n+1)/u(n+1);

for k=n:-1:1

p(k)=(q(k)-r(k)*p(k+1))/u(k);

end

fprintf('三对角线性方程组的解为：');

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);

for i=1:n1+1

for j=1:n

if x1(i)>=x(j)&&x1(i)<=x(j+1) %利用书本111页算法EV ASPLINE求解拟合曲线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);

end

end

end

plot(x,-y,'x') %画出插值点

hold on

plot(x1,-s) %画出三次样条插值拟合曲线

hold on

title('三次样条插值法拟合电缆曲线');

xlabel('河流宽度/m');

ylabel('河流深度/m');

Length=0;

for i=1:n1

L=sqrt((x1(i+1)-x1(i))^2+(s(i+1)-s(i))^2); %计算电缆长度

Length=Length+L;

end

fprintf('电缆长度(m)=');

disp(Length);