%% Solving Linear Ordinary Differential Equations with Matlab's symbolic toolbox: spring-mass-damping example (forced response)
% © 2022, Antonio Sala Piqueras, Universitat Politècnica de València. All rights reserved.
% This code executed in Matlab R2022a
% Video-presentation: http://personales.upv.es/asala/YT/V/masmusymForzEN.html
% Objectives: analyzing the solution of linear ODEs provided by Matlab (symbolic toolbox dsolve), understanding its meaning, and graphically representing it  in a spring-damper model, in this case with a forced response.
%
% Table of Contents
% EDO in non-normalized (second derivatives) form
%     General solution
%     Particular solution for given initial conditions (position and speed) and input
%
%% 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}+F(t)

syms y(t) F(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+F(t)
%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=2;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 and some integrals (convolutions) of force input:

sol=simplify(dsolve(ODE_mass_spring))

%% Particular solution for given initial conditions (position and speed) and input

ODE_mass_spring=subs(ODE_mass_spring,F(t),-2.5+6*sin(3*t));
position=simplify(dsolve(ODE_mass_spring, y(0)==3, vel(0)==0))

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

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

% Transients take more or less e^{-\frac{t_{est}}{5}}\approx 0.02, so:

t_est=log(0.02)*(-5)

% The system's proper oscillation frequency is:

omega_p=3*sqrt(11)/5

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

Tf=30; %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,
xline(t_est,'g')
legend(["Position","Velocity"], location='best')

%%