STMint package

STMint module

class STMint.STMint.STMint(vars=None, dynamics=None, preset='', preset_mult=1, variational_order=1)[source]

Bases: object

State Transition Matrix Integrator

A tool for numerical integration of variational equations associated with a symbolically specified dynamical system.

Constructor Parameters:

vars (1-dimensional sympy matrix):

The variables used in the symbolic integration.

dynamics (sympy expression(s)):

The dynamics to be symbolically integrated.

preset (string):

Dynamic and Variational equation preset. Current presets are:

“twoBody”:: Two body motion.
“twoBodyEarth”:: Two body motion around Earth.
“twoBodySun”:: Two body motion around the Sun.
“threeBody”:: Three body motion.
“threeBodySunEarth”:: Three body motion around the Sun and Earth.
“threeBodyEarthMoon”:: Three body motion around the Earth and Moon.

preset_mult (float):

Constant multiple of potential V for 2-body motion.

Note: Only needed if preset = “”.

variational_order (int):

Order of variational equations to be computed:

0 - for no variational equations.
1 - for first order variational equations.
2 - for first and second order variational equations.

vars

The variables used in the symbolic integration.

Type:: 1-d sympy matrix

dynamics

The dynamics to be symbolically integrated.

Type:: sympy expression

lambda_dynamics

The lambdified dynamic equations.

Type:: lambdafied sympy expression

lambda_dynamics_and_variational

The lambdified dynamic and variational equations.

Type:: lambdafied sympy expression

dynVar_int(t_span, y0, output='raw', method='DOP853', t_eval=None, dense_output=False, events=None, vectorized=False, args=None, **options)[source]

Clone of scipy.solve_ivp

Method uses _dynamicsVariationalSolver to solve an initial value problem with given dynamics and variational equations. This method has the same optional arguments as Scipy’s solve_ivp function.

Non-optional arguments are listed below. See documentation of solve_ivp for a full list and description of arguments and returns https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html

Parameters:

t_span (2-tuple of floats) – Interval of integration (t0, tf). The solver starts with t=t0 and integrates until it reaches t=tf.
y0 (array_like, shape (n,)) – Initial state. For problems in the complex domain, pass y0 with a complex data type (even if the initial value is purely real).
output (str) –
Output of dynVar_int, options include:
- ”raw”:
  Raw bunch object from scipy.solve_ivp.
- ”final”:
  The state vector and STM at the final time only.
- ”all”:
  The state vector and STM at all times.

Returns:

If output is ‘raw’:

Bunch object with multiple defined fields, such as:

t (ndarray, shape (n_points,)):
Time points.

y (ndarray, shape (n, n_points)):
Values of the solution at t.

sol (OdeSolution or None):
Found solution as OdeSolution instance; None if dense_output was set to False.

If output is ‘final’:

state (n-array):: The state vector.
stm (ndarray):: The state transition matrix.

If output is ‘all’:

states (n-array):: The state vectors.
stms (ndarray):: The state transition matricies.
ts (n-array):: The time steps of integration.

Return type:

varies

dynVar_int2(t_span, y0, output='raw', method='DOP853', t_eval=None, dense_output=False, events=None, vectorized=False, args=None, **options)[source]

Clone of scipy.solve_ivp

Method uses _dynamicsVariationalSolver to solve an initial value problem with given dynamics and variational equations. This method has the same optional arguments as Scipy’s solve_ivp function. Note that this method also integrates second order variational equations to obtain a second order state transition tensor.

Non-optional arguments are listed below. See documentation of solve_ivp for a full list and description of arguments and returns https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html

Parameters:

t_span (2-tuple of floats) – Interval of integration (t0, tf). The solver starts with t=t0 and integrates until it reaches t=tf.
y0 (array_like, shape (n,)) – Initial state. For problems in the complex domain, pass y0 with a complex data type (even if the initial value is purely real).
output (str) –

Output of dynVar_int, options include:

”raw”:
Raw bunch object from solve_ivp.

”final”:
The state vector, STM, and STT at the final time only.

”all”:
The state vector, STM, and STT at all times.

Returns:

If output is ‘raw’:

Bunch object with multiple defined fields, such as:

t (ndarray, shape (n_points,)):: Time points.
y (ndarray, shape (n, n_points)):: Values of the solution at t.
sol (OdeSolution or None):: Found solution as OdeSolution instance; None if dense_output was set to False.

If output is ‘final’:

state (n-array):: The state vector.
stm (ndarray):: The state transition matrix.
stt (ndarray):: The state transition tensor.

