计算方法上机报告姚威.docx
- 文档编号:3717222
- 上传时间:2022-11-24
- 格式:DOCX
- 页数:59
- 大小:474.90KB
计算方法上机报告姚威.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等分,然后采用如下算法结构:
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 计算方法 上机 报告