Solve_ivp события - PullRequest
       10

Solve_ivp события

0 голосов
/ 30 апреля 2020

Я пытаюсь реализовать событие в моем коде так, чтобы решатель ode останавливался, когда переменная ep становится больше единицы, но я не уверен, как это сделать, поскольку в справочном руководстве scipy не так много литературы , Вот мой текущий вывод и код Вывод из кода.

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

plt.rcParams["font.size"] = 6

plt.autoscale(enable=True, axis='both', tight=None)

#Define Parameters.
n = float(input('Enter the value of n: '))


##############################################################
#I.C
H0 = (20.2,)
H_0 = 20.2
dH0 = 0.3       

#Define the potential
def V_func(t, n):
    return t**n
##############################################################
#Range of Phi.
t_span = np.linspace(start = (3*H_0**2 - 2*dH0**2)**(1/n) , stop = 0 , num = 800) 

##############################################################
#Define H-J equation
def f(t, H, n):
    dh = np.sqrt((3/2)*H**2 - (1/2)*V_func(t ,n))
    return dh

##############################################################
#Find H for each value of phi by solving ODE.
def I_end(t, dh, n):
    return dh
I_end.terminal = True
I_end.direction = -1


sol = integrate.solve_ivp(f, (((3*H_0**2 - 2*dH0**2)**(1/n)),0), H0, t_eval = t_span, method = 'Radau' , rtol = 1e-9, events = I_end, args = (n,),)

H1 = sol.y
t = sol.t

H = H1[0,:]

dH = np.sqrt((3/2)*H**2 - (1/2)*V_func(t ,n))
##############################################################
#New
ep = 2*(dH/H)**2    
##############################################################
#Plot epsilon

print(H)

plt.plot(t, dH)

#plt.yscale('log')

plt.title('$\epsilon$ Vs $\phi$ for n = ' +str(n))

plt.xlabel('$\phi$')

plt.ylabel('$\epsilon$')

plt.show()

Спасибо.

Добро пожаловать на сайт PullRequest, где вы можете задавать вопросы и получать ответы от других членов сообщества.
...