diff options
author | Alberto Ruiz <aruiz@um.es> | 2010-02-07 09:29:50 +0000 |
---|---|---|
committer | Alberto Ruiz <aruiz@um.es> | 2010-02-07 09:29:50 +0000 |
commit | aef0333b5180ea79e539bd53194f1dfed20b7db5 (patch) | |
tree | 6ece2434ecacab194331120bd47d09ab04ae0f4f /lib/Numeric/GSL/gsl-aux.c | |
parent | 634683fcfab73a0bd830fef03fb9f4603ba837b6 (diff) |
added odeSolve
Diffstat (limited to 'lib/Numeric/GSL/gsl-aux.c')
-rw-r--r-- | lib/Numeric/GSL/gsl-aux.c | 89 |
1 files changed, 88 insertions, 1 deletions
diff --git a/lib/Numeric/GSL/gsl-aux.c b/lib/Numeric/GSL/gsl-aux.c index 0c71ca1..0fcde1c 100644 --- a/lib/Numeric/GSL/gsl-aux.c +++ b/lib/Numeric/GSL/gsl-aux.c | |||
@@ -22,6 +22,7 @@ | |||
22 | #include <gsl/gsl_complex_math.h> | 22 | #include <gsl/gsl_complex_math.h> |
23 | #include <gsl/gsl_rng.h> | 23 | #include <gsl/gsl_rng.h> |
24 | #include <gsl/gsl_randist.h> | 24 | #include <gsl/gsl_randist.h> |
25 | #include <gsl/gsl_odeiv.h> | ||
25 | #include <string.h> | 26 | #include <string.h> |
26 | #include <stdio.h> | 27 | #include <stdio.h> |
27 | 28 | ||
@@ -507,7 +508,7 @@ double f_aux_min(const gsl_vector*x, void *pars) { | |||
507 | 508 | ||
508 | 509 | ||
509 | void df_aux_min(const gsl_vector * x, void * pars, gsl_vector * g) { | 510 | void df_aux_min(const gsl_vector * x, void * pars, gsl_vector * g) { |
510 | Tfdf * fdf = ((Tfdf*) pars); | 511 | Tfdf * fdf = ((Tfdf*) pars); |
511 | double* p = (double*)calloc(x->size,sizeof(double)); | 512 | double* p = (double*)calloc(x->size,sizeof(double)); |
512 | double* q = (double*)calloc(g->size,sizeof(double)); | 513 | double* q = (double*)calloc(g->size,sizeof(double)); |
513 | int k; | 514 | int k; |
@@ -797,3 +798,89 @@ int random_vector(int seed, int code, RVEC(r)) { | |||
797 | } | 798 | } |
798 | } | 799 | } |
799 | #undef RAN | 800 | #undef RAN |
801 | |||
802 | ////////////////////////////////////////////////////// | ||
803 | // ODE | ||
804 | |||
805 | typedef struct {int n; int (*f)(double,int, const double*, int, double *); int (*j)(double,int, const double*, int, int, double*);} Tode; | ||
806 | |||
807 | int odefunc (double t, const double y[], double f[], void *params) { | ||
808 | Tode * P = (Tode*) params; | ||
809 | (P->f)(t,P->n,y,P->n,f); | ||
810 | return GSL_SUCCESS; | ||
811 | } | ||
812 | |||
813 | int odejac (double t, const double y[], double *dfdy, double dfdt[], void *params) { | ||
814 | Tode * P = ((Tode*) params); | ||
815 | (P->j)(t,P->n,y,P->n,P->n,dfdy); | ||
816 | int j; | ||
817 | for (j=0; j< P->n; j++) | ||
818 | dfdt[j] = 0.0; | ||
819 | return GSL_SUCCESS; | ||
820 | } | ||
821 | |||
822 | |||
823 | int ode(int method, double h, double eps_abs, double eps_rel, | ||
824 | int f(double, int, const double*, int, double*), | ||
825 | int jac(double, int, const double*, int, int, double*), | ||
826 | KRVEC(xi), KRVEC(ts), RMAT(sol)) { | ||
827 | |||
828 | const gsl_odeiv_step_type * T; | ||
829 | |||
830 | switch(method) { | ||
831 | case 0 : {T = gsl_odeiv_step_rk2; break; } | ||
832 | case 1 : {T = gsl_odeiv_step_rk4; break; } | ||
833 | case 2 : {T = gsl_odeiv_step_rkf45; break; } | ||
834 | case 3 : {T = gsl_odeiv_step_rkck; break; } | ||
835 | case 4 : {T = gsl_odeiv_step_rk8pd; break; } | ||
836 | case 5 : {T = gsl_odeiv_step_rk2imp; break; } | ||
837 | case 6 : {T = gsl_odeiv_step_rk4imp; break; } | ||
838 | case 7 : {T = gsl_odeiv_step_bsimp; break; } | ||
839 | case 8 : {T = gsl_odeiv_step_gear1; break; } | ||
840 | case 9 : {T = gsl_odeiv_step_gear2; break; } | ||
841 | default: ERROR(BAD_CODE); | ||
842 | } | ||
843 | |||
844 | |||
845 | gsl_odeiv_step * s = gsl_odeiv_step_alloc (T, xin); | ||
846 | gsl_odeiv_control * c = gsl_odeiv_control_y_new (eps_abs, eps_rel); | ||
847 | gsl_odeiv_evolve * e = gsl_odeiv_evolve_alloc (xin); | ||
848 | |||
849 | Tode P; | ||
850 | P.f = f; | ||
851 | P.j = jac; | ||
852 | P.n = xin; | ||
853 | |||
854 | gsl_odeiv_system sys = {odefunc, odejac, xin, &P}; | ||
855 | |||
856 | double t = tsp[0]; | ||
857 | |||
858 | double* y = (double*)calloc(xin,sizeof(double)); | ||
859 | int i,j; | ||
860 | for(i=0; i< xin; i++) { | ||
861 | y[i] = xip[i]; | ||
862 | solp[i] = xip[i]; | ||
863 | } | ||
864 | |||
865 | for (i = 1; i < tsn ; i++) | ||
866 | { | ||
867 | double ti = tsp[i]; | ||
868 | while (t < ti) | ||
869 | { | ||
870 | gsl_odeiv_evolve_apply (e, c, s, | ||
871 | &sys, | ||
872 | &t, ti, &h, | ||
873 | y); | ||
874 | // if (h < hmin) h = hmin; | ||
875 | } | ||
876 | for(j=0; j<xin; j++) { | ||
877 | solp[i*xin + j] = y[j]; | ||
878 | } | ||
879 | } | ||
880 | |||
881 | free(y); | ||
882 | gsl_odeiv_evolve_free (e); | ||
883 | gsl_odeiv_control_free (c); | ||
884 | gsl_odeiv_step_free (s); | ||
885 | return 0; | ||
886 | } | ||