1

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

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

4

1 回答 1

0

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

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

它应该可以正常工作。vode是Livermore ODEPACK系列的另一个出色求解器。另一种方法是使用外部强制或事件。

于 2021-11-30T23:49:56.607 回答