数值大作业Word文档格式.docx
《数值大作业Word文档格式.docx》由会员分享，可在线阅读，更多相关《数值大作业Word文档格式.docx（16页珍藏版）》请在冰豆网上搜索。
假设一个无阻尼的理想单摆，初始摆角为
，摆长为l，小球的质量为m，重力加速度为g
模型求解：
根据公式：
转动惯量·
角加速度=力矩，得
1.1小角度模型
在小角度（一般认为
）的情况下，作近似
则运动方程为
这是一个典型的常微分方程，可以解得周期为
1.2大角度模型
当初始摆角
较大（一般认为
）时，上述的近似处理
误差很大，不能使用该方法来求解周期。
此时需要精确求解。
根据微积分和微分方程的知识，可以解得运动周期的精确表达式为
设k为精确周期
与近似周期
之比，即
不妨设
则当摆角为
时，周期比为
为了为之后的计算值提供对比，首先用MATLAB内置函数ellipke计算周期比，将其作为标准值，代码如下所示：
theta=1:
90;
%bydegree
thm=theta*pi/180;
%byradians
k=ellipke（sin（thm/2）.^2）*2/pi;
%result
xlswrite（'
elli.xls'
[theta;
to]）%writeresultstoEXCELdocument
计算结果如下表所示
作出摆角与周期比的函数图像，如下所示。
可以看到，当角度很小时，周期比约等于1；
当角度增大时，周期比偏离1。
可见，当单摆作大角度运动时，近似周期公式是不可靠的。
2.函数编写
为了更好地体现数值计算，在本文的函数计算中，均使用四则运算，所以需要编自行写正弦函数mysin和开方函数mysqrt。
然后在此基础上，用多种数值积分方法求函数
的值。
2.1mysin函数编写（计算正弦值）
主要思想是利用泰勒公式来求sin（x）的多项式逼近。
根据泰勒公式，sin（x）的泰勒展开式如下所示：
在本函数中，控制计算值达到小数点后10位，函数代码如下所示：
%Thisfunctioncalculatethevalueofsin（x）
%UsingTaylor'
sexpansion
function[y]=mysin（x）
y=0;
tol=0.5e-10;
i=0;
er=1;
whileer>
tol
temp=y;
y=y+（-1）^i*x^（2*i+1）/factorial（2*i+1）;
er=abs（temp-y）;
i=i+1;
end
2.2mysqrt函数编写（计算开方值）
%Thisfunctioncalculatethesquarerootofinpoutx
%Argumentxmustbetweenzeroandone
function[y]=mysqrt（x）
low=0;
high=1;
%theinitiallowerandupperbounds
%theprincipleof10significantdigits
mid=（high+low）/2;
%theestimateofzeropointvalue
symst
f（t）=t^2-x;
while（er>
tol）
iff（mid）==0
break;
elseiff（mid）*f（low）<
high=mid;
elselow=mid;
temp=mid;
%storethelastestimatedvalue
mid=（high+low）/2;
er=abs（temp-mid）;
y=mid;
2.3f（x）函数编写
令
其中x是初始摆角，即最大摆角。
则周期比可表示为
下面分别用复合梯形积分法、Romberg积分法、Gauss积分法，来构造
函数。
函数的输入变量是摆角（角度值），输出变量是周期比。
2.3.1复合梯形积分法
%Youshouldinputthemaximumswingangle（bydegree）astheArgument
%ThisfunctionreturntherelativecycleofSimplePendulum
function[y2]=tra（x）
formatlong
x=x/180*pi;
g=@（t）1/mysqrt（1-（mysin（x/2）^2）*（mysin（t））^2）;
%compositetrapzoidalformula
a=0;
b=pi/2;
%integratingbounds
n=50;
h=（b-a）/n;
%stepsize
s=0;
fort=a+h:
h:
b-h
s=s+g（t）;
I=h/2*（g（a）+g（b））+h*s;
y2=I*2/pi;
2.3.2Romberg积分法
function[y3]=rom（x）
%Rombergmethod
%aandbareupperandlowerlimitsofintegration
m=1;
h=b-a;
j=0;
R（1,1）=h*（g（a）+g（b））/2;
while（j<
4）
j=j+1;
h=h/2;
s=0;
forp=1:
m
t=a+h*（2*p-1）;
R（j+1,1）=R（j,1）/2+h*s;
m=m*2;
fork=1:
j
R（j+1,k+1）=R（j+1,k）+（R（j+1,k）-R（j,k））/（4^k-1）;
2/pi*R
y3=2/pi*R（j+1,j+1）;
2.3.3Gauss积分法（两点）
function[y1]=gau（x）
%Gaussmethod
I=pi/4*（5/9*g（pi/4-pi/4*sqrt（15）/5）+8/9*g（pi/4+pi/4*0）+5/9*g（pi/4+pi/4*sqrt（15）/5））;
y1=I*2/pi;
2.4myspline函数编写（计算三次样条函数）
为了能够通过少许插值点来知道任意摆角的周期比，可以通过对插值点进行三次样条插值，具体说明见后文。
该函数的输入量是摆角，输出量是周期比，计算代码如下所示：
%xisindependentvariableofthisfunction
%tandvisthevectorofsamplingpoint
functiony=myspline（x,t,v）
n=length（t）;
%Lengthoftheinput
%ThecoefficientsinAmatrix.
lambda=zeros（n-1,1）;
u=zeros（n-1,1）;
h=zeros（n-1,1）;
g=zeros（n-1,1）;
b=zeros（n,1）;
%A*M=bforcalculatecoefficentmatrixM.
%S（i）istheresultofyineachsplinesection.
A=eye（n）*2;
S=zeros（n-1,1）;
fori=1:
n-1
h（i）=t（i+1）-t（i）;
fori=2:
lambda（i）=h（i）/（h（i-1）+h（i））;
A（i,i+1）=lambda（i）;
u（i-1）=1-lambda（i）;
A（i,i-1）=u（i-1）;
n-2
g（i）=6/（h（i）+h（i+1））*...
（（v（i+1）-v（i+2））/（t（i+1）-t（i+2））-...
（v（i）-v（i+1））/（t（i）-t（i+1）））;
b（i+1）=g（i）;
%CalculatethematrixMforlatersplineinterpolation
M=A\b;
S（i）=1/（6*h（i））*（M（i）*（t（i+1）-x）^3+M（i+1）*（x-t（i））^3）+...
（v（i）-（M（i）*h（i）^2）/6）*（t（i+1）-x）/h（i）+...
（v（i+1）-（M（i+1）*h（i）^2）/6）*（x-t（i））/h（i）;
ifx<
=t（i+1）
y=S（i）;
break
ifx>
t（i+1）
3.运动周期求解
在f（x）函数编写模块中，已经用复合梯形、Romberg和Gauss积分法三种方法写了对
的求解算法。
由于Gauss积分法的精度不高，故不采用该方法。
复合梯形积分法可以通过增加小区间的数量来提高精度，Romberg积分法亦可以通过增加外推次数来提高精度。
又由于Romberg方法可以通过控制外推终止条件来控制计算结果的精度，所以下面主要使用Romberg方法来计算（精确到小数点后6位）：
为了达到需要的精度，将代码做一些修改（修改部分加粗）
tol=0.5e-6;
err=1;
%Theinitialvalueoferror
4&
&
err>
err=2/pi*abs（（R（j+1,j+1）-R（j,j）））;
通过以下代码计算各摆角对应的周期，并将数据写入Excel表中。
（注：
为了方便，以3°
为步长从0°
到90°
取值）
step=3;
i=0;
forx=0:
step:
90
i=i+1;
y=rom（x）;
result（i,1）=x;
result（i,2）=y;
rom.xls'
result）;
计算结果与标准值的对比结果如下：
可以看出，在精确到小数点后6位的情况下，计算值与标准值之间几乎没有误差。
4.函数拟合
为了能够方便计算任意一摆角的周期比，可以在标准值中取一些插值点，然后利用这些插值点计算三次样条函数。
有了三次样条函数后，就可以很方便地计算任意摆角的周期比。
4.1插值点的选取
每隔
取一个插值点，代码如下：
theta=0:
5:
to=ellipke（sin（thm/2）.^2）*2/pi;
4.2插值点的选取
选好插值点后，再调用之前写好的myspline函数即可计算各个样条函数，输出摆角对应周期比的估计值，以计算
对应的周期比为例：
myspline（87,theta,to）;
输出结果为1.166704
4.3结果分析
首先画出插值点的散点图和样条函数的图像：
%thescatterplot
set（0,'
defaultfigurecolor'
'
w'
）;
%setthebackgroundcolortowhite
0.01:
y1（i）=myspline（x,theta,to）;
x=0:
scatter（x,y1,'
.r'
holdon
scatter（theta,to,'
*k'
图像如下所示，可见吻合地很好：
为了更加精确地分析样条函数的误差情况，下面作出各点的样条函数值和标准值的差的函数
以
为步长）
%errorplot
1:
y1=myspline（x,theta,to）;
t=x/180*pi;
y2=ellipke（sin（t/2）.^2）*2/pi;
line（[x,x],[0,y1-y2]）;
图像如下所示：
可以看出，误差已经很小了，达到了10-5数量级，但是在两端的误差相对来说比较大。
此时可以增加相应位置的插值点来改善：
theta=[0:
4,5:
85,86:
90];
再次作出误差图像，如下：
可以看出，误差已经达到了10-8数量级，精度已经很高了。
从图中也可以看出，误差都是为负的，也就是说样条函数的估计值要小于标准值。
所以，各个区间的样条函数均是下凹程度要大于实际的函数。
如果还需要提高精度，可以继续增加插值点。
但是采用样条函数的意义本在于为了更加方便计算任意摆角的周期比；
若是插值点过多，不同区间用不同的样条函数，导致样条函数也更多，反而繁琐。
所以在实际使用中应当在两者中权衡。
总结
本文中提供了复合梯形积分法、Romberg积分法、Gauss积分法三种数值积分法来计算周期比，其中重点分析了Romberg积分法，因为其比较方便控制精度；
梯形积分法也很适合，篇幅有限，就没有重复赘述。
为了能够得到周期比与摆角的近似函数，以方便计算任意角度的周期比，本文中通过在标准值中取了若干插值点来求三次样条函数。
若是三次样条函数在某些值处误差较大，可以通过在该处添加插值点，从而大大提高精确度。
升级会员 