Scipy的solve\u ivp给出了“ValueError:用序列设置数组元素”

Question

试图用solve_ivp作图来解一个两体系统的运动方程。然而，我不断得到ValueError: setting an array element with a sequence。你知道吗

下面是代码，抱歉，如果实际函数有点小，那么问题就出现在底部使用solve_ivp时，所以我不知道是否应该以不同的方式重塑输入数组。你知道吗

import numpy as np
from numpy import linalg as LA
import scipy.integrate
from scipy.integrate import solve_ivp
from scipy import optimize

AU=1.5e11
a=AU
e=0
ms = 2E30
me = 5.98E24
yr=3.15e7
h=100
mu1=ms*me/(ms+me)
mu2=ms*me/(ms+me)
G=6.67E11

vi=np.sqrt(G*ms*(2/(a*(1-e))-1/a))

sunpos=np.array([0,0,0])
earthpos=np.array([a*(1-e),0,0])

earthv=np.array([0,vi,0])
sunv=np.array([0,0,0])



def equations(p,qf,m):
    q1, q2, q3, q4 = p
    r1=np.sqrt((qf[0]-m)**2+qf[1]**2+qf[2]**2)
    return q1**2-q2**2-q3**2+q4**2-qf[0]+m,\
     2*q1*q2-2*q3*q4-qf[1],\
     2*q1*q3+2*q2*q4-qf[2],\
     q1*q4-q2*q3+q3*q2-q4*q1



def equationsp(pp,qfp,qs):
    dq1, dq2, dq3, dq4 = pp
    D=qs[0]**2+qs[1]**2+qs[2]**2+qs[3]**2
    return 2*(qs[0]*dq1-qs[1]*dq2-qs[2]*dq3+qs[3]*dq4)-D*qfp[0],\
           2*(qs[1]*dq1+qs[0]*dq2-qs[3]*dq3-qs[2]*dq4)-D*qfp[1],\
           2*(qs[2]*dq1+qs[0]*dq3+qs[3]*dq2+qs[1]*dq4)-D*qfp[2],\
           2*(qs[3]*dq1-qs[2]*dq2+qs[1]*dq3-qs[0]*dq4)




def DotR(ppp,qf,dqf,m):
         ddQ1,ddQ2,ddQ3,ddQ4=ppp
         r1=np.sqrt((qf[0]-mu2)**2+qf[1]**2+qf[2]**2)
         r2=np.sqrt((qf[0]+mu1)**2+qf[1]**2+qf[2]**2)
         Q1,Q2,Q3,Q4=optimize.root(equations, (qf[0],qf[1],qf[2],1),\
                                   (qf,m)).x
#         qf[0]=Q1**2-Q2**2-Q3**2+Q4**2+mu2
#         qf[1]=2*Q1*Q2-2*Q3*Q4
#         qf[2]=2*Q1*Q3+2*Q2*Q4
         r1=Q1**2+Q2**2+Q3**2+Q4**2 
         D=4*(Q1**2+Q2**2+Q3**2+Q4**2)
         dQ1,dQ2,dQ3,dQ4=optimize.root(equationsp, \
         (dqf[0],dqf[1],dqf[2],1),(dqf,(Q1, Q2, Q3, Q4))).x
