我一直在尝试使用solve_ivp求解方程组，但得到一个错误说：
值错误：使用序列设置数组元素。
我一直在尝试不同的事情，但没有任何效果，我想我会在这里分享代码，希望有人可以帮助我。不要太担心下面的geta33中的代码，这是一种正则化方法，我正在尝试实现问题在于一阶和solve_ivp。
import numpy as np
from scipy import optimize
from numpy import linalg as LA
from scipy.integrate import solve_ivp
G = 1
mu1=3.0616236117028484e+20
mu2=3.0616236117028484e+20
def geta33(Qf1, Qf2, dQf, m2):
Q1, Q2, Q3, Q4 = Qf1
Q5, Q6, Q7, Q8 = Qf2
dQ1, dQ2, dQ3, dQ4 = dQf
q1 = Q1 ** 2 - Q2 ** 2 - Q3 ** 2 + Q4 ** 2
q2 = 2 * Q1 * Q2 - 2 * Q3 * Q4
q3 = 2 * Q1 * Q3 + 2 * Q2 * Q4
q4 = Q5 ** 2 - Q6 ** 2 - Q7 ** 2 + Q8 ** 2
q5 = 2 * Q5 * Q6 - 2 * Q7 * Q8
q6 = 2 * Q5 * Q7 + 2 * Q6 * Q8
qf1 = np.array([q1, q2, q3])
qf2 = np.array([q4, q5, q6])
norm1 = sum((qf1 - qf2) ** 2) ** 0.5
a0 = -G * m2 * (qf1[0] - qf2[0]) / norm1 ** 3
a1 = -G * m2 * (qf1[1] - qf2[1]) / norm1 ** 3
a2 = -G * m2 * (qf1[2] - qf2[2]) / norm1 ** 3
D = 4 * (Q1 ** 2 + Q2 ** 2 + Q3 ** 2 + Q4 ** 2)
r1 = np.sqrt((q1 - mu2) ** 2 + q2 ** 2 + q3 ** 2)
r2 = np.sqrt((q1 + mu1) ** 2 + q2 ** 2 + q3 ** 2)
D = 4 * (Q1 ** 2 + Q2 ** 2 + Q3 ** 2 + Q4 ** 2)
dD = 8 * (Q1 * dQ1 + Q2 * dQ2 + Q3 * dQ3 + Q4 * dQ4)
dq1 = 2 * (Q1 * dQ1 - Q2 * dQ2 - Q3 * dQ3 + Q4 * dQ4)
dq2 = 2 * (Q2 * dQ1 + Q1 * dQ2 - Q4 * dQ3 - Q3 * dQ4)
dq3 = 2 * (Q3 * dQ1 + Q1 * dQ3 + Q4 * dQ2 + Q2 * dQ4)
ddq0 = a0
ddq1 = a1
ddq2 = a2
qpp0 = D ** 2 * (a0 + (dD / D ** 3) * dq1)
qpp1 = D ** 2 * (a1 + (dD / D ** 3) * dq2)
qpp2 = D ** 2 * (a2 + (dD / D ** 3) * dq3)
C1 = 2 * (dQ1 * dQ1 - dQ2 * dQ2 - dQ3 * dQ3 + dQ4 * dQ4)
C2 = 2 * (dQ2 * dQ1 + dQ1 * dQ2 - dQ4 * dQ3 - dQ3 * dQ4)
C3 = 2 * (dQ3 * dQ1 + dQ1 * dQ3 + dQ4 * dQ2 + dQ2 * dQ4)
C4 = 2 * (dQ4 * dQ1 - dQ3 * dQ2 + dQ2 * dQ3 - dQ1 * dQ4)
c = np.array(
[
[2 * Q1, -2 * Q2, -2 * Q3, 2 * Q4],
[2 * Q2, 2 * Q1, -2 * Q4, -2 * Q3],
[2 * Q3, 2 * Q4, 2 * Q1, 2 * Q2],
[2 * Q4, -2 * Q3, 2 * Q2, -2 * Q1],
]
)
b = np.array([qpp0 - C1, qpp1 - C2, qpp2 - C3, -C4])
return [*dQf, *np.linalg.solve(c, b)]
state0 = np.array(
[
np.array([22.33824111, 0.0, 0.0, 0.0]),
np.array([0.03862643, -0.0, 0.0, -0.0]),
np.array([0.00000000e00, -3.50534105e13, 0.00000000e00, 0.00000000e00]),
3.0616327659574474e20,
]
)
me = 3.0616327659574474e20
t = 0
T = 10 ** 7
AU = 499
a = AU
def firstorder(t, state):
pos, vel = state.reshape(2, -1)
return [
vel[0][0],
vel[0][1],
vel[0][2],
vel[0][3],
vel[1],
*geta33(pos[0], pos[1], vel[0], me),
]
sol = solve_ivp(firstorder, [0, T], state0, first_step=1e5, atol=1e-6 * a)
您可以通过将state0的定义替换为 1D 数组来解决此问题：
state0 = np.array([ # We define input as a 1D array
*[ 22.33824111, 0. , 0. , 0. ],
*[ 0.03862643, -0. , 0. , -0. ],
*[ 0.00000000e+00, -3.50534105e+13, 0.00000000e+00, 0.00000000e+00],
3.0616327659574474e+20
])
并在firstorder函数的开头添加一些代码，以您喜欢的格式转换数组：
def firstorder(t, state_vec):
state = np.array([ # We convert the 1D array to an array of arrays
np.array(state_vec[0:4]),
np.array(state_vec[4:8]),
np.array(state_vec[8:12]),
state_vec[12]
])
pos, vel = state.reshape(2,-1)
return [vel[0][0],vel[0][1],vel[0][2],vel[0][3],vel[1], *geta33(pos[0],pos[1],vel[0],me)]
代码正在运行，但您现在遇到溢出问题，在我看来这是另一个问题。你应该在你的程序上工作，以避免使用像3.0e+20这样的非常大的数字，也避免采用这么大的数字的三次幂。这是输出：
% python3 script.py
script.py:23: RuntimeWarning: overflow encountered in double_scalars
a0=-G*m2*(qf1[0]-qf2[0])/norm1**3
script.py:24: RuntimeWarning: overflow encountered in double_scalars
a1=-G*m2*(qf1[1]-qf2[1])/norm1**3
script.py:25: RuntimeWarning: overflow encountered in double_scalars
a2=-G*m2*(qf1[2]-qf2[2])/norm1**3
script.py:37: RuntimeWarning: overflow encountered in double_scalars
qpp0=D**2*(a0+(dD/D**3)*dq1)
script.py:38: RuntimeWarning: overflow encountered in double_scalars
qpp1=D**2*(a1+(dD/D**3)*dq2)
script.py:39: RuntimeWarning: overflow encountered in double_scalars
qpp2=D**2*(a2+(dD/D**3)*dq3)
script.py:22: RuntimeWarning: overflow encountered in square
norm1=sum( (qf1-qf2)**2 )**0.5
script.py:27: RuntimeWarning: overflow encountered in double_scalars
r1=np.sqrt((q1-mu2)**2+q2**2+q3**2)
script.py:28: RuntimeWarning: overflow encountered in double_scalars
r2=np.sqrt((q1+mu1)**2+q2**2+q3**2)
script.py:37: RuntimeWarning: invalid value encountered in double_scalars
qpp0=D**2*(a0+(dD/D**3)*dq1)
script.py:38: RuntimeWarning: invalid value encountered in double_scalars
qpp1=D**2*(a1+(dD/D**3)*dq2)
script.py:39: RuntimeWarning: invalid value encountered in double_scalars
qpp2=D**2*(a2+(dD/D**3)*dq3)
祝你好运，你快到了！