CODEWITHSUNDEEP
python-3.xscipyodeint
Mohammad Hani
Mohammad Hani

Reputation: 11

Solving simple ODE using scipy odeint gives straight line at 0

I am trying to solve a simple ODE: dN/dt = N*(rho(t)-beta)/lambda

Rho is a function of time and I've generated it using linspace. The code is working for other equations but somehow gives a flat straight line at 0. (You can see it in the graph). Any guidelines about how to correct it?

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

def model2(N, t, rho):
    beta_val = 0.0065
    lambda_val = 0.00002
    k = (rho - beta_val) / lambda_val
    dNdt = k*N
    print(rho)
    return dNdt

# initial condition
N0 = [0]

# number of time points
n = 200

# time points
t = np.linspace(0,200,n)

rho = np.linspace(6,9,n)
#rho =np.array([6,6.1,6.2,6.3,6.4,6.5,6.6,6.7,6.8,6.9,7.0,7.1,7.2,7.3,7.4,7.5,7.6,7.7,7.8,7.9])  # Array of constants

# store solution
NSol = np.empty_like(t)
# record initial conditions
NSol[0] = N0[0]

# solve ODE
for i in range(1,n):
    # span for next time step
    tspan = [t[i-1],t[i]]
    # solve for next step
    N = odeint(model2,N0,tspan,args=(rho[i],))
    print(N)
    # store solution for plotting
    NSol[i] = N[0][0]
    # next initial condition
    #z0 = N0[0]

# plot results
plt.plot(t,rho,'g:',label='rho(t)')
plt.plot(t,NSol,'b-',label='NSol(t)')
plt.ylabel('values')
plt.xlabel('time')
plt.legend(loc='best')
plt.show()

This is the graph I get after running this code

Upvotes: 1

Views: 337

Answers (1)

xdze2
xdze2

Reputation: 4151

I modified your code (and the coefficients) to make it work. When coefficients are also dependent of t, they have to be python functions called by the derivative function:

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

# Define
def model2(N, t, rho):
    beta_val = 0.0065
    lambda_val = 0.02
    k = ( rho(t) - beta_val )/lambda_val
    dNdt = k*N
    return dNdt

def rho(t):
    return .001 + .003/20*t


# Solve
tspan = np.linspace(0, 20, 10)
N0 = .01
N = odeint(model2, N0 , tspan, args=(rho,))

# Plot
plt.plot(tspan, N, label='NS;ol(t)');
plt.ylabel('N');
plt.xlabel('time'); plt.legend(loc='best');

Upvotes: 0

Related Questions