求解 ODE - 只有正解答案

【问题标题】：Solving ODEs - only positive solutions求解 ODE - 只有正解
【发布时间】：2017-11-20 21:48:24
【问题描述】：

我正在尝试解决仅限于正解的 ODE，即：

dx/dt=f(x)

x>=0。

在 MATLAB 中，这很容易实现。 R 是否有任何解决方法或包将解决方案空间限制为仅正值？

这对我来说非常重要，不幸的是没有其他选择。我现在搜索了一段时间，但没有任何成功。 :-(

【问题讨论】：

您能说一下（技术上）它是如何在 MATLAB 中完成的吗？最简单的方法是将系统更改为 d(log(x))/dt = f(log(x))。
没有任何示例代码和示例数据，很难更具体。因此，在回答您的问题时，是的，这是可能的。看看R包deSolve，你可以在里面subset满足一定约束的变量。
@MauritsEvers，我不确定这是否容易，但我同意需要更多细节。
很难解释，但这里是用于 MATLAB 的 ode45 radford.edu/~thompson/RP/nonnegative.pdf 的原始出版物。 R 的 deSolve 也使用类似的代码，但它不能强制执行非负性。我只是想用 x>=0 来解决所有 x 的系统 dx/dt=f(x)。假设（仅举一个例子）我有一个生物系统 f(x,p)，其测量浓度为 y，参数为 p。显然，系统的统计数据和观察结果必须是正数或等于零。现在我想推断参数，但 y 和 y 必须保持正数。
我的意思是 x 和 y 必须保持积极态度。我知道这是一个非常普遍的问题，但我真的需要帮助。由于我有一个给定的网络结构，因此很难对系统进行日志转换。我现在搜索了半年的解决方案......但没有任何成功......不幸的是......

标签： r constraints ode

【解决方案1】：

这里还不够。对于我熟悉的各种问题，修改系统以在对数转换后的状态变量的规模上运行效果很好（您始终可以对结果进行反向转换，例如将它们与数据进行比较）。例如，我已将其与SIR model in epidemiology 一起使用。我将尝试使用@MauritsEver 的示例来说明如何将系统转换为在对数刻度上运行：

library(deSolve)
model <- function (time, y, parms) {
   with(as.list(c(y, parms)), {
       dlogN <-   r * (1 - exp(logN) / K)
       list(dlogN)
   })
}

# Starting conditions
y <- c(logN = log(0.1))
parms <- c(r = 0.1, K = 10)
times <- seq(0, 100, 1)
out <- as.data.frame(ode(y, times, model, parms))
out_backtran <- transform(out,N=exp(logN))
plot(N~time,data=out_backtran)

这种方法有以下缺点：

它不会处理恰好在边界上的解，并且会遇到“过快”接近边界的解（即状态变量在有限时间内收敛到零）的问题
如其所写，需要手动翻译。完全有可能编写一个允许用户输入一组方程和一组变换并自动应用变换的系统，但这需要一些努力。
它可能会稍微增加计算量（任何时候我们必须使用状态变量的原始尺度值来求幂）

【讨论】：

只是为了我的理解...您是否错过了 r * N * (1 - N / K) 中的 N？我想避免这个问题，但当然这会迫使变量保持积极。我真的不明白为什么有一个 ocatve，MATLAB 和 MAPLE 能够在 ODE 系统中强制非负性但 R 不是......
不，我没有错过。它消失了，因为 d(log(N))/dt = (dN/dt)/N ...

【解决方案2】：

没有关于 ODE 的任何具体示例代码或详细信息，很难更具体。可能很简单，具体取决于问题。

这是一个使用deSolve 及其函数deSolve::subset 的简单示例。

# Example straight from the deSolve manual
library(deSolve);
model <- function (time, y, parms) {
    with(as.list(c(y, parms)), {
        dN <-   r * N * (1 - N / K);
        list(dN)
    })
}

# Starting conditions
y <- c(N = 0.1);
parms <- c(r = 0.1, K = 10);
times <- seq(0, 100, 1);

# Solve ODE and plot
out <- ode(y, times, model, parms);
plot(out, type = "l", xlim = c(0, 100));

我们现在对解决方案 time 和 subset 施加约束。

# Constrain: time > 20 and plot
out.constrained <- subset(out, select = c("time", "N"), subset = time > 20);
plot(out.constrained, type = "l", xlim = c(0, 100));

【讨论】：

我明白了。但不幸的是，这个想法不仅仅是删除所有与 y 值低于零匹配的 x 值。例如，函数 y=sqrt(x) 可以求解 y>0 和 yradford.edu/~thompson/vodef90web/vodef90source/misc.html 具有此功能。我还与 odeSolve 的作者争论了一段时间。但没有成功......我只想解决限制为x>0的约束ODE系统。我真的不确定如何构建一个有意义的示例，因为这个约束基本上可以应用于任何 ODE 系统。
我没有关注你。您是在谈论耦合 ODE 的约束系统吗？在单个 ODE 的情况下，施加约束 x>=0 只是在约束 x(t) 的图像，不是吗？
问题在于精确解为正但收敛到零。数值求解器，尤其是显式方法，往往会越过零线。正如 cmets 中所建议的，可以使用替换 x=exp(u) 来强制执行积极性。

【解决方案3】：

@Anto 添加一个 if 语句将其设置为零

dA_solid = -k1*(A_sat - A_bulk)
if((A_solid+dA_solid)<0 {dA_solid = -1*(A_solid)}
dA_bulk = -k2*A_bulk*B + k1*(A_sat - A_bulk)

【讨论】：

虽然这会有所帮助，但数字噪声仍有可能将 A_solid 推为负值。

【解决方案4】：

我想我也有同样的问题，可能会给你一个例子。

考虑一个溶解过程，其中产品 A_solid 在溶剂中以速率常数 k_1 溶解到 A_bulk 中（该反应可以在两个方向上进行）。 A_solid 溶解在溶剂中，直到 A_bulk 达到饱和 A_sat。此外，A_bulk 与产物 B 以速率常数 k_2 反应生成 C。这是反应的图片：

Dissolution process

这是我为反应写的代码：

library(deSolve)

# inputs 
T = 0  + 273.15 # K (Kelvin ) / tTemperature 
V = 50 # mL / Volume
A_solid = 125/V # mmol/mL = mol/L / initial concentration of product A_solid
B = 100/V # mol/L / initial concentration of product B

# parameters
R = 8.314 # J/(K*mol) / gas constant
expfact_sat = 2
E_a_sat = 10^3

params <- c(k1 = 0.1, # rate constant of dissolution
            k2 = 3*10^(-3), # rate constant of reaction
            A_sat = expfact_sat*exp(-E_a_sat/(R*T))) # saturation of the A_bulk into the solvent

# initial values
state <- c(A_solid = A_solid, A_bulk = 0, B = B, C = 0)

# system of differential equations
derivs <- function(t, y, parms) {
  with(as.list(c(y, parms)), {
    dA_solid = -k1*(A_sat - A_bulk)
    dA_bulk = -k2*A_bulk*B + k1*(A_sat - A_bulk)
    dB = -k2*A_bulk*B
    dC = k2*A_bulk*B
    return(list(c(dA_solid, dA_bulk, dB, dC)))
  })                                                        
}

times = seq(0, 500, by = 0.01)
init <- ode(y = state, func = derivs, time = times, parms = params)

l = dim(init)[1]-1
matplot(init[,1], init[,-1], type = "l", lty = 1:1, lwd = c(2),
        col = 1:l, xlab = "time [min]", ylab = "concn [mol/L]")
legend("topright", colnames(init)[-1], col = 1:l, lwd = c(2))

这里的问题是，如果你做图你会看到A_solid低于0，这意味着A_solid的浓度是负的，这是无意义的。

最终的情节应该是这样的：

Dissolution

如果您对如何处理此问题有任何建议或解决方案，那就太好了。

@BenBolker：对数转换的问题是，正如你在缺点中所说，A_solid 的浓度在有限时间内变为 0，所以我不确定我们是否仍然可以应用这种技术。

P.S.：我是这个论坛的新手，所以无法直接显示图片。

【讨论】：