小编kil*_*roc的帖子

我目前正在致力于实施 Adams Bashforth Moulton 方法来解决钟摆问题。

我当前的代码如下：

import numpy as np

#####################################################################################
# DEF func
#####################################################################################
def func(y,t):

    n=y.size
    f=np.zeros(6)

    xp=np.array([y[3],y[4],y[5]])[np.newaxis]

    mass=1.
    F1=0.
    F2=0.
    F3=-mass*9.8
    F=np.array([F1,F2,F3])[np.newaxis]

    phix=2.*y[0]
    phiy=2.*y[1]
    phiz=2.*y[2]
    G=np.array([phix,phiy,phiz])[np.newaxis]

    H=2.*np.eye(3)
    lambd=(mass*np.dot(xp,np.dot(H,xp.T))+np.dot(F,G.T))/np.dot(G,G.T)

    f[0]=y[3]
    f[1]=y[4]
    f[2]=y[5]
    for k in range(0,3):
        f[k+3]=(F[0,k]-lambd*G[0,k])/mass

    return f

def dF_matrix(y):

    n=y.size
    dF=np.zeros((6,6))

    xp=np.array([y[3],y[4],y[5]])[np.newaxis]

    mass=1.
    F1=0.
    F2=0.
    F3=-mass*9.8
    F=np.array([F1,F2,F3])[np.newaxis]

    phix=2.*y[0]
    phiy=2.*y[1]
    phiz=2.*y[2]
    G=np.array([phix,phiy,phiz])[np.newaxis]

    H=2.*np.eye(3)
    lambd=(mass*np.dot(xp,np.dot(H,xp.T))+np.dot(F,G.T))/np.dot(G,G.T)

    dF[0,3]=1
    dF[1,4]=1
    dF[2,5]=1
    dF[3,0]=(y[0]*F1+2*lambd)/mass
    dF[3,1]=(y[0]*F2)/mass
    dF[3,2]=(y[0]*F3)/mass
    dF[3,3]=phix*y[3]
    dF[3,4]=phix*y[4]
    dF[3,5]=phix*y[5]
    dF[4,0]=(y[1]*F1)/mass
    dF[4,1]=(y[1]*F2+2*lambd)/mass
    dF[4,2]=(y[1]*F3)/mass
    dF[4,3]=phiy*y[3]
    dF[4,4]=phiy*y[4]
    dF[4,5]=phiy*y[5]
    dF[5,0]=(y[2]*F1)/mass
    dF[5,1]=(y[2]*F2)/mass
    dF[5,2]=(y[2]*F3+2*lambd)/mass
    dF[5,3]=phiz*y[3]
    dF[5,4]=phiz*y[4]
    dF[5,5]=phiz*y[5]

    return dF …