娘核種1つ 娘核種２つ 核種３つ以上 Cs-137 ルンゲ=クッタで解く
（画像を押すと拡大） 参考： housya_sample.zip (マクロを使うので自己責任でお願いします）
dN1/dt = －λ1・N1 —親核種
dN2/dt = 0.8＊λ1・N1 -λ2・N2 —娘1核種
dN3/dt = 0.2＊λ1・N1 -λ3・N3 —娘2核種
dN4/dt = ＋λ2・N2 ＋λ3・N3 -λ4・N4 —孫核種
dN5/dt = ＋λ4・N4 —ひ孫核種
Option Explicit
Sub 放射能_娘ex2() 'ExelVBA
Dim t As Double, dt As Double '時間
Dim KK As Long
Dim N(10) As Double 'mol数 10核種
Dim λ(10) 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) = Range("F4") '1mol 親
N(2) = Range("F5") '0mol 娘
N(3) = Range("F6") '0mol 娘2
N(4) = Range("F7") '0mol 娘3
N(5) = Range("F8") '0mol 娘4
λ(1) = Range("I4") '崩壊定数 親
λ(2) = Range("I5") '崩壊定数 娘1
λ(3) = Range("I6") '崩壊定数 娘2
λ(4) = Range("I7") '崩壊定数 娘3
λ(5) = Range("I8") '崩壊定数 娘4
cnt1 = 0
終了時刻 = 5
'------------------------------------------
Call クリア
Do
With sh
.Cells(cnt1 + 12, 27).Value = t 't
.Cells(cnt1 + 12, 28).Value = N(1) 'N1 mol
.Cells(cnt1 + 12, 29).Value = N(2) 'N2 mol
.Cells(cnt1 + 12, 30).Value = N(3) 'N3 mol
.Cells(cnt1 + 12, 31).Value = N(4) 'N4 mol
.Cells(cnt1 + 12, 32).Value = N(5) 'N5 mol
End With
cnt1 = cnt1 + 1
ans = RK4(t, dt, λ, N) 'ルンゲ＝クッタ
For KK = 1 To 10
N(KK) = ans(KK)
Next KK
t = t + dt
Loop While (t - dt) < 終了時刻
End Sub
Function RK4(ByVal t, ByVal dt, ByVal λ, ByVal N) As Variant
'下準備-----------------------------------------------------
Dim LL(10) 'λ
Dim NN1(10), NN2(10), NN3(10), NN4(10) 'N
Dim k1NN(10), k2NN(10), k3NN(10), k4NN(10) 'k1NN--k4NN(識別番号）
Dim KK As Long
For KK = 1 To 10
LL(KK) = λ(KK)
Next KK
'下準備-----------------------------------------------------
'k1
For KK = 1 To 10
NN1(KK) = N(KK) '前処理
Next KK
For KK = 1 To 10
k1NN(KK) = dt * 娘Gen関数(KK, t, LL, NN1)
Next KK
'k2
For KK = 1 To 10
NN2(KK) = NN1(KK) + k1NN(KK) / 2 '前処理
Next KK
For KK = 1 To 10
k2NN(KK) = dt * 娘Gen関数(KK, t + dt / 2, LL, NN2)
Next KK
'k3
For KK = 1 To 10
NN3(KK) = NN1(KK) + k2NN(KK) / 2 '前処理
Next KK
For KK = 1 To 10
k3NN(KK) = dt * 娘Gen関数(KK, t + dt / 2, LL, NN3)
Next KK
'k4
For KK = 1 To 10
NN4(KK) = NN1(KK) + k3NN(KK) '前処理
Next KK
For KK = 1 To 10
k4NN(KK) = dt * 娘Gen関数(KK, t + dt, LL, NN4)
Next KK
'計
For KK = 1 To 10
NN1(KK) = NN1(KK) + 1 / 6 * (k1NN(KK) + 2 * k2NN(KK) + 2 * k3NN(KK) + k4NN(KK))
Next KK
RK4 = NN1
End Function
Function 娘Gen関数(ByVal Gen, ByVal t, ByVal LL, ByVal NN) As Double
Dim ans As Variant
If Gen = 1 Then '親
'dN1/dt = －λ1・N1
ans = -LL(1) * NN(1)
ElseIf Gen = 2 Then '娘
'd/dt・N2（t)＝0.8*λ1・Ｎ1 －λ2・Ｎ2
ans = 0.8 * LL(1) * NN(1) - LL(2) * NN(2)
ElseIf Gen = 3 Then '娘2
'd/dt・N3（t)＝0.2*λ1・Ｎ1 －λ3・Ｎ3
ans = 0.2 * LL(1) * NN(1) - LL(3) * NN(3)
ElseIf Gen = 4 Then '娘3
'd/dt・N4（t)＝ +λ2・Ｎ2 +λ3・Ｎ3 -λ4・Ｎ4
ans = LL(2) * NN(2) + LL(3) * NN(3) - LL(4) * NN(4)
ElseIf Gen = 5 Then '娘4
'd/dt・N5（t)＝ +λ4・Ｎ4
ans = LL(4) * NN(4)
End If
娘Gen関数 = ans
End Function
Sub クリア()
Dim sh As Worksheet
Set sh = ThisWorkbook.Worksheets(1)
sh.Range("AA12:AK500000").ClearContents
End Sub