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