due 09/12/2018 at 11:59pm
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
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.
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.
import math
from math import sqrt as sqrt
import numpy as np
import matplotlib.pyplot as plt
#to display inside the notebook!
%matplotlib inline
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
.
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.
We define our functions. comp_E
compute the electric field created at the location ro
by a charge qp
located at a point rp
.
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
.
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$.
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
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.ylabel("Deflection angle (rad)")
plt.show()
We see that the fit tells us that the angle varies as $\frac{1}{v^2}$
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.
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.
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.
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()