If output is ‘all’:

states (n-array):: The state vectors.
stms (ndarray):: The state transition matricies.
stts (ndarray):: The state transition tensors.
ts (n-array):: The time steps of integration.

Return type:

varies

dyn_int(t_span, y0, method='DOP853', t_eval=None, dense_output=False, events=None, vectorized=False, args=None, **options)[source]

Clone of scipy.solve_ivp

Method uses _dynamicsSolver to solve an initial value problem with given dynamics. This method has the same arguments and Scipy’s solve_ivp function.

Non-optional arguments are listed below. See documentation of solve_ivp for a full list and description of arguments and returns https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html

Parameters:

t_span (2-tuple of floats) – Interval of integration (t0, tf). The solver starts with t=t0 and integrates until it reaches t=tf.
y0 (array_like, shape (n,)) – Initial state. For problems in the complex domain, pass y0 with a complex data type (even if the initial value is purely real).

Returns:

t (ndarray, shape (n_points,)):: Time points.
y (ndarray, shape (n, n_points)):: Values of the solution at t.
sol (OdeSolution or None):: Found solution as OdeSolution instance; None if dense_output was set to False.

Return type:

Bunch object with multiple defined fields

TensorNormUtilities module

STMint.TensorNormUtilities.MM_iterate(stringEin, tensOrder, tens, K, vec)[source]

Function to perform one step of polynomial optimization on a scalar valued polynomial on unit sphere

Single step

Parameters:

stringEin (string) – String to instruct einsum to perform contractions
tensOrder (int) – Order of the tensor
tens (np array) – Tensor
K (double) – Damping constant
vec (np array) – Vector

Returns:

vecNew (np array)

vecNorm (float)

STMint.TensorNormUtilities.MM_iteration(tens, vecGuess, maxIter, tol)[source]

Function to perform polynomial optimization on a scalar valued polynomial on unit sphere

Parameters:

tens (np array) – Tensor
vec (np array) – Vector
maxIter (int) – Max number of iterations to perform
tol (float) – Tolerance for difference and iterates

Returns:

eigVec (np array)

eigValue (np array)

STMint.TensorNormUtilities.cocycle1(stm10, stm21)[source]

Function to find STM along two combined subintervals

The cocycle conditon equation is used to find Phi(t2,t_0)=Phi(t2,t_1)*Phi(t1,t_0)

Args:

stm10 (np array)
State transition matrix from time 0 to 1

stm21 (np array)
State transition matrix from time 1 to 2

Returns:

stm20 (np array)
State transition matrix from time 0 to 2

STMint.TensorNormUtilities.cocycle2(stm10, stt10, stm21, stt21)[source]

Function to find STM and STT along two combined subintervals

The cocycle conditon equation is used to find Phi(t2,t0)=Phi(t2,t1)*Phi(t1,t0)

and the generalized cocycle condition is used to find Psi(t2,t0)

Args:

stm10 (np array): State transition matrix from time 0 to 1
stt10 (np array): State transition tensor from time 0 to 1
stm21 (np array): State transition matrix from time 1 to 2
stt21 (np array): State transition tensor from time 1 to 2

Returns:

stm20 (np array): State transition matrix from time 0 to 2
stt20 (np array): State transition tensor from time 0 to 2

STMint.TensorNormUtilities.get_polynomial_bound(tens)[source]

Function to find a bound on the value of a sclar valued polynomial on the unit sphere

Parameters:

tens (np array) – Tensor

Returns:

K (double): Bound on the polynomial on the unit sphere

STMint.TensorNormUtilities.nonlin_index_2(stm, stt)[source]

Function to calculate the nonlinearity index

Using tensor eigenvalues, the quotient of the induced 2-norm of the STT with the 2-norm of the STM

Parameters:

stm (np array) – State transition matrix (used to generate guess)
stt (np array) – Arbitrary order state transition tensor

Returns:

nonlinearity_index (float)

STMint.TensorNormUtilities.nonlin_index_DEMoN2(stm, stt)[source]

Function to calculate the nonlinearity index

Using tensor eigenvalues, the quotient of the induced 2-norm of the STT with the 2-norm of the STM

Parameters:

stm (np array) – State transition matrix (used to generate guess)
stt (np array) – Arbitrary order state transition tensor

Returns:

nonlinearity_index (float)

STMint.TensorNormUtilities.nonlin_index_TEMoN3(stm, stt)[source]

Function to calculate the nonlinearity index

