%% second_order_sys.m Analysis of second order systems
% BJ Furman
% 12SEP2015
%% The 2nd Order System - Symbolic Solution
%
% The unforced equation of motion is:
%
%  Mx''(t) + cx'(t) + kx(t) = 0            (1)
%   or
%  x''(t) + 2*z*wn*x'(t) + wn^2*x(t) = 0    (2)   
%  where z == c/(2*M) the damping ratio
%  and wn == sqrt(k/m) the natural frequency
%
% The following solves the differential equation symbolically, and then
% uses the resulting time function, x(t) to generate a numeric solution 
clear all
syms z wn x(t)   % define the symbolic variables
xdot = diff(x,t);
xddot = diff(xdot,t);
x(t) = dsolve(xddot + 2*z*wn*xdot + wn^2*x(t) == 0, x(0)==2, xdot(0)==0)
z = 0.1;         % damping ratio
wn = 5;          % natural frequency
t = linspace(0,20,500);
res = eval(vectorize(x));  % vectorize insert .ops for elem-by-elem math
plot(t,res)
xlabel('Time,s')
ylabel('x(t)')
title('Unforced Response From Symbolic Result')
grid
%% The 2nd Order System - Numeric Solution 
%
% The unforced equation of motion is:
%
%  Mx''(t) + cx'(t) + kx(t) = 0            (1)
%   or
%  x''(t) + 2*z*wn*x'(t) + wn^2*x(t) = 0    (2)  
%
%    where z == c/(2*M) the damping ratio
%    and wn == sqrt(k/m) the natural frequency
%
% Note that:
%
%  x'(t) = dx/dt                           (3)
%  x''(t) = -(2*z*wn*x'(t) + wn^2*x(t))    (4)
%
% Matlab has a suite of numerical solvers for solving ODEs with initial 
% conditions. To use them for 2nd or higher order ODEs, one MUST first
% express the ODE as a set of first order ODEs in a Matlab function M-file.
%
% For example, given the ODE: x''' + ax'' + bx' + cx = d, let:
%
%  x   = x(1)                                  (5)
%  x'  = dxdts1 = x(2)                         (6)
%  x'' = dxdts2 = x(3)                         (7)
%  x''' = -( a*x(3) + b*x(2) + c*x(1)) + d     (8)
%
% then the set of 1st order ODEs that are needed will be eqs 6-8
% See the Matlab function M-file sec_ord.m that implements equations 6-8:
%
%   function dxdts = sec_ord(t,x,params)
%    dxdts is the column vector of the first order differential equations
%    We pass in the parameters for z and wn in the vector params
%      z = params(1);
%      wn = params(2);
%      dxdts = zeros(2,1);
%      dxdts(1) = x(2);
%      dxdts(2) = -(2*z*wn*x(2) + wn^2*x(1));
%
clear all
z = 0.1;
wn = 5;
params = [z wn];
t_start = 0;
t_end = 20;
t_span = [t_start t_end];
x0 = 2;
xdot0 = 0;
init_cond = [x0 xdot0];
% The call to the ode45 solver is below these comments. Note the argument list:
%
%   1. @name-of-function M-file containing the 1st order odes
%   2. vector consisting of the start time and end time for the integration
%   (alternatively you could pass a time vector with the times you want the
%   integration performed at. Must be in increasing or decreasing order)
%   3. vector consisting of the initial conditions x0, xdot0,..., x_n-1dot0
%   4. vector of strings defining options to control the integration. Just use
%   empty brackets, [], unless you have good reason to do otherwise.
%   5. Any parameters you wish to pass into the function with the 1st order
%   odes.
[t,x] = ode45(@sec_ord,t_span,init_cond,[],params);
plot(t,x)   % note that x is a matrix with x(:,1) position and x(:,2) velocity
xlabel('Time, s')
ylabel('Response')
title('Numerical Solution to Unforced Second Order ODE')
legend('x(t)','xdot(t)')
grid
%% Numerical Solution for Forced Response
clear all
close all
z = 0.1;
wn = 5;
params = [z wn];
t_start = 0;
t_end = 20;
% Let's define the integration times explicitly
t=linspace(t_start,t_end,500);
x0 = 0;
xdot0 = 0;
init_cond = [x0 xdot0];
F_t = 5;
[t,x] = ode45(@sec_ord_forced,t,init_cond,[],params,F_t);
plot(t,x)   % note that x is a matrix with x(:,1) position and x(:,2) velocity
xlabel('Time, s')
ylabel('Response')
title('Numerical Solution to Forced Second Order ODE')
legend('x(t)','xdot(t)')
grid
%% Numerical Solution for Forced Response That Is A Function of Time
% Reference: http://www.mathworks.com/matlabcentral/answers/97074-how-do-i-solve-an-ode-with-time-dependent-parameters-in-matlab
clear all
close all
z = 0.1;
wn = 5;
params = [z wn];
t_start = 0;
t_end = 20;
% Let's define the integration times explicitly
t=linspace(t_start,t_end,500);
x0 = 0;
xdot0 = 0;
init_cond = [x0 xdot0];
% define the time vector for the forcing function
ft=linspace(0,20,100);
% define the forcing function
%F_t = 5*ones(1,length(ft));
F_t = 5*sin(ft);
%
%   function dxdts = sec_ord_forced_F_t(t,x,ft,F_t,params)
%   dxdts is the column vector of the first order differential equations
%   We pass in the parameters for z and wn in the vector params
%   z = params(1);
%   wn = params(2);
%   F_t = interp1(ft,F_t,t);
%   dxdts = zeros(2,1);
%   dxdts(1) = x(2);
%   dxdts(2) = -(2*z*wn*x(2) + wn^2*x(1)) + F_t;
%
% call ode45 passing it the time vector for the integration, the time
% vector for the forcing function, the values of the forcing function, and
% the parameters for the odes
[t,x] = ode45(@(t,x) sec_ord_forced_F_t(t,x,ft,F_t,params),t,init_cond);
plot(t,x)   % note that x is a matrix with x(:,1) position and x(:,2) velocity
xlabel('Time, s')
ylabel('Response')
title('Numerical Solution to Second Order ODE With Time Dependent Forcing')
legend('x(t)','xdot(t)')
grid
%%