Page MenuHomec4science
Files F93764298

osc1d.py
No OneTemporary
Actions

File Metadata

Created
Sun, Dec 1, 07:40

osc1d.py
View Options

import numpy as np
import matplotlib.pylab as plt
def exact(t):
# exact solution based on the initial value
x_ex = np.sin(t)
v_ex = np.cos(t)
return x_ex, v_ex
def euler(x,v,dt):
for i in range(x.size-1):
x[i+1] = x[i] + dt * v[i]
v[i+1] = v[i] - dt * x[i]
return x, v
def euler_pred_corr(x,v,dt):
for i in range(x.size-1):
x[i+1] = x[i] + dt * v[i]
v[i+1] = v[i] - dt * x[i]
x[i+1] = x[i] + dt * (v[i] + v[i+1])/2
v[i+1] = v[i] - dt * (x[i] + x[i+1])/2
return x, v
def euler2nd(x,v,dt):
#implement here 2nd order Euler integrator
for i in range(x.size-1):
x[i+1] = x[i] + dt * v[i] - dt**2 * x[i] / 2.
v[i+1] = v[i] - dt * x[i]
return x, v
def verlet(x,v,dt):
#implement here Verlet integrator (Leap-frog)
for i in range(x.size-1):
v[i+1] = v[i] - dt * x[i]
x[i+1] = x[i] + dt * v[i+1]
# transpose velocity
for i in range(v.size-1):
v[i] = (v[i] + v[i+1]) / 2.
return x, v
def verletv(x,v,dt):
#implement here velocity Verlet integrator (velocity Crank Nicholson)
for i in range(x.size-1):
x[i+1] = x[i] + dt * v[i] - 0.5 * dt**2 * x[i]
v[i+1] = v[i] - ( x[i] + x[i+1] ) * dt / 2.
return x, v
def fit(t,E):
poly = np.polyfit(t,E,1)
g_E = lambda x: poly[0]*x + poly[1]
return g_E(t),poly[0] #return fitted E and slope
def plotxv(x_ex, v_ex, x, v, t):
font = {'fontname':'sans-serif'}
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(t, x_ex, 'r-', t, x, 'r.', t, v_ex, 'g-', t, v, 'g.')
ax.set_xlabel('time', **font)
ax.autoscale(axis='x', tight=True)
plt.show()
def plotE(E, t):
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(t, E, 'b.')
ax.set_xlabel('time')
ax.set_ylabel('energy')
ax.autoscale(axis='x', tight=True)
plt.show()

Event Timeline

c4science · Help