Lorenz.py 文件源码

def Jacobian(ssp, dt, nstp):
    """
    Jacobian function for the trajectory started on ssp, evolved for time t

    Inputs:
    ssp: Initial state space point. dx1 NumPy array: ssp = [x, y, z]
    t: Integration time
    Outputs:
    J: Jacobian of trajectory f^t(ssp). dxd NumPy array
    """
    #CONSTRUCT THIS FUNCTION
    #Hint: See the Jacobian calculation in CycleStability.py
    #J = None
    Jacobian0 = np.identity(3)  # COMPLETE THIS LINE. HINT: Use np.identity(DIMENSION)
    #Initial condition for Jacobian integral is a d+d^2 dimensional matrix
    #formed by concatenation of initial condition for state space and the
    #Jacobian:
    sspJacobian0 = np.zeros(3 + 3 ** 2)  # Initiate
    sspJacobian0[0:3] = ssp  # First 3 elemenets
    sspJacobian0[3:] = np.reshape(Jacobian0, 9)  # Remaining 9 elements
    tInitial = 0  # Initial time
    tFinal = dt*nstp  # Final time
    Nt = nstp  # Number of time points to be used in the integration

    tArray = np.linspace(tInitial, tFinal, Nt)  # Time array for solution

    sspJacobianSolution = odeint(JacobianVelocity, sspJacobian0, tArray)

    xt = sspJacobianSolution[:, 0]  # Read x(t)
    yt = sspJacobianSolution[:, 1]  # Read y(t)
    zt = sspJacobianSolution[:, 2]  # Read z(t)

    #Read the Jacobian for the periodic orbit:
    J = sspJacobianSolution[-1, 3:].reshape((3, 3))

    return J