|
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 |
|
|
30 ## $d = {1 \over 4} \left( {d^2 \over dx^2} M(x,y) + {d^2 \over dy^2} M(x,y) \right)$ |
|
|
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 |
|
|
47 ## scalar defining the spacing in all dimensions. Or alternative, the spacing |
|
|
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 |
