function poisson_2d ( nx, ny, rhs, u_exact, uxxyy_exact )

%*****************************************************************************80
%
%% poisson_2d() computes a solution to the Poisson equation on a rectangle.
%
%  Discussion:
%
%    This program runs serially.  Its output is used as a benchmark for
%    comparison with similar programs run in a parallel environment.
%
%    The Poisson equation
%
%      - DEL^2 U(x,y) = F(x,y)
%
%    is solved on the unit square [0,1] x [0,1] using a grid of NX by
%    NX evenly spaced points.  The first and last points in each direction
%    are boundary points.
%
%    The boundary conditions and F are set so that the exact solution is
%
%      U(x,y) = sin ( pi * x * y )
%
%    so that
%
%      - DEL^2 U(x,y) = pi^2 * ( x^2 + y^2 ) * sin ( pi * x * y )
%
%    The Jacobi iteration is repeatedly applied until convergence is detected.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    03 October 2024
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    integer NX, NY, the number of nodes in the X and Y directions.
%
%    f = rhs ( x, y ): evaluates the right hand side function.
%
%    u = u_exact ( x, y ): evaluates the exact solution.
%
%    uxxyy = uxxyy_exact ( x, y ): evaluates the exact uxxyy.
%
%  Local:
%
%    real DX, DY, the X and Y spacing.
%
%    integer IT_MAX: the maximum number of Jacobi iterations allowed.
%
%    real TOLERANCE: convergence is declared if successive iterates differ
%    by less than this value in norm.
%
  a = 0.0;
  b = 1.0;
  c = 0.0;
  d = 1.0;
  dx = ( b - a ) / ( nx - 1 );
  dy = ( d - c ) / ( ny - 1 );
  it_max = 1000;
  tolerance = 0.000001;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'poisson_2d():\n' );
  fprintf ( 1, '  A program for solving the Poisson equation:\n' );
  fprintf ( 1, '    -DEL^2 U = F(X,Y)\n' );
  fprintf ( 1, '  on the rectangle 0 <= X <= 1, 0 <= Y <= 1.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  F(X,Y) = pi^2 * ( x^2 + y^2 ) * sin ( pi * x * y )\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The number of X grid points is %d\n', nx );
  fprintf ( 1, '  The number of Y grid points is %d\n', ny );
  fprintf ( 1, '  The X grid spacing is %f\n', dx );
  fprintf ( 1, '  The Y grid spacing is %f\n', dy );
%
%  Set the right hand side.
%
  f = rhs ( nx, ny );
  fnorm = norm ( f, 'fro' ) / nx / ny;
  fprintf ( 1, '  RMS of F = %g\n', fnorm );
%
%  Set the initial solution estimate.
%  We are "allowed" to pick up the boundary conditions exactly.
%
  unew = zeros ( nx, ny );

  unew(1,   1:ny) = f(1,   1:ny);
  unew(  nx,1:ny) = f(  nx,1:ny);
  unew(1:nx,1)    = f(1:nx,1);
  unew(1:nx,  ny) = f(1:nx,  ny);

  unew_norm = norm ( unew, 'fro' ) / nx / ny;;
%
%  Set up the exact solution.
%
  x = linspace ( a, b, nx );
  y = linspace ( c, d, ny );
  [ X, Y ] = meshgrid ( x, y );
  uexact = u_exact ( X, Y );
  unorm = norm ( uexact, 'fro' ) / nx / ny;
  fprintf ( 1, '  RMS of exact solution U = %g\n', unorm );
%
%  Do the iteration.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'Step    |Unew|     |Unew-U|     |Unew-Exact|\n' );
  fprintf ( 1, '\n' );

  error = norm ( unew - uexact, 'fro' ) / nx / ny;
  fprintf ( 1, '%4d  %12g                %12g\n', 0, unew_norm, error );

  for it = 1 : it_max

    u = unew;
%
%  Perform one Jacobi sweep.
%
    unew = jacobi ( nx, ny, dx, dy, f, u );
%
%  Check for convergence.
%
    u_norm = unew_norm;
    unew_norm = norm ( unew, 'fro' ) / nx / ny;
    diff = norm ( unew - u, 'fro' ) / nx / ny;
    error = norm ( unew - uexact, 'fro' ) / nx / ny;

    fprintf ( 1, '%3d  %12g  %12g  %12g\n', it, unew_norm, diff, error );

    if ( diff <= tolerance )
      fprintf ( 1, '  The iteration has converged.\n' );
      return
    end

  end

  fprintf ( 1, '  The iteration has NOT converged.\n' );

  return
end
function unew = jacobi ( nx, ny, dx, dy, f, u )

%*****************************************************************************80
%
%% jacobi() carries out one step of the Jacobi iteration.
%
%  Discussion:
%
%    Assuming DX = DY, we can approximate
%
%      - ( d/dx d/dx + d/dy d/dy ) U(X,Y) 
%
%    by
%
%      ( U(i-1,j) + U(i+1,j) + U(i,j-1) + U(i,j+1) - 4*U(i,j) ) / dx / dy
%
%    The discretization employed below will not be correct in the general
%    case where DX and DY are not equal.  It's only a little more complicated
%    to allow DX and DY to be different, but we're not going to worry about 
%    that right now.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    25 October 2011
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    integer NX, NY, the X and Y grid dimensions.
%
%    real DX, DY, the spacing between grid points.
%
%    real F(NX,NY), the right hand side data.
%
%    real U(NX,NY), the previous solution estimate.
%
%  Output:
%
%    real UNEW(NX,NY), the updated solution estimate.
%

%
%  The boundary values are stored in F.  Copy them.
%
  unew = f;
%
%  The interior values are set by a Jacobi iteration.
%
  unew(2:nx-1,2:ny-1) = 0.25 * ( ...
      u(1:nx-2,2:ny-1) ...
    + u(2:nx-1,3:ny) ...
    + u(2:nx-1,1:ny-2) ...
    + u(3:nx,  2:ny-1) ...
    + f(2:nx-1,2:ny-1) * dx * dy );

  return
end

