双摆模型matlab,科学网—双摆的MATLAB程序模拟 - 王福昌的博文

g = 9.81;

theta1 = y(1);       % angle 1

theta2 = y(2);       % angle 2

p1 = y(3);           % momentum

p2 = y(4);           % momentum

% The derivatives

dy(1) = 6/(m*L^2) * (2*p1-3*cos(theta1-theta2)*p2) / ...

(16-9*cos(theta1-theta2)^2);

dy(2) = 6/(m*L^2)*(8*p2-3*cos(theta1-theta2)*p1) / ...

(16-9*cos(theta1-theta2)^2);

dy(3) = -1/2*m*L^2*(dy(1)*dy(2)*sin(theta1-theta2)+3*g/L*sin(theta1));

dy(4) = -1/2*m*L^2*(-dy(1)*dy(2)*sin(theta1-theta2)+g/L*sin(theta2));

dy = dy(:);

主程序

m = 1;

L = 1;

theta1 = 3*pi/4;

theta2 = 3*pi/8;

t = linspace(0, 10, 200);

% Solving ODE of a double pendulum

[T,Y] = ode45(@(t, x) double_pendulum(t, x, m, L), ...

t, [theta1, theta2, 0, 0]);

% Calculating join