deSolve：具有两个连续动力学的微分方程

Question

我正在使用 deSolve::ode() 模拟具有流动水和温度梯度的环形管。环被建模为一个向量，其中每个元素都有一个温度值和位置。

我正在对热扩散公式建模：

1)

但我也在努力沿着圆环移动水。理论上，就是将tube vector中i元素的温度替换为更早的s元素的温度。由于s不一定是整数，所以可以分为整数部分（n）和小数部分（p ): s=n+p。因此，由于水的移动导致的温度变化变为：

2)

问题是 s 等于 v 由 dt 评估的水流速度 v在 ode 求解器的每次迭代中。

我的想法是把现象看做加法，即先算(1)，再算(2)，最后加起来。我担心时间的影响。带有隐式方法的 ode 求解器自动确定时间步长并线性缩小单一变化增量。

我的问题是在导数函数中只返回 (1) + (2) 是否正确，或者我是否应该将这两个过程分开并分别计算导数。在第二种情况下，建议的方法是什么？

编辑：根据@tpetzoldt 的建议，我尝试使用 ReacTran::advection.1D() 实现水流。我的模型有多个温度变化源：自发对称热扩散；水流；一个热源，如果传感器附近（放置在热源之前）的温度低于下限阈值，则开启该热源，如果高于上限阈值则关闭该热源；由周期性外部温度决定的恒定散热。

在“流水”部分下面还有我之前版本的代码，现在被ReacTran::advection.1D()代替了。 plot_type 参数允许可视化水管（“管道”）中温度的时间序列，或传感器（加热器前后）的温度序列。

library(deSolve)
library(dplyr)
library(ggplot2)
library(tidyr)
library(ReacTran)

test <- function(simTime = 5000, vel = 1, L = 500, thresh = c(16, 25), heatT = 25,
                                    heatDisp = .0025,   baseTemp = 15, alpha = .025,
                                    adv_method = 'up', plot_type = c('pipe', 'sensors')) {
    
    plot_type <- match.arg(plot_type)

    thresh <- c(16, 25)

    sensorP <- round(L/2)

    vec <- c(rep(baseTemp, L), 0)

    eventfun <- function(t, y, pars) {

        heat <- y[L + 1] > 0

        if (y[sensorP] < thresh[1] & heat == FALSE) { # if heat is FALSE -> T was above the threshold
            #browser()
            y[L + 1] <- heatT
        }

        if (y[sensorP] > thresh[2] & heat == TRUE) { # if heat is TRUE -> T was below the threshold
            #browser()
            y[L + 1] <- 0
        }

        return(y)
    }

    rootfun <- function (t, y, pars) {

        heat <- y[L + 1] > 0

        trigger_root <- 1

        if (y[sensorP] < thresh[1] & heat == FALSE & t > 1) { # if heat is FALSE -> T was above the threshold
            #browser()
            trigger_root <- 0
        }

        if (y[sensorP] > thresh[2] & heat == TRUE & t > 1) { # if heat is TRUE -> T was below the threshold
            #browser()
            trigger_root <- 0
        }


        return(trigger_root)
    }

    roll <- function(x, n) {
        x[((1:length(x)) - (n + 1)) %% length(x) + 1]
    }

    fun <- function(t, y, pars) {

        v <- y[1:L]

        # Heat diffusion: dT/dt = alpha * d2T/d2X
        d2Td2X <- c(v[2:L], v[1]) + c(v[L], v[1:(L - 1)]) - 2 * v

        dT_diff <- pars * d2Td2X

        # Moving water
        # nS <- floor(vel)
        # pS <- vel - nS
        #
        # v_shifted <- roll(v, nS)
        # nS1 <- nS + 1
        # v_shifted1 <- roll(v, nS + 1)
        #
        # dT_flow <- v_shifted + pS * (v_shifted1 - v_shifted) - v
        dT_flow <- advection.1D(v, v = vel, dx = 1, C.up = v[L], C.down = v[1],
                                                        adv.method = adv_method)$dC

        dT <- dT_flow + dT_diff

        # heating of the ring after the sensor
        dT[sensorP + 1] <- dT[sensorP  + 1] + y[L + 1]

        # heat dispersion
        dT <- dT - heatDisp * (v - baseTemp + 2.5 * sin(t/(60*24) * pi * 2))

        return(list(c(dT, 0)))
    }

    out <- ode.1D(y = vec, times = 1:simTime, func = fun, parms = alpha, nspec = 1,
                                    events = list(func = eventfun, root = T),
                                    rootfunc = rootfun)


    if (plot_type == 'sensors') {

        ## Trend of the temperature at the sensors levels
        out %>%
            {.[,c(1, sensorP + 1, sensorP + 3, L + 2)]} %>%
            as.data.frame() %>%
            setNames(c('time', 'pre', 'post', 'heat')) %>%
            mutate(Amb = baseTemp + 2.5 * sin(time/(60*24) * pi * 2)) %>%
            pivot_longer(-time, values_to = "val", names_to = "trend") %>%
            ggplot(aes(time, val)) +
            geom_hline(yintercept = thresh) +
            geom_line(aes(color = trend)) +
            theme_minimal() +
            theme(panel.spacing=unit(0, "lines")) +
            labs(x = 'time', y = 'T°', color = 'sensor')
    } else {

    ## Trend of the temperature in the whole pipe
    out %>%
        as.data.frame() %>%
        pivot_longer(-time, values_to = "val", names_to = "x") %>%
        filter(time %in% round(seq.int(1, simTime, length.out = 40))) %>%
        ggplot(aes(as.numeric(x), val)) +
        geom_hline(yintercept = thresh) +
        geom_line(alpha = .5, show.legend = FALSE) +
        geom_point(aes(color = val)) +
        scale_color_gradient(low = "#56B1F7", high = "red") +
        facet_wrap(~ time) +
        theme_minimal() +
        theme(panel.spacing=unit(0, "lines")) +
        labs(x = 'x', y = 'T°', color = 'T°')
    }
}

有趣的是，设置更多的段数 (L = 500) 和高速 (vel = 2) 可以在 post 加热传感器中观察到尖峰序列。此外，处理时间急剧增加，但更多的是速度增加的影响，而不是管道分辨率增加的影响。

我现在最大的疑问是 ReacTran::advection.1D() 在我的上下文中是否有意义，因为我正在模拟水温，而这个函数似乎与流水中溶质的浓度更相关。

Answer 1

该问题看起来像是一个具有移动相和固定相的 PDE 示例。可以在 Soetaert 和 Meysman (2012) doi.org/10.1016/j.envsoft.2011.08.011.

关于 ReachTran 的论文中找到关于使用 R/deSolve 的“线法”(MOL) 方法的一个很好的介绍

可以在 slide 55 of some workshop slides, more in the teaching package RTM.

找到 PDE 示例

R/deSolve/ReacTran 试图让 ODEs/PDEs 变得容易，但陷阱仍然存在。如果出现数值离散或振荡，则可能是由于违反了 Courant–Friedrichs–Lewy condition.

deSolve：具有两个连续动力学的微分方程

deSolve: differential equations with two consecutive dynamics

r

ode

desolve