octave_packages/m/general/del2.m

   1 ## Copyright (C) 2000-2012 Kai Habel
   2 ## Copyright (C) 2007  David Bateman
   3 ##
   4 ## This file is part of Octave.
   5 ##
   6 ## Octave is free software; you can redistribute it and/or modify it
   7 ## under the terms of the GNU General Public License as published by
   8 ## the Free Software Foundation; either version 3 of the License, or (at
   9 ## your option) any later version.
  10 ##
  11 ## Octave is distributed in the hope that it will be useful, but
  12 ## WITHOUT ANY WARRANTY; without even the implied warranty of
  13 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  14 ## General Public License for more details.
  15 ##
  16 ## You should have received a copy of the GNU General Public License
  17 ## along with Octave; see the file COPYING.  If not, see
  18 ## <http://www.gnu.org/licenses/>.
  19
  20 ## -*- texinfo -*-
  21 ## @deftypefn  {Function File} {@var{d} =} del2 (@var{M})
  22 ## @deftypefnx {Function File} {@var{d} =} del2 (@var{M}, @var{h})
  23 ## @deftypefnx {Function File} {@var{d} =} del2 (@var{M}, @var{dx}, @var{dy}, @dots{})
  24 ##
  25 ## Calculate the discrete Laplace
  26 ## @tex
  27 ## operator $( \nabla^2 )$.
  28 ## @end tex
  29 ## @ifnottex
  30 ## operator.
  31 ## @end ifnottex
  32 ## For a 2-dimensional matrix @var{M} this is defined as
  33 ## @tex
  34 ## $$d = {1 \over 4} \left( {d^2 \over dx^2} M(x,y) + {d^2 \over dy^2} M(x,y) \right)$$
  35 ## @end tex
  36 ## @ifnottex
  37 ##
  38 ## @example
  39 ## @group
  40 ##       1    / d^2            d^2         \
  41 ## D  = --- * | ---  M(x,y) +  ---  M(x,y) |
  42 ##       4    \ dx^2           dy^2        /
  43 ## @end group
  44 ## @end example
  45 ##
  46 ## @end ifnottex
  47 ## For N-dimensional arrays the sum in parentheses is expanded to include second
  48 ## derivatives over the additional higher dimensions.
  49 ##
  50 ## The spacing between evaluation points may be defined by @var{h}, which is a
  51 ## scalar defining the equidistant spacing in all dimensions.  Alternatively,
  52 ## the spacing in each dimension may be defined separately by @var{dx},
  53 ## @var{dy}, etc.  A scalar spacing argument defines equidistant spacing,
  54 ## whereas a vector argument can be used to specify variable spacing.  The
  55 ## length of the spacing vectors must match the respective dimension of
  56 ## @var{M}.  The default spacing value is 1.
  57 ##
  58 ## At least 3 data points are needed for each dimension.  Boundary points are
  59 ## calculated from the linear extrapolation of interior points.
  60 ##
  61 ## @seealso{gradient, diff}
  62 ## @end deftypefn
  63
  64 ## Author:  Kai Habel <kai.habel@gmx.de>
  65
  66 function D = del2 (M, varargin)
  67
  68   if (nargin < 1)
  69     print_usage ();
  70   endif
  71
  72   nd = ndims (M);
  73   sz = size (M);
  74   dx = cell (1, nd);
  75   if (nargin == 2 || nargin == 1)
  76     if (nargin == 1)
  77       h = 1;
  78     else
  79       h = varargin{1};
  80     endif
  81     for i = 1 : nd
  82       if (isscalar (h))
  83         dx{i} = h * ones (sz (i), 1);
  84       else
  85         if (length (h) == sz (i))
  86           dx{i} = diff (h)(:);
  87         else
  88           error ("del2: dimensionality mismatch in %d-th spacing vector", i);
  89         endif
  90       endif
  91     endfor
  92   elseif (nargin - 1 == nd)
  93     ## Reverse dx{1} and dx{2} as the X-dim is the 2nd dim of the ND array
  94     tmp = varargin{1};
  95     varargin{1} = varargin{2};
  96     varargin{2} = tmp;
  97
  98     for i = 1 : nd
  99       if (isscalar (varargin{i}))
 100         dx{i} = varargin{i} * ones (sz (i), 1);
 101       else
 102         if (length (varargin{i}) == sz (i))
 103           dx{i} = diff (varargin{i})(:);
 104         else
 105           error ("del2: dimensionality mismatch in %d-th spacing vector", i);
 106         endif
 107       endif
 108     endfor
 109   else
 110     print_usage ();
 111   endif
 112
 113   idx = cell (1, nd);
 114   for i = 1: nd
 115     idx{i} = ":";
 116   endfor
 117
 118   D = zeros (sz);
 119   for i = 1: nd
 120     if (sz(i) >= 3)
 121       DD = zeros (sz);
 122       idx1 = idx2 = idx3 = idx;
 123
 124       ## interior points
 125       idx1{i} = 1 : sz(i) - 2;
 126       idx2{i} = 2 : sz(i) - 1;
 127       idx3{i} = 3 : sz(i);
 128       szi = sz;
 129       szi (i) = 1;
 130
 131       h1 = repmat (shiftdim (dx{i}(1 : sz(i) - 2), 1 - i), szi);
 132       h2 = repmat (shiftdim (dx{i}(2 : sz(i) - 1), 1 - i), szi);
 133       DD(idx2{:}) = ((M(idx1{:}) - M(idx2{:})) ./ h1 + ...
 134                      (M(idx3{:}) - M(idx2{:})) ./ h2) ./ (h1 + h2);
 135
 136       ## left and right boundary
 137       if (sz(i) == 3)
 138         DD(idx1{:}) = DD(idx3{:}) = DD(idx2{:});
 139       else
 140         idx1{i} = 1;
 141         idx2{i} = 2;
 142         idx3{i} = 3;
 143         DD(idx1{:}) = (dx{i}(1) + dx{i}(2)) / dx{i}(2) * DD (idx2{:}) - ...
 144             dx{i}(1) / dx{i}(2) * DD (idx3{:});
 145
 146         idx1{i} = sz(i);
 147         idx2{i} = sz(i) - 1;
 148         idx3{i} = sz(i) - 2;
 149         DD(idx1{:}) =  (dx{i}(sz(i) - 1) + dx{i}(sz(i) - 2)) / ...
 150             dx{i}(sz(i) - 2) * DD (idx2{:}) - ...
 151             dx{i}(sz(i) - 1) / dx{i}(sz(i) - 2) * DD (idx3{:});
 152       endif
 153
 154       D += DD;
 155     endif
 156   endfor
 157
 158   D = D ./ nd;
 159 endfunction