cylin4.m

matlab中用于解决高等动力学hamilton系统的相关程序

function cylin3
clc
clear all
[t,Y] = ode23(@equ3,[0 100],[-10 0 0 0]);
global t
subplot(2,2,1);
plot(t,Y(:,1))
title('x')
grid;
subplot(2,2,2);
plot(t,Y(:,2))
title('x dot')
grid;

subplot(2,2,3);
plot(t,Y(:,3)) 
title('teta')
grid;
subplot(2,2,4);
plot(t,Y(:,4)) 
title('teta dot')
grid;

function dy = equ3(t,y)
dy = zeros(4,1);

g=9.81;
T=100;
I=50;
muk=0;

r=0.5;
R=1;
m=100;
mic=0;

h=2;



dy(1) = y(2);
dy(2) = (((abs(R*y(4)-y(2))<=mic)*(m*R*T*(R*cos(atan2((r-1/2*h./(y(1).^2+h^2).*(2*h*r - 2*sign(y(1)).*(y(1).^4-y(1).^2*r^2+h^2*y(1).^2).^(1/2)))./y(1), 1/2./(y(1).^2+h^2).*(2*h*r - 2*sign(y(1)).*(y(1).^4-y(1).^2*r^2+h^2*y(1).^2).^(1/2)))) + r)/(m*R^2 + I) - T*cos(atan2((r-1/2*h./(y(1).^2+h^2).*(2*h*r - 2*sign(y(1)).*(y(1).^4-y(1).^2*r^2+h^2*y(1).^2).^(1/2)))./y(1), 1/2./(y(1).^2+h^2).*(2*h*r - 2*sign(y(1)).*(y(1).^4-y(1).^2*r^2+h^2*y(1).^2).^(1/2))))) + (abs(R*y(4)-y(2))>mic)*(-muk*(m*g - T*sin(atan2((r-1/2*h./(y(1).^2+h^2).*(2*h*r - 2*sign(y(1)).*(y(1).^4-y(1).^2*r^2+h^2*y(1).^2).^(1/2)))./y(1), 1/2./(y(1).^2+h^2).*(2*h*r - 2*sign(y(1)).*(y(1).^4-y(1).^2*r^2+h^2*y(1).^2).^(1/2)))))*sign(y(2)-R*y(4)))) + T*cos(atan2((r-1/2*h./(y(1).^2+h^2).*(2*h*r - 2*sign(y(1)).*(y(1).^4-y(1).^2*r^2+h^2*y(1).^2).^(1/2)))./y(1), 1/2./(y(1).^2+h^2).*(2*h*r - 2*sign(y(1)).*(y(1).^4-y(1).^2*r^2+h^2*y(1).^2).^(1/2)))))/m;
dy(3) = y(4);
dy(4) = (T*r - ((abs(R*y(4)-y(2))<=mic)*(m*R*T*(R*cos(atan2((r-1/2*h./(y(1).^2+h^2).*(2*h*r - 2*sign(y(1)).*(y(1).^4-y(1).^2*r^2+h^2*y(1).^2).^(1/2)))./y(1), 1/2./(y(1).^2+h^2).*(2*h*r - 2*sign(y(1)).*(y(1).^4-y(1).^2*r^2+h^2*y(1).^2).^(1/2)))) + r)/(m*R^2 + I) - T*cos(atan2((r-1/2*h./(y(1).^2+h^2).*(2*h*r - 2*sign(y(1)).*(y(1).^4-y(1).^2*r^2+h^2*y(1).^2).^(1/2)))./y(1), 1/2./(y(1).^2+h^2).*(2*h*r - 2*sign(y(1)).*(y(1).^4-y(1).^2*r^2+h^2*y(1).^2).^(1/2))))) + (abs(R*y(4)-y(2))>mic)*(-muk*(m*g - T*sin(atan2((r-1/2*h./(y(1).^2+h^2).*(2*h*r - 2*sign(y(1)).*(y(1).^4-y(1).^2*r^2+h^2*y(1).^2).^(1/2)))./y(1), 1/2./(y(1).^2+h^2).*(2*h*r - 2*sign(y(1)).*(y(1).^4-y(1).^2*r^2+h^2*y(1).^2).^(1/2)))))*sign(y(2)-R*y(4))))*R)/I;

%f = ((abs(R*y(4)-y(2))<=mic)*(m*R*T*(R*cos(atan2((r-1/2*h./(y(1).^2+h^2).*(2*h*r - 2*sign(y(1)).*(y(1).^4-y(1).^2*r^2+h^2*y(1).^2).^(1/2)))./y(1), 1/2./(y(1).^2+h^2).*(2*h*r - 2*sign(y(1)).*(y(1).^4-y(1).^2*r^2+h^2*y(1).^2).^(1/2)))*180/pi) + r)/(m*R^2 + I) - T*cos(atan2((r-1/2*h./(y(1).^2+h^2).*(2*h*r - 2*sign(y(1)).*(y(1).^4-y(1).^2*r^2+h^2*y(1).^2).^(1/2)))./y(1), 1/2./(y(1).^2+h^2).*(2*h*r - 2*sign(y(1)).*(y(1).^4-y(1).^2*r^2+h^2*y(1).^2).^(1/2)))*180/pi)) + (abs(R*y(4)-y(2))>mic)*(-muk*(m*g - T*sin(atan2((r-1/2*h./(y(1).^2+h^2).*(2*h*r - 2*sign(y(1)).*(y(1).^4-y(1).^2*r^2+h^2*y(1).^2).^(1/2)))./y(1), 1/2./(y(1).^2+h^2).*(2*h*r - 2*sign(y(1)).*(y(1).^4-y(1).^2*r^2+h^2*y(1).^2).^(1/2)))*180/pi))*sign(y(2)-R*y(4))))
%f1 = (-muk*(m*g - T*sin(atan2((r-1/2*h./(y(1).^2+h^2).*(2*h*r - 2*sign(y(1)).*(y(1).^4-y(1).^2*r^2+h^2*y(1).^2).^(1/2)))./y(1), 1/2./(y(1).^2+h^2).*(2*h*r - 2*sign(y(1)).*(y(1).^4-y(1).^2*r^2+h^2*y(1).^2).^(1/2)))*180/pi))*sign(y(2)-R*y(4)))
%f2 = (m*R*T*(R*cos(atan2((r-1/2*h./(y(1).^2+h^2).*(2*h*r - 2*sign(y(1)).*(y(1).^4-y(1).^2*r^2+h^2*y(1).^2).^(1/2)))./y(1), 1/2./(y(1).^2+h^2).*(2*h*r - 2*sign(y(1)).*(y(1).^4-y(1).^2*r^2+h^2*y(1).^2).^(1/2)))*180/pi) + r)/(m*R^2 + I) - T*cos(atan2((r-1/2*h./(y(1).^2+h^2).*(2*h*r - 2*sign(y(1)).*(y(1).^4-y(1).^2*r^2+h^2*y(1).^2).^(1/2)))./y(1), 1/2./(y(1).^2+h^2).*(2*h*r - 2*sign(y(1)).*(y(1).^4-y(1).^2*r^2+h^2*y(1).^2).^(1/2)))*180/pi))