%% Comparación de solvers (ode45 vs ode15s) según el problema sea stiff/non-stiff
% © 2022, Antonio Sala Piqueras, Universitat Politècnica de València. Todos los derechos reservados.
% Este código ejecutó en Matlab R2021b
% Presentación en vídeo: http://personales.upv.es/asala/YT/V/ode45vs15s.html
% Objetivo: Los sistemas "stiff" (rígidos) son aquéllos con constantes de tiempo (polos de linealización) muy diferentes entre sí. En algunos casos, puede convenir solver diferente.
% Tabla de Contenidos
% Planteamiento del problema
% Integración numérica, comparación de solvers
%     Con ode45
%     Con ode15s
% Conclusiones
%% Planteamiento del problema
% Modelo (lineal), pero la idea aplica a no-lineal

A=[-1 0.05 0;0 -1.25*20 0.05;0.0 0 -1.5*4000]
eig(A)
odefun=@(t,x) A*x+[1; 2+2*cos(t/10); 3]*sin(t)^3;

% Condición inicial, tiempo final y opciones de tolerancia:

x0=[1;-1;1];
tf=500;
opts=odeset(AbsTol=1e-4,RelTol=1e-5);

%% Integración numérica, comparación de solvers
%% Con ode45

tic
[T,X]=ode45(odefun,[0 tf],x0,opts);
toc
size(T)

%% Con ode15s

tic
[T2,X2]=ode15s(odefun,[0 tf],x0,opts);
toc
size(T2)
plot(T,X), hold on
plot(T2,X2), hold off, grid on

%% Conclusiones
% ode45 es el "básico", que Matlab recomienda como primera opción a probar en problemas "suaves" y con constantes de tiempo todas "parecidas"... Pero con constantes de tiempo muy diferentes, tarda mucho. El solver ode15s tarda más en el caso "non-stiff" pero menos en el "stiff".
% En cálculo numérico con modelos de orden elevado (incluso con decenas de miles de variables de estado, elementos finitos), diferentes constantes de tiempo, discontinuidades en las ecuaciones de estado (impactos, singularidades, entradas o condiciones iniciales impulsivas...) elegir el solver más adecuado podría requerir un estudio cuidadoso por un experto en cálculo numérico y software asociado antes de correr una simulación que tarde "dos semanas" en dar un "mal resultado".
% Matlab tiene todas estas alternativas para resolver ecuaciones diferenciales ordinarias:
% Solver Accuracy When to Use
% ode45 Nonstiff Medium
% Most of the time. ode45 should be the first solver you try.
% ode23 Low
% ode23 can be more efficient than ode45 at problems with crude tolerances, or in the presence of moderate stiffness.
% ode113 Low to High
% ode113 can be more efficient than ode45 at problems with stringent error tolerances, or when the ODE function is expensive to evaluate.
% ode78 High
% ode78 can be more efficient than ode45 at problems with smooth solutions that have high accuracy requirements.
% ode89 High
% ode89 can be more efficient than ode78 on very smooth problems, when integrating over long time intervals, or when tolerances are especially tight.
%
% ode15s Stiff Low to Medium
% Try ode15s when ode45 fails or is inefficient and you suspect that the problem is stiff. Also use ode15s when solving differential algebraic equations (DAEs).
% ode23s Low
% ode23s can be more efficient than ode15s at problems with crude error tolerances. It can solve some stiff problems for which ode15s is not effective.
% ode23s computes the Jacobian in each step, so it is beneficial to provide the Jacobian via odeset to maximize efficiency and accuracy.
% If there is a mass matrix, it must be constant.
% ode23t Low
% Use ode23t if the problem is only moderately stiff and you need a solution without numerical damping.
% ode23t can solve differential algebraic equations (DAEs).
% ode23tb Low
% Like ode23s, the ode23tb solver might be more efficient than ode15s at problems with crude error tolerances.
% ode15i Fully implici tLow
% Use ode15i for fully implicit problems f(t,y,y’) = 0 and for differential algebraic equations (DAEs) of index 1.
%