/* 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 #include #include #include #include #include #include #include #include #include #include #include #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 */ /* 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 ()); }