|
Packit |
67cb25 |
/* deriv/deriv.c
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* Copyright (C) 2004, 2007 Brian Gough
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* This program is free software; you can redistribute it and/or modify
|
|
Packit |
67cb25 |
* it under the terms of the GNU General Public License as published by
|
|
Packit |
67cb25 |
* the Free Software Foundation; either version 3 of the License, or (at
|
|
Packit |
67cb25 |
* your option) any later version.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* This program is distributed in the hope that it will be useful, but
|
|
Packit |
67cb25 |
* WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
Packit |
67cb25 |
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
|
|
Packit |
67cb25 |
* General Public License for more details.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* You should have received a copy of the GNU General Public License
|
|
Packit |
67cb25 |
* along with this program; if not, write to the Free Software
|
|
Packit |
67cb25 |
* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
#include <config.h>
|
|
Packit |
67cb25 |
#include <stdlib.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_math.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_errno.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_deriv.h>
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static void
|
|
Packit |
67cb25 |
central_deriv (const gsl_function * f, double x, double h,
|
|
Packit |
67cb25 |
double *result, double *abserr_round, double *abserr_trunc)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* Compute the derivative using the 5-point rule (x-h, x-h/2, x,
|
|
Packit |
67cb25 |
x+h/2, x+h). Note that the central point is not used.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Compute the error using the difference between the 5-point and
|
|
Packit |
67cb25 |
the 3-point rule (x-h,x,x+h). Again the central point is not
|
|
Packit |
67cb25 |
used. */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double fm1 = GSL_FN_EVAL (f, x - h);
|
|
Packit |
67cb25 |
double fp1 = GSL_FN_EVAL (f, x + h);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double fmh = GSL_FN_EVAL (f, x - h / 2);
|
|
Packit |
67cb25 |
double fph = GSL_FN_EVAL (f, x + h / 2);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double r3 = 0.5 * (fp1 - fm1);
|
|
Packit |
67cb25 |
double r5 = (4.0 / 3.0) * (fph - fmh) - (1.0 / 3.0) * r3;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double e3 = (fabs (fp1) + fabs (fm1)) * GSL_DBL_EPSILON;
|
|
Packit |
67cb25 |
double e5 = 2.0 * (fabs (fph) + fabs (fmh)) * GSL_DBL_EPSILON + e3;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* The next term is due to finite precision in x+h = O (eps * x) */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double dy = GSL_MAX (fabs (r3 / h), fabs (r5 / h)) *(fabs (x) / h) * GSL_DBL_EPSILON;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* The truncation error in the r5 approximation itself is O(h^4).
|
|
Packit |
67cb25 |
However, for safety, we estimate the error from r5-r3, which is
|
|
Packit |
67cb25 |
O(h^2). By scaling h we will minimise this estimated error, not
|
|
Packit |
67cb25 |
the actual truncation error in r5. */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
*result = r5 / h;
|
|
Packit |
67cb25 |
*abserr_trunc = fabs ((r5 - r3) / h); /* Estimated truncation error O(h^2) */
|
|
Packit |
67cb25 |
*abserr_round = fabs (e5 / h) + dy; /* Rounding error (cancellations) */
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_deriv_central (const gsl_function * f, double x, double h,
|
|
Packit |
67cb25 |
double *result, double *abserr)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double r_0, round, trunc, error;
|
|
Packit |
67cb25 |
central_deriv (f, x, h, &r_0, &round, &trunc);
|
|
Packit |
67cb25 |
error = round + trunc;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (round < trunc && (round > 0 && trunc > 0))
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double r_opt, round_opt, trunc_opt, error_opt;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Compute an optimised stepsize to minimize the total error,
|
|
Packit |
67cb25 |
using the scaling of the truncation error (O(h^2)) and
|
|
Packit |
67cb25 |
rounding error (O(1/h)). */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double h_opt = h * pow (round / (2.0 * trunc), 1.0 / 3.0);
|
|
Packit |
67cb25 |
central_deriv (f, x, h_opt, &r_opt, &round_opt, &trunc_opt);
|
|
Packit |
67cb25 |
error_opt = round_opt + trunc_opt;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Check that the new error is smaller, and that the new derivative
|
|
Packit |
67cb25 |
is consistent with the error bounds of the original estimate. */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (error_opt < error && fabs (r_opt - r_0) < 4.0 * error)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
r_0 = r_opt;
|
|
Packit |
67cb25 |
error = error_opt;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
*result = r_0;
|
|
Packit |
67cb25 |
*abserr = error;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static void
|
|
Packit |
67cb25 |
forward_deriv (const gsl_function * f, double x, double h,
|
|
Packit |
67cb25 |
double *result, double *abserr_round, double *abserr_trunc)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* Compute the derivative using the 4-point rule (x+h/4, x+h/2,
|
|
Packit |
67cb25 |
x+3h/4, x+h).
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Compute the error using the difference between the 4-point and
|
|
Packit |
67cb25 |
the 2-point rule (x+h/2,x+h). */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double f1 = GSL_FN_EVAL (f, x + h / 4.0);
|
|
Packit |
67cb25 |
double f2 = GSL_FN_EVAL (f, x + h / 2.0);
|
|
Packit |
67cb25 |
double f3 = GSL_FN_EVAL (f, x + (3.0 / 4.0) * h);
|
|
Packit |
67cb25 |
double f4 = GSL_FN_EVAL (f, x + h);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double r2 = 2.0*(f4 - f2);
|
|
Packit |
67cb25 |
double r4 = (22.0 / 3.0) * (f4 - f3) - (62.0 / 3.0) * (f3 - f2) +
|
|
Packit |
67cb25 |
(52.0 / 3.0) * (f2 - f1);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Estimate the rounding error for r4 */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double e4 = 2 * 20.67 * (fabs (f4) + fabs (f3) + fabs (f2) + fabs (f1)) * GSL_DBL_EPSILON;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* The next term is due to finite precision in x+h = O (eps * x) */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double dy = GSL_MAX (fabs (r2 / h), fabs (r4 / h)) * fabs (x / h) * GSL_DBL_EPSILON;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* The truncation error in the r4 approximation itself is O(h^3).
|
|
Packit |
67cb25 |
However, for safety, we estimate the error from r4-r2, which is
|
|
Packit |
67cb25 |
O(h). By scaling h we will minimise this estimated error, not
|
|
Packit |
67cb25 |
the actual truncation error in r4. */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
*result = r4 / h;
|
|
Packit |
67cb25 |
*abserr_trunc = fabs ((r4 - r2) / h); /* Estimated truncation error O(h) */
|
|
Packit |
67cb25 |
*abserr_round = fabs (e4 / h) + dy;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_deriv_forward (const gsl_function * f, double x, double h,
|
|
Packit |
67cb25 |
double *result, double *abserr)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double r_0, round, trunc, error;
|
|
Packit |
67cb25 |
forward_deriv (f, x, h, &r_0, &round, &trunc);
|
|
Packit |
67cb25 |
error = round + trunc;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (round < trunc && (round > 0 && trunc > 0))
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double r_opt, round_opt, trunc_opt, error_opt;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Compute an optimised stepsize to minimize the total error,
|
|
Packit |
67cb25 |
using the scaling of the estimated truncation error (O(h)) and
|
|
Packit |
67cb25 |
rounding error (O(1/h)). */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double h_opt = h * pow (round / (trunc), 1.0 / 2.0);
|
|
Packit |
67cb25 |
forward_deriv (f, x, h_opt, &r_opt, &round_opt, &trunc_opt);
|
|
Packit |
67cb25 |
error_opt = round_opt + trunc_opt;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Check that the new error is smaller, and that the new derivative
|
|
Packit |
67cb25 |
is consistent with the error bounds of the original estimate. */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (error_opt < error && fabs (r_opt - r_0) < 4.0 * error)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
r_0 = r_opt;
|
|
Packit |
67cb25 |
error = error_opt;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
*result = r_0;
|
|
Packit |
67cb25 |
*abserr = error;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_deriv_backward (const gsl_function * f, double x, double h,
|
|
Packit |
67cb25 |
double *result, double *abserr)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
return gsl_deriv_forward (f, x, -h, result, abserr);
|
|
Packit |
67cb25 |
}
|