Haskell 中的自动雅可比矩阵

Question

我正在尝试求解 ODE 集，但遇到了一些问题，因为我正在研究的集非常僵硬。我的相关经验很低，但我知道我的下一步应该是在我正在使用的求解器中包含雅可比矩阵。

我当前的设置使用 Numeric.GSL.Ode 中的 OdeSolveV 和

ODEMethod RKf45

。我的目标是将

ODEMethod

升级为

BSimp

。 ODE 示例中给出了包含雅可比矩阵的示例，其中

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

据我了解，此函数

jac

返回雅可比矩阵，其中包含

x

和

v

的当前值。到目前为止，这非常简单，除非我错过了一些东西。然而，我正在使用的系统没有 2 个参数，而是 35 个参数，并且我不想手动执行此矩阵，因此我正在尝试找出一种自动执行此操作的方法。

Numeric.AD.Jacobian似乎有一种方法可以准确地给我这个，但我不明白如何实现它。我认为在这种情况下我想要的是一个函数

jac

，它根据

xdot

和传入的参数返回雅可比矩阵。这可行吗？

Answer 1

最简单（尽管不一定是最优雅）的解决方案是将

xdot

定义为多态数值参数的列表函数，因为它们可以直接实例化为用于 ODE 求解器的数字类型，或者自动微分值 -无穷小对。然后，您只需根据求解器的需要将每个版本包装在向量/矩阵中：

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)