Using tensor eigenvalues, the quotient of the induced 2-norm of the 3rd order term in the CGT series with the 2-norm of the second order term in the CGT series

Parameters:

stm (np array) – State transition matrix (used to generate guess)
stt (np array) – Arbitrary order state transition tensor

Returns:

nonlinearity_index (float)

STMint.TensorNormUtilities.nonlin_index_inf_2(stm, stt)[source]

Function to calculate the nonlinearity index

The induced infinity-2 norm is used in this calculation

Args:

stm (np array)
State transition matrix

stt (np array)
Second order state transition tensor

Returns:
nonlinearity_index (float)

STMint.TensorNormUtilities.nonlin_index_junkins_scale_free(stm, stt)[source]

Function to calculate the nonlinearity index

The induced 2 norm of the unfolded STT is used in this calculation This gives the quotient of the (Frobenius, 2)-norm of the second order STT with the Frobenius norm of the STM.

Args:

stm (np array)
State transition matrix

stt (np array)
Second order state transition tensor

Returns:
nonlinearity_index (float)

STMint.TensorNormUtilities.nonlin_index_unfold_bound(stm, stt)[source]

Function to calculate the nonlinearity index

The induced 2 norm of the unfolded STT is used in this calculation. This is a bound on the 2-norm of the second order STT quotiented with the 2-norm of the STM

Args:

stm (np array)
State transition matrix

stt (np array)
Second order state transition tensor

Returns:
nonlinearity_index (float)

STMint.TensorNormUtilities.power_iterate(stringEin, tensOrder, tens, vec)[source]

Function to perform one higher order power iteration on a symmetric tensor

Single step

Parameters:

stringEin (string) – String to instruct einsum to perform contractions
tensOrder (int) – Order of the tensor
tens (np array) – Tensor
vec (np array) – Vector

Returns:

vecNew (np array)

vecNorm (float)

STMint.TensorNormUtilities.power_iterate_string(tens)[source]

Function to calculate the index string for einsum (up to 26 dimensional tensor)

Parameters:: tens (np array) – Tensor
Returns:: einsum string to perform power iteration (string)

STMint.TensorNormUtilities.power_iterate_symmetrizing(stringEin, tensOrder, tens, vec)[source]

Function to perform one higher order power iteration on a non-symmetric tensor

Parameters:

stringEin (string) – String to instruct einsum to perform contractions
tensOrder (int) – Order of the tensor
tens (np array) – Tensor
vec (np array) – Vector

Returns:

vecNew (np array)

vecNorm (float)

STMint.TensorNormUtilities.power_iteration(tens, vecGuess, maxIter, tol)[source]

Function to perform higher order power iteration on a symmetric tensor

Parameters:

tens (np array) – Tensor
vec (np array) – Vector
maxIter (int) – Max number of iterations to perform
tol (float) – Tolerance for difference and iterates

Returns:

eigVec (np array)

eigValue (np array)

STMint.TensorNormUtilities.power_iteration_symmetrizing(tens, vecGuess, maxIter, tol)[source]

Function to perform higher order power iteration on a non-symmetric tensor

Parameters:

tens (np array) – Tensor
vec (np array) – Vector
maxIter (int) – Max number of iterations to perform
tol (float) – Tolerance for difference and iterates

Returns:

eigVec (np array)

eigValue (np array)

STMint.TensorNormUtilities.stt_2_norm(stm, stt)[source]

Function to calculate the norm of the state transition tensor, and the input unit vector that leads to that norm.

The maximum eigenvalue of the tensor squared computed with symmetrization along the way

Parameters:

stm (np array) – State transition matrix
stt (np array) – Second order state transition tensor

Returns:

sttArgMax (np array): Input unit vector that maximizes the STT
sqrt(sttNorm) (float): Square root of the norm of the STT

STMint.TensorNormUtilities.symmetrize_tensor(tens)[source]

Symmetrize a tensor

Parameters:: tens (np array) – Tensor
Returns:: symTens (np array)

STMint.TensorNormUtilities.tensor_2_norm(tens, guessVec)[source]

Function to calculate the norm of a state transition tensor

The square root of the maximum eigenvalue of the tensor squared

Parameters:

tens (np array) – Arbitrary 1-m tensor
guessVec (np array) – Guess vector for input that maximizes the tensor

Returns:

nonlinearity_index (float)

STMint.TensorNormUtilities.tensor_square_string(tens)[source]

Function to calculate the index string for einsum (up to 1-13 dimensional tensor) :param tens: Tensor :type tens: np array

Returns:: einsum string to perform tensor squaring (string)