Mercurial > hg > octave-nkf
annotate scripts/general/del2.m @ 8499:e935a5158b0f
[docs] use $$ for display
| author | Brian Gough <bjg@gnu.org> |
|---|---|
| date | Tue, 13 Jan 2009 00:27:29 -0500 |
| parents | fa78cb8d8a5c |
| children | eb63fbe60fab |
| rev | line source |
|---|---|
| 7017 | 1 ## Copyright (C) 2000, 2007 Kai Habel |
| 6788 | 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 | |
| 7016 | 8 ## the Free Software Foundation; either version 3 of the License, or (at |
| 9 ## your option) any later version. | |
| 6788 | 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 | |
| 7016 | 17 ## along with Octave; see the file COPYING. If not, see |
| 18 ## <http://www.gnu.org/licenses/>. | |
| 6788 | 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 ## Calculates the discrete Laplace operator. If @var{m} is a matrix this is | |
| 26 ## defined as | |
| 27 ## | |
| 28 ## @iftex | |
| 29 ## @tex | |
| 8499 | 30 ## $$d = {1 \over 4} \left( {d^2 \over dx^2} M(x,y) + {d^2 \over dy^2} M(x,y) \right)$$ |
| 6788 | 31 ## @end tex |
| 32 ## @end iftex | |
| 33 ## @ifnottex | |
| 34 ## @example | |
| 35 ## @group | |
| 36 ## 1 / d^2 d^2 \ | |
| 37 ## D = --- * | --- M(x,y) + --- M(x,y) | | |
| 38 ## 4 \ dx^2 dy^2 / | |
| 39 ## @end group | |
| 40 ## @end example | |
| 41 ## @end ifnottex | |
| 42 ## | |
| 43 ## The above to continued to N-dimensional arrays calculating the second | |
| 44 ## derivative over the higher dimensions. | |
| 45 ## | |
| 46 ## The spacing between evaluation points may be defined by @var{h}, which is a | |
|
8347
fa78cb8d8a5c
corrections for typos
Brian Gough<bjg@network-theory.co.uk>
parents:
7669
diff
changeset
|
47 ## scalar defining the spacing in all dimensions. Or alternatively, the spacing |
| 6788 | 48 ## in each dimension may be defined separately by @var{dx}, @var{dy}, etc. |
| 49 ## Scalar spacing values give equidistant spacing, whereas vector spacing | |
| 50 ## values can be used to specify variable spacing. The length of the vectors | |
| 51 ## must match the respective dimension of @var{m}. The default spacing value | |
| 52 ## is 1. | |
| 53 ## | |
| 54 ## You need at least 3 data points for each dimension. Boundary points are | |
| 55 ## calculated as the linear extrapolation of the interior points. | |
| 56 ## | |
| 57 ## @seealso{gradient, diff} | |
| 58 ## @end deftypefn | |
| 59 | |
| 60 ## Author: Kai Habel <kai.habel@gmx.de> | |
| 61 | |
| 62 function D = del2 (M, varargin) | |
| 63 | |
| 64 if (nargin < 1) | |
| 65 print_usage (); | |
| 66 endif | |
| 67 | |
| 68 nd = ndims (M); | |
| 69 sz = size (M); | |
| 70 dx = cell (1, nd); | |
| 71 if (nargin == 2 || nargin == 1) | |
| 72 if (nargin == 1) | |
| 73 h = 1; | |
| 74 else | |
| 7669 | 75 h = varargin{1}; |
| 6788 | 76 endif |
| 77 for i = 1 : nd | |
| 78 if (isscalar (h)) | |
| 79 dx{i} = h * ones (sz (i), 1); | |
| 80 else | |
| 81 if (length (h) == sz (i)) | |
| 82 dx{i} = diff (h)(:); | |
| 83 else | |
| 84 error ("dimensionality mismatch in %d-th spacing vector", i); | |
| 85 endif | |
| 86 endif | |
| 87 endfor | |
| 88 elseif (nargin - 1 == nd) | |
| 89 ## Reverse dx{1} and dx{2} as the X-dim is the 2nd dim of the ND array | |
| 90 tmp = varargin{1}; | |
| 91 varargin{1} = varargin{2}; | |
| 92 varargin{2} = tmp; | |
| 93 | |
| 94 for i = 1 : nd | |
| 95 if (isscalar (varargin{i})) | |
| 96 dx{i} = varargin{i} * ones (sz (i), 1); | |
| 97 else | |
| 98 if (length (varargin{i}) == sz (i)) | |
| 99 dx{i} = diff (varargin{i})(:); | |
| 100 else | |
| 101 error ("dimensionality mismatch in %d-th spacing vector", i); | |
| 102 endif | |
| 103 endif | |
| 104 endfor | |
| 105 else | |
| 106 print_usage (); | |
| 107 endif | |
| 108 | |
| 109 idx = cell (1, nd); | |
| 110 for i = 1: nd | |
| 111 idx{i} = ":"; | |
| 112 endfor | |
| 113 | |
| 114 D = zeros (sz); | |
| 115 for i = 1: nd | |
| 116 if (sz(i) >= 3) | |
| 117 DD = zeros (sz); | |
| 118 idx1 = idx2 = idx3 = idx; | |
| 119 | |
| 120 ## interior points | |
| 121 idx1{i} = 1 : sz(i) - 2; | |
| 122 idx2{i} = 2 : sz(i) - 1; | |
| 123 idx3{i} = 3 : sz(i); | |
| 124 szi = sz; | |
| 125 szi (i) = 1; | |
| 126 | |
| 127 h1 = repmat (shiftdim (dx{i}(1 : sz(i) - 2), 1 - i), szi); | |
| 128 h2 = repmat (shiftdim (dx{i}(2 : sz(i) - 1), 1 - i), szi); | |
| 129 DD(idx2{:}) = ((M(idx1{:}) - M(idx2{:})) ./ h1 + ... | |
| 130 (M(idx3{:}) - M(idx2{:})) ./ h2) ./ (h1 + h2); | |
| 131 | |
| 132 ## left and right boundary | |
| 133 if (sz(i) == 3) | |
| 134 DD(idx1{:}) = DD(idx3{:}) = DD(idx2{:}); | |
| 135 else | |
| 136 idx1{i} = 1; | |
| 137 idx2{i} = 2; | |
| 138 idx3{i} = 3; | |
| 139 DD(idx1{:}) = (dx{i}(1) + dx{i}(2)) / dx{i}(2) * DD (idx2{:}) - ... | |
| 140 dx{i}(1) / dx{i}(2) * DD (idx3{:}); | |
| 141 | |
| 142 idx1{i} = sz(i); | |
| 143 idx2{i} = sz(i) - 1; | |
| 144 idx3{i} = sz(i) - 2; | |
| 145 DD(idx1{:}) = (dx{i}(sz(i) - 1) + dx{i}(sz(i) - 2)) / ... | |
| 146 dx{i}(sz(i) - 2) * DD (idx2{:}) - ... | |
| 147 dx{i}(sz(i) - 1) / dx{i}(sz(i) - 2) * DD (idx3{:}); | |
| 148 endif | |
| 149 | |
| 150 D += DD; | |
| 151 endif | |
| 152 endfor | |
| 153 | |
| 154 D = D ./ nd; | |
| 155 endfunction |