３体問題

温度シミュレーション (下部のみ１００℃　周りは０℃の例）

Option Explicit
Option Base 0
Private Const XX = 41
Private Const YY = 41
Private T(YY, XX)          '　　Tはこのプロシージャ内で共通化している

Sub 境界条件(T)
Dim x As Integer
Dim y As Integer
Erase T
For x = 0 To XX - 1      '　　 上面
  T(0, x) = 0
Next x
For x = 0 To XX - 1      '　 　下面だけ100度
    T(YY - 1, x) = 100
Next x
For y = 0 To YY - 1      '　　 左面
  T(y, 0) = 0
Next y
For y = 0 To YY - 1      '　　 右面
  T(y, XX - 1) = 0
Next y
End Sub

Sub メイン()
Dim i As Integer
Dim j As Integer
Dim Max誤差 As Double
Dim 旧data As Double
Call 境界条件(T)
Do
Max誤差 = 0
For i = 1 To YY - 2
    For j = 1 To XX - 2
        旧data = T(i, j)      '　前データを保管
        T(i, j) = (T(i + 1, j) + T(i - 1, j) + T(i, j + 1) + T(i, j - 1)) / 4    '　まわりの平均値　２次元なら１／４、３次元なら１／６
         '　前データとの差を見て差が一番大きいものをMax誤差とする
         If Max誤差 < Abs(T(i, j) - 旧data) Then　
            Max誤差 = Abs(T(i, j) - 旧data)
         End If
    Next j
Next i
Loop While Max誤差 > 0.001　　　'　差が収束したら終了。

'描画する
For i = 0 To YY
    For j = 0 To XX
      Cells(i + 5, j + 5) = T(i, j)  
    Next j
Next i
End Sub

熱拡散方程式（全海洋蒸発後、 １年間　地表面で４０００℃が続くとしたモデルの例）

　初期値：地中４０℃

（温度拡散率は砂の係数を使った　：　χ＝２×１０^(-7)）

Option Explicit
Option Base 0
Sub 一次元拡散方程式()

Range("E5:P1040000").ClearContents
 
Const Xmax = 21                 '　メートル
Const tmax = 350001           '　秒

'　砂の係数　2 * 10 ^ (-7)
Const χ = 2 * 10 ^ (-7) * 100  '温度拡散率　　１００秒を１単位とした
Const 描画圧縮率 = 10000

Dim i As Long, j As Long
Dim T(Xmax, tmax)

Erase T

For j = 0 To (tmax - 1)      '地表表面を4000℃
    T(0, j) = 4000
Next j
   
For j = 0 To (tmax - 1)      '地中を40℃
   For i = 1 To Xmax
    T(i, j) = 40
   Next i
Next j

For j = 0 To (tmax - 1)
  For i = 1 To (Xmax - 1)
    '両端は除いて計算
    T(i, j + 1) = T(i, j) + χ * (T(i + 1, j) + T(i - 1, j) - 2 * T(i, j))
  Next i
Next j

For j = 0 To (tmax - 1)
  If Int(j / 描画圧縮率) = j / 描画圧縮率 Then
            Cells(j / 描画圧縮率 + 5, 5) = j
        For i = 0 To (Xmax - 1)
            Cells(j / 描画圧縮率 + 5, i + 6) = T(i, j)   　　 ' T(x,t)
        Next i
  End If

Next j

End Sub