使用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((直到该点，它仍会输出结果(。如果求解方法切换到"0"，则模型可以成功运行；ode45"；，但是随后非常慢，或者如果切换到显式方法；euler"；但它只能工作到alpha足够低。

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

Answer 1

似乎if线使模型变得非常僵硬。这并不奇怪，因为常微分方程在定义上是连续的和可微的。在实际情况中违反这一点并不罕见，但幸运的是，求解器非常健壮。然而，总是可以"；将解算器靠墙驱动"；，这里的情况似乎就是这样。在这种情况下有几种可能性：调整公差，通过使用具有圆角边缘的矩形信号使信号更平滑，改变网格。有时，一个更健壮的求解器就可以了。默认的lsoda对大多数应用程序来说都很好，但在这种情况下vode会更好。将对ode的调用替换为以下行：

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

而且它应该可以正常工作。vode是Livermore ODEPACK家族的另一个优秀的求解器。另一种方法是使用外部强制或事件。