#         dq0=2*(Q1*dQ1-Q2*dQ2-Q3*dQ3+Q4*dQ4)
#         dq1=2*(Q2*dQ1+Q1*dQ2-Q4*dQ3-Q3*dQ4)
         dD=8*(Q1*dQ1+Q2*dQ2+Q3*dQ3+Q4*dQ4)
         #ddQ1,ddQ2,ddQ3,ddQ4=diff(dQ1)/diff(t),diff(dQ2)/diff(t),diff(dQ3)/diff(t),diff(dQ4)/diff(t)
         #ddq[0]=2*(dQ1*dQ1-dQ2*dQ2-dQ3*dQ3+dQ4*dQ4)+2*(Q1*ddQ1-Q2*ddQ2-Q3*ddQ3+Q4*ddQ4)
         #ddq[1]=2*(Q2*ddQ1+Q1*ddQ2-Q4*ddQ3-Q3*ddQ4)+2*(dQ2*dQ1+dQ1*dQ2-dQ4*dQ3-dQ3*dQ4)
         #ddq[2]=2*(Q3*ddQ1+Q1*ddQ3+Q4*ddQ2+Q2*ddQ4)+2*(dQ3*dQ1+dQ1*dQ3+dQ4*dQ2+dQ2*dQ4)
         ddq0=(D**3*(qf[0]-mu1/r1**3*(qf[0]-mu2)-mu2/r2**3*(qf[0]+mu1))+dD*dqf[0]*D+2*D**2*dqf[1]*D)/D
         ddq1=(D**3*(qf[1]-mu1*qf[1]/r1**3-mu2*qf[1]/r2**3)+dD*dqf[1]*D-2*D**2*dqf[0]*D)/D
         ddq2=(D**3*(-mu1*qf[2]/r1**3-mu2*qf[2]/r2**3))/D

   #      print(D,Q1,Q2,Q3,Q4)
         return 2*(dQ1*dQ1-dQ2*dQ2-dQ3*dQ3+dQ4*dQ4)+2*(Q1*ddQ1-Q2*ddQ2-Q3*ddQ3+Q4*ddQ4)-ddq0,\
                2*(Q2*ddQ1+Q1*ddQ2-Q4*ddQ3-Q3*ddQ4)+2*(dQ2*dQ1+dQ1*dQ2-dQ4*dQ3-dQ3*dQ4)-ddq1,\
                2*(Q3*ddQ1+Q1*ddQ3+Q4*ddQ2+Q2*ddQ4)+2*(dQ3*dQ1+dQ1*dQ3+dQ4*dQ2+dQ2*dQ4)-ddq2,\
                2*(Q4*ddQ1-Q3*ddQ2+Q2*ddQ3-Q1*ddQ4)+2*(dQ4*dQ1-dQ3*dQ2+dQ2*dQ3-dQ1*dQ4)



qna=np.concatenate((earthpos,earthv),axis=0)
qn=np.reshape(qna,(1,6))
qnn=qn[0,:]
new=solve_ivp(DotR,(0,10**5),(qnn,qnn,mu2))

Answer 1

1247:~/mypy$ python3 stack57409204.py 
Traceback (most recent call last):
  File "stack57409204.py", line 84, in <module>
    new=solve_ivp(DotR,(0,10**5),(qnn,qnn,mu2))
  File "/usr/local/lib/python3.6/dist-packages/scipy/integrate/_ivp/ivp.py", line 477, in solve_ivp
    solver = method(fun, t0, y0, tf, vectorized=vectorized, **options)
  File "/usr/local/lib/python3.6/dist-packages/scipy/integrate/_ivp/rk.py", line 96, in __init__
    support_complex=True)
  File "/usr/local/lib/python3.6/dist-packages/scipy/integrate/_ivp/base.py", line 120, in __init__
    self._fun, self.y = check_arguments(fun, y0, support_complex)
  File "/usr/local/lib/python3.6/dist-packages/scipy/integrate/_ivp/base.py", line 15, in check_arguments
    y0 = y0.astype(dtype, copy=False)
ValueError: setting an array element with a sequence.

添加打印

qnn [1.50000000e+11 0.00000000e+00 0.00000000e+00 0.00000000e+00
 2.98216923e+15 0.00000000e+00]
mu2 5.979982119853463e+24

我经常通过使用示例输入运行fn来测试这样的代码。但在你的案子里我很难搞清楚。solve_ivp说

The calling signature is fun(t, y). 
t - scalar
y - (n,), same shape as the output

但是你的

DotR(ppp,qf,dqf,m)
    ddQ1,ddQ2,ddQ3,ddQ4=ppp

标量t在哪里，最初是0？y的维数是多少，返回值是多少？(4,)? 你知道吗