In [1]:
from hypore import *

In [2]:
# create state vector with 4 components and derivatives up to order 2
# and make available in global namespace
xx = gen_state(2,3)
xx.T


Out[2]:
$$\left[\begin{matrix}x_{1} & x_{2} & \dot{x}_{1} & \dot{x}_{2} & \ddot{x}_{1} & \ddot{x}_{2} & \dddot{x}_{1} & \dddot{x}_{2}\end{matrix}\right]$$

In [3]:
A0 = sp.Matrix([
    [0,0,0],
    [xdot2,0,-1],
    [0,1,0]])
A1 = sp.Matrix([
    [-xdot2,-xdot1,1],
    [0,0,0],
    [0,0,0]])
A0,A1 # A(d/dt) = A0 + A1*d/dt


Out[3]:
$$\left ( \left[\begin{matrix}0 & 0 & 0\\\dot{x}_{2} & 0 & -1\\0 & 1 & 0\end{matrix}\right], \quad \left[\begin{matrix}- \dot{x}_{2} & - \dot{x}_{1} & 1\\0 & 0 & 0\\0 & 0 & 0\end{matrix}\right]\right )$$

In [4]:
# manually:
# st.rnd_number_rank(spread(A0,A1,beta=0)) \
# == st.rnd_number_rank(spread_aug(A0,A1,beta=0)) \
# and st.rnd_number_rank(spread(A0,A1,beta=0)) == spread(A0,A1,beta=0).shape[0]

In [5]:
# st.rnd_number_rank(spread(A0,A1,beta=1)) \
# == st.rnd_number_rank(spread_aug(A0,A1,beta=1)) \
# and st.rnd_number_rank(spread(A0,A1,beta=1)) == spread(A0,A1,beta=1).shape[0]

In [6]:
# automated:
# returns degree beta of inverse (-1 => not unimodular)
beta_ = is_unimodular(A0,A1)
if (beta_==-1):
    print("A is not unimodular")
else:
    print("A is unimodular")


A is unimodular

In [7]:
sA1 = Top(A0,A1,beta=beta_)
sA1


Out[7]:
$$\left[\begin{matrix}0 & 0 & 0 & - 1.0 \dot{x}_{2} & - 1.0 \dot{x}_{1} & 1.0 & 0 & 0 & 0\\1.0 \dot{x}_{2} & 0 & -1.0 & 0 & 0 & 0 & 0 & 0 & 0\\0 & 1.0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & - 1.0 \ddot{x}_{2} & - 1.0 \ddot{x}_{1} & 0 & - 1.0 \dot{x}_{2} & - 1.0 \dot{x}_{1} & 1.0\\1.0 \ddot{x}_{2} & 0 & 0 & 1.0 \dot{x}_{2} & 0 & -1.0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 1.0 & 0 & 0 & 0 & 0\end{matrix}\right]$$

In [8]:
sA1_rpinv = sp.simplify(right_pseudo_inverse(sA1))
sA1_rpinv


Out[8]:
$$\left[\begin{matrix}\frac{1.0}{\ddot{x}_{2}} & 0 & 0 & 0 & \frac{1.0}{\ddot{x}_{2}} & \frac{1.0 \dot{x}_{1}}{\ddot{x}_{2}}\\0 & 0 & 1.0 & 0 & 0 & 0\\\frac{1.0 \dot{x}_{2}}{\ddot{x}_{2}} & -1.0 & 0 & 0 & \frac{1.0 \dot{x}_{2}}{\ddot{x}_{2}} & \frac{1.0 \dot{x}_{1}}{\ddot{x}_{2}} \dot{x}_{2}\\0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 1.0\\1.0 & 0 & 0 & 0 & 0 & 1.0 \dot{x}_{1}\\0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 1.0 & 0 & 1.0 \ddot{x}_{1}\end{matrix}\right]$$

In [9]:
I, O = sp.eye(3), sp.zeros(3,6)
IO = st.col_stack(I,O)
B_ = IO*sA1_rpinv
B_


Out[9]:
$$\left[\begin{matrix}\frac{1.0}{\ddot{x}_{2}} & 0 & 0 & 0 & \frac{1.0}{\ddot{x}_{2}} & \frac{1.0 \dot{x}_{1}}{\ddot{x}_{2}}\\0 & 0 & 1.0 & 0 & 0 & 0\\\frac{1.0 \dot{x}_{2}}{\ddot{x}_{2}} & -1.0 & 0 & 0 & \frac{1.0 \dot{x}_{2}}{\ddot{x}_{2}} & \frac{1.0 \dot{x}_{1}}{\ddot{x}_{2}} \dot{x}_{2}\end{matrix}\right]$$

In [10]:
B0 = st.col_select(B_, 0,1,2)
B1 = st.col_select(B_, 3,4,5)
# unimodular inverse of A(d/dt)
B0, B1 # B(d/dt) = B0 + B1*d/dt


Out[10]:
$$\left ( \left[\begin{matrix}\frac{1.0}{\ddot{x}_{2}} & 0 & 0\\0 & 0 & 1.0\\\frac{1.0 \dot{x}_{2}}{\ddot{x}_{2}} & -1.0 & 0\end{matrix}\right], \quad \left[\begin{matrix}0 & \frac{1.0}{\ddot{x}_{2}} & \frac{1.0 \dot{x}_{1}}{\ddot{x}_{2}}\\0 & 0 & 0\\0 & \frac{1.0 \dot{x}_{2}}{\ddot{x}_{2}} & \frac{1.0 \dot{x}_{1}}{\ddot{x}_{2}} \dot{x}_{2}\end{matrix}\right]\right )$$