How to solve the varying length of strings of double pendulum in the simulation?

1 次查看（过去 30 天）
reema shrestha
reema shrestha 2017-8-2
评论： reema shrestha 2017-8-3
In my code,while the simulation is running,the length of the strings vary though the lengths are constant in the code. How can I solve it?
function double_pendulum(m1,m2,l1,l2,theta1,theta2,tmax)
%value of the constants
m1=0.1;
m2=0.1;
l1=0.5;
l2=0.5;
theta1=86;
theta2=86;
g=9.8;
tmax=20;
t(1)=0;
%RK matrices
c=[0 1/2 1/2 1];
a=[0 0 0 0;0 1/2 0 0;0 0 1/2 0;0 0 0 1];
w=[1/6 2/6 2/6 1/6];
theta1=theta1*pi/180;
theta2=theta2*pi/180;
s(:,1)=[theta1 theta2 0 0]';
F=@(t,s)[s(3) s(4) (-m2*l1*(s(3))^2*sin(s(1)-s(2))*cos(s(1)-s(2))+g*m2*sin(s(2))*cos(s(1)-s(2))-m2*l2*(s(4)^2)*sin(s(1)-s(2))-(m1+m2)*g*sin(s(1)))/l1*(m1+m2)-m2*l1*(cos(s(1)-s(2))^2) (m2*l2*(s(4)^2)*sin(s(1)-s(2))*cos(s(1)-s(2))+g*sin(s(1))*cos(s(1)-s(2))*(m1+m2)+l1*(s(3))^2*sin(s(1)-s(2))*(m1+m2)-g*sin(s(2))*(m1+m2))/l2*(m1+m2)-m2*l2*cos(s(1)-s(2))^2]';
h=0.1;
i=1;
while t(i)<tmax %==================================================
%Computation of position and angular Velocity
k=zeros (length(s(:,1)),length(c));
for j=1 : length(c)
k(:,j)=h*F(t(i)+h*c(j),s(:,i)+k*a(j,:)')
end
s(:,i+1)=s(:,i)+k*w';
t(i+1)=t(i)+h;
i=i+1;
end
x1=l1*sin(s(1,:));
x2=l1*sin(s(1,:))+l2*sin(s(2,:)); % x- axis position of pendulum...... s(1,:)= All theta values
y1=-l1*cos(s(1,:)); %y-axis position of pendulum
y2=-l1*cos(s(1,:))-l2*cos(s(2,:));
arraysize=size(t);
axissize = [min(min(x1),min(x2)) max(max(x1),max(x2)) min(min(y1),min(y2)) max(max(y1),max(y2))]; % Sets axis size for animation.
max(arraysize);
pendulumtopx = sum(x1)/max(arraysize);
pendulumtopy = sum(y1)/max(arraysize);
for i = 1 : max(arraysize)
plotarrayy = [pendulumtopy y1(i)]; % Plots solution at each time interval
plotarrayx = [pendulumtopx x1(i)];
plotarray2x = [x1(i) x2(i)];
plotarray2y = [y1(i) y2(i)];
plot(x1(i),y1(i),'o',x2(i),y2(i),'o','markersize',10,'markerfacecolor','g')
hold on
plot(x1(1:i),y1(1:i),'r');
plot(x2(1:i),y2(1:i),'b');
plot(plotarrayx,plotarrayy)
plot(plotarray2x,plotarray2y)
title(['Double pendulum simulation'])
hold off
xlabel('x','fontsize',12)
ylabel('y','fontsize',12)
axis(axissize);
pause(0.001); % Pause time between animations.
i=i+1;
end
figure
subplot(4,1,1)
hold on
plot(s(1,:),s(2,:),'b')
title('Double pendulum phase portrait','fontsize',12)
xlabel('theta1','fontsize',12)
ylabel('theta2','fontsize',12)
subplot(4,1,2)
hold on
plot(s(1,:),s(3,:),'r')
title('Double pendulum phase portrait','fontsize',12)
xlabel('theta1','fontsize',12)
ylabel('velocity1','fontsize',12)
% axis([min(s1) max(s1) min(s2) max(s2)])
subplot(4,1,3)% Plotting angle velocity vs time
plot(t,s(3,:),t,0)
title('Angle velocity vs Time Graph');
xlabel('time');
ylabel('\theta1(t)')
subplot(4,1,4)
plot(t,s(4,:),t,0)
set(gca,'XLim',[0 tmax])
title('Angle velocity vs Time Graph')
xlabel('time')
ylabel('\theta2(t)')

请先登录，再回答此问题。

回答（1 个）

James Tursa
James Tursa 2017-8-2
Do you have a profile of how the string lengths change as a function of time? E.g., can you do something like this:
Have functions defined
function L = l1(t)
L = % put code here for length as function of t
end
function L = l2(t)
L = % put code here for length as function of t
end
And then make your function handle call l1 and l2 as functions of t:
F=@(t,s)[s(3) s(4) (-m2*l1(t)*(s(3))^2*sin(s(1)-s(2))*cos(s(1)-s(2))+g*m2*sin(s(2))*cos(s(1)-s(2))-m2*l2*(s(4)^2)*sin(s(1)-s(2))-(m1+m2)*g*sin(s(1)))/l1(t)*(m1+m2)-m2*l1(t)*(cos(s(1)-s(2))^2) (m2*l2(t)*(s(4)^2)*sin(s(1)-s(2))*cos(s(1)-s(2))+g*sin(s(1))*cos(s(1)-s(2))*(m1+m2)+l1(t)*(s(3))^2*sin(s(1)-s(2))*(m1+m2)-g*sin(s(2))*(m1+m2))/l2(t)*(m1+m2)-m2*l2(t)*cos(s(1)-s(2))^2]';
Of if the lengths are a function of something else, put that into the l1 and l2 function code.

类别

在 Help Center 和 File Exchange 中查找有关 Programming 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by