Differential correction

This method is based on small (\(\approx 10^{-5}\)) corrections to the initial condition of an unstable periodic orbit at a desired total energy. The procedure is started from the linear solutions of the Hamilton’s equations at the index-1 saddle equilibrium point and which generates a small amplitude (\(\approx 10^{-5}\)) periodic orbit. This is then fed into an iterative procedure based on state transition matrix that calculates correction to the initial condition based on error in the terminal coordinates of the periodic orbit. This leads to convergence within 3 steps in the sense of the trajectory returning to the initial condition. Once a small amplitude periodic orbit is obtained, numerical continuation increases the amplitude and correspondingly the total energy of the unstable periodic orbit. Then, a combination of bracketing and bisection method computes the UPO at the desired energy for a specified tolerance.

differential_correction.clean_up_matrix(A)

Returns any complex matrix A with entries where the real or complex part has absolute value smaller than TOL set to 0, where TOL is set inside this function.

differential_correction.eig_get(A, discrete)

Returns the eigenvalues and eigenvectors of the matrix A spanning the three local subspaces <Es,Eu,Ec> where A is MxM and s+u+c=M, which locally approximate the invariant manifolds <Ws,Wu,Wc>

This function is designed for continuous dynamical systems only.

Parameters:

• A (2d numpy array) – Input matrix
• discrete (int) – flag to set if the dynamical system is discrete or continuous

Returns:

sn,un,cn,Ws,Wu,Wc – eigenvalues and corresponding eigenvectors

Return type:

1d numpy array

differential_correction.get_eq_pts(eqNum, init_guess_eqpt_model, grad_pot_model, par)

Returns configuration space coordinates of the equilibrium points.

get_eq_pts(eqNum, init_guess_eqpt_model, grad_pot_model, par) solves the coordinates of the equilibrium point for a Hamiltonian of the form kinetic energy (KE) + potential energy (PE).

Parameters:

• eqNum (int) – 1 is the saddle-center equilibrium point 2 or 3 is the center-center equilibrium point
• init_guess_eqpt_model (function name) – function that returns the initial guess for the equilibrium point
• grad_pot_model (function name) – function that defines the vector of potential gradient
• par (float (list)) – model parameters

Returns:

configuration space coordinates

Return type:

float (list)

differential_correction.get_po_diffcorr(x0, diffcorr_setup_model, conv_coord_model, diffcorr_acc_corr_model, ham2dof_model, half_period_model, pot_energy_model, variational_eqns_model, plot_iter_orbit_model, par)

Returns an initial condition and time period of an unstable periodic orbit using differential correction

The initial condition is generated using a small (first-order derived using state transition matrix) correction of the initial guess where the assumption is the unstable periodic orbit projects as a function nameaight line on the configuration space.

Parameters:

• x0 (float) – guess of the initial condition for the unstable periodic orbit
• diffcorr_setup_model (function name) – function that returns the combination of coordinates for applying terminal and periodic orbit conditions
• conv_coord_model (function name) – function that returns the coordinate for convergence criteria
• diffcorr_acc_corr_model (function name) – function that returns the corrected phase space coordinate and where the correction term is derived
• ham2dof_model (function name) – function that returns the Hamiltonian vector field at an input phase space coordinate and time
• half_period_model (function name) – function that returns the event criteria in terms of the coordinate that is set to zero for half-period of the unstable periodic orbit
• pot_energy_model (function name) – function that returns the potential energy of Hamiltonian
• variational_eqns_model (function name) – function that returns the variational equations of the dynamical system
• plot_iter_orbit_model (function name) – function to plot the computed orbit in the 3D phase space of 2 position and 1 momentum coordinate
• par (float (list)) – model parameters

Returns:

• x0po (1d numpy array) – initial condition for the periodic orbit
• t1 (float) – time period for the periodic orbit

differential_correction.get_po_fam(eqNum, Ax1, Ax2, nFam, po_fam_file, init_guess_eqpt_model, grad_pot_model, jacobian_model, guess_lin_model, diffcorr_setup_model, conv_coord_model, diffcorr_acc_corr_model, ham2dof_model, half_period_model, pot_energy_model, variational_eqns_model, plot_iter_orbit_model, par)

Returns a family of unstable periodic orbits given two seed initial conditions and time periods

