娘核種1つ 娘核種２つ 核種３つ以上 Cs-137 ルンゲ=クッタで解く
dN2/dt = λ1・N1 - λ2・N2 —娘核種
dN1/dt = – λ1・N1 —親核種
(画像を押すと拡大)
Option Explicit
Sub 放射性物質_減衰()
Dim t As Double, dt As Double '時間 Exel VBA
Dim N(2) As Double 'mol数
Dim λ(2) As Double '崩壊定数
Dim ans As Variant
Dim cnt1 As Long
Dim sh As Worksheet
Dim 終了時刻 As Double
Set sh = ThisWorkbook.Worksheets(1)
dt = 0.001
'初期条件----------------------------------
N(1) = 1 '1mol 親
N(2) = 0 '0mol 娘
λ(1) = 1 '崩壊定数 親
λ(2) = 0.5 '崩壊定数 娘
cnt1 = 0
終了時刻 = 5
'------------------------------------------
With sh '表記
.Cells(9, 4).Value = "t" 't
.Cells(9, 5).Value = "N1" 'N1 mol
.Cells(9, 6).Value = "N2" 'N2 mol
End With
sh.Range("D10:F500000").ClearContents
Do While t <= 終了時刻
With sh
.Cells(cnt1 + 10, 4).Value = t 't
.Cells(cnt1 + 10, 5).Value = N(1) 'N1 mol
.Cells(cnt1 + 10, 6).Value = N(2) 'N2 mol
End With
cnt1 = cnt1 + 1
ans = RK4(t, dt, λ, N) 'ルンゲ＝クッタ
N(2) = ans(2)
N(1) = ans(1)
t = t + dt
Loop
sh.Range("N3") = Time
End Sub
Function RK4(ByVal t, ByVal dt, ByVal λ, ByVal N) As Variant
'下準備-----------------------------------------------------
Dim LL(2) 'λ
Dim NN1(2), NN2(2), NN3(2), NN4(2) 'N
Dim k1NN(2), k2NN(2), k3NN(2), k4NN(2) 'k1NN--k4NN(識別番号）
LL(1) = λ(1)
LL(2) = λ(2)
'下準備-----------------------------------------------------
'k1
NN1(1) = N(1) '前処理
NN1(2) = N(2) '前処理
k1NN(2) = dt * 娘関数(t, LL, NN1)
k1NN(1) = dt * 親関数(t, LL, NN1)
'k2
NN2(2) = NN1(2) + k1NN(2) / 2 '前処理
NN2(1) = NN1(1) + k1NN(1) / 2 '前処理
k2NN(2) = dt * 娘関数(t + dt / 2, LL, NN2)
k2NN(1) = dt * 親関数(t + dt / 2, LL, NN2)
'k3
NN3(2) = NN1(2) + k2NN(2) / 2 '前処理
NN3(1) = NN1(1) + k2NN(1) / 2 '前処理
k3NN(2) = dt * 娘関数(t + dt / 2, LL, NN3)
k3NN(1) = dt * 親関数(t + dt / 2, LL, NN3)
'k4
NN4(2) = NN1(2) + k3NN(2) '前処理
NN4(1) = NN1(1) + k3NN(1) '前処理
k4NN(2) = dt * 娘関数(t + dt, LL, NN4)
k4NN(1) = dt * 親関数(t + dt, LL, NN4)
NN1(2) = NN1(2) + 1 / 6 * (k1NN(2) + 2 * k2NN(2) + 2 * k3NN(2) + k4NN(2))
NN1(1) = NN1(1) + 1 / 6 * (k1NN(1) + 2 * k2NN(1) + 2 * k3NN(1) + k4NN(1))
RK4 = NN1
End Function
Function 娘関数(ByVal t, ByVal LL, ByVal NN) As Double
'dN2/dt＝λ1・Ｎ1 －λ2・Ｎ2
娘関数 = LL(1) * NN(1) - LL(2) * NN(2)
End Function
Function 親関数(ByVal t, ByVal LL, ByVal NN) As Double
'dN1/dt = －λ1・N1
親関数 = -LL(1) * NN(1)
End Function