数学实验报告
日期：20 年 月 日
第二篇：数学实验教程_综合实验报告范例
单摆运动
这是一个我们熟悉的物理模型，可看作工程技术中一些振动问题的简化。
图1中一根长L的（无弹性的）细线，一端固定，另一端悬挂一质量为m的小球，在重力作用下小球处于竖直的平衡位置。使小球偏离平衡位置一个初始角度，然后让它自由落下，在不考虑空气阻力的情况下，小球就会沿圆弧摆动。
问题1 建立该物理系统的数学模型；
问题2 当初始角度较小时，求解数学模型；
问题3 当初始角度较大时，还可用问题2中的方法吗？
图1中以θ=0位平衡位置，以右边为正方向建立摆角θ的坐标系。在小球摆动过程中的任一位置θ，小球所受重力沿运动轨迹方向的分力为-mgsinθ（负号表示力的方向与θ的正方向相反），利用牛顿第二定律即得微分方程
设小球初始偏离角度为θ0，且无初速，则方程的初始条件为
求解(1)、(2)时，在θ0不大的条件下，可将方程（1）中的sinθ近似为θ，于是得到线性常系数微分方程
容易算出方程（3）在初始条件（2）下的解为
由解（4）显然可知，简谐运动的周期为
解：
（1） 理论说明部分
描述单摆运动规律的微分方程（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时的数值解与近似解计算表
初始条件为 θ=30o时的数值解与近似解计算表
B．数值解图形
C．近似解图形
（3）简要结论
从上面的数据可以看出，处事角度为10o时精确（数值）解与近似解相差不大，而初始角度为30o时，随着时间的增加差别就很大了。
（4）计算程序
① Matlab 5.0 程序如下:
% 定义微分方程组
function xdot=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的取值范围确定在[0 10]
% 微分方程组的初始解为：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',[t0 tf],x0);
g=9.8;l=0.25;w=sqrt(g/l);
y=a*cos(w*t);
% 输出数值解 x(:,1)和近似解 y
data=[t,x(:,1),y];
% 描出数值解 x(:,1)和近似解 y 的图形
hold on
plot(t,x(:,1),'-k*')
plot(t,y,'-r*')
plot(t,y-x(:,1),'-.bs')
xlabel('t')
ylabel('theta')
title('theta=0.1745 的数值解')
hold off