In [2]:
from numpy import *

In [3]:
import matplotlib.pyplot as plt
%config InlineBackend.figure_format = 'retina'
%matplotlib inline
#%matplotlib notebook

In [4]:
def f(X):
x,y = X
return array([-y,x])

def eulerstep(f,X,t,h):
return X + f(X)*h

X = array([1.,0.])
T = 2*pi  # final time we want to get to
n = 50  # number of steps
h = T/n # time-step ("Delta-t")

#print(0,X)
plt.subplot(111,aspect=1)

# draw exact solution
t = linspace(0,2*pi,300)
plt.plot(cos(t),sin(t),'k',alpha=0.5)
for i in range(n):
newX =  eulerstep(f,X,t,h)  # Euler step
plt.plot([X[0],newX[0]],[X[1],newX[1]],'bo-',ms=2,alpha=0.3)
X = newX
t = (i+1)*h
#print(t,X)
plt.show()

In [5]:
def f(X,t):
x,y = X
return array([-y,x])

def eulerstep(f,X,t,h):
return X + f(X,t)*h

def improved_euler_step(f,X,t,h):
f1 = f(X,t)
X1 =  X + f1*h  # provisional Euler step
f2 = f(X1,t+h)
fav = (f1+f2)/2
return X + h*fav

def rk4(f,X,t,h):
# See Runge-Kutta in Wikipedia
pass

X = array([1.,0.])
T = 2*pi  # final time we want to get to
n = 50  # number of steps
h = T/n # time-step ("Delta-t")

#print(0,X)
plt.subplot(111,aspect=1)

# draw exact solution
t = linspace(0,2*pi,300)
plt.plot(cos(t),sin(t),'k',alpha=0.5)
for i in range(n):
newX =  improved_euler_step(f,X,t,h)  # Euler step
plt.plot([X[0],newX[0]],[X[1],newX[1]],'bo-',ms=2,alpha=0.3)
X = newX
t = (i+1)*h
#print(t,X)
plt.show()


# Newton's cannon¶

In [9]:
def f(X,t):
G = 1
M = 1
x,y,u,v = X
r = sqrt(x*x + y*y)
r32 = r**3
return array([u,v,-G*M*x/r32,-G*M*y/r32])

def eulerstep(f,X,t,h):
return X + f(X,t)*h

def improved_euler_step(f,X,t,h):
f1 = f(X,t)
X1 =  X + f1*h  # provisional Euler step
f2 = f(X1,t+h)
fav = (f1+f2)/2
return X + h*fav

def rk4(f,X,t,h):
# See Runge-Kutta in Wikipedia
pass

T = 30  # final time we want to get to
#n = 50  # number of steps
#h = T/n # time-step ("Delta-t")
h = .05
n = int(T/h)
#print(0,X)
plt.figure(figsize=(12,12))
plt.subplot(111,aspect=1)

# draw exact solution
t = linspace(0,2*pi,300)
plt.plot(cos(t),sin(t),'k',alpha=0.5)

for u0 in [1.,1.1,1.2]:
X = array([0,1.2,u0,0])
for i in range(n):
newX =  improved_euler_step(f,X,t,h)  # Euler step
plt.plot([X[0],newX[0]],[X[1],newX[1]],'b-')#,ms=2,alpha=0.3)
X = newX
t = (i+1)*h
#print(t,X)
plt.show()


# Animation¶

In [15]:
M1 = 0.8
M2 = 1 - M1
R1 = M2
R2 = M1
n = 100
t = linspace(0,2*pi,n,endpoint=False)
x1 = R1*cos(t)
y1 = R1*sin(t)
x2 =-R2*cos(t)
y2 =-R2*sin(t)
fig = plt.figure(figsize=(4,4))
ax = plt.subplot(111,aspect=1)
ax.plot(x1,y1,'b')
ax.plot(x2,y2,'r');

In [16]:
from matplotlib import animation
from IPython.display import HTML

In [17]:
def init():
#star1.set_data = []
#star2.set_data = []
return star1,star2

def animate(i): # plot ith frame
star1.set_data(x1[i],y1[i])
star2.set_data(x2[i],y2[i])
return star1,star2

In [23]:
fig = plt.figure(figsize=(4,4))
ax = plt.subplot(111,aspect=1)
#ax.plot(0,0,'ko')
#star1, = ax.plot(x1[0],y1[0],'bo',ms=6)
#star2, = ax.plot(x2[0],y2[0],'ro',ms=6)
star1, = ax.plot([],[],'bo',ms=10*sqrt(M1))
star2, = ax.plot([],[],'ro',ms=10*sqrt(M2))
ax.set_xlim(-2,2)
ax.set_ylim(-2,2)
#,-2,2)
anim = animation.FuncAnimation(fig,animate,init_func=init,frames=n,interval=5,blit=True)
HTML(anim.to_jshtml())

Out[23]:

Once Loop Reflect