Converting GSL ODE solver to Python

Question

I solve a set of couples ODEs which I solve using the GSL ODE solver similar to this example. Currently this is automates by writing a file in python e.g

text = """
#include <stdio.h>
#include <gsl/gsl_errno.h>

...

"""

Then replacing strings in text with the relevant words and writing to a file script.c. I then using os.system to run it e.g

os.system("gcc -Wall -I/usr/local/include -c script.c")
os.system("gcc -static script.o -lgsl -lgslcblas -lm")
os.system("./a.out >  %s" % (file_name) )

All of this very inelegant and so I have been reading about alternatives and have stumbled PyGSL, CythonGSL so far. These seem to lack proper documentation and I'm not really smart enough to work out how they work!

Any advise would be appreciated.

Answer 1

PyGSL might be a reasonable option. (CythonGSL, too, but I haven't tried it.)

I'm running Ubuntu 12.04 (64 bit), and for this installation, I am using the Anaconda python distribution. (It might work fine with the default Python in Ubuntu, but I didn't try.)

I have GSL 1.15 installed from source, and I just installed PyGSL 0.9.5 from source by following the instructions in the README. That is, I ran

$ python setup.py build
$ python setup.py install

in the top directory of the extracted tar file. Like many python packages that build C libraries, a lot of warnings are printed during the build step, but it completed without errors.

Here's some code that demonstrates the use of the ODE solver in PyGSL:

#
# pendulum_demo.py
#
# Use PyGSL to solve the differential equations for a pendulum with
# friction.  Plot the solution using matplotlib.
#

import numpy as np
import matplotlib.pyplot as plt
from pygsl import odeiv


def pendsys(t, y, args):
    """
    The right-hand side of the first order system of differential equations.
    """
    theta, v = y
    g, b, L, m = args

    f = np.zeros((2,))
    f[0] = v
    f[1] = -b * v / (m * L ** 2) - np.sin(theta) * g / L

    return f


def pendjac(t, y, args):
    """
    The Jacobian of the system.
    """
    theta, v = y
    g, b, L, m = args

    jac = np.zeros((2, 2))
    jac[0, 1] = 1.0
    jac[1, 0] = -g / L * np.cos(theta)
    jac[1, 1] = -b / (m * L ** 2)

    dfdt = np.zeros((2,))

    return jac, dfdt


# Pendulum parameter values
#
# Gravitational constant
g = 9.81
# Friction coefficient
b = 0.5
# Pendulum length
L = 1.0
# Pendulum bob mass
m = 2.0

# Initial conditions
theta0 = np.pi - 0.01
v0 = 0.0
y = [theta0, v0]
t = 0

# Solver control parameters.
abserr = 1e-11
relerr = 1e-9

stoptime = 12.0

# Create the GSL ODE solver
N = 2
step    = odeiv.step_rk8pd(N, pendsys, pendjac, args=(g, b, L, m))
control = odeiv.control_y_new(step, abserr, relerr)
evolve  = odeiv.evolve(step, control, N)

# h is the initial step size.
h = stoptime / 500.0

# The time values and points in the solution are saved in the lists
# tsol and ysol, respectively.  The lists are initialized with
# the initial conditions.
tsol = [t]
ysol = [y]

# Call evolve.apply(...) until the solution reaches stoptime
while t < stoptime:
    t, h, y = evolve.apply(t, stoptime, h, y)
    tsol.append(t)
    ysol.append(y)

tsol = np.array(tsol)
ysol = np.array(ysol)

plt.plot(tsol, ysol[:, 0], label=r'$\theta$')
plt.plot(tsol, ysol[:, 1], label='v')
plt.xlabel('t')
plt.grid(True)
plt.legend()
plt.title('Pendulum with Friction')
plt.show()

The resulting plot: