Automatic Jacobian matrix in Haskell

I'm trying to solve an ODE set, but I've run into some problems, as the set I'm looking into is very stiff. My related experience is low, but I've understood that my next step should be to include a jacobian matrix in the solver I'm using.

My current setup uses OdeSolveV from Numeric.GSL.Ode with ODEMethod RKf45. My goal is to upgrade ODEMethod to BSimp. An example with an included jacobian is given in the ODE examples for hmatrix:

vanderpol' mu = do
    let xdot mu t (toList->[x,v]) = fromList [v, -x + mu * v * (1-x^2)]
        jac t (toList->[x,v]) = (2><2) [ 0          ,          1
                                       , -1-2*x*v*mu, mu*(1-x**2) ]
        ts = linspace 1000 (0,50)
        hi = (ts!1 - ts!0)/100
        sol = toColumns $ odeSolveV (MSBDF jac) hi 1E-8 1E-8 (xdot mu) (fromList [1,0]) ts
    mplot sol

To my understanding this function jac returns the jacobian matrix with the current values of x and v. This is quite straightforward thus far, unless I've missed something. However, the system I'm working with has not 2, but 35 parameters, and I would like to not do this matrix by hand, so I'm trying to figure out a way to do this automatically.

Numeric.AD.Jacobian seems to have a way of giving me exactly this, but I do not understand how to implement it. What I think I want in this scenario is a function jac that returns the jacobian matrix based on xdot and the passed in parameters. Is this viable?

Solution

The easiest (albeit not necessarily most elegant) solution is to define xdot as a list function over polymorphic numerical arguments, as those can be instantiated to either directly the number type used for the ODE solver, or the automatic-differentiation value-infinitesimal pairs. Then you just wrap each version in vectors/matrices as required for the solver:

import Numeric.GSL.ODE
import Numeric.LinearAlgebra
import Numeric.AD

vanderpol :: Vector Double  -- ^ Time points
         -> [Vector Double]
vanderpol ts = toColumns $
         odeSolveV (BSimp $ \t -> fromLists . jac t . toList)
                   0.1 1e-8 1e-8
                   (\t -> fromList . xdot t . toList)
                   (fromList [0.5,0]) ts
 where xdot :: Num a => a -> [a] -> [a]
       xdot t [x,v] = [v, -x*(1-x^2)]
       jac :: Double -> [Double] -> [[Double]]
       jac t = jacobian (xdot $ realToFrac t)

The realToFrac is required because from the perspective of jac, xdot does not accept a Double argument but rather a Reverse s Double argument.