下面的代码有问题。我尝试了很多网页,以找到一些解决方案,但我无法设法。问题是,在只打印4行之后运行此脚本时,Python将重新启动shell。在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发布于 2019-01-11 17:28:34
现在完全的免责声明,我完全不知道这意味着什么,数学和意义的东西是我无法理解的。然而,我试图最好地模仿solve_ivp的行为,因为它不接受args。向救援队致敬。
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为什么会这样崩溃。
编辑:
对于op,下面是代码的更新方式。
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您提供的代码并不总是生成负值。我不认为您理解您为什么期望或希望它们这样做,但是结果应该与原始代码相匹配。
发布于 2019-01-11 22:28:03
您的代码将像预期的那样工作,只需做一个小的更改:避免零长度的集成间隔(如果遇到这种情况),直接返回初始值。
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 https://stackoverflow.com/questions/54150112
复制相似问题