%% Exponencial de Matrices, respuesta temporal dx/dt=Ax
% © 2019, Antonio Sala Piqueras. Universitat Politècnica de València. Todos los derechos reservados.
% Vídeo-presentación en http://personales.upv.es/asala/YT/V/rtssml.html
%
%% 1. Cálculo numérico de la exponencial: representación gráfica de respuesta libre

%matriz del sistema dxdt=Ax:
A=[-4.5 -2 3;-4 -5 -0.5;-6 -7 -8];
%condición inicial para simulaciones (arbitraria):
x0=[1;0.2;0.6];

% Simulación exacta (NO hay integración numérica Euler/Runge-Kutta)

tiempos=0:0.02:3;
Nmuestras=length(tiempos);
estados=zeros(3,Nmuestras);
for muestra=1:Nmuestras
    estados(:,muestra)=expm(A*tiempos(muestra))*x0;
end
plot(tiempos,estados,'LineWidth',2)
hold on
plot(0,x0,'o','LineWidth',5)
hold off
grid on
title('Respuesta temporal (variables estado internas)')

% Simulación de una ecuación de salida: y=Cx

C=[1 0 -1];
plot(tiempos,C*estados)
grid on
title('Respuesta temporal (variable de salida)')

%% 2. Cálculo simbólico: fórmula explícita de la respuesta temporal

syms lambda omega alpha t real
%syms t real
jojo=[alpha -omega;omega alpha]
eig(jojo)
simplify(expm(jojo*t),'Criterion','preferReal','Steps', 40)
J=[lambda 1 0;0 lambda 1;0 0 lambda]
expm(J*t)

%


%Matlab sabe hacerla esta exponencial... o casi
tic,otracosa=simplify(expm(A*t)),toc

%

%con un poco de ayuda tarda mucho menos y está todo más claro:
[V,D]=eig(A);
D
Vs=vpa(V,6);Ds=vpa(D,6);
tic,exponencialdeAt=simplify(Vs*expm(Ds*t)*inv(Vs)),toc
vpa(exponencialdeAt*x0,4) %salen todo números reales excepto errores de redondeo:

% Ecuación de salida:

vpa(C*exponencialdeAt*x0,4)