%% Solving Linear Ordinary Differential Equations with Matlab's symbolic toolbox: spring-mass-damping example
%%
% © 2022, Antonio Sala Piqueras, Universitat Politècnica de València. All rights reserved.
% This code executed in Matlab R2022a
% Video-presentations:
% http://personales.upv.es/asala/YT/V/masmusym1EN.html , http://personales.upv.es/asala/YT/V/masmusym2EN.html
% Objectives: analyzing the solution of linear ODEs provided by Matlab (symbolic toolbox), understanding its meaning, and representing it graphically in a spring-damper (free response) model, and understanding a very simple animation code.
%
% Table of Contents
% EDO in non-normalized (second derivatives) form
%     General solution
%     Particular solution for given intial conditions (position and speed)
% ODE in normalised form (dx/dt=Ax)
%     General solution x=[pos; vel];
%     Particular solution for given initial conditions
%
%% EDO in non-normalized (second derivatives) form
% Consider the following spring (with some damping, too, by friction with walls and spring losses):
%
% Let us solve the equation of a mass-spring-damper system (linearized, zero equilibrium; linear viscous friction assumed), starting at initial conditions out of equilibrium, written as:
% M\frac{d^2 y}{dt^2}=-ky-b\frac{dy}{dt}

syms y(t) %symbolic function of time
syms M k b real %constant parameters
vel=diff(y); % notation for velocity
accel=diff(y,2); % notation for acceleration

ODE_SYM= M*accel==-k*y-b*vel
%if it were a nonlinear ODE, a solution could not, in general, be found in the form of a symbolic expression
M_num=0.5;k_num=1;b_num=0.2; %numerical values for constant parameters
%let's replace symbols by actual numbers
ODE_mass_spring=subs(ODE_SYM, {M,k,b}, {M_num,k_num,b_num})

%% General solution
% With no fixed initial conditions, there are infinitely many solutions, formally related to two integration constants:

sol=simplify(dsolve(ODE_mass_spring))

% This is the so-called general solution, where C_1 and C_2 influence initial conditions.
% Decay rate: 1/5 s^{-1}, free oscillation frequency: 7/5 rad/s.

sol_velocity=simplify(diff(sol))

% Let's test that it verifies the ODE:

simplify( M_num*diff(sol,2) == -k_num*sol-b_num*diff(sol) )
subs(sol,t,0)%initial position
subs(sol_velocity,t,0)%initial velocity

%% Particular solution for given intial conditions (position and speed)

position=dsolve(ODE_mass_spring, y(0)==3, vel(0)==2)

% Let's check that the prescribed initial conditions are indeed met:

subs(position,t,0)
vel=simplify(diff(position))
subs(vel,t,0)

% If we draw position (thick blue) and speed (dashed red), we get:

Tf=25; %the final time at which we wish to stop plotting
fplot(position,[0 Tf],LineWidth=3), hold on
fplot(vel,[0 Tf],LineStyle='-.'), hold off, grid on,
legend(["Position","Velocity"], location='best')

%
%% ODE in normalised form (dx/dt=Ax)
% We pursue a model in the form \frac{dx}{dt}=Ax (no input, so no need of \frac{dx}{dt}=Ax+Bu for the moment being).
% We check that we get coincident results.

syms pos(t) vel(t) %position and velocity

% We express the mass-spring-damping model as a pair of first-order differential equations:

ODE_SYM= [ diff(pos) == vel;   diff(vel) == -k/M*pos-b/M*vel ]

% We can express the same thing in matrix form (because the equations are linear):

x=[pos; vel]; %state fector
A=[0 1;-k/M -b/M];
ODE_SYM=    diff(x) == A*x %normalised matrix form representation of the ODE

% Matrix form is convenient for subsequent theoretical analysis (not in the scope of this material).
%% General solution x=[pos; vel];
% We will solve with dsolve command, once the numerical values of the constant parameters have been replaced:

ODE_mass_spring=subs(ODE_SYM, {M,k,b}, {M_num,k_num,b_num})
A_numeric=subs(A, {M,k,b}, {M_num,k_num,b_num}) %this is the matrix A whose eigenvalues, etc. are theoretically important
sol=dsolve(ODE_mass_spring)
simplify(sol.pos)
simplify(sol.vel)

% It's the general solution. Initial conditions would be

subs(sol,t,0)

%
%% Particular solution for given initial conditions

spring_state=dsolve(ODE_mass_spring, pos(0)==3, vel(0)==2);
simplify(spring_state.pos)
simplify(spring_state.vel)

% Prescribed initial conditions do hold:

subs(spring_state,t,0)

% Let us plot the resulting time response (of course, identical to the one at first section):

fplot(spring_state.pos,[0 Tf],LineWidth=3), hold on
fplot(spring_state.vel,[0 Tf],LineStyle='-.'), hold off, grid on, legend(["Position","Velocity"], location='best')