/* ode-initval/test.c
*
* Copyright (C) 2009, 2010 Tuomo Keskitalo
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 3 of the License, or (at
* your option) any later version.
*
* This program is distributed in the hope that it will be useful, but
* WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
* General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program; if not, write to the Free Software
* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
*/
/* Some functions and tests based on test.c by G. Jungman. */
#include <config.h>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <gsl/gsl_test.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_ieee_utils.h>
#include <gsl/gsl_odeiv2.h>
#include "odeiv_util.h"
/* Maximum number of ODE equations */
#define MAXEQ 15
/* Maximum number of ODE solvers */
#define MAXNS 20
/* Track number of function and jacobian evaluations in tests
with global variables
*/
int nfe;
int nje;
/**********************************************************/
/* ODE test system definitions */
/**********************************************************/
/* RHS for f=2. Solution y = 2 * t + t0 */
int
rhs_linear (double t, const double y[], double f[], void *params)
{
extern int nfe;
nfe += 1;
f[0] = 2.0;
return GSL_SUCCESS;
}
int
jac_linear (double t, const double y[], double *dfdy, double dfdt[],
void *params)
{
extern int nje;
nje += 1;
dfdy[0] = 0.0;
dfdt[0] = 0.0;
return GSL_SUCCESS;
}
gsl_odeiv2_system rhs_func_lin = {
rhs_linear,
jac_linear,
1,
0
};
/* RHS for f=y. Equals y=exp(t) with initial value y(0)=1.0 */
int
rhs_exp (double t, const double y[], double f[], void *params)
{
extern int nfe;
nfe += 1;
f[0] = y[0];
return GSL_SUCCESS;
}
int
jac_exp (double t, const double y[], double *dfdy, double dfdt[],
void *params)
{
extern int nje;
nje += 1;
dfdy[0] = y[0];
dfdt[0] = 0.0;
return GSL_SUCCESS;
}
gsl_odeiv2_system rhs_func_exp = {
rhs_exp,
jac_exp,
1,
0
};
int
rhs_sin (double t, const double y[], double f[], void *params)
{
extern int nfe;
nfe += 1;
f[0] = -y[1];
f[1] = y[0];
return GSL_SUCCESS;
}
int
jac_sin (double t, const double y[], double *dfdy, double dfdt[],
void *params)
{
extern int nje;
nje += 1;
dfdy[0] = 0.0;
dfdy[1] = -1.0;
dfdy[2] = 1.0;
dfdy[3] = 0.0;
dfdt[0] = 0.0;
dfdt[1] = 0.0;
return GSL_SUCCESS;
}
gsl_odeiv2_system rhs_func_sin = {
rhs_sin,
jac_sin,
2,
0
};
/* Sine/cosine with random failures */
static int rhs_xsin_reset = 0;
static int jac_xsin_reset = 0;
int
rhs_xsin (double t, const double y[], double f[], void *params)
{
static int n = 0, m = 0;
extern int nfe;
nfe += 1;
if (rhs_xsin_reset)
{
rhs_xsin_reset = 0;
n = 0;
m = 1;
}
n++;
if (n >= m)
{
m = n * 1.3;
return GSL_EFAILED;
}
if (n > 40 && n < 65)
{
f[0] = GSL_NAN;
f[1] = GSL_NAN;
return GSL_EFAILED;
}
f[0] = -y[1];
f[1] = y[0];
return GSL_SUCCESS;
}
int
jac_xsin (double t, const double y[], double *dfdy, double dfdt[],
void *params)
{
static int n = 0;
extern int nje;
nje += 1;
if (jac_xsin_reset)
{
jac_xsin_reset = 0;
n = 0;
}
n++;
if (n > 50 && n < 55)
{
dfdy[0] = GSL_NAN;
dfdy[1] = GSL_NAN;
dfdy[2] = GSL_NAN;
dfdy[3] = GSL_NAN;
dfdt[0] = GSL_NAN;
dfdt[1] = GSL_NAN;
return GSL_EFAILED;
}
dfdy[0] = 0.0;
dfdy[1] = -1.0;
dfdy[2] = 1.0;
dfdy[3] = 0.0;
dfdt[0] = 0.0;
dfdt[1] = 0.0;
return GSL_SUCCESS;
}
gsl_odeiv2_system rhs_func_xsin = {
rhs_xsin,
jac_xsin,
2,
0
};
/* RHS for classic stiff example
dy0 / dt = 998 * y0 + 1998 * y1 y0(0) = 1.0
dy1 / dt = -999 * y0 - 1999 * y1 y1(0) = 0.0
solution is
y0 = 2 * exp(-t) - exp(-1000 * t)
y1 = - exp(-t) + exp(-1000 * t)
*/
int
rhs_stiff (double t, const double y[], double f[], void *params)
{
extern int nfe;
nfe += 1;
f[0] = 998.0 * y[0] + 1998.0 * y[1];
f[1] = -999.0 * y[0] - 1999.0 * y[1];
return GSL_SUCCESS;
}
int
jac_stiff (double t, const double y[], double *dfdy, double dfdt[],
void *params)
{
extern int nje;
nje += 1;
dfdy[0] = 998.0;
dfdy[1] = 1998.0;
dfdy[2] = -999.0;
dfdy[3] = -1999.0;
dfdt[0] = 0.0;
dfdt[1] = 0.0;
return GSL_SUCCESS;
}
gsl_odeiv2_system rhs_func_stiff = {
rhs_stiff,
jac_stiff,
2,
0
};
/* Cosine function */
int
rhs_cos (double t, const double *y, double *dydt, void *params)
{
dydt[0] = cos (t);
return GSL_SUCCESS;
}
int
jac_cos (double t, const double y[], double *dfdy, double dfdt[],
void *params)
{
dfdy[0] = 0.0;
dfdt[0] = -sin (t);
return GSL_SUCCESS;
}
gsl_odeiv2_system rhs_func_cos = {
rhs_cos,
jac_cos,
1,
0
};
/* Broken problem for testing numerical problems in user function that
leads to decrease of step size in gsl_odeiv2_evolve below machine
precision.
*/
int
rhs_broken (double t, const double y[], double f[], void *params)
{
if (t < 10.0)
{
f[0] = 1.0;
}
else
{
f[0] = GSL_NAN;
return 123;
}
return GSL_SUCCESS;
}
int
jac_broken (double t, const double y[], double *dfdy, double dfdt[],
void *params)
{
if (t < 10.0)
{
dfdy[0] = 0.0;
dfdt[0] = 0.0;
}
else
{
dfdy[0] = GSL_NAN;
dfdt[0] = GSL_NAN;
return 123;
}
return GSL_SUCCESS;
}
gsl_odeiv2_system rhs_func_broken = {
rhs_broken,
jac_broken,
1,
0
};
/* Immediate user break (at t > 1.5) test sine system */
int
rhs_sin_ub (double t, const double y[], double f[], void *params)
{
extern int nfe;
nfe += 1;
f[0] = -y[1];
f[1] = y[0];
if (t > 1.5)
{
return GSL_EBADFUNC;
}
return GSL_SUCCESS;
}
int
jac_sin_ub (double t, const double y[], double *dfdy, double dfdt[],
void *params)
{
extern int nje;
nje += 1;
dfdy[0] = 0.0;
dfdy[1] = -1.0;
dfdy[2] = 1.0;
dfdy[3] = 0.0;
dfdt[0] = 0.0;
dfdt[1] = 0.0;
if (t > 1.5)
{
return GSL_EBADFUNC;
}
return GSL_SUCCESS;
}
gsl_odeiv2_system rhs_func_sin_ub = {
rhs_sin_ub,
jac_sin_ub,
2,
0
};
/* Badly scaled random function */
int
rhs_br (double t, const double *y, double *dydt, void *params)
{
dydt[0] = (rand () - RAND_MAX / 2) * 2e100;
return GSL_SUCCESS;
}
int
jac_br (double t, const double y[], double *dfdy, double dfdt[], void *params)
{
dfdy[0] = (rand () - RAND_MAX / 2) * 2e100;
dfdt[0] = (rand () - RAND_MAX / 2) * 2e100;
return GSL_SUCCESS;
}
gsl_odeiv2_system rhs_func_br = {
rhs_br,
jac_br,
1,
0
};
/* stepfn and stepfn2 based on testcases from Frank Reininghaus
<frank78ac@googlemail.com>
*/
/* Derivative change at t=0, small derivative */
int
rhs_stepfn (double t, const double *y, double *dydt, void *params)
{
if (t >= 1.0)
dydt[0] = 1;
else
dydt[0] = 0;
return GSL_SUCCESS;
}
int
jac_stepfn (double t, const double y[], double *dfdy, double dfdt[],
void *params)
{
dfdy[0] = 0.0;
dfdt[0] = 0.0;
return GSL_SUCCESS;
}
gsl_odeiv2_system rhs_func_stepfn = {
rhs_stepfn,
jac_stepfn,
1,
0
};
/* Derivative change at t=0, large derivative */
int
rhs_stepfn2 (double t, const double *y, double *dydt, void *params)
{
if (t >= 0.0)
dydt[0] = 1e300;
else
dydt[0] = 0;
return GSL_SUCCESS;
}
gsl_odeiv2_system rhs_func_stepfn2 = {
rhs_stepfn2,
jac_stepfn,
1,
0
};
/* Volterra-Lotka predator-prey model
f0 = (a - b * y1) * y0 y0(0) = 2.725
f1 = (-c + d * y0) * y1 y1(0) = 1.0
*/
int
rhs_vl (double t, const double y[], double f[], void *params)
{
const double a = -1.0;
const double b = -1.0;
const double c = -2.0;
const double d = -1.0;
extern int nfe;
nfe += 1;
f[0] = (a - b * y[1]) * y[0];
f[1] = (-c + d * y[0]) * y[1];
return GSL_SUCCESS;
}
int
jac_vl (double t, const double y[], double *dfdy, double dfdt[], void *params)
{
const double a = -1.0;
const double b = -1.0;
const double c = -2.0;
const double d = -1.0;
extern int nje;
nje += 1;
dfdy[0] = a - b * y[1];
dfdy[1] = -b * y[0];
dfdy[2] = d * y[1];
dfdy[3] = -c + d * y[0];
dfdt[0] = 0.0;
dfdt[1] = 0.0;
return GSL_SUCCESS;
}
gsl_odeiv2_system rhs_func_vl = {
rhs_vl,
jac_vl,
2,
0
};
/* van Der Pol oscillator
f0 = y1 y0(0) = 1.0
f1 = -y0 + mu * y1 * (1 - y0^2) y1(0) = 0.0
*/
int
rhs_vanderpol (double t, const double y[], double f[], void *params)
{
const double mu = 10.0;
extern int nfe;
nfe += 1;
f[0] = y[1];
f[1] = -y[0] + mu * y[1] * (1.0 - y[0] * y[0]);
return GSL_SUCCESS;
}
int
jac_vanderpol (double t, const double y[], double *dfdy, double dfdt[],
void *params)
{
const double mu = 10.0;
extern int nje;
nje += 1;
dfdy[0] = 0.0;
dfdy[1] = 1.0;
dfdy[2] = -2.0 * mu * y[0] * y[1] - 1.0;
dfdy[3] = mu * (1.0 - y[0] * y[0]);
dfdt[0] = 0.0;
dfdt[1] = 0.0;
return GSL_SUCCESS;
}
gsl_odeiv2_system rhs_func_vanderpol = {
rhs_vanderpol,
jac_vanderpol,
2,
0
};
/* Stiff trigonometric example
f0 = -50 * (y0 - cos(t)) y0(0) = 0.0
*/
int
rhs_stifftrig (double t, const double y[], double f[], void *params)
{
extern int nfe;
nfe += 1;
f[0] = -50 * (y[0] - cos (t));
return GSL_SUCCESS;
}
int
jac_stifftrig (double t, const double y[], double *dfdy, double dfdt[],
void *params)
{
extern int nje;
nje += 1;
dfdy[0] = -50;
dfdt[0] = -50 * sin (t);
return GSL_SUCCESS;
}
gsl_odeiv2_system rhs_func_stifftrig = {
rhs_stifftrig,
jac_stifftrig,
1,
0
};
/* The Oregonator - chemical Belusov-Zhabotinskii reaction
y0(0) = 1.0, y1(0) = 2.0, y2(0) = 3.0
*/
int
rhs_oregonator (double t, const double y[], double f[], void *params)
{
const double c1 = 77.27;
const double c2 = 8.375e-6;
const double c3 = 0.161;
extern int nfe;
nfe += 1;
f[0] = c1 * (y[1] + y[0] * (1 - c2 * y[0] - y[1]));
f[1] = 1 / c1 * (y[2] - y[1] * (1 + y[0]));
f[2] = c3 * (y[0] - y[2]);
return GSL_SUCCESS;
}
int
jac_oregonator (double t, const double y[], double *dfdy, double dfdt[],
void *params)
{
const double c1 = 77.27;
const double c2 = 8.375e-6;
const double c3 = 0.161;
extern int nje;
nje += 1;
dfdy[0] = c1 * (1 - 2 * c2 * y[0] - y[1]);
dfdy[1] = c1 * (1 - y[0]);
dfdy[2] = 0.0;
dfdy[3] = 1 / c1 * (-y[1]);
dfdy[4] = 1 / c1 * (-1 - y[0]);
dfdy[5] = 1 / c1;
dfdy[6] = c3;
dfdy[7] = 0.0;
dfdy[8] = -c3;
dfdt[0] = 0.0;
dfdt[1] = 0.0;
dfdt[2] = 0.0;
return GSL_SUCCESS;
}
gsl_odeiv2_system rhs_func_oregonator = {
rhs_oregonator,
jac_oregonator,
3,
0
};
/* E5 - a stiff badly scaled chemical problem by Enright, Hull &
Lindberg (1975): Comparing numerical methods for stiff systems of
ODEs. BIT, vol. 15, pp. 10-48.
f0 = -a * y0 - b * y0 * y2 y0(0) = 1.76e-3
f1 = a * y0 - m * c * y1 * y2 y1(0) = 0.0
f2 = a * y0 - b * y0 * y2 - m * c * y1 * y2 + c * y3 y2(0) = 0.0
f3 = b * y0 * y2 - c * y3 y3(0) = 0.0
*/
int
rhs_e5 (double t, const double y[], double f[], void *params)
{
const double a = 7.89e-10;
const double b = 1.1e7;
const double c = 1.13e3;
const double m = 1.0e6;
extern int nfe;
nfe += 1;
f[0] = -a * y[0] - b * y[0] * y[2];
f[1] = a * y[0] - m * c * y[1] * y[2];
f[3] = b * y[0] * y[2] - c * y[3];
f[2] = f[1] - f[3];
return GSL_SUCCESS;
}
int
jac_e5 (double t, const double y[], double *dfdy, double dfdt[], void *params)
{
const double a = 7.89e-10;
const double b = 1.1e7;
const double c = 1.13e3;
const double m = 1.0e6;
extern int nje;
nje += 1;
dfdy[0] = -a - b * y[2];
dfdy[1] = 0.0;
dfdy[2] = -b * y[0];
dfdy[3] = 0.0;
dfdy[4] = a;
dfdy[5] = -m * c * y[2];
dfdy[6] = -m * c * y[1];
dfdy[7] = 0.0;
dfdy[8] = a - b * y[2];
dfdy[9] = -m * c * y[2];
dfdy[10] = -b * y[0] - m * c * y[1];
dfdy[11] = c;
dfdy[12] = b * y[2];
dfdy[13] = 0.0;
dfdy[14] = b * y[0];
dfdy[15] = -c;
dfdt[0] = 0.0;
dfdt[1] = 0.0;
dfdt[2] = 0.0;
dfdt[3] = 0.0;
return GSL_SUCCESS;
}
gsl_odeiv2_system rhs_func_e5 = {
rhs_e5,
jac_e5,
4,
0
};
/* Chemical reaction system of H.H. Robertson (1966): The solution of
a set of reaction rate equations. In: J. Walsh, ed.: Numer.
Anal., an Introduction, Academ. Press, pp. 178-182.
f0 = -a * y0 + b * y1 * y2 y0(0) = 1.0
f1 = a * y0 - b * y1 * y2 - c * y1^2 y1(0) = 0.0
f2 = c * y1^2 y2(0) = 0.0
*/
int
rhs_robertson (double t, const double y[], double f[], void *params)
{
const double a = 0.04;
const double b = 1.0e4;
const double c = 3.0e7;
extern int nfe;
nfe += 1;
f[0] = -a * y[0] + b * y[1] * y[2];
f[2] = c * y[1] * y[1];
f[1] = -f[0] - f[2];
return GSL_SUCCESS;
}
int
jac_robertson (double t, const double y[], double *dfdy, double dfdt[],
void *params)
{
const double a = 0.04;
const double b = 1.0e4;
const double c = 3.0e7;
extern int nje;
nje += 1;
dfdy[0] = -a;
dfdy[1] = b * y[2];
dfdy[2] = b * y[1];
dfdy[3] = a;
dfdy[4] = -b * y[2] - 2 * c * y[1];
dfdy[5] = -b * y[1];
dfdy[6] = 0.0;
dfdy[7] = 2 * c * y[1];
dfdy[8] = 0.0;
dfdt[0] = 0.0;
dfdt[1] = 0.0;
dfdt[2] = 0.0;
return GSL_SUCCESS;
}
gsl_odeiv2_system rhs_func_robertson = {
rhs_robertson,
jac_robertson,
3,
0
};
/* A two-dimensional oscillating Brusselator system.
f0 = a + y0^2 * y1 - (b + 1) * y0 y0(0) = 1.5
f1 = b * y0 - y0^2 * y1 y1(0) = 3.0
*/
int
rhs_brusselator (double t, const double y[], double f[], void *params)
{
const double a = 1.0;
const double b = 3.0;
extern int nfe;
nfe += 1;
f[0] = a + y[0] * y[0] * y[1] - (b + 1.0) * y[0];
f[1] = b * y[0] - y[0] * y[0] * y[1];
return GSL_SUCCESS;
}
int
jac_brusselator (double t, const double y[], double *dfdy, double dfdt[],
void *params)
{
const double b = 3.0;
extern int nje;
nje += 1;
dfdy[0] = 2 * y[0] * y[1] - (b + 1.0);
dfdy[1] = y[0] * y[0];
dfdy[2] = b - 2 * y[0] * y[1];
dfdy[3] = -y[0] * y[0];
dfdt[0] = 0;
dfdt[1] = 0;
return GSL_SUCCESS;
}
gsl_odeiv2_system rhs_func_brusselator = {
rhs_brusselator,
jac_brusselator,
2,
0
};
/* Ring Modulator, stiff ODE of dimension 15.
Reference: Walter M. Lioen, Jacques J.B. de Swart, Test Set for
Initial Value Problem Solvers, Release 2.1 September 1999,
http://ftp.cwi.nl/IVPtestset/software.htm
*/
#define NRINGMOD 15
int
rhs_ringmod (double t, const double y[], double f[], void *params)
{
const double c = 1.6e-8;
const double cs = 2e-12;
const double cp = 1e-8;
const double r = 25e3;
const double rp = 50e0;
const double lh = 4.45e0;
const double ls1 = 2e-3;
const double ls2 = 5e-4;
const double ls3 = 5e-4;
const double rg1 = 36.3;
const double rg2 = 17.3;
const double rg3 = 17.3;
const double ri = 5e1;
const double rc = 6e2;
const double gamma = 40.67286402e-9;
const double delta = 17.7493332;
const double pi = 3.141592653589793238462643383;
const double uin1 = 0.5 * sin (2e3 * pi * t);
const double uin2 = 2 * sin (2e4 * pi * t);
const double ud1 = +y[2] - y[4] - y[6] - uin2;
const double ud2 = -y[3] + y[5] - y[6] - uin2;
const double ud3 = +y[3] + y[4] + y[6] + uin2;
const double ud4 = -y[2] - y[5] + y[6] + uin2;
const double qud1 = gamma * (exp (delta * ud1) - 1.0);
const double qud2 = gamma * (exp (delta * ud2) - 1.0);
const double qud3 = gamma * (exp (delta * ud3) - 1.0);
const double qud4 = gamma * (exp (delta * ud4) - 1.0);
extern int nfe;
nfe += 1;
f[0] = (y[7] - 0.5 * y[9] + 0.5 * y[10] + y[13] - y[0] / r) / c;
f[1] = (y[8] - 0.5 * y[11] + 0.5 * y[12] + y[14] - y[1] / r) / c;
f[2] = (y[9] - qud1 + qud4) / cs;
f[3] = (-y[10] + qud2 - qud3) / cs;
f[4] = (y[11] + qud1 - qud3) / cs;
f[5] = (-y[12] - qud2 + qud4) / cs;
f[6] = (-y[6] / rp + qud1 + qud2 - qud3 - qud4) / cp;
f[7] = -y[0] / lh;
f[8] = -y[1] / lh;
f[9] = (0.5 * y[0] - y[2] - rg2 * y[9]) / ls2;
f[10] = (-0.5 * y[0] + y[3] - rg3 * y[10]) / ls3;
f[11] = (0.5 * y[1] - y[4] - rg2 * y[11]) / ls2;
f[12] = (-0.5 * y[1] + y[5] - rg3 * y[12]) / ls3;
f[13] = (-y[0] + uin1 - (ri + rg1) * y[13]) / ls1;
f[14] = (-y[1] - (rc + rg1) * y[14]) / ls1;
return GSL_SUCCESS;
}
int
jac_ringmod (double t, const double y[], double *dfdy, double dfdt[],
void *params)
{
const double c = 1.6e-8;
const double cs = 2e-12;
const double cp = 1e-8;
const double r = 25e3;
const double rp = 50e0;
const double lh = 4.45e0;
const double ls1 = 2e-3;
const double ls2 = 5e-4;
const double ls3 = 5e-4;
const double rg1 = 36.3;
const double rg2 = 17.3;
const double rg3 = 17.3;
const double ri = 5e1;
const double rc = 6e2;
const double gamma = 40.67286402e-9;
const double delta = 17.7493332;
const double pi = 3.141592653589793238462643383;
const double uin2 = 2 * sin (2e4 * pi * t);
const double ud1 = +y[2] - y[4] - y[6] - uin2;
const double ud2 = -y[3] + y[5] - y[6] - uin2;
const double ud3 = +y[3] + y[4] + y[6] + uin2;
const double ud4 = -y[2] - y[5] + y[6] + uin2;
const double qpud1 = gamma * delta * exp (delta * ud1);
const double qpud2 = gamma * delta * exp (delta * ud2);
const double qpud3 = gamma * delta * exp (delta * ud3);
const double qpud4 = gamma * delta * exp (delta * ud4);
extern int nje;
size_t i;
nje += 1;
for (i = 0; i < NRINGMOD * NRINGMOD; i++)
{
dfdy[i] = 0.0;
}
dfdy[0 * NRINGMOD + 0] = -1 / (c * r);
dfdy[0 * NRINGMOD + 7] = 1 / c;
dfdy[0 * NRINGMOD + 9] = -0.5 / c;
dfdy[0 * NRINGMOD + 10] = -dfdy[0 * NRINGMOD + 9];
dfdy[0 * NRINGMOD + 13] = dfdy[0 * NRINGMOD + 7];
dfdy[1 * NRINGMOD + 1] = dfdy[0 * NRINGMOD + 0];
dfdy[1 * NRINGMOD + 8] = dfdy[0 * NRINGMOD + 7];
dfdy[1 * NRINGMOD + 11] = dfdy[0 * NRINGMOD + 9];
dfdy[1 * NRINGMOD + 12] = dfdy[0 * NRINGMOD + 10];
dfdy[1 * NRINGMOD + 14] = dfdy[0 * NRINGMOD + 13];
dfdy[2 * NRINGMOD + 2] = (-qpud1 - qpud4) / cs;
dfdy[2 * NRINGMOD + 4] = qpud1 / cs;
dfdy[2 * NRINGMOD + 5] = -qpud4 / cs;
dfdy[2 * NRINGMOD + 6] = (qpud1 + qpud4) / cs;
dfdy[2 * NRINGMOD + 9] = 1 / cs;
dfdy[3 * NRINGMOD + 3] = (-qpud2 - qpud3) / cs;
dfdy[3 * NRINGMOD + 4] = -qpud3 / cs;
dfdy[3 * NRINGMOD + 5] = qpud2 / cs;
dfdy[3 * NRINGMOD + 6] = (-qpud2 - qpud3) / cs;
dfdy[3 * NRINGMOD + 10] = -1 / cs;
dfdy[4 * NRINGMOD + 2] = qpud1 / cs;
dfdy[4 * NRINGMOD + 3] = -qpud3 / cs;
dfdy[4 * NRINGMOD + 4] = (-qpud1 - qpud3) / cs;
dfdy[4 * NRINGMOD + 6] = (-qpud1 - qpud3) / cs;
dfdy[4 * NRINGMOD + 11] = 1 / cs;
dfdy[5 * NRINGMOD + 2] = -qpud4 / cs;
dfdy[5 * NRINGMOD + 3] = qpud2 / cs;
dfdy[5 * NRINGMOD + 5] = (-qpud2 - qpud4) / cs;
dfdy[5 * NRINGMOD + 6] = (qpud2 + qpud4) / cs;
dfdy[5 * NRINGMOD + 12] = -1 / cs;
dfdy[6 * NRINGMOD + 2] = (qpud1 + qpud4) / cp;
dfdy[6 * NRINGMOD + 3] = (-qpud2 - qpud3) / cp;
dfdy[6 * NRINGMOD + 4] = (-qpud1 - qpud3) / cp;
dfdy[6 * NRINGMOD + 5] = (qpud2 + qpud4) / cp;
dfdy[6 * NRINGMOD + 6] = (-qpud1 - qpud2 - qpud3 - qpud4 - 1 / rp) / cp;
dfdy[7 * NRINGMOD + 0] = -1 / lh;
dfdy[8 * NRINGMOD + 1] = dfdy[7 * NRINGMOD + 0];
dfdy[9 * NRINGMOD + 0] = 0.5 / ls2;
dfdy[9 * NRINGMOD + 2] = -1 / ls2;
dfdy[9 * NRINGMOD + 9] = -rg2 / ls2;
dfdy[10 * NRINGMOD + 0] = -0.5 / ls3;
dfdy[10 * NRINGMOD + 3] = 1 / ls3;
dfdy[10 * NRINGMOD + 10] = -rg3 / ls3;
dfdy[11 * NRINGMOD + 1] = dfdy[9 * NRINGMOD + 0];
dfdy[11 * NRINGMOD + 4] = dfdy[9 * NRINGMOD + 2];
dfdy[11 * NRINGMOD + 11] = dfdy[9 * NRINGMOD + 9];
dfdy[12 * NRINGMOD + 1] = dfdy[10 * NRINGMOD + 0];
dfdy[12 * NRINGMOD + 5] = dfdy[10 * NRINGMOD + 3];
dfdy[12 * NRINGMOD + 12] = dfdy[10 * NRINGMOD + 10];
dfdy[13 * NRINGMOD + 0] = -1 / ls1;
dfdy[13 * NRINGMOD + 13] = -(ri + rg1) / ls1;
dfdy[14 * NRINGMOD + 1] = dfdy[13 * NRINGMOD + 0];
dfdy[14 * NRINGMOD + 14] = -(rc + rg1) / ls1;
for (i = 0; i < NRINGMOD; i++)
{
dfdt[i] = 0.0;
}
return GSL_SUCCESS;
}
gsl_odeiv2_system rhs_func_ringmod = {
rhs_ringmod,
jac_ringmod,
NRINGMOD,
NULL
};
/**********************************************************/
/* Functions for carrying out tests */
/**********************************************************/
void
test_odeiv_stepper (const gsl_odeiv2_step_type * T,
const gsl_odeiv2_system * sys, const double h,
const double t, const char desc[], const double ystart[],
const double yfin[], const double relerr)
{
/* tests stepper T with one fixed length step advance of system sys
and compares the result with given values yfin
*/
double y[MAXEQ] = { 0.0 };
double yerr[MAXEQ] = { 0.0 };
double scale_abs[MAXEQ];
size_t ne = sys->dimension;
size_t i;
gsl_odeiv2_driver *d;
for (i = 0; i < MAXEQ; i++)
{
scale_abs[i] = 1.0;
}
d = gsl_odeiv2_driver_alloc_scaled_new (sys, T, h, relerr, relerr,
1.0, 0.0, scale_abs);
DBL_MEMCPY (y, ystart, MAXEQ);
{
int s = gsl_odeiv2_step_apply (d->s, t, h, y, yerr, 0, 0, sys);
if (s != GSL_SUCCESS)
{
gsl_test (s, "test_odeiv_stepper: %s step_apply returned %d", desc,
s);
}
}
for (i = 0; i < ne; i++)
{
gsl_test_rel (y[i], yfin[i], relerr,
"%s %s step(%d)", gsl_odeiv2_step_name (d->s), desc, i);
}
gsl_odeiv2_driver_free (d);
}
void
test_stepper (const gsl_odeiv2_step_type * T)
{
/* Tests stepper T with a step of selected systems */
double y[MAXEQ] = { 0.0 };
double yfin[MAXEQ] = { 0.0 };
/* Step length */
double h;
/* Required tolerance */
double err_target;
/* classic stiff */
h = 1e-7;
err_target = 1e-4;
y[0] = 1.0;
y[1] = 0.0;
{
const double e1 = exp (-h);
const double e2 = exp (-1000.0 * h);
yfin[0] = 2.0 * e1 - e2;
yfin[1] = -e1 + e2;
}
test_odeiv_stepper (T, &rhs_func_stiff, h, 0.0, "classic_stiff",
y, yfin, err_target);
/* linear */
h = 1e-1;
err_target = 1e-10;
y[0] = 0.58;
yfin[0] = y[0] + 2 * h;
test_odeiv_stepper (T, &rhs_func_lin, h, 0.0, "linear",
y, yfin, err_target);
/* exponential */
h = 1e-4;
err_target = 1e-8;
y[0] = exp (2.7);
yfin[0] = exp (2.7 + h);
test_odeiv_stepper (T, &rhs_func_exp, h, 2.7, "exponential",
y, yfin, err_target);
/* cosine-sine */
h = 1e-3;
err_target = 1e-6;
y[0] = cos (1.2);
y[1] = sin (1.2);
yfin[0] = cos (1.2 + h);
yfin[1] = sin (1.2 + h);
test_odeiv_stepper (T, &rhs_func_sin, h, 1.2, "cosine-sine",
y, yfin, err_target);
}
void
test_evolve_system (const gsl_odeiv2_step_type * T,
const gsl_odeiv2_system * sys,
double t0, double t1, double hstart,
double y[], double yfin[],
double err_target, const char *desc)
{
/* Tests system sys with stepper T. Step length is controlled by
error estimation from the stepper.
*/
int steps = 0;
size_t i;
double t = t0;
double h = hstart;
/* Tolerance factor in testing errors */
const double factor = 10;
gsl_odeiv2_driver *d =
gsl_odeiv2_driver_alloc_standard_new (sys, T, h, err_target, err_target,
1.0, 0.0);
double *y_orig = (double *) malloc (sys->dimension * sizeof (double));
int s = 0;
extern int nfe, nje;
nfe = 0;
nje = 0;
while (t < t1)
{
double t_orig = t;
memcpy (y_orig, y, sys->dimension * sizeof (double));
s = gsl_odeiv2_evolve_apply (d->e, d->c, d->s, sys, &t, t1, &h, y);
#ifdef DEBUG
printf ("test_evolve_system at: %.5e %.5e %.5e %d\n", t, y[0], y[1],
gsl_odeiv2_step_order (d->s));
#endif
if (s != GSL_SUCCESS)
{
/* check that t and y are unchanged */
gsl_test_abs (t, t_orig, 0.0, "%s t must be restored on failure",
gsl_odeiv2_step_name (d->s));
for (i = 0; i < sys->dimension; i++)
{
gsl_test_abs (y[i], y_orig[i], 0.0,
"%s y must be restored on failure",
gsl_odeiv2_step_name (d->s), desc, i);
}
if (sys != &rhs_func_xsin)
{
/* apart from xsin, other functions should not return errors */
gsl_test (s, "%s evolve_apply returned %d",
gsl_odeiv2_step_name (d->s), s);
break;
}
}
if (steps > 100000)
{
gsl_test (GSL_EFAILED,
"%s evolve_apply reached maxiter",
gsl_odeiv2_step_name (d->s));
break;
}
steps++;
}
gsl_test (s, "%s %s [%g,%g], %d steps (nfe %d, nje %d) completed",
gsl_odeiv2_step_name (d->s), desc, t0, t1, steps, nfe, nje);
/* err_target is target error of one step. Test if stepper has made
larger error than (tolerance factor times) the number of steps
times the err_target.
*/
for (i = 0; i < sys->dimension; i++)
{
gsl_test_abs (y[i], yfin[i], factor * d->e->count * err_target,
"%s %s evolve(%d)", gsl_odeiv2_step_name (d->s), desc, i);
}
free (y_orig);
gsl_odeiv2_driver_free (d);
}
int
sys_driver (const gsl_odeiv2_step_type * T,
const gsl_odeiv2_system * sys,
double t0, double t1, double hstart,
double y[], double epsabs, double epsrel, const char desc[])
{
/* This function evolves a system sys with stepper T from t0 to t1.
Step length is varied via error control with possibly different
absolute and relative error tolerances.
*/
int s = 0;
int steps = 0;
double t = t0;
double h = hstart;
gsl_odeiv2_driver *d =
gsl_odeiv2_driver_alloc_standard_new (sys, T, h, epsabs, epsrel,
1.0, 0.0);
extern int nfe, nje;
nfe = 0;
nje = 0;
while (t < t1)
{
s = gsl_odeiv2_evolve_apply (d->e, d->c, d->s, sys, &t, t1, &h, y);
#ifdef DEBUG
printf ("sys_driver at: %.5e %.5e %.5e %d\n", t, y[0], y[1],
gsl_odeiv2_step_order (d->s));
#endif
if (s != GSL_SUCCESS)
{
gsl_test (s, "sys_driver: %s evolve_apply returned %d",
gsl_odeiv2_step_name (d->s), s);
break;
}
if (steps > 1e7)
{
gsl_test (GSL_EMAXITER,
"sys_driver: %s evolve_apply reached maxiter at t=%g",
gsl_odeiv2_step_name (d->s), t);
s = GSL_EMAXITER;
break;
}
steps++;
}
gsl_test (s, "%s %s [%g,%g], %d steps (nfe %d, nje %d) completed",
gsl_odeiv2_step_name (d->s), desc, t0, t1, steps, nfe, nje);
gsl_odeiv2_driver_free (d);
return s;
}
void
test_evolve_linear (const gsl_odeiv2_step_type * T, double h, double err)
{
/* Test linear evolve */
double y[1];
double yfin[1];
y[0] = 1.0;
yfin[0] = 9.0;
test_evolve_system (T, &rhs_func_lin, 0.0, 4.0, h, y, yfin, err,
"linear[0,4]");
}
void
test_evolve_exp (const gsl_odeiv2_step_type * T, double h, double err)
{
/* Test exponential evolve */
double y[1];
double yfin[1];
y[0] = 1.0;
yfin[0] = exp (2.0);
test_evolve_system (T, &rhs_func_exp, 0.0, 2.0, h, y, yfin, err,
"exp[0,2]");
}
void
test_evolve_sin (const gsl_odeiv2_step_type * T, double h, double err)
{
/* Test sinusoidal evolve */
double y[2];
double yfin[2];
y[0] = 1.0;
y[1] = 0.0;
yfin[0] = cos (2.0);
yfin[1] = sin (2.0);
test_evolve_system (T, &rhs_func_sin, 0.0, 2.0, h, y, yfin, err,
"sine[0,2]");
}
void
test_evolve_xsin (const gsl_odeiv2_step_type * T, double h, double err)
{
/* Test sinusoidal evolve including a failing window */
double y[2];
double yfin[2];
y[0] = 1.0;
y[1] = 0.0;
yfin[0] = cos (2.0);
yfin[1] = sin (2.0);
rhs_xsin_reset = 1;
jac_xsin_reset = 1;
test_evolve_system (T, &rhs_func_xsin, 0.0, 2.0, h, y, yfin, err,
"sine[0,2] w/errors");
}
void
test_evolve_stiff1 (const gsl_odeiv2_step_type * T, double h, double err)
{
/* Test classical stiff problem evolve, t=[0,1] */
double y[2];
double yfin[2];
y[0] = 1.0;
y[1] = 0.0;
{
double arg = 1.0;
double e1 = exp (-arg);
double e2 = exp (-1000.0 * arg);
yfin[0] = 2.0 * e1 - e2;
yfin[1] = -e1 + e2;
}
test_evolve_system (T, &rhs_func_stiff, 0.0, 1.0, h, y, yfin, err,
"stiff[0,1]");
}
void
test_evolve_stiff5 (const gsl_odeiv2_step_type * T, double h, double err)
{
/* Test classical stiff problem evolve, t=[0,5] */
double y[2];
double yfin[2];
y[0] = 1.0;
y[1] = 0.0;
{
double arg = 5.0;
double e1 = exp (-arg);
double e2 = exp (-1000.0 * arg);
yfin[0] = 2.0 * e1 - e2;
yfin[1] = -e1 + e2;
}
test_evolve_system (T, &rhs_func_stiff, 0.0, 5.0, h, y, yfin, err,
"stiff[0,5]");
}
void
test_evolve_negative_h (const gsl_odeiv2_step_type * T, double h, double err)
{
/* Test evolution in negative direction */
/* Tolerance factor in testing errors */
const double factor = 10;
const gsl_odeiv2_system sys = rhs_func_cos;
double t = 0;
double t1 = -4.0;
double y = 0.0;
double yfin = sin (t1);
gsl_odeiv2_driver *d;
/* Make initial h negative */
h = -fabs (h);
d = gsl_odeiv2_driver_alloc_standard_new (&sys, T, h, err, err, 1.0, 0.0);
while (t > t1)
{
int status = gsl_odeiv2_evolve_apply (d->e, d->c, d->s,
&sys, &t, t1, &h, &y);
if (status != GSL_SUCCESS)
{
gsl_test (status, "%s evolve_apply returned %d for negative h",
gsl_odeiv2_step_name (d->s), status);
break;
}
}
gsl_test_abs (y, yfin, factor * d->e->count * err,
"%s evolution with negative h", gsl_odeiv2_step_name (d->s));
gsl_odeiv2_driver_free (d);
}
void
test_broken (const gsl_odeiv2_step_type * T, double h, double err)
{
/* Check for gsl_odeiv2_evolve_apply. The user function fails at
t>=10, which in this test causes step size to decrease below
machine precision and return with a failure code.
*/
/* Tolerance factor in testing errors */
const double factor = 10;
const gsl_odeiv2_system sys = rhs_func_broken;
gsl_odeiv2_driver *d = gsl_odeiv2_driver_alloc_y_new (&sys, T, h, err, err);
double t = 0;
double t1 = 100.0;
double y = 0.0;
const double final_max_h = GSL_DBL_EPSILON;
int status;
while (t < t1)
{
status =
gsl_odeiv2_evolve_apply (d->e, d->c, d->s, &sys, &t, t1, &h, &y);
if (status != GSL_SUCCESS)
{
break;
}
}
gsl_test_abs (h + final_max_h, final_max_h, factor * final_max_h,
"%s test_broken: step size at break point",
gsl_odeiv2_step_name (d->s));
gsl_test_abs (t, 10.0, factor * err,
"%s test_broken: point of break",
gsl_odeiv2_step_name (d->s));
/* GSL_FAILURE results from stepper failure, 123 from user function */
if (status != GSL_FAILURE && status != 123)
{
gsl_test (status, "%s test_broken: evolve return value %d",
gsl_odeiv2_step_name (d->s), status);
}
else
{
gsl_test (GSL_SUCCESS, "%s test_broken: evolve return value %d",
gsl_odeiv2_step_name (d->s), status);
}
gsl_odeiv2_driver_free (d);
}
void
test_stepsize_fail (const gsl_odeiv2_step_type * T, double h)
{
/* Check for gsl_odeiv2_evolve_apply. The user function works
apparently fine (returns GSL_SUCCESS) but is faulty in this case.
gsl_odeiv2_evolve_apply therefore decreases the step-size below
machine precision and therefore stops with GSL_FAILURE.
*/
/* Error tolerance */
const double epsabs = 1e-16;
const double epsrel = 1e-6;
/* Tolerance factor in testing errors */
const double factor = 10;
const double final_max_h = GSL_DBL_EPSILON;
const gsl_odeiv2_system sys = rhs_func_br;
gsl_odeiv2_driver *d =
gsl_odeiv2_driver_alloc_y_new (&sys, T, h, epsabs, epsrel);
double t = 1.0;
double t1 = 1e5;
double y = 0.0;
int status;
while (t < t1)
{
status =
gsl_odeiv2_evolve_apply (d->e, d->c, d->s, &sys, &t, t1, &h, &y);
if (status != GSL_SUCCESS)
{
break;
}
}
gsl_test_abs (h + final_max_h, final_max_h, factor * final_max_h,
"%s test_stepsize_fail: step size at break point",
gsl_odeiv2_step_name (d->s));
gsl_test_abs (t, 1.0, 1e-6,
"%s test_stepsize_fail: point of break",
gsl_odeiv2_step_name (d->s));
gsl_test_int (status, GSL_FAILURE,
"%s test_stepsize_fail: evolve return value",
gsl_odeiv2_step_name (d->s));
gsl_odeiv2_driver_free (d);
}
void
test_user_break (const gsl_odeiv2_step_type * T, double h)
{
/* Tests for user signaled immediate break */
const double tol = 1e-8;
gsl_odeiv2_driver *d = gsl_odeiv2_driver_alloc_y_new (&rhs_func_sin_ub, T,
h, tol, tol);
double y[] = { 1.0, 0.0 };
double t = 0.0;
const double t1 = 8.25;
int s = gsl_odeiv2_driver_apply (d, &t, t1, y);
gsl_test ((s - GSL_EBADFUNC), "%s test_user_break returned %d",
gsl_odeiv2_step_name (d->s), s);
gsl_odeiv2_driver_free (d);
}
void
test_stepfn (const gsl_odeiv2_step_type * T)
{
/* Test evolve on a small derivative change at t=0 */
double epsabs = 1e-16;
double epsrel = 1e-6;
const gsl_odeiv2_system sys = rhs_func_stepfn;
double t = 0.0;
double h = 1e-6;
double y = 0.0;
int i = 0;
int status = 0;
gsl_odeiv2_driver *d =
gsl_odeiv2_driver_alloc_y_new (&sys, T, h, epsabs, epsrel);
while (t < 2 && i < 1000000)
{
status =
gsl_odeiv2_evolve_apply (d->e, d->c, d->s, &sys, &t, 2, &h, &y);
#ifdef DEBUG
printf ("test_stepfn at: i=%d status=%d t=%g h=%g y=%g\n",
i, status, t, h, y);
#endif
if (status != GSL_SUCCESS)
break;
i++;
}
gsl_test (status, "evolve step function, return value (stepfn/%s): %d",
gsl_odeiv2_step_name (d->s), status);
gsl_test_abs (t, 2, 1e-16, "evolve step function, t (stepfn/%s)",
gsl_odeiv2_step_name (d->s));
gsl_test_rel (y, 1, epsrel, "evolve step function, y (stepfn/%s)",
gsl_odeiv2_step_name (d->s));
gsl_odeiv2_driver_free (d);
}
void
test_stepfn2 (const gsl_odeiv2_step_type * T)
{
/* Test evolve on a large derivative change at t=0 */
double epsabs = 1e-16;
double epsrel = 1e-6;
const gsl_odeiv2_system sys = rhs_func_stepfn2;
double t = -1.0;
double h = 1e-6;
double y = 0.0;
int i = 0;
int status;
const int maxiter = 100000;
gsl_odeiv2_driver *d =
gsl_odeiv2_driver_alloc_yp_new (&sys, T, h, epsabs, epsrel);
while (t < 1.0 && i < maxiter)
{
status =
gsl_odeiv2_evolve_apply (d->e, d->c, d->s, &sys, &t, 1.0, &h, &y);
#ifdef DEBUG
printf ("test_stepfn2 at: i=%d status=%d t=%g h=%g y=%g\n",
i, status, t, h, y);
#endif
if (status != GSL_SUCCESS)
break;
i++;
}
if (i >= maxiter)
printf ("FAIL: evolve big step function, (stepfn2/%s) reached maxiter\n",
gsl_odeiv2_step_name (d->s));
gsl_test_abs (t, 1.0, 1e-16, "evolve big step function, t (stepfn2/%s)",
gsl_odeiv2_step_name (d->s));
gsl_test_rel (y, 1e300, epsrel, "evolve big step function, y (stepfn2/%s)",
gsl_odeiv2_step_name (d->s));
gsl_odeiv2_driver_free (d);
}
void
test_nonstiff_problems (void)
{
/* Compares output of non-stiff (or only moderately stiff) problems
with several steppers
*/
const gsl_odeiv2_step_type *steppers[MAXNS];
/* Required error tolerance for each stepper. */
double err_target[MAXNS];
/* initial values for each ode-solver */
double y[MAXEQ * MAXNS];
size_t i, k, p;
/* Number of problems to test */
#define CONST_NONSTIFF_NPROB 4
/* Problems, their names and number of equations */
const gsl_odeiv2_system *prob[] =
{ &rhs_func_vl, &rhs_func_vanderpol, &rhs_func_stifftrig,
&rhs_func_brusselator
};
const char *probname[] =
{ "volterra-lotka", "vanderpol", "stifftrig", "brusselator" };
const size_t sd[] = { 2, 2, 1, 2 };
/* Integration interval for problems */
const double start[CONST_NONSTIFF_NPROB] = { 0.0 };
const double end[] = { 9.0, 100.0, 1.5, 20.0 };
const double epsabs = 1e-8;
const double epsrel = 1e-8;
const double initstepsize = 1e-5;
/* Steppers */
steppers[0] = gsl_odeiv2_step_rk4;
err_target[0] = 1e-6;
steppers[1] = gsl_odeiv2_step_rk2;
err_target[1] = 1e-6;
steppers[2] = gsl_odeiv2_step_rkf45;
err_target[2] = 1e-6;
steppers[3] = gsl_odeiv2_step_rkck;
err_target[3] = 1e-6;
steppers[4] = gsl_odeiv2_step_rk8pd;
err_target[4] = 1e-6;
steppers[5] = gsl_odeiv2_step_rk1imp;
err_target[5] = 1e-3;
steppers[6] = gsl_odeiv2_step_rk2imp;
err_target[6] = 1e-5;
steppers[7] = gsl_odeiv2_step_rk4imp;
err_target[7] = 1e-6;
steppers[8] = gsl_odeiv2_step_bsimp;
err_target[8] = 1e-6;
steppers[9] = gsl_odeiv2_step_msadams;
err_target[9] = 1e-5;
steppers[10] = gsl_odeiv2_step_msbdf;
err_target[10] = 1e-5;
steppers[11] = 0;
/* Loop over problems */
for (p = 0; p < CONST_NONSTIFF_NPROB; p++)
{
/* Initialize */
for (i = 0; i < MAXNS * MAXEQ; i++)
{
y[i] = 0.0;
}
for (i = 0; i < MAXNS; i++)
{
switch (p)
{
case 0:
y[i * sd[p]] = 2.725;
y[i * sd[p] + 1] = 1.0;
break;
case 1:
y[i * sd[p]] = 1.0;
y[i * sd[p] + 1] = 0.0;
break;
case 2:
y[i * sd[p]] = 0.0;
break;
case 3:
y[i * sd[p]] = 1.5;
y[i * sd[p] + 1] = 3.0;
break;
default:
gsl_test (GSL_EFAILED,
"test_nonstiff_problems: initialization error\n");
return;
}
}
/* Call each solver for the problem */
for (i = 0; steppers[i] != 0; i++)
{
int s = sys_driver (steppers[i], prob[p], start[p], end[p],
initstepsize, &y[sd[p] * i],
epsabs, epsrel, probname[p]);
if (s != GSL_SUCCESS)
{
gsl_test (s, "test_nonstiff_problems %s %s",
steppers[i]->name, probname[p]);
}
}
/* Compare results */
for (i = 1; steppers[i] != 0; i++)
for (k = 0; k < sd[p]; k++)
{
const double val1 = y[k];
const double val2 = y[sd[p] * i + k];
gsl_test_rel (val1, val2,
(GSL_MAX (err_target[0], err_target[i])),
"%s/%s %s [%d]",
steppers[0]->name, steppers[i]->name,
probname[p], k);
}
}
}
void
test_stiff_problems (void)
{
/* Compares output of stiff problems with several steppers */
const gsl_odeiv2_step_type *steppers[MAXNS];
/* Required error tolerance for each stepper. */
double err_target[MAXNS];
/* initial values for each ode-solver */
double y[MAXEQ * MAXNS];
size_t i, k, p;
/* Number of problems to test */
#define CONST_STIFF_NPROB 3
/* Problems, their names and number of equations */
const gsl_odeiv2_system *prob[] =
{ &rhs_func_oregonator, &rhs_func_e5, &rhs_func_robertson };
const char *probname[] = { "oregonator", "e5", "robertson" };
const size_t sd[] = { 3, 4, 3 };
/* Integration interval for problems */
const double start[CONST_STIFF_NPROB] = { 0.0 };
const double end[] = { 360.0, 1e4, 1e4 };
const double epsabs = 1e-40;
const double epsrel = 1e-7;
const double initstepsize = 1e-5;
/* Steppers */
steppers[0] = gsl_odeiv2_step_bsimp;
err_target[0] = 1e-6;
steppers[1] = gsl_odeiv2_step_rk1imp;
err_target[1] = 5e-3;
steppers[2] = gsl_odeiv2_step_rk2imp;
err_target[2] = 5e-5;
steppers[3] = gsl_odeiv2_step_rk4imp;
err_target[3] = 5e-5;
steppers[4] = gsl_odeiv2_step_msbdf;
err_target[4] = 1e-4;
steppers[5] = 0;
/* Loop over problems */
for (p = 0; p < CONST_STIFF_NPROB; p++)
{
/* Initialize */
for (i = 0; i < MAXNS * MAXEQ; i++)
{
y[i] = 0.0;
}
for (i = 0; i < MAXNS; i++)
{
switch (p)
{
case 0:
y[i * sd[p]] = 1.0;
y[i * sd[p] + 1] = 2.0;
y[i * sd[p] + 2] = 3.0;
break;
case 1:
y[i * sd[p]] = 1.76e-3;
y[i * sd[p] + 1] = 0.0;
y[i * sd[p] + 2] = 0.0;
y[i * sd[p] + 3] = 0.0;
break;
case 2:
y[i * sd[p]] = 1.0;
y[i * sd[p] + 1] = 0.0;
y[i * sd[p] + 2] = 0.0;
break;
default:
gsl_test (GSL_EFAILED,
"test_stiff_problems: initialization error\n");
return;
}
}
/* Call each solver for the problem */
for (i = 0; steppers[i] != 0; i++)
{
int s = sys_driver (steppers[i], prob[p], start[p], end[p],
initstepsize, &y[sd[p] * i],
epsabs, epsrel, probname[p]);
if (s != GSL_SUCCESS)
{
gsl_test (s, "test_stiff_problems %s %s",
steppers[i]->name, probname[p]);
}
}
/* Compare results */
for (i = 1; steppers[i] != 0; i++)
for (k = 0; k < sd[p]; k++)
{
const double val1 = y[k];
const double val2 = y[sd[p] * i + k];
gsl_test_rel (val1, val2,
(GSL_MAX (err_target[0], err_target[i])),
"%s/%s %s [%d]",
steppers[0]->name, steppers[i]->name,
probname[p], k);
}
}
}
void
test_extreme_problems (void)
{
/* Compares output of numerically particularly demanding problems
with several steppers */
const gsl_odeiv2_step_type *steppers[MAXNS];
/* Required error tolerance for each stepper. */
double err_target[MAXNS];
/* initial values for each ode-solver */
double y[MAXEQ * MAXNS];
size_t i, k, p;
/* Number of problems to test */
#define CONST_EXTREME_NPROB 3
/* Problems, their names and number of equations */
const gsl_odeiv2_system *prob[] =
{ &rhs_func_e5, &rhs_func_robertson, &rhs_func_ringmod };
const char *probname[] = { "e5_bigt", "robertson_bigt", "ringmod" };
const size_t sd[] = { 4, 3, NRINGMOD };
/* Integration interval for problems */
const double start[CONST_EXTREME_NPROB] = { 0.0 };
const double end[CONST_EXTREME_NPROB] = { 1e11, 1e11, 1e-5 };
const double epsabs[CONST_EXTREME_NPROB] =
{ 1e1 * GSL_DBL_MIN, 1e1 * GSL_DBL_MIN, 1e-12 };
const double epsrel[CONST_EXTREME_NPROB] = { 1e-12, 1e-12, 1e-12 };
const double initstepsize[CONST_EXTREME_NPROB] = { 1e-5, 1e-5, 1e-10 };
/* Steppers */
steppers[0] = gsl_odeiv2_step_bsimp;
err_target[0] = 1e-3;
steppers[1] = gsl_odeiv2_step_msbdf;
err_target[1] = 1e-3;
steppers[2] = 0;
/* Loop over problems */
for (p = 0; p < CONST_EXTREME_NPROB; p++)
{
/* Initialize */
for (i = 0; i < MAXNS * MAXEQ; i++)
{
y[i] = 0.0;
}
for (i = 0; i < MAXNS; i++)
{
switch (p)
{
case 0:
y[i * sd[p]] = 1.76e-3;
y[i * sd[p] + 1] = 0.0;
y[i * sd[p] + 2] = 0.0;
y[i * sd[p] + 3] = 0.0;
break;
case 1:
y[i * sd[p]] = 1.0;
y[i * sd[p] + 1] = 0.0;
y[i * sd[p] + 2] = 0.0;
break;
case 2:
{
size_t j;
for (j = 0; j < NRINGMOD; j++)
{
y[i * sd[p] + j] = 0.0;
}
}
break;
default:
gsl_test (GSL_EFAILED,
"test_extreme_problems: initialization error\n");
return;
}
}
/* Call each solver for the problem */
for (i = 0; steppers[i] != 0; i++)
{
int s = sys_driver (steppers[i], prob[p], start[p], end[p],
initstepsize[p], &y[sd[p] * i],
epsabs[p], epsrel[p], probname[p]);
if (s != GSL_SUCCESS)
{
printf ("start=%.5e, initstepsize=%.5e\n", start[p],
initstepsize[p]);
gsl_test (s, "test_extreme_problems %s %s", steppers[i]->name,
probname[p]);
}
}
/* Compare results */
for (i = 1; steppers[i] != 0; i++)
for (k = 0; k < sd[p]; k++)
{
const double val1 = y[k];
const double val2 = y[sd[p] * i + k];
gsl_test_rel (val1, val2,
(GSL_MAX (err_target[0], err_target[i])),
"%s/%s %s [%d]",
steppers[0]->name, steppers[i]->name,
probname[p], k);
}
}
}
void
test_driver (const gsl_odeiv2_step_type * T)
{
/* Tests for gsl_odeiv2_driver object */
int s;
const double tol = 1e-8;
const double hmax = 1e-2;
double y[] = { 1.0, 0.0 };
double t = 0.0;
const double tinit = 0.0;
const double t1 = 8.25;
const double t2 = 100;
const double t3 = -2.5;
const size_t minsteps = ceil (t1 / hmax);
const size_t maxsteps = 20;
const double hmin = 1e-10;
const unsigned long int nfsteps = 100000;
const double hfixed = 0.000025;
gsl_odeiv2_driver *d = gsl_odeiv2_driver_alloc_y_new (&rhs_func_sin, T,
1e-3, tol, tol);
gsl_odeiv2_driver_set_hmax (d, hmax);
s = gsl_odeiv2_driver_apply (d, &t, t1, y);
if (s != GSL_SUCCESS)
{
gsl_test (s, "%s test_driver apply returned %d",
gsl_odeiv2_step_name (d->s), s);
}
/* Test that result is correct */
gsl_test_rel (y[0], cos (t1), d->n * tol, "%s test_driver y0",
gsl_odeiv2_step_name (d->s));
gsl_test_rel (y[1], sin (t1), d->n * tol, "%s test_driver y1",
gsl_odeiv2_step_name (d->s));
/* Test that maximum step length is obeyed */
if (d->n < minsteps)
{
gsl_test (1, "%s test_driver steps %d < minsteps %d \n",
gsl_odeiv2_step_name (d->s), d->n, minsteps);
}
else
{
gsl_test (0, "%s test_driver max step length test",
gsl_odeiv2_step_name (d->s));
}
/* Test changing integration direction from forward to backward */
gsl_odeiv2_driver_reset_hstart (d, -1e-3);
s = gsl_odeiv2_driver_apply (d, &t, tinit, y);
if (s != GSL_SUCCESS)
{
gsl_test (s, "%s test_driver apply returned %d",
gsl_odeiv2_step_name (d->s), s);
}
gsl_test_rel (y[0], cos (tinit), d->n * tol,
"%s test_driver y0", gsl_odeiv2_step_name (d->s));
gsl_test_rel (y[1], sin (tinit), d->n * tol,
"%s test_driver y1", gsl_odeiv2_step_name (d->s));
gsl_odeiv2_driver_free (d);
/* Test that maximum number of steps is obeyed */
d = gsl_odeiv2_driver_alloc_y_new (&rhs_func_sin, T, 1e-3, tol, tol);
gsl_odeiv2_driver_set_hmax (d, hmax);
gsl_odeiv2_driver_set_nmax (d, maxsteps);
s = gsl_odeiv2_driver_apply (d, &t, t2, y);
if (d->n != maxsteps + 1)
{
gsl_test (1, "%s test_driver steps %d, expected %d",
gsl_odeiv2_step_name (d->s), d->n, maxsteps + 1);
}
else
{
gsl_test (0, "%s test_driver max steps test",
gsl_odeiv2_step_name (d->s));
}
gsl_odeiv2_driver_free (d);
/* Test that minimum step length is obeyed */
d = gsl_odeiv2_driver_alloc_y_new (&rhs_func_broken, T, 1e-3, tol, tol);
gsl_odeiv2_driver_set_hmin (d, hmin);
y[0] = 0.0;
t = 0.0;
s = gsl_odeiv2_driver_apply (d, &t, t2, y);
if (s != GSL_ENOPROG)
{
gsl_test (1, "%s test_driver min step test returned %d",
gsl_odeiv2_step_name (d->s), s);
}
else
{
gsl_test (0, "%s test_driver min step test",
gsl_odeiv2_step_name (d->s));
}
gsl_odeiv2_driver_free (d);
/* Test negative integration direction */
d = gsl_odeiv2_driver_alloc_y_new (&rhs_func_sin, T, -1e-3, tol, tol);
y[0] = 1.0;
y[1] = 0.0;
t = 2.5;
s = gsl_odeiv2_driver_apply (d, &t, t3, y);
{
const double tol = 1e-3;
const double test = fabs (t - t3);
const double val = fabs (sin (-5.0) - y[1]);
if (test > GSL_DBL_EPSILON)
{
gsl_test (1,
"%s test_driver negative dir diff %e, expected less than %e",
gsl_odeiv2_step_name (d->s), test, GSL_DBL_EPSILON);
}
else if (val > tol)
{
gsl_test (1,
"%s test_driver negative dir val %e, expected less than %e",
gsl_odeiv2_step_name (d->s), val, tol);
}
else
{
gsl_test (s, "%s test_driver negative direction test",
gsl_odeiv2_step_name (d->s));
}
}
/* Test driver_apply_fixed_step */
gsl_odeiv2_driver_reset_hstart (d, 1e-3);
s = gsl_odeiv2_driver_apply_fixed_step (d, &t, hfixed, nfsteps, y);
if (s != GSL_SUCCESS)
{
gsl_test (s, "%s test_driver apply_fixed_step returned %d",
gsl_odeiv2_step_name (d->s), s);
}
{
const double tol = 1e-3;
const double val = fabs (sin (-2.5) - y[1]);
if (fabs (t) > nfsteps * GSL_DBL_EPSILON)
{
gsl_test (1,
"%s test_driver apply_fixed_step t %e, expected less than %e",
gsl_odeiv2_step_name (d->s), fabs (t),
nfsteps * GSL_DBL_EPSILON);
}
else if (val > tol)
{
gsl_test (1,
"%s test_driver apply_fixed_step val %e, expected less than %e",
gsl_odeiv2_step_name (d->s), val, tol);
}
else
{
gsl_test (s, "%s test_driver apply_fixed_step test",
gsl_odeiv2_step_name (d->s));
}
}
gsl_odeiv2_driver_free (d);
}
void
benchmark_precision (void)
{
/* Tests steppers with several error tolerances and evaluates their
performance and precision.
*/
/* Maximum number of tolerances to be tested */
#define MAXNT 12
const gsl_odeiv2_step_type *steppers[MAXNS];
/* Required error tolerance for each stepper. */
double err_target[MAXNS];
/* initial values for each ode-solver */
double y[MAXEQ * MAXNS * MAXNT];
/* precise results from e.g. analytical solution */
double yres[MAXEQ];
size_t i, j, k, p;
/* Number of problems to test */
#define CONST_BENCHMARK_PRECISION_NPROB 3
/* Problems, their names and number of equations */
const gsl_odeiv2_system *prob[] =
{ &rhs_func_sin, &rhs_func_exp, &rhs_func_stiff };
const char *probname[] =
{ "rhs_func_sin", "rhs_func_exp", "rhs_func_stiff" };
const size_t sd[] = { 2, 1, 2 };
/* Integration interval for problems */
const double start[] = { 0.0, 0.0, 0.0 };
const double end[] = { 2.4, 8.4, 1.2 };
const double epsabs[] = { 1e-4, 1e-6, 1e-8, 1e-10, 1e-12, 1e-14, 0 };
const double epsrel[] = { 1e-4, 1e-6, 1e-8, 1e-10, 1e-12, 1e-14, 0 };
const double initstepsize = 1e-5;
/* number of function and jacobian evaluations */
extern int nfe, nje;
int nnfe[MAXNS * MAXNT] = { 0.0 };
int nnje[MAXNS * MAXNT] = { 0.0 };
/* Steppers */
steppers[0] = gsl_odeiv2_step_rk4;
err_target[0] = 1e-6;
steppers[1] = gsl_odeiv2_step_rk2;
err_target[1] = 1e-6;
steppers[2] = gsl_odeiv2_step_rkf45;
err_target[2] = 1e-6;
steppers[3] = gsl_odeiv2_step_rkck;
err_target[3] = 1e-6;
steppers[4] = gsl_odeiv2_step_rk8pd;
err_target[4] = 1e-6;
steppers[5] = gsl_odeiv2_step_rk1imp;
err_target[5] = 1e-3;
steppers[6] = gsl_odeiv2_step_rk2imp;
err_target[6] = 1e-5;
steppers[7] = gsl_odeiv2_step_rk4imp;
err_target[7] = 1e-6;
steppers[8] = gsl_odeiv2_step_bsimp;
err_target[8] = 1e-6;
steppers[9] = gsl_odeiv2_step_msadams;
err_target[9] = 1e-5;
steppers[10] = gsl_odeiv2_step_msbdf;
err_target[10] = 1e-5;
steppers[11] = 0;
/* Loop over problems */
for (p = 0; p < CONST_BENCHMARK_PRECISION_NPROB; p++)
{
/* Initialize */
for (i = 0; i < MAXNS * MAXEQ * MAXNT; i++)
{
y[i] = 0.0;
}
for (i = 0; i < MAXNS * MAXNT; i++)
{
switch (p)
{
case 0:
y[i * sd[p]] = 1.0;
y[i * sd[p] + 1] = 0.0;
yres[0] = cos (2.4);
yres[1] = sin (2.4);
break;
case 1:
y[i * sd[p]] = 1.0;
yres[0] = exp (8.4);
break;
case 2:
y[i * sd[p]] = 1.0;
y[i * sd[p] + 1] = 0.0;
yres[0] = 2 * exp (-1.2) - exp (-1000 * 1.2);
yres[1] = -exp (-1.2) + exp (-1000 * 1.2);
break;
default:
gsl_test (GSL_EFAILED,
"benchmark_precision: initialization error\n");
return;
}
}
/* Call each solver for the problem */
for (i = 0; steppers[i] != 0; i++)
for (j = 0; epsrel[j] != 0; j++)
{
int s = sys_driver (steppers[i], prob[p], start[p], end[p],
initstepsize, &y[sd[p] * (j + i * MAXNT)],
epsabs[j], epsrel[j], probname[p]);
if (s != GSL_SUCCESS)
{
for (k = 0; k < sd[p]; k++)
{
y[sd[p] * (j + i * MAXNT) + k] = GSL_NAN;
}
}
else
{
nnfe[j + i * MAXNT] = nfe;
nnje[j + i * MAXNT] = nje;
}
}
/* Print results */
printf
("benchmark_precision: diff = (y_true - y) / y_true with %s\n epsrel: ",
probname[p]);
for (i = 0; epsrel[i] != 0.0; i++)
{
printf ("%12.0e ", epsrel[i]);
}
printf ("\n");
for (i = 0; steppers[i] != 0; i++)
{
for (k = 0; k < sd[p]; k++)
{
printf ("%8s diff[%d]: ", steppers[i]->name, (int) k);
for (j = 0; epsrel[j] != 0; j++)
{
const double diff =
(yres[k] - y[sd[p] * (j + i * MAXNT) + k]) / yres[k];
printf ("%12.5e ", diff);
}
printf ("\n");
}
}
printf ("\n");
/* Print number of function/jacobian evaluations */
printf
("benchmark_precision: number of function/jacobian evaluations with %s\n epsrel: ",
probname[p]);
for (i = 0; epsrel[i] != 0.0; i++)
{
printf ("%12.0e ", epsrel[i]);
}
printf ("\n");
for (i = 0; steppers[i] != 0; i++)
{
printf ("%8s nfe/nje: ", steppers[i]->name);
for (j = 0; epsrel[j] != 0; j++)
{
printf ("%9d %9d | ", nnfe[j + i * MAXNT], nnje[j + i * MAXNT]);
}
printf ("\n");
}
printf ("\n");
}
}
void
test_evolve_temp (const gsl_odeiv2_step_type * T, double h, double err)
{
/* Temporary test */
double y[3];
double yfin[3];
y[0] = 1.0;
y[1] = 0.0;
y[2] = 0.0;
yfin[0] = 0;
yfin[1] = 0;
yfin[2] = 0;
test_evolve_system (T, &rhs_func_stiff, 0.0, 1.0, h, y, yfin, err, "temp");
}
/**********************************************************/
/* Main function */
/**********************************************************/
int
main (void)
{
/* Benchmark routine to compare stepper performance */
/* benchmark_precision(); return 0; */
/* Test single problem, for debugging purposes */
/* test_evolve_temp (gsl_odeiv2_step_msadams, 1e-3, 1e-8); return 0; */
int i;
struct ptype
{
const gsl_odeiv2_step_type *type;
double h;
}
p[MAXNS];
struct ptype explicit_stepper[MAXNS];
p[0].type = gsl_odeiv2_step_rk4;
p[0].h = 1.0e-3;
p[1].type = gsl_odeiv2_step_rk2;
p[1].h = 1.0e-3;
p[2].type = gsl_odeiv2_step_rkf45;
p[2].h = 1.0e-3;
p[3].type = gsl_odeiv2_step_rkck;
p[3].h = 1.0e-3;
p[4].type = gsl_odeiv2_step_rk8pd;
p[4].h = 1.0e-3;
p[5].type = gsl_odeiv2_step_rk1imp;
p[5].h = 1.0e-3;
p[6].type = gsl_odeiv2_step_rk2imp;
p[6].h = 1.0e-3;
p[7].type = gsl_odeiv2_step_rk4imp;
p[7].h = 1.0e-3;
p[8].type = gsl_odeiv2_step_bsimp;
p[8].h = 1.0e-3;
p[9].type = gsl_odeiv2_step_msadams;
p[9].h = 1.0e-3;
p[10].type = gsl_odeiv2_step_msbdf;
p[10].h = 1.0e-3;
p[11].type = 0;
gsl_ieee_env_setup ();
/* Basic tests for all steppers */
for (i = 0; p[i].type != 0; i++)
{
test_stepper (p[i].type);
test_evolve_linear (p[i].type, p[i].h, 1e-10);
test_evolve_exp (p[i].type, p[i].h, 1e-5);
test_evolve_sin (p[i].type, p[i].h, 1e-8);
test_evolve_xsin (p[i].type, p[i].h, 1e-8);
test_evolve_xsin (p[i].type, 0.1, 1e-8);
test_evolve_stiff1 (p[i].type, p[i].h, 1e-7);
test_evolve_stiff5 (p[i].type, p[i].h, 1e-7);
test_evolve_negative_h (p[i].type, p[i].h, 1e-7);
test_broken (p[i].type, p[i].h, 1e-8);
test_stepsize_fail (p[i].type, p[i].h);
test_user_break (p[i].type, p[i].h);
test_driver (p[i].type);
}
/* Derivative test for explicit steppers */
explicit_stepper[0].type = gsl_odeiv2_step_rk4;
explicit_stepper[0].h = 1.0e-3;
explicit_stepper[1].type = gsl_odeiv2_step_rk2;
explicit_stepper[1].h = 1.0e-3;
explicit_stepper[2].type = gsl_odeiv2_step_rkf45;
explicit_stepper[2].h = 1.0e-3;
explicit_stepper[3].type = gsl_odeiv2_step_rkck;
explicit_stepper[3].h = 1.0e-3;
explicit_stepper[4].type = gsl_odeiv2_step_rk8pd;
explicit_stepper[4].h = 1.0e-3;
explicit_stepper[5].type = gsl_odeiv2_step_msadams;
explicit_stepper[5].h = 1.0e-3;
explicit_stepper[6].type = 0;
for (i = 0; explicit_stepper[i].type != 0; i++)
{
test_stepfn (explicit_stepper[i].type);
test_stepfn2 (explicit_stepper[i].type);
}
/* Special tests */
test_nonstiff_problems ();
test_stiff_problems ();
test_extreme_problems ();
exit (gsl_test_summary ());
}