Mercurial > octave
view scripts/general/del2.m @ 31058:12f8fb75fc30
quadgk.m: Overhaul function to add new "ArrayValued" option (bug #62468)
* quadgk.m: Add "ArrayValued" option to docstring. Re-phrase bits of existing
documentation for clarity. Rename input validation variable "str" to "prop"
for clarity. Add decode for "arrayvalued" property. If limits of integration
are high to low then reverse them in function rather than recursively calling
quadgk(). Add FIXME note about missing input validation for several properties.
Delete variable 'h0' which was never used. New code branch if "ArrayValued"
is true which uses same code strategy as for normal quadgk, but subfunction
__quadgk_eval_array__ for evaluating the function in a vectorized manner.
Add BIST tests for "ArrayValued" input.
* quadgk.m (__quadgk_eval__): Eliminate "too_close" output and code to
calculate "too_close" which was never used.
* quadgk.m (__quadgk_eval_array__): New subfunction.
author | Michael Leitner <michael.leitner@frm2.tum.de>, Rik <rik@octave.org> |
---|---|
date | Thu, 02 Jun 2022 08:45:10 -0700 |
parents | 796f54d4ddbf |
children | a40c0b7aa376 |
line wrap: on
line source
######################################################################## ## ## Copyright (C) 2000-2022 The Octave Project Developers ## ## See the file COPYRIGHT.md in the top-level directory of this ## distribution or <https://octave.org/copyright/>. ## ## 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 3 of the License, 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, see ## <https://www.gnu.org/licenses/>. ## ######################################################################## ## -*- texinfo -*- ## @deftypefn {} {@var{L} =} del2 (@var{M}) ## @deftypefnx {} {@var{L} =} del2 (@var{M}, @var{h}) ## @deftypefnx {} {@var{L} =} del2 (@var{M}, @var{dx}, @var{dy}, @dots{}) ## ## Calculate the discrete Laplace ## @tex ## operator $( \nabla^2 )$. ## @end tex ## @ifnottex ## operator. ## @end ifnottex ## ## For a 2-dimensional matrix @var{M} this is defined as ## @tex ## $$L = {1 \over 4} \left( {d^2 \over dx^2} M(x,y) + {d^2 \over dy^2} M(x,y) \right)$$ ## @end tex ## @ifnottex ## ## @example ## @group ## 1 / d^2 d^2 \ ## L = --- * | --- M(x,y) + --- M(x,y) | ## 4 \ dx^2 dy^2 / ## @end group ## @end example ## ## @end ifnottex ## For N-dimensional arrays the sum in parentheses is expanded to include ## second derivatives over the additional higher dimensions. ## ## The spacing between evaluation points may be defined by @var{h}, which is a ## scalar defining the equidistant spacing in all dimensions. Alternatively, ## the spacing in each dimension may be defined separately by @var{dx}, ## @var{dy}, etc. A scalar spacing argument defines equidistant spacing, ## whereas a vector argument can be used to specify variable spacing. The ## length of the spacing vectors must match the respective dimension of ## @var{M}. The default spacing value is 1. ## ## Dimensions with fewer than 3 data points are skipped. Boundary points are ## calculated from the linear extrapolation of interior points. ## ## Example: Second derivative of 2*x^3 ## ## @example ## @group ## f = @@(x) 2*x.^3; ## dd = @@(x) 12*x; ## x = 1:6; ## L = 4*del2 (f(x)); ## assert (L, dd (x)); ## @end group ## @end example ## ## @seealso{gradient, diff} ## @end deftypefn function L = del2 (M, varargin) if (nargin < 1) print_usage (); endif nd = ndims (M); sz = size (M); dx = cell (1, nd); if (nargin == 1) for i = 1 : nd dx(i) = ones (sz(i), 1); endfor elseif (nargin == 2 && isscalar (varargin{1})) h = varargin{1}; for i = 1 : nd dx(i) = h * ones (sz(i), 1); endfor elseif (numel (varargin) <= nd) ndx = numel (varargin); varargin(ndx+1:nd) = 1; # Fill missing dims with 1. ## Reverse dx{1} and dx{2} as the X-dim is the 2nd dim of a meshgrid array varargin([1, 2]) = varargin([2, 1]); for i = 1 : nd arg = varargin{i}; if (isscalar (arg)) dx(i) = arg * ones (sz(i), 1); elseif (isvector (arg)) if (length (arg) != sz(i)) error ("del2: number of elements in spacing vector %d does not match dimension %d of M", i, i); endif dx(i) = diff (varargin{i})(:); else error ("del2: spacing element %d must be a scalar or vector", i); endif endfor else print_usage (); endif idx = cell (1, nd); idx(:) = ":"; L = 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 L += DD; endif endfor L ./= nd; endfunction ## 3x3 constant test %!test %! a = ones (3,3); %! b = del2 (a); %! assert (b(:,1), [0.00;0.00;0.00]); %! assert (b(:,2), [0.00;0.00;0.00]); %! assert (b(:,3), [0.00;0.00;0.00]); ## 3x3 planar test %!test %! a = [1,2,3;2,3,4;3,4,5]; %! b = del2 (a); %! assert (b(:,1), [0.00;0.00;0.00]); %! assert (b(:,2), [0.00;0.00;0.00]); %! assert (b(:,3), [0.00;0.00;0.00]); ## 3x3 corner test %!test %! a = zeros (3,3); %! a(1,1) = 1.0; %! b = 2*del2 (a); %! assert (b(:,1), [1.00;0.50;0.50]); %! assert (b(:,2), [0.50;0.00;0.00]); %! assert (b(:,3), [0.50;0.00;0.00]); %! assert (b, flipud (2*del2 (flipud (a)))); %! assert (b, fliplr (2*del2 (fliplr (a)))); %! assert (b, flipud (fliplr (2*del2 (fliplr (flipud (a)))))); ## 3x3 boundary test %!test %! a = zeros (3,3); %! a(2,1)=1.0; %! b = 2*del2 (a); %! assert (b(:,1), [-1.00;-0.50;-1.00]); %! assert (b(:,2), [0.00;0.50;0.00]); %! assert (b(:,3), [0.00;0.50;0.00]); %! assert (b, flipud (2*del2 (flipud (a)))); %! assert (b, fliplr (2*del2 (fliplr (a)))); %! assert (b, flipud (fliplr (2*del2 (fliplr (flipud (a)))))); ## 3x3 center test %!test %! a = zeros (3,3); %! a(2,2) = 1.0; %! b = del2 (a); %! assert (b(:,1), [0.00;-0.50;0.00]); %! assert (b(:,2), [-0.50;-1.00;-0.50]); %! assert (b(:,3), [0.00;-0.50;0.00]); ## 4x4 constant test %!test %! a = ones (4,4); %! b = del2 (a); %! assert (b(:,1), [0.00;0.00;0.00;0.00]); %! assert (b(:,2), [0.00;0.00;0.00;0.00]); %! assert (b(:,3), [0.00;0.00;0.00;0.00]); %! assert (b(:,4), [0.00;0.00;0.00;0.00]); ## 4x4 planar test %!test %! a = [1,2,3,4;2,3,4,5;3,4,5,6;4,5,6,7]; %! b = del2 (a); %! assert (b(:,1), [0.00;0.00;0.00;0.00]); %! assert (b(:,2), [0.00;0.00;0.00;0.00]); %! assert (b(:,3), [0.00;0.00;0.00;0.00]); %! assert (b(:,4), [0.00;0.00;0.00;0.00]); ## 4x4 corner test %!test %! a = zeros (4,4); %! a(1,1) = 1.0; %! b = 2*del2 (a); %! assert (b(:,1), [2.00;0.50;0.00;-0.50]); %! assert (b(:,2), [0.50;0.00;0.00;0.00]); %! assert (b(:,3), [0.00;0.00;0.00;0.00]); %! assert (b(:,4), [-0.50;0.00;0.00;0.00]); %! assert (b, flipud (2*del2 (flipud (a)))); %! assert (b, fliplr (2*del2 (fliplr (a)))); %! assert (b, flipud (fliplr (2*del2 (fliplr (flipud (a)))))); ## 9x9 center test %!test %! a = zeros (9,9); %! a(5,5) = 1.0; %! b = 2*del2 (a); %! assert (b(:,1), [0.00;0.00;0.00;0.00;0.00;0.00;0.00;0.00;0.00]); %! assert (b(:,2), [0.00;0.00;0.00;0.00;0.00;0.00;0.00;0.00;0.00]); %! assert (b(:,3), [0.00;0.00;0.00;0.00;0.00;0.00;0.00;0.00;0.00]); %! assert (b(:,4), [0.00;0.00;0.00;0.00;0.50;0.00;0.00;0.00;0.00]); %! assert (b(:,5), [0.00;0.00;0.00;0.50;-2.00;0.50;0.00;0.00;0.00]); %! assert (b(:,6), b(:,4)); %! assert (b(:,7), b(:,3)); %! assert (b(:,8), b(:,2)); %! assert (b(:,9), b(:,1)); ## 9x9 boundary test %!test %! a = zeros (9,9); %! a(1,5) = 1.0; %! b = 2*del2 (a); %! assert (b(1,:), [0.00,0.00,0.00,0.50,0.00,0.50,0.00,0.00,0.00]); %! assert (b(2,:), [0.00,0.00,0.00,0.00,0.50,0.00,0.00,0.00,0.00]); %! assert (b(3:9,:), zeros (7,9)); %! a(1,5) = 0.0; %! a(5,1) = 1.0; %! b = 2*del2 (a); %! assert (b(:,1), [0.00;0.00;0.00;0.50;0.00;0.50;0.00;0.00;0.00]); %! assert (b(:,2), [0.00;0.00;0.00;0.00;0.50;0.00;0.00;0.00;0.00]); %! assert (b(:,3:9), zeros (9,7)); ## 9x9 dh center test %!test %! a = zeros (9,9); %! a(5,5) = 1.0; %! b = 8*del2 (a,2); %! assert (b(:,1:3), zeros (9,3)); %! assert (b(:,4), [0.00;0.00;0.00;0.00;0.50;0.00;0.00;0.00;0.00]); %! assert (b(:,5), [0.00;0.00;0.00;0.50;-2.00;0.50;0.00;0.00;0.00]); %! assert (b(:,6), b(:,4)); %! assert (b(:,7:9), zeros (9,3)); ## 9x9 dx test %!test %! a = zeros (9,9); %! a(5,5) = 1.0; %! b = 4*del2 (a,2,1); %! assert (b(1:3,:), zeros (3,9)); %! assert (b(4,:), [0.00;0.00;0.00;0.00;1.00;0.00;0.00;0.00;0.00]'); %! assert (b(5,:), [0.00;0.00;0.00;0.25;-2.5;0.25;0.00;0.00;0.00]'); %! assert (b(6,:), b(4,:)); %! assert (b(7:9,:), zeros (3,9)); ## 9x9 dy test %!test %! a = zeros (9,9); %! a(5,5) = 1.0; %! b = 4*del2 (a,1,2); %! assert (b(:,1:3), zeros (9,3)); %! assert (b(:,4), [0.00;0.00;0.00;0.00;1.00;0.00;0.00;0.00;0.00]); %! assert (b(:,5), [0.00;0.00;0.00;0.25;-2.5;0.25;0.00;0.00;0.00]); %! assert (b(:,6), b(:,4)); %! assert (b(:,7:9), zeros (9,3)); ## 3D test %!test %! a = zeros (9,9,9); %! a(5,5,5) = 1.0; %! b = 8*3*del2 (a,2); %! assert (b(:,:,1:3), zeros (9,9,3)); %! assert (b(:,1:3,:), zeros (9,3,9)); %! assert (b(1:3,:,:), zeros (3,9,9)); %! assert (b(4:5,4,4), [0.0,0.0]'); %! assert (b(5,5,4), 1.00); %! assert (b(4,4,5), 0.00); %! assert (b(5,4,5), 1.00); %! assert (b(5,5,5),-6.00); %! assert (b, flip (b,1)); %! assert (b, flip (b,2)); %! assert (b, flip (b,3)); %!test <*51728> %! x = linspace (-2*pi, 2*pi); %! U = cos (x); %! L = 4*del2 (U, x); ## Test input validation %!error <Invalid call> del2 () %!error <Invalid call> del2 (1, 1, 2, 3) %!error <in spacing vector 1> del2 (1, 2, [1 1]) %!error <in spacing vector 2> del2 (1, [1 1], 2) %!error <must be a scalar or vector> del2 (1, ones (2,2), 2)