% msd_fd.m
%
% Used to plot the solution of a mass-spring-damper system using
% a finite difference approach.
%
% Gary P. Scavone
% McGill University, 2004.

% Simulation Parameters

fs = 44100;      % sample rate
T = 1/fs;        % sample period
N = 1000;        % iterations

% System Parameters

m = 0.5;          % mass in kg
k = 5000000;      % spring constant (N/m)
w0 = sqrt(k/m);   % natural angular frequency of system
f0 = w0/(2*pi);   % natural frequency of system (no damping)
mu = 50;          % damping constant
x0 = 0.01;        % initial displacement
v0 = 0;           % initial velocity

% System Equation: m(x'') + mu(x') + kx = 0, where x' = dx/dt

% Resulting difference equation after substituting the finite
% difference approximation to the derivative:
%
% x[n] = a x[n-1] - b x[n-2]
%
% where:
d = m + mu*T + k*T^2;
a = (2*m + mu*T) / d;
b = -m / d;

% Displacement Initialization
x = zeros(1,N);
xm1 = x0;      % x[n-1]
xm2 = xm1;     % x[n-2]

% Now do simulation:
for n = 1:N
  
  x(n) = a*xm1 + b*xm2;
  xm2 = xm1;
  xm1 = x(n);
  
end

plot(x)
grid
xlabel('Time Steps (samples)');
ylabel('Mass Displacement (m)');
title('Finite Difference');
