##================================================================
## PROBLEM DESCRIPTION
##================================================================
## Mass particles are positioned at points 1, 2, ..., n of the vertical axis,
## pointing down.
## Assume m(i) = is the mass of the i-th particle.
## The masses are connected with springs. The top mass is attached to
## the origin with a spring. The bottom mass is dangling freely.
## Hookes law is assumed to be in effect, and the hook constant is hook=2.0. 
## It is assumed that the system is in equilibrium.
##
## Now, let us suppose that the initial velocity of the masses is
## 0 and the displacements from the equilibrium are x(0), x(1), ..., x(n).
## Find the position and velocity of the masses after time t.
##================================================================
## END OF PROBLEM DESCRIPTION
##================================================================

##================================================================
## ANALYSIS
##================================================================
## Let us assume that each x(i) depends on t, but we will not
## explicitly carry this dependence in our formulas.
##
## The force acting on the i-th particle is Hooke's constant multiplied by:
##
##   (x(i-1)-x(i)) + (x(i+1)-x(i)) = x(i+1)+x(i-1)-2*x(i)
##
## The exceptions are:
## i = 1;  -x(i) + (x(i+1)-x(i)) = x(i+1) - 2 * x(i)
## i = n;  x(i-1)-x(i) = 
##
## Newton's law for the i-th particle is:
##
##   d^2 x(i) / dt^2 = hooke*(x(i+1)+x(i-1)-2*x(i))
##
## with the exceptions mentioned before.
## We write the system in the matrix form:
##     d^2 x / dt^2 = A * x
##================================================================
## END OF ANALYSIS
##================================================================

##
## Explanation of args:
## t      - end time of the simulation
## x0     - vector of initial displacements
## v0     - vector of initial velocities (defaults to zeros(n))
## m      - vector of masses (defaults to ones(n))
## h      - Hooke's constant (defaults to 1)
## 
function [x,v] = springs(t, x0, v0, m, h)
  if(nargin < 2)
    error("Must provide time and initial displacement.");
  endif
  n = length(x0);
  if(nargin < 3)
    v0 = zeros(n, 1);
  endif
  if(nargin < 4)
    m = ones(n, 1);
  endif
  if(nargin < 5)
    h = 1;
  endif
  A = diag(ones(n - 1, 1), 1) + diag(ones(n - 1, 1), -1) \
      + diag([-2*ones(n-1, 1); -1]);
  A = h * A;
  A = inv(diag(m)) * A;
  B = sqrtm(-A);
  x = cosm(t * B) * x0 + inv(B) * sinm(t * B) * v0;
  v = - B * sinm(t * B) * x0 + cosm(t * B) * v0;
endfunction

% ;;; Local Variables: ***
% ;;; mode: octave ***
% ;;; End: ***
