如何将F#函数传入RProvider调用R deSolve包的ode求解器

你需要把deSolve包示例里的R版LVmod函数改成F#编写，再传入R的ode函数中，下面是具体实现方案：

原R示例代码 #

LVmod <- function(Time, State, Pars) {
    with(as.list(c(State, Pars)), {
        Ingestion    <- rIng  * Prey * Predator
        GrowthPrey   <- rGrow * Prey * (1 - Prey/K)
        MortPredator <- rMort * Predator
        dPrey        <- GrowthPrey - Ingestion
        dPredator    <- Ingestion * assEff - MortPredator
        return(list(c(dPrey, dPredator)))
    })
}
pars  <- c(rIng   = 0.2,    # /day, rate of ingestion
rGrow  = 1.0,    # /day, growth rate of prey
rMort  = 0.2 ,   # /day, mortality rate of predator
assEff = 0.5,    # -, assimilation efficiency
K      = 10)     # mmol/m3, carrying capacity
yini  <- c(Prey = 1, Predator = 2)
times <- seq(0, 200, by = 1)
out   <- ode(yini, times, LVmod, pars)
summary(out)
plot(out)

F#实现代码 #

先确保项目已引用RProvider和RProvider.deSolve包，再编写如下代码：

open RProvider
open RProvider.deSolve
open RProvider.base

// 定义F#版的Lotka-Volterra模型函数
let lvmod (time: float) (state: float[]) (pars: float[]) =
    // 解析状态变量与参数
    let prey = state.[0]
    let predator = state.[1]
    let rIng = pars.[0]
    let rGrow = pars.[1]
    let rMort = pars.[2]
    let assEff = pars.[3]
    let K = pars.[4]

    // 计算各变量值
    let ingestion = rIng * prey * predator
    let growthPrey = rGrow * prey * (1.0 - prey / K)
    let mortPredator = rMort * predator
    let dPrey = growthPrey - ingestion
    let dPredator = ingestion * assEff - mortPredator

    // 返回符合R要求的列表格式
    R.list([|dPrey; dPredator|])

// 定义参数、初始值和时间序列
let pars = [|0.2; 1.0; 0.2; 0.5; 10.0|] // 对应rIng, rGrow, rMort, assEff, K
let yini = [|1.0; 2.0|] // Prey, Predator
let times = R.seq(0.0, 200.0, by = 1.0).AsNumeric()

// 调用R的ode求解器
let result = R.ode(yini, times, lvmod, pars)

// 查看结果摘要并绘图
R.summary(result) |> ignore
R.plot(result) |> ignore

关键注意点 #

• F#函数参数要和R的ode要求匹配：依次为单浮点数时间、状态数组、参数数组
• 最终必须返回R的list类型，否则deSolve无法识别导数结果
• 参数和初始值用数组传递，对应R中向量的结构
• 调用R函数时用R.xxx的方式，结果可通过AsNumeric()等方法转换为F#类型

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