%% Ecuaciones Euler-Lagrange / Hamilton de sistema de 2 masas unidas por 1 muelle
% © 2021, Antonio Sala Piqueras. Universitat Politécnica de València. CC-BY-SA 4.0
% Este código funcionó correctamente en Matlab R2021b
% Presentación en vídeo: http://personales.upv.es/asala/YT/V/twomELH.html
% Objetivo: comprender el significado de las distintas ecuaciones Euler-Lagrange y Hamilton en un sistema mecánico sencillo.
% Tabla de Contenidos
% Energías potencial y cinética
% Ecuaciones Euler-Lagrange
% Hamiltoniano, energía mecánica total
%
%

syms x1 x2 v1 v2 a1 a2 real
syms M1 M2 K real

%% Energías potencial y cinética

T=0.5*(M1*v1^2+M2*v2^2)%cinética
V=0.5*K*(x2-x1)^2 %potencial

%% Ecuaciones Euler-Lagrange

Lagrangian=T-V; %escalar
q=[x1;x2];
dotq=[v1;v2];
ddotq=[a1;a2];

p=jacobian(Lagrangian,dotq) %momentum, cant. movimiento

% Para cualquier \phi(q,\dot q)= \phi(x_1,x_2,v_1,v_2), tenemos: \frac{d\phi}{dt}=\frac{\partial \phi}{\partial q}\dot q+\frac{\partial \phi}{\partial \dot q}\ddot q

dpdt=jacobian(p,[q; dotq])*[dotq; ddotq]
syms F1 F2 real %fuerzas externas, actuan alineadas con cada GL.

% Ecuaciones Euler-Lagrange:    \frac{dp}{dt}=\frac{d}{dt}\left(\frac{\partial L}{\partial \dot q}\right)=\frac{\partial L}{\partial q} +Q

EcsMovEL=  dpdt == simplify( jacobian(Lagrangian,q)'+[F1; F2] ,50) 

% En simulación, se despejan las aceleraciones:

acels=solve(EcsMovEL,a1,a2);
acels.a1
acels.a2

%% Hamiltoniano, energía mecánica total
% *El enfoque Hamiltoniano no suele ser abordado en cursos de ingeniería, puede ser "saltado" por estudiantes que comienzan en la temática de modelado de sistemas físicos.

syms p1 p2 real %variables canónicas momentum
H0= -Lagrangian + [p1 p2]*dotq
v_from_p=solve(p==[p1 p2],dotq)
H=subs(H0,{v1,v2},{v_from_p.v1,v_from_p.v2}) %coords. canónicas (x,p)

% Ecuaciones de Hamilton
%        \frac{dp}{dt}=-\frac{\partial H}{\partial q} +Q,      \frac{dq}{dt}=\frac{\partial H}{\partial p}

dpdt=simplify(-jacobian(H,q)'+ [F1; F2],50)
dqdt=jacobian(H,[p1;p2])'

% Estas ecuaciones ya se podrían simular directamente: son una representación interna de primer orden con estados (x_1,x_2,p_1,p_2).
% Obviamente, equivalen a las de Euler Lagrange, dado que resultan en las mismas aceleraciones:

d2qdt_Hamilton=simplify( jacobian(dqdt,[p1;p2])*dpdt , 200)

% Aparte, Hamiltoniano y Energía son "lo mismo" en sistemas escleronómicos (sin grados de libertad "forzados" a funciones del tiempo):

Energia=T+V %x v
subs(H,{p1,p2},{p(1),p(2)})

% De modo que se conserva en ausencia de fuerzas externas:

dHdt=simplify(jacobian(H,[x1 x2 p1 p2])*[dqdt;dpdt],50) %en (x,p)
subs(dHdt,{p1,p2},{p(1),p(2)}) %en (x,v)