0

After a previous question I am trying to solve some ODEs using GSL. The complicating factor is I want it written in python and so am trying to get scipy.weave to help.

Using the standard GSL ODE example I have the following which is based on this example

import scipy.weave as weave

time =0.0
y1 = [0.0,0.0]
# Define the function 

support_code=r'''

 int
 func (double t, const double y[], double f[],
       void *params)
 {
   double mu = *(double *)params;
   f[0] = y[1];
   f[1] = -y[0] - mu*y[1]*(y[0]*y[0] - 1);
   return GSL_SUCCESS;
 }

 int
 jac (double t, const double y[], double *dfdy, 
      double dfdt[], void *params)
 {
   double mu = *(double *)params;
   gsl_matrix_view dfdy_mat 
     = gsl_matrix_view_array (dfdy, 2, 2);
   gsl_matrix * m = &dfdy_mat.matrix; 
   gsl_matrix_set (m, 0, 0, 0.0);
   gsl_matrix_set (m, 0, 1, 1.0);
   gsl_matrix_set (m, 1, 0, -2.0*mu*y[0]*y[1] - 1.0);
   gsl_matrix_set (m, 1, 1, -mu*(y[0]*y[0] - 1.0));
   dfdt[0] = 0.0;
   dfdt[1] = 0.0;
   return GSL_SUCCESS;
 }
'''


# Define the main

c_main = r'''
 int
 main (void)
 {
   double mu = 10;
   gsl_odeiv2_system sys = {func, jac, 2, &mu};

   gsl_odeiv2_driver * d = 
     gsl_odeiv2_driver_alloc_y_new (&sys, gsl_odeiv2_step_rk8pd,
                  1e-6, 1e-6, 0.0);
   int i;
   double t1 = 100.0;


   for (i = 1; i <= 100; i++)
     {
       double ti = i * t1 / 100.0;
       int status = gsl_odeiv2_driver_apply (d, &t, ti, y);

       if (status != GSL_SUCCESS)
    {
      printf ("error, return value=%d\n", status);
      break;
    }

       time = t;
     }

   gsl_odeiv2_driver_free (d);
   return 0;
 }
'''


compiler = 'gcc'
vars = ['time','y']
libs = ['gsl','gslcblas','m']
headers = ['<stdio.h>','<gsl/gsl_errno.h>','<gsl/gsl_odeiv2.h>','<gsl/gsl_matrix.h>']


res = weave.inline(c_main, vars, compiler=compiler,
                   libraries=libs, headers=headers,
                   support_code = support_code)

print time

which returns this horrible spew. I really don't have a clue what I am doing but know this should be fairly easy! The trouble is all of the examples of weave on the web are based on printf or short iterations.

Community
  • 1
  • 1
Greg
  • 11,654
  • 3
  • 44
  • 50
  • 1
    Your other question mentions pygsl (http://pygsl.sourceforge.net/). I would recommend using that before trying to create your own wrapper. – Warren Weckesser May 03 '13 at 13:09
  • Yes but could not manage to get it installed properly, perhaps I will have another go at that. – Greg May 03 '13 at 13:17
  • 1
    It has worked well in the past, but I haven't tried it in a while. The date of the pygsl tar file on sourceforge is 2010, and GSL 1.15 was released in May, 2011, so pygsl might not be up to date. – Warren Weckesser May 03 '13 at 13:27
  • OK, I just tried PyGSL, and it worked fine. I added a new answer to your other question. Further comments about using PyGSL should be made over there; sorry for going a bit off-topic in this question. – Warren Weckesser May 03 '13 at 19:54

0 Answers0