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.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 |
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' ) |