Revision history [back]

I guess you mean the example at the end of the Solving ODE numerically by GSL section.

You can add a constructor to your class which takes params as an argument and stores it in the instance created. The code to use this would be (mostly copied from the reference manual)

sage: params = [<value>, <value>, <value>]
sage: T = ode_solver()
sage: T.algorithm = "bsimp"
sage: osc = duffing_osc(params) # note the argument
sage: T.function = osc
sage: T.ode_solve(...)

The new class would look like:

cdef class duffing_osc(sage.gsl.ode.ode_system):
    def __init__(self, params):
        self.params = params

    cdef int c_f(self,double t, double *y,double *dydt):
        dydt[0]=y[1]
        dydt[1]=y[0]-self.params[0]*y[1]-y[0]*y[0]*y[0]+\
                self.params[1]*cos(self.params[2]*t)
        return GSL_SUCCESS
    cdef int c_j(self, double t,double *y,double *dfdy,double *dfdt):
        dfdy[0]=0
        dfdy[1]=1.0
        dfdy[2]=-3*y[0]*y[0]+1
        dfdy[3]=-self.params[0]*y[1]
        dfdt[0]=0
        dfdt[1]=-self.params[1]*self.params[2]*sin(self.params[2]*t)
        return GSL_SUCCESS

In terms of performance, there is a problem with this code.

The elements of params are Sage types, either Python floats, RDF elements, or even multiprecision floats (backed by MPFR). Whenever these are multiplied by a double, Sage will automatically convert this value to a Sage/Python object, do the arithmetic, and convert the result back to a double. This is very slow, compared to a simple float multiplication.

To work around this, you can convert params to a vector of doubles in the constructor. Assuming params has 3 elements, then the first lines of the class would be

cdef class duffing_osc(sage.gsl.ode.ode_system):
    cdef double params[3]
    def __init__(self, params):
        for i in range(3):
            self.params[i] = params[i] # type conversion is automatic