Solution to function of a DE is not plotting correctly for all values of time

Python Help
1

Hello. I am trying to write the following code to have the solution of lam_solution=lamb_i_star for t<=0, and have lam_solution = lam_func(time, dsigma_dt_solution) for t>=0. It keeps giving me a non-existent jump (I am trying to replicate a figure in a paper). Any help would be great. I shouldn’t have to manually draw the line for lamb_i_star when t<0.

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

gamma = 11.2
lamb_i_star = 0.46
lamb_f_star = 0.95
lamb_c = 1.17
phi_rb_i = 0.049
mu_f = lamb_f_star/(1-(lamb_f_star/lamb_c))
phi_rb_i_star = phi_rb_i+(lamb_i_star/gamma)
phi_rb_f_star = phi_rb_i+(lamb_f_star/gamma)
sigma_i = lamb_i_star/phi_rb_i_star
#print(sigma_i)
sigma_f = lamb_f_star/phi_rb_f_star
#print(sigma_f)

def func(t,y):
    sigma = y[0]
    
    dsigma_dt = (mu_f*y[0])*((1-(y[0]/sigma_f))/(1-(y[0]/gamma)))
    
    return [dsigma_dt]

def solve_func(func, t_eval):
    y0 = [lamb_i_star / phi_rb_i_star]
    solution = solve_ivp(func, [t_eval[0], t_eval[-1]], y0, method='LSODA', t_eval=t_eval)
    return solution.y[0]

def lam_func(t, dsigma_dt):
    a = (lamb_f_star * dsigma_dt) / sigma_f
    b = (lamb_f_star * sigma_i) / (lamb_i_star * sigma_f)
    c = (lamb_c / lamb_f_star) - 1
    d = sigma_i / dsigma_dt
    g = (lamb_f_star / lamb_c) - 1
    h = 1 + (c * (sigma_i / gamma))
    i = 1 + (c * (sigma_f / gamma))
    j = 1 + (c * (dsigma_dt / gamma))
    
    return a/(((b-(h/i))*(d**g)*(np.exp(-lamb_f_star*t)))+(j/i))

# Define the time range
time = np.linspace(-2, 5, 100)

# Solve dsigma_dt for all time points
dsigma_dt_solution = solve_func(func, time)

lam_solution = lam_func(time, dsigma_dt_solution)
lam_solution[time < 0] = lamb_i_star

print(lam_solution)

plt.plot(time, lam_solution)
plt.xlabel('t')
plt.ylabel(r'$\lambda$(t)')
plt.yticks(np.arange(0, 1.2, 0.2))
plt.plot([0, 0], [0, 1], 'k', linewidth=2)
plt.plot([0, time[-1]], [lamb_f_star, lamb_f_star], 'r--')
plt.grid(True)
plt.show()
2

I see that you try to check the computed values before plotting (print(lam_solution)). Do these values look right? If something is incorrect, specifically where?

3

The values are correct, yes. It turns out that the issue is with how my np.linspace() is defined, it is circling around t=0 whereas I need to redefine my time span to include t=0 (maybe by using np.arange()). Other than that, there are no issues. Problem solved. Thank you!

4

np.linspace(-2, 5, 100) makes 100 values between -2 and 5, (including both endpoints), that are equally spaced. The spacing between them is therefore 7/99 (rather, the closest approximation of that in floating-point numbers), and none of the resulting values will be equal to zero, yes. Of course, repeatedly adding 7/100 (even in exact fractions rather than floating-point numbers) from a starting point of -2 will also skip over 0 - in the same way that enumerating all the multiples of 7 will skip over 200 (because that isn’t a multiple).

If the X coordinate values need to include zero, then this is really a math/logic/reasoning question, not a Python or Numpy question.