提供一些数据时计算常微分方程中的参数

Question

代码：

from scipy.integrate import odeint
import numpy as np
import matplotlib.pyplot as plt

# parameters

S = 0.0001
M = 30.03
K = 113.6561
Vr = 58
R = 8.3145
T = 298.15
Q = 0.000133
Vp = 0.000022
Mr = 36
Pvap = 1400
wf = 0.001
tr = 1200
mass = 40000

# define t
time = 14400
t = np.arange(0, time + 1, 1)

# define initial state
Cv0 = (mass / Vp) * wf  # Cv(0)
Cr0 = (mass / Vp) * (1 - wf)
Cair0 = 0  # Cair(0)


# define function and solve ode
def model(x, t):
    C = x[0]  # C is Cair(t)
    c = x[1]  # c is Cv(t)
    a = Q + (K * S / Vr)
    b = (K * S * M) / (Vr * R * T)
    s = (K * S * M) / (Vp * R * T)
    w = (1 - wf) * 1000
    Peq = (c * Pvap) / (c + w * c * M / Mr)
    Pair = (C * R * T) / M
    dcdt = -s * (Peq - Pair)
    if t <= tr:
        dCdt = -a * C + b * Peq
    else:
        dCdt = -a * C
    return [dCdt, dcdt]

x = odeint(model, [Cair0, Cv0], t)

C = x[:, 0]
c = x[:, 1]

现在，当我知道 C(0)(当 t 为 0)和 C(tr)(当 t 为 tr)(因此我知道两种 t 和 C(t ))。

我找到了一些与此相关的链接（Curve Fit Parameters in Multiple ODE Function、Solving ODE with Python reversely、https://medium.com/analytics-vidhya/coronavirus-in-italy-ode-model-an-parameter-optimization-forecast-with-python-c1769cf7a511、https://kitchingroup.cheme.cmu.edu/blog/2013/02/18/Fitting-a-numerical-ODE-solution-to-data/），尽管我无法掌握主题。

我可以用两个data((0, C(0)), (tr, C(tr))和ode来细化参数wf吗？

Answer 1

首先，ODE 求解器采用平滑的右手边函数。因此，您代码中的 if t <= tr:... 语句将不起作用。必须进行两个单独的集成来处理不连续性。积分到 tf，然后使用 tf 处的解作为初始条件，在 tf 之外积分以用于新的 ODE 函数。

但您的主要问题（解决wf）似乎只涉及集成到tf（而不是超出），因此我们在解决wf 时可以忽略该问题

现在，当我知道 C(0)(当 t 为 0)和 C(tr)(当 t 为 tr)时，我想计算 wf 值(因此我知道两种 t 和 C(t)) .

您可以对wf 进行非线性求解：

from scipy.integrate import odeint
import numpy as np
import matplotlib.pyplot as plt

# parameters
S = 0.0001
M = 30.03
K = 113.6561
Vr = 58
R = 8.3145
T = 298.15
Q = 0.000133
Vp = 0.000022
Mr = 36
Pvap = 1400
mass = 40000

# initial condition for wf
wf_initial = 0.02

# define t
tr = 1200
t_eval = np.array([0, tr], np.float)

# define initial state. This is C(t = 0)
Cv0 = (mass / Vp) * wf_initial  # Cv(0)
Cair0 = 0  # Cair(0)
init_cond = np.array([Cair0, Cv0],np.float)

# Definte the final state. This is C(t = tr)
final_state = 3.94926615e-03

# define function and solve ode
def model(x, t, wf):
    C = x[0]  # C is Cair(t)
    c = x[1]  # c is Cv(t)
    a = Q + (K * S / Vr)
    b = (K * S * M) / (Vr * R * T)
    s = (K * S * M) / (Vp * R * T)
    w = (1 - wf) * 1000
    Peq = (c * Pvap) / (c + w * c * M / Mr)
    Pair = (C * R * T) / M
    dcdt = -s * (Peq - Pair)
    dCdt = -a * C + b * Peq
    return [dCdt, dcdt]

# define non-linear system to solve
def function(x):
    wf = x[0]
    x = odeint(model, init_cond, t_eval, args = (wf,), rtol = 1e-10, atol = 1e-10)    
    return x[-1,0] - final_state

from scipy.optimize import root
sol = root(function, np.array([wf_initial]),  method='lm')

print(sol.success)

wf_solution = sol.x[0]
x = odeint(model, init_cond, t_eval, args = (wf_solution,), rtol = 1e-10, atol = 1e-10)

print(wf_solution)
print(x[-1])
print(final_state)