在 R 中使用 deSolve::ode 热扩散具有温度激活火焰的环
Heat diffusion a ring with a temperature activated flame using deSolve::ode in R
我正在尝试模拟一个环,如果温度低于某个值,该环会在某一点被加热。这是我的 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 无法管理的东西?
似乎 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 代码:
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 无法管理的东西?
似乎 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 的另一个优秀求解器。