Python перезапускает оболочку, когда я пытался сделать математику

Question

У меня проблема со следующим кодом. Я пробовал много веб-страниц, чтобы найти какие-то решения, но мне не удалось. Проблема в том, что когда я запускаю этот скрипт после 4 строк печати, Python перезапускает оболочку. В jupyter отправляется сообщение: The kernel appears to have died. It will restart automatically.

Я пытаюсь прочитать файл первого скрипта, содержащий (300,000) строки данных. затем после вычисления ODE и других функций я ожидаю распечатать результаты для сравнения

import numpy as np
from scipy.integrate import odeint
from math import *
from scipy.integrate import quad  
import math
import pandas as pd

hr, dr, cr, br = np.genfromtxt('outputs/new.txt',unpack=True)

def OD_H(od, z, c, b):
   Omegai = (1-od)
   div1 = np.divide(1, (1 + z), where = (1 + z)!= 0)     
   dMdt = -(div1) * (2 *(1-od)* (-2 + (od/(6 * c))) + 3 - 3 * b**2 * Omegai - 3 * od)
return dMdt

def ant(H0, z, od0, c, b):
    z1 = 0
    od = odeint(OD_H, od0, [z1, z], args=(c, b))[-1]                   
    return od   

def dec(H0, z, od0, c, b):
    od = ant(H0, z, od0, c, b)
    q = -1 - (-2 + od/(6 * c))
    return q

for i in range(len(hr)):
    for z in range (0,1):
        print(dec(hr[i], z, dr[i], cr[i], br[i]),hr[i], dr[i], cr[i], br[i])

Это простой код, но я не знаю, что это за проблема. Я действительно ценю любую помощь.

входной файл (new.txt) может быть

71.076588184266 0.40147988209522 0.080396967668756 0.050302016457046
71.02284157687 0.39756707964421 0.080918035449145 0.050501956013259
71.102923163306 0.41587392748136 0.07823452108922 0.049336707395359
70.860444589498 0.46748446539443 0.072392464271658 0.046667808684486
70.181278149341 0.44888833570037 0.077917371645449 0.04777288009128
70.588452406351 0.49035265611716 0.072303154996487 0.045942096884044
70.588452406351 0.49035265611716 0.072303154996487 0.045942096884044
70.011812869146 0.44210637315163 0.07871914246357 0.048393990901086
69.807956729005 0.41349634394633 0.082020266421564 0.049900569076028
69.807956729005 0.41349634394633 0.082020266421564 0.049900569076028
69.807956729005 0.41349634394633 0.082020266421564 0.049900569076028
70.123419349447 0.43961350279409 0.07862300319627 0.048607832896286
70.361666430312 0.41397677666087 0.080502527828865 0.049745843125116
70.361666430312 0.41397677666087 0.080502527828865 0.049745843125116
70.357430153315 0.41485946940097 0.08042642593323 0.049703105696664
70.357430153315 0.41485946940097 0.08042642593323 0.049703105696664
71.551080656041 0.51047305096688 0.066682530098241 0.0446474321235

Answer 1

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

def ant(H0, z, od0, c, b):
    z1 = 0
    if type(od0) is np.float64: od0 = np.array([od0]); # for uniform output
    od = od0 if abs(z-z1) < 1e-15 else odeint(OD_H, od0, [z1, z], args=(c, b))[-1]                   
    return od   

Answer 2

Хорошо, так что проблема с odeint. Вместо docs рекомендуется переключиться на solve_ivp .

Теперь полный отказ от ответственности, я абсолютно не знаю, что это значит, математика и значения вещей были за пределами меня. Тем не менее, я попытался лучше всего подражать поведению odeint с solve_ivp, который из коробки не принимает аргументы. Представляя, спасение гетто лямбда.

od = odeint(OD_H, od0, [z1, z], args=(c, b))[-1] #before
od = solve_ivp(lambda od, z: OD_H(od, z, c, b), t_span = [z1, z], y0 = [od0])['y'][-1][-1] #after

обратите внимание, что это не ПОЛНОСТЬЮ замена, результат решения ivp представляет плавающее число таким образом, и вы захотите обернуть его в [od] вместо точного соответствия старому результату. Что касается минимального кода, который я использовал, чтобы сузить odeint, вы можете увидеть поле битвы ниже.

import numpy as np
from scipy.integrate import odeint, solve_ivp
#from math import *
from scipy.integrate import quad  
#import math
#import pandas as pd

hr, dr, cr, br = np.genfromtxt('new.txt',unpack=True)




def OD_H(od, z, c, b):
   Omegai = (1-od)
   div1 = np.divide(1, (1 + z), where = (1 + z)!= 0)     
   dMdt = -(div1) * (2 *(1-od)* (-2 + (od/(6 * c))) + 3 - 3 * b**2 * Omegai - 3 * od)
   return dMdt

def ant(H0, z, od0, c, b):
    z1 = 0
    #od = odeint(OD_H, od0, [z1, z], args=(c, b))[-1] #this solver crashed
    od = solve_ivp(lambda od, z: OD_H(od, z, c, b), t_span = [z1, z], y0 = [od0])['y'][-1][-1] #this worked out. perhaps wrap in square brackets [od] if needed.
    #od = OD_H(od0,z,c,b)  #this alone without the solvers worked fine              
    return od   

#def dec(H0, z, od0, c, b): #remove the middleman
#    od = ant(H0, z, od0, c, b)
#    q = -1 - (-2 + od/(6 * c))
#    return od

