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

   1 /* ode-initval/rkf45.c
   2  *
   3  * Copyright (C) 2001, 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 /* Runge-Kutta-Fehlberg 4(5)*/
  21
  22 /* Reference eg. Hairer, E., Norsett S.P., Wanner, G. Solving ordinary
  23    differential equations I, Nonstiff Problems, 2nd revised edition,
  24    Springer, 2000.
  25 */
  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 /* Runge-Kutta-Fehlberg coefficients. Zero elements left out */
  36
  37 static const double ah[] = { 1.0/4.0, 3.0/8.0, 12.0/13.0, 1.0, 1.0/2.0 };
  38 static const double b3[] = { 3.0/32.0, 9.0/32.0 };
  39 static const double b4[] = { 1932.0/2197.0, -7200.0/2197.0, 7296.0/2197.0};
  40 static const double b5[] = { 8341.0/4104.0, -32832.0/4104.0, 29440.0/4104.0, -845.0/4104.0};
  41 static const double b6[] = { -6080.0/20520.0, 41040.0/20520.0, -28352.0/20520.0, 9295.0/20520.0, -5643.0/20520.0};
  42
  43 static const double c1 = 902880.0/7618050.0;
  44 static const double c3 = 3953664.0/7618050.0;
  45 static const double c4 = 3855735.0/7618050.0;
  46 static const double c5 = -1371249.0/7618050.0;
  47 static const double c6 = 277020.0/7618050.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   1.0 / 360.0,
  54   0.0,
  55   -128.0 / 4275.0,
  56   -2197.0 / 75240.0,
  57   1.0 / 50.0,
  58   2.0 / 55.0
  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 rkf45_state_t;
  73
  74 static void *
  75 rkf45_alloc (size_t dim)
  76 {
  77   rkf45_state_t *state = (rkf45_state_t *) malloc (sizeof (rkf45_state_t));
  78
  79   if (state == 0)
  80     {
  81       GSL_ERROR_NULL ("failed to allocate space for rkf45_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 rkf45_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   rkf45_state_t *state = (rkf45_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] +  ah[0] * 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] + c5 * k5[i] + c6 * k6[i];
 292       y[i] += h * d_i;
 293     }
 294
 295   /* Derivatives at output */
 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
 306           return s;
 307         }
 308     }
 309
 310   /* difference between 4th and 5th order */
 311   for (i = 0; i < dim; i++)
 312     {
 313       yerr[i] = h * (ec[1] * k1[i] + ec[3] * k3[i] + ec[4] * k4[i]
 314                      + ec[5] * k5[i] + ec[6] * k6[i]);
 315     }
 316
 317   return GSL_SUCCESS;
 318 }
 319
 320
 321 static int
 322 rkf45_reset (void *vstate, size_t dim)
 323 {
 324   rkf45_state_t *state = (rkf45_state_t *) vstate;
 325
 326   DBL_ZERO_MEMSET (state->k1, dim);
 327   DBL_ZERO_MEMSET (state->k2, dim);
 328   DBL_ZERO_MEMSET (state->k3, dim);
 329   DBL_ZERO_MEMSET (state->k4, dim);
 330   DBL_ZERO_MEMSET (state->k5, dim);
 331   DBL_ZERO_MEMSET (state->k6, dim);
 332   DBL_ZERO_MEMSET (state->ytmp, dim);
 333   DBL_ZERO_MEMSET (state->y0, dim);
 334
 335   return GSL_SUCCESS;
 336 }
 337
 338 static unsigned int
 339 rkf45_order (void *vstate)
 340 {
 341   rkf45_state_t *state = (rkf45_state_t *) vstate;
 342   state = 0; /* prevent warnings about unused parameters */
 343   return 5;
 344 }
 345
 346 static void
 347 rkf45_free (void *vstate)
 348 {
 349   rkf45_state_t *state = (rkf45_state_t *) vstate;
 350
 351   free (state->ytmp);
 352   free (state->y0);
 353   free (state->k6);
 354   free (state->k5);
 355   free (state->k4);
 356   free (state->k3);
 357   free (state->k2);
 358   free (state->k1);
 359   free (state);
 360 }
 361
 362 static const gsl_odeiv_step_type rkf45_type = { "rkf45",        /* name */
 363   1,                            /* can use dydt_in */
 364   0,                            /* gives exact dydt_out */
 365   &rkf45_alloc,
 366   &rkf45_apply,
 367   &rkf45_reset,
 368   &rkf45_order,
 369   &rkf45_free
 370 };
 371
 372 const gsl_odeiv_step_type *gsl_odeiv_step_rkf45 = &rkf45_type;