使用序列设置数组元素（Python）

Question

我一直在尝试使用solve_ivp求解一个方程组，但我得到一个错误提示：

ValueError:使用序列设置数组元素

我一直在尝试不同的东西，但没有任何效果，我想我会在这里分享代码，希望有人能帮助我。不要太担心下面geta33中的代码，这是一种正则化方法，我正在尝试实现一阶问题并解决_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)

Answer 1

您可以通过用1D数组替换state0的定义来解决此问题：

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)

祝你好运，你就快到了