The 1D Wave Equation: Finite Difference Scheme

% matlab script waveeq1dfd.m
% finite difference scheme for the 1D wave equation 
% *fixed boundary conditions
% *raised cosine initial conditions

%%%%%% begin global parameters

SR = 44100;                                     % sample rate (Hz)
f0 = 200;                                       % fundamental frequency (Hz)
TF = 1;                                         % duration of simulation (s)
ctr = 0.4;                                      % center location of excitation (0-1)
wid = 0.1;                                      % width of excitation 
u0 = 0;                                         % maximum initial displacement
v0 = 1;                                         % maximum initial velocity
lambda = 1;                                     % Courant number
rp = 0.5;                                       % position of readout (0-1)

%%%%%% end global parameters

%%%%%% begin derived parameters

gamma = 2*f0;                                   % wave equation free parameter
k = 1/SR;                                       % time step
NF = floor(SR*TF);                              % duration of simulation (samples)
h = gamma*k/lambda;                             % grid spacing
N = floor(1/h);                                 % number of subdivisions of spatial domain
h = 1/N;                                        % reset h 
lambda = gamma*k/h;                             % reset Courant number
s0 = 2*(1-lambda^2); s1 = lambda^2;             % scheme parameters  
rp_int = 1+floor(N*rp);                         % rounded grid index for readout
rp_frac = 1+rp/h-rp_int;                        % fractional part of readout location

%%%%%% initialize grid functions and output

u = zeros(N+1,1); u1 = zeros(N+1,1); u2 = zeros(N+1,1); 
out = zeros(NF,1); 

%%%%%% create raised cosine

rc = zeros(N+1,1);
for qq=1:N+1
    pos = (qq-1)*h; dist = ctr-pos; 
    if(abs(dist)<=wid/2)
        rc(qq) = 0.5*(1+cos(2*pi*dist/wid));
    end
end

%%%%%% set initial conditions

u2 = u0*rc; u1 = (u0+k*v0)*rc; 

%%%%%% start main loop

for n=3:NF
    u(2:N) = -u2(2:N)+s0*u1(2:N)+s1*(u1(1:N-1)+u1(3:N+1));  % scheme calculation
    out(n) = (1-rp_frac)*u(rp_int)+rp_frac*u(rp_int+1);     % readout
    u2 = u1; u1 = u;                                        % update of grid variables
end

%%%%% end main loop

% plot output waveform

plot([0:NF-1]*k, out, 'k'); 
xlabel('t'); ylabel('u'); title('1D Wave Equation: Difference Scheme Output');
axis tight

% play sound

soundsc(out,SR);