当我尝试做数学时，Python 重新启动了 shell答案

【问题标题】：Python restarts the shell when I tried to do math当我尝试做数学时，Python 重新启动了 shell
【发布时间】：2019-01-11 16:06:13
【问题描述】：

我对以下代码有疑问。我尝试了很多网页来找到一些解决方案，但我无法做到。问题是当我在打印 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

【问题讨论】：

用 30 行数据或其他东西创建一个不同的 new.txt 副本。如果它加载，您可以确认我的怀疑：这很可能是您没有足够的内存来保存所有数据。只是一个典型的内存问题。
@ParitoshSingh 是的，在添加 4 GB 额外的 ddr3 RAM 后，我尝试了一些先前的计算有这个问题，我的问题得到了解决。但现在我有 8 GB RAM。但现在让我来做你的建议。等一下……
@ParitoshSingh 不，文本文件又有 17 行我有同样的问题
您能否提供仍然允许复制问题的最小输入？
@ParitoshSingh 什么意思，我听不懂

标签： python shell kernel jupyter ode

【解决方案1】：

好的，所以问题出在odeint。 docs 建议改用 solve_ivp。

现在完全免责声明，我完全不知道这意味着什么，事物的数学和含义超出了我的范围。但是，我尝试使用solve_ivp 最好地模仿odeint 行为，它开箱即用不接受args。展示，ghetto lambda 来救援。

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

请注意，这不是完全替代，solve 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

您提供的代码并不总是产生负值。我想不明白您为什么期望或希望他们这样做，但结果应该与原始代码匹配。

【讨论】：

在做任何事情之前，我必须先问一个问题。如果使用def dec()函数，则打印结果为positive，而全部为负数。当然以这种方式我不知道最终的答案。我的意思是他们是错误的结果
如果结果错误，那么它们是您提供的代码的结果。你必须检查方程的逻辑等等。检查更新，看看是否给你任何线索。 @MaY 恐怕我无法帮助您了解所有数学背后的逻辑。
不，不，代码是正确的，即使在之前打印了 5 行的代码中，我们也可以看到负值。但是这段代码的答案完全是肯定的。此代码用于许多先前的工作。
我已经发布了新代码的答案。它也有负值。所以我不确定你指的是什么。
别担心，很高兴我能帮上忙。

【解决方案2】：

您的代码只需稍作改动即可按预期工作：避免零长度的积分区间，如果遇到这种情况，直接返回初始值。

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

【讨论】：