Page MenuHomec4science
Files F102417458

serie1.py
No OneTemporary
Actions

File Metadata

Created
Thu, Feb 20, 12:36

serie1.py
View Options

#!/usr/bin/python
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
# parameters
e = 0.6
h = 2e-4
N = 200
def gHq(q):
return q / np.linalg.norm(q)**3
def euler_step(p, q, h, Hp, Hq):
qold = np.copy(q)
q += h * Hp(p)
p -= h * Hq(qold)
def euler_symplectic(p, q, h, Hp, Hq):
q += h * Hp(p)
p -= h * Hq(q)
def euler_implicit(p, q, h, Hp, Hq):
f = lambda x: x - q - h * Hp(p - h * Hq(x))
q = scipy.optimize.fsolve(f)
p -= h * Hq(q)
def implicit_midpoint(p, q, h, Hp, Hq):
pnext = lambda x: p - h * Hq((q + x)/ 2.)
f = lambda x: x - (q + h * Hp( (p + pnext(x)) / 2.))
q = scipy.optimize.fsolve(f)
p = pnext(q)
def plot_energy(E, T):
for i in range(1,N+1):
E[i] = np.abs(E[i] - E[0])
E[0] = 0
fig, axis = plt.subplots()
axis.set_title("Energy conservation")
axis.set_xlabel("t")
axis.set_ylabel("|E - E_0|")
a,b = np.polyfit(np.log(T[1:-1]),np.log(E[1:-1]),1)
print("Energy error fit: %.5g * t^%.2f" % (np.exp(b), a))
axis.loglog(T[1:-1], E[1:-1], 'b*', label='Integrator')
axis.legend()
def plot_L(L, T):
for i in range(1,N+1):
L[i] = np.abs(L[i] - L[0])
L[0] = 0
fig, axis = plt.subplots()
axis.set_title("Angular momentum conservation")
axis.set_xlabel("t")
axis.set_ylabel("|L - L_0|")
a,b = np.polyfit(np.log(T[1:-1]),np.log(L[1:-1]),1)
print("L error fit: %.5g * t^%.2f" % (np.exp(b), a))
axis.loglog(T[1:-1], L[1:-1], 'b*', label='Integrator')
axis.legend()
# calls
def collect_energy(E, p, q):
E.append(np.dot(p, p) / 2.- 1 / np.linalg.norm(q))
def collect_L(L, p, q):
L.append(q[0] * p[1] - q[1] * p[0])
def main():
# initial conditions
q = np.array([1. - e, 0])
p = np.array([0, np.sqrt((1. + e) / (1. - e))])
# set integrator
#integrate = euler_step
integrate = euler_symplectic
# energy collection setup
E = []
collect_energy(E, p, q)
L = []
collect_L(L, p, q)
T = [h * n for n in range(N+1)]
# simulate
for n in range(N):
integrate(p, q, h, lambda x: x, gHq)
collect_energy(E, p, q)
collect_L(L, p, q)
# analysis and plots
plot_energy(E, T)
plot_L(L, T)
plt.show()
main()

Event Timeline

c4science · Help