Python har en randomgenerator som del af numpy pakken
import numpy as np
Ovenstående linje importerer en funktion
np.random som har mange funktioner til at generere tilfældige tal på. Det er denne funktion vi bruger her.
Links til programmer og hjælpesider
I dette program har jeg 500 kerner (n_0) som kan henfalde med en sandsynlighed (henfaldskonstant) k (her valgt til 0,01).
Tiden er sat til 800s.
For hver sekund spørger jeg alle resterende kerner om de er henfaldet. Det gøres ved at jeg genererer et tilfældigt tal.
tal= np.random.random()
n_kerner er det antal kerner som IKKE er henfaldet.
j-loopet kører 800 gange for 800 sekunder. i-loopet kører n_kerner gange dvs at den for hvert sekund spørger ALLE resterende (n_kerner) kerner om de er henfaldet.
Kernen henfalder HVIS if tal < k . Random funktionen genererer et tal mellem 0 og 1 med mange decimaler så det er ligepræcis 1% af den (k = 0.01) som opfylder kriteriet til henfald.
Grafen ser ud som vist her til højre. Læg mærke til at man kan aflæse en halveringstid
Find siderne om plotfunktioner og tilføj en legende som viser antal kerner og henfaldskonstanten k
Design n tilføjelse til programmet sådan at halveringstiden bliver beregnet af programmet.
import numpy as np
import matplotlib.pyplot as plt
n_0 =500
n_kerner = n_0
antalgange = 800
k = 0.01
nlist=[]
tlist=[]
for j in range(antalgange):
for i in range(n_kerner):
tal = np.random.random()
if tal < k:
n_kerner=n_kerner-1
nlist.append(n_kerner)
tlist.append(j)
plt.xlabel("tid(s)")
plt.ylabel("antal henfald")
plt.axis([-10,antalgange+20,-10,n_0+50])
plt.plot(tlist,nlist,linestyle='dotted',color='blue')
Jeg synes selv at det er en ret finurlig måde at gøre det her på. Og jeg har det med for at vise hvor forskelligt man kan gøre ting i python.
I begge løsninger kan man tilføje en regression som enten kan være lineær (så skulle man tage ln(x) af alle y-værdier. Eller så kunne man tilføje en eksponentiel regression.
import numpy as np
import matplotlib.pyplot as plt
#definitionerne
N = 1000 # tusinde kerner til at starte med
k = 0.05 # henfaldskonstanten
tid = 200 # antal gange vi vil spørge
#interne definitioner
Nlist = [] #definerer en tom liste som bliver fyldt op med 1 2 3 4 ....1000
Nantal = [] #bliver fyld op med de tilbageværende kerner 1000 998 995 990 .. 0
tidlist = [] #x-akse til plot af tid 1 2 3 4 .... 200
# producere en liste med N tal 1 .. N
i = 1
while i <= N: #opfyldning af Nlist med tal 1 2 3 4 .... N
Nlist.append(i)
i +=1
i = 1
while i <= tid: #opfyldning af tidlist som er x-aksen med 1 2 3 .. tid
tidlist.append(i)
i +=1
#nu løber vi listen igennem og spørger hver kerne om den er henfaldet
j = 1
while j <= tid: #vi spørger tid (200) antal gange
i = 1
while i < len(Nlist): #det er længden af Nlist som angiver antal ikke henfaldne kernerne
rand = np.random.randint(100000) #jeg genererer et stort random tal
if rand < (100000*k): #k er 0,05 så kun hvis det genererede tal er mindre kører "if"
Nlist.pop() #vi skærer (popfunktion) en af listen så len(Nlist) bliver en mindre
i += 1 #i loopet kører indtil alle tilbægeværende kerner er spurgt
Nantal.append(len(Nlist)) #vi skriver antal af kerner tilbage til listen
j +=1 #j-loopet skal køre tid gange. (200)
minlegende = ("N = {0} k = {1} tid = {2}s ".format(N,k,tid))
plt.plot(tidlist,Nantal,'r+',label=minlegende) #ro betyder rød og cirkel fx gx virker også
plt.xlabel("tid (s)")
plt.ylabel("N (antal)")
plt.legend(loc='upper right')
plt.title('Radioaktiv henfald med random generator')
plt.show()
Her kan du se indholdet af Nantal som jo indeholder alle y-værdier til
[954, 913, 870, 824, 792, 760, 737, 705, 678, 642, 610, 584, 557, 524, 499, 474, 452, 432, 404, 386, 362, 336, 325, 309, 298, 279, 262, 249, 241, 228, 217, 208, 200, 187, 180, 169, 158, 151, 146, 143, 140, 134, 126, 120, 115, 114, 109, 104, 102, 98, 91, 86, 79, 74, 73, 70, 66, 63, 56, 56, 52, 52, 49, 47, 44, 42, 36, 36, 35, 35, 32, 30, 27, 26, 26, 23, 20, 19, 17, 17, 15, 15, 15, 13, 11, 11, 11, 11, 9, 9, 9, 9, 9, 9, 8, 7, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 4, 4, 3, 3, 3, 3, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]
Her kan du se indholdet af tidlist som indeholder alle x-værdier (kommer ikke bag på os :-) )
[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, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100, 101, 102, 103, 104, 105, 106, 107, 108, 109, 110, 111, 112, 113, 114, 115, 116, 117, 118, 119, 120, 121, 122, 123, 124, 125, 126, 127, 128, 129, 130, 131, 132, 133, 134, 135, 136, 137, 138, 139, 140, 141, 142, 143, 144, 145, 146, 147, 148, 149, 150, 151, 152, 153, 154, 155, 156, 157, 158, 159, 160, 161, 162, 163, 164, 165, 166, 167, 168, 169, 170, 171, 172, 173, 174, 175, 176, 177, 178, 179, 180, 181, 182, 183, 184, 185, 186, 187, 188, 189, 190, 191, 192, 193, 194, 195, 196, 197, 198, 199, 200]