优化根会收敛吗？

Question

尝试使用ks正则化和runge-kutte四阶方法绘制双体轨道。问题似乎是DotR函数中的optimize.root似乎不收敛，我只得到一条直线。在下面的代码中，我在第一个根查找之后添加了一个exit（），并且我还打印了等式函数的返回，以确定是否实际找到了零

我不知道我是否应该继续尝试根查找的不同起始值，或者是否有其他方法更适合这个问题

import numpy as np
from scipy import optimize
from numpy import linalg as LA
import matplotlib.pyplot as plt

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

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])

norme=sum( (earthpos-sunpos)**2 )**0.5
norma=sum( (sunpos-earthpos)**2 )**0.5

def equations(p,qf,m):
    q1, q2, q3, q4 = p
    r1=np.sqrt((qf[0]-m)**2+qf[1]**2+qf[2]**2)
#    print(q1,q2,q3,q4,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],\
     q4*q1-q2*q3+q3*q2-q1*q4)
    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],\
     q4*q1-q2*q3+q3*q2-q1*q4


def jaceq(p,qf):
    q1, q2, q3, q4 = p
    return 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],[q4,-q3,q2,-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 jaceqp(pp,qfp,qs):
    dq1, dq2, dq3, dq4 = pp
    D=qs[0]**2+qs[1]**2+qs[2]**2+qs[3]**2
    return np.array([[2*qs[0],-2*qs[1],-2*qs[2],2*qs[3]],\
                     [2*qs[1],2*qs[1],-2*qs[3],-2*qs[2]],\
                     [2*qs[2],2*qs[3],2*qs[0],2*qs[1]],\
                     [2*qs[3],-2*qs[2],2*qs[1],-2*qs[0]]])


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]/2,1),\
                                   (qf,m)).x
         exit()                          
         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
         dD=8*(Q1*dQ1+Q2*dQ2+Q3*dQ3+Q4*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

         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)


eartha=G*ms*(earthpos-sunpos)/norme**3
suna=-G*me*(sunpos-earthpos)/norma**3

xarray=[]
yarray=[]
zarray=[]
xarray.append(earthpos[0])
yarray.append(earthpos[1])
zarray.append(earthpos[2])

t=0
T=10**4
while t<T:
    k1v1=h*earthv
    k1v2=h*sunv
    k1a1=h*eartha
    k1a2=h*suna
    earthpos=earthpos+.5*k1v1
    sunpos=sunpos+.5*k1v2
    earthv=earthv+.5*k1a1
    sunv=sunv+.5*k1a2
    eartha=optimize.root(DotR, (1, 1, 1, 1),(earthv,earthpos,ms)).x[0:3]
    suna=optimize.root(DotR, (1, 1, 1, 1),(sunpos,sunv,-me)).x[0:3]

    k2v1=h*earthv
    k2v2=h*sunv
    k2a1=h*eartha
    k2a2=h*suna
    earthpos=earthpos+.5*k2v1
    sunpos=sunpos+.5*k2v2
    earthv=earthv+.5*k2a1
    sunv=sunv+.5*k2a2
    eartha=optimize.root(DotR, (1, 1, 1, 1),(earthv,earthpos,ms)).x[0:3]
    suna=optimize.root(DotR, (1, 1, 1, 1),(sunpos,sunv,-me)).x[0:3]

    k3v1=h*earthv
    k3v2=h*sunv
    k3a1=h*eartha
    k3a2=h*suna
    earthpos=earthpos+k3v1
    sunpos=sunpos+k3v2
    earthv=earthv+k3a1
    sunv=sunv+k3a2
    eartha=optimize.root(DotR, (1, 1, 1, 1),(earthv,earthpos,ms)).x[0:3]
    suna=optimize.root(DotR, (1, 1, 1, 1),(sunpos,sunv,-me)).x[0:3]

    k4v1=h*earthv
    k4v2=h*sunv
    k4a1=h*eartha
    k4a2=h*suna
    earthpos=earthpos+h/6.0*(k1v1+2*k2v1+2*k3v1+k4v1)
    sunpos=sunpos+h/6.0*(k1v2+2*k2v2+2*k3v2+k4v2)
    earthv=earthv+h/6.0*(k1a1+2*k2a1+2*k3a1+k4a1)
    sunv=sunv+h/6.0*(k1a2+2*k2a2+2*k3a2+k4a2)
    eartha=optimize.root(DotR, (1, 1, 1, 1),(earthv,earthpos,ms)).x[0:3]
    suna=optimize.root(DotR, (1, 1, 1, 1),(sunpos,sunv,-me)).x[0:3]
    xarray.append(earthpos[0])
    yarray.append(earthpos[1])

    t=t+h

plt.plot(xarray,yarray)
plt.savefig('orbit.png')