# eigenvector example
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
a = np.array([[-1,-10,-10],[1,0,0],[0,1,0]],dtype=float)
b = np.array([[0],[0],[0]],dtype=float)
c = np.eye(3)
d = np.zeros((3,1))
sys = signal.StateSpace(a, b, c, d)
t, yout, xout = signal.lsim(sys, U=0, T=np.linspace(0,10,100), X0=[0,-1,1])
plt.figure(figsize=(8,4), dpi=100)
plt.plot(t,xout)
plt.grid()
plt.show()
# (right) eigenvectors
D,V = np.linalg.eig(a)
print(f'eval_r:\n{D}')
print(f'evec_r:\n{V}')
# left eigenvectors
D,W = np.linalg.eig(a.T)
print(f'eval_l:\n{D}')
print(f'evec_l:\n{W}')
eval_r: [-8.32667268e-16+3.16227766j -8.32667268e-16-3.16227766j -1.00000000e+00+0.j ] evec_r: [[ 9.49157996e-01+0.00000000e+00j 9.49157996e-01-0.00000000e+00j 5.77350269e-01+0.00000000e+00j] [-1.72167284e-17-3.00150113e-01j -1.72167284e-17+3.00150113e-01j -5.77350269e-01+0.00000000e+00j] [-9.49157996e-02-2.32611182e-17j -9.49157996e-02+2.32611182e-17j 5.77350269e-01+0.00000000e+00j]] eval_l: [-3.33066907e-16+3.16227766j -3.33066907e-16-3.16227766j -1.00000000e+00+0.j ] evec_l: [[-6.42824347e-02+0.20327891j -6.42824347e-02-0.20327891j 9.95037190e-02+0.j ] [-7.07106781e-01+0.j -7.07106781e-01-0.j 1.08016402e-16+0.j ] [-6.42824347e-01-0.20327891j -6.42824347e-01+0.20327891j 9.95037190e-01+0.j ]]
# w3: left eigenvector associated with lambda=-1
w3 = W[:,2]
print(f'w3:\n{w3}')
w3: [9.95037190e-02+0.j 1.08016402e-16+0.j 9.95037190e-01+0.j]
# (w_3^T)x history
plt.figure(figsize=(8,4), dpi=100)
plt.plot(t,xout.dot(w3))
plt.grid()
plt.show()
/usr/local/lib/python3.10/dist-packages/matplotlib/cbook/__init__.py:1335: ComplexWarning: Casting complex values to real discards the imaginary part return np.asarray(x, float)
# time history for x_0 = V[:,0].real
t, yout, xout = signal.lsim(sys, U=0, T=np.linspace(0,10,100), X0=V[:,0].real)
plt.figure(figsize=(8,4), dpi=100)
plt.plot(t,xout)
plt.grid()
plt.show()