matlab微分方程数值求解（MATLAB推导钟摆的运动方程）

钟摆模型主要物理参量如下图所示：

钟摆模型

钟摆物理方程

syms m a g theta(t)

eqn = m*a == -m*g*sin(theta)

对于有长度的摆锤，摆锤的加速度：

syms r

eqn = subs(eqn,a,r*diff(theta,2))

使用隔离，隔离方程eqn中的角加速度。得到：

eqn = isolate(eqn,diff(theta,2))

将常数和频率收集到一个参数中，该参数也称为固有频率：

syms omega_0

eqn = subs(eqn,g/r,omega_0^2)

运动方程是非线性的，所以很难解析求解。假设角度很小，用泰勒展开式将方程线性化：

syms x

approx = taylor(sin(x),x,'Order',2);

approx = subs(approx,x,theta(t))

eqnLinear = subs(eqn,sin(theta(t)),approx)

用dsolve解方程组。将初始条件指定为第二个参数。通过使用假设来简化解决方案：

syms theta_0 theta_t0

theta_t = diff(theta);

cond = [theta(0) == theta_0, theta_t(0) == theta_t0];

assume(omega_0,'real')

thetaSol(t) = dsolve(eqnLinear,cond)

钟摆周期：

gValue = 9.81;

rValue = 1;

omega_0Value = sqrt(gValue/rValue);

T = 2*pi/omega_0Value;

theta_0Value = 0.1*pi; % Solution only valid for small angles.

theta_t0Value = 0; % Initially at rest.

vars = [omega_0 theta_0 theta_t0];

values = [omega_0Value theta_0Value theta_t0Value];

thetaSolPlot = subs(thetaSol,vars,values);

画出钟摆运动：

fplot(thetaSolPlot(t*T)/pi, [0 5]);

grid on;

title('Harmonic Pendulum Motion');

xlabel('t/T');

ylabel('\theta/\pi');

钟摆动画：

x_pos = sin(thetaSolPlot);

y_pos = -cos(thetaSolPlot);

fanimator(@fplot,x_pos,y_pos,'ko','MarkerFaceColor','k','AnimationRange',[0 5*T]);

hold on;

fanimator(@(t) plot([0 x_pos(t)],[0 y_pos(t)],'k-'),'AnimationRange',[0 5*T]);

fanimator(@(t) text(-0.3,0.3,"Timer: " num2str(t,2) " s"),'AnimationRange',[0 5*T]);

,