%
%   Solution of some linear dispersive wave PDEs
%
%   Andre Weideman                             5/10/95
%                                  Revised:    2/20/98
%

%  Input:

disp('*** Solution of Linear Dispersive Wave PDEs ***')
disp(' ')
disp('    (1) Klein-Gordon  ')
disp('    (2) Korteweg-de Vries  ')
disp('    (3) Benjamin-Ono ')
disp('    (4) Schroedinger ')
disp('    (5) Regularized Long Wave ')
PDE = input(' Which one? (1)--(5) ');


u = input(' Initial condition (string function of x, in single quotes) = ? ');
                                                 % Example:  'sech(x)'
if PDE == 1; disp(' Note: For KG equation u_t(x,0) = 0 is assumed'); end
L = input(' Space interval [-L,L], L = ? ');     % 2L = interval length
t = input(' Time t = ? ');                       % t = time 

%  Set up the x-grid and wave numbers

N = 2^10;                 % N = number of grid points; the less smooth the 
                          % initial condition, the larger N needs to be.
n = [-N:N-1];             
h = L/N;                  % h = grid spacing
x = n*h;                  % x = grid
k = n*pi/L;               % k = wave numbers representable on this grid

%  Define dispersion relation

if PDE == 1
w = sqrt(k.^2+1);         % Klein-Gordon
elseif PDE == 2
w = k.^3;                 % Korteweg-de Vries 
elseif PDE == 3
w = sign(k).*k.^2;        % Benjamin-Ono
elseif PDE == 4 
w = k.^2;                 % Schroedinger
elseif PDE == 5
w = k./(1+k.^2);          % Regularized Long Wave
end

u0 = eval(u);

% Take Fourier transform of initial data

a0 = fftshift(fft(u0));

% Advance in Fourier space

if PDE == 1
at = a0.*cos(w*t);        % Klein-Gordon is second order; u_t(x,0)=0 is assumed
else
at = a0.*exp(-i*w*t);     % The others are first order w.r.t. time.
end

% Return to physical space via inverse transform

ut = ifft(fftshift(at));

% Plot the solution

plot(x,real(ut));











