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]:
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]:
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")
In [7]:
sA1 = Top(A0,A1,beta=beta_)
sA1
Out[7]:
In [8]:
sA1_rpinv = sp.simplify(right_pseudo_inverse(sA1))
sA1_rpinv
Out[8]:
In [9]:
I, O = sp.eye(3), sp.zeros(3,6)
IO = st.col_stack(I,O)
B_ = IO*sA1_rpinv
B_
Out[9]:
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]: