solve_osc1d.py
e760d1e4250bmaster

solve_osc1d.py
View Options

			#!/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))