syms theta1(t) theta2(t) theta3(t)
m1=1; m2=1; m3=1; l1=1; l2=1;
l3=1; g=9.81; tau1=1; tau2=1; tau3=1;
e1= diff(theta1,2)*(l1^2*(m1+m2+m3)+l1)+...
diff(theta1,2)*((m2+m3)*(l1*l2*cos(theta1-theta2)))+...
diff(theta3,2)*(m3*l1*l3*cos(theta1-theta3))==...
-(m2+m3)*(l1*l2*diff(theta2)^2*sin(theta1-theta2))-...
m3*l1*l3*diff(theta3)^2*sin(theta1-theta3)-...
(m1+m2+m3)*g*l1*cos(theta1)+tau1-tau2;
e2= diff(theta1,2)*((m2+m3)*l1*l2*cos(theta1-theta2))+...
diff(theta2,2)*(l2^2*(m2+m3)+l2)+...
diff(theta3,2)*(m3*l2*l3*cos(theta2-theta3))==...
(m2+m3)*l1*l2*diff(theta1)^2*sin(theta1-theta2)-...
m3*l2*l3*diff(theta3)^2*sin(theta2-theta3)-...
(m2+m3)*g*l2*cos(theta2)+tau2-tau3;
e3= diff(theta1,2)*(m3*l1*l3*cos(theta1-theta3))+...
diff(theta2,2)*(m3*l2*l3*cos(theta2-theta3))+...
diff(theta3,2)*(l3^2*m3+l3)==...
m3*l1*l3*diff(theta1)^2*sin(theta1-theta3)+...
m3*l2*l3*diff(theta2)^2*sin(theta2-theta3)-...
m3*g*l3*cos(theta3)+tau3;
vars = [theta1(t); theta2(t); theta3(t)]
V = odeToVectorField([e1,e2,e3])
M = matlabFunction(V,'vars', {'t','Y'})
interval = [0 20];
y0 = [0; pi/2; 0; pi/2; 0; pi/2];
ySol = ode45(M,interval,y0);
tValues = linspace(interval(1),interval(2),1000);
for i = 1:6
yValues = deval(ySol,tValues,i);
plot(tValues,yValues)
hold on
end
Best Answer