数学实验教程-综合实验报告范例

单摆运动这是一个我们熟悉的物理模型，可看作工程技术中一些振动问题的简化。图1中一根长L的（无弹性的）细线，一端固定，另一端悬挂一质量为m的小球，在重力作用下小球处于竖直的平衡位置。使小球偏离平衡位置一个初始角度，然后让它自由落下，在不考虑空气阻力的情况下，小球就会沿圆弧摆动。问题1建立该物理系统的数学模型；问题2当初始角度较小时，求解数学模型；问题3当初始角度较大时，还可用问题2中的方法吗？图1中以θ=0位平衡位置，以右边为正方向建立摆角θ的坐标系。在小球摆动过程中的任一位置θ，小球所受重力沿运动轨迹方向的分力为-mgsinθ（负号表示力的方向与θ的正方向相反），利用牛顿第二定律即得微分方程设小球初始偏离角度为θ0，且无初速，则方程的初始条件为求解(1)、(2)时，在θ0不大的条件下，可将方程（1）中的sinθ近似为θ，于是得到线性常系数微分方程容易算出方程（3）在初始条件（2）下的解为由解（4）显然可知，简谐运动的周期为解：理论说明部分描述单摆运动规律的微分方程（1）是2阶微分方程，无解析解，但可用Matlab或其它软件编程求其数值解，但都需要先将它化成方程组的形式。令则微分方程(1)化为初始条件转化为在前面的两式中，g=9.8，l=25，x10为10o=0.1745(弧度)及30o=0.5236(弧度)两种情况.对于近似解(4)式，周期（2）方法描述令y=x1，z=x2，步长h=0.1，则原方程组可改写为①向前欧拉公式对n=0，1，2，…，向前欧拉计算公式如下：②改进的欧拉公式对n=0，1，2，…，改进的欧拉计算公式如下：③二阶龙格—库塔公式对n=0，1，2，…，二阶龙格—库塔公式计算公式如下：④四阶龙格—库塔公式对n=0，1，2，…，四阶龙格—库塔公式计算公式如下：（3）计算数据及数据图形①Matlab库函数ode45(…)计算结果及图形A.计算数据下面的数据表分别是对应于两种情况下的计算结果比较表。初始条件为θ=10o时的数值解与近似解计算表t数值解近似解t数值解近似解00.17450.17453.4606-0.0974-0.09790.00070.17450.17453.7106-0.1187-0.11930.00150.17450.17453.9606-0.1372-0.13770.00220.17450.17454.2106-0.1523-0.15270.00300.17450.17454.4606-0.1637-0.16400.00660.17450.17454.7106-0.1711-0.17130.01030.17450.17454.9606-0.1743-0.17440.01400.17450.17455.2106-0.1734-0.17320.01770.17450.17455.4606-0.1681-0.16780.03620.17450.17455.7106-0.1588-0.15830.05460.17440.17445.9606-0.1457-0.14500.07310.17430.17436.2106-0.1289-0.12810.09150.17420.17426.4606-0.1091-0.10800.18380.17340.17336.7106-0.0866-0.08530.27610.17190.17196.9606-0.0619-0.06050.36830.16990.16997.2106-0.0358-0.03430.46060.16730.16737.4606-0.0087-0.00720.71060.15760.15757.71060.01850.02010.96060.14400.14397.96060.04530.04681.21060.12690.12678.21060.07100.07251.46060.10680.10658.46060.09490.09631.71060.08400.08368.71060.11650.11781.96060.05910.05878.96060.13530.13652.21060.03290.03249.21060.15080.15172.46060.00580.00539.46060.16260.16332.7106-0.0214-0.02209.59550.16740.16792.9606-0.0482-0.04879.73030.17090.17133.2106-0.0737-0.07429.86520.17330.173510.00000.17440.1745初始条件为θ=30o时的数值解与近似解计算表t数值解近似解t数值解近似解00.52360.52363.5507-0.3030-0.31780.00030.52360.52363.8007-0.3649-0.37880.00050.52360.52364.0507-0.4181-0.43050.00080.52360.52364.3007-0.4614-0.47170.00100.52360.52364.5507-0.4937-0.50140.00230.52360.52364.8007-0.5145-0.51880.00360.52360.52365.0507-0.5233-0.52350.00490.52360.52365.3007-0.5198-0.51540.00620.52360.52365.5507-0.5042-0.49470.01260.52360.52365.8007-0.4768-0.46190.01900.52360.52366.0507-0.4381-0.41790.02540.52350.52356.3007-0.3891-0.36360.03180.52350.52356.5507-0.3308-0.30040.06380.52320.52326.8007-0.2646-0.22980.09590.52270.52277.0507-0.1919-0.15370.12790.52200.52197.3007-0.1146-0.07380.15990.52110.52107.5507-0.03450.00790.32010.51360.51317.80070.04640.08940.48030.50110.50018.05070.12620.16880.64050.48390.48218.30070.20290.24400.80070.46190.45928.55070.27470.31321.05070.41880.41438.80070.33990.37481.30070.36570.35949.05070.39690.42721.55070.30390.29569.30070.44450.46921.80070.23480.22469.55070.48150.49972.05070.16000.14819.80070.50720.51802.30070.08130.06809.85050.51090.52012.55070.0006-0.01379.90040.51420.52172.8007-0.0801-0.09519.95020.51690.52293.0507-0.1589-0.174210.00000.51920.52353.3007-0.2338-0.2491B．数值解图形0012345678910-0.2-0.15-0.1-0.0500.050.10.150.2tthetaTheta=10o的数值解0012345678910-0.8-0.6-0.4-0.200.20.40.6tthetatheta=30的数值解C．近似解图形0012345678910-0.2-0.15-0.1-0.0500.050.10.150.2tthetaTheta=10o的近似解0012345678910-0.8-0.6-0.4-0.200.20.40.6tthetatheta=30o的近似解（3）简要结论从上面的数据可以看出，处事角度为10o时精确（数值）解与近似解相差不大，而初始角度为30o时，随着时间的增加差别就很大了。（4）计算程序①Matlab5.0程序如下:%定义微分方程组functionxdot=danbai(t,x)g=9.8;l=0.25;xdot=zeros(2,1);xdot(1)=x(2);xdot(2)=-g/l*sin(x(1));%调用ode45(…)解此微分方程组%周期近似为10s,因此，时间t的取值范围确定在[010]%微分方程组的初始解为：x0=[0.1745,0](即10o)，x0=[0.5236,0](即30o)，程序中只有第一种初始解；要得到第二种初始解，只需将a=0.1745改成0.5236重新运行即可。t0=0;tf=1.0;a=0.1745;x0=[a,0];[t,x]=ode45('danbai',[t0tf],x0);g=9.8;l=0.25;w=sqrt(g/l)