多斑块ODE模型多强迫项适配：R deSolve与编译C代码问题

我正在用R的deSolve包构建多斑块SEIR模型，每个时间步长内不同斑块间存在个体流动，会引发跨斑块感染；同时每个斑块配有独立的外部强迫参数（代表不同环境条件）。模型已在基础R中实现，但因斑块、 compartments数量多且模拟时长较长，希望转换为编译代码提升运行速度。

目前多斑块功能已实现，但为每个斑块添加独立强迫参数时遇到障碍：传入强迫项时，自动执行的checkforcings函数不允许超过两列的矩阵。不确定最优方案是重写ode和checkforcings函数覆盖原有逻辑，还是在C代码中重构强迫项数据？最终模型包含195个斑块，希望实现自动化，避免编写数千个方程或数百个函数。

也可直接建议更换语言，想了解语言选择建议，比如Julia是否合适？

极简示例代码 #

R代码 #

# Packages #########################################################

library(deSolve)
library(ggplot2); theme_set(theme_bw())
library(tidyr)
library(dplyr)


# Initial Parameters and things ####################################

times <- 1:500
n_patch <- 2
patch_ind <- 100
state_names <- (c("S", "I"))
n_state <- length(state_names)

x <- rep(0, n_patch*n_state)
names(x) <- unlist(lapply(state_names, function(x) paste(x, 
                                                         stringr::str_pad(seq(n_patch), width = 3, side = "left", pad =0), 
                                                         sep = "_")))

#start with infected individuals in patch 1
x[startsWith(names(x), "S")] <- patch_ind
x['S_001'] <- x['S_001'] - 5
x['I_001'] <- x['I_001'] + 5
x['I_002'] <- x['I_002'] + 20

params <- c(gamma = 0.1, betam = 0.2)

#seasonality
forcing <- data.frame(times = times,
                      rain = rep(rep(c(0.95,1.05), each = 50), 5))
new_approx_fun <- function(rain.column, t){
  approx_col <- approxfun(rain.column, rule = 2)
      return(approx_col(t))
    }
    
rainfall2 <- data.frame(P1 = forcing$rain, 
                        P2 = forcing$rain+0.01)

# model in R

r.mod2 <- function(t,x,params){
  # turn state.vec into matrix
  # columns are different states, rows are different patches
  states <- matrix(x,
                   nrow = n_patch,
                   ncol = n_state, byrow = F)
  
  S <- states[,1]
  I <- states[,2]
  N <- rowSums(states[,1:2])
  
  with(as.list(params),{
    #seasonal forcing
    rain <- as.numeric(apply(as.matrix(rainfall2), MARGIN = 2, FUN = new_approx_fun, t = t))
    
    dS <- gamma*I - rain*betam*S*I/N 
    dI <- rain*betam*S*I/N  - gamma*I
    
    return(list(c(dS, dI), rain))
  })
}

out.R2 <- data.frame(ode(y = x, times =times, func = r.mod2,
                        parms = params))

#create seasonality for C
ftime <- seq(0, max(times), by = 0.1)
rain.ft <- approx(times, rainfall2$P1, xout = ftime, rule = 2)$y
forcings2 <- cbind(ftime, rain.ft, rain.ft +0.01)

# C model
system("R CMD SHLIB ex-patch-season-multi.c")
dyn.load(paste("ex-patch-season-multi", .Platform$dynlib.ext, sep = ""))
out.dll <- data.frame(ode(y = x, times = times, func = "derivsc",
                          dllname = "ex-patch-season-multi", initfunc = "parmsc",
                          parms = params, forcings = forcings2, 
                          initforc = "forcc", nout = 1, outnames = "rain"))

C代码 #

#include <R.h>
#include <math.h>
#include <Rmath.h>
// this is for testing to try and get different forcing for each patch //

/*define parameters, pay attention to order */

static double parms[2];
static double forc[1];

#define gamma parms[0]
#define betam parms[1]
//define forcing
#define rain forc[0]



