Sample Scripts from GB Books
GB005:  1次元/2次元過渡問題

今度はそれぞれ200度K、400度Kに保った鉄棒A, Bを時刻 t = 0 で接触させたときの温度分布の変化を調べてみることにします。ただし鉄棒の周囲はすべて熱的に遮断されているものとします。この例でキーとなるのは初期状態の設定に際し ustep というステップ関数（unit step function、Heavi-side function）を用いる点です。

例えば ustep(x - 0.5) という関数は x < 0.5 のとき0、x > 0.5 のときに1という値を取る右図のような関数となります。

基本形は transient01a.pde と変わりません。
  TITLE
    'Temperature Step in Bar'    { transient01b.pde }

  SELECT
    Errlim = 1e-4

  VARIABLES
    temp                  { Temperature }

  DEFINITIONS
    Lx = 1.0  Ly = 0.1
    heat = 0                  { Heat source }
    k = 82  rcp = 7.87e3*449  { Iron parameters }
    fluxd_x = -k*dx(temp)  fluxd_y = -k*dy(temp)
    fluxd = vector(fluxd_x, fluxd_y)  fluxd_m = magnitude(fluxd)

初期状態をステップ関数を用いて設定します。
  INITIAL VALUES
    temp = 200 + 200*ustep(x - Lx/2)  { Unit step function }

  EQUATIONS
    div(fluxd) + rcp*dt(temp) = heat

境界条件としてはすべて断熱型を使用します。境界条件がすべてNatural型だった場合には関数値がユニークに決まらないことがありますが（GB001参照）、今回は初期値を指定しているため問題は起こりません。
  BOUNDARIES
    Region 1
      Start(0,0)
        Natural(temp) = 0 
        Line to (Lx, 0) to (Lx, Ly) to (0, Ly) to Close

  TIME
    from 0 to 100000

プロットを出力するタイミングは対数スケールを意識して設定します。
  PLOTS
    for t = 10, 30, 100, 300, 1000, 3000, 10000, 30000, 100000
    Elevation(temp) from (0, Ly/2) to (Lx, Ly/2)
    Contour(temp) painted

  END

今回、プロットを出力するタイミングは for ... by ... to クローズを使用せず、すべてを明示的な形で設定しています。それぞれの時点におけるElevationプロット、Contourプロットが計18枚出力されるわけですが、ここではその一部を示すに留めます。

(1) Elevation(temp) from (0, Ly/2) to (Lx, Ly/2) at t = 100
2本の鉄棒を接触させてから100秒経過した時点における鉄棒軸上での温度分布をプロットしたものです。鉄棒の左側の部分と右側の部分とでは200度の温度差があります。

(2) Contour(temp) painted at t = 100
t = 100 の時点における温度分布を等温線図の形でプロットしたものです。

(3) Elevation(temp) from (0, Ly/2) to (Lx, Ly/2) at t = 100000
t = 100000 における温度分布のグラフです。十分な時間の経過後であるため、一様な温度分布となっています。

(4) Elevation(temp) from (0, Ly/2) to (Lx, Ly/2) at t = 30000
(3)より一つ前の時点である t = 30000 における温度分布を示すグラフです。Movie再生すると最後のフレームである(3)とその一つ前のフレーム(4)の間で、グラフ形状に突然の変化が生じるように見えますが、表示されている temp のスケールに注意してください。例えば下図の場合、299.62から299.91の範囲だけが拡大されて表示されているためS字状の曲線に見えますが、これを(1)と同じ200から400のスケールでプロットすれば全くの直線に見えるはずです。

なお、プロット下部に表示されている温度の積分値（この場合は線積分）が t の値によらずほぼ一定値を保っている点に注意してください。これはエネルギー保存則に基づく帰結と言えます。