使用R deSolve::ode实现带温度触发加热的环形热扩散问题

问题根源 #

你遇到的报错是因为在微分方程右端函数中引入了硬阈值导致的导数不连续点：if (v[50] < 25)的判断会让第50个节点的温度导数在阈值处发生跳变。deSolve默认的隐式求解器（如lsoda）的自适应步长算法假设右端函数连续可微，碰到不连续点时会不断缩小步长尝试拟合，直到超过maxsteps上限报错，或步长过小导致运行速度骤降。
这类带条件触发的微分方程完全可以用ODE求解，只要调整实现方式避开硬不连续点即可，以下是两种可行方案：

方案1：使用deSolve内置事件机制（最准确） #

deSolve专门提供了事件处理框架，用于处理阈值触发这类不连续逻辑，求解器会自动识别事件点、重启步长计算，不会出现卡步长的问题，代码修改如下：

library(deSolve)
library(dplyr)
library(ggplot2)
library(tidyr)

local({
  heatT <- 100
  v <- c(rep(1, 49), heatT, rep(1, 50))
  alpha <- .02
  
  # 根函数：返回值过0时触发事件，这里监测v[50]是否降到25
  rootfun <- function(t, v, pars) {
    return(v[50] - 25)
  }
  
  # 事件函数：触发事件时执行的操作，可根据实际加热需求调整逻辑
  eventfun <- function(t, v, pars) {
    # 对应持续加热逻辑，若只需单次加热到100℃可直接写v[50] <- 100
    v[50] <- v[50] + (100 - v[50]) * 0.1 
    return(v)
  }
  
  fun <- function(t, v, pars) {
    L <- length(v)
    d2T <- c(v[2:L], v[1]) + c(v[L], v[1:(L - 1)]) - 2 * v
    dt <- pars * d2T
    return(list(dt - .005 * (v - 1)))
  }
  
  # 调用ode时传入根函数和事件函数
  ode(v, 1:200, fun, parms = alpha, 
      rootfun = rootfun, 
      events = list(func = eventfun, root = TRUE))
}) %>% as.data.frame() %>% 
  pivot_longer(-time, values_to = "val", names_to = "x") %>% 
  filter(time %in% round(seq.int(1, 200, length.out = 40))) %>%
  ggplot(aes(as.numeric(x), val)) +
  geom_line(alpha = .5, show.legend = FALSE) +
  geom_point(aes(color = val)) +
  scale_color_gradient(low = "#56B1F7", high = "red") +
  facet_wrap(~ time) +
  theme_minimal() +
  scale_y_continuous(limits = c(0, 100)) +
  labs(x = 'x', y = 'T', color = 'T')

方案2：硬阈值光滑近似（实现最简单） #

如果你不需要完全精准的硬阈值，也可以用sigmoid类的光滑函数替换if判断，让导数在阈值附近平滑过渡，求解器就可以正常运行，仅需替换原if行即可：

# 把原if注释行替换为下面的代码
# scale参数越小，过渡区间越窄，越接近硬阈值效果
dt[50] <- dt[50] + (100 - v[50]) * plogis(25 - v[50], scale = 0.1)

这个方案不需要修改其他逻辑，运行速度和原无阈值版本几乎一致，适合对阈值精度要求不高的场景。

两种方案都可以在默认隐式lsoda求解器下快速运行，不需要切换求解方法或者调大maxsteps参数。

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