Mercurial > octave
diff scripts/general/del2.m @ 6788:c81a0f3f5a82
[project @ 2007-07-23 22:05:29 by dbateman]
author | dbateman |
---|---|
date | Mon, 23 Jul 2007 22:05:30 +0000 |
parents | |
children | 93c65f2a5668 |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/general/del2.m Mon Jul 23 22:05:30 2007 +0000 @@ -0,0 +1,156 @@ +## Copyright (C) 2000 Kai Habel +## Copyright (C) 2007 David Bateman +## +## 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} {@var{d} =} del2 (@var{m}) +## @deftypefnx {Function File} {@var{d} =} del2 (@var{m}, @var{h}) +## @deftypefnx {Function File} {@var{d} =} del2 (@var{m}, @var{dx}, @var{dy}, @dots{}) +## +## Calculates the discrete 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 dy^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 as the linear extrapolation of the interior points. +## +## @seealso{gradient, diff} +## @end deftypefn + +## Author: Kai Habel <kai.habel@gmx.de> + +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) + DD = zeros (sz); + idx1 = idx2 = idx3 = idx; + + ## interior points + idx1{i} = 1 : sz(i) - 2; + idx2{i} = 2 : sz(i) - 1; + idx3{i} = 3 : sz(i); + szi = sz; + szi (i) = 1; + + h1 = repmat (shiftdim (dx{i}(1 : sz(i) - 2), 1 - i), szi); + h2 = repmat (shiftdim (dx{i}(2 : sz(i) - 1), 1 - i), szi); + DD(idx2{:}) = ((M(idx1{:}) - M(idx2{:})) ./ h1 + ... + (M(idx3{:}) - M(idx2{:})) ./ h2) ./ (h1 + h2); + + ## left and right boundary + if (sz(i) == 3) + DD(idx1{:}) = DD(idx3{:}) = DD(idx2{:}); + else + idx1{i} = 1; + idx2{i} = 2; + idx3{i} = 3; + DD(idx1{:}) = (dx{i}(1) + dx{i}(2)) / dx{i}(2) * DD (idx2{:}) - ... + dx{i}(1) / dx{i}(2) * DD (idx3{:}); + + idx1{i} = sz(i); + idx2{i} = sz(i) - 1; + idx3{i} = sz(i) - 2; + DD(idx1{:}) = (dx{i}(sz(i) - 1) + dx{i}(sz(i) - 2)) / ... + dx{i}(sz(i) - 2) * DD (idx2{:}) - ... + dx{i}(sz(i) - 1) / dx{i}(sz(i) - 2) * DD (idx3{:}); + endif + + D += DD; + endif + endfor + + D = D ./ nd; +endfunction