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

   1 /* ode-initval/rk8pd.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 8(9), Prince-Dormand
  21  *
  22  * High Order Embedded Runge-Kutta Formulae
  23  * P.J. Prince and J.R. Dormand
  24  * J. Comp. Appl. Math.,7, pp. 67-75, 1981
  25  */
  26
  27 /* Author:  G. Jungman
  28  */
  29 #include <config.h>
  30 #include <stdlib.h>
  31 #include <string.h>
  32 #include <gsl/gsl_errno.h>
  33 #include <gsl/gsl_odeiv.h>
  34
  35 #include "odeiv_util.h"
  36
  37 /* Prince-Dormand constants */
  38
  39 static const double Abar[] = {
  40   14005451.0 / 335480064.0,
  41   0.0,
  42   0.0,
  43   0.0,
  44   0.0,
  45   -59238493.0 / 1068277825.0,
  46   181606767.0 / 758867731.0,
  47   561292985.0 / 797845732.0,
  48   -1041891430.0 / 1371343529.0,
  49   760417239.0 / 1151165299.0,
  50   118820643.0 / 751138087.0,
  51   -528747749.0 / 2220607170.0,
  52   1.0 / 4.0
  53 };
  54
  55 static const double A[] = {
  56   13451932.0 / 455176623.0,
  57   0.0,
  58   0.0,
  59   0.0,
  60   0.0,
  61   -808719846.0 / 976000145.0,
  62   1757004468.0 / 5645159321.0,
  63   656045339.0 / 265891186.0,
  64   -3867574721.0 / 1518517206.0,
  65   465885868.0 / 322736535.0,
  66   53011238.0 / 667516719.0,
  67   2.0 / 45.0
  68 };
  69
  70 static const double ah[] = {
  71   1.0 / 18.0,
  72   1.0 / 12.0,
  73   1.0 / 8.0,
  74   5.0 / 16.0,
  75   3.0 / 8.0,
  76   59.0 / 400.0,
  77   93.0 / 200.0,
  78   5490023248.0 / 9719169821.0,
  79   13.0 / 20.0,
  80   1201146811.0 / 1299019798.0
  81 };
  82
  83 static const double b21 = 1.0 / 18.0;
  84 static const double b3[] = { 1.0 / 48.0, 1.0 / 16.0 };
  85 static const double b4[] = { 1.0 / 32.0, 0.0, 3.0 / 32.0 };
  86 static const double b5[] = { 5.0 / 16.0, 0.0, -75.0 / 64.0, 75.0 / 64.0 };
  87 static const double b6[] = { 3.0 / 80.0, 0.0, 0.0, 3.0 / 16.0, 3.0 / 20.0 };
  88 static const double b7[] = {
  89   29443841.0 / 614563906.0,
  90   0.0,
  91   0.0,
  92   77736538.0 / 692538347.0,
  93   -28693883.0 / 1125000000.0,
  94   23124283.0 / 1800000000.0
  95 };
  96 static const double b8[] = {
  97   16016141.0 / 946692911.0,
  98   0.0,
  99   0.0,
 100   61564180.0 / 158732637.0,
 101   22789713.0 / 633445777.0,
 102   545815736.0 / 2771057229.0,
 103   -180193667.0 / 1043307555.0
 104 };
 105 static const double b9[] = {
 106   39632708.0 / 573591083.0,
 107   0.0,
 108   0.0,
 109   -433636366.0 / 683701615.0,
 110   -421739975.0 / 2616292301.0,
 111   100302831.0 / 723423059.0,
 112   790204164.0 / 839813087.0,
 113   800635310.0 / 3783071287.0
 114 };
 115 static const double b10[] = {
 116   246121993.0 / 1340847787.0,
 117   0.0,
 118   0.0,
 119   -37695042795.0 / 15268766246.0,
 120   -309121744.0 / 1061227803.0,
 121   -12992083.0 / 490766935.0,
 122   6005943493.0 / 2108947869.0,
 123   393006217.0 / 1396673457.0,
 124   123872331.0 / 1001029789.0
 125 };
 126 static const double b11[] = {
 127   -1028468189.0 / 846180014.0,
 128   0.0,
 129   0.0,
 130   8478235783.0 / 508512852.0,
 131   1311729495.0 / 1432422823.0,
 132   -10304129995.0 / 1701304382.0,
 133   -48777925059.0 / 3047939560.0,
 134   15336726248.0 / 1032824649.0,
 135   -45442868181.0 / 3398467696.0,
 136   3065993473.0 / 597172653.0
 137 };
 138 static const double b12[] = {
 139   185892177.0 / 718116043.0,
 140   0.0,
 141   0.0,
 142   -3185094517.0 / 667107341.0,
 143   -477755414.0 / 1098053517.0,
 144   -703635378.0 / 230739211.0,
 145   5731566787.0 / 1027545527.0,
 146   5232866602.0 / 850066563.0,
 147   -4093664535.0 / 808688257.0,
 148   3962137247.0 / 1805957418.0,
 149   65686358.0 / 487910083.0
 150 };
 151 static const double b13[] = {
 152   403863854.0 / 491063109.0,
 153   0.0,
 154   0.0,
 155   -5068492393.0 / 434740067.0,
 156   -411421997.0 / 543043805.0,
 157   652783627.0 / 914296604.0,
 158   11173962825.0 / 925320556.0,
 159   -13158990841.0 / 6184727034.0,
 160   3936647629.0 / 1978049680.0,
 161   -160528059.0 / 685178525.0,
 162   248638103.0 / 1413531060.0,
 163   0.0
 164 };
 165
 166 typedef struct
 167 {
 168   double *k[13];
 169   double *ytmp;
 170   double *y0;
 171 }
 172 rk8pd_state_t;
 173
 174 static void *
 175 rk8pd_alloc (size_t dim)
 176 {
 177   rk8pd_state_t *state = (rk8pd_state_t *) malloc (sizeof (rk8pd_state_t));
 178   int i, j;
 179
 180   if (state == 0)
 181     {
 182       GSL_ERROR_NULL ("failed to allocate space for rk8pd_state", GSL_ENOMEM);
 183     }
 184
 185   state->ytmp = (double *) malloc (dim * sizeof (double));
 186
 187   if (state->ytmp == 0)
 188     {
 189       free (state);
 190       GSL_ERROR_NULL ("failed to allocate space for ytmp", GSL_ENOMEM);
 191     }
 192
 193   state->y0 = (double *) malloc (dim * sizeof (double));
 194
 195   if (state->y0 == 0)
 196     {
 197       free (state->ytmp);
 198       free (state);
 199       GSL_ERROR_NULL ("failed to allocate space for y0", GSL_ENOMEM);
 200     }
 201
 202   for (i = 0; i < 13; i++)
 203     {
 204       state->k[i] = (double *) malloc (dim * sizeof (double));
 205
 206       if (state->k[i] == 0)
 207         {
 208           for (j = 0; j < i; j++)
 209             {
 210               free (state->k[j]);
 211             }
 212           free (state->y0);
 213           free (state->ytmp);
 214           free (state);
 215           GSL_ERROR_NULL ("failed to allocate space for k's", GSL_ENOMEM);
 216         }
 217     }
 218
 219   return state;
 220 }
 221
 222
 223 static int
 224 rk8pd_apply (void *vstate,
 225              size_t dim,
 226              double t,
 227              double h,
 228              double y[],
 229              double yerr[],
 230              const double dydt_in[],
 231              double dydt_out[], const gsl_odeiv_system * sys)
 232 {
 233   rk8pd_state_t *state = (rk8pd_state_t *) vstate;
 234
 235   size_t i;
 236
 237   double *const ytmp = state->ytmp;
 238   double *const y0 = state->y0;
 239   /* Note that k1 is stored in state->k[0] due to zero-based indexing */
 240   double *const k1 = state->k[0];
 241   double *const k2 = state->k[1];
 242   double *const k3 = state->k[2];
 243   double *const k4 = state->k[3];
 244   double *const k5 = state->k[4];
 245   double *const k6 = state->k[5];
 246   double *const k7 = state->k[6];
 247   double *const k8 = state->k[7];
 248   double *const k9 = state->k[8];
 249   double *const k10 = state->k[9];
 250   double *const k11 = state->k[10];
 251   double *const k12 = state->k[11];
 252   double *const k13 = state->k[12];
 253
 254   DBL_MEMCPY (y0, y, dim);
 255
 256   /* k1 step */
 257   if (dydt_in != NULL)
 258     {
 259       DBL_MEMCPY (k1, dydt_in, dim);
 260     }
 261   else
 262     {
 263       int s = GSL_ODEIV_FN_EVAL (sys, t, y, k1);
 264
 265       if (s != GSL_SUCCESS)
 266         {
 267           return s;
 268         }
 269     }
 270
 271   for (i = 0; i < dim; i++)
 272     ytmp[i] = y[i] + b21 * h * k1[i];
 273
 274   /* k2 step */
 275   {
 276     int s = GSL_ODEIV_FN_EVAL (sys, t + ah[0] * h, ytmp, k2);
 277
 278     if (s != GSL_SUCCESS)
 279       {
 280         return s;
 281       }
 282   }
 283
 284   for (i = 0; i < dim; i++)
 285     ytmp[i] = y[i] + h * (b3[0] * k1[i] + b3[1] * k2[i]);
 286
 287   /* k3 step */
 288   {
 289     int s = GSL_ODEIV_FN_EVAL (sys, t + ah[1] * h, ytmp, k3);
 290
 291     if (s != GSL_SUCCESS)
 292       {
 293         return s;
 294       }
 295   }
 296
 297   for (i = 0; i < dim; i++)
 298     ytmp[i] = y[i] + h * (b4[0] * k1[i] + b4[2] * k3[i]);
 299
 300   /* k4 step */
 301   {
 302     int s = GSL_ODEIV_FN_EVAL (sys, t + ah[2] * h, ytmp, k4);
 303
 304     if (s != GSL_SUCCESS)
 305       {
 306         return s;
 307       }
 308   }
 309
 310   for (i = 0; i < dim; i++)
 311     ytmp[i] = y[i] + h * (b5[0] * k1[i] + b5[2] * k3[i] + b5[3] * k4[i]);
 312
 313   /* k5 step */
 314   {
 315     int s = GSL_ODEIV_FN_EVAL (sys, t + ah[3] * h, ytmp, k5);
 316
 317     if (s != GSL_SUCCESS)
 318       {
 319         return s;
 320       }
 321   }
 322
 323   for (i = 0; i < dim; i++)
 324     ytmp[i] = y[i] + h * (b6[0] * k1[i] + b6[3] * k4[i] + b6[4] * k5[i]);
 325
 326   /* k6 step */
 327   {
 328     int s = GSL_ODEIV_FN_EVAL (sys, t + ah[4] * h, ytmp, k6);
 329
 330     if (s != GSL_SUCCESS)
 331       {
 332         return s;
 333       }
 334   }
 335
 336   for (i = 0; i < dim; i++)
 337     ytmp[i] =
 338       y[i] + h * (b7[0] * k1[i] + b7[3] * k4[i] + b7[4] * k5[i] +
 339                   b7[5] * k6[i]);
 340
 341   /* k7 step */
 342   {
 343     int s = GSL_ODEIV_FN_EVAL (sys, t + ah[5] * h, ytmp, k7);
 344
 345     if (s != GSL_SUCCESS)
 346       {
 347         return s;
 348       }
 349   }
 350
 351   for (i = 0; i < dim; i++)
 352     ytmp[i] =
 353       y[i] + h * (b8[0] * k1[i] + b8[3] * k4[i] + b8[4] * k5[i] +
 354                   b8[5] * k6[i] + b8[6] * k7[i]);
 355
 356   /* k8 step */
 357   {
 358     int s = GSL_ODEIV_FN_EVAL (sys, t + ah[6] * h, ytmp, k8);
 359
 360     if (s != GSL_SUCCESS)
 361       {
 362         return s;
 363       }
 364   }
 365
 366   for (i = 0; i < dim; i++)
 367     ytmp[i] =
 368       y[i] + h * (b9[0] * k1[i] + b9[3] * k4[i] + b9[4] * k5[i] +
 369                   b9[5] * k6[i] + b9[6] * k7[i] + b9[7] * k8[i]);
 370
 371   /* k9 step */
 372   {
 373     int s = GSL_ODEIV_FN_EVAL (sys, t + ah[7] * h, ytmp, k9);
 374
 375     if (s != GSL_SUCCESS)
 376       {
 377         return s;
 378       }
 379   }
 380
 381   for (i = 0; i < dim; i++)
 382     ytmp[i] =
 383       y[i] + h * (b10[0] * k1[i] + b10[3] * k4[i] + b10[4] * k5[i] +
 384                   b10[5] * k6[i] + b10[6] * k7[i] + b10[7] * k8[i] +
 385                   b10[8] * k9[i]);
 386
 387   /* k10 step */
 388   {
 389     int s = GSL_ODEIV_FN_EVAL (sys, t + ah[8] * h, ytmp, k10);
 390
 391     if (s != GSL_SUCCESS)
 392       {
 393         return s;
 394       }
 395   }
 396
 397   for (i = 0; i < dim; i++)
 398     ytmp[i] =
 399       y[i] + h * (b11[0] * k1[i] + b11[3] * k4[i] + b11[4] * k5[i] +
 400                   b11[5] * k6[i] + b11[6] * k7[i] + b11[7] * k8[i] +
 401                   b11[8] * k9[i] + b11[9] * k10[i]);
 402
 403   /* k11 step */
 404   {
 405     int s = GSL_ODEIV_FN_EVAL (sys, t + ah[9] * h, ytmp, k11);
 406
 407     if (s != GSL_SUCCESS)
 408       {
 409         return s;
 410       }
 411   }
 412
 413   for (i = 0; i < dim; i++)
 414     ytmp[i] =
 415       y[i] + h * (b12[0] * k1[i] + b12[3] * k4[i] + b12[4] * k5[i] +
 416                   b12[5] * k6[i] + b12[6] * k7[i] + b12[7] * k8[i] +
 417                   b12[8] * k9[i] + b12[9] * k10[i] + b12[10] * k11[i]);
 418
 419   /* k12 step */
 420   {
 421     int s = GSL_ODEIV_FN_EVAL (sys, t + h, ytmp, k12);
 422
 423     if (s != GSL_SUCCESS)
 424       {
 425         return s;
 426       }
 427   }
 428
 429   for (i = 0; i < dim; i++)
 430     ytmp[i] =
 431       y[i] + h * (b13[0] * k1[i] + b13[3] * k4[i] + b13[4] * k5[i] +
 432                   b13[5] * k6[i] + b13[6] * k7[i] + b13[7] * k8[i] +
 433                   b13[8] * k9[i] + b13[9] * k10[i] + b13[10] * k11[i] +
 434                   b13[11] * k12[i]);
 435
 436   /* k13 step */
 437   {
 438     int s = GSL_ODEIV_FN_EVAL (sys, t + h, ytmp, k13);
 439
 440     if (s != GSL_SUCCESS)
 441       {
 442         return s;
 443       }
 444   }
 445
 446   /* final sum  */
 447   for (i = 0; i < dim; i++)
 448     {
 449       const double ksum8 =
 450         Abar[0] * k1[i] + Abar[5] * k6[i] + Abar[6] * k7[i] +
 451         Abar[7] * k8[i] + Abar[8] * k9[i] + Abar[9] * k10[i] +
 452         Abar[10] * k11[i] + Abar[11] * k12[i] + Abar[12] * k13[i];
 453       y[i] += h * ksum8;
 454     }
 455
 456   /* Evaluate dydt_out[]. */
 457
 458   if (dydt_out != NULL)
 459     {
 460       int s = GSL_ODEIV_FN_EVAL (sys, t + h, y, dydt_out);
 461
 462       if (s != GSL_SUCCESS)
 463         {
 464           /* Restore initial values */
 465           DBL_MEMCPY (y, y0, dim);
 466           return s;
 467         }
 468     }
 469
 470   /* error estimate */
 471   for (i = 0; i < dim; i++)
 472     {
 473       const double ksum8 =
 474         Abar[0] * k1[i] + Abar[5] * k6[i] + Abar[6] * k7[i] +
 475         Abar[7] * k8[i] + Abar[8] * k9[i] + Abar[9] * k10[i] +
 476         Abar[10] * k11[i] + Abar[11] * k12[i] + Abar[12] * k13[i];
 477       const double ksum7 =
 478         A[0] * k1[i] + A[5] * k6[i] + A[6] * k7[i] + A[7] * k8[i] +
 479         A[8] * k9[i] + A[9] * k10[i] + A[10] * k11[i] + A[11] * k12[i];
 480       yerr[i] = h * (ksum7 - ksum8);
 481     }
 482
 483   return GSL_SUCCESS;
 484 }
 485
 486 static int
 487 rk8pd_reset (void *vstate, size_t dim)
 488 {
 489   rk8pd_state_t *state = (rk8pd_state_t *) vstate;
 490
 491   int i;
 492
 493   for (i = 0; i < 13; i++)
 494     {
 495       DBL_ZERO_MEMSET (state->k[i], dim);
 496     }
 497
 498   DBL_ZERO_MEMSET (state->y0, dim);
 499   DBL_ZERO_MEMSET (state->ytmp, dim);
 500
 501   return GSL_SUCCESS;
 502 }
 503
 504 static unsigned int
 505 rk8pd_order (void *vstate)
 506 {
 507   rk8pd_state_t *state = (rk8pd_state_t *) vstate;
 508   state = 0; /* prevent warnings about unused parameters */
 509   return 8;
 510 }
 511
 512 static void
 513 rk8pd_free (void *vstate)
 514 {
 515   rk8pd_state_t *state = (rk8pd_state_t *) vstate;
 516   int i;
 517
 518   for (i = 0; i < 13; i++)
 519     {
 520       free (state->k[i]);
 521     }
 522   free (state->y0);
 523   free (state->ytmp);
 524   free (state);
 525 }
 526
 527 static const gsl_odeiv_step_type rk8pd_type = { "rk8pd",        /* name */
 528   1,                            /* can use dydt_in */
 529   1,                            /* gives exact dydt_out */
 530   &rk8pd_alloc,
 531   &rk8pd_apply,
 532   &rk8pd_reset,
 533   &rk8pd_order,
 534   &rk8pd_free
 535 };
 536
 537 const gsl_odeiv_step_type *gsl_odeiv_step_rk8pd = &rk8pd_type;