数值大作业Word文档格式.docx

1、假设一个无阻尼的理想单摆，初始摆角为，摆长为l，小球的质量为m，重力加速度为g模型求解：根据公式：转动惯量角加速度 =力矩，得1.1小角度模型在小角度（一般认为）的情况下，作近似则运动方程为这是一个典型的常微分方程，可以解得周期为1.2 大角度模型当初始摆角较大（一般认为）时，上述的近似处理误差很大，不能使用该方法来求解周期。此时需要精确求解。根据微积分和微分方程的知识，可以解得运动周期的精确表达式为设k为精确周期与近似周期之比，即不妨设则当摆角为时，周期比为为了为之后的计算值提供对比，首先用MATLAB内置函数ellipke计算周期比，将其作为标准值，代码如下所示：theta=1:90;%b

2、y degreethm=theta*pi/180;%by radiansk=ellipke（sin（thm/2）.2）*2/pi;%resultxlswrite（elli.xls,theta;to） %write results to EXCEL document计算结果如下表所示作出摆角与周期比的函数图像，如下所示。可以看到，当角度很小时，周期比约等于1；当角度增大时，周期比偏离1。可见，当单摆作大角度运动时，近似周期公式是不可靠的。2. 函数编写为了更好地体现数值计算，在本文的函数计算中，均使用四则运算，所以需要编自行写正弦函数mysin和开方函数mysqrt。然后在此基础上，用多种数值积

3、分方法求函数的值。2.1 mysin函数编写（计算正弦值）主要思想是利用泰勒公式来求sin（x）的多项式逼近。根据泰勒公式，sin（x）的泰勒展开式如下所示：在本函数中，控制计算值达到小数点后10位，函数代码如下所示：%This function calculate the value of sin（x）%Using Taylors expansionfunction y = mysin（x）y = 0;tol = 0.5e-10;i = 0;er = 1;whileertol temp=y; y = y+（-1）i*x（2*i+1）/factorial（2*i+1）;er = abs（tem

4、p-y）;i = i+1;end2.2 mysqrt函数编写（计算开方值）%This function calculate the square root of inpout x%Argument x must between zero and onefunction y = mysqrt（x）low = 0;high = 1; %the initial lower and upper bounds %the principle of 10 significant digitsmid = （high+low）/2; %the estimate of zero point valuesymstf

5、（t） = t2-x;while（ertol）if f（mid）=0break;elseif f（mid）*f（low） high=mid;else low=mid; temp=mid; %store the last estimated value mid=（high+low）/2;er=abs（temp-mid）;y = mid;2.3 f（x）函数编写令其中x是初始摆角，即最大摆角。则周期比可表示为下面分别用复合梯形积分法、Romberg积分法、Gauss积分法，来构造函数。函数的输入变量是摆角（角度值），输出变量是周期比。2.3.1 复合梯形积分法%You should input t

6、he maximum swing angle（by degree） as the Argument%This function return the relative cycle of Simple Pendulumfunction y2 = tra（x）format longx=x/180*pi;g=（t）1/mysqrt（1-（mysin（x/2）2）*（mysin（t）2）;%composite trapzoidal formulaa=0;b=pi/2; %integrating boundsn=50;h=（b-a）/n; %step sizes=0;for t=a+h:h:b-h s=

7、s+g（t）;I=h/2*（g（a）+g（b）+h*s;y2=I*2/pi;2.3.2 Romberg积分法function y3 = rom（x）%Romberg method %a and b are upper and lower limits of integrationm=1;h=b-a;j=0;R（1,1）=h*（g（a）+g（b）/2;while（j4） j=j+1; h=h/2; s=0;for p=1:m t=a+h*（2*p-1）; R（j+1,1）=R（j,1）/2+h*s; m=m*2;for k=1:j R（j+1,k+1）=R（j+1,k）+（R（j+1,k）-R（

8、j,k）/（4k-1）;2/pi*Ry3=2/pi*R（j+1,j+1）;2.3.3 Gauss积分法（两点）function y1 = gau（x）%Gauss methodI=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.4 myspline函数编写（计算三次样条函数）为了能够通过少许插值点来知道任意摆角的周期比，可以通过对插值点进行三次样条插值，具体说明见后文。该函数的输入量是摆角，输出量是周期比，计算代码如下所示：%x is independent

9、 variable of this function%t and v is the vector of sampling pointfunction y = myspline（ x,t,v ）n = length（t）; % Length of the input%The coefficients in A matrix.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=b for calculate coefficent matrix M .% S（i）

10、is the result of y in each spline section.A = eye（n）*2;S = zeros（n-1,1）;for i = 1:n-1 h（i） = t（i+1） - t（i）;for i = 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（

11、i+1） = g（i）;% Calculate the matrix M for later spline interpolationM = Ab; 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）;if x t（i+1）3. 运动周期求解在f（x）函数编写模块中，已经用复合梯形、Romberg和Gauss积分法三种方法写了对的求解算法。由于Gauss积分法的精度不高，故不采用该方

12、法。复合梯形积分法可以通过增加小区间的数量来提高精度，Romberg积分法亦可以通过增加外推次数来提高精度。又由于Romberg方法可以通过控制外推终止条件来控制计算结果的精度，所以下面主要使用Romberg方法来计算（精确到小数点后6位）：为了达到需要的精度，将代码做一些修改（修改部分加粗）tol=0.5e-6;err=1; %The initial value of error4&err err=2/pi*abs（R（j+1,j+1）-R（j,j）;通过以下代码计算各摆角对应的周期，并将数据写入Excel表中。（注：为了方便，以3为步长从0到90取值）step=3;i=0;for x=0:

13、step:90i=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插值点的选取选好插值点后，再调用之前写好的my

14、spline函数即可计算各个样条函数，输出摆角对应周期比的估计值，以计算对应的周期比为例：myspline（87,theta,to）;输出结果为1.1667044.3 结果分析首先画出插值点的散点图和样条函数的图像：%the scatter plotset（0,defaultfigurecolor,w）; %set the background color to white0.01: y1（i）=myspline（x,theta,to）;x=0:scatter（x,y1,.rhold onscatter（theta,to,*k图像如下所示，可见吻合地很好：为了更加精确地分析样条函数的误差情况，

15、下面作出各点的样条函数值和标准值的差的函数以为步长）%error plot1: 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数量级，精度已经很高了。从图中也可以看出，误差都是为负的，也就是说样条函数的估计值要小于标准值。所以，各个区间的样条函数均是

16、下凹程度要大于实际的函数。如果还需要提高精度，可以继续增加插值点。但是采用样条函数的意义本在于为了更加方便计算任意摆角的周期比；若是插值点过多，不同区间用不同的样条函数，导致样条函数也更多，反而繁琐。所以在实际使用中应当在两者中权衡。总结本文中提供了复合梯形积分法、Romberg积分法、Gauss积分法三种数值积分法来计算周期比，其中重点分析了Romberg积分法，因为其比较方便控制精度；梯形积分法也很适合，篇幅有限，就没有重复赘述。为了能够得到周期比与摆角的近似函数，以方便计算任意角度的周期比，本文中通过在标准值中取了若干插值点来求三次样条函数。若是三次样条函数在某些值处误差较大，可以通过在该处添加插值点，从而大大提高精确度。