Попытка заставить ODEint работать для связанного ODE с изменяющимися во времени параметрами - PullRequest
Новая
       9

Попытка заставить ODEint работать для связанного ODE с изменяющимися во времени параметрами

0 голосов
/ 09 мая 2019

Я пытаюсь использовать ODEint для решения связанного набора плотностей нуклидов в ODE во время нуклеосинтеза Большого взрыва. Я начинаю с очень простой версии, просто отслеживающей плотности нейтронов, протонов и дейтерия на основе одной реакции n + p -> D, для которой требуется объединенный набор из 3 уравнений, для каждого из которых требуются параметры (скорости реакции), которые меняются во времени. Более сложная версия будет иметь до 20 параметров скорости реакции (все меняющиеся со временем), которые нужно будет передать в ODE, но я упростил это до трех функций здесь, для тестирования. К сожалению, код (ниже) выполняется, но плотность нуклидов, которую он мне дает, - мусор. Я все еще новичок в программировании и Python в частности, и помощь будет очень признательна.

Кажется, что он полностью интегрируется, но выводит мусор и выдает предупреждение во время выполнения:

... ODEintWarning: избыточная работа, выполненная для этого вызова (возможно, неправильный тип Dfun). Запустите с full_output = 1, чтобы получить количественную информацию. warnings.warn (warning_msg, 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])
...