3．ルンゲ・クッタ法のサンプル

　　　（プログラム（マクロ）の提供を主としている）

 　ルンゲ・クッタ法の微係数のテイラー展開の

　4次の項まで考慮の式は、以下の書籍より出典。

　丸善　木村欽一著

　エクセルで解く線形・非線形方程式の数値計算

　木村先生ありがとうございました。

3－1.熱（定係数1階線形方程式）

実験の系の微分方程式

Temp_end / R = C・dθ/dt + θ / R

dθ/dt = Temp_end / C / R - θ / C / R

ここで　θ：温度[K]     C:熱容量[Ws/K]     R : 熱抵抗[K/W]

Excel macro        sub expCR()

       Sub expCR()

       tend = Range("B5")

       step = Range("B6")

       V = Range("B7")

       c = Range("B8")

       rhou = Range("B9")

       A = Range("B10")

       alpha = Range("B11")

       Temp_end= Range("B12")

       temp0 = Range("B13")

       [F5].Select

       cc = V * c * rhou

       ActiveCell.Offset(0, 0) = cc

       R = 1# / A / alpha

       ActiveCell.Offset(1, 0) = R

       dt = tend / step

       ActiveCell.Offset(2, 0) = dt

       [O1].Select

       t = 0#

       temp = temp0

       For i = 1 To step

       Call rk2(dt, t, temp, tempnew, cc, R, Temp_end)

       ActiveCell.Offset(i, 0) = t

       ActiveCell.Offset(i, 1) = temp

       t = t + dt

       temp = tempnew

       Next i

       End Sub

      Sub rk2(dt, t, temp, tempnew, cc, R, Temp_end)

      k1 = dt * f(t, temp, cc, R, Temp_end)

      k2 = dt * f(t + dt / 2, temp + k1 / 2, cc, R, Temp_end)

      k3 = dt * f(t + dt / 2, temp + k2 / 2, cc, R, Temp_end)

      k4 = dt * f(t + dt, temp + k3, cc, R, Temp_end)

      tempnew = temp + (k1 + 2 * (k2 + k3) + k4) / 6

      End Sub

      Function f(t, temp, cc, R, Temp_end)

      f = Temp_end / cc / R - temp / cc / R

      End Function