solve_osc1d.py
No OneTemporary
Actions

Subscribers

None

File Metadata

	#!/usr/bin/env python

	import osc1d as solver
	import numpy as np

	# dt: step
	# x[0], v[0]: boundary conditions
	# j: number of cycles
	n = 100; j=2
	dt = 2np.pij/n #step size
	t = np.zeros(n)

	for i in range(n):
	t[i] = i*dt
	x = np.zeros(n); v = np.zeros(n)
	x[0] = 0; v[0] = 1

	X_ex, V_ex = solver.exact(t)

	#simple Euler integrator
	#X, V = solver.euler(x,v,dt)

	#uncomment to use predictor-corrector Euler
	#X, V = solver.euler_pred_corr(x,v,dt)

	#second order Euler integrator
	#X, V = solver.euler2nd(x,v,dt)

	#verlet integrator
	#X, V = solver.verlet(x,v,dt)

	#verlet velocity integrator
	X, V = solver.verletv(x,v,dt)

	# maximum deviation
	maxerr = max(np.abs(X - X_ex))
	print("#INFO: Maximum deviation: ", maxerr)

	E = x2 + v2

	# fit energy, get slope
	fitted, slope = solver.fit(t, E)
	print("#INFO: Slope of fitted energy: ", slope)

	# computation of noise
	n = len(E)
	noise = np.sqrt( np.sum((E - fitted)**2) / (n-2))
	print("#INFO: Noise in the energy: ", noise)

	print("#INFO: number of steps: %f" %n)
	#print("#INFO: number of cycles: %f" %j)
	print("#INFO: time step: %f" %dt)
	print("# time x v x-analytical v-analytic energy")
	for t,x,v,x_ex,v_ex,E in zip(t,X,V,X_ex,V_ex,E):
	print (" {:5.4f} {:5.4f} {:5.4f} {:5.4f} {:5.4f} {:5.4f}"\
	.format(t,x,v,x_ex,v_ex,E))