% Explicit finite difference scheme for vibrating string equation
% with damping:
%                 mu y_tt = T y_xx - alpha y_t.
% Note that this scheme becomes unstable if CFL condition is
% violated, i.e. if dt/hx is too large.

clear all

% Parameters
lx = 1;                   % length of string (m)
nx = 128;                 % number of grid points
tau = 150;                % string tension (N)
mu = 0.002;               % string mass (kg/m)
alpha = 0.05;             % velocity damping term
t_final = 0.1;            % time to stop simulation (s)
dt = 1/44100;             % length of time step (s)

% Compute some constants
hx = lx/(nx-1);           % length discretization
xj = [0:hx:lx]; 
dt2 = dt^2;
c2  = tau/mu;             % frequency constant
c1  = alpha/mu;           % damping constant

% Initial condition: string displaced 2 mm at x = lx/4.
y(1:nx/4) = 0.002*[0:nx/4-1]/(nx/4);
y(nx/4+1:nx) = 0.002*(1 - 4*[1:3*nx/4]/(3*nx));
yt = zeros(1,length(y));    % initial velocity zero

ym = y;
y = ym + dt*yt;

% Start stepping
for j = 1:floor(t_final/dt),
  t = j*dt;
  yp = (2*y - ym + dt2*c2*diff([-y(2), y, -y(nx-1)], 2)/(hx^2) + (c1*dt/2)*ym)/(1+(c1*dt)/2);
  ym = y; y = yp;
  if mod(j,10) == 0
    plot(xj,y); axis([0 lx -0.002 0.002]); drawnow;
  end
  ytrace(j) = y(nx/4);
  tj(j) = t;
end

plot(tj,ytrace)