如何为Haskell中Van der Pol方程的xdot函数添加参数？

问题背景 #

我想给一段求解Van der Pol方程的Haskell代码添加参数，原代码如下：

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)

我尝试传递参数mu给xdot，但不清楚正确写法，尝试代码如下：

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 mu . toList)
                   0.1 1e-8 1e-8
                   (\t -> fromList . xdot t mu . toList)
                   (fromList [0.5,0]) ts
 where xdot :: Num a => a -> [a] -> [a]
       xdot t mu [x,v] = [v, -x*mu*(1-x^2)]
       jac :: Double -> [Double] -> [[Double]]
       jac t = jacobian (xdot $ realToFrac t mu)
       mu = 2 :: Double

请问如何修改代码才能将Double或(Double, Double)类型的参数传递给xdot？这是否需要修改jac函数？

解决方案 #

核心思路 #

要给xdot传递额外参数，需调整参数顺序以适配自动微分的jacobian函数，同时jac函数必须同步修改——因为它依赖xdot生成雅可比矩阵。

单参数（Double）的正确实现 #

我们可以把参数作为vanderpol的输入，让函数更灵活：

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

vanderpol :: Double          -- ^ Van der Pol参数mu
          -> Vector Double   -- ^ 时间点
          -> [Vector Double]
vanderpol mu ts = toColumns $
         odeSolveV (BSimp $ \t -> fromLists . jac t mu . toList)
                   0.1 1e-8 1e-8
                   (\t -> fromList . xdot t mu . toList)
                   (fromList [0.5,0]) ts
 where 
   -- 调整参数顺序：先固定mu和t，再接收状态变量[x,v]
   xdot :: Num a => a -> a -> [a] -> [a]
   xdot t mu [x,v] = [v, -x * mu * (1 - x^2)]
   
   -- jac函数接收mu，传递给xdot后生成雅可比矩阵
   jac :: Double -> Double -> [Double] -> [[Double]]
   jac t mu = jacobian (xdot (realToFrac t) (realToFrac mu))

多参数（(Double, Double)）的实现 #

如果需要传递多个参数（比如mu和k），只需打包成元组调整参数即可：

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

vanderpol :: (Double, Double)    -- ^ 参数元组(mu, k)
          -> Vector Double       -- ^ 时间点
          -> [Vector Double]
vanderpol (mu, k) ts = toColumns $
         odeSolveV (BSimp $ \t -> fromLists . jac t (mu, k) . toList)
                   0.1 1e-8 1e-8
                   (\t -> fromList . xdot t (mu, k) . toList)
                   (fromList [0.5,0]) ts
 where 
   xdot :: Num a => a -> (a, a) -> [a] -> [a]
   xdot t (mu, k) [x,v] = [v, -k * x * mu * (1 - x^2)]
   
   jac :: Double -> (Double, Double) -> [Double] -> [[Double]]
   jac t (mu, k) = jacobian (xdot (realToFrac t) (realToFrac mu, realToFrac k))

关键说明 #

1. 参数顺序调整：把固定参数（mu、t）放在状态变量前面，确保jacobian仅对状态变量（x、v）求导，固定参数会被当作常量处理。
2. jac函数的修改：必须同步修改jac函数，让它接收额外参数并传递给xdot，否则雅可比矩阵无法正确反映含参数的微分方程。
3. 灵活性提升：将参数作为vanderpol的输入参数，而非在where块硬编码，方便调用时动态调整参数值。

内容的提问来源于stack exchange，提问作者Pompan