import numpy as np
import matplotlib.pyplot as plt
from numpy.linalg import inv

# Define the system parameters
m1 = 2
m2 = 3
k1 = 100
k2 = 200
c = 1

# Define the system matrices
A = np.array([[0, 0, 1, 0],
[(k2 - k1) / m2, -k2 / m1, 0, 0],
[0, 0, 0, 1],
[k2 / m2, -(k2 + k1) / m2, 0, -c / m2]])

B = np.array([[0], [0], [0], [1 / m2]])
C = np.array([[1, 0, 0, 0]])

# Define an initial state for simulation
x0 = np.zeros(shape=(4, 1))

# Define the number of time-samples used for the simulation
time = 30000
# Discretization constant
h = 0.001
# Define an input sequence for the simulation
input_seq = 10 * np.ones(shape=(1, time))
#plt.plot(input_sequence)

I = np.identity(A.shape[0]) #this is an identy matrix
A_D = inv(I - h * A)
B_D = h * np.matmul(A_D, B)
C_D = C

# Check the eigenvalues
# Continuous-time system
eigen_A = np.linalg.eig(A)[0]
# Discrete-time system
print(eigen_A)
eigen_A_D = np.linalg.eig(A_D)[0]

def simulate(Am, Bm, Cm, initial_state, input_sequence, time_steps):
Xm = np.zeros(shape=(Am.shape[0], time_steps + 1))
Ym = np.zeros(shape=(Cm.shape[0], time_steps + 1))

for i in range(0, time_steps):
if i == 0:
Xm[:, [i]] = initial_state
Ym[:, [i]] = np.matmul(Cm, initial_state)
input_tmp = input_sequence[0, i].reshape(Bm.shape[1], 1)
x = np.matmul(Am, initial_state) + np.matmul(Bm, input_tmp)
else:
Xm[:, [i]] = x
Ym[:, [i]] = np.matmul(Cm, x)
input_tmp = input_sequence[0, i].reshape(Bm.shape[1], 1)
x = np.matmul(Am, x) + np.matmul(Bm, input_tmp)

Xm[:, [-1]] = x
Ym[:, [-1]] = np.matmul(Cm, x)
return Xm, Ym

states, output = simulate(A_D, B_D, C_D, x0, input_seq, time)

plt.plot(states[1, :], color='blue', linewidth=2)
plt.xlabel('Discrete time instant-k')
plt.ylabel('State')
plt.title('System Response')
plt.savefig('step_response_Euler.png', dpi=600)
plt.show()






여기에서 시뮬을 돌리고 싶은데 혹시 어떻게 돌려야 하는지 알려주실 수 있을까요...?