数值大作业.docx

1、数值大作业目录大角度单摆的运动周期求解 31. 单摆运动模型建立 31.1 小角度模型 32.2 大角度模型 32. 函数编写 52.1 mysin函数编写（计算正弦值） 52.2 mysqrt函数编写（计算开方值） 62.3 f(x)函数编写 72.3.1 复合梯形积分法 72.3.2 Romberg积分法 82.3.3 Gauss积分法（两点） 92.4 myspline函数编写（计算三次样条函数） 93. 运动周期求解 104. 函数拟合 134.1 插值点的选取 134.2 插值点的选取 134.3 结果分析 13总结 15大角度单摆的运动周期求解摘要：无阻尼条件下单摆的周期是第一类椭

2、圆积分，在小角度下可以通过取近似来计算周期。但是在大角度下上式不成立，在本文中利用数值积分方法对各个摆角下的周期进行求解，同时亦通过求三次样条插值函数来近似原函数。1. 单摆运动模型建立模型假设：1. 单摆作无阻尼运动2. 摆绳的质量可忽略，且在运动中无形变3. 悬挂的小球可以看做质点参数设计：假设一个无阻尼的理想单摆，初始摆角为，摆长为l，小球的质量为m，重力加速度为g模型求解：根据公式：转动惯量角加速度 =力矩，得1.1小角度模型在小角度（一般认为）的情况下，作近似则运动方程为这是一个典型的常微分方程，可以解得周期为1.2 大角度模型当初始摆角较大（一般认为）时，上述的近似处理误差很大，不

3、能使用该方法来求解周期。此时需要精确求解。根据微积分和微分方程的知识，可以解得运动周期的精确表达式为设k为精确周期与近似周期之比，即不妨设则当摆角为时，周期比为为了为之后的计算值提供对比，首先用MATLAB内置函数ellipke计算周期比，将其作为标准值，代码如下所示：theta=1:90;%by degreethm=theta*pi/180;%by radiansk=ellipke(sin(thm/2).2)*2/pi;%resultxlswrite(elli.xls,theta;to) %write results to EXCEL document计算结果如下表所示作出摆角与周期比的函数

