I'm quite new to programming, and I'm trying to compare the decay of atoms with a Monte Carlo random.random() function against an analytical solution. I can't seem to figure out how to print the graph for the Monte Carlo solution though.
import numpy as np
import matplotlib.pyplot as plt
import random
"""
Define Constants
"""
Lambda = 0.4 #This probability that an atom will decay, where 0 means no decays and 1 means 100% chance of decay
dt = 1 #This is the step size, time should be adjusted depending on time unit of lambda
NumberOfAtoms = 1000 #This is the starting amount of atoms in our sample
"""
Introduce arrays to hold values for the time and amount of atoms we have for each step iteration
"""
Time = []
Atoms = []
while NumberOfAtoms >=0:
"""
While giving the statement that 'WHILE' the number of atoms is above or equal to 0, we perform
the following
"""
NumberOfDecays = 0 #The atom will stop decaying (decay is 0) when the NumberOfAtoms = 0
for Atom in range(NumberOfAtoms):
"""
Here we are using the random number generator to oppose the chance of a decay, lambda
"""
r = random.random()
if (r < Lambda): #If random generated number is not over 0.4, the atom will not decay
NumberOfDecays += 1 #this action is repeated until the 'WHILE' statement is fulfilled
NumberOfAtoms = NumberOfDecays
if NumberOfAtoms < 1:
print 'The atoms have all decayed'
break
print 'The number of atoms left is', NumberOfAtoms
'''Here we are defining variables for a new function, the Analytical solution to the decay equation'''
t1 = 0
tEnd = 20
xdata =[]
yana= []
xxdata =[]
Nana = []
NewAtoms = NumberOfAtoms
AtomsInitial = 1000 #Setting same initial atom count to match the Monte Carlo Solution
for t in np.arange(t1, tEnd, dt):
xdata.append(t)
yana.append(AtomsInitial*np.exp(-Lambda*t) )
#Atoms.append(NumberOfAtoms)
plt.plot(xdata,yana, label = 'Analytical Solution')
plt.title('Number of Atoms vs Time')
plt.xlabel('Time')
plt.ylabel('Number of Atoms')
plt.legend(loc='upper right')
