西安交通大学计算方法B上机报告

计算方法上机报告

姓名：

学号：

班级：能动上课班级：

题目及求解：

一、对以下和式计算：

∑

∞

?

??

??+-+-+-+=0681581482184161n n n n S n ，要求：

① 若只需保留11个有效数字，该如何进行计算； ② 若要保留30个有效数字，则又将如何进行计算；

1 算法思想

（1）根据精度要求估计所加的项数，可以使用后验误差估计，通项为： 1421114

16818485861681

n n n

a n n n n n ε??=

---<< ?+++++??； （2）为了保证计算结果的准确性，写程序时，从后向前计算； （3）使用Matlab 时，可以使用以下函数控制位数： digits(位数)或vpa(变量，精度为数)

2 算法结构

;0=s

??

?

??+-+-+-+=

681581482184161n n n n t n ； for 0,1,2,,n i =??? if 10m t -≤

end;

for ,1,2,,0n i i i =--???

;s s t =+

3 Matlab 源程序

clear; %清除工作空间变量 clc; %清除命令窗口命令 m=input('请输入有效数字的位数m='); %输入有效数字的位数 s=0;

for n=0:50

t=(1/16^n)*(4/(8*n+1)-2/(8*n+4)-1/(8*n+5)-1/(8*n+6));

if t<=10^(-m) %判断通项与精度的关系break;

end

end;

fprintf('需要将n值加到n=%d\n',n-1); %需要将n值加到的数值

for i=n-1:-1:0

t=(1/16^i)*(4/(8*i+1)-2/(8*i+4)-1/(8*i+5)-1/(8*i+6));

s=s+t; %求和运算

end

s=vpa(s,m) %控制s的精度

4 结果与分析

若保留11位有效数字，则n=7，此时求解得：

s =3.1415926536；

若保留30位有效数字时，则n=22, 此时求解得：

s =3.8。

通过上面的实验结果可以看出，通过从后往前计算，这种算法很好的保证了计算结果要求保留的准确数字位数的要求。

二、某通信公司在一次施工中，需要在水面宽度为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所示，光缆从A

h。

点铺至B点，在某点处的深度为i

图2-1 海底光缆线的长度预测模型

计算光缆长度时，用如下公式：

20

()L f x ds =

?

20

0(f x =? 19

1

(k k

k f x +==∑?

=2 算法结构

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）获取M 的矩阵元素个数，存入m 8）For m k ,,3,2 =

8.1 k k k l a ?-1/μ

8.2 k k k k c l b μ??-1- 8.3 k k k k l d γγ??-1- 9）m m m M ?μγ/

10）For 1,,2,1 --=m m k

10.1 k k k k k M M c ??-+μγ/)(1 11）获取x 的元素个数存入s 12）k ?1

13） For 1,,2,1-=s i

13.1 if i x x ≤~

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

14）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

?-+-++---

3 Matlab 源程序

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