计算方法上机报告姚威资料.docx
《计算方法上机报告姚威资料.docx》由会员分享，可在线阅读，更多相关《计算方法上机报告姚威资料.docx（59页珍藏版）》请在冰豆网上搜索。
计算方法上机报告姚威资料
计算方法B上机作业报告
姓名：
姚威
班级：
硕5081班
学号：
3115312022
学科专业：
电气工程
任课教师：
苏剑老师
2015年12月
1.上机题一
对以下和式计算：
，要求：
（1）若只需保留11个有效数字，该如何进行计算；
（2）若要保留30个有效数字，则又将如何进行计算；
计算方法
任何一个β进制、具有t位有效数字的实数
总可以表示成
，其中，
是满足
的整数。
有效数字是指一个近似数的“有意义”的数字的数位，通常在十进制中讨论，设
其中
，近似数
，若
，则称
是
的具有t位有效数字的近似数，或称
具有t位有效数字。
根据有效数字的要求可以得到计算公式的精度要求。
再利用后验误差估计式可以根据精度要求得到和式累加的项数。
首先满足题目要求的情况下，利用后验误差估计计算出算法中的最大运行次数。
即以一下条件作为迭代终止的准则：
，其中
为误差数值。
该实验第一问利用后验误差估计截取的项数，使误差小于1e-10；第二问利用后验误差估计截取的项数，使误差小于1e-30。
得出运行次数，作为累加的终止条件。
特别地，为了减小舍入误差在计算S时所采用的方法是逆序相加，其依据是：
S中各项的绝对值为递减的正数，而两个数量级相差较大的数字相加减时，较小数的有效数字会丢失，从而导致最后的运算结果失真。
为避免“大数吃小数”现象的发生，采用逆序相加；同时为了避免相近数相减和减小计算步骤，我们将
化简为
进行计算。
另外，要限定计算结果的有效数字的个数，在MATLAB软件中可直接调用函数vpa（变量名,位数）来控制计算结果的位数，最终完成精度要求较高的计算。
解题源程序
functionsolve1
clc;
clear;
sum1=0;
sum2=0;
n=0;
N=0;
%第一小问
%根据精度要求利用后验误差式，估计要使误差小于
，需要累加的项数n，所以要保证计算精度第n+1项要小于
，才能保证计算的精度
whileabs（（960*n^2+1208*n+376）/（（8*n+1）*（8*n+4）*（8*n+5）*（8*n+6））/（16.^n））>10^（-11）
%为了减小误差，我们对计算公式进行变形，一是防止相近数相减，二是为了减少计算步骤；
n=n+1;
end
n;
%累加进行的项数
disp（'累加的次数:
'）
n
%为减少舍入误差的影响，按绝对值递增的顺序求和。
因此按n按从小到大的顺序相加,步距取为-1
fori=n:
-1:
0
sum1=sum1+（（960*i.^2+1208*i+376）/（（8*i+1）*（8*i+4）*（8*i+5）*（8*i+6）））/（16.^i）;
end
disp（'保留11位有效数字时的近似值:
'）
Sum1=vpa（sum1,11）
%第二小问
%根据精度要求利用后验误差式，估计要使误差小于0.5e-30需要累加的项数N
While
abs（（960*N^2+1208*N+376）/（（8*N+1）*（8*N+4）*（8*N+5）*（8*N+6））/（16.^N））>10^（-30）
N=N+1;
end
N;
%累加进行的项数
disp（'累加的次数:
'）
N
%为减少舍入误差的影响，按绝对值递增的顺序求和。
因此按N按从小到大的顺序相加,步距取为-1
forj=N:
-1:
0
sum2=sum2+（（960*j.^2+1208*j+376）/（（8*j+1）*（8*j+4）*（8*j+5）*（8*j+6）））/（16.^j）;
end
disp（'保留30位有效数字时的近似值:
'）
Sum2=vpa（sum2,30）
End
计算结果
累加的次数:
n=8
保留11位有效数字时的近似值:
sum1=3.1415926536
累加的次数:
N=23
保留30位有效数字时的近似值:
sum2=3.14159265358979323846264338328
实验结果分析及总结
通过运行上述程序可知：
本题通过计算可以发现，计算值和
的值非常接近，由
，可知这里的
当t=10时，
，我们为了得到有效数位为11，所以计算的最后一项t应该等于t=-11。
同理，为了得到有效数位为30，所以计算的最后一项t应该等于t=-30。
第一问中，在满足误差小于1e-10时，截取的项数为8，，保留11位有效数字估算值:
3.1415926536
。
第二问中，在满足误差小于1e-30时，截取的项数为23，保留30位有效数字估算值:
3.14159265358979323846264338328。
从上述的算法思想中可以看出，浮点数运算中不仅要满足误差要求，还要尽可能地减少计算量，而其计算量的大小很大程度上是由算式要求的精度来决定的；此外还要考虑舍入误差的影响，这时就要对所运算数据的性质进行分析，设置合适的算法，从而提高运算的精度。
2.上机题二
某通信公司在一次施工中，需要在水面宽度为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）预测所需光缆长度的近似值，并作出铺设河底光缆的曲线图；
三次样条插值拟合数据点
2.1.1算法组织
使用分段插值多项式拟合数据的曲线没有足够的光滑性，所以此题采用三次样条函数拟合所测数据点，即在每一个子区间
上是一个三次多项式，这些多项式的曲线在节点
处“光滑”地连接起来。
这样的分段三次多项式称为三次样条函数。
样条函数为：
记
则上式可写成：
这是一个有n-1个方程，n+1个未知数
的线性方程组，还需要两个边界条件才能确定唯一解。
（1）自然样条
（2）指定端点的一阶导数
。
此时得
（3）若两端点的导数值均不知道时，可以通过推导得到
最后可以得到这些方程组的统一形式
本题采用边界条件三求解，其中
2.1.2算法
SPLINEM（
）
[本算法计算三次样条的参数Mi，存放于数组M中，参数
根据不同边界条件给定。
]
1.1
2.1
2.1.1
3.
4.
4.1
4.2
4.3
用算法TSS（a,b,c,M,n,M）
算法TTS（a，b，c，d，n，x）
[本算法解对角线性方程组，系数存于数组a、b和c中，d是有段向量。
]
2.1.3算法的Matlab实现
三次样条插值及复化Simpson公式主程序
clc
clear
%%构造数据矩阵B
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（y）;
x=0:
20/（N-1）:
20;
Data=zeros（2,N）;
fori=1:
N
Data（1,i）=x（i）;
Data（2,i）=y（i）;
end
Data;
%构造步长矩阵h
h=zeros（1,N-1）;
fork=1:
N-1
h（k）=x（k+1）-x（k）;
end
%利用差分表计算d1，dN及di
x=Data（1,:
）;
y=Data（2,:
）;
fori=1:
N
D（i,1）=x（i）;
D（i,2）=y（i）;
end
fori=2:
N
D（i,3）=（y（i）-y（i-1））/（x（i）-x（i-1））;
end
fori=3:
N
D（i,4）=（D（i,3）-D（i-1,3））/（x（i）-x（i-2））;
end
fori=4:
N
D（i,5）=（D（i,4）-D（i-1,4））/（x（i）-x（i-3））;
end
D;
%构造M的系数矩阵A
A=zeros（N）;
%由第三类边界条件可得
A（1,2）=-2;
%由第三类边界条件可得
A（N,N-1）=-2;
fori=1:
N
A（i,i）=2;
end
fori=2:
N-1
A（i,i+1）=h（i）/（h（i-1）+h（i））;
A（i,i-1）=h（i-1）/（h（i-1）+h（i））;
end
%构造M的结果矩阵d
d=zeros（N,1）;
d
（1）=-12*h
（1）*D（4,5）;
d（N）=12*h（N-1）*D（N,5）;
forj=2:
N-1
d（j）=6*D（j,4）;
end
%利用TSS方法解方程组
M=TSS（A,d）;
%得到各段三次样条插值的函数表达式，并作出曲线图
y=Data（2,:
）;
[P,N]=size（Data）;
fori=1:
N-1
xx=x;
Sy=getS（i,xx,M,h,y）;
sx=xx（i）:
0.001:
xx（i+1）;
sy=Sy（sx）;
str=[repmat（'（',21,1）num2str（x'）repmat（',',21,1）num2str（y'）];
plot（x,y,'o'）%原探测数据点及数值
text（x,-y,cellstr（str））
stem（x,-y,'filled','LineWidth',3）%花火柴杆图
plot（sx,sy,'r','LineWidth',3）%按插值函数表达式作图
holdon
end
%添加水平线
l=line（[0,20],[0,0]）;
%设置直线的宽度和颜色
set（l,'Linewidth',1.5）
set（l,'color','k'）
%设置坐标轴刻度
set（gca,'Xtick',[0:
1:
20]）
%设置坐标轴刻度
set（gca,'Ytick',[-15:
2:
15]）
%添加图形标题和坐标轴名称
title（'作出的铺设河底光缆的曲线图'）
grid
xlabel（'探测的一组等分点位置'）
ylabel（'相应探测点对应的深度数据'）
%利用复化Simpson公式计算光缆的长度
y=Data（2,:
）;
[P,N]=size（Data）;
sl=0;
fori=1:
N-1
Sy=getS（i,xx,M,h,y）;
symst
dsy=diff（Sy（t）,t）;
g=sqrt（1+dsy^2）;
u=i-1;
v=i;
sl=sl+simpson（g,u,v,20）;
end
Sy;
g;
disp（'所需光缆长度的近似值为:
'）;
sl
TSS解方程组的子程序
%%追赶法解方程组
functionx=TSS（A,d）
[M,N]=size（A）;
b（N）=A（N,N）;
fori=1:
N-1
a（i+1）=A（i+1,i）;
b（i）=A（i,i）;
c（i）=A（i,i+1）;
end
u
（1）=b
（1）;
y
（1）=d
（1）;
fori=2:
N
l（i）=a（i）/u（i-1）;
u（i）=b（i）-l（i）*c（i-1）;
y（i）=d（i）-l（i）*y（i-1）;
end
x（N）=y（N）/u（N）;
fori=N-1:
-1:
1
x（i）=（y（i）-c（i）*x（i+1））/u（i）;
end
x;
end
得到子区间函数表达式的子程序
functionSy=getS（k,xx,m,h,y）
symsSxSy
i=k+1;
S=（xx（i）-x）^3*m（i-1）/（6*h（i-1））+（x-xx（i-1））^3*m（i）/（6*h（i-1））+（y（i-1）-h（i-1）^2*m（i-1）/6）*（xx（i）-x）/h（i-1）+（y（i）-h（i-1）^2*m（i）/6）*（x-xx（i-1））/h（i-1）
Sy=inline（（xx（i）-x）^3*m（i-1）/（6*h（i-1））+（x-xx（i-1））^3*m（i）/（6*h（i-1））+（y（i-1）-h（i-1）^2*m（i-1）/6）*（xx（i）-x）/h（i-1）+（y（i）-h（i-1）^2*m（i）/6）*（x-xx（i-1））/h（i-1））;
end
2.1.4实验结果及分析
由上面程序可以得到如下结果：
在
时的函数表达式为
S=（845002892516739*（x-1）^3）/4503599627370496-（3313195038365*x）/8796093022208+（7514075829091499*x^3）/54043195528445952+5177804441890613/562949953421312
在
时的函数表达式为
S=2621488579161321/281474976710656-（2489533575770641*（x-1）^3）/6755399441055744-（7514075829091499*（x-2）^3）/54043195528445952-（277217198887389*x）/562949953421312
在
时的函数表达式为
S=（650311338860789*（x-2）^3）/1688849860263936-（418602078066937*x）/562949953421312+（2489533575770641*（x-3）^3）/6755399441055744+2762873458340869/281474976710656
在
时的函数表达式为
S=（2624680819031*x）/4398046511104-（1092494410535677*（x-3）^3）/6755399441055744-（650311338860789*（x-4）^3）/1688849860263936+3262063247973023/562949953421312
在
时的函数表达式为
S=（79721223946041*x）/140737488355328+（8155793057367127*（x-4）^3）/27021597764222976+（1092494410535677*（x-5）^3）/6755399441055744+3330360244180239/562949953421312
在
时的函数表达式为
S=（814815571273233*x）/562949953421312-（3544057556774495*（x-5）^3）/54043195528445952-（8155793057367127*（x-6）^3）/27021597764222976+425353433367447/281474976710656
在
时的函数表达式为
S=（274137988850561*x）/281474976710656+（566803888785965*（x-6）^3）/54043195528445952+（3544057556774495*（x-7）^3）/54043195528445952+306243053520945/70368744177664
在
时的函数表达式为
S=（617478217955525*x）/562949953421312-（137********91219*（x-7）^3）/216172782113783808-（566803888785965*（x-8）^3）/54043195528445952+1965528746386739/562949953421312
在
时的函数表达式为
（545284355633681*x）/562949953421312+（405384572326449*（x-8）^3）/9007199254740992+（137********91219*（x-9）^3）/216172782113783808+2543079644961491/562949953421312
在
时的函数表达式为
S=（179389613108493*x）/562949953421312-（3156841900829689*（x-9）^3）/135********111488-（405384572326449*（x-10）^3）/9007199254740992+5836132327688183/562949953421312
在
时的函数表达式为
S=（3156841900829689*（x-11）^3）/135********111488-（814194773093437*（x-10）^3）/9007199254740992-（480342372812029*x）/562949953421312+12433452186893403/562949953421312
在
时的函数表达式为
S=（814194773093437*（x-12）^3）/9007199254740992-（139********62219*（x-11）^3）/9007199254740992-（706852419243083*x）/562949953421312+14925062697634997/562949953421312
在
时的函数表达式为
S=（5376237321301877*（x-12）^3）/54043195528445952-（745487614849087*x）/562949953421312+（139********62219*（x-13）^3）/9007199254740992+153********907045/562949953421312
在
时的函数表达式为
S=（493605447569857*（x-13）^3）/72057594037927936-（550210270623049*x）/562949953421312-（5376237321301877*（x-14）^3）/54043195528445952+12850079569968551/562949953421312
在
时的函数表达式为
S=139********622509/562949953421312-（857131708001443*（x-14）^3）/6755399441055744-（493605447569857*（x-15）^3）/72057594037927936-（314201753442023*x）/281474976710656
在
时的函数表达式为
S=（2888372562541421*（x-15）^3）/9007199254740992-（447606859642455*x）/1125899906842624+（857131708001443*（x-16）^3）/6755399441055744+157********360463/1125899906842624
在
时的函数表达式为
S=（122327*********7*x）/1125899906842624+（3896133174998587*（x-16）^3）/27021597764222976-（2888372562541421*（x-17）^3）/9007199254740992-10982567962964529/1125899906842624
在
时的函数表达式为
S=（637035986695777*x）/562949953421312-（134********82487*（x-17）^3）/36028797018963968-（3896133174998587*（x-18）^3）/27021597764222976-5923044627661189/562949953421312
在
时的函数表达式为
S=（255325723948709*x）/281474976710656+（614517988558921*（x-18）^3）/135********111488+（134********82487*（x-19）^3）/36028797018963968-3648122929290727/562949953421312
在
时的函数表达式为
S=（157********7967*x）/562949953421312-（1412101843757803*（x-19）^3）/135********111488-（614517988558921*（x-20）^3）/135********111488+152********79421/281474976710656
所需光缆长度的近似值为:
sl=
25.4533
最后得到用三次样条插值拟合所测数据点的曲线如下图所示。
由此可以看出三次样条插值在节点处的光滑性得到了改善，但是要使用方程组计算增大了计算量，同时还要附加边界条件，分段三次样条插值对图形的控制能力还不够灵活。
图中的圆圈是取的基点值，从图中可以看到三次样条插值达到较好的效果。
复化Simpson公式近似计算光缆长度
铺设河底光缆的曲线图如上图所示。
本实验采用复化Simpson公式近似计算光缆长度。
2.2.1算法组织
1、算法依据及思想
由以上三次样条插值结果可以得到在20个区间上的是三次插值函数，结果如上述运行结果所示。
本题要求近似计算光缆长度，即在每个区间上采用弧长积分公式，便可得到所在区间的近似光缆长度。
各个区间的近似光缆长度之和即为所求近似计算光缆长度。
其中
为对应区间上的三次插值函数。
利用复化Simpson公式近似计算光缆长度的基本思想如下。
利用定积分对区间的区间可加性，将求积公式分别应用于每一个小区间上，设
是区间上[a,b]上满足
的n+1个点，由于
。
取
等距地分布在区间[a,b]上，此时
，
=
.。
得到复化梯形公式为：
误差为：
得到复化simpson公式为：
误差为：
本题采用在每个区间上20等分的方法，近似求解各个区间的定积分。
2、算法结构
首先本实验将所在的每个区间分别进行20等分，然后采用如下算法结构：
升级会员 