Figure out parameter in ordinary differential equation when some data provided

Question

Code:

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]

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

I found some links(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/) related to this, although I cannot get the hang of subject.

Can I fine parameter wf with two data((0, C(0)), (tr, C(tr)) and ode?

Answer 1

First, ODE solvers assume smooth right-hand-side functions. So the if t <= tr:... statement in your code isn't going to work. Two separate integrations must be done to deal with the discontinuity. Integrate to tf, then use the solution at tf as initial conditions to integrate beyond tf for the new ODE function.

But it seems like your main problem (solving for wf) only involves integrating to tf (not beyond), so we can ignore that issue when solving for wf

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

You can do a non-linear solve for 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)