Last Update: Oct. 4, 2007
/***********************************************/
/* T.Kouya's GSL sample program collection */
/* Stiff and Non-stiff Linear Ordinary */
/* Differential Equations */
/* Written by Tomonori Kouya */
/* */
/* Version 0.01: 2007-08-13 */
/***********************************************/
#include <stdio.h>
#include <gsl/gsl_errno.h>// GSL_SUCCESS ...
#include <gsl/gsl_math.h> // cos(x), sin(x) ...
#include <gsl/gsl_odeiv.h>// ODE solver
/* Return relative and absolute errors */
void vec_rel_error(double relerr[], double approx_val[], double true_val[], int dimension)
{
static int i;
static double abs_error, rel_error;
for(i = 0; i < dimension; i++)
{
abs_error = fabs(approx_val[i] - true_val[i]);
rel_error = abs_error;
if(true_val[i] >= 1.0e-15)
rel_error /= fabs(true_val[i]);
relerr[i] = rel_error;
}
return;
}
/* Dimension of ODEs */
#define DIM 2
/* Common True Solution */
void true_solution(double ret_y[], double x)
{
ret_y[0] = exp(-x);
ret_y[1] = exp(-x) + cos(x);
}
/* Definition of ODE system */
/* Non-stiff Linear ODE (cf. http://na-inet.jp/nasoft/chap16.pdf) */
/* dydx[] := func(x, y[]) */
int nonstiff_lode_func(double x, const double y[], double dydx[], void *params)
{
dydx[0] = -2.0 * y[0] + y[1] - cos(x);
dydx[1] = 2.0 * y[0] - 3.0 * y[1] + 3.0 * cos(x) - sin(x);
return GSL_SUCCESS;
}
/* Jacobian matrix of nonstiff_lode_func */
int jac_nonstiff_lode_func(double x, const double y[], double *dfdy, double dfdx[], void *params)
{
/* Jacobian matrix -> dfdy */
*(dfdy + 0 * DIM + 0) = -2.0;
*(dfdy + 0 * DIM + 1) = 1.0;
*(dfdy + 1 * DIM + 0) = 2.0;
*(dfdy + 1 * DIM + 1) = -3.0;
/* df/dt */
dfdx[0] = sin(x);
dfdx[1] = -3.0 * sin(x) - cos(x);
return GSL_SUCCESS;
}
/* Definition of ODE system */
/* Stiff Linear ODE (cf. http://na-inet.jp/nasoft/chap16.pdf) */
/* dydx[] := func(x, y[]) */
int stiff_lode_func(double x, const double y[], double dydx[], void *params)
{
dydx[0] = -2.0 * y[0] + y[1] - cos(x);
dydx[1] = 1998.0 * y[0] - 1999.0 * y[1] + 1999.0 * cos(x) - sin(x);
return GSL_SUCCESS;
}
/* Jacobian matrix of stiff_lode_func */
int jac_stiff_lode_func(double x, const double y[], double *dfdy, double dfdx[], void *params)
{
/* Jacobian matrix -> dfdy */
*(dfdy + 0 * DIM + 0) = -2.0;
*(dfdy + 0 * DIM + 1) = 1.0;
*(dfdy + 1 * DIM + 0) = 1998.0;
*(dfdy + 1 * DIM + 1) = -1999.0;
/* df/dt */
dfdx[0] = sin(x);
dfdx[1] = -1999.0 * sin(x) - cos(x);
return GSL_SUCCESS;
}
int main(void)
{
int i;
/* Integration Interval: [0, 10] */
/* Initial min stepsize: 1.0e-5 */
/* Initial value : y(0) = [1, 2]^T */
double x_start = 0.0, x_end = 10.0;
double h_init = 1.0e-5;
double y_init[DIM] = {1.0, 2.0};
double relerr_y[DIM], true_y[DIM];
/* Variables used for evolving */
int status_nonstiff, status_stiff;
double h_nonstiff, x_nonstiff, y_nonstiff[DIM];
double h_stiff, x_stiff, y_stiff[DIM];
/* Definitions to determine ODE solvers */
const gsl_odeiv_step_type *solver_nonstiff = gsl_odeiv_step_rkf45; // Runge-Kutta Felberg (4, 5)
const gsl_odeiv_step_type *solver_stiff = gsl_odeiv_step_rk4imp; // Fully implicit RK Gauss 4th order
/* Memories of stepsizes */
gsl_odeiv_step *step_nonstiff, *step_stiff;
/* Constants to control errors */
gsl_odeiv_control *tol_nonstiff, *tol_stiff;
/* Memories needed to evolve ODE solvers */
gsl_odeiv_evolve *evol_nonstiff, *evol_stiff;
/* Definitions of ODE systems */
gsl_odeiv_system sys_nonstiff, sys_stiff;
/* Non-stiff Problem */
/* ODE system */
sys_nonstiff.function = nonstiff_lode_func;
sys_nonstiff.jacobian = jac_nonstiff_lode_func;
sys_nonstiff.dimension = (size_t)(DIM);
sys_nonstiff.params = NULL;
/* ODE solver */
/* Determine constans for error controling */
/* Preparing evolution */
step_nonstiff = gsl_odeiv_step_alloc(solver_nonstiff, DIM);
tol_nonstiff = gsl_odeiv_control_standard_new(1.0e-10, 1.0e-5, 1.0, 0.0);
evol_nonstiff = gsl_odeiv_evolve_alloc(DIM);
/* Initialize for integration */
h_nonstiff = h_init;
x_nonstiff = x_start;
for(i = 0; i < DIM; i++)
y_nonstiff[i] = y_init[i];
/* Integration with stepsize control */
printf("Non-stiff Problem:\n");
printf(" Solver: %s, Controling: %s\n", gsl_odeiv_step_name(step_nonstiff), gsl_odeiv_control_name(tol_nonstiff));
printf(" x y[0] y[1] Relerr y[0] Relerr y[1]\n");
while(x_nonstiff < x_end)
{
status_nonstiff = gsl_odeiv_evolve_apply(evol_nonstiff, tol_nonstiff, step_nonstiff, &sys_nonstiff, &x_nonstiff, x_end, &h_nonstiff, y_nonstiff);
if(status_nonstiff != GSL_SUCCESS)
{
fprintf(stderr, "ERROR: Non-stiff Problem at x = %25.17e\n", x_nonstiff);
break;
}
true_solution(true_y, x_nonstiff);
vec_rel_error(relerr_y, y_nonstiff, true_y, DIM);
printf("%15.7e %25.17e %25.17e %10.3e %10.3e\n", x_nonstiff, y_nonstiff[0], y_nonstiff[1], relerr_y[0], relerr_y[1]);
}
/* Free */
gsl_odeiv_evolve_free(evol_nonstiff);
gsl_odeiv_control_free(tol_nonstiff);
gsl_odeiv_step_free(step_nonstiff);
/* Stiff Problem */
/* ODE system */
sys_stiff.function = stiff_lode_func;
sys_stiff.jacobian = jac_stiff_lode_func;
sys_stiff.dimension = (size_t)(DIM);
sys_stiff.params = NULL;
/* ODE solver */
/* Determine constans for error controling */
/* Preparing evolution */
step_stiff = gsl_odeiv_step_alloc(solver_stiff, DIM);
tol_stiff = gsl_odeiv_control_standard_new(1.0e-10, 1.0e-5, 1.0, 0.0);
evol_stiff = gsl_odeiv_evolve_alloc(DIM);
/* Initialize for integration */
h_stiff = h_init;
x_stiff = x_start;
for(i = 0; i < DIM; i++)
y_stiff[i] = y_init[i];
/* Integration with stepsize control */
printf("Stiff Problem:\n");
printf(" Solver: %s, Controling: %s\n", gsl_odeiv_step_name(step_stiff), gsl_odeiv_control_name(tol_stiff));
printf(" x y[0] y[1] Relerr y[0] Relerr y[1]\n");
while(x_stiff < x_end)
{
status_stiff = gsl_odeiv_evolve_apply(evol_stiff, tol_stiff, step_stiff, &sys_stiff, &x_stiff, x_end, &h_stiff, y_stiff);
if(status_stiff != GSL_SUCCESS)
{
fprintf(stderr, "ERROR: Stiff Problem at x = %25.17e\n", x_stiff);
break;
}
true_solution(true_y, x_stiff);
vec_rel_error(relerr_y, y_stiff, true_y, DIM);
printf("%15.7e %25.17e %25.17e %10.3e %10.3e\n", x_stiff, y_stiff[0], y_stiff[1], relerr_y[0], relerr_y[1]);
}
/* Free */
gsl_odeiv_evolve_free(evol_stiff);
gsl_odeiv_control_free(tol_stiff);
gsl_odeiv_step_free(step_stiff);
return 0;
}
$ ./linear_ode
Non-stiff Problem:
Solver: rkf45, Controling: standard
x y[0] y[1] Relerr y[0] Relerr y[1]
1.0000000e-05 9.99990000049999828e-01 1.99998999999999993e+00 0.000e+00 0.000e+00
6.0000000e-05 9.99940001799964007e-01 1.99993999999996408e+00 0.000e+00 0.000e+00
3.1000000e-04 9.99690048045035251e-01 1.99968999999503572e+00 0.000e+00 0.000e+00
1.5600000e-03 9.98441216167510692e-01 1.99843999936775751e+00 0.000e+00 0.000e+00
7.8100000e-03 9.92220418808186233e-01 1.99218992091321856e+00 3.692e-15 3.567e-15
3.9060000e-02 9.61693005778894716e-01 1.96093026112794799e+00 5.827e-11 5.610e-11
(OMIT)
9.7572849e+00 5.78719406467666030e-05 -9.45169137882558519e-01 7.013e-06 8.551e-10
9.7874581e+00 5.61518489106037504e-05 -9.34893173693603030e-01 7.137e-06 8.440e-10
9.8176313e+00 5.44828825914175032e-05 -9.23766026470360435e-01 7.265e-06 8.332e-10
9.8478044e+00 5.28635220633079945e-05 -9.11797828859521320e-01 7.396e-06 8.227e-10
9.8779776e+00 5.12922928725537602e-05 -8.98999479068398610e-01 7.529e-06 8.122e-10
9.9081507e+00 4.97677643944452827e-05 -8.85382630946418447e-01 7.664e-06 8.018e-10
9.9383239e+00 4.82885485301170624e-05 -8.70959683378632055e-01 7.798e-06 7.913e-10
9.9684971e+00 4.68532984421406228e-05 -8.55743769000906318e-01 7.932e-06 7.807e-10
9.9986702e+00 4.54607073277795066e-05 -8.39748742247069080e-01 8.065e-06 7.699e-10
1.0000000e+01 4.54002943932760116e-05 -8.39026129912536334e-01 8.032e-06 7.658e-10
Stiff Problem:
Solver: rk4imp, Controling: standard
x y[0] y[1] Relerr y[0] Relerr y[1]
1.0000000e-05 9.99990000049999828e-01 1.99998999999999971e+00 0.000e+00 1.110e-16
6.0000000e-05 9.99940001799964007e-01 1.99993999999996386e+00 0.000e+00 1.110e-16
3.1000000e-04 9.99690048045035362e-01 1.99968999999464692e+00 1.111e-16 1.944e-13
1.5600000e-03 9.98441216170407153e-01 1.99843999358035185e+00 2.901e-12 2.896e-09
4.3639627e-03 9.95645547812364029e-01 1.99563145641214712e+00 2.296e-09 2.289e-06
5.8180308e-03 9.94198861343948082e-01 1.99418160531753896e+00 1.667e-10 1.661e-07
8.0221772e-03 9.92009914828045125e-01 1.99197728977807276e+00 2.257e-10 2.246e-07
(OMIT)
9.9804552e+00 4.62958761739525062e-05 -8.49496859624095046e-01 2.480e-06 2.294e-07
9.9825539e+00 4.61988174664507395e-05 -8.48387990757782595e-01 2.482e-06 2.291e-07
9.9846526e+00 4.61019622419043659e-05 -8.47275384776054685e-01 2.484e-06 2.288e-07
9.9867512e+00 4.60053100737143624e-05 -8.46159046580215457e-01 2.486e-06 2.285e-07
9.9888499e+00 4.59088605361768437e-05 -8.45038981088006569e-01 2.488e-06 2.282e-07
9.9909486e+00 4.58126132044794909e-05 -8.43915193233583549e-01 2.490e-06 2.279e-07
9.9930473e+00 4.57165676547016873e-05 -8.42787687967497146e-01 2.492e-06 2.276e-07
9.9951460e+00 4.56207234638109064e-05 -8.41656470256668010e-01 2.494e-06 2.273e-07
9.9972447e+00 4.55250802096614788e-05 -8.40521545084366828e-01 2.496e-06 2.270e-07
9.9993434e+00 4.54296374709929523e-05 -8.39382917450192223e-01 2.498e-06 2.267e-07
1.0000000e+01 4.53998991727971499e-05 -8.39026068028467908e-01 6.738e-07 6.112e-08