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

   1 /* ode-initval/test_odeiv.c
   2  *
   3  * Copyright (C) 2004 Tuomo Keskitalo
   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 /* Some functions and tests based on test.c by G. Jungman.
  21 */
  22
  23 #include <config.h>
  24 #include <stdlib.h>
  25 #include <stdio.h>
  26 #include <string.h>
  27 #include <math.h>
  28 #include <gsl/gsl_test.h>
  29 #include <gsl/gsl_errno.h>
  30 #include <gsl/gsl_math.h>
  31 #include <gsl/gsl_matrix.h>
  32 #include <gsl/gsl_linalg.h>
  33 #include <gsl/gsl_ieee_utils.h>
  34 #include <gsl/gsl_odeiv.h>
  35 #include "odeiv_util.h"
  36
  37 /* Maximum number of ODE equations */
  38 #define MAXEQ 4
  39
  40 /* RHS for f=2. Solution y = 2 * t + t0 */
  41
  42 int
  43 rhs_linear (double t, const double y[], double f[], void *params)
  44 {
  45   f[0] = 2.0;
  46
  47   return GSL_SUCCESS;
  48 }
  49
  50 int
  51 jac_linear (double t, const double y[], double *dfdy, double dfdt[],
  52             void *params)
  53 {
  54   dfdy[0] = 0.0;
  55   dfdt[0] = 0.0;
  56
  57   return GSL_SUCCESS;
  58 }
  59
  60 gsl_odeiv_system rhs_func_lin = {
  61   rhs_linear,
  62   jac_linear,
  63   1,
  64   0
  65 };
  66
  67 /* RHS for f=y. Equals y=exp(t) with initial value y(0)=1.0 */
  68
  69 int
  70 rhs_exp (double t, const double y[], double f[], void *params)
  71 {
  72   f[0] = y[0];
  73
  74   return GSL_SUCCESS;
  75 }
  76
  77 int
  78 jac_exp (double t, const double y[], double *dfdy, double dfdt[],
  79          void *params)
  80 {
  81   dfdy[0] = y[0];
  82   dfdt[0] = 0.0;
  83
  84   return GSL_SUCCESS;
  85 }
  86
  87 gsl_odeiv_system rhs_func_exp = {
  88   rhs_exp,
  89   jac_exp,
  90   1,
  91   0
  92 };
  93
  94 /* RHS for f0 = -y1, f1 = y0
  95    equals y = [cos(t), sin(t)] with initial values [1, 0]
  96 */
  97
  98 int
  99 rhs_sin (double t, const double y[], double f[], void *params)
 100 {
 101   f[0] = -y[1];
 102   f[1] = y[0];
 103
 104   return GSL_SUCCESS;
 105 }
 106
 107 int
 108 jac_sin (double t, const double y[], double *dfdy, double dfdt[],
 109          void *params)
 110 {
 111   dfdy[0] = 0.0;
 112   dfdy[1] = -1.0;
 113   dfdy[2] = 1.0;
 114   dfdy[3] = 0.0;
 115
 116   dfdt[0] = 0.0;
 117   dfdt[1] = 0.0;
 118
 119   return GSL_SUCCESS;
 120 }
 121
 122 gsl_odeiv_system rhs_func_sin = {
 123   rhs_sin,
 124   jac_sin,
 125   2,
 126   0
 127 };
 128
 129 /*  Sine/cosine with random failures */
 130
 131 int
 132 rhs_xsin (double t, const double y[], double f[], void *params)
 133 {
 134   static int n = 0;
 135
 136   n++;
 137
 138   if (n > 40 && n < 65) {
 139     f[0] = GSL_NAN;
 140     f[1] = GSL_NAN;
 141     return GSL_EFAILED;
 142   }
 143
 144   f[0] = -y[1];
 145   f[1] = y[0];
 146
 147   return GSL_SUCCESS;
 148 }
 149
 150 int
 151 jac_xsin (double t, const double y[], double *dfdy, double dfdt[],
 152          void *params)
 153 {
 154   static int n = 0;
 155
 156   n++;
 157
 158   if (n > 50 && n < 55) {
 159     dfdy[0] = GSL_NAN;
 160     dfdy[1] = GSL_NAN;
 161     dfdy[2] = GSL_NAN;
 162     dfdy[3] = GSL_NAN;
 163
 164     dfdt[0] = GSL_NAN;
 165     dfdt[1] = GSL_NAN;
 166     return GSL_EFAILED;
 167   }
 168
 169   dfdy[0] = 0.0;
 170   dfdy[1] = -1.0;
 171   dfdy[2] = 1.0;
 172   dfdy[3] = 0.0;
 173
 174   dfdt[0] = 0.0;
 175   dfdt[1] = 0.0;
 176
 177   return GSL_SUCCESS;
 178 }
 179
 180 gsl_odeiv_system rhs_func_xsin = {
 181   rhs_xsin,
 182   jac_xsin,
 183   2,
 184   0
 185 };
 186
 187
 188 /* RHS for classic stiff example
 189    dy0 / dt =  998 * y0 + 1998 * y1    y0(0) = 1.0
 190    dy1 / dt = -999 * y0 - 1999 * y1    y1(0) = 0.0
 191
 192    solution is
 193    y0 = 2 * exp(-t) - exp(-1000 * t)
 194    y1 = - exp(-t) + exp(-1000 * t)
 195 */
 196
 197 int
 198 rhs_stiff (double t, const double y[], double f[], void *params)
 199 {
 200   f[0] = 998.0 * y[0] + 1998.0 * y[1];
 201   f[1] = -999.0 * y[0] - 1999.0 * y[1];
 202
 203   return GSL_SUCCESS;
 204 }
 205
 206 int
 207 jac_stiff (double t, const double y[], double *dfdy, double dfdt[],
 208            void *params)
 209 {
 210   dfdy[0] = 998.0;
 211   dfdy[1] = 1998.0;
 212   dfdy[2] = -999.0;
 213   dfdy[3] = -1999.0;
 214
 215   dfdt[0] = 0.0;
 216   dfdt[1] = 0.0;
 217
 218   return GSL_SUCCESS;
 219 }
 220
 221 gsl_odeiv_system rhs_func_stiff = {
 222   rhs_stiff,
 223   jac_stiff,
 224   2,
 225   0
 226 };
 227
 228 /* van Der Pol oscillator:
 229    f0 = y1                           y0(0) = 1.0
 230    f1 = -y0 + mu * y1 * (1 - y0^2)   y1(0) = 0.0
 231 */
 232
 233 int
 234 rhs_vanderpol (double t, const double y[], double f[], void *params)
 235 {
 236   const double mu = 10.0;
 237
 238   f[0] = y[1];
 239   f[1] = -y[0] + mu * y[1] * (1.0 - y[0]*y[0]);
 240
 241   return GSL_SUCCESS;
 242 }
 243
 244 int
 245 jac_vanderpol (double t, const double y[], double *dfdy, double dfdt[],
 246                void *params)
 247 {
 248   const double mu = 10.0;
 249
 250   dfdy[0] = 0.0;
 251   dfdy[1] = 1.0;
 252   dfdy[2] = -2.0 * mu * y[0] * y[1] - 1.0;
 253   dfdy[3] = mu * (1.0 - y[0] * y[0]);
 254
 255   dfdt[0] = 0.0;
 256   dfdt[1] = 0.0;
 257
 258   return GSL_SUCCESS;
 259 }
 260
 261 gsl_odeiv_system rhs_func_vanderpol = {
 262   rhs_vanderpol,
 263   jac_vanderpol,
 264   2,
 265   0
 266 };
 267
 268 /* The Oregonator - chemical Belusov-Zhabotinskii reaction
 269    y0(0) = 1.0, y1(0) = 2.0, y2(0) = 3.0
 270 */
 271
 272 int
 273 rhs_oregonator (double t, const double y[], double f[], void *params)
 274 {
 275   const double c1=77.27;
 276   const double c2=8.375e-6;
 277   const double c3=0.161;
 278
 279   f[0] = c1 * (y[1] + y[0] * (1 - c2 * y[0] - y[1]));
 280   f[1] = 1/c1 * (y[2] - y[1] * (1 + y[0]));
 281   f[2] = c3 * (y[0] - y[2]);
 282
 283   return GSL_SUCCESS;
 284 }
 285
 286 int
 287 jac_oregonator (double t, const double y[], double *dfdy, double dfdt[],
 288                 void *params)
 289 {
 290   const double c1=77.27;
 291   const double c2=8.375e-6;
 292   const double c3=0.161;
 293
 294   dfdy[0] = c1 * (1 - 2 * c2 * y[0] - y[1]);
 295   dfdy[1] = c1 * (1 - y[0]);
 296   dfdy[2] = 0.0;
 297
 298   dfdy[3] = 1/c1 * (-y[1]);
 299   dfdy[4] = 1/c1 * (-1 - y[0]);
 300   dfdy[5] = 1/c1;
 301
 302   dfdy[6] = c3;
 303   dfdy[7] = 0.0;
 304   dfdy[8] = -c3;
 305
 306   dfdt[0] = 0.0;
 307   dfdt[1] = 0.0;
 308   dfdt[2] = 0.0;
 309
 310   return GSL_SUCCESS;
 311 }
 312
 313 gsl_odeiv_system rhs_func_oregonator = {
 314   rhs_oregonator,
 315   jac_oregonator,
 316   3,
 317   0
 318 };
 319
 320 /* Volterra-Lotka predator-prey model
 321
 322    f0 = (a - b * y1) * y0     y0(0) = 3.0
 323    f1 = (-c + d * y0) * y1    y1(0) = 1.0
 324  */
 325
 326 int
 327 rhs_vl (double t, const double y[], double f[], void *params)
 328 {
 329   const double a = 1.0;
 330   const double b = 1.0;
 331   const double c = 1.0;
 332   const double d = 1.0;
 333
 334   f[0] = (a - b * y[1]) * y[0];
 335   f[1] = (-c + d * y[0]) * y[1];
 336
 337   return GSL_SUCCESS;
 338 }
 339
 340 int
 341 jac_vl (double t, const double y[], double *dfdy, double dfdt[],
 342             void *params)
 343 {
 344   const double a = 1.0;
 345   const double b = 1.0;
 346   const double c = 1.0;
 347   const double d = 1.0;
 348
 349   dfdy[0] = a - b * y[1];
 350   dfdy[1] = -b * y[0];
 351   dfdy[2] = d * y[1];
 352   dfdy[3] = -c + d * y[0];
 353
 354   dfdt[0] = 0.0;
 355   dfdt[1] = 0.0;
 356
 357   return GSL_SUCCESS;
 358 }
 359
 360 gsl_odeiv_system rhs_func_vl = {
 361   rhs_vl,
 362   jac_vl,
 363   2,
 364   0
 365 };
 366
 367 /* Stiff trigonometric example
 368
 369    f0 = -50 * (y0 - cos(t))    y0(0) = 0.0
 370  */
 371
 372 int
 373 rhs_stifftrig (double t, const double y[], double f[], void *params)
 374 {
 375   f[0] = -50 * (y[0] - cos(t));
 376
 377   return GSL_SUCCESS;
 378 }
 379
 380 int
 381 jac_stifftrig (double t, const double y[], double *dfdy, double dfdt[],
 382             void *params)
 383 {
 384   dfdy[0] = -50;
 385
 386   dfdt[0] = -50 * sin(t);
 387
 388   return GSL_SUCCESS;
 389 }
 390
 391 gsl_odeiv_system rhs_func_stifftrig = {
 392   rhs_stifftrig,
 393   jac_stifftrig,
 394   1,
 395   0
 396 };
 397
 398 /* E5 - a stiff badly scaled chemical problem by Enright, Hull &
 399    Lindberg (1975): Comparing numerical methods for stiff systems of
 400    ODEs. BIT, vol. 15, pp. 10-48.
 401
 402    f0 = -a * y0 - b * y0 * y2                            y0(0) = 1.76e-3
 403    f1 = a * y0 - m * c * y1 * y2                         y1(0) = 0.0
 404    f2 = a * y0 - b * y0 * y2 - m * c * y1 * y2 + c * y3  y2(0) = 0.0
 405    f3 = b * y0 * y2 - c * y3                             y3(0) = 0.0
 406  */
 407
 408 int
 409 rhs_e5 (double t, const double y[], double f[], void *params)
 410 {
 411   const double a = 7.89e-10;
 412   const double b = 1.1e7;
 413   const double c = 1.13e3;
 414   const double m = 1.0e6;
 415
 416   f[0] = -a * y[0] - b * y[0] * y[2];
 417   f[1] = a * y[0] - m * c * y[1] * y[2];
 418   f[3] = b * y[0] * y[2] - c * y[3];
 419   f[2] = f[1] - f[3];
 420
 421   return GSL_SUCCESS;
 422 }
 423
 424 int
 425 jac_e5 (double t, const double y[], double *dfdy, double dfdt[],
 426             void *params)
 427 {
 428   const double a = 7.89e-10;
 429   const double b = 1.1e7;
 430   const double c = 1.13e3;
 431   const double m = 1.0e6;
 432
 433   dfdy[0] = -a - b * y[2];
 434   dfdy[1] = 0.0;
 435   dfdy[2] = -b * y[0];
 436   dfdy[3] = 0.0;
 437
 438   dfdy[4] = a;
 439   dfdy[5] = -m * c * y[2];
 440   dfdy[6] = -m * c * y[1];
 441   dfdy[7] = 0.0;
 442
 443   dfdy[8] = a - b * y[2];
 444   dfdy[9] = -m * c * y[2];
 445   dfdy[10] = -b * y[0] - m * c * y[1];
 446   dfdy[11] = c;
 447
 448   dfdy[12] = b * y[2];
 449   dfdy[13] = 0.0;
 450   dfdy[14] = b * y[0];
 451   dfdy[15] = -c;
 452
 453   dfdt[0] = 0.0;
 454   dfdt[1] = 0.0;
 455   dfdt[2] = 0.0;
 456   dfdt[3] = 0.0;
 457
 458   return GSL_SUCCESS;
 459 }
 460
 461 gsl_odeiv_system rhs_func_e5 = {
 462   rhs_e5,
 463   jac_e5,
 464   4,
 465   0
 466 };
 467
 468 void
 469 test_odeiv_stepper (const gsl_odeiv_step_type *T, const gsl_odeiv_system *sys,
 470                     const double h, const double t, const char desc[],
 471                     const double ystart[], const double yfin[],
 472                     const double relerr)
 473 {
 474   /* tests stepper T with one fixed length step advance of system sys
 475      and compares with given values yfin
 476   */
 477
 478   double y[MAXEQ] = {0.0};
 479   double yerr[MAXEQ] = {0.0};
 480   size_t ne = sys->dimension;
 481   size_t i;
 482
 483   gsl_odeiv_step *step = gsl_odeiv_step_alloc (T, ne);
 484
 485   DBL_MEMCPY (y, ystart, MAXEQ);
 486
 487   {
 488     int s = gsl_odeiv_step_apply (step, t, h, y, yerr, 0, 0, sys);
 489     if (s != GSL_SUCCESS)
 490       {
 491         gsl_test(s, "test_odeiv_stepper: %s step_apply returned %d", desc, s);
 492       }
 493   }
 494
 495   for (i = 0; i < ne; i++)
 496     {
 497       gsl_test_rel (y[i], yfin[i], relerr,
 498                     "%s %s step(%d)",
 499                     gsl_odeiv_step_name (step), desc,i);
 500     }
 501
 502   gsl_odeiv_step_free (step);
 503 }
 504
 505 void
 506 test_stepper (const gsl_odeiv_step_type *T)
 507 {
 508   /* Tests stepper T with a step of selected systems */
 509
 510   double y[MAXEQ] = {0.0};
 511   double yfin[MAXEQ] = {0.0};
 512
 513   /* Step length */
 514   double h;
 515
 516   /* Required tolerance */
 517   double err_target;
 518
 519   /* linear */
 520   h = 1e-1;
 521   err_target = 1e-10;
 522   y[0] = 0.58;
 523   yfin[0] = y[0] + 2 * h;
 524   test_odeiv_stepper (T, &rhs_func_lin, h, 0.0, "linear",
 525                       y, yfin, err_target);
 526
 527   /* exponential */
 528   h = 1e-4;
 529   err_target = 1e-8;
 530   y[0] = exp(2.7);
 531   yfin[0] = exp(2.7 + h);
 532   test_odeiv_stepper (T, &rhs_func_exp, h, 2.7, "exponential",
 533                       y, yfin, err_target);
 534   /* cosine-sine */
 535   h = 1e-3;
 536   err_target = 1e-6;
 537   y[0] = cos(1.2);
 538   y[1] = sin(1.2);
 539   yfin[0] = cos(1.2 + h);
 540   yfin[1] = sin(1.2 + h);
 541   test_odeiv_stepper (T, &rhs_func_sin, h, 1.2, "cosine-sine",
 542                       y, yfin, err_target);
 543
 544   /* classic stiff */
 545   h = 1e-7;
 546   err_target = 1e-4;
 547   y[0] = 1.0;
 548   y[1] = 0.0;
 549
 550   {
 551     const double e1 = exp (-h);
 552     const double e2 = exp (-1000.0 * h);
 553     yfin[0] = 2.0 * e1 - e2;
 554     yfin[1] = -e1 + e2;
 555   }
 556
 557   test_odeiv_stepper (T, &rhs_func_stiff, h, 0.0, "classic_stiff",
 558                       y, yfin, err_target);
 559 }
 560
 561 void
 562 test_evolve_system (const gsl_odeiv_step_type * T,
 563                     const gsl_odeiv_system * sys,
 564                     double t0, double t1, double hstart,
 565                     double y[], double yfin[],
 566                     double err_target, const char *desc)
 567 {
 568   /* Tests system sys with stepper T. Step length is controlled by
 569      error estimation from the stepper.
 570   */
 571
 572   int steps = 0;
 573   size_t i;
 574
 575   double t = t0;
 576   double h = hstart;
 577
 578   /* Tolerance factor in testing errors */
 579   const double factor = 10;
 580
 581   gsl_odeiv_step * step = gsl_odeiv_step_alloc (T, sys->dimension);
 582
 583   gsl_odeiv_control *c =
 584     gsl_odeiv_control_standard_new (err_target, err_target, 1.0, 0.0);
 585
 586   gsl_odeiv_evolve *e = gsl_odeiv_evolve_alloc (sys->dimension);
 587
 588   while (t < t1)
 589     {
 590       int s = gsl_odeiv_evolve_apply (e, c, step, sys, &t, t1, &h, y);
 591
 592       if (s != GSL_SUCCESS && sys != &rhs_func_xsin)
 593         {
 594           gsl_test(s, "%s evolve_apply returned %d",
 595                    gsl_odeiv_step_name (step), s);
 596           break;
 597         }
 598
 599       if (steps > 100000)
 600         {
 601           gsl_test(GSL_EFAILED,
 602                    "%s evolve_apply reached maxiter",
 603                    gsl_odeiv_step_name (step));
 604           break;
 605         }
 606
 607       steps++;
 608     }
 609
 610   /* err_target is target error of one step. Test if stepper has made
 611      larger error than (tolerance factor times) the number of steps
 612      times the err_target */
 613
 614   for (i = 0; i < sys->dimension; i++)
 615     {
 616       gsl_test_abs (y[i], yfin[i], factor * e->count * err_target,
 617                     "%s %s evolve(%d)",
 618                     gsl_odeiv_step_name (step), desc, i);
 619     }
 620
 621   gsl_odeiv_evolve_free (e);
 622   gsl_odeiv_control_free (c);
 623   gsl_odeiv_step_free (step);
 624 }
 625
 626 int
 627 sys_driver (const gsl_odeiv_step_type * T,
 628             const gsl_odeiv_system * sys,
 629             double t0, double t1, double hstart,
 630             double y[], double epsabs, double epsrel,
 631             const char desc[])
 632 {
 633   /* This function evolves a system sys with stepper T from t0 to t1.
 634      Step length is varied via error control with possibly different
 635      absolute and relative error tolerances.
 636   */
 637
 638   int s = 0;
 639   int steps = 0;
 640
 641   double t = t0;
 642   double h = hstart;
 643
 644   gsl_odeiv_step * step = gsl_odeiv_step_alloc (T, sys->dimension);
 645
 646   gsl_odeiv_control *c =
 647     gsl_odeiv_control_standard_new (epsabs, epsrel, 1.0, 0.0);
 648   gsl_odeiv_evolve *e = gsl_odeiv_evolve_alloc (sys->dimension);
 649
 650   while (t < t1)
 651     {
 652       s = gsl_odeiv_evolve_apply (e, c, step, sys, &t, t1, &h, y);
 653
 654       if (s != GSL_SUCCESS)
 655         {
 656           gsl_test(s, "sys_driver: %s evolve_apply returned %d",
 657                    gsl_odeiv_step_name (step), s);
 658           break;
 659         }
 660
 661       if (steps > 1e7)
 662         {
 663           gsl_test(GSL_EMAXITER,
 664                    "sys_driver: %s evolve_apply reached maxiter at t=%g",
 665                    gsl_odeiv_step_name (step), t);
 666           s = GSL_EMAXITER;
 667           break;
 668         }
 669
 670       steps++;
 671     }
 672
 673   gsl_test(s, "%s %s [%g,%g], %d steps completed",
 674            gsl_odeiv_step_name (step), desc, t0, t1, steps);
 675
 676   gsl_odeiv_evolve_free (e);
 677   gsl_odeiv_control_free (c);
 678   gsl_odeiv_step_free (step);
 679
 680   return s;
 681 }
 682
 683 void
 684 test_compare_vanderpol (void)
 685 {
 686   /* Compares output of van Der Pol oscillator with several steppers */
 687
 688   /* system dimension */
 689   const size_t sd = 2;
 690
 691   const gsl_odeiv_step_type *steppers[20];
 692   const gsl_odeiv_step_type **T;
 693
 694   /* Required error tolerance for each stepper. */
 695   double err_target[20];
 696
 697   /* number of ODE solvers */
 698   const size_t ns = 11;
 699
 700   /* initial values for each ode-solver */
 701   double y[11][2];
 702   double *yp = &y[0][0];
 703
 704   size_t i, j, k;
 705   int status = 0;
 706
 707   /* Parameters for the problem and stepper  */
 708   const double start = 0.0;
 709   const double end = 100.0;
 710   const double epsabs = 1e-8;
 711   const double epsrel = 1e-8;
 712   const double initstepsize = 1e-5;
 713
 714   /* Initialize */
 715
 716   steppers[0] = gsl_odeiv_step_rk2;
 717   err_target[0] = 1e-6;
 718   steppers[1] = gsl_odeiv_step_rk4;
 719   err_target[1] = 1e-6;
 720   steppers[2] = gsl_odeiv_step_rkf45;
 721   err_target[2] = 1e-6;
 722   steppers[3] = gsl_odeiv_step_rkck;
 723   err_target[3] = 1e-6;
 724   steppers[4] = gsl_odeiv_step_rk8pd;
 725   err_target[4] = 1e-6;
 726   steppers[5] = gsl_odeiv_step_rk2imp;
 727   err_target[5] = 1e-5;
 728   steppers[6] = gsl_odeiv_step_rk2simp;
 729   err_target[6] = 1e-5;
 730   steppers[7] = gsl_odeiv_step_rk4imp;
 731   err_target[7] = 1e-6;
 732   steppers[8] = gsl_odeiv_step_bsimp;
 733   err_target[8] = 1e-7;
 734   steppers[9] = gsl_odeiv_step_gear1;
 735   err_target[9] = 1e-2;
 736   steppers[10] = gsl_odeiv_step_gear2;
 737   err_target[10] = 1e-6;
 738   steppers[11] = 0;
 739
 740   T = steppers;
 741
 742   for (i = 0; i < ns; i++)
 743     {
 744       y[i][0] = 1.0;
 745       y[i][1] = 0.0;
 746     }
 747
 748   /* Call each solver for the problem */
 749
 750   i = 0;
 751   while (*T != 0)
 752     {
 753       {
 754         int s = sys_driver (*T, &rhs_func_vanderpol,
 755                             start, end, initstepsize, &yp[i],
 756                             epsabs, epsrel, "vanderpol");
 757         if (s != GSL_SUCCESS)
 758           {
 759             status++;
 760           }
 761       }
 762
 763       T++;
 764       i += sd;
 765     }
 766
 767   if (status != GSL_SUCCESS)
 768     {
 769       return;
 770     }
 771
 772   /* Compare results */
 773
 774   T = steppers;
 775
 776   for (i = 0; i < ns; i++)
 777     for (j = i+1; j < ns; j++)
 778       for (k = 0; k < sd; k++)
 779         {
 780           const double val1 = yp[sd * i + k];
 781           const double val2 = yp[sd * j + k];
 782           gsl_test_abs (val1, val2,
 783                         ( GSL_MAX(err_target[i], err_target[j]) ),
 784                         "%s/%s vanderpol",
 785                         T[i]->name, T[j]->name);
 786         }
 787
 788 }
 789
 790 void
 791 test_compare_oregonator (void)
 792 {
 793   /* Compares output of the Oregonator with several steppers */
 794
 795   /* system dimension */
 796   const size_t sd = 3;
 797
 798   const gsl_odeiv_step_type *steppers[20];
 799   const gsl_odeiv_step_type **T;
 800
 801   /* Required error tolerance for each stepper. */
 802   double err_target[20];
 803
 804   /* number of ODE solvers */
 805   const size_t ns = 2;
 806
 807   /* initial values for each ode-solver */
 808   double y[2][3];
 809   double *yp = &y[0][0];
 810
 811   size_t i, j, k;
 812   int status = 0;
 813
 814   /* Parameters for the problem and stepper  */
 815   const double start = 0.0;
 816   const double end = 360.0;
 817   const double epsabs = 1e-8;
 818   const double epsrel = 1e-8;
 819   const double initstepsize = 1e-5;
 820
 821   /* Initialize */
 822
 823   steppers[0] = gsl_odeiv_step_rk2simp;
 824   err_target[0] = 1e-6;
 825   steppers[1] = gsl_odeiv_step_bsimp;
 826   err_target[1] = 1e-6;
 827   steppers[2] = 0;
 828
 829   T = steppers;
 830
 831   for (i = 0; i < ns; i++)
 832     {
 833       y[i][0] = 1.0;
 834       y[i][1] = 2.0;
 835       y[i][2] = 3.0;
 836     }
 837
 838   /* Call each solver for the problem */
 839
 840   i = 0;
 841   while (*T != 0)
 842     {
 843       {
 844         int s = sys_driver (*T, &rhs_func_oregonator,
 845                             start, end, initstepsize, &yp[i],
 846                             epsabs, epsrel, "oregonator");
 847
 848         if (s != GSL_SUCCESS)
 849           {
 850             status++;
 851           }
 852       }
 853
 854       T++;
 855       i += sd;
 856     }
 857
 858   if (status != GSL_SUCCESS)
 859     {
 860       return;
 861     }
 862
 863   /* Compare results */
 864
 865   T = steppers;
 866
 867   for (i = 0; i < ns; i++)
 868     for (j = i+1; j < ns; j++)
 869       for (k = 0; k < sd; k++)
 870         {
 871           const double val1 = yp[sd * i + k];
 872           const double val2 = yp[sd * j + k];
 873           gsl_test_rel (val1, val2,
 874                         ( GSL_MAX(err_target[i], err_target[j]) ),
 875                         "%s/%s oregonator",
 876                         T[i]->name, T[j]->name);
 877         }
 878
 879 }
 880
 881 void
 882 test_evolve_linear (const gsl_odeiv_step_type * T, double h, double err)
 883 {
 884   double y[1];
 885   double yfin[1];
 886
 887   y[0] = 1.0;
 888   yfin[0] = 9.0;
 889   test_evolve_system (T, &rhs_func_lin, 0.0, 4.0, h, y, yfin, err,
 890                       "linear[0,4]");
 891 }
 892
 893 void
 894 test_evolve_exp (const gsl_odeiv_step_type * T, double h, double err)
 895 {
 896   double y[1];
 897   double yfin[1];
 898
 899   y[0] = 1.0;
 900   yfin[0] = exp (2.0);
 901   test_evolve_system (T, &rhs_func_exp, 0.0, 2.0, h, y, yfin, err,
 902                       "exp[0,2]");
 903 }
 904
 905 void
 906 test_evolve_sin (const gsl_odeiv_step_type * T, double h, double err)
 907 {
 908   double y[2];
 909   double yfin[2];
 910
 911   y[0] = 1.0;
 912   y[1] = 0.0;
 913   yfin[0] = cos (2.0);
 914   yfin[1] = sin (2.0);
 915   test_evolve_system (T, &rhs_func_sin, 0.0, 2.0, h, y, yfin, err,
 916                       "sine[0,2]");
 917 }
 918
 919 void
 920 test_evolve_xsin (const gsl_odeiv_step_type * T, double h, double err)
 921 {
 922   double y[2];
 923   double yfin[2];
 924
 925   y[0] = 1.0;
 926   y[1] = 0.0;
 927   yfin[0] = cos (2.0);
 928   yfin[1] = sin (2.0);
 929   test_evolve_system (T, &rhs_func_xsin, 0.0, 2.0, h, y, yfin, err,
 930                       "sine[0,2] w/errors");
 931 }
 932
 933
 934 void
 935 test_evolve_stiff1 (const gsl_odeiv_step_type * T, double h, double err)
 936 {
 937   double y[2];
 938   double yfin[2];
 939
 940   y[0] = 1.0;
 941   y[1] = 0.0;
 942   {
 943     double arg = 1.0;
 944     double e1 = exp (-arg);
 945     double e2 = exp (-1000.0 * arg);
 946     yfin[0] = 2.0 * e1 - e2;
 947     yfin[1] = -e1 + e2;
 948   }
 949   test_evolve_system (T, &rhs_func_stiff, 0.0, 1.0, h, y, yfin, err,
 950                       "stiff[0,1]");
 951 }
 952
 953 void
 954 test_evolve_stiff5 (const gsl_odeiv_step_type * T, double h, double err)
 955 {
 956   double y[2];
 957   double yfin[2];
 958
 959   y[0] = 1.0;
 960   y[1] = 0.0;
 961   {
 962     double arg = 5.0;
 963     double e1 = exp (-arg);
 964     double e2 = exp (-1000.0 * arg);
 965     yfin[0] = 2.0 * e1 - e2;
 966     yfin[1] = -e1 + e2;
 967   }
 968   test_evolve_system (T, &rhs_func_stiff, 0.0, 5.0, h, y, yfin, err,
 969                       "stiff[0,5]");
 970 }
 971
 972 /* Test cases from Frank Reininghaus <frank78ac@googlemail.com> */
 973
 974 int rhs_stepfn (double t, const double * y, double * dydt, void * params) {
 975   if (t >= 1.0)
 976     dydt [0] = 1;
 977   else
 978     dydt [0] = 0;
 979
 980   return GSL_SUCCESS;
 981 }
 982
 983 void test_stepfn (void) {
 984   /* infinite loop for epsabs = 1e-18, but not for 1e-17 */
 985   double epsabs = 1e-18;
 986   double epsrel = 1e-6;
 987
 988   const gsl_odeiv_step_type * T = gsl_odeiv_step_rk2;
 989   gsl_odeiv_step * s = gsl_odeiv_step_alloc (T, 1);
 990   gsl_odeiv_control * c = gsl_odeiv_control_y_new (epsabs, epsrel);
 991   gsl_odeiv_evolve * e = gsl_odeiv_evolve_alloc (1);
 992   gsl_odeiv_system sys = {rhs_stepfn, 0, 1, 0};
 993
 994   double t = 0.0;
 995   double h = 1e-6;
 996   double y = 0.0;
 997   int i = 0;
 998   int status;
 999
1000   while (t < 2.0 && i < 1000000) {
1001     status = gsl_odeiv_evolve_apply (e, c, s, &sys, &t, 2, &h, &y);
1002 #ifdef DEBUG
1003     printf("i=%d status=%d t=%g h=%g y=%g\n", i, status, t, h, y);
1004 #endif
1005     if (status != GSL_SUCCESS)
1006       break;
1007
1008     i++;
1009   }
1010
1011   gsl_test_abs(t, 2.0, 1e-16, "evolve step function, t (stepfn/rk2)");
1012   gsl_test_rel(y, 1.0, epsrel, "evolve step function, y (stepfn/rk2)");
1013
1014   gsl_odeiv_evolve_free (e);
1015   gsl_odeiv_control_free (c);
1016   gsl_odeiv_step_free (s);
1017 }
1018
1019 int rhs_stepfn2 (double t, const double * y, double * dydt, void * params) {
1020   if (t >= 0.0)
1021     dydt [0] = 1e300;
1022   else
1023     dydt [0] = 0;
1024
1025   return GSL_SUCCESS;
1026 }
1027
1028 void test_stepfn2 (void) {
1029   /* infinite loop for epsabs = 1e-25, but not for 1e-24 */
1030   double epsabs = 1e-25;
1031   double epsrel = 1e-6;
1032
1033   const gsl_odeiv_step_type * T = gsl_odeiv_step_rk2;
1034   gsl_odeiv_step * s = gsl_odeiv_step_alloc (T, 1);
1035   gsl_odeiv_control * c = gsl_odeiv_control_y_new (epsabs, epsrel);
1036   gsl_odeiv_evolve * e = gsl_odeiv_evolve_alloc (1);
1037   gsl_odeiv_system sys = {rhs_stepfn2, 0, 1, 0};
1038
1039   double t = -1.0;
1040   double h = 1e-6;
1041   double y = 0.0;
1042
1043   int i = 0;
1044   int status;
1045
1046   while (t < 1.0 && i < 10000) {
1047     status = gsl_odeiv_evolve_apply (e, c, s, &sys, &t, 1.0, &h, &y);
1048 #ifdef DEBUG
1049     printf("i=%d status=%d t=%g h=%g y=%g\n", i, status, t, h, y);
1050 #endif
1051     if (status != GSL_SUCCESS)
1052       break;
1053
1054     i++;
1055   }
1056
1057   gsl_test_abs(t, 1.0, 1e-16, "evolve big step function, t (stepfn2/rk2)");
1058   gsl_test_rel(y, 1e300, epsrel, "evolve big step function, y (stepfn2/rk2)");
1059
1060   gsl_odeiv_evolve_free (e);
1061   gsl_odeiv_control_free (c);
1062   gsl_odeiv_step_free (s);
1063 }
1064
1065 int rhs_stepfn3 (double t, const double * y, double * dydt, void * params) {
1066
1067   static int calls = 0;
1068
1069   if (t >= 0.0)
1070     dydt [0] = 1e300;
1071   else
1072     dydt [0] = 0;
1073
1074   calls++;
1075
1076   return (calls < 100000) ? GSL_SUCCESS : -999;
1077 }
1078
1079
1080 void test_stepfn3 (void) {
1081   /* infinite loop for epsabs = 1e-26, but not for 1e-25 */
1082   double epsabs = 1e-26;
1083   double epsrel = 1e-6;
1084
1085   const gsl_odeiv_step_type * T = gsl_odeiv_step_rkf45;
1086   gsl_odeiv_step * s = gsl_odeiv_step_alloc (T, 1);
1087   gsl_odeiv_control * c = gsl_odeiv_control_y_new (epsabs, epsrel);
1088   gsl_odeiv_evolve * e = gsl_odeiv_evolve_alloc (1);
1089   gsl_odeiv_system sys = {rhs_stepfn3, 0, 1, 0};
1090
1091   double t = -1.0;
1092   double h = 1e-6;
1093   double y = 0.0;
1094
1095   int i = 0;
1096   int status;
1097
1098   while (t < 1.0 && i < 10000) {
1099     status = gsl_odeiv_evolve_apply (e, c, s, &sys, &t, 1.0, &h, &y);
1100 #ifdef DEBUG
1101     printf("i=%d status=%d t=%g h=%g y=%g\n", i, status, t, h, y);
1102 #endif
1103     if (status != GSL_SUCCESS)
1104       break;
1105
1106     i++;
1107   }
1108
1109   gsl_test_abs(t, 1.0, 1e-16, "evolve big step function, t (stepfn3/rkf45)");
1110   gsl_test_rel(y, 1e300, epsrel, "evolve big step function, y (stepfn3/rkf45)");
1111
1112   gsl_odeiv_evolve_free (e);
1113   gsl_odeiv_control_free (c);
1114   gsl_odeiv_step_free (s);
1115 }
1116
1117 int rhs_cos (double t, const double * y, double * dydt, void * params) {
1118   dydt [0] = cos (t);
1119   return GSL_SUCCESS;
1120 }
1121
1122 int jac_cos (double t, const double y[], double *dfdy, double dfdt[],
1123              void *params)
1124 {
1125   dfdy[0] = 0.0;
1126   dfdt[0] = -sin(t);
1127
1128   return GSL_SUCCESS;
1129 }
1130
1131 /* Test evolution in negative direction */
1132
1133 void
1134 test_evolve_negative_h (const gsl_odeiv_step_type * T, double h, double err)
1135 {
1136   /* Tolerance factor in testing errors */
1137   const double factor = 10;
1138
1139   gsl_odeiv_step * step = gsl_odeiv_step_alloc (T, 1);
1140   gsl_odeiv_control * c = gsl_odeiv_control_standard_new (err, err, 1.0, 0.0);
1141   gsl_odeiv_evolve * e = gsl_odeiv_evolve_alloc (1);
1142   gsl_odeiv_system sys = {rhs_cos, jac_cos, 1, 0};
1143
1144   double t = 0;
1145   double t1 = -4.0;
1146
1147   double y = 0.0;
1148   double yfin = sin (t1);
1149
1150   /* Make initial h negative */
1151   h = -fabs(h);
1152
1153   while (t > t1) {
1154     int status = gsl_odeiv_evolve_apply (e, c, step, &sys, &t, t1, &h, &y);
1155
1156     if (status != GSL_SUCCESS)
1157       {
1158         gsl_test(status, "%s evolve_apply returned %d for negative h",
1159                  gsl_odeiv_step_name (step), status);
1160         break;
1161       }
1162   }
1163
1164   gsl_test_abs (y, yfin, factor * e->count * err,
1165                 "evolution with negative h (using %s)",
1166                 gsl_odeiv_step_name (step));
1167
1168   gsl_odeiv_evolve_free (e);
1169   gsl_odeiv_control_free (c);
1170   gsl_odeiv_step_free (step);
1171 }
1172
1173 int
1174 main (void)
1175 {
1176   int i;
1177
1178   struct ptype
1179   {
1180     const gsl_odeiv_step_type *type;
1181     double h;
1182   }
1183   p[20];
1184
1185   p[0].type = gsl_odeiv_step_rk2;
1186   p[0].h = 1.0e-3;
1187   p[1].type = gsl_odeiv_step_rk4;
1188   p[1].h = 1.0e-3;
1189   p[2].type = gsl_odeiv_step_rkf45;
1190   p[2].h = 1.0e-3;
1191   p[3].type = gsl_odeiv_step_rkck;
1192   p[3].h = 1.0e-3;
1193   p[4].type = gsl_odeiv_step_rk8pd;
1194   p[4].h = 1.0e-3;
1195   p[5].type = gsl_odeiv_step_rk2imp;
1196   p[5].h = 1.0e-3;
1197   p[6].type = gsl_odeiv_step_rk2simp;
1198   p[6].h = 1.0e-3;
1199   p[7].type = gsl_odeiv_step_rk4imp;
1200   p[7].h = 1.0e-3;
1201   p[8].type = gsl_odeiv_step_bsimp;
1202   p[8].h = 1.0e-3;
1203   p[9].type = gsl_odeiv_step_gear1;
1204   p[9].h = 1.0e-3;
1205   p[10].type = gsl_odeiv_step_gear2;
1206   p[10].h = 1.0e-3;
1207   p[11].type = 0;
1208
1209   gsl_ieee_env_setup ();
1210
1211   for (i = 0; p[i].type != 0; i++)
1212     {
1213       test_stepper(p[i].type);
1214     }
1215
1216   for (i = 0; p[i].type != 0; i++)
1217     {
1218       test_evolve_linear (p[i].type, p[i].h, 1e-10);
1219       test_evolve_exp (p[i].type, p[i].h, 1e-6);
1220       test_evolve_sin (p[i].type, p[i].h, 1e-8);
1221       test_evolve_xsin (p[i].type, p[i].h, 1e-8);
1222       test_evolve_stiff1 (p[i].type, p[i].h, 1e-7);
1223       test_evolve_stiff5 (p[i].type, p[i].h, 1e-7);
1224       test_evolve_negative_h (p[i].type, p[i].h, 1e-7);
1225     }
1226
1227   test_compare_vanderpol();
1228   test_compare_oregonator();
1229
1230   test_stepfn();
1231   test_stepfn2();
1232   test_stepfn3();
1233
1234   exit (gsl_test_summary ());
1235 }