cylin5.m

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

function cylin5
clc
clear all
[t,Y] = ode45(@equ5,[0 10],[-10 0 0 0]);

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 dx = equ5(t,x)
dx = zeros(4,1);

u(1)=9.81;
u(2)=100;
u(3)=50;
u(4)=0;

u(5)=0.5;
u(6)=1;
u(7)=100;
u(8)=0;
u(9)=2;

dx(1) = x(2);
dx(2) = (((abs(u(6)*x(4)-x(2))<=u(8))*(u(7)*u(6)*u(2)*(u(6)*cos(atan2((u(5)-1/2*u(9)./(x(1).^2+u(9)^2).*(2*u(9)*u(5) - 2*sign(x(1)).*(x(1).^4-x(1).^2*u(5)^2+u(9)^2*x(1).^2).^(1/2)))./x(1), 1/2./(x(1).^2+u(9)^2).*(2*u(9)*u(5) - 2*sign(x(1)).*(x(1).^4-x(1).^2*u(5)^2+u(9)^2*x(1).^2).^(1/2)))) + u(5))/(u(7)*u(6)^2 + u(3)) - u(2)*cos(atan2((u(5)-1/2*u(9)./(x(1).^2+u(9)^2).*(2*u(9)*u(5) - 2*sign(x(1)).*(x(1).^4-x(1).^2*u(5)^2+u(9)^2*x(1).^2).^(1/2)))./x(1), 1/2./(x(1).^2+u(9)^2).*(2*u(9)*u(5) - 2*sign(x(1)).*(x(1).^4-x(1).^2*u(5)^2+u(9)^2*x(1).^2).^(1/2))))) + (abs(u(6)*x(4)-x(2))>u(8))*(-u(4)*(u(7)*u(1) - u(2)*sin(atan2((u(5)-1/2*u(9)./(x(1).^2+u(9)^2).*(2*u(9)*u(5) - 2*sign(x(1)).*(x(1).^4-x(1).^2*u(5)^2+u(9)^2*x(1).^2).^(1/2)))./x(1), 1/2./(x(1).^2+u(9)^2).*(2*u(9)*u(5) - 2*sign(x(1)).*(x(1).^4-x(1).^2*u(5)^2+u(9)^2*x(1).^2).^(1/2)))))*sign(x(2)-u(6)*x(4)))) + u(2)*cos(atan2((u(5)-1/2*u(9)./(x(1).^2+u(9)^2).*(2*u(9)*u(5) - 2*sign(x(1)).*(x(1).^4-x(1).^2*u(5)^2+u(9)^2*x(1).^2).^(1/2)))./x(1), 1/2./(x(1).^2+u(9)^2).*(2*u(9)*u(5) - 2*sign(x(1)).*(x(1).^4-x(1).^2*u(5)^2+u(9)^2*x(1).^2).^(1/2)))))/u(7);
dx(3) = x(4);
dx(4) = (u(2)*u(5) - ((abs(u(6)*x(4)-x(2))<=u(8))*(u(7)*u(6)*u(2)*(u(6)*cos(atan2((u(5)-1/2*u(9)./(x(1).^2+u(9)^2).*(2*u(9)*u(5) - 2*sign(x(1)).*(x(1).^4-x(1).^2*u(5)^2+u(9)^2*x(1).^2).^(1/2)))./x(1), 1/2./(x(1).^2+u(9)^2).*(2*u(9)*u(5) - 2*sign(x(1)).*(x(1).^4-x(1).^2*u(5)^2+u(9)^2*x(1).^2).^(1/2)))) + u(5))/(u(7)*u(6)^2 + u(3)) - u(2)*cos(atan2((u(5)-1/2*u(9)./(x(1).^2+u(9)^2).*(2*u(9)*u(5) - 2*sign(x(1)).*(x(1).^4-x(1).^2*u(5)^2+u(9)^2*x(1).^2).^(1/2)))./x(1), 1/2./(x(1).^2+u(9)^2).*(2*u(9)*u(5) - 2*sign(x(1)).*(x(1).^4-x(1).^2*u(5)^2+u(9)^2*x(1).^2).^(1/2))))) + (abs(u(6)*x(4)-x(2))>u(8))*(-u(4)*(u(7)*u(1) - u(2)*sin(atan2((u(5)-1/2*u(9)./(x(1).^2+u(9)^2).*(2*u(9)*u(5) - 2*sign(x(1)).*(x(1).^4-x(1).^2*u(5)^2+u(9)^2*x(1).^2).^(1/2)))./x(1), 1/2./(x(1).^2+u(9)^2).*(2*u(9)*u(5) - 2*sign(x(1)).*(x(1).^4-x(1).^2*u(5)^2+u(9)^2*x(1).^2).^(1/2)))))*sign(x(2)-u(6)*x(4))))*u(6))/u(3);