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 joint coordinates for animation purposes

x = [L*sin(Y(:,1)),  L*sin(Y(:,1))+L*sin(Y(:,2))];

y = [-L*cos(Y(:,1)), -L*cos(Y(:,1))-L*cos(Y(:,2))];

% Convert radians to degrees

ang = Y(:,1:2)*180/pi;

% Set up first frame

figure('Color', 'white')

subplot(2,1,1)

plot(T, ang, 'LineWidth', 2)

hh1(1) = line(T(1), ang(1,1), 'Marker', '.', 'MarkerSize', 20, ...

'Color', 'b');

hh1(2) = line(T(1), ang(1,2), 'Marker', '.', 'MarkerSize', 20, ...

'Color', 'r');

xlabel('time (sec)')

ylabel('angle (deg)')

subplot(2,1,2)

hh2 = plot([0, x(1,1);x(1,1), x(1,2)], [0, y(1,1);y(1,1), y(1,2)], ...

'.-', 'MarkerSize', 20, 'LineWidth', 2);

axis equal

axis([-2*L 2*L -2*L 2*L])

ht = title(sprintf('Time: %0.2f sec', T(1)));

% Get figure size

pos = get(gcf, 'Position');

width = pos(3);

height = pos(4);

% Preallocate data (for storing frame data)

mov = zeros(height, width, 1, length(T));%, 'uint8');

% Loop through by changing XData and YData

for id = 1:length(T)

% Update graphics data. This is more efficient than recreating plots.

set(hh1(1), 'XData', T(id), 'YData', ang(id, 1))

set(hh1(2), 'XData', T(id), 'YData', ang(id, 2))

set(hh2(1), 'XData', [0, x(id, 1)], 'YData', [0, y(id, 1)])

set(hh2(2), 'XData', x(id, :), 'YData', y(id, :))

set(ht, 'String', sprintf('Time: %0.2f sec', T(id)))

%下面输出 gif 动画

frame=getframe(gcf);

im=frame2im(frame);%制作gif文件，图像必须是index索引图像

[I,map]=rgb2ind(im,20);

if id ==1

imwrite(I,map,'doublependulumanimation.gif','gif', 'Loopcount',inf,'DelayTime',0.5);%第一次必须创建！

else

imwrite(I,map,'doublependulumanimation.gif','gif','WriteMode','append','DelayTime',0.5);

end

end