htswanalysis/MACS/lib/gsl/gsl-1.11/ode-initval/rkck.c

   1 /* ode-initval/rkck.c
   2  *
   3  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
   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 /* Runge-Kutta 4(5), Cash-Karp */
  21
  22 /* Reference: Cash, J.R., Karp, A.H., ACM Transactions of Mathematical
  23    Software, vol. 16 (1990) 201-222. */
  24
  25 /* Author:  G. Jungman
  26  */
  27 #include <config.h>
  28 #include <stdlib.h>
  29 #include <string.h>
  30 #include <gsl/gsl_errno.h>
  31 #include <gsl/gsl_odeiv.h>
  32
  33 #include "odeiv_util.h"
  34
  35 /* Cash-Karp constants */
  36 static const double ah[] = { 1.0 / 5.0, 0.3, 3.0 / 5.0, 1.0, 7.0 / 8.0 };
  37 static const double b21 = 1.0 / 5.0;
  38 static const double b3[] = { 3.0 / 40.0, 9.0 / 40.0 };
  39 static const double b4[] = { 0.3, -0.9, 1.2 };
  40 static const double b5[] = { -11.0 / 54.0, 2.5, -70.0 / 27.0, 35.0 / 27.0 };
  41 static const double b6[] =
  42   { 1631.0 / 55296.0, 175.0 / 512.0, 575.0 / 13824.0, 44275.0 / 110592.0,
  43     253.0 / 4096.0 };
  44 static const double c1 = 37.0 / 378.0;
  45 static const double c3 = 250.0 / 621.0;
  46 static const double c4 = 125.0 / 594.0;
  47 static const double c6 = 512.0 / 1771.0;
  48
  49 /* These are the differences of fifth and fourth order coefficients
  50    for error estimation */
  51
  52 static const double ec[] = { 0.0,
  53   37.0 / 378.0 - 2825.0 / 27648.0,
  54   0.0,
  55   250.0 / 621.0 - 18575.0 / 48384.0,
  56   125.0 / 594.0 - 13525.0 / 55296.0,
  57   -277.0 / 14336.0,
  58   512.0 / 1771.0 - 0.25
  59 };
  60
  61 typedef struct
  62 {
  63   double *k1;
  64   double *k2;
  65   double *k3;
  66   double *k4;
  67   double *k5;
  68   double *k6;
  69   double *y0;
  70   double *ytmp;
  71 }
  72 rkck_state_t;
  73
  74 static void *
  75 rkck_alloc (size_t dim)
  76 {
  77   rkck_state_t *state = (rkck_state_t *) malloc (sizeof (rkck_state_t));
  78
  79   if (state == 0)
  80     {
  81       GSL_ERROR_NULL ("failed to allocate space for rkck_state", GSL_ENOMEM);
  82     }
  83
  84   state->k1 = (double *) malloc (dim * sizeof (double));
  85
  86   if (state->k1 == 0)
  87     {
  88       free (state);
  89       GSL_ERROR_NULL ("failed to allocate space for k1", GSL_ENOMEM);
  90     }
  91
  92   state->k2 = (double *) malloc (dim * sizeof (double));
  93
  94   if (state->k2 == 0)
  95     {
  96       free (state->k1);
  97       free (state);
  98       GSL_ERROR_NULL ("failed to allocate space for k2", GSL_ENOMEM);
  99     }
 100
 101   state->k3 = (double *) malloc (dim * sizeof (double));
 102
 103   if (state->k3 == 0)
 104     {
 105       free (state->k2);
 106       free (state->k1);
 107       free (state);
 108       GSL_ERROR_NULL ("failed to allocate space for k3", GSL_ENOMEM);
 109     }
 110
 111   state->k4 = (double *) malloc (dim * sizeof (double));
 112
 113   if (state->k4 == 0)
 114     {
 115       free (state->k3);
 116       free (state->k2);
 117       free (state->k1);
 118       free (state);
 119       GSL_ERROR_NULL ("failed to allocate space for k4", GSL_ENOMEM);
 120     }
 121
 122   state->k5 = (double *) malloc (dim * sizeof (double));
 123
 124   if (state->k5 == 0)
 125     {
 126       free (state->k4);
 127       free (state->k3);
 128       free (state->k2);
 129       free (state->k1);
 130       free (state);
 131       GSL_ERROR_NULL ("failed to allocate space for k5", GSL_ENOMEM);
 132     }
 133
 134   state->k6 = (double *) malloc (dim * sizeof (double));
 135
 136   if (state->k6 == 0)
 137     {
 138       free (state->k5);
 139       free (state->k4);
 140       free (state->k3);
 141       free (state->k2);
 142       free (state->k1);
 143       free (state);
 144       GSL_ERROR_NULL ("failed to allocate space for k6", GSL_ENOMEM);
 145     }
 146
 147   state->y0 = (double *) malloc (dim * sizeof (double));
 148
 149   if (state->y0 == 0)
 150     {
 151       free (state->k6);
 152       free (state->k5);
 153       free (state->k4);
 154       free (state->k3);
 155       free (state->k2);
 156       free (state->k1);
 157       free (state);
 158       GSL_ERROR_NULL ("failed to allocate space for y0", GSL_ENOMEM);
 159     }
 160
 161   state->ytmp = (double *) malloc (dim * sizeof (double));
 162
 163   if (state->ytmp == 0)
 164     {
 165       free (state->y0);
 166       free (state->k6);
 167       free (state->k5);
 168       free (state->k4);
 169       free (state->k3);
 170       free (state->k2);
 171       free (state->k1);
 172       free (state);
 173       GSL_ERROR_NULL ("failed to allocate space for ytmp", GSL_ENOMEM);
 174     }
 175
 176   return state;
 177 }
 178
 179
 180 static int
 181 rkck_apply (void *vstate,
 182             size_t dim,
 183             double t,
 184             double h,
 185             double y[],
 186             double yerr[],
 187             const double dydt_in[],
 188             double dydt_out[], const gsl_odeiv_system * sys)
 189 {
 190   rkck_state_t *state = (rkck_state_t *) vstate;
 191
 192   size_t i;
 193
 194   double *const k1 = state->k1;
 195   double *const k2 = state->k2;
 196   double *const k3 = state->k3;
 197   double *const k4 = state->k4;
 198   double *const k5 = state->k5;
 199   double *const k6 = state->k6;
 200   double *const ytmp = state->ytmp;
 201   double *const y0 = state->y0;
 202
 203   DBL_MEMCPY (y0, y, dim);
 204
 205   /* k1 step */
 206   if (dydt_in != NULL)
 207     {
 208       DBL_MEMCPY (k1, dydt_in, dim);
 209     }
 210   else
 211     {
 212       int s = GSL_ODEIV_FN_EVAL (sys, t, y, k1);
 213
 214       if (s != GSL_SUCCESS)
 215         {
 216           return s;
 217         }
 218     }
 219
 220   for (i = 0; i < dim; i++)
 221     ytmp[i] = y[i] + b21 * h * k1[i];
 222
 223   /* k2 step */
 224   {
 225     int s = GSL_ODEIV_FN_EVAL (sys, t + ah[0] * h, ytmp, k2);
 226
 227     if (s != GSL_SUCCESS)
 228       {
 229         return s;
 230       }
 231   }
 232
 233   for (i = 0; i < dim; i++)
 234     ytmp[i] = y[i] + h * (b3[0] * k1[i] + b3[1] * k2[i]);
 235
 236   /* k3 step */
 237   {
 238     int s = GSL_ODEIV_FN_EVAL (sys, t + ah[1] * h, ytmp, k3);
 239
 240     if (s != GSL_SUCCESS)
 241       {
 242         return s;
 243       }
 244   }
 245
 246   for (i = 0; i < dim; i++)
 247     ytmp[i] = y[i] + h * (b4[0] * k1[i] + b4[1] * k2[i] + b4[2] * k3[i]);
 248
 249   /* k4 step */
 250   {
 251     int s = GSL_ODEIV_FN_EVAL (sys, t + ah[2] * h, ytmp, k4);
 252
 253     if (s != GSL_SUCCESS)
 254       {
 255         return s;
 256       }
 257   }
 258
 259   for (i = 0; i < dim; i++)
 260     ytmp[i] =
 261       y[i] + h * (b5[0] * k1[i] + b5[1] * k2[i] + b5[2] * k3[i] +
 262                   b5[3] * k4[i]);
 263
 264   /* k5 step */
 265   {
 266     int s = GSL_ODEIV_FN_EVAL (sys, t + ah[3] * h, ytmp, k5);
 267
 268     if (s != GSL_SUCCESS)
 269       {
 270         return s;
 271       }
 272   }
 273
 274   for (i = 0; i < dim; i++)
 275     ytmp[i] =
 276       y[i] + h * (b6[0] * k1[i] + b6[1] * k2[i] + b6[2] * k3[i] +
 277                   b6[3] * k4[i] + b6[4] * k5[i]);
 278
 279   /* k6 step and final sum */
 280   {
 281     int s = GSL_ODEIV_FN_EVAL (sys, t + ah[4] * h, ytmp, k6);
 282
 283     if (s != GSL_SUCCESS)
 284       {
 285         return s;
 286       }
 287   }
 288
 289   for (i = 0; i < dim; i++)
 290     {
 291       const double d_i = c1 * k1[i] + c3 * k3[i] + c4 * k4[i] + c6 * k6[i];
 292       y[i] += h * d_i;
 293     }
 294
 295   /* Evaluate dydt_out[]. */
 296
 297   if (dydt_out != NULL)
 298     {
 299       int s = GSL_ODEIV_FN_EVAL (sys, t + h, y, dydt_out);
 300
 301       if (s != GSL_SUCCESS)
 302         {
 303           /* Restore initial values */
 304           DBL_MEMCPY (y, y0, dim);
 305           return s;
 306         }
 307     }
 308
 309   /* difference between 4th and 5th order */
 310   for (i = 0; i < dim; i++)
 311     {
 312       yerr[i] = h * (ec[1] * k1[i] + ec[3] * k3[i] + ec[4] * k4[i]
 313                      + ec[5] * k5[i] + ec[6] * k6[i]);
 314     }
 315
 316   return GSL_SUCCESS;
 317 }
 318
 319
 320 static int
 321 rkck_reset (void *vstate, size_t dim)
 322 {
 323   rkck_state_t *state = (rkck_state_t *) vstate;
 324
 325   DBL_ZERO_MEMSET (state->k1, dim);
 326   DBL_ZERO_MEMSET (state->k2, dim);
 327   DBL_ZERO_MEMSET (state->k3, dim);
 328   DBL_ZERO_MEMSET (state->k4, dim);
 329   DBL_ZERO_MEMSET (state->k5, dim);
 330   DBL_ZERO_MEMSET (state->k6, dim);
 331   DBL_ZERO_MEMSET (state->ytmp, dim);
 332   DBL_ZERO_MEMSET (state->y0, dim);
 333
 334   return GSL_SUCCESS;
 335 }
 336
 337 static unsigned int
 338 rkck_order (void *vstate)
 339 {
 340   rkck_state_t *state = (rkck_state_t *) vstate;
 341   state = 0; /* prevent warnings about unused parameters */
 342   return 5; /* FIXME: should this be 4? */
 343 }
 344
 345 static void
 346 rkck_free (void *vstate)
 347 {
 348   rkck_state_t *state = (rkck_state_t *) vstate;
 349
 350   free (state->ytmp);
 351   free (state->y0);
 352   free (state->k6);
 353   free (state->k5);
 354   free (state->k4);
 355   free (state->k3);
 356   free (state->k2);
 357   free (state->k1);
 358   free (state);
 359 }
 360
 361 static const gsl_odeiv_step_type rkck_type = { "rkck",  /* name */
 362   1,                            /* can use dydt_in */
 363   1,                            /* gives exact dydt_out */
 364   &rkck_alloc,
 365   &rkck_apply,
 366   &rkck_reset,
 367   &rkck_order,
 368   &rkck_free
 369 };
 370
 371 const gsl_odeiv_step_type *gsl_odeiv_step_rkck = &rkck_type;