0

I have a reaction system whose mechanism is defined in the form of a nonlinear first-order ODE as:

dx(t)/dt = - gaf(1 - rho)(1 - exp(-kt))(x(t)^2 + bx(t))/(cx(t) + p)

where, g, a, f, b, c, rho, k and p are all constants.

I’m trying to compile a code in python to obtain the time trajectory of the variable x(t) from time t1 to time tf where the kinetic range of the profile is always preceded by a waiting period of 2 min before the reaction starts with constant x0 value in that period starting a time t0 = 0. My aim is to export the whole profile to excel (as .xlsx file) in regular spaced intervals of 2.5 secs as shown in the link to the figure below.

Overall Trajectory Profile - Plotted Curve

I wish the system be defined in the form of a piecewised function, but when I run the code below in python it gives me the error: "float() argument must be a string or a number, not 'OdeResult'". I’m stuck here! Any suggestions on how to remedy this problem or alternatives?

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

f = 1
rho = 0.1
r = 0.5
a = 0.08
x0 = 2
b = x0*(1/r - 1)
k = 1.0e-2
g = 0.138
b = 2
c = 1.11
p = 2.174

tau = (1/k)*np.log(f/(f-rho))    # Inhibition time (induction period)

t0 = 0
ti = 120
t1 = ti + tau
tf = 1000
delta_t = 2.5
t_end = tf + delta_t
dt = 1.0e-5
t_eval=np.arange(t0, t_end, dt)

def ode(t, x):
    return -g*a*f*(1-rho)*(1-np.exp(-k*t))*(x**2 + b*x)/(c*x + p)
sol = solve_ivp(ode, [t0, t_end], [x0], t_eval=t_eval, method='RK45')

plt.figure(figsize = (12, 4))
plt.subplot(121)
curve, = plt.plot(sol.t, sol.y[0])
plt.xlabel("t")
plt.ylabel("x(t)")
plt.show()

t = np.arange(t0, tf, delta_t)

def x(t):
    if(t0 <= t <= ti): return x0    # Waiting period
    elif(ti < t < t1): return x0    # Lag period associated to 'tau' (inhibition period)
    else: return sol                # Dynamic range (problem here!!). TypeError: float() argument must be a string or a number, not 'OdeResult'

y = []
for i in range(len(t)):
   y.append(x(t[i]))

plt.plot(t, y, c='red', ls='', ms=5, marker='.')
ax = plt.gca()
ax.set_ylim([0, x0])
plt.show()
CCMJ
  • 1
  • 1

1 Answers1

0

Replace t_eval=t_eval with dense_output=True. Then the returned sol object contains an interpolation function sol so that you could complete your x function as

    return sol.sol(t)[0]

You do not need an else branch if the branch before is concluded with a return statement.

You might want to include the derivative zero on the constant segment inside the ODE function, so that the x(t) function is not needed. Then you would get the desired values also direct with the instant evaluation using t_eval as done originally.

Lutz Lehmann
  • 25,219
  • 2
  • 22
  • 51
  • I'm not sure how your last suggestion could be implemented, but it is definitely the best choice to the piecewised representation. When is replaced by the resulting solution is numerically less accurate than the original. I've changed the integration method to which gave a much better agreement with the original profile. Yet, there are obvious differences. I would therefore be more interested on how to include the zero derivative inside the ODE function. If you could exemplify I'd be much appreciated. I'm new to python! – CCMJ Mar 18 '22 at 08:53
  • You would just insert `if t – Lutz Lehmann Mar 18 '22 at 09:32