octave-maintainers
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

NDArray version of dicrete laplace operator


From: David Bateman
Subject: NDArray version of dicrete laplace operator
Date: Tue, 17 Jul 2007 19:07:41 +0200
User-agent: Thunderbird 1.5.0.7 (X11/20060921)

I have made an NDArray version of the discrete laplace operator "del2",
and have attached it here. As this is a matlab core function. I wanted
to propose this for inclusion in Octave. However, although I can get the
same behavior as the 2D version in octave-forge, the edge points of the
operator differ from the matlab version. It appears that the matlab
version of "del2" calculates the interior points first and then
extrapolates the edge points from this calculation, whereas the
octave-forge version of del2 uses a nearest value approximation.

Should I keep the octave-forge behavior or try and duplicate the Matlab
behavior of del2?

Regards
David


-- 
David Bateman                                address@hidden
Motorola Labs - Paris                        +33 1 69 35 48 04 (Ph) 
Parc Les Algorithmes, Commune de St Aubin    +33 6 72 01 06 33 (Mob) 
91193 Gif-Sur-Yvette FRANCE                  +33 1 69 35 77 01 (Fax) 

The information contained in this communication has been classified as: 

[x] General Business Information 
[ ] Motorola Internal Use Only 
[ ] Motorola Confidential Proprietary

## Copyright (C) 2000  Kai Habel
##
## This file is part of Octave.
##
## Octave is free software; you can redistribute it and/or modify it
## under the terms of the GNU General Public License as published by
## the Free Software Foundation; either version 2, or (at your option)
## any later version.
##
## Octave is distributed in the hope that it will be useful, but
## WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
## General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with Octave; see the file COPYING.  If not, write to the Free
## Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
## 02110-1301, USA.

## -*- texinfo -*-
## @deftypefn {Function File} address@hidden =} del2 (@var{m})
## @deftypefnx {Function File} address@hidden =} del2 (@var{m}, @var{h})
## @deftypefnx {Function File} address@hidden =} del2 (@var{m}, @var{dx}, 
@var{dy}, @dots{})
##
## Calculates the dicrete Laplace operator. If @var{m} is a matrix this is
## defined as
##
## @iftex
## @tex
## $d = {1 \over 4} \left( {d^2 \over dx^2} M(x,y) + {d^2 \over dy^2} M(x,y) 
\right)$
## @end tex
## @end iftex
## @ifnottex
## @example
## @group
##       1    / d^2            d^2         \
## D  = --- * | ---  M(x,y) +  ---  M(x,y) | 
##       4    \ dx^2           dx^2        /
## @end group
## @end example
## @end ifnottex
##
## The above to continued to N-dimensional arrays calculating the second
## derivative over the higher dimensions.
##
## The spacing between evaluation points may be defined by @var{h}, which is a
## scalar defining the spacing in all dimensions. Or alternative, the spacing
## in each dimension may be defined separately by @var{dx}, @var{dy}, etc. 
## Scalar spacing values give equidistant spacing, whereas vector spacing 
## values can be used to specify variable spacing. The length of the vectors
## must match the respective dimension of @var{m}. The default spacing value
## is 1.
##
## You need at least 3 data points for each dimension.
## Boundary points are calculated with y0'' and y2'' respectively.
## For interior point y1'' is taken. 
##
## @example
## y0''(i) = 1/(dy^2) *(y(i)-2*y(i+1)+y(i+2)).
## y1''(i) = 1/(dy^2) *(y(i-1)-2*y(i)+y(i+1)).
## y2''(i) = 1/(dy^2) *(y(i-2)-2*y(i-1)+y(i)).
## @end example
##
## @end deftypefn

## Author:  Kai Habel <address@hidden>

function D = del2 (M, varargin)
  
  if (nargin < 1)
    print_usage ();
  endif

  nd = ndims (M);
  sz = size (M);
  dx = cell (1, nd);
  if (nargin == 2 || nargin == 1)
    if (nargin == 1)
      h = 1;
    else
      h = varargin{1}
    endif
    for i = 1 : nd
      if (isscalar (h))
        dx{i} = h * ones (sz (i), 1);
      else
        if (length (h) == sz (i))
          dx{i} = diff (h)(:);
        else
          error ("dimensionality mismatch in %d-th spacing vector", i);
        endif
      endif
    endfor
  elseif (nargin - 1 == nd)
    ## Reverse dx{1} and dx{2} as the X-dim is the 2nd dim of the ND array
    tmp = varargin{1};
    varargin{1} = varargin{2};
    varargin{2} = tmp;

    for i = 1 : nd
      if (isscalar (varargin{i}))
        dx{i} = varargin{i} * ones (sz (i), 1);
      else
        if (length (varargin{i}) == sz (i))
          dx{i} = diff (varargin{i})(:);
        else
          error ("dimensionality mismatch in %d-th spacing vector", i);
        endif
      endif
    endfor
  else
    print_usage ();
  endif

  idx = cell (1, nd);
  for i = 1: nd
    idx{i} = ":";
  endfor

  D = zeros (sz);
  for i = 1: nd
    if (sz(i) >= 3)
      idx1 = idx2 = idx3 = idx;

      ## left and right boundary
      idx1{i} = 1;
      idx2{i} = 2;
      idx3{i} = 3;
      D(idx1{:}) = D(idx1{:}) + ... 
          (M(idx1{:}) - 2 * M(idx2{:}) + M(idx3{:})) / ...
          (dx{i}(1) * dx{i}(2));

      idx1{i} = sz(i);
      idx2{i} = sz(i) - 1;
      idx3{i} = sz(i) - 2;
      D(idx1{:}) = D(idx1{:}) + ...
          (M(idx3{:}) - 2 * M(idx2{:}) + M(idx1{:})) / ...
          (dx{i}(sz(i) - 2) * dx{i}(sz(i) - 1));

      ## interior points
      idx1{i} = 1 : sz(i) - 2;
      idx2{i} = 2 : sz(i) - 1;
      idx3{i} = 3 : sz(i);
      szi = sz;
      szi (i) = 1;
      D(idx2{:}) = D(idx2{:}) ...
          + (M(idx3{:}) - 2 * M(idx2{:}) + M(idx1{:})) ./ ...
          repmat (shiftdim (dx{i}(1 : sz(i) - 2) .* ...
                            dx{i}(2 : sz(i) - 1), 1 - i), szi);
    endif
  endfor

  D = D ./ (2 .* nd);
endfunction

reply via email to

[Prev in Thread] Current Thread [Next in Thread]