Parameters:

• eqNum (int) – 1 is the saddle-center equilibrium point 2 or 3 is the center-center equilibrium point
• Ax2 (Ax1,) – amplitude of the guess periodic orbit (typically 1e-4 or 1e-5)
• nFam (int) – number of members in the family (monotonically increasing energy) of the unstable periodic orbits
• po_fam_file (function name) – file name to save the members in the family of the unstable periodic orbits
• init_guess_eqpt_model (function name) – function that returns the initial guess for the equilibrium point
• grad_pot_model (function name) – function that returns gradient of the potential energy function.
• jacobian_model (function name) – function that returns Jacobian of the vector field.
• guess_lin_model (function name) – function that returns the initial guess based on linearization around the equilibrium point.
• diffcorr_setup_model (function name) – function that returns the combination of coordinates for applying terminal and periodic orbit conditions
• conv_coord_model (function name) – function that returns the coordinate for convergence criteria
• diffcorr_acc_corr_model (function name) – function that returns the corrected phase space coordinate and where the correction term is derived
• ham2dof_model (function name) – function that returns the Hamiltonian vector field at an input phase space coordinate and time
• half_period_model (function name) – function that returns the event criteria in terms of the coordinate that is set to zero for half-period of the unstable periodic orbit
• pot_energy_model (function name) – function that returns the potential energy of Hamiltonian
• variational_eqns_model (function name) – function that returns the variational equations of the dynamical system
• plot_iter_orbit_model (function name) – function to plot the computed orbit in the 3D phase space of 2 position and 1 momentum coordinate

par : float (list)

model parameters

Returns:

• x0po (2d numpy array) – array of initial conditions for the members in the family of the unstable periodic orbits
• T (1d numpy array) – column vector of time period for the members in the family of the unstable periodic orbits

differential_correction.get_po_guess_linear(eqNum, Ax, init_guess_eqpt_model, grad_pot_model, jacobian_model, guess_lin_model, par)

Returns an initial guess for the differential correction method.

Uses a small displacement from the equilibrium point (in a direction on the collinear point’s center manifold) as a first guess for a planar periodic orbit (called a Lyapunov orbit). This initial condition and period are to be used as a first guess for the differential correction routine.

Parameters:

• eqNum (int) – index of the equilibrium point (typically 1 or 2 or 3)
• Ax (float) – small nondimensional amplitude of the periodic orbit (<< 1)
• init_guess_eqpt_model (function name) – function that returns an initial guess to solve the equilibrium point using numerical solvers.
• grad_pot_model (function name) – function that returns gradient of the potential energy function.
• jacobian_model (function name) – function that returns Jacobian of the vector field.
• guess_lin_model (function name) – function that returns the initial guess based on linearization around the equilibrium point.

Returns:

• x0poGuess (1d numpy array) – guess initial condition for the periodic orbit (first guess)
• TGuess (float) – gues time period for the periodic orbit (first guess)

differential_correction.get_pot_surf_proj(xVec, yVec, pot_energy_model, par)

Returns projection of the potential energy (PE) surface on the configuration space

Parameters:

• yVec (xVec,) – x,y-coordinates that discretizes the x, y domain of the configuration space
• pot_energy_model (function name) – function that returns the potential energy of Hamiltonian
• parameters (float (list)) – model parameters

Returns:

values of the PE

Return type:

2d numpy array

differential_correction.get_total_energy(orbit, pot_energy_model, parameters)

Returns total energy (value of Hamiltonian) of a phase space point on an orbit

get_total_energy(orbit, pot_energy_model, parameters) returns the total energy for a Hamiltonian of the form KE + PE.

Parameters:

• orbit (float (list)) – phase space coordinates (x,y,px,py) of a point on an orbit
• pot_energy_model (function name) – function that returns the potential energy of Hamiltonian
• parameters (float (list)) – model parameters

Returns:

total energy (value of Hamiltonian)

Return type:

scalar

differential_correction.po_bracket_energy(energyTarget, x0podata, po_brac_file, diffcorr_setup_model, conv_coord_model, diffcorr_acc_corr_model, ham2dof_model, half_period_model, pot_energy_model, variational_eqns_model, plot_iter_orbit_model, par)

Returns two unstable periodic orbits that bracket (bound in energy values) the unstable periodic orbit at the desired energy

