试图让ODEint为具有时变参数的耦合ODE工作

Question

我试着用ODEint来解大爆炸核合成过程中ODE控制核素密度的耦合集。我从一个非常简单的版本开始，只是根据单个反应n+p跟踪中子、质子和氘的密度-->；D、 这需要一组耦合的3个方程，所有这些方程都需要随时间变化的参数（反应速率）。更复杂的版本将有多达20个反应速率参数（都随时间而变化），这些参数必须传递到ODE中，但我在这里将其简化为三个函数，用于测试。不幸的是代码（如下）运行，但它给我的核素密度是垃圾。我对编码，尤其是python还是新手，非常感谢您的帮助

它似乎完全集成，但输出垃圾，并发出运行时警告：

…ODEintWarning:在此调用上完成了过多的工作（可能是错误的Dfun类型）。以full_output=1运行以获取定量信息。 警告。警告（警告消息，ODEintWarning）

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Fri May  3 12:42:24 2019
"""

import scipy.integrate as integ
import numpy as np
import matplotlib.pyplot as plt


#Initate time array
n=200
t0=0.1
tf=100
time_vec = np.linspace(t0,tf,n)

#Initiate Nuclide Abundances 
X_n0 = 0.5
X_p0 = 0.5
X_D0 = 0.0
X0_vec = (X_n0, X_p0, X_D0)

#define functions of changing parameters
def a(t):
    return np.sqrt(t)

def rho_b(t, a):
    return 1 / (a(t) ** 3)

#define functions of changing reaction rates
def brac_pn(t, rho_b):
    return 2.5E4 * rho_b(t, a)


#Define function to plug into ODEint
def EvolveDensity(X, t, a, rho_b, brac_pn):
    dX_n_dt = (X[0] * X[1]) + (brac_pn(t, rho_b) * X[2] )
    dX_p_dt = (X[0] * X[1]) + (brac_pn(t, rho_b) * X[2] )
    dX_D_dt = - (X[0] * X[1]) + (brac_pn(t, rho_b) * X[2] )
    return [dX_n_dt, dX_p_dt, dX_D_dt]


#Solve Ode through time
Densities = integ.odeint(EvolveDensity, X0_vec, time_vec, args=(a, rho_b, brac_pn))


plt.plot(time_vec,Densities[:,0])
plt.plot(time_vec,Densities[:,1])
plt.plot(time_vec,Densities[:,2])