Homework 1¶

due 09/12/2018 at 11:59pm

Problem 1¶

Compare the final positions of the particle as computed by the numerical integration of the ordinary differential equation using ODE_integration_E and compare it to the analytical (i.e. correct) integration given by Eq. $(2.5)$. To do this, plot the final distance obtained by ODE_integration_E along the y-axis as a function of the time step and show that it converges towards the value of the analytical solution

Problem 2¶

Compute the angle of deflection of an electron "flying by" a fixed ion as a function of the initial velocity of the electron. Infer the dependence of the angle of deflection with velocity. Comment on how the discretization impacts the angle at low velocities.

Problem 3¶

Plot the strength (color) and direction (arrows) of the electric field produced by one single charge in 2D. Try to do it in 3D after that. While it is not that difficult, 3D plot tend to be difficult to read if not done properly.

In [29]:
import math
from math import sqrt as sqrt
import numpy as np
import matplotlib.pyplot as plt
#to display inside the notebook!
%matplotlib inline


Solutions¶

Problem 1¶

We have found that $$y=\frac{q}{2m}E_y t^2+y_0=\frac{q}{2m}E_y t^2$$ and we want to compare the analytical result to the numerical result from ODE_integration_E given below.

In [35]:
def ODE_integration_E(x_initial,y_initial,vx_initial,vy_initial,t_initial,t_final,n_steps):
x=np.full(1,x_initial) # we do this to initialize the x variable as a numpy array
y=np.full(1,y_initial)
px_old=vx_initial*m
py_old=vy_initial*m
dt=(t_final-t_initial)/n_steps
vx=0
vy=0
for i in range(1,n_steps,1):
fx=q*Ex #this is the force along x
fy=q*Ey #this is the force along y
px_new=fx*dt+px_old #new momentum
py_new=fy*dt+py_old
vx=px_new/m #new velocity
vy=py_new/m
x_new=vx*dt+x[i-1] #new position
y_new=vy*dt+y[i-1] # in python arrays start at index 0
x=np.append(x,x_new) #the positions are appended at the end of the x and y arrays using numpy.
y=np.append(y,y_new)
px_old=px_new
py_old=py_new
return x,y


We now write the code to compare both results as a function of the time step dt.

In [57]:
n_max=600
err=np.zeros(n_max)
m=9.10938356e-31
q=-1.60217662e-19
Ex=0
Ey=1
t_f=1
vx_final=0
vy_final=0
y_real=q/(2*m)*Ey*t_f**2
for i in range(n_max):
x,y=ODE_integration_E(0,0,0,1000,0,t_f,i+1)
err[i]=math.log10(abs((y[i]-y_real)/y_real))
plt.plot(np.linspace(0,n_max,len(err)),err)
plt.show()


We find that the relative error converges towards $10^{-2.5}$ after 600 steps.

Problem 2¶

We define our functions. comp_E compute the electric field created at the location ro by a charge qp located at a point rp.

In [3]:
def comp_E(rp,qp,ro):
r=ro-rp
rn=np.linalg.norm(r)
return qp/(4*np.pi*8.85e-12)*r/rn**3


We the integrate the ODE using the integration function used in Ch.02. We used here a modified version of ODE_integration_B.

In [4]:
def ODE_integration_E(rp,initial_position,initial_velocity,time_frame,n_steps=200,q=-1.6e-19,m=9.1e-31,qp=1.6e-19):
loc=np.zeros((3,n_steps)) #3 for x,y,z
loc[:,0]=initial_position
v_old=initial_velocity
dt=(time_frame[1]-time_frame[0])/n_steps
f=np.zeros(3)
v_new=np.zeros(3)
for i in range(1,n_steps,1):
f=q*comp_E(rp,qp,loc[:,i-1])
v_new=f/m*dt+v_old #new momentum
loc[:,i]=v_new*dt+loc[:,i-1] #new position
v_old=v_new
return loc