Generates a family of periodic orbits (po) given a pair of seed initial conditions from a data file, while targeting a specific periodic orbit. This is performed using a scaling factor of the numerical continuation step size which is used to obtain initial guess for the periodic orbit. The energy of target periodic orbit should be higher than the input periodic orbit.

Parameters:

• energyTarget (float) – energy of the target unstable periodic orbit
• x0podata (2d numpy array) – array of input members of the family of the unstable periodic orbit
• po_brac_file (str) – file name to save the 2 bounds of the target unstable periodic orbit
• diffcorr_setup_model (function name) – function that returns the combination of coordinates for applying terminal and periodic orbit conditions
• conv_coord_model (function name) – function that returns the coordinate for convergence criteria
• diffcorr_acc_corr_model (function name) – function that returns the corrected phase space coordinate and where the correction term is derived
• ham2dof_model (function name) – function that returns the Hamiltonian vector field at an input phase space coordinate and time
• half_period_model (function name) – function that returns the event criteria in terms of the coordinate that is set to zero for half-period of the unstable periodic orbit
• pot_energy_model (function name) – function that returns the potential energy of Hamiltonian
• variational_eqns_model (function name) – function that returns the variational equations of the dynamical system
• plot_iter_orbit_model (function name) – function to plot the computed orbit in the 3D phase space of 2 position and 1 momentum coordinate
• par (float (list)) – model parameters

Returns:

• x0po (2d numpy array) – 2 bracketing unstable periodic orbits that bound the desired unstable periodic orbit
• T (1d numpy array) – time period of the 2 bracketing unstable periodic orbits

differential_correction.po_target_energy(x0po, energyTarget, po_target_file, diffcorr_setup_model, conv_coord_model, diffcorr_acc_corr_model, ham2dof_model, half_period_model, pot_energy_model, variational_eqns_model, plot_iter_orbit_model, par)

po_target_energy computes the periodic orbit of target energy using bisection method.

Using bisection method on the lower and higher energy values of the POs to find the PO with the target energy. Use this condition to integrate with event function of half-period defined by maximum distance from the initial point on the PO

Parameters:

• x0po (1d numpy array) – Initial conditions for the periodic orbit with the last two initial conditions bracketing (lower and higher than) the target energy
• diffcorr_setup_model (function name) – function that returns the combination of coordinates for applying terminal and periodic orbit conditions
• conv_coord_model (function name) – function that returns the coordinate for convergence criteria
• diffcorr_acc_corr_model (function name) – function that returns the corrected phase space coordinate and where the correction term is derived
• ham2dof_model (function name) – function that returns the Hamiltonian vector field at an input phase space coordinate and time
• half_period_model (function name) – function that returns the event criteria in terms of the coordinate that is set to zero for half-period of the unstable periodic orbit
• pot_energy_model (function name) – function that returns the potential energy of Hamiltonian
• variational_eqns_model (function name) – function that returns the variational equations of the dynamical system
• plot_iter_orbit_model (function name) – function to plot the computed orbit in the 3D phase space of 2 position and 1 momentum coordinate
• par (float (list)) – model parameters

Returns:

• x0_PO (1d numpy array) – Initial condition of the target unstable periodic orbit
• T_PO (float) – Time period of the target unstable periodic orbit
• ePO (float) – Energy of the target unstable periodic orbit.

differential_correction.remove_infinitesimals(A, TOL=1e-14)

Returns any complex matrix A with entries where the real or complex part has absolute value smaller than TOL set to 0

differential_correction.state_transit_matrix(tf, x0, par, variational_eqns_model, fixed_step=0)

Returns state transition matrix, the trajectory, and the solution of the variational equations over a length of time

In particular, for periodic solutions of % period tf=T, one can obtain the monodromy matrix, PHI(0,T).

Parameters:

• tf (float) – final time for the integration
• x0 (float) – initial condition
• par (float (list)) – model parameters
• variational_eqns_model (function name) – function that returns the variational equations of the dynamical system

Returns:

• t (1d numpy array) – solution time
• x (2d numpy array) – solution of the phase space coordinates
• phi_tf (2d numpy array) – state transition matrix at the final time, tf
• PHI (1d numpy array,) – solution of phase space coordinates and corresponding tangent space coordinates