function cof = r8vec_dif ( n, h )
%*****************************************************************************80
%
%% r8vec_dif() computes coefficients for estimating the N-th derivative.
%
% Discussion:
%
% The routine computes the N+1 coefficients for a centered finite difference
% estimate of the N-th derivative of a function.
%
% The estimate has the form
%
% FDIF(N,X) = Sum (I = 0 to N) COF(I) * F ( X(I) )
%
% To understand the computation of the coefficients, it is enough
% to realize that the first difference approximation is
%
% FDIF(1,X) = F(X+DX) - F(X-DX) ) / (2*DX)
%
% and that the second difference approximation can be regarded as
% the first difference approximation repeated:
%
% FDIF(2,X) = FDIF(1,X+DX) - FDIF(1,X-DX) / (2*DX)
% = F(X+2*DX) - 2 F(X) + F(X-2*DX) / (4*DX)
%
% and so on for higher order differences.
%
% Thus, the next thing to consider is the integer coefficients of
% the sampled values of F, which are clearly the Pascal coefficients,
% but with an alternating negative sign. In particular, if we
% consider row I of Pascal's triangle to have entries j = 0 through I,
% then P(I,J) = P(I-1,J-1) - P(I-1,J), where P(*,-1) is taken to be 0,
% and P(0,0) = 1.
%
% 1
% -1 1
% 1 -2 1
% -1 3 -3 1
% 1 -4 6 -4 1
% -1 5 -10 10 -5 1
% 1 -6 15 -20 15 -6 1
%
% Next, note that the denominator of the approximation for the
% N-th derivative will be (2*DX)**N.
%
% And finally, consider the location of the N+1 sampling
% points for F:
%
% X-N*DX, X-(N-2)*DX, X-(N-4)*DX, ..., X+(N-4)*DX, X+(N-2*DX), X+N*DX.
%
% Thus, a formula for evaluating FDIF(N,X) is
%
% fdif = 0.0
% do i = 0, n
% xi = x + (2*i-n) * h
% fdif = fdif + cof(i) * f(xi)
% end do
%
% Licensing:
%
% This code is distributed under the GNU LGPL license.
%
% Modified:
%
% 26 April 2005
%
% Author:
%
% John Burkardt
%
% Input:
%
% integer N, the order of the derivative to be approximated.
% N must be 0 or greater.
%
% real H, the half spacing between points.
% H must be positive.
%
% Output:
%
% real COF(1:N+1), the coefficients needed to approximate
% the N-th derivative of a function F.
%
if ( n < 0 )
fprintf ( 1, '\n' );
fprintf ( 1, 'R8VEC_DIF - Fatal error!\n' );
fprintf ( 1, ' Derivative order N = %d\n', n );
fprintf ( 1, ' but N must be at least 0.\n' );
error ( 'R8VEC_DIF - Fatal error!' );
end
if ( h <= 0.0 )
fprintf ( 1, '\n' );
fprintf ( 1, 'R8VEC_DIF - Fatal error!\n' );
fprintf ( 1, ' The half sampling spacing is H = %f\n', h );
fprintf ( 1, ' but H must be positive.\n' );
error ( 'R8VEC_DIF - Fatal error!' );
end
for i = 0 : n
cof(i+1) = 1.0;
for j = i-1 : -1 : 1
cof(j+1) = -cof(j+1) + cof(j);
end
if ( 0 < i )
cof(1) = -cof(1);
end
end
cof(1:n+1) = cof(1:n+1) / ( 2.0 * h ).^n;
return
end