我正在尝试用 RK4 方法求解非线性微分方程组。我在该方程式中的一个参数是一个非常小的数字，大约为 10^{-121} 因此为了更精确地工作，我使用 Decimal 包。但我被“decimal.InvalidOperation”警告困住了。
def V(x,y,N):
v = np.divide(y**2, (1-x**2-y**2)) * 3* H02 * (Om1 * np.exp(-3*N) + Om2 * np.exp(-4*N))
return v
def f1 (x,y,l,q):
return -3*x + l * np.sqrt(3/2) * y**2+ 3/2 * x * ( 2*x**2 + q*(1 - x**2 y**2))
def f2 (x,y,l,q):
return - l * np.sqrt(3/2) * y * x + 3/2 * y * (2 * x**2 + q*(1 - x**2 - y**2))
def f3 (N, x, y, l, g, M, al):
exp1, G, m = decimal.Decimal(-2*g).exp(), decimal.Decimal(g), decimal.Decimal(M)
xx, a, ll = decimal.Decimal(x), decimal.Decimal(al), decimal.Decimal(l)
VV = decimal.Decimal(V(x,y,N))
log = decimal.Decimal(VV/(m*exp1) + 1).ln()
Tanh = 1/G * log - 1
# gamm = -(G*(1 - Tanh**2) - 2 * Tanh)/(np.sqrt(6*a)*ll)
return float(ll *xx/np.sqrt(a)* (G*(1 - Tanh**2) - 2 * Tanh + np.sqrt(6*a)*ll))
def R_K_4(xi,yi,li, g, M, al):
N = int(round((lnaf-lnai)/dlna))
lna = np.linspace(lnai, lnaf, N)
x = np.zeros(N)
y = np.zeros(N)
l = np.zeros(N)
x[0],y[0],l[0] = xi,yi,li
for i in range(0,N-1):
q = 1 + 1/(3*(1 + 5780*np.exp(lna[i])))
kx1 = f1( x[i], y[i], l[i], q )
ky1 = f2( x[i], y[i], l[i], q )
kl1 = f3( lna[i], x[i], y[i], l[i], g, M, al)
kx2 = f1( x[i]+kx1*dlna/2, y[i]+ky1*dlna/2, l[i]+kl1*dlna/2, q )
ky2 = f2( x[i]+kx1*dlna/2, y[i]+ky1*dlna/2, l[i]+kl1*dlna/2, q )
kl2 = f3( lna[i] + dlna/2, x[i]+kx1*dlna/2, y[i]+ky1*dlna/2, l[i]+kl1*dlna/2, g, M, al )
kx3 = f1( x[i]+kx2*dlna/2, y[i]+ky2*dlna/2, l[i]+kl2*dlna/2, q )
ky3 = f2( x[i]+kx2*dlna/2, y[i]+ky2*dlna/2, l[i]+kl2*dlna/2, q )
kl3 = f3( lna[i] + dlna/2, x[i]+kx2*dlna/2, y[i]+ky2*dlna/2, l[i]+kl2*dlna/2, g, M, al )
kx4 = f1( x[i]+kx3*dlna, y[i]+ky3*dlna, l[i]+kl3*dlna, q )
ky4 = f2( x[i]+kx3*dlna, y[i]+ky3*dlna, l[i]+kl3*dlna, q )
kl4 = f3( lna[i] + dlna, x[i]+kx3*dlna, y[i]+ky3*dlna, l[i]+kl3*dlna, g, M, al )
x[i+1] = x[i] +dlna/6*(kx1 + 2*kx2 + 2*kx3 + kx4)
y[i+1] = y[i] +dlna/6*(ky1 + 2*ky2 + 2*ky3 + ky4)
l[i+1] = l[i] +dlna/6*(kl1 + 2*kl2 + 2*kl3 + kl4)
return x,y,l
当我将求解器称为
H02 = 1.44*1e-122
Om1, Om2 = 0.31, 5.45*1e-5
lnai, lnaf, dlna = -15, 10, 1e-2
x11, y11, l11 = R_K_4(0, 2.4e-11, -0.919, 128, 1e-10, 7/3)
我收到上述错误。但我认为错误发生在该步骤之前，但我无法确定到底哪里出了问题。请帮助我解决这个问题！