首页
学习
活动
专区
圈层
工具
发布
首页
学习
活动
专区
圈层
工具
MCP广场
社区首页 >问答首页 >当我尝试做数学时,Python重新启动shell。

当我尝试做数学时,Python重新启动shell。
EN

Stack Overflow用户
提问于 2019-01-11 16:06:13
回答 2查看 123关注 0票数 0

下面的代码有问题。我尝试了很多网页,以找到一些解决方案,但我无法设法。问题是,在只打印4行之后运行此脚本时,Python将重新启动shell。在jupyter中,它发送消息:The kernel appears to have died. It will restart automatically.

我尝试在包含(300,000)行数据的脚本的第一个位置读取一个文件。然后,在计算ODE和其他函数之后,我希望打印结果进行比较。

代码语言:javascript
运行
复制
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)可以是

代码语言:javascript
运行
复制
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
EN

回答 2

Stack Overflow用户

回答已采纳

发布于 2019-01-11 17:28:34

好的,所以问题在于odeint文档建议改用ivp

现在完全的免责声明,我完全不知道这意味着什么,数学和意义的东西是我无法理解的。然而,我试图最好地模仿solve_ivp的行为,因为它不接受args。向救援队致敬。

代码语言:javascript
运行
复制
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范围的最小代码,您可以看到下面的战场。

代码语言:javascript
运行
复制
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,下面是代码的更新方式。

代码语言:javascript
运行
复制
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])

上述数据的输出如下所示:

代码语言:javascript
运行
复制
[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个值打印。

代码语言:javascript
运行
复制
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])

以及产出:

代码语言:javascript
运行
复制
[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

您提供的代码并不总是生成负值。我不认为您理解您为什么期望或希望它们这样做,但是结果应该与原始代码相匹配。

票数 1
EN

Stack Overflow用户

发布于 2019-01-11 22:28:03

您的代码将像预期的那样工作,只需做一个小的更改:避免零长度的集成间隔(如果遇到这种情况),直接返回初始值。

代码语言:javascript
运行
复制
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   
票数 1
EN
页面原文内容由Stack Overflow提供。腾讯云小微IT领域专用引擎提供翻译支持
原文链接:

https://stackoverflow.com/questions/54150112

复制
相关文章

相似问题

领券
问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档