
import numpy as np
import matplotlib.pyplot as plt

f = 1.
w = 2*np.pi*f

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
	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')

Oscilador Harmônico