数学建模单摆问题

1数学建模实验第一次实验：单摆摆动问题分析结构9109175004陈建勇2单摆摆动问题分析摘要根据平常接触到的摆钟、秋千等实物中，我们可以抽象出单摆的模型。针对本题题目要求，我们抽象出最为简单的单摆模型，在此基础上加以分析。在理想条件下，单摆的摆动规律大致分为两种情况：小角度摆动和大角度摆动，分别针对这两种情况，从摆动微分方程出发，之后采取不同的方法分析。小角度摆动时，可做三角近似代替，将非线性微分方程转化为线性微分方程，进而求出其解析解，得到小摆角时单摆运动规律。通过matlab软件的验证，我们可以明显的看出结果与实际相符的很好。针对第二种情况，根据文献由相图的角度得到单摆的准确周期公式，但是实际计算比较困难，借助数学软件可以方便的求出实际准确周期。同样我们也可以用近似的方法计算单摆周期，如格林函数法、余弦函数法得到的近似计算公式，可根据误差要求选取不同的公式。本题的难点主要在于摆动微分方程的求解，在不同的条件下对近似解与精确解的选择。关键字：单摆周期小角度摆动大角度变动摆动微分方程3一：问题重述针对理想条件下的单摆，分析在小摆角和大摆角两种不同情况下的运动规律。二：模型假设1.悬挂小球的细线伸缩和质量均忽略不记，线长比小球的直径大得多；2.装置严格水平；3.不受空气阻力，且无驱动力。三：符号说明符号含义（rad）单摆偏离平衡位置的角位移0(rad)单摆的最大摆角g重力加速度，取9.8m/s2l细线长，取1mt(s)单摆摆动时间四：模型建立与求解1.最简单的单摆模型(如图1)图1简单单摆模型2.小角度时单摆运动规律（5o）单摆的运动微分方程1为：22ddt+sinlg=0(1)当摆角很小时,sin,故方程1可简化为:422ddt+lg=0（2）这是一个简单的谐振动方程，其解析解为：=Acos(00)(3)其固定角频率为：0=lg（4）得其周期为：T0=g220l（5）可以利用matlab软件在[0,5o]分别作出方程（1）和方程（2）的解得图像，如图图2小角度单摆摆动规律（—方程（1）的解，**方程（2）的解）由图像可以看出两方程的解的图像几乎吻合，可以说明当较小时（5o），两方程的解几乎相等，故周期公式此时较为准确。上述结论仅仅适用于摆角很小时（5o），当摆角很大时，方程sin不5再成立，方程（1）和方程（2）的解不再相近，故周期公式（5）不再成立。下面我们继续讨论摆角比较大时的单摆运动规律。3.大摆角时单摆运动规律（1）文献【2】从相图的角度得出单摆运动周期的精确解为：T=202020sin)2/(sin1dT2（6）为研究大摆角时单摆运动周期准确解，我们利用matlab软件在[0,2]区间上做出T0、0TT的图像，得到周期的准确解。图3单摆大摆角周期准确解如图可以看出，随着摆角的增大，单摆的运动周期逐渐增大，0TT也随之增大。（2）其他参考公式为求出单摆的周期，有时我们仅仅需要比较简单的近似公式还计算，T2为文献【3】利用格林函数法得到的近似公式：T2=T0（1+384-164020）（7）经计算，T2的相对误差随着0的增加而迅速增加，当摆角达到360（540）时6误差达到0.1%（0.5%），这对于精确计算已经不满足要求。3T为文献【4】给出另外一个简单近似计算公式：3T=T0)2/cos(10（8）该公式简单实用，由公式可以计算，当摆角为570是相对误差为0.1%；当摆角900时误差未达到0.75%。T4为文献【5】利用另外一种近似求解的方法——余弦函数法求的周期：T4=T0{1+2sin42312sin410402+……}（9）由公式显而易见可得，周期随着摆角的增大而增大。五．结论对于类似单摆的运动问题，我们很希望求得微分方程的严格的解析表达式进而得到完全准确的运动规律，但很多时候都无法直接求出这类非线性微分方程准确解，那么我们则可以考虑采取近似求解的方式，之后再采取一定的方式验证结果；或者可以直接借助数学软件接触需要的数值解。Matlab软件为解决这些数学问题提供了很大的方便，我们可以方便的利用它为我们排忧解惑。7参考文献：【1】万明理，何金娜，基于matlab下对单摆实验中大摆角问题的讨论[J]，大学物理实验，第23卷6期，2010年12月；【2】金亚平，单摆周期的相图求解[J]，大学物理，2000，19（10）；【3】孙春峰，非线性单摆的格林函数解法[J]，大学物理，2004,23（1）；【4】RRKiddandS.LFogg,”Asimpleformularforthelarge-anglepen-dulumperiod”[J],Phys,Teacher2002;【5】韦德泉，王秋芳，单摆角振幅对其周期的影响[J]，洛阳大学学报，2003,18（2）.附件：程序1：functionxp=pp1(t,x)xp=zeros(2,1);xp(1)=x(2);xp(2)=-9.8*sin(x(1));endfunctionwp=pp2(t,x)wp=zeros(2,1);wp(1)=x(2);wp(2)=-9.8*x(1);endt0=0;tf=10;[t,x]=ode45('pp1',[t0,tf],[0,pi/36]);[t1,x1]=ode45('pp2',[t0,tf],[0,pi/36]);plot(t,x(:,1),'k-');holdonplot(t1,x1(:,1),'k*')axis([010-0.040.04])title('С½Ç¶Èµ¥°Ú°Ú¶¯¹æÂÉ')xlabel('t/s')ylabel('¦È/rad')text(3,-0.032,'¡ª¡ª·½³Ì£¨1£©µÄ½â')text(3,-0.036,'*****·½³Ì£¨2£©µÄ½â')程序2：a=0;b=pi/2;n=1000;h=(b-a)/n;l=1;g=9.8;8h1=pi/(2*n);c=a:h1:b;m=2*pi*sqrt(l/g);x=a;s=0;fori1=1:(n+1)f0=4*sqrt(l/g)/sqrt(1-(sin(c(i1)/2))^2*(sin(x))^2);fori2=1:nx=x+h;f1=4*sqrt(l/g)/sqrt(1-(sin(c(i1)/2))^2*(sin(x))^2);s=s+(f0+f1)*h/2;f0=f1;enddisp(s);s1(i1)=s;s=0;endq=s1/m;plot(c,s1,'linewidth',3)gridholdonplot(c,m,'linewidth',3)holdonplot(c,q,'linewidth',3)title('µ¥°Ú´ó°Ú½ÇʱÖÜÆÚ')xlabel('¦È/rad')ylabel('´ó°Ú½Çʱµ¥°ÚÖÜÆÚµÄ׼ȷ½â')text(0.8,1.9,'¹ÌÓÐÖÜÆÚT0')text(1.2,2.15,'´ó½Ç¶ÈÖÜÆÚT')text(0.8,1.1,'ÖÜÆÚ±ÈT/T0')