计算方法上机作业

计算方法上机作业

上机实习题目

1.某通信公司在一次施工中，需要在水面宽度为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)请用合适的曲线拟合所测数据点；

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

（1）算法思想

分段多项式是由一些在相互连接的区间上的不同多项式连接而成的一条连续曲线，其中三次样条插值方法是一种具有较好“光滑性”的分段插值方法。在本题中，假设所铺设的光缆足够柔软，在铺设过程中光缆触地走势光滑，紧贴地面，并且忽略水流对光缆的冲击。

计算光缆长度时，用如下公式： 200()Lfxds

20'20()1()fxfxdx

191'20()1()kkkfxfxdx

22()()xy

本题采取三次样条插值的方法，因为三次样条插值方法是一种具有较好“光滑性”的分段插值方法。根据提供的数据，只用x,y值，不包含导数值，因此采用第三类三次插值多项式进行插值编程。

设计算法如下：

1. For ni,,2,1,0

1.1 iiMy

2. For 2,1k

2.1 For knni,,1,

2.1.1 ikiiiiMxxMM)/()(1

3. 101hxx

4. For 1-,,2,1ni

4.1 11iiihxx

4.2 bacchhhiiiiii2;1;)/(11

4.3 iidM16

5. 0000;;cMdMdnn

nnnbab2;;20

6. 1111,db

7. 获取M的矩阵元素个数，存入m 8. For mk,,3,2

8.1 kkkla1/

8.2 kkkkclb1-

8.3 kkkkld1-

9. mmmM/

10. For 1,,2,1mmk

10.1 kkkkkMMc/)(1

11. 获取x的元素个数存入s

12. k1

13. For 1,,2,1si

13.1 if ixx~ then ki；break

else ki1

14. xxxxxxhxxkkkkˆ~;~;11

yhxhMyxhMyxMxMkkkkkk~/]ˆ)6()6(6ˆ6[2211331

（3）源程序

clear;

clc;

x=0:1:20; %产生从0到20含21个等分点的数组

X=0:0.2: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.80,10.93];

%等分点位置的深度数据

n=length(x); %等分点的数目

N=length(X);

%% 求三次样条插值函数s(x) M=y;

for k=2:3; %计算二阶差商并存放在M中

for i=n:-1:k;

M(i)=(M(i)-M(i-1))/(x(i)-x(i-k+1));

end

end

h(1)=x(2)-x(1); %计算三对角阵系数a,b,c及右端向量d

for i=2:n-1;

h(i)=x(i+1)-x(i);

c(i)=h(i)/(h(i)+h(i-1));

a(i)=1-c(i);

b(i)=2;

d(i)=6*M(i+1);

end

M(1)=0; %选择自然边界条件

M(n)=0;

b(1)=2;

b(n)=2;

c(1)=0;

a(n)=0;

d(1)=0;

d(n)=0;

u(1)=b(1); %对三对角阵进行LU分解

y1(1)=d(1);

for k=2:n;

l(k)=a(k)/u(k-1);

u(k)=b(k)-l(k)*c(k-1);

y1(k)=d(k)-l(k)*y1(k-1);

end

M(n)=y1(n)/u(n); %追赶法求解样条参数M(i)

for k=n-1:-1:1;

M(k)=(y1(k)-c(k)*M(k+1))/u(k);

end s=zeros(1,N);

for m=1:N;

k=1;

for i=2:n-1

if X(m)<=x(i);

k=i-1;

break;

else

k=i;

end

end

H=x(k+1)-x(k); %在各区间用三次样条插值函数计算X点处的值

x1=x(k+1)-X(m);

x2=X(m)-x(k); s(m)=(M(k)*(x1^3)/6+M(k+1)*(x2^3)/6+(y(k)-(M(k)*(H^2)/6))*x1+(y(k+1)-(M(k+1)*(H^2)/6))*x2)/H;

end

%% 计算所需光缆长度

L=0; %计算所需光缆长度

for i=2:N

L=L+sqrt((X(i)-X(i-1))^2+(s(i)-s(i-1))^2);

end

disp('所需光缆长度为 L=');

disp(L);

figure

plot(x,y,'*',X,s,'-') %绘制铺设河底光缆的曲线图

xlabel('位置X/测量点','fontsize',16); %标注坐标轴含义

ylabel('深度Y/m','fontsize',16);

title('铺设河底光缆的曲线图','fontsize',16);

grid;

（4）结果与分析

铺设海底光缆的曲线图如下图所示：

拟合结果表明，运用分段三次样条插值所得的拟合曲线能较准确地反映铺设光缆的走势图。 计算出所需光缆的长度为 L=26.4844m。可以用Newton法进行拟合，求得拟合曲线和光缆长度，也可以用三次样条法来拟合，则精度会更高，根据实际光缆的现实特性，三次样条法来拟合。

2.假定某天的气温变化记录如下表所示，试用数据拟合的方法找出这一天的气温变化的规律;试计算这一天的平均气温,并试估计误差。

时刻 0 1 2 3 4 5 6 7 8 9 10 11 12

平均气温 15 14 14 14 14 15 16 18 20 20 23 25 28

时刻 13 14 15 16 17 18 19 20 21 22 23 24

平均气温 31 34 31 29 27 25 24 22 20 18 17 16

（1）解题思想和数学原理：

对于具体实验时，通常不是先给出函数的解析式，再进行实验，而是通过实验的观察和测量给出离散的一些点，再来求出具体的函数解析式。又因为测量误差的存在，实际真实的解析式曲线并不一定通过测量给出的所有点。最小二乘法，形成法方程是求解这一问题的很好的方法，本实验运用这一方法实现对给定数据的拟合。

（2）matlab源程序：

x=[0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24]; %给出题目数据（时间）

y=[15 14 14 14 14 15 16 18 20 20 23 25 28 31 32 31 29 27 25 24 22 20 18 17

16]; %给出题目数据（温度） plot(x,y, 'm*') %画出各个离散数据点

hold on

for n=2:4; %2、3、4代表拟合函数最高项次数

alltemp=25; % alltemp代表数据点总共有25个

A=zeros(n+1,n+1); %定义初始正规方程组的系数矩阵A

C=ones(n+1,1); %定义初始正规方程组的系数向量C

D=zeros(n+1,1); %定义初始正规方程组的向量D

for i=1:n+1

for j=1:n+1

for k1=1: alltemp

A(i,j)=A(i,j)+(x(k1).^(i-1+j-1));

end

end

for k2=1: alltemp

D(i,1)=D(i,1)+(x(k2).^(i-1)).*(y(k2));

end

end %以上为计算出正规方程组矩阵A、D的所有元素的程序

tol=1.0e-12;

maxit=1000;

C=bicg(A,D,tol,maxit); %使用bicg迭代算出正规方程组的系数向量C

p=0; %误差分量

E=0; %误差总量

if n==2

b=0:24;

f=C(1)+C(2).*b+C(3).*(b.^2);

p=y(b+1)-f;

for v=1:25

E=E+(p(v)).^2;

end

plot(b,f, 'r-')