Background

In this lab we will be using ODE or ordinary differential equation, is when we use a function and also apply its derivatives, to show the relationship between a function and its rate of change. In this lab we will also use a method called Taylor expansion is used for an approximation of a function, this is done by using an infinite sum of terms, and its derivatives for said function of x. This method is useful, when the equations or functions cannot be solved analytically. In this lab we will also be using Euler method, which will come in handy when we are trying to solve ODE equations.

Euler method in radioactive decay

we will use the equation that describes how the number of surviving nuclei changes with time is:

dN / dt=−N/τ where N = number of nuclei present τ = half live of the nuclei

next step will be to solve the equation by separating variables N(t)=N0e^−t/τ this being the analytic solution with respect to time.

Using Euler method N(t+Δt) ≈ N(t) − (N/τ)(Δt)

we will first have to create a function called initialize this function will create a nstep which is tmax/dt

def intialize(N0,Tao,tmax,dt):
  n_steps = int(tmax/dt)
  t = np.zeros(n_steps)
  n_nuclei = np.zeros(n_steps)
  n_nuclei[0] = N0
  return t,n_nuclei

then will need a calculate function will loop over arrays and calculate
n_nuclei[i] = n_nuclei[i-1] - n_nuclei[i-1]/tau * dt t[i] = t[i-1] + dt and will return t,n_nuclei arrays

'''n_nuclei[i] = n_nuclei[i-1] - n_nuclei[i-1]/tau * dt
t[i] = t[i-1] + dt
'''
def calcuate(t,n_nuclei,dt,tau):
  for i in range(1,len(t)):
    n_nuclei[i] =n_nuclei[i-1]-n_nuclei[i-1]/tau*dt
    t[i] = t[i-1]+dt
  return t,n_nuclei

Calc with 100 nuclei, decay-life 1 second, tmax 5 seconds and timestep .05

# define initial values
N0 = 100
tau =1
tmax =5
dt =.05

# call initialize function
t,n_nuclei= intialize(N0,tau,tmax,dt)

# call your calculate function
t,n_nuclei = calcuate(t,n_nuclei,dt,tau)
print(t)
print(n_nuclei)

plt.plot(b[0],b[1],'ks')
plt.plot(a,N0*np.exp(-a/tau),'k-')
plt.show()

Changing the timestamps

When viewing the code below you will see that we change the time steps from dt=.5 to dt=.05 and you will see howt he accuracy begins to increase when we change the time steps, in comparison to the analytic answer we have created with the python code

##########################################################
## code to get number of nuclei and time for dt = 0.05 s
##########################################################
N0 = 100
tau =1
tmax =5
dt=.05

# initialize the arrays
tp05,np05 = intialize(N0,tau,tmax,dt)

# call calculate, and save results to tp05, np05 (time for dt=.05s, number for dt=.05s)
newcalc_1 = calcuate(tp05,np05,dt,tau)
print(tp05,np05)
##########################################################
## code to get number of nuclei and time for dt = 0.2 s
##########################################################
dt_2=.2

# initialize the arrays
tp2,np2 = intialize(N0,tau,tmax,dt_2)

# call calculate, and save results to tp2, np2 (time for dt=.2s, number for dt=.2s)
newcalc_2 = calcuate(tp2,np2,dt_2,tau)

##########################################################
## code to get number of nuclei and time for dt = 0.5 s
##########################################################

dt_3 = .5
# initialize the arrays
tp3,np3 = intialize(N0,tau,tmax,dt_3)
# call calculate, and save results to tp5, np5 (time for dt=.5s, number for dt=.5s)
newcalc_3 = calcuate(tp3,np3,dt_3,tau)

# code to plot the 3 different results
plt.plot(tp05,np05,'ks')
plt.plot(tp2,np2,'rs')
plt.plot(tp3,np3,'bs')
# code to plot the analytic solution
plot_analytic = N0*np.exp(-tp05/tau)
plt.plot(tp05,plot_analytic,'k-')
# include a legend
plt.legend(['dt=.05','dt=.2','dt=.5','Analytic'])
# label your axes
plt.xlabel('Time')
plt.ylabel('Number of Nuclei')
plt.title('Number of Nuclei vs Time')

Another example Euler method on air resistance

def acceleration(P,m,v):
  a = P/(m*v)
  return a

def initialize(v0,tmax,dt):
  nsteps = int(tmax/dt)
  t = np.zeros(nsteps)
  v = np.zeros(nsteps)
  v[0] = v0
  return t,v

def calculate(t,v,dt,P,m):
  # create i in range of length
  for i in range(1,len(t)):
    v[i] = v[i-1] + acceleration(P,m,v[i-1]) * dt
    t[i] = t[i-1] + dt
  return t,v


P = 400
m = 70
v0 = 4
dt = .1
tmax = 200

# call initialize and save output as t,v
t,v= initialize(v0,tmax,dt)

# call calculate, and save output as tnoair,vnoair
tnoair,vnoair = calculate(t,v,dt,P,m)


###############################################################
# PLOT YOUR RESULTS
###############################################################
# open a figure, and make your figure size (6,6)

plt.figure(figsize=(6,6))
plt.plot(t,v,tnoair,vnoair,label='no air')
plt.xlabel("time")
plt.ylabel('velocity')
plt.title('Velocity vs Time')
plt.legend(['no air'])
plt.show()

# plot velocity vs time with a line
# be sure to label the x and y axes
# include a legend and figure title

B2=0.5*0.5*1.225*0.33

def accel_resist(p,m,v,B2):
  a = p/(m*v) - (B2/(m))*v**2
  return a

def calcuate_resist(t,v,dt,P,m,B2):
  for i in range(1,len(t)):
    v[i] = v[i-1] + accel_resist(P,m,v[i-1],B2)*dt
    t[i] = t[i-1] + dt
  return t,v

P = 400
m = 70
v0 =4
dt =.1
tmax =200
rho = 1.225
C = .5
A = .33

# calculate B2
B2 = .5*C*rho*A

# call intialize from Part 1, save output as tnoair,vnoair
tnoair,vnoair = initialize(v0,tmax,dt)
# calculate the Euler method with no air resistance
tnoair,vnoair = calculate(tnoair,vnoair,dt,P,m)
# save output as tnoair, vnoair

# call intialize, save output as tair,vair
tair,vair = initialize(v0,tmax,dt)
# calculate the Euler method WITH air resistance
tair,vair = calcuate_resist(tair,vair,dt,P,m,B2)
# save output as tair, vair


###############################################################
# PLOT YOUR RESULTS
###############################################################

# open a figure, and make your figure size (6,6)

# plot velocity with no air resistance vs time with a solid line

# plot velocity with air resistance vs time with a dashed line

# be sure to label the x and y axes

# include a legend and figure title
plt.figure(figsize=(6,6))
plt.plot(tnoair,vnoair,label='no air')
plt.xlabel("time")
plt.ylabel('velocity')
plt.title('Velocity vs Time')
plt.legend(['no air'])
plt.plot(tair,vair,label='air')
plt.xlabel("time")
plt.ylabel('velocity')
plt.title('Velocity vs Time')
plt.legend(['no air'])
plt.show()