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

Posted 2023-02-18

【中文标题】deSolve：具有两个连续动力学的微分方程

【英文标题】：deSolve: differential equations with two consecutive dynamics

【发布时间】：2022-01-09 19:03:22

【问题描述】：

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

我正在模拟热扩散公式：

1)

但我也在努力让水沿着环移动。从理论上讲，它只是将管矢量中元素 i 处的温度替换为前面元素 s 处的温度。由于s可能不是整数，所以可以分为整数部分（n）和小数部分（p）：s=n+p。因此，由于水的移动引起的温度变化变为：

2)

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

我的想法是将这些现象视为加法，即首先计算 (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) 可以观察到后加热传感器中的尖峰序列。此外，处理时间急剧增加，但更多的是速度增加的影响，而不是管道分辨率的增加。

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

【问题讨论】：

这不是一个真正的编程问题，但是可以使用包 deSolve 中的函数 ode.1D 来解决此类问题，并且在需要时还可以使用包 ReacTran。网上已经有一些资料（包小插曲、论文），你也可以看看dynamic-r.github.io

Soetaert 和 Meysman (2012) @987654325 关于 ReachTran 的论文中可以找到关于 R/deSolve 的“线法” (MOL) 方法的一个很好的介绍@

所以这意味着简单地将热扩散和水流的 dT 相加是不正确的？

我建议阅读一些材料，然后找到一个好的解决方案。也可以在特定的代码示例中讨论它。

这是一个有趣的案例，但也有些复杂。我做了一些小的技术编辑：库，FALSE/TRUE 而不是 F/T，但我看不到你提供的绘图函数的振荡。无论如何，我建议将问题重新编辑回原始版本（没有代码的版本），然后用代码开始新的版本。

【参考方案1】：

这个问题看起来像一个具有移动和固定相位的 PDE 示例。可以在 Soetaert 和 Meysman (2012) doi.org/10.1016/j.envsoft.2011.08.011 的关于 ReachTran 的论文中找到关于 R/deSolve 的“线方法”(MOL) 方法的一个很好的介绍。

可以在slide 55 of some workshop slides 找到示例 PDE，更多信息请参见教学包RTM。

R/deSolve/ReacTran 试图使 ODE/PDE 变得简单，但仍然存在缺陷。如果出现数值离散或振荡，可能是由于违反了Courant–Friedrichs–Lewy condition。

【讨论】：

最终我使用 ReacTran::advection.1D 来制作水流效果。我不是 100% 确定我在做什么，因为在我的情况下，我没有分散在水中的具有给定浓度的物质，但我正在模拟水温本身。如果您想看一下，我会发布代码。