Return to call_gsl_ode.c CVS log

Up to [local] / OpenXM / src / ox_gsl

Annotation of OpenXM/src/ox_gsl/call_gsl_ode.c, Revision 1.1

1.1     ! takayama    1: /* $OpenXM$
        !             2: */
        !             3: //#include <gsl/gsl_types.h>
        !             4: //#include <gsl/gsl_sys.h>
        !             5: #include <unistd.h>
        !             6: #include <gsl/gsl_sf_result.h>
        !             7: #include <gsl/gsl_errno.h>
        !             8: #include <gsl/gsl_matrix.h>
        !             9: #include <gsl/gsl_odeiv.h>
        !            10: #include "ox_gsl.h"
        !            11: extern int Debug;
        !            12: // local prototype declarations
        !            13:
        !            14: int GSL_ode_dim=1;
        !            15: #define VEC_MAX 30
        !            16: cmo *GSL_ode_func[VEC_MAX];
        !            17: int f_ode(double x, const double y[], double f[], void *params)
        !            18: {
        !            19:   int debug_f_ode=0;
        !            20:   int dim;
        !            21:   int i;
        !            22:   double d;
        !            23:   char *yy[VEC_MAX]={"y0","y1","y2","y3","y4","y5","y6","y7","y8","y9",
        !            24:                 "y10","y11","y12","y13","y14","y15","y16","y17","y18","y19",
        !            25:                 "y20","y21","y22","y23","y24","y25","y26","y27","y28","y29"
        !            26:                };
        !            27:   /* result is stored in f */
        !            28:   dim = GSL_ode_dim;
        !            29:   if (dim > VEC_MAX) GSL_ERROR("f_n supports functions with args <= VEC_MAX",GSL_ETOL);
        !            30:   init_dic();
        !            31:   register_entry("x",x);
        !            32:   for (i=0; i<dim; i++) {
        !            33:     if (debug_f_ode) ox_printf("y%d=%lg, ",i,y[i]);
        !            34:     register_entry(yy[i],y[i]);
        !            35:   }
        !            36:   if (debug_f_ode) { ox_printf("\n"); fflush(NULL); }
        !            37:   for (i=0; i<dim; i++) {
        !            38:     if (eval_cmo(GSL_ode_func[i],&d)==0) GSL_ERROR("eval_cmo fails in f_ode",GSL_ETOL);
        !            39:     if (debug_f_ode) ox_printf("\nGSL_ode_func[%d](...) -> d=%lg\n",i,d);
        !            40:     f[i] = d;
        !            41:   }
        !            42:   return GSL_SUCCESS;
        !            43: }
        !            44:
        !            45: void  call_gsl_odeiv_step(char *type) {
        !            46:   int argc;
        !            47:   int dim;
        !            48:   cmo *f;
        !            49:   cmo *y;
        !            50:   double x0,x1,h;
        !            51:   double abs_err = 1e-6;
        !            52:   double rel_err = 1e-5;
        !            53:   cmo *err;
        !            54:   char *method;
        !            55:   cmo **cr;
        !            56:   int i;
        !            57:   cmo *ans;
        !            58:   const gsl_odeiv_step_type *T;
        !            59:   method = type;
        !            60:   //  gsl_set_error_handler_off();
        !            61:   gsl_set_error_handler((gsl_error_handler_t *)myhandler);
        !            62:   argc = get_i(); // number of args
        !            63:   if (argc < 5) {
        !            64:     pops(argc);
        !            65:     push(make_error2("The argc must be >=5 for gsl_odeiv_step.",NULL,0,-1));
        !            66:     return;
        !            67:   }
        !            68:   f = pop(); /* ff is a list of cmo's y'=f(x,y); f is a vector valued function */
        !            69:   dim = get_length(f);
        !            70:   GSL_ode_dim=dim;
        !            71:   if (dim < 1) {
        !            72:     push(make_error2("gsl_odeiv_step: the first argument must be a list",NULL,0,-1));
        !            73:     return;
        !            74:   }
        !            75:   for (i=0; i<dim; i++) {
        !            76:     GSL_ode_func[i] = element_of_at(f,i);
        !            77:   }
        !            78:   y = pop(); /* initial value of y (vector) */
        !            79:   if (get_length(y) != dim) {
        !            80:     push(make_error2("gsl_odeiv_step: dim f != dim y (initial value).",NULL,0,-1));
        !            81:     return;
        !            82:   }
        !            83:   x0 = get_double(); // start
        !            84:   x1 = get_double(); // end
        !            85:   h  = get_double(); // step
        !            86:   if (argc >= 6) {
        !            87:     err = pop();
        !            88:     // not implemented. set
        !            89:   }
        !            90:   if (argc >= 7) {
        !            91:     method = get_string();
        !            92:   }
        !            93:
        !            94:   if (strcmp(method,"rk4")==0) {
        !            95:     T = gsl_odeiv_step_rk4;
        !            96:   }else{ // default
        !            97:     T = gsl_odeiv_step_rk4;
        !            98:   }
        !            99:   gsl_odeiv_step *s = gsl_odeiv_step_alloc(T, 2);
        !           100:
        !           101:   gsl_odeiv_control *c = gsl_odeiv_control_y_new(abs_err, rel_err);
        !           102:
        !           103:   gsl_odeiv_evolve *e = gsl_odeiv_evolve_alloc(dim);
        !           104:   gsl_odeiv_system sys = {f_ode, NULL, dim, NULL};
        !           105:   /* x : start, x1 : goal */
        !           106:   double *yy;
        !           107:   int tmp_len;
        !           108:   yy = cmo2double_list(&tmp_len,y);
        !           109:
        !           110:   while (x0 < x1) {
        !           111:     int status = gsl_odeiv_evolve_apply(e, c, s, &sys, &x0, x1, &h, yy);
        !           112:     if (status != GSL_SUCCESS) break;
        !           113:     ox_printf("x=%.5e, yy[0]=%.5e, h=%lg,\n", x0, yy[0], h);
        !           114:     // Todo, store intermediate values.
        !           115:   }
        !           116:
        !           117:   gsl_odeiv_evolve_free(e);
        !           118:   gsl_odeiv_control_free(c);
        !           119:   gsl_odeiv_step_free(s);
        !           120:
        !           121:   cr = (cmo **) GC_malloc(sizeof(cmo *)*(dim+1));
        !           122:   cr[0] = (cmo *)new_cmo_double(x0);
        !           123:   for (i=1; i<=dim; i++) {
        !           124:     cr[i] = (cmo *)new_cmo_double(yy[i-1]);
        !           125:   }
        !           126:   ans = (cmo *)new_cmo_list_array((void *)cr,dim+1);
        !           127:   push(ans);
        !           128: }