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

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

我找到了一些与此相关的链接(, , 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[ 时忽略该问题=19=]

Now, I want to figure out wf value when I know C(0)(when t is 0) and C(tr)(when t is tr)(Therefore I know two kind of t and 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)

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

Figure out parameter in ordinary differential equation when some data provided

python

parameters

ode