西皮的solve_ivp给了"ValueError: setting an array element with a sequence"

Question

尝试使用solve_ivp绘制图解 2 个体系统运动方程。但是，我一直在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，(？