## Copyright (C) 2006, 2007 G. D. McBain ## ## This program 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 of the License, or (at ## your option) any later version. ## ## This program 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 this program; if not, write to the Free Software ## Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301USA ## -*- texinfo -*- ## @deftypefn {Function File} {} finidiff (@var{x}, @var{s}, @var{k}, @var{b}, @var{z}) ## Compute the finite difference coefficients for the derivative of ## order @var{k} at the point @var{s} using the abscissae @var{x}; exact ## for polynomials of degree less than @address@hidden (where @code{n = ## length (@var{x})}) with a zero of order @var{b} at @var{z}, or, ## equivalently, exact for (x - z)^p for p = b:(b+n). ## ## The @var{z} is irrelevant if @var{b} is zero, and the last two ## arguments default to zero. Setting @var{b} greater than zero is ## useful for homogeneous Dirichlet boundary conditions. ## ## An interpolation stencil is returned if @address@hidden ## @end deftypefn ## ### BACKGROUND THEORY ## ## The idea is that we want the discrete operator (a row of ## coefficients) to act on a discrete function (column of ordinates) to ## give the value of the derivative. We determine the n = length (x) ## coefficients by requiring the row to work exactly for n test ## functions, the monomials (x-z).^(p-i) for i=1:n where p = length (x) ## + b. ## ## ### ## The shapes ## ## We want 1*n weights to act on n*1 ordinates to give 1 answer (or ## perhaps m*n to act on n*1 to give m*1?). ## ## The condition determining d is d * V = g, which is (1*n) * (n*n) = ## (1*n), which we transpose to get V' * d' = g', which is (n*n) * (n*1) ## = (n*1), so d' = inverse (V') * g', or d = (inverse (V') * g')', ## which is exactly what is computed by Octave's `right division' ## operator; i.e. d = g / V . ## ## Now the columns of V represent the test functions (by their ## ordinates) and the elements of the row g are the exact answers. ## ## ### ## The elements ## ## Now, the q-th derivative of x.^r is ## ## r*(r-1)*...*(r-(q-1))*x.^(r-q) ## ## = prod (r-(1:q-1)) * x .^ (r-q) ## ## = gamma (r+1) / gamma (r-q+1) * x .^ (r-q) , ## ## unless of course q > r, in which case it's zero; thus that of ## (x-z).^(p-i) at x = s is ## ## gamma (p-i+1) ./ gamma (p-i-q+1) * (x-z) .^ (p-i-q) ## ## for i = 1:(p-q) and zero for i = (p-q+1):n . ## ### REFERENCES ## ## P. J. Davis & I. Polonsky 1965 Numerical interpolation, ## differentiation, and integration. In M. Abramowitz & I. Stegun, ## Handbook of Mathematical Functions, Dover, New York, ch. 25, ## pp.~875--924. ## Author: G. D. McBain ## Created: 2006-11-27 ## Keywords: finite difference, uneven grid, boundary conditions function d = finidiff (x, s, q, b, z) switch nargin case num2cell (0:2) print_usage (); case 3 b=0; z=s; endswitch n = length (x); p = n + b; i = 1:(p-q); d = postpad (gamma (p-i+1) ./ gamma (p-i-q+1) .* (s-z) .^ (p-i-q), n) (:) ' \ / vander (postpad (x-z, p)) (1:n,1:n); endfunction %!test# central three-point 2nd derivative on evenly spaced grid %! assert (finidiff (-1:1, 0, 2), [1, -2, 1]); %!test# central five-point 4th derivative on evenly spaced grid %! assert (finidiff (-2:2, 0, 4), [1, -4, 6, -4, 1], 1e-12); %!test# differentiate x.^(0:8) 8 times on 9 random points; %! ## the answers should all be zero except the last, 8! = 40320 %! n = 8; %! x = rand (1, n+1); # note: not even necessarily sorted %! d = finidiff (x, 0.5, n); %! for i = 0:n %! kthderiv(i+1) = d*(x.^i).'; %! endfor %! assert (kthderiv', prepad (gamma (n+1), n+1), 1e-4); %!demo# derivatives on even grid, no boundary conditions %! # reproduce table 25.2 of Davis & Polonsky (1965, p. 914) %! for k = 1:5 %! printf ("k=%d\n", k); %! for m = (max (2, k)):5 %! printf ("\t\t m=%d\n\n", m); %! for j = 0:m %! printf ("%1d ", j); %! d = gamma (m+1) / gamma (k+1) * finidiff (0:m, j, k); %! printf ("% 5.0f ", d'); printf ("\n"); %! endfor %! printf ("\n"); %! endfor %! printf ("\n"); %! endfor %!demo# ditto, but with f (0) = 0 %! # ditto minus coefficient of f (0) %! for k = 1:5 %! printf ("k=%d\n", k); %! for m = (max (2, k)):5 %! printf ("\t\t m=%d\n\n", m); %! for j = 0:m %! printf ("%1d ", j); %! d = gamma (m+1) / gamma (k+1) * finidiff (1:m, j, k, 1, 0); %! printf ("% 5.0f ", d'); printf ("\n"); %! endfor %! printf ("\n"); %! endfor %! printf ("\n"); %! endfor %!demo# three-point Lagrangian interpolation coefficients %! # reproduce p. 900 of Davis & Polonsky (1965) %! x = -1:1; %! for p = 0:0.01:1 %! printf ("%4.2f % 7.5f % 7.5f % 7.5f\n", p, finidiff (x, p, 0)); %! endfor %!demo# four-point Lagrangian interpolation coefficients %! # reproduce p. 901 of Davis & Polonsky (1965) %! x = -1:2; %! for p = 0:0.01:0.5 %! printf ("%4.2f % 7.7f % 7.7f % 7.7f % 7.7f\n", p, finidiff (x, p, 0)); %! endfor %!demo# eight-point Lagrangian interpolation coefficients %! # reproduce data of Davis & Polonsky (1965, table 25.1, p. 913) %! x = -3:4; %! for p = 0:0.1:4 %! printf ("%3.1f", p); %! printf (" % 7.5f", finidiff (x, p, 0)); printf ("\n"); %! endfor %!demo# Lagrangian integration coefficients---CURRENTLY BROKEN %! # reproduce data of Davis & Polonsky (1965, table 25.3, p. 915) %! printf ("%d ", finidiff (-4:4, 1, -1) * 3628800); printf ("\n"); %! %! ## It is true that the formula works correctly giving the definite %! ## integral from z.