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

Question

我正在尝试模拟一个环，如果温度低于某个值，该环会在某一点被加热。这是我的 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] 告诉模型如果低于 25°，则增加段 50 的温度。如果注释此行，则模型工作正常。如果该线处于活动状态，则一旦达到 25°，模型就会失败（要求增加 maxsteps）（直到该点仍会输出结果）。如果将求解方法切换为“ode45”，模型可以运行成功，但随后会非常慢，或者如果切换为“euler”等显式方法，但它只能在 alpha 足够低之前工作。

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

Answer 1

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

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

它应该可以正常工作。 vode 是 Livermore ODEPACK family. Another approach is to use an external forcing or an event 的另一个优秀求解器。

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

Heat diffusion a ring with a temperature activated flame using deSolve::ode in R

r

simulation

ode

differential-equations

desolve