4、图像，如下所示。可以看到，当角度很小时，周期比约等于1；当角度增大时，周期比偏离1。可见，当单摆作大角度运动时，近似周期公式是不可靠的。2. 函数编写为了更好地体现数值计算，在本文的函数计算中，均使用四则运算，所以需要编自行写正弦函数mysin和开方函数mysqrt。然后在此基础上，用多种数值积分方法求函数的值。2.1 mysin函数编写（计算正弦值）主要思想是利用泰勒公式来求sin(x)的多项式逼近。根据泰勒公式，sin(x)的泰勒展开式如下所示：在本函数中，控制计算值达到小数点后10位，函数代码如下所示：%This function calculate the value of sin(x

5、)%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(temp-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

6、= 0;high = 1; %the initial lower and upper boundstol = 0.5e-10; %the principle of 10 significant digitsmid = (high+low)/2; %the estimate of zero point valueer = 1;symstf(t) = t2-x;while(ertol)if f(mid)=0break;elseif f(mid)*f(low)0 high=mid;else low=mid;end temp=mid; %store the last estimated value m

7、id=(high+low)/2;er=abs(temp-mid);endy = mid;2.3 f(x)函数编写令其中x是初始摆角，即最大摆角。则周期比可表示为下面分别用复合梯形积分法、Romberg积分法、Gauss积分法，来构造函数。函数的输入变量是摆角（角度值），输出变量是周期比。2.3.1 复合梯形积分法%You should input the maximum swing angle(by degree) as the Argument%This function return the relative cycle of Simple Pendulumfunction y2 = tr

8、a(x)format longx=x/180*pi;symstg=(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=s+g(t);endI=h/2*(g(a)+g(b)+h*s;y2=I*2/pi;2.3.2 Romberg积分法%You should input the maximum swing angle(by degree) as the Ar

9、gument%This function return the relative cycle of Simple Pendulumfunction y3 = rom(x)format longx=x/180*pi;symstg=(t)1/mysqrt(1-(mysin(x/2)2)*(mysin(t)2);%Romberg methoda=0;b=pi/2; %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

10、=1:m t=a+h*(2*p-1); s=s+g(t);end 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(j,k)/(4k-1);endend2/pi*Ry3=2/pi*R(j+1,j+1);2.3.3 Gauss积分法（两点）%You should input the maximum swing angle(by degree) as the Argument%This function return the relative cycle of Simple Pendulumfunction

11、 y1 = gau(x)format longx=x/180*pi;symstg=(t)1/mysqrt(1-(mysin(x/2)2)*(mysin(t)2);%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 indep

12、endent 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 .%

13、 S(i) 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);endfor i = 2:n-1 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);endfor 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

14、)-v(i+1)/(t(i)-t(i+1); b(i+1) = g(i);end% Calculate the matrix M for later spline interpolationM = Ab;for i = 1:n-1 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) y = S(i);endend3. 运动周期求解在f(x)函数编写模块中，已经

15、用复合梯形、Romberg和Gauss积分法三种方法写了对的求解算法。由于Gauss积分法的精度不高，故不采用该方法。复合梯形积分法可以通过增加小区间的数量来提高精度，Romberg积分法亦可以通过增加外推次数来提高精度。又由于Romberg方法可以通过控制外推终止条件来控制计算结果的精度，所以下面主要使用Romberg方法来计算（精确到小数点后6位）：为了达到需要的精度，将代码做一些修改（修改部分加粗）%You should input the maximum swing angle(by degree) as the Argument%This function return the re

16、lative cycle of Simple Pendulumfunction y3 = rom(x)format longx=x/180*pi;symstg=(t)1/mysqrt(1-(mysin(x/2)2)*(mysin(t)2);%Romberg methoda=0;b=pi/2; %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;tol=0.5e-6;err=1; %The initial value of errorwhile(jtol) j=j+1; h=h

17、/2; s=0;for p=1:m t=a+h*(2*p-1); s=s+g(t);end 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(j,k)/(4k-1);end err=2/pi*abs(R(j+1,j+1)-R(j,j);end2/pi*Ry3=2/pi*R(j+1,j+1);通过以下代码计算各摆角对应的周期，并将数据写入Excel表中。（注：为了方便，以3为步长从0到90取值）step=3;i=0;for x=0:step:90i=i+1; y=rom(x); result(i,1) =

18、 x; result(i,2) = y;endxlswrite(rom.xls,result);计算结果与标准值的对比结果如下：可以看出，在精确到小数点后6位的情况下，计算值与标准值之间几乎没有误差。4. 函数拟合为了能够方便计算任意一摆角的周期比，可以在标准值中取一些插值点，然后利用这些插值点计算三次样条函数。有了三次样条函数后，就可以很方便地计算任意摆角的周期比。4.1 插值点的选取每隔取一个插值点，代码如下：theta=0:5:90;thm=theta*pi/180;to=ellipke(sin(thm/2).2)*2/pi;4.2插值点的选取选好插值点后，再调用之前写好的mysplin

19、e函数即可计算各个样条函数，输出摆角对应周期比的估计值，以计算对应的周期比为例：myspline(87,theta,to);输出结果为1.1667044.3 结果分析首先画出插值点的散点图和样条函数的图像：%the scatter plotset(0,defaultfigurecolor,w); %set the background color to whitei=0;for x=0:0.01:90i=i+1; y1(i)=myspline(x,theta,to);endx=0:0.01:90;scatter(x,y1,.r);hold onscatter(theta,to,*k);图像如下

20、所示，可见吻合地很好：为了更加精确地分析样条函数的误差情况，下面作出各点的样条函数值和标准值的差的函数（注：以为步长）%error plotset(0,defaultfigurecolor,w); %set the background color to whitefor x=0:1:90 y1=myspline(x,theta,to); t=x/180*pi; y2=ellipke(sin(t/2).2)*2/pi; line(x,x,0,y1-y2);end图像如下所示：可以看出，误差已经很小了，达到了10-5数量级，但是在两端的误差相对来说比较大。此时可以增加相应位置的插值点来改善：th

21、eta=0:1:4,5:5:85,86:1:90;thm=theta*pi/180;to=ellipke(sin(thm/2).2)*2/pi;再次作出误差图像，如下：可以看出，误差已经达到了10-8数量级，精度已经很高了。从图中也可以看出，误差都是为负的，也就是说样条函数的估计值要小于标准值。所以，各个区间的样条函数均是下凹程度要大于实际的函数。如果还需要提高精度，可以继续增加插值点。但是采用样条函数的意义本在于为了更加方便计算任意摆角的周期比；若是插值点过多，不同区间用不同的样条函数，导致样条函数也更多，反而繁琐。所以在实际使用中应当在两者中权衡。总结本文中提供了复合梯形积分法、Romberg积分法、Gauss积分法三种数值积分法来计算周期比，其中重点分析了Romberg积分法，因为其比较方便控制精度；梯形积分法也很适合，篇幅有限，就没有重复赘述。为了能够得到周期比与摆角的近似函数，以方便计算任意角度的周期比，本文中通过在标准值中取了若干插值点来求三次样条函数。若是三次样条函数在某些值处误差较大，可以通过在该处添加插值点，从而大大提高精确度。