wave.m

function wave(bc,L,c,n,dt,figno)

%
%    wave(bc,L,c,n,dt,figno)
%
%  Solve one-dimensional wave equation
%
%      u_{tt} = c^2 u_{xx}
%
%   using finite difference numerics.
%
%  bc =   'f'  -- Dirichlet boundary conditions: both ends fixed  (default if omitted)
%         'n'  -- Neumann boundary conditions: both ends free
%         'p'  -- periodic boundary conditions
%
%  L  = length of interval 0 <= x <= L  (default = pi)
%
%  c  = wave speed  (default = f)
%
%  n  = number of space mesh points  (default = 65)
%
%  dt = time step size  (default = .005)
%
%  figno = Figure number  (default = 1)
%
%  The default initial conditions are
%
%     u(0,x)   = exp(-40*(x-1).^2)
%     u_t(0,x) = 0
%
%  These can be changed by clicking 
%     the "init. data" button on the graphics.

if nargin < 1, bc = 'f'; end
if nargin < 2, L = pi; end
if nargin < 3, c = 1; end
if nargin < 4, n = 64; end
if nargin < 5, dt = .005; end
if nargin < 6, figno = 1; end

global buttonid;

buttonid=0;

icons = ['text(.5,.5,''pause''     ,''Horiz'',''center'')'
         'text(.5,.5,''show time'' ,''Horiz'',''center'')'
	 	 'text(.5,.5,''restart''   ,''Horiz'',''center'')'
         'text(.5,.5,''init. data'',''Horiz'',''center'')'
         'text(.5,.5,''input''     ,''Horiz'',''center'')'
	 	 'text(.5,.5,''stop''      ,''Horiz'',''center'')'];
callbacks = ['button(1)'
             'button(2)'
             'button(3)'
   	    	 'button(4)'
             'button(5)'
             'button(6)'];
presstype =  ['toggle'
              'flash '
              'flash '
              'flash '
              'flash '
	    	  'flash '];

hf = figure(figno); 

bg = btngroup(hf,'GroupID','b','ButtonId',['p';'s';'r';'0';'d';'i'],...
            'IconFunctions',icons,'GroupSize',[1 6],...
	    'CallBack',callbacks,'Position',[.1 0 .8 .07],...
            'PressType',presstype);

h = L/n;

l = dt*c/h; l2 = l^2;

M = l2 *( diag(ones(n,1),1) + diag(ones(n,1),-1));
switch bc 
	case 'n'
		M(1,2) = 2*l2; M(n+1,n) = 2*l2;
	case 'p'
		M(1,n+1) = l2; M(n+1,1) = l2;
end

x = linspace(0,L,n+1)';

x0 = 1;
f0 = 0*x;
del = f0; del(25) = 1;
f1 = exp(-40*(x-x0).^2);
f2 = (x < 1) .*x + (x>1 & x < 2.8) .* (3 - 2*x) + (x > 2.8) .* (x - pi) * 2.6 / (pi - 2.8); 
f3 = 1 - abs(2*x/L - 1) + .1 * sin(15*x);

f = f1; 
g = f0;
      
clf;
ax = axes('Position',[.15 .17 .75 .75]);
A = 2;
t=0; 

sol= plot(x,f,'EraseMode','background');
axis([0 L -A A]);
w = f;
z = (1 - l2) * f + M * f/2  + dt * g;
set(sol,'ydata',z);
t = t + dt;

stopit = 0;

while stopit == 0 
  switch buttonid
   case 1  
      title(['Time: ',num2str(t)]);
	  buttonid = 0;
      while buttonid == 0, pause(.2); end 
	  buttonid = 0;
   case 2
      title(['Time: ',num2str(t)]);
	  buttonid = 0;
   case 3
      clf;
	  sol= plot(x,f,'EraseMode','background');
	axis([0 L -A A]);
	w = f;
	z = (1 - l2) * f + M * f/2  + dt * g;
	set(sol,'ydata',z);
	t = t + dt;
	buttonid = 0;
   case 4
      title(['Time: ',num2str(t)]);
	  f = input('Initial displacement: ');
      g = input('Initial velocity: ');
      clf;
	  sol= plot(x,f,'EraseMode','background');
	axis([0 L -A A]);
	w = f;
	z = (1 - l2) * f + M * f/2  + dt * g;
	set(sol,'ydata',z);
	t = t + dt;
	buttonid = 0;
   case 5
      title(['Time: ',num2str(t)]);
	  keyboard
	  buttonid = 0;
   case 6
	  stopit = 1;
      title(['Time: ',num2str(t)]);
	  buttonid = 0;
  end

	u = 2 * (1 - l2) * z + M * z - w;
	w = z;
	z = u;
	set(sol,'ydata',z);drawnow; 

	t=t+dt;

end