We now build the for loop to compute the total deflection. We have place the stationary charge at the origin so the last position coming out of the for loop can be used to compute the angle $\alpha$ using $$\alpha=\arctan \frac{y_{n_{steps-1}}}{x_{n_{steps-1}}}.$$ Note that $\arctan$ gives angle values such that $\alpha \in [-\frac{\pi}{2},\frac{\pi}{2}]$. So for negative angles which are in $[-\frac{3\pi}{2},-\frac{\pi}{2}]$ we need to subtract $\pi$ to $\arctan$.

In [5]:
ri=np.array([-9e-3,1e-3,0])
rp=np.zeros(3)
vi=np.linspace(1000,5000,40)
angle=np.zeros(len(vi))
time_frame=np.array([0,50e-6])
n_steps=20000
for i in range(len(vi)):
v=np.array([vi[i],0,0])
loc=ODE_integration_E(rp,ri,v,time_frame,n_steps=n_steps)
angle[i]=np.arctan(loc[1,n_steps-1]/loc[0,n_steps-1])
if (loc[0,n_steps-1]<0):
angle[i]=-np.pi+angle[i]
angle=abs(angle)


We now use the optimize package from scipy to do our curve fitting. It is important to realize that the fitting has to be done properly. This is why we use high velocities for vi above. Using low velocities would

In [6]:
from scipy.optimize import curve_fit
def func(x,a,b):
return a*x**b
fig, ax = plt.subplots(1, 1, figsize=(8, 8))
popt, pcov = curve_fit(func, vi, angle)
ax.plot(vi,angle,'r+',label='computed data')
ax.plot(vi,func(vi, *popt), 'g--',label='fit: a=%5.3f, b=%5.3f' % tuple(popt))
plt.legend()
plt.xlabel("Initial velocity (m/s)")
plt.show()


We see that the fit tells us that the angle varies as $\frac{1}{v^2}$

Problem 3¶

We define our function using numpy arrays. Clearly this method is most readable. It is not the most efficient though. However it is important to try to depart from the C++ or Java formalism to really benefit from Python.

In [7]:
def comp_E(rp,box,n,q=-1.6e-19):
X,Y=np.meshgrid(np.linspace(box[0],box[1],n[0]),np.linspace(box[2],box[3],n[1]))
RX=X-rp[0]
RY=Y-rp[1]
R=np.sqrt(RX**2+RY**2)
Ex=1./(4*np.pi*8.85e-12)*q/R**2*RX/R
Ey=1./(4*np.pi*8.85e-12)*q/R**2*RY/R
E=np.hypot(Ex,Ey)
return X,Y,Ex,Ey,E

X,Y,Ex,Ey,E=comp_E((0,0),(-1,1,-1,1),(80,80))


Note that we had to clip the electric field strength for our color plot so that we can appreciate the color scale.

In [8]:
fig, ax = plt.subplots(1, 1, figsize=(8, 8))
strm=ax.streamplot(X,Y,Ex,Ey,color=np.clip(E,0,E.max()/30),linewidth=1, cmap=plt.cm.nipy_spectral,
density=2, arrowstyle='->', arrowsize=1.5)
ax.axis('equal')
cbar=fig.colorbar(strm.lines)
cbar.ax.set_ylabel('E(V/m)', rotation=0)
plt.show()


If you do not want to clip it, then use the $\log_{10}$ scale. In case you are not familiar with it $\log_{10}(10^x)=x$. Also note that the color scale below is a really good choice when the color spread is rather flat, as it is often the case with logarithms, because each shade has been desaturated.

In [9]:
fig, ax = plt.subplots(1, 1, figsize=(8, 8))
strm=ax.streamplot(X,Y,Ex,Ey,color=np.log10(E),linewidth=1, cmap=plt.cm.nipy_spectral,
density=2, arrowstyle='->', arrowsize=1.5)
ax.axis('equal')
cbar=fig.colorbar(strm.lines)
cbar.ax.set_ylabel('$\log_{10}$E (V/m)', rotation=0)
plt.show()

In [ ]: