Vi vil nu prøve at lave et eksponentielt fit til et randomgenereret radioaktiv henfaldsspektrum.
Vi tager udangspunkt i programmerne
Programmet radioaktiv henfald er vist her til højre.
For at kombinere de to programmer skal det OUTPUT fra programmet til højre kunne læses af programmet "eksponentiel fit".
Prøv selv at kombinere de to programmer.
Løsningen er her nedenunder hvis du går istå. Husk at der er mange måder at gøre det her på.
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')
Det første problem er at outputtet fra henfaldsprogrammet af af typen "list" og input af fitprogrammet er af typen "array"
Et hurtigt opslag på nettet giver et løsningsforslag som vist her til højre.
import numpy as np
my_list = [2,4,6,8,10]
my_array = np.array(my_list)
# printing my_array
print my_array
Gennemgå koden og slet irrelevante plotfunktioner. Sørg for at variablerne er det samme i begge programmer.
x_vaerdier = tlist_array
y_vaerdier = nlist_array
import numpy as np
import matplotlib.pyplot as plt
n_0 =500
n_kerner = n_0
antalgange = 300
k = 0.02
nlist=[]
tlist=[]
#tall = random.random()
#print(tall)
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')
#plt.axis([-10, 500, 0.5 * np.min(tlist), 1.1 * np.max(tlist)])
#konvertere list til array sådan at fitteprogrammet kan læse dem.
nlist_array = np.array(nlist)
tlist_array = np.array(tlist)
#nye definitioner til fitteprogrammet
from scipy.optimize import curve_fit
def func(x, a, b, c):
return a * np.exp(-b * x) + c
x_vaerdier = tlist_array
y_vaerdier = nlist_array
popt, pcov = curve_fit(func, x_vaerdier, y_vaerdier, p0=[250, 0.13, 0])
a=popt[0].round(3)
b=popt[1].round(5)
c=popt[2].round(3)
minlegende=('y = {0} * exp(-{1}x) + {2}'.format(a,b,c))
x_plot=np.linspace(0,antalgange,100)
plt.plot(x_plot, func(x_plot, *popt), 'r-',label=minlegende)
plt.scatter(x_vaerdier,y_vaerdier,label='datapunkter')
plt.title("x vs y = a*exp(-bx) +c")
plt.xlabel('x-værdier')
plt.ylabel('y-værdier')
plt.legend()
plt.show()
# equation popt indeholder a b og c i en liste med mange decimaler
# round runder til 2 decimaler og print(f'.....') fungerer som .format
print(f'The equation of regression line is y={a}exp(-{b}x)+{c}')
halveringstid=np.log(2)/b
t_en_halv=halveringstid.round(3)
print(f'halveringstid = {t_en_halv} s')
Programmets output er her
The equation of regression line is y=499.834exp(-0.02143x)+0.236
halveringstid = 32.345 s
Her kan du se at det er en rigtig simulering.
a-værdien som er antal kerner til at begynde med er 499.834 meget tæt på
n_0 = 500
b-værdien er 0.02143 som er tæt på de
k = 0.02
c-værdien har ikke rigtig nogen fysisk betydning.
Vi ved at der gælder :
Og siden henfaldskonstanten k = 0.02 i forsøget og genfindes som 0.02143 efter simuleringen kan vi beregne halveringstiden til
t(1/2) = 32,345s
Prøv at aflæse grafen og gør dig klar at det er et realistisk bud.
Ændr startværdierne for "n_0" og "antalgange" diskutér inflydelsen på resultaterne.
Her til højre ser du
n_0 = 30
n_kerner = n_0
antalgange = 60
k = 0.1
sammenlign med ovenstående.
Forstil dig at du skal lave et rigtigt forsøg med en radioaktiv kilde og et Geiger Müller rør. Hvilke fysiske parametre svaerer n_0 og antalgange til ?