diff options
Diffstat (limited to 'gsl-1.9/ode-initval/test.c')
| -rw-r--r-- | gsl-1.9/ode-initval/test.c | 1030 |
1 files changed, 1030 insertions, 0 deletions
diff --git a/gsl-1.9/ode-initval/test.c b/gsl-1.9/ode-initval/test.c new file mode 100644 index 0000000..7fb6001 --- /dev/null +++ b/gsl-1.9/ode-initval/test.c @@ -0,0 +1,1030 @@ +/* ode-initval/test_odeiv.c + * + * Copyright (C) 2004 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 2 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_odeiv.h> +#include "odeiv_util.h" + +/* Maximum number of ODE equations */ +#define MAXEQ 4 + +/* RHS for f=2. Solution y = 2 * t + t0 */ + +int +rhs_linear (double t, const double y[], double f[], void *params) +{ + f[0] = 2.0; + + return GSL_SUCCESS; +} + +int +jac_linear (double t, const double y[], double *dfdy, double dfdt[], + void *params) +{ + dfdy[0] = 0.0; + dfdt[0] = 0.0; + + return GSL_SUCCESS; +} + +gsl_odeiv_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) +{ + f[0] = y[0]; + + return GSL_SUCCESS; +} + +int +jac_exp (double t, const double y[], double *dfdy, double dfdt[], + void *params) +{ + dfdy[0] = y[0]; + dfdt[0] = 0.0; + + return GSL_SUCCESS; +} + +gsl_odeiv_system rhs_func_exp = { + rhs_exp, + jac_exp, + 1, + 0 +}; + +/* RHS for f0 = -y1, f1 = y0 + equals y = [cos(t), sin(t)] with initial values [1, 0] +*/ + +int +rhs_sin (double t, const double y[], double f[], void *params) +{ + 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) +{ + 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_odeiv_system rhs_func_sin = { + rhs_sin, + jac_sin, + 2, + 0 +}; + +/* Sine/cosine with random failures */ + +int +rhs_xsin (double t, const double y[], double f[], void *params) +{ + static int n = 0; + + n++; + + 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; + + 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_odeiv_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) +{ + 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) +{ + 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_odeiv_system rhs_func_stiff = { + rhs_stiff, + jac_stiff, + 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; + + 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; + + 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_odeiv_system rhs_func_vanderpol = { + rhs_vanderpol, + jac_vanderpol, + 2, + 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; + + 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; + + 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_odeiv_system rhs_func_oregonator = { + rhs_oregonator, + jac_oregonator, + 3, + 0 +}; + +/* Volterra-Lotka predator-prey model + + f0 = (a - b * y1) * y0 y0(0) = 3.0 + 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 = 1.0; + const double d = 1.0; + + 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 = 1.0; + const double d = 1.0; + + 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_odeiv_system rhs_func_vl = { + rhs_vl, + jac_vl, + 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) +{ + f[0] = -50 * (y[0] - cos(t)); + + return GSL_SUCCESS; +} + +int +jac_stifftrig (double t, const double y[], double *dfdy, double dfdt[], + void *params) +{ + dfdy[0] = -50; + + dfdt[0] = -50 * sin(t); + + return GSL_SUCCESS; +} + +gsl_odeiv_system rhs_func_stifftrig = { + rhs_stifftrig, + jac_stifftrig, + 1, + 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; + + 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; + + 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_odeiv_system rhs_func_e5 = { + rhs_e5, + jac_e5, + 4, + 0 +}; + +void +test_odeiv_stepper (const gsl_odeiv_step_type *T, const gsl_odeiv_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 with given values yfin + */ + + double y[MAXEQ] = {0.0}; + double yerr[MAXEQ] = {0.0}; + size_t ne = sys->dimension; + size_t i; + + gsl_odeiv_step *step = gsl_odeiv_step_alloc (T, ne); + + DBL_MEMCPY (y, ystart, MAXEQ); + + { + int s = gsl_odeiv_step_apply (step, 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_odeiv_step_name (step), desc,i); + } + + gsl_odeiv_step_free (step); +} + +void +test_stepper (const gsl_odeiv_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; + + /* 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); + + /* 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); +} + +void +test_evolve_system (const gsl_odeiv_step_type * T, + const gsl_odeiv_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_odeiv_step * step = gsl_odeiv_step_alloc (T, sys->dimension); + + gsl_odeiv_control *c = + gsl_odeiv_control_standard_new (err_target, err_target, 1.0, 0.0); + + gsl_odeiv_evolve *e = gsl_odeiv_evolve_alloc (sys->dimension); + + while (t < t1) + { + int s = gsl_odeiv_evolve_apply (e, c, step, sys, &t, t1, &h, y); + + if (s != GSL_SUCCESS && sys != &rhs_func_xsin) + { + gsl_test(s, "%s evolve_apply returned %d", + gsl_odeiv_step_name (step), s); + break; + } + + if (steps > 100000) + { + gsl_test(GSL_EFAILED, + "%s evolve_apply reached maxiter", + gsl_odeiv_step_name (step)); + break; + } + + steps++; + } + + /* 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 * e->count * err_target, + "%s %s evolve(%d)", + gsl_odeiv_step_name (step), desc, i); + } + + gsl_odeiv_evolve_free (e); + gsl_odeiv_control_free (c); + gsl_odeiv_step_free (step); +} + +int +sys_driver (const gsl_odeiv_step_type * T, + const gsl_odeiv_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_odeiv_step * step = gsl_odeiv_step_alloc (T, sys->dimension); + + gsl_odeiv_control *c = + gsl_odeiv_control_standard_new (epsabs, epsrel, 1.0, 0.0); + gsl_odeiv_evolve *e = gsl_odeiv_evolve_alloc (sys->dimension); + + while (t < t1) + { + s = gsl_odeiv_evolve_apply (e, c, step, sys, &t, t1, &h, y); + + if (s != GSL_SUCCESS) + { + gsl_test(s, "sys_driver: %s evolve_apply returned %d", + gsl_odeiv_step_name (step), s); + break; + } + + if (steps > 1e7) + { + gsl_test(GSL_EMAXITER, + "sys_driver: %s evolve_apply reached maxiter at t=%g", + gsl_odeiv_step_name (step), t); + s = GSL_EMAXITER; + break; + } + + steps++; + } + + gsl_test(s, "%s %s [%g,%g], %d steps completed", + gsl_odeiv_step_name (step), desc, t0, t1, steps); + + gsl_odeiv_evolve_free (e); + gsl_odeiv_control_free (c); + gsl_odeiv_step_free (step); + + return s; +} + +void +test_compare_vanderpol (void) +{ + /* Compares output of van Der Pol oscillator with several steppers */ + + /* system dimension */ + const size_t sd = 2; + + const gsl_odeiv_step_type *steppers[20]; + const gsl_odeiv_step_type **T; + + /* Required error tolerance for each stepper. */ + double err_target[20]; + + /* number of ODE solvers */ + const size_t ns = 11; + + /* initial values for each ode-solver */ + double y[11][2]; + double *yp = &y[0][0]; + + size_t i, j, k; + int status = 0; + + /* Parameters for the problem and stepper */ + const double start = 0.0; + const double end = 100.0; + const double epsabs = 1e-8; + const double epsrel = 1e-8; + const double initstepsize = 1e-5; + + /* Initialize */ + + steppers[0] = gsl_odeiv_step_rk2; + err_target[0] = 1e-6; + steppers[1] = gsl_odeiv_step_rk4; + err_target[1] = 1e-6; + steppers[2] = gsl_odeiv_step_rkf45; + err_target[2] = 1e-6; + steppers[3] = gsl_odeiv_step_rkck; + err_target[3] = 1e-6; + steppers[4] = gsl_odeiv_step_rk8pd; + err_target[4] = 1e-6; + steppers[5] = gsl_odeiv_step_rk2imp; + err_target[5] = 1e-5; + steppers[6] = gsl_odeiv_step_rk2simp; + err_target[6] = 1e-5; + steppers[7] = gsl_odeiv_step_rk4imp; + err_target[7] = 1e-6; + steppers[8] = gsl_odeiv_step_bsimp; + err_target[8] = 1e-7; + steppers[9] = gsl_odeiv_step_gear1; + err_target[9] = 1e-2; + steppers[10] = gsl_odeiv_step_gear2; + err_target[10] = 1e-6; + steppers[11] = 0; + + T = steppers; + + for (i = 0; i < ns; i++) + { + y[i][0] = 1.0; + y[i][1] = 0.0; + } + + /* Call each solver for the problem */ + + i = 0; + while (*T != 0) + { + { + int s = sys_driver (*T, &rhs_func_vanderpol, + start, end, initstepsize, &yp[i], + epsabs, epsrel, "vanderpol"); + if (s != GSL_SUCCESS) + { + status++; + } + } + + T++; + i += sd; + } + + if (status != GSL_SUCCESS) + { + return; + } + + /* Compare results */ + + T = steppers; + + for (i = 0; i < ns; i++) + for (j = i+1; j < ns; j++) + for (k = 0; k < sd; k++) + { + const double val1 = yp[sd * i + k]; + const double val2 = yp[sd * j + k]; + gsl_test_abs (val1, val2, + ( GSL_MAX(err_target[i], err_target[j]) ), + "%s/%s vanderpol", + T[i]->name, T[j]->name); + } + +} + +void +test_compare_oregonator (void) +{ + /* Compares output of the Oregonator with several steppers */ + + /* system dimension */ + const size_t sd = 3; + + const gsl_odeiv_step_type *steppers[20]; + const gsl_odeiv_step_type **T; + + /* Required error tolerance for each stepper. */ + double err_target[20]; + + /* number of ODE solvers */ + const size_t ns = 2; + + /* initial values for each ode-solver */ + double y[2][3]; + double *yp = &y[0][0]; + + size_t i, j, k; + int status = 0; + + /* Parameters for the problem and stepper */ + const double start = 0.0; + const double end = 360.0; + const double epsabs = 1e-8; + const double epsrel = 1e-8; + const double initstepsize = 1e-5; + + /* Initialize */ + + steppers[0] = gsl_odeiv_step_rk2simp; + err_target[0] = 1e-6; + steppers[1] = gsl_odeiv_step_bsimp; + err_target[1] = 1e-6; + steppers[2] = 0; + + T = steppers; + + for (i = 0; i < ns; i++) + { + y[i][0] = 1.0; + y[i][1] = 2.0; + y[i][2] = 3.0; + } + + /* Call each solver for the problem */ + + i = 0; + while (*T != 0) + { + { + int s = sys_driver (*T, &rhs_func_oregonator, + start, end, initstepsize, &yp[i], + epsabs, epsrel, "oregonator"); + + if (s != GSL_SUCCESS) + { + status++; + } + } + + T++; + i += sd; + } + + if (status != GSL_SUCCESS) + { + return; + } + + /* Compare results */ + + T = steppers; + + for (i = 0; i < ns; i++) + for (j = i+1; j < ns; j++) + for (k = 0; k < sd; k++) + { + const double val1 = yp[sd * i + k]; + const double val2 = yp[sd * j + k]; + gsl_test_rel (val1, val2, + ( GSL_MAX(err_target[i], err_target[j]) ), + "%s/%s oregonator", + T[i]->name, T[j]->name); + } + +} + +void +test_evolve_linear (const gsl_odeiv_step_type * T, double h, double err) +{ + 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_odeiv_step_type * T, double h, double err) +{ + 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_odeiv_step_type * T, double h, double err) +{ + 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_odeiv_step_type * T, double h, double err) +{ + 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_xsin, 0.0, 2.0, h, y, yfin, err, + "sine[0,2] w/errors"); +} + + +void +test_evolve_stiff1 (const gsl_odeiv_step_type * T, double h, double err) +{ + 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_odeiv_step_type * T, double h, double err) +{ + 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]"); +} + + +int +main (void) +{ + int i; + + struct ptype + { + const gsl_odeiv_step_type *type; + double h; + } + p[20]; + + p[0].type = gsl_odeiv_step_rk2; + p[0].h = 1.0e-3; + p[1].type = gsl_odeiv_step_rk4; + p[1].h = 1.0e-3; + p[2].type = gsl_odeiv_step_rkf45; + p[2].h = 1.0e-3; + p[3].type = gsl_odeiv_step_rkck; + p[3].h = 1.0e-3; + p[4].type = gsl_odeiv_step_rk8pd; + p[4].h = 1.0e-3; + p[5].type = gsl_odeiv_step_rk2imp; + p[5].h = 1.0e-3; + p[6].type = gsl_odeiv_step_rk2simp; + p[6].h = 1.0e-3; + p[7].type = gsl_odeiv_step_rk4imp; + p[7].h = 1.0e-3; + p[8].type = gsl_odeiv_step_bsimp; + p[8].h = 1.0e-3; + p[9].type = gsl_odeiv_step_gear1; + p[9].h = 1.0e-3; + p[10].type = gsl_odeiv_step_gear2; + p[10].h = 1.0e-3; + p[11].type = 0; + + gsl_ieee_env_setup (); + + for (i = 0; p[i].type != 0; i++) + { + test_stepper(p[i].type); + } + + for (i = 0; p[i].type != 0; i++) + { + test_evolve_linear (p[i].type, p[i].h, 1e-10); + test_evolve_exp (p[i].type, p[i].h, 1e-6); + test_evolve_sin (p[i].type, p[i].h, 1e-8); + test_evolve_xsin (p[i].type, p[i].h, 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_compare_vanderpol(); + test_compare_oregonator(); + + exit (gsl_test_summary ()); +} |