function exercise2 ( )

%*****************************************************************************80
%
%% exercise2 estimates the Laplacian on a 2D periodic domain.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 June 2020
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'exercise2\n' );
  fprintf ( 1, '  Estimate the Laplacian on a 2D periodic region.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  nx  ny  || error ||\n' );
  fprintf ( 1, '\n' );

  for n = [ 7, 13, 25, 49 ]

    nx = n;
    xmin = 0.0;
    xmax = 2.0;
    x = linspace ( xmin, xmax, nx );
    dx = ( xmax - xmin ) / ( nx - 1 );

    ny = n;
    ymin = -1.0;
    ymax =  1.0;
    y = linspace ( ymin, ymax, ny );
    dy = ( ymax - ymin ) / ( ny - 1 );
%
%  Create periodic grid [ XP, YP ], by reducing nx and ny and removing last points.
%
    nxp = nx - 1;
    nyp = ny - 1;
    xp = x(1:nxp);
    yp = y(1:nyp);
    [ XP, YP ] = ndgrid ( xp, yp );
%
%  Compute function and exact Laplacian.
%
    UP = sin ( 2.0 * pi * XP ) .* sin ( 2.0 * pi * YP );
    UXXP = - 8.0 * pi^2 * sin ( 2.0 * pi * XP ) .* sin ( 2.0 * pi * YP );
%
%  Compute discrete Laplacian and compare.
%
    LP = laplacian5_torus ( UP, dx, dy );
    diff = rms_norm ( UXXP - LP );
    fprintf ( 1, '  %2d  %2d  %g\n', nx, ny, diff );

  end

  return
end
function L = laplacian5_torus ( A, dx, dy )

%*****************************************************************************80
%
%% laplacian5_torus uses a 5 point Laplacian stencil on a torus.
%
%  Modified:
%
%    28 June 2020
%

%
%  Vectorized version:
%
  if ( true )

    f_xx = (        circshift ( A, [  1,  0 ] ) ... 
           - 2.0 .* circshift ( A, [  0,  0 ] ) ...
                  + circshift ( A, [ -1,  0 ] ) ) / dx^2;

    f_yy = (        circshift ( A, [  0,  1 ] ) ...
           - 2.0 .* circshift ( A, [  0,  0 ] ) ...
                  + circshift ( A, [  0, -1 ] ) ) / dy^2;

    L = f_xx + f_yy;
%
%  Scalar version.
%
  else

    [ nxp, nyp ] = size ( A );
    L = zeros ( nxp, nyp );

    for j = 1 : nyp

      jp1 = j + 1;
      if ( nyp < jp1 )
        jp1 = 1;
      end
      jm1 = j - 1;
      if ( jm1 < 1 )
        jm1 = nyp;
      end

      for i = 1 : nxp

        ip1 = i + 1;
        if ( nxp < ip1 )
          ip1 = 1;
        end
        im1 = i - 1;
        if ( im1 < 1 )
          im1 = nxp;
        end

        L(i,j) = ( A(ip1,j) - 2.0 * A(i,j) + A(im1,j) ) / dx^2 ...
               + ( A(i,jp1) - 2.0 * A(i,j) + A(i,jm1) ) / dy^2;
      end
    end

  end

  return
end

