#!/usr/bin/octave -q

%-----------------------------------------------------------------------------
% Jonathan R. Senning <jonathan.senning@gordon.edu>
% Gordon College
% December 1, 2000
% May 1, 2003  Revised to include SOR code and provide proper comments.
% October 21, 2008 Revised to correct PDE problem
% November 25, 2008 Revised to add optional graphics commands
%
% Use finite differences to solve the elliptic Helmholtz problem on a square
% domain.  The equation solved is
%
%       u_xx + u_yy + f(x,y)u = g(x,y)
%
% with f(x,y) = -1/25 and g(x,y) = 0 and boundary conditions given by the
% known solution u(x,y) = cosh(x/5) + cosh(y/5).  This problem and solution
% is given in "Numerical Mathematics and Computing", 6th edition, by Cheney
% and Kincaid, Thompson Brooks/Cole, 2008.
%-----------------------------------------------------------------------------

% Get grid size from command line
arg_list = argv ();
if ( nargin == 0 )
  n = 20;   % default
  m = 20;
elseif ( nargin == 1 )
  n = str2double( arg_list{1} );
  m = n;
else
  n = str2double( arg_list{1} );
  m = str2double( arg_list{2} );
end

% Set the size of the problem.  (Domain is the unit square)
xmin = 0.0;
xmax = 1.0;
ymin = 0.0;
ymax = 1.0;

h = ( xmax - xmin ) / n;
h2 = h * h;

% Set the SOR parameter omega and a tolerance
w = 2.0 / ( 1.0 + 0.5 * ( sin( pi / n ) + sin( pi / m ) ) );
tol = 1e-6;

% Make sure that the step size is the same in both x and y directions
if ( abs( m * h - ( ymax - ymin ) ) > tol )
  fprintf( stderr, 'Step size in x and y directions must be equal\n' );
  break
end

% Allocate array for grid of u(x,y) values and determine the x and y values
% at each grid point.
u = ones( n+1, m+1 );
x = linspace( xmin, xmax, n+1 );
y = linspace( ymin, ymax, m+1 );

% Assign the boundary values.  In this case we do this by using the known
% solution to this problem.
solution = @( x, y ) ( cosh( x / 5 ) + cosh( y / 5 ) );

u(:,1)   = solution( x, ymin )';
u(:,m+1) = solution( x, ymax )';
u(1,:)   = solution( xmin, y );
u(n+1,:) = solution( xmax, y );

% Fill the arrays holding the values of the functions f(x,y) and g(x,y) from
% the Helmholtz equation
f = -0.04 * ones( n+1, m+1 );
g = zeros( n+1, m+1 );

% Begin SOR iterations.  Continue until the norm of the difference of
% consecutive estimates is less than tol.
u0 = u .+ tol;
k = 1;
while ( norm( u - u0, inf ) >= tol )
  u0 = u;
  for j = 2 : m
    for i = 2 : n
      u(i,j) = ( 1 - w ) * u(i,j) ...
          + w * ( u(i,j-1) + u(i-1,j) + u(i+1,j) + u(i,j+1) - h2 * g(i,j) ) ...
                / ( 4 - h2 * f(i,j) );
    end
  end
  fprintf( '*' );
  fflush( stdout );
  k = k + 1;
end
fprintf( '\n' );

% Verify solution
sum = 0.0;
for j = 2 : m
  for i = 2 : n
    sum = sum + ( u(i,j) - solution( (i-1)*h, (j-1)*h ) )^2;
  end
end
fprintf( 'Exact Error Norm = %e\n', sqrt( sum ) );

display_graphics = 1;
if ( display_graphics == 1 )
   % Display an image showing the solution.
   colormap( 'jet' )
   ncolor = size( colormap )(1)
   umin = min( min( u ) );
   umax = max( max( u ) );

   image( x, y, ncolor * ( u' - umin ) / ( umax - umin ) );
   xlabel( 'x' );
   ylabel( 'y' );
   drawnow();
   input( 'Press Enter to continue... ' );

   % Display a surface plot showing the solution.
   surf( x, y, u' );
   xlabel( 'x' );
   ylabel( 'y' );
   drawnow();
   input( 'Press Enter to continue... ' );

   % Display a contour plot showing the solution.
   levels = linspace( umin, umax, 21 );
   contour( x, y, u', levels );
   xlabel( 'x' );
   ylabel( 'y' );
   drawnow();
   input( 'Press Enter to continue... ' );
end

% End of file