使用 R 中的 deSolve::ode 对具有温度激活火焰的环进行热扩散

Posted 2023-02-18

【中文标题】使用 R 中的 deSolve::ode 对具有温度激活火焰的环进行热扩散

【英文标题】：Heat diffusion a ring with a temperature activated flame using deSolve::ode in R

【发布时间】：2022-01-07 06:56:39

【问题描述】：

我正在尝试模拟一个在温度低于某个值时会在某一点加热的环。这是我的 R 代码：

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

local(
    heatT <- 100
    v <- c(rep(1, 49), heatT, rep(1, 50))
    alpha <- .02
    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
        
        # Uncomment to trigger the problem
        #if (v[50] < 25) dt[50] <- 100 - v[50]
        
        return(list(dt - .005 * (v - 1)))
    
    
    ode(v, 1:200, fun, parms = alpha)
) %>% 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')

这条线：if (v[50] < 25) dt[50] <- 100 - v[50] 告诉模型如果段 50 的温度低于 25°，则增加其温度。 如果注释了此行，则模型可以正常工作。如果该线处于活动状态，则模型会在达到 25° 时失败（要求增加maxsteps）（它仍然会输出结果直到该点）。 如果求解方法切换到“ode45”，模型可以成功运行，但速度很慢，或者如果切换到像“euler”这样的显式方法，但它只能在 alpha 足够低时运行。

是否有正确的方法来实现它以便使用默认的隐式方法快速运行它，或者它只是 ode 无法管理的东西？

【参考方案1】：

if-line 似乎使模型非常僵硬。这并不奇怪，因为 ODE 在定义上是连续且可微的。在实际案例中违反这一点的情况并不少见，但幸运的是，求解器非常健壮。但是，总是有可能“将求解器推到墙上”，这似乎就是这种情况。在这种情况下有几种可能性：调整容差，通过使用具有圆形边缘的较小矩形信号使信号更平滑一点，更改网格。有时，更强大的求解器会做。默认的lsoda 适用于大多数应用程序，但在这种情况下vode 会更好。将对ode 的调用替换为以下行：

ode(v, 1:200, fun, parms = alpha, method = "vode")

它应该可以正常工作。 vode 是 Livermore ODEPACK 系列的另一个出色求解器。另一种方法是使用外部forcing or an event。

【讨论】：

谢谢，投票成功！！！您能否详细说明或提供其他建议的参考，例如如何调整容差、使信号不那么矩形是什么意思，或者如何更改网格？

再看一遍，vode 产生了奇怪的结果。即使达到 25° 后的新传热是一次性的，峰值温度也会不断升高，直到达到 ~65°，然后再次下降。相反，使用 ode45 它开始在 25° 左右振荡，我相信这是预期的一次热爆发，不是吗？是什么推动了这两种算法之间的差异？当然我更喜欢使用vode，因为它更快。

不管怎样，使用root事件完全解决了问题，谢谢！

没错。这是最直接的方法。我刚刚开始向您指出这种方法。另一种方法是扩展方程，而不是 if，例如使用饱和函数（在生命科学和化学中称为 Michelis-Menten 或 Monod）。关于ode45，我怀疑振荡是不正确的。这听起来有道理，但可能是数值不稳定的迹象。