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

Question

I'm trying to model a ring that is heated at one point if the temperature goes below a certain value. Here's my R code:

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')

The line: if (v[50] < 25) dt[50] <- 100 - v[50] tells the model to increase temperature on segment 50 if it goes below 25°. If this line is commented the model works fine. If the line is active the model fails (asking to increase maxsteps) as soon 25° are reached (it still outputs the results until that point). The model can run successfully if the solving method is switched to "ode45", but then is very slow, or if switched to an explicit method like "euler" but then it works only until alpha is low enough.

Is there a correct way to implement this in order to run it fast with the default implicit methods or it is simply something that ode cannot manage?

Answer 1

It seems that the if-line makes the model very stiff. This is not surprising as ODEs are continuous and differentiable by definition. It is not uncommon that this is violated in practical cases but the solvers are, fortunately, quite robust. However, it is always possible to "drive the solvers against a wall", that seems to be the case here. There are several possibilities ins such cases: tune the tolerances, make the signal a little bit smoother by using a less rectangular signal with rounded edges, change the grid. Sometimes, a more robust solver will do. The default lsoda is fine for most applications, but in this case vode would be better. Replace the call to ode with the following line:

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

and it should work without error. vode is another excellent solver of the Livermore ODEPACK family. Another approach is to use an external forcing or an event.