htswanalysis/MACS/lib/gsl/gsl-1.11/deriv/deriv.c

   1 /* deriv/deriv.c
   2  *
   3  * Copyright (C) 2004, 2007 Brian Gough
   4  *
   5  * This program is free software; you can redistribute it and/or modify
   6  * it under the terms of the GNU General Public License as published by
   7  * the Free Software Foundation; either version 3 of the License, or (at
   8  * your option) any later version.
   9  *
  10  * This program is distributed in the hope that it will be useful, but
  11  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13  * General Public License for more details.
  14  *
  15  * You should have received a copy of the GNU General Public License
  16  * along with this program; if not, write to the Free Software
  17  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
  18  */
  19
  20 #include <config.h>
  21 #include <stdlib.h>
  22 #include <gsl/gsl_math.h>
  23 #include <gsl/gsl_errno.h>
  24 #include <gsl/gsl_deriv.h>
  25
  26 static void
  27 central_deriv (const gsl_function * f, double x, double h,
  28                double *result, double *abserr_round, double *abserr_trunc)
  29 {
  30   /* Compute the derivative using the 5-point rule (x-h, x-h/2, x,
  31      x+h/2, x+h). Note that the central point is not used.
  32
  33      Compute the error using the difference between the 5-point and
  34      the 3-point rule (x-h,x,x+h). Again the central point is not
  35      used. */
  36
  37   double fm1 = GSL_FN_EVAL (f, x - h);
  38   double fp1 = GSL_FN_EVAL (f, x + h);
  39
  40   double fmh = GSL_FN_EVAL (f, x - h / 2);
  41   double fph = GSL_FN_EVAL (f, x + h / 2);
  42
  43   double r3 = 0.5 * (fp1 - fm1);
  44   double r5 = (4.0 / 3.0) * (fph - fmh) - (1.0 / 3.0) * r3;
  45
  46   double e3 = (fabs (fp1) + fabs (fm1)) * GSL_DBL_EPSILON;
  47   double e5 = 2.0 * (fabs (fph) + fabs (fmh)) * GSL_DBL_EPSILON + e3;
  48
  49   /* The next term is due to finite precision in x+h = O (eps * x) */
  50
  51   double dy = GSL_MAX (fabs (r3 / h), fabs (r5 / h)) *(fabs (x) / h) * GSL_DBL_EPSILON;
  52
  53   /* The truncation error in the r5 approximation itself is O(h^4).
  54      However, for safety, we estimate the error from r5-r3, which is
  55      O(h^2).  By scaling h we will minimise this estimated error, not
  56      the actual truncation error in r5. */
  57
  58   *result = r5 / h;
  59   *abserr_trunc = fabs ((r5 - r3) / h); /* Estimated truncation error O(h^2) */
  60   *abserr_round = fabs (e5 / h) + dy;   /* Rounding error (cancellations) */
  61 }
  62
  63 int
  64 gsl_deriv_central (const gsl_function * f, double x, double h,
  65                    double *result, double *abserr)
  66 {
  67   double r_0, round, trunc, error;
  68   central_deriv (f, x, h, &r_0, &round, &trunc);
  69   error = round + trunc;
  70
  71   if (round < trunc && (round > 0 && trunc > 0))
  72     {
  73       double r_opt, round_opt, trunc_opt, error_opt;
  74
  75       /* Compute an optimised stepsize to minimize the total error,
  76          using the scaling of the truncation error (O(h^2)) and
  77          rounding error (O(1/h)). */
  78
  79       double h_opt = h * pow (round / (2.0 * trunc), 1.0 / 3.0);
  80       central_deriv (f, x, h_opt, &r_opt, &round_opt, &trunc_opt);
  81       error_opt = round_opt + trunc_opt;
  82
  83       /* Check that the new error is smaller, and that the new derivative
  84          is consistent with the error bounds of the original estimate. */
  85
  86       if (error_opt < error && fabs (r_opt - r_0) < 4.0 * error)
  87         {
  88           r_0 = r_opt;
  89           error = error_opt;
  90         }
  91     }
  92
  93   *result = r_0;
  94   *abserr = error;
  95
  96   return GSL_SUCCESS;
  97 }
  98
  99
 100 static void
 101 forward_deriv (const gsl_function * f, double x, double h,
 102                double *result, double *abserr_round, double *abserr_trunc)
 103 {
 104   /* Compute the derivative using the 4-point rule (x+h/4, x+h/2,
 105      x+3h/4, x+h).
 106
 107      Compute the error using the difference between the 4-point and
 108      the 2-point rule (x+h/2,x+h).  */
 109
 110   double f1 = GSL_FN_EVAL (f, x + h / 4.0);
 111   double f2 = GSL_FN_EVAL (f, x + h / 2.0);
 112   double f3 = GSL_FN_EVAL (f, x + (3.0 / 4.0) * h);
 113   double f4 = GSL_FN_EVAL (f, x + h);
 114
 115   double r2 = 2.0*(f4 - f2);
 116   double r4 = (22.0 / 3.0) * (f4 - f3) - (62.0 / 3.0) * (f3 - f2) +
 117     (52.0 / 3.0) * (f2 - f1);
 118
 119   /* Estimate the rounding error for r4 */
 120
 121   double e4 = 2 * 20.67 * (fabs (f4) + fabs (f3) + fabs (f2) + fabs (f1)) * GSL_DBL_EPSILON;
 122
 123   /* The next term is due to finite precision in x+h = O (eps * x) */
 124
 125   double dy = GSL_MAX (fabs (r2 / h), fabs (r4 / h)) * fabs (x / h) * GSL_DBL_EPSILON;
 126
 127   /* The truncation error in the r4 approximation itself is O(h^3).
 128      However, for safety, we estimate the error from r4-r2, which is
 129      O(h).  By scaling h we will minimise this estimated error, not
 130      the actual truncation error in r4. */
 131
 132   *result = r4 / h;
 133   *abserr_trunc = fabs ((r4 - r2) / h); /* Estimated truncation error O(h) */
 134   *abserr_round = fabs (e4 / h) + dy;
 135 }
 136
 137 int
 138 gsl_deriv_forward (const gsl_function * f, double x, double h,
 139                    double *result, double *abserr)
 140 {
 141   double r_0, round, trunc, error;
 142   forward_deriv (f, x, h, &r_0, &round, &trunc);
 143   error = round + trunc;
 144
 145   if (round < trunc && (round > 0 && trunc > 0))
 146     {
 147       double r_opt, round_opt, trunc_opt, error_opt;
 148
 149       /* Compute an optimised stepsize to minimize the total error,
 150          using the scaling of the estimated truncation error (O(h)) and
 151          rounding error (O(1/h)). */
 152
 153       double h_opt = h * pow (round / (trunc), 1.0 / 2.0);
 154       forward_deriv (f, x, h_opt, &r_opt, &round_opt, &trunc_opt);
 155       error_opt = round_opt + trunc_opt;
 156
 157       /* Check that the new error is smaller, and that the new derivative
 158          is consistent with the error bounds of the original estimate. */
 159
 160       if (error_opt < error && fabs (r_opt - r_0) < 4.0 * error)
 161         {
 162           r_0 = r_opt;
 163           error = error_opt;
 164         }
 165     }
 166
 167   *result = r_0;
 168   *abserr = error;
 169
 170   return GSL_SUCCESS;
 171 }
 172
 173 int
 174 gsl_deriv_backward (const gsl_function * f, double x, double h,
 175                     double *result, double *abserr)
 176 {
 177   return gsl_deriv_forward (f, x, -h, result, abserr);
 178 }