Radio-random

Simulering af radioaktiv henfald

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

Den nemme løsning

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

OPGAVE 1:

Find siderne om plotfunktioner og tilføj en legende som viser antal kerner og henfaldskonstanten k

OPGAVE 2:

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')

En spændende men lidt mere kompliceret løsning med lister

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]

Page updated

Google Sites

Report abuse