使用odeint求解不同阶数的多极子ode

Question

我正在尝试解决三个问题，其中一个是一阶，其他是二阶。 这是方程式：

这是代码：

#note that:
p[0]=y
p[1]=dy/dt
p[2]=q
p[3]=u
p[4]=du/dt

#initial values and requirements :

import numpy as np
import time
from scipy.integrate import odeint

g=1e-25
k=4.14*(10**-5)
t=(10**31)*np.logspace(0.247237,2.35443,10**7)
t=t.tolist()
z0=[2.435*1e26,0,1/3400]

#equations:

def my_eqs1(p,t):
     return[p[1],-2*p[2]*1.4441*(10**-33)*np.sqrt(0.31*(p[2]**-4)*(p[2]+1/3400)+0.69)*p[1]-
       (r*p[2])**(2)*p[0],p[2]*p[2]*1.4441*(10**-33)*np.sqrt(0.31*(p[2]**-4)*(p[2]+1/3400)+0.69)]
       
       

#solving:
start_time=time.perf_counter()
r=1e-29
z1=odeint(my_eqs1,z0,t)
end_time=time.perf_counter()
print(end_time-start_time,"seconds")

在下一步中，我必须使用

p[1]

来解决另一个二阶颂歌，所以我尝试将

my_eqs1

重新定义为休闲以同时解决所有这些问题：

def my_eqs2(p,t):
    return[p[1],-2*p[2]*1.4441*(10**-2)*np.sqrt(0.31*(p[2]**-4)*(p[2]+1/3400)+0.69)*p[1]-
           ((1e31*r*p[2])**2)*p[0],p[2]*p[2]*1.4441*(10**-2)*np.sqrt(0.31*(p[2]**-4)*(p[2]+1/3400)+0.69)
           ,1e-62*p[4],-((k**2)+(g*p[1]/k))*p[3]]

和新的

z0

：

z0=[2.435*1e26,0,1/3400,1/np.sqrt(2*k),0.35927359523]

但是这次

odeint

给出的答案一点都不好。

这段代码有什么问题？

通过重新定义函数和缩放

t

，（如t->1e31*t）我得到了一些答案：

t=np.logspace(0.247237,2.35443,10**7)

def my_eqs3(p,t):
        return[p[1],-2*p[2]*1.4441*(10**-2)*np.sqrt(0.31*(p[2]**-4)*(p[2]+1/3400)+0.69)*p[1]-
               ((1e31*r*p[2])**2)*p[0],p[2]*p[2]*1.4441*(10**-2)*np.sqrt(0.31*(p[2]**-4)*(p[2]+1/3400)+0.69)
               ,1e-62*p[4],-((k**2)+(g*p[1]/k))*p[3]]