import numpy as np
import matplotlib.pyplot as plt
f = 4.
w = 2*np.pi*f
xi = 2.
N = 10000
t = np.linspace(0,5,N)
h = t[1] - t[0]
def func(X):
x, vx = X
dx = vx
dvx = -w**2*x - xi*vx
return np.array([dx, dvx])
def RK4(X,h):
X = np.array(X)
K1 = h*func(X)
K2 = h*func(X+K1/2)
K3 = h*func(X+K2/2)
K4 = h*func(X+K3)
return X + (K1+2*K2+2*K3+K4)/6
x = np.zeros(N)
vx = np.zeros(N)
x[0] = 1.3
vx[0] = 0
for i in range(1,N):
X0 = x[i-1], vx[i-1]
x[i], vx[i] = RK4(X0,h)
fig = plt.figure()
plt.plot(t,x,'r-')
plt.grid(ls='--',alpha=0.5)
plt.xlabel('time [sec]')
plt.ylabel('x [meters]')
plt.show()
fig.savefig('fig.png')
