/***********************************************/ /* 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 #include // GSL_SUCCESS ... #include // cos(x), sin(x) ... #include // 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; }