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([
    [1,xdot1],
    [x2,x1*x2]])
A1 = sp.Matrix([
    [1,x1],
    [x2,0]])
A2 = sp.Matrix([
    [1,0],
    [0,0]])
A0,A1,A2 # A(d/dt) = A0 + A1*d/dt + A2*(d/dt)**2


    

    
Out[3]:





$$\left ( \left[\begin{matrix}1 & \dot{x}_{1}\\x_{2} & x_{1} x_{2}\end{matrix}\right], \quad \left[\begin{matrix}1 & x_{1}\\x_{2} & 0\end{matrix}\right], \quad \left[\begin{matrix}1 & 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]:

    

# st.rnd_number_rank(spread(A0,A1,beta=2)) \
# == st.rnd_number_rank(spread_aug(A0,A1,beta=2)) \
# and st.rnd_number_rank(spread(A0,A1,beta=2)) == spread(A0,A1,beta=2).shape[0]


    

In [7]:

    

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


    

    




A is unimodular

In [8]:

    

beta_


    

    
Out[8]:





$$2$$

In [9]:

    

sA2 = Top(A0,A1,A2,beta=beta_)
sA2


    

    
Out[9]:





$$\left[\begin{matrix}1.0 & 1.0 \dot{x}_{1} & 1.0 & 1.0 x_{1} & 1.0 & 0 & 0 & 0 & 0 & 0\\1.0 x_{2} & 1.0 x_{1} x_{2} & 1.0 x_{2} & 0 & 0 & 0 & 0 & 0 & 0 & 0\\0 & 1.0 \ddot{x}_{1} & 1.0 & 2.0 \dot{x}_{1} & 1.0 & 1.0 x_{1} & 1.0 & 0 & 0 & 0\\1.0 \dot{x}_{2} & 1.0 x_{1} \dot{x}_{2} + 1.0 x_{2} \dot{x}_{1} & 1.0 x_{2} + 1.0 \dot{x}_{2} & 1.0 x_{1} x_{2} & 1.0 x_{2} & 0 & 0 & 0 & 0 & 0\\0 & 1.0 \dddot{x}_{1} & 0 & 3.0 \ddot{x}_{1} & 1.0 & 3.0 \dot{x}_{1} & 1.0 & 1.0 x_{1} & 1.0 & 0\\1.0 \ddot{x}_{2} & 1.0 x_{1} \ddot{x}_{2} + 1.0 x_{2} \ddot{x}_{1} + 2.0 \dot{x}_{1} \dot{x}_{2} & 1.0 \ddot{x}_{2} + 2.0 \dot{x}_{2} & 2.0 x_{1} \dot{x}_{2} + 2.0 x_{2} \dot{x}_{1} & 1.0 x_{2} + 2.0 \dot{x}_{2} & 1.0 x_{1} x_{2} & 1.0 x_{2} & 0 & 0 & 0\end{matrix}\right]$$

In [10]:

    

sA2_rpinv = right_pseudo_inverse(sA2)
sA2_rpinv.shape


    

    
Out[10]:





$$\left ( 10, \quad 6\right )$$

In [11]:

    

I, O = sp.eye(2), sp.zeros(2,8)
IO = st.col_stack(I,O)
B_ = IO*sA2_rpinv
B_ = sp.simplify(B_)


    

In [12]:

    

B_


    

    
Out[12]:





$$\left[\begin{matrix}1.0 & \frac{1.0 \dot{x}_{2}}{x_{2}^{2}} & 0 & - \frac{1.0}{x_{2}} & 0 & 0\\- \frac{1.0}{x_{1}} & \frac{1}{x_{1} x_{2}^{3}} \left(1.0 x_{2}^{2} - 1.0 x_{2} \ddot{x}_{2} - 1.0 x_{2} \dot{x}_{2} + 2.0 \dot{x}_{2}^{2}\right) & - \frac{1.0}{x_{1}} & \frac{1}{x_{1} x_{2}^{2}} \left(1.0 x_{2} - 2.0 \dot{x}_{2}\right) & 0 & \frac{1.0}{x_{1} x_{2}}\end{matrix}\right]$$

In [13]:

    

B_.shape


    

    
Out[13]:





$$\left ( 2, \quad 6\right )$$

In [14]:

    

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


    

    
Out[14]:





$$\left ( \left[\begin{matrix}1.0 & \frac{1.0 \dot{x}_{2}}{x_{2}^{2}}\\- \frac{1.0}{x_{1}} & \frac{1}{x_{1} x_{2}^{3}} \left(1.0 x_{2}^{2} - 1.0 x_{2} \ddot{x}_{2} - 1.0 x_{2} \dot{x}_{2} + 2.0 \dot{x}_{2}^{2}\right)\end{matrix}\right], \quad \left[\begin{matrix}0 & - \frac{1.0}{x_{2}}\\- \frac{1.0}{x_{1}} & \frac{1}{x_{1} x_{2}^{2}} \left(1.0 x_{2} - 2.0 \dot{x}_{2}\right)\end{matrix}\right], \quad \left[\begin{matrix}0 & 0\\0 & \frac{1.0}{x_{1} x_{2}}\end{matrix}\right]\right )$$