Решение простого ODE с использованием scipy odeint дает прямую линию в 0

Я пытаюсь решить простую ОДУ: dN/dt = N*(rho(t)-бета)/лямбда

Rho — это функция времени, и я сгенерировал ее с помощью linspace. Код работает для других уравнений, но каким-то образом дает ровную прямую линию в точке 0. (Вы можете видеть это на графике). Любые рекомендации о том, как это исправить?

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

Это график, который я получаю после запуска этого кода


python-3.x scipy odeint
person Mohammad Hani    schedule 15.08.2018    source источник

Ответы (1)


arrow_upward
0
arrow_downward

Я изменил ваш код (и коэффициенты), чтобы он работал. Когда коэффициенты также зависят от t, они должны быть функциями Python, вызываемыми производной функцией:

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');
person xdze2    schedule 16.08.2018
comment
Это работает отлично, но я хочу использовать rho из диапазона, который я определил в своем коде, то есть 6-9, и я не смог найти способ его использовать. - person Mohammad Hani; 16.08.2018
comment
Я думаю, что ошибка возникает из-за используемых значений, решение dN/dt = k*N с константой k равно exp(k*t), а если rho ~ 6, k составляет около 3e6, таким образом, N(200)~ exp(1e8)... - person xdze2; 16.08.2018