Implementation of piecewised function in python with numerical ODE solution. TypeError: float() argument must be a string or a number, not 'OdeResult'

Question

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

Answer 1

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.