for i in range(5): #simple check instead
    z = 0 #you do not need a loop here
    res = ant(hr[i], z, dr[i], cr[i], br[i])
    print(res)
        #print(dec(hr[i], z, dr[i], cr[i], br[i]),hr[i], dr[i], cr[i], br[i]) print was not the culprit

Суть в том, что я не совсем знаю (пока), почему odeint разбился таким образом.

EDIT:

Для операции, вот как будет выглядеть обновленный код.

import numpy as np
from scipy.integrate import odeint, solve_ivp
#from math import *
from scipy.integrate import quad  
#import math
#import pandas as pd

hr, dr, cr, br = np.genfromtxt('new.txt',unpack=True)

def OD_H(od, z, c, b):
   Omegai = (1-od)
   div1 = np.divide(1, (1 + z), where = (1 + z)!= 0)     
   dMdt = -(div1) * (2 *(1-od)* (-2 + (od/(6 * c))) + 3 - 3 * b**2 * Omegai - 3 * od)
   return dMdt

def ant(H0, z, od0, c, b):
    z1 = 0
    #od = odeint(OD_H, od0, [z1, z], args=(c, b))[-1]  
    od = [solve_ivp(lambda od, z: OD_H(od, z, c, b), t_span = [z1, z], y0 = [od0])['y'][-1][-1]]
    return od   

def dec(H0, z, od0, c, b):
    od = ant(H0, z, od0, c, b)
    q = -1 - (-2 + od/(6 * c))
    return q

for i in range(len(hr)):
    for z in range (0,1):
        print(dec(hr[i], z, dr[i], cr[i], br[i]),hr[i], dr[i], cr[i], br[i])

Вывод вышеуказанных данных выглядит следующим образом:

[0.16771346] 71.076588184266 0.40147988209522 0.080396967668756 0.050302016457046
[0.18113212] 71.02284157687 0.39756707964421 0.080918035449145 0.050501956013259
[0.11404428] 71.102923163306 0.41587392748136 0.07823452108922 0.049336707395359
[-0.07627332] 70.860444589498 0.46748446539443 0.072392464271658 0.046667808684486
[0.03981973] 70.181278149341 0.44888833570037 0.077917371645449 0.04777288009128
[-0.13031641] 70.588452406351 0.49035265611716 0.072303154996487 0.045942096884044
[-0.13031641] 70.588452406351 0.49035265611716 0.072303154996487 0.045942096884044
[0.06395836] 70.011812869146 0.44210637315163 0.07871914246357 0.048393990901086
[0.15976794] 69.807956729005 0.41349634394633 0.082020266421564 0.049900569076028
[0.15976794] 69.807956729005 0.41349634394633 0.082020266421564 0.049900569076028
[0.15976794] 69.807956729005 0.41349634394633 0.082020266421564 0.049900569076028
[0.06809821] 70.123419349447 0.43961350279409 0.07862300319627 0.048607832896286
[0.14293214] 70.361666430312 0.41397677666087 0.080502527828865 0.049745843125116
[0.14293214] 70.361666430312 0.41397677666087 0.080502527828865 0.049745843125116
[0.14029196] 70.357430153315 0.41485946940097 0.08042642593323 0.049703105696664
[0.14029196] 70.357430153315 0.41485946940097 0.08042642593323 0.049703105696664
[-0.27587903] 71.551080656041 0.51047305096688 0.066682530098241 0.0446474321235

Код, который я могу использовать для сравнения, таков. Обратите внимание на старый решатель odeint, но только с 4 напечатанными значениями.

import numpy as np
from scipy.integrate import odeint, solve_ivp
#from math import *
from scipy.integrate import quad  
#import math
#import pandas as pd

hr, dr, cr, br = np.genfromtxt('new.txt',unpack=True)

def OD_H(od, z, c, b):
   Omegai = (1-od)
   div1 = np.divide(1, (1 + z), where = (1 + z)!= 0)     
   dMdt = -(div1) * (2 *(1-od)* (-2 + (od/(6 * c))) + 3 - 3 * b**2 * Omegai - 3 * od)
   return dMdt

def ant(H0, z, od0, c, b):
    z1 = 0
    od = odeint(OD_H, od0, [z1, z], args=(c, b))[-1]  
    #od = [solve_ivp(lambda od, z: OD_H(od, z, c, b), t_span = [z1, z], y0 = [od0])['y'][-1][-1]]
    return od   

def dec(H0, z, od0, c, b):
    od = ant(H0, z, od0, c, b)
    q = -1 - (-2 + od/(6 * c))
    return q

for i in range(4): #simplified to avoid crash
    for z in range (0,1):
        print(dec(hr[i], z, dr[i], cr[i], br[i]),hr[i], dr[i], cr[i], br[i])

И вывод:

[0.16771346] 71.076588184266 0.40147988209522 0.080396967668756 0.050302016457046
[0.18113212] 71.02284157687 0.39756707964421 0.080918035449145 0.050501956013259
[0.11404428] 71.102923163306 0.41587392748136 0.07823452108922 0.049336707395359
[-0.07627332] 70.860444589498 0.46748446539443 0.072392464271658 0.046667808684486

Код, который вы указали, не всегда генерирует отрицательные значения. Я не берусь понять, почему вы ожидаете или хотите, чтобы они это делали, но результаты должны совпадать с исходным кодом.