From 197e88c3b56d28840217010a2871c6ea3a4dd1a4 Mon Sep 17 00:00:00 2001 From: Alberto Ruiz Date: Wed, 21 May 2014 10:30:55 +0200 Subject: update dependencies, move examples etc --- packages/gsl/src/Numeric/GSL/gsl-ode.c | 182 +++++++++++++++++++++++++++++++++ 1 file changed, 182 insertions(+) create mode 100644 packages/gsl/src/Numeric/GSL/gsl-ode.c (limited to 'packages/gsl/src/Numeric/GSL/gsl-ode.c') diff --git a/packages/gsl/src/Numeric/GSL/gsl-ode.c b/packages/gsl/src/Numeric/GSL/gsl-ode.c new file mode 100644 index 0000000..3f2771b --- /dev/null +++ b/packages/gsl/src/Numeric/GSL/gsl-ode.c @@ -0,0 +1,182 @@ + +#ifdef GSLODE1 + +////////////////////////////// ODE V1 ////////////////////////////////////////// + +#include + +typedef struct {int n; int (*f)(double,int, const double*, int, double *); int (*j)(double,int, const double*, int, int, double*);} Tode; + +int odefunc (double t, const double y[], double f[], void *params) { + Tode * P = (Tode*) params; + (P->f)(t,P->n,y,P->n,f); + return GSL_SUCCESS; +} + +int odejac (double t, const double y[], double *dfdy, double dfdt[], void *params) { + Tode * P = ((Tode*) params); + (P->j)(t,P->n,y,P->n,P->n,dfdy); + int j; + for (j=0; j< P->n; j++) + dfdt[j] = 0.0; + return GSL_SUCCESS; +} + + +int ode(int method, double h, double eps_abs, double eps_rel, + int f(double, int, const double*, int, double*), + int jac(double, int, const double*, int, int, double*), + KRVEC(xi), KRVEC(ts), RMAT(sol)) { + + const gsl_odeiv_step_type * T; + + switch(method) { + case 0 : {T = gsl_odeiv_step_rk2; break; } + case 1 : {T = gsl_odeiv_step_rk4; break; } + case 2 : {T = gsl_odeiv_step_rkf45; break; } + case 3 : {T = gsl_odeiv_step_rkck; break; } + case 4 : {T = gsl_odeiv_step_rk8pd; break; } + case 5 : {T = gsl_odeiv_step_rk2imp; break; } + case 6 : {T = gsl_odeiv_step_rk4imp; break; } + case 7 : {T = gsl_odeiv_step_bsimp; break; } + case 8 : { printf("Sorry: ODE rk1imp not available in this GSL version\n"); exit(0); } + case 9 : { printf("Sorry: ODE msadams not available in this GSL version\n"); exit(0); } + case 10: { printf("Sorry: ODE msbdf not available in this GSL version\n"); exit(0); } + default: ERROR(BAD_CODE); + } + + gsl_odeiv_step * s = gsl_odeiv_step_alloc (T, xin); + gsl_odeiv_control * c = gsl_odeiv_control_y_new (eps_abs, eps_rel); + gsl_odeiv_evolve * e = gsl_odeiv_evolve_alloc (xin); + + Tode P; + P.f = f; + P.j = jac; + P.n = xin; + + gsl_odeiv_system sys = {odefunc, odejac, xin, &P}; + + double t = tsp[0]; + + double* y = (double*)calloc(xin,sizeof(double)); + int i,j; + for(i=0; i< xin; i++) { + y[i] = xip[i]; + solp[i] = xip[i]; + } + + for (i = 1; i < tsn ; i++) + { + double ti = tsp[i]; + while (t < ti) + { + gsl_odeiv_evolve_apply (e, c, s, + &sys, + &t, ti, &h, + y); + // if (h < hmin) h = hmin; + } + for(j=0; j + +typedef struct {int n; int (*f)(double,int, const double*, int, double *); int (*j)(double,int, const double*, int, int, double*);} Tode; + +int odefunc (double t, const double y[], double f[], void *params) { + Tode * P = (Tode*) params; + (P->f)(t,P->n,y,P->n,f); + return GSL_SUCCESS; +} + +int odejac (double t, const double y[], double *dfdy, double dfdt[], void *params) { + Tode * P = ((Tode*) params); + (P->j)(t,P->n,y,P->n,P->n,dfdy); + int j; + for (j=0; j< P->n; j++) + dfdt[j] = 0.0; + return GSL_SUCCESS; +} + + +int ode(int method, double h, double eps_abs, double eps_rel, + int f(double, int, const double*, int, double*), + int jac(double, int, const double*, int, int, double*), + KRVEC(xi), KRVEC(ts), RMAT(sol)) { + + const gsl_odeiv2_step_type * T; + + switch(method) { + case 0 : {T = gsl_odeiv2_step_rk2; break; } + case 1 : {T = gsl_odeiv2_step_rk4; break; } + case 2 : {T = gsl_odeiv2_step_rkf45; break; } + case 3 : {T = gsl_odeiv2_step_rkck; break; } + case 4 : {T = gsl_odeiv2_step_rk8pd; break; } + case 5 : {T = gsl_odeiv2_step_rk2imp; break; } + case 6 : {T = gsl_odeiv2_step_rk4imp; break; } + case 7 : {T = gsl_odeiv2_step_bsimp; break; } + case 8 : {T = gsl_odeiv2_step_rk1imp; break; } + case 9 : {T = gsl_odeiv2_step_msadams; break; } + case 10: {T = gsl_odeiv2_step_msbdf; break; } + default: ERROR(BAD_CODE); + } + + Tode P; + P.f = f; + P.j = jac; + P.n = xin; + + gsl_odeiv2_system sys = {odefunc, odejac, xin, &P}; + + gsl_odeiv2_driver * d = + gsl_odeiv2_driver_alloc_y_new (&sys, T, h, eps_abs, eps_rel); + + double t = tsp[0]; + + double* y = (double*)calloc(xin,sizeof(double)); + int i,j; + int status=0; + for(i=0; i< xin; i++) { + y[i] = xip[i]; + solp[i] = xip[i]; + } + + for (i = 1; i < tsn ; i++) + { + double ti = tsp[i]; + + status = gsl_odeiv2_driver_apply (d, &t, ti, y); + + if (status != GSL_SUCCESS) { + printf ("error in ode, return value=%d\n", status); + break; + } + +// printf ("%.5e %.5e %.5e\n", t, y[0], y[1]); + + for(j=0; j