熱方程式のFTCS法、格子を細かくすると時間刻みが厳しくなる理由

熱の広がりを細かく計算しようと、位置の間隔を半分にする。ところが時間の刻みをそのままにすると、温度が一つおきに上下し、計算を進めるほど振れ幅が大きくなることがあります。熱が発生したのではなく、計算手順が小さな凹凸を増幅している状態です。

熱方程式を解くFTCS法は、今ある温度から次の時刻の温度を直接求める、比較的単純な方法です。その扱いやすさと引き換えに、位置と時間の刻みを別々には決められません。なぜ格子を細かくすると時間刻みも厳しくなるのかを、更新式と短い計算例から見ていきます。

次の温度を、今の三つの温度から求める

ここでは内部に熱源がなく、熱拡散率αが一定で正の値を持つ、一次元の熱方程式を考えます。温度をT、位置をx、時間をtとすると、式は∂T/∂t = α∂²T/∂x²です。αの単位はm²/sで、温度の凹凸がならされる速さに関係します。

FTCSはForward Time, Centered Spaceの略です。時間方向は今から一歩先への変化、空間方向は左右の点を使って微分を近似します。等間隔Δxで並んだ格子点のうち、i番目の温度をT_iとし、上付きのnを時間の段階とすると、更新式は次の形になります。

Tin+1 = r Ti−1n + (1−2r) Tin + r Ti+1n
r = αΔt/Δx²

右辺にあるのは、すべて今の時刻の温度です。未知の次時刻の温度をまとめて解く必要がないため、この方法は陽解法に分類されます。

rが0から1/2までなら、三つの重みr、1−2r、rはすべて0以上で、足すと1になります。次の温度は、今の自分と左右の温度の範囲に収まります。一方、rが1/2を超えると中央の重みが負になり、この「三点の平均として更新する」という性質が失われます。

交互に並ぶ小さな凹凸が、条件を教えてくれる

増幅の仕方を見やすくするため、8個の格子点を輪のようにつなぎ、右端の次を左端とする周期境界を使います。ある基準温度からの差を無次元化し、初めの値を+1、−1、+1、−1……と交互に置きます。実際の材料の測定値ではなく、数値計算の挙動を調べるための入力です。

この並びでは、ある点の値がuなら左右はどちらも−uです。更新式に入れると、次の値は(1−4r)uになります。r = 0.4では毎回−0.6倍、r = 0.5では−1倍、r = 0.6では−1.4倍です。どの場合も符号は交代しますが、振れ幅は順に小さくなる、変わらない、大きくなるという違いが出ます。

10回更新すると、振れ幅はそれぞれ約0.00605、1、28.9になります。r = 0.6のケースでは、ごく小さな誤差として混じった交互の模様でも、更新するたびに育ってしまいます。これが数値的不安定性の具体例です。

この結果は、特別な並びだけの偶然ではありません。周期境界または無限に続く一様格子で、凹凸を波の成分に分けて調べると、一回の更新での倍率はG = 1 − 4r sin²(kΔx/2)となります。kは空間方向の波数です。どの波数でも|G| ≤ 1を満たす条件が0 ≤ r ≤ 1/2で、格子点ごとに交互になる成分はsin²の部分が1となる端のケースに当たります。

格子幅を半分にしたら、時間刻みは4分の1へ

比rの分母にはΔxの二乗があります。同じ熱拡散率のまま格子幅だけを半分にすると、rは4倍になります。もともとr = 0.4で計算していても、時間刻みを変えなければr = 1.6となり、安定範囲を大きく超えます。同じrに戻すには、Δtを4分の1にする必要があります。

たとえば説明用にα = 1 × 10⁻⁵ m²/sと置くと、Δx = 1 mmでの境界はΔt = Δx²/(2α) = 0.05秒です。Δx = 0.5 mmなら0.0125秒になります。特定の材料を推奨する値ではなく、長さの単位をmにそろえて計算した例です。実際にr = 0.4を選ぶなら、それぞれ0.04秒と0.01秒になります。

細かな格子には、短い距離にある急な温度差も表せるという利点があります。ただし、その細かな凹凸を陽解法で追うには、時間方向も小刻みに進める必要があります。空間の解像度を上げたときに、計算時間まで増える理由の一つです。

安定範囲に入った後で、精度を確かめる

r = 0.5の計算では、交互の模様は発散しませんでした。しかし振れ幅が1のまま残るので、熱によって凹凸がならされる様子を適切に表しているとは言えません。安定条件は、誤差を際限なく増幅させないための条件です。それだけで、欲しい温度分布が正確に得られたことにはなりません。

計算を確認するときは、まず採用した方程式・格子・境界条件に合う安定条件を調べ、その範囲内で格子幅と時間刻みを小さくした結果を比較します。温度の最大値や注目する位置の変化が、細分化しても十分変わらなくなったかを見るためです。熱拡散率が場所で変わる場合、多次元の場合、別の境界処理や時間積分法を使う場合には、ここで示した一次元FTCS法の式をそのまま使うことはできません。

出典

• Gerry Recktenwald, FTCS for the 1D Heat Equation, in a Nutshell、Portland State University、Winter 2020、1〜2ページ。更新式、安定条件、刻みと誤差の関係。
• John S Butler, The Explicit Forward Time Centered Space (FTCS) Difference Equation for the Heat Equation、Numerical Analysis with Applications in Python。「Stability Analysis」の増幅因子。版日付の記載がない公開HTMLを2026年10月25日に確認。本文の8点周期格子と数値・図は、明示した更新式による独自計算。