/* initialize parameters */
void parmsc(void (* odeparms)(int *, double *)){
  int N=2;
  odeparms(&N, parms);
}

/* forcing */
void forcc(void (* odeforcs)(int *, double *))
{
  int N=1;
  odeforcs(&N, forc);
}

/* model function */
void derivsc(int *neq, double *t, double *y, double *ydot, double *yout, int *ip){
  
  //use for-loops for patches
  //define all variables at start of block
  int npatch=2;
  double S[npatch]; double I[npatch]; double N[npatch];
  int i; 
  for(i=0; i<npatch; i++){
    S[i] = y[i];
  };
  for(i=0; i <npatch; i++){
    int ind = npatch+i;
    I[i] = y[ind];
  };
  for(i=0; i<npatch; i++){
    N[i] = S[i] + I[i];
  };
  
  //use for loops for equations 
  {
    // Susceptible
    for(i=0; i<npatch; i++){
      ydot[i] = gamma*I[i] - rain*betam*I[i]*S[i]/N[i] ;
    };
    //infected
    for(i=0; i<npatch; i++){
      int ind=npatch+i;
      ydot[ind] = rain*betam*I[i]*S[i]/N[i] - gamma*I[i];
    };
  };
  
  yout[0] = rain;
  
}

一、在deSolve框架内解决多强迫项问题 #

1. 手动传入强迫数据，绕过checkforcings限制 #

无需重写核心函数，可将多斑块强迫数据打包为参数的一部分传入，在C代码中自行处理插值：

• R端：将所有斑块的强迫时间序列合并为大向量，同时传入时间点、斑块数等元信息作为参数。
• C端：根据当前时间t，对每个斑块的强迫序列做线性插值，获取对应值。

2. 使用函数式强迫项替代矩阵 #

deSolve支持传入返回多值的强迫函数，替代矩阵格式：

• R端定义强迫函数，输入时间t返回所有斑块的强迫值向量：

  forcing_fun <- function(t) {
    sapply(rainfall2, function(col) approxfun(col, rule=2)(t))
  }
  

• 调用ode时指定forcings = forcing_fun，同时修改C代码中forc数组的大小为n_patch，调整forcc函数接收所有强迫值。

二、语言选择：Julia是理想替代方案 #

对于195个斑块的大规模模型，Julia具备显著优势：

• 原生高性能：代码速度接近C，无需额外编译步骤，直接编写模型即可获得高效计算能力。
• ODE生态成熟：DifferentialEquations.jl包支持大规模系统、自定义强迫项、并行计算，处理多斑块逻辑极为便捷。
• 代码简洁自动化：用数组和循环即可轻松处理所有斑块，无需手动编写数千个方程，代码结构与R类似但效率更高。
• 强迫项处理灵活：可直接在模型函数中为每个斑块调用独立的插值函数，无列数限制。

简化版Julia示例框架：

using DifferentialEquations, Plots

n_patch = 2
timespan = (1.0, 500.0)
# 为每个斑块创建插值函数
rain_interps = [approx(1:500, rainfall2[:,i], rule=2) for i in 1:n_patch]

# 模型函数
function seir!(du, u, p, t)
    gamma, betam = p
    S = @view u[1:n_patch]
    I = @view u[n_patch+1:2n_patch]
    N = S .+ I
    # 获取每个斑块的强迫值
    rain = [interp(t) for interp in rain_interps]
    
    du[1:n_patch] = gamma .* I .- rain .* betam .* S .* I ./ N
    du[n_patch+1:2n_patch] = rain .* betam .* S .* I ./ N .- gamma .* I
end

# 初始条件与参数
u0 = [fill(100, n_patch); fill(0, n_patch)]
u0[1] -=5; u0[n_patch+1] +=5; u0[n_patch+2] +=20
p = (0.1, 0.2)

# 求解模型
prob = ODEProblem(seir!, u0, timespan, p)
sol = solve(prob, Tsit5())

内容的提问来源于stack exchange，提问作者m.evans