拡散方程式は、物質の濃度や温度などが時間とともに空間へ広がっていく現象を表す偏微分方程式です。
熱伝導、化学物質の移動、半導体内部のキャリア輸送、生体組織への薬剤浸透など、幅広い分野で利用されています。
しかし、常微分方程式とは異なり、時間と位置という複数の変数を含むため、初めて学ぶ方には解き方が複雑に見えるかもしれません。
基本的な境界値問題では、変数分離法によって時間と空間を分け、固有値問題を解いた後にフーリエ級数を用いて初期条件を満たすという流れが中心になります。
一方、無限領域や外力を含む問題ではグリーン関数が便利であり、解析解を得にくい複雑な形状では差分法や有限要素法などの数値解法が必要です。
本記事では、拡散方程式の基本形から変数分離法の具体的な計算手順、境界条件の扱い方、フーリエ級数、グリーン関数、数値計算まで順を追って解説します。
拡散方程式は変数分離法と初期条件の展開によって解けます
それではまず、拡散方程式の基本的な解き方について解説していきます。
拡散方程式の標準形を確認する
一次元の拡散現象を表す標準的な方程式は、濃度や温度を表す関数を位置と時間の関数として記述します。
代表的な形は次のとおりです。
∂u／∂t＝D∂²u／∂x²
ここでuは位置xと時間tに依存する物理量であり、物質拡散では濃度、熱伝導では温度に相当します。
Dは拡散係数で、物質や媒質によって決まる正の定数です。
拡散係数が大きいほど空間的な偏りが速く解消され、分布が短時間でなだらかになります。
右辺にある二階の空間微分は、分布の曲がり具合を示す量です。
ある位置の値が周囲より高い場合は分布が上向きに膨らみ、その値が時間とともに低下します。
反対に周囲より低い位置では、近くから物質や熱が流れ込むため値が上昇する仕組みです。
拡散方程式を解くとは、方程式だけを変形することではなく、初期条件と境界条件を同時に満たす関数uを求めることです。
同じ偏微分方程式でも、端点の条件や初めの分布が異なれば得られる解も変わります。
したがって、問題文では方程式に加えて初期条件と境界条件を丁寧に確認する必要があります。
変数分離法で空間と時間を分ける
変数分離法では、解uを位置だけの関数Xと時間だけの関数Tの積として仮定します。
u（x，t）＝X（x）T（t）
この形を拡散方程式へ代入すると、左辺はXとTの一階時間微分の積になり、右辺はTとXの二階空間微分の積になります。
X（x）dT／dt＝DT（t）d²X／dx²
1／D・1／T・dT／dt＝1／X・d²X／dx²
左辺は時間tだけに依存し、右辺は位置xだけに依存しています。
異なる変数だけで構成された両辺がすべてのxとtに対して等しくなるには、両方が同じ定数でなければなりません。
そこで分離定数をマイナスλと置きます。
1／D・1／T・dT／dt＝1／X・d²X／dx²＝－λ
マイナスを付ける理由は、両端の値をゼロに固定する境界条件では、正弦関数による減衰解が得られるためです。
これにより、一つの偏微分方程式を二つの常微分方程式へ分解できます。
d²X／dx²＋λX＝0
dT／dt＋DλT＝0
空間側の方程式は境界条件から許されるλを決める固有値問題です。
時間側の方程式は指数関数で解けるため、空間側よりも比較的簡単に処理できます。
固有関数を重ね合わせて初期条件を満たす
区間0以上L以下で両端の値が常にゼロとなる境界条件を考えます。
u（0，t）＝0
u（L，t）＝0
変数分離形へ代入すると、空間関数はX（0）＝0とX（L）＝0を満たさなければなりません。
空間方程式の一般解は正弦関数と余弦関数の組み合わせです。
X（x）＝A cos（√λx）＋B sin（√λx）
X（0）＝0からA＝0となります。
さらにX（L）＝0を満たすためには、sin（√λL）＝0が必要です。
したがって、√λLはπの整数倍でなければなりません。
λn＝（nπ／L）²
Xn（x）＝sin（nπx／L）
時間側の解は、それぞれの固有値に対応して指数関数的に減衰します。
Tn（t）＝exp｛－D（nπ／L）²t｝
一つのnに対応する解だけでは一般的な初期分布を表せないため、すべての固有関数を重ね合わせます。
u（x，t）＝ΣBn sin（nπx／L）exp｛－D（nπ／L）²t｝
係数Bnは初期条件u（x，0）＝f（x）をフーリエ正弦級数へ展開して決定します。
これが、有限区間における拡散方程式の代表的な解析解です。
境界条件を整理して固有値問題を解く手順
続いては、境界条件の種類と固有値問題の解き方を確認していきます。
ディリクレ境界条件では端の値を指定する
ディリクレ境界条件とは、領域の端における関数uの値を指定する条件です。
熱伝導問題では端点の温度を一定に保つ状況、物質拡散では境界の濃度を一定にする状況に対応します。
両端をゼロにする同次境界条件は、最も基本的な例です。
| 境界条件 | 物理的な意味 | 代表的な固有関数 |
|---|---|---|
| u（0，t）＝u（L，t）＝0 | 両端の値をゼロに固定 | 正弦関数 |
| u（0，t）＝A | 左端の値を一定に固定 | 定常解を引いた後に正弦関数 |
| u（0，t）＝A、u（L，t）＝B | 両端を異なる一定値に固定 | 線形定常解と過渡解 |
境界値がゼロでない場合、そのまま変数分離法を適用すると計算しにくくなります。
この場合は、境界条件を満たす時間に依存しない定常解を先に求める方法が有効です。
例えば左端がA、右端がBなら、定常解vは位置に対する一次関数になります。
v（x）＝A＋（B－A）x／L
元の解をu＝v＋wと置けば、wは両端でゼロとなります。
その結果、wについて通常のフーリエ正弦級数を用いた解法を適用できます。
非同次境界条件を同次化する操作は、計算ミスを減らすうえで重要です。
ノイマン境界条件では流束を指定する
ノイマン境界条件は、関数そのものではなく空間微分の値を指定する条件です。
熱伝導では熱流束、物質移動では物質流束の指定に対応します。
端から熱や物質が出入りしない断熱境界では、空間微分がゼロになります。
∂u／∂x（0，t）＝0
∂u／∂x（L，t）＝0
この条件を空間関数Xに適用すると、代表的な固有関数は余弦関数です。
Xn（x）＝cos（nπx／L）
ノイマン境界条件ではn＝0のモードも含まれます。
n＝0の固有関数は空間的に一定であり、時間が経過しても減衰しません。
これは、境界から物質が出入りしない閉じた系では全物質量が保存され、最終的な平均濃度がゼロになるとは限らないためです。
ディリクレ条件では正弦関数、ノイマン条件では余弦関数が現れやすいという対応を覚えると、解法の見通しが立てやすくなります。
ただし、片側がディリクレ条件で反対側がノイマン条件となる混合境界条件では、半整数倍の固有値が現れる場合があります。
条件を暗記するだけでなく、一般解へ実際に代入して固有値を求める姿勢が大切です。
ロビン境界条件では値と流束を組み合わせる
ロビン境界条件は、関数uとその空間微分を組み合わせた条件です。
物体表面から周囲の流体へ熱が伝わる対流熱伝達などで用いられます。
－k∂u／∂x＝h（u－u∞）
ここでkは熱伝導率、hは熱伝達率、u∞は周囲の温度です。
境界の値と流束が互いに関係するため、固有値を単純な整数倍として書けないことがあります。
三角関数を含む超越方程式が得られ、数値的に固有値を求める必要が生じるでしょう。
例えば、無次元化した空間関数に対してtanを含む固有値方程式が現れます。
固有値が求まれば、それぞれに対応する固有関数を構成し、直交性を利用して係数を決定する流れは同じです。
境界条件が複雑でも、空間固有値問題と時間減衰問題へ分けるという基本構造は変わりません。
フーリエ級数を使って初期分布を表す方法
続いては、初期条件からフーリエ係数を求める方法を確認していきます。
初期条件を固有関数の和として展開する
変数分離法で得られる一つ一つの固有関数は、特定の波形しか表せません。
ところが実際の初期分布は、三角形、矩形、区分的な関数など、さまざまな形を取ります。
そこで、固有関数を無限に重ね合わせることで初期分布を表現します。
両端がゼロとなるディリクレ境界条件では、初期関数fを正弦級数で展開します。
f（x）＝ΣBn sin（nπx／L）
正弦関数同士には直交性があります。
∫0からL sin（nπx／L）sin（mπx／L）dx＝0
ただしnとmは異なる整数です。
この性質を利用すると、それぞれの係数Bnだけを取り出せます。
Bn＝2／L∫0からL f（x）sin（nπx／L）dx
初期関数を積分へ代入し、nごとに係数を計算すればよいでしょう。
区分的に定義された関数では、積分区間を分割して計算します。
対称性がある場合は、偶数番目または奇数番目の係数がゼロになり、式を簡略化できることもあります。
具体例として一定の初期分布を解く
区間の内部で初期濃度が一定値C0であり、時刻tが正になった後は両端の濃度がゼロに保たれる問題を考えます。
初期条件は次のように表せます。
u（x，0）＝C0
ただし0より大きくLより小さいxです。
フーリエ係数は一定値C0を積分へ代入して求めます。
Bn＝2C0／L∫0からL sin（nπx／L）dx
Bn＝2C0／（nπ）｛1－（－1）n乗｝
nが偶数なら係数はゼロになります。
nが奇数なら係数は4C0／nπです。
したがって、解は奇数番目のモードだけを含みます。
u（x，t）＝4C0／πΣnが奇数 1／n sin（nπx／L）exp｛－D（nπ／L）²t｝
時刻tがゼロに近いときは、多数のモードが重なって一定に近い初期分布を形成します。
時間が経過すると高いnの項ほど速く減衰するため、細かな空間変化から先に消えていきます。
十分に時間が経つと最も減衰の遅いn＝1の項が支配的になるでしょう。
フーリエ級数解の物理的な意味を理解する
各固有関数は、拡散系に存在できる空間的なモードを表します。
nが大きいモードほど空間内で細かく振動し、二階微分の大きさも増加します。
その結果、指数関数に含まれる減衰率D（nπ／L）²が大きくなり、短時間で消滅します。
| モード | 空間的な特徴 | 減衰の速さ |
|---|---|---|
| n＝1 | 最も緩やかな分布 | 最も遅い |
| n＝2 | 区間内に細かな変化が増える | n＝1より速い |
| 大きなn | 短い波長の変動 | 非常に速い |
拡散によって分布が滑らかになるのは、高周波成分が急速に減衰するためです。
初期分布に不連続点があっても、正の時間が経過すると解は滑らかな形へ変化します。
この性質は熱方程式の平滑化効果とも呼ばれます。
フーリエ級数は単なる計算手段ではなく、複雑な初期分布を減衰速度の異なる空間モードへ分解する考え方です。
どのモードが長時間残るのかを確認すれば、厳密な無限級数をすべて計算しなくても挙動を予測できます。
グリーン関数を用いて無限領域や外力付き問題を解く方法
続いては、グリーン関数による拡散方程式の解法を確認していきます。
無限領域の基本解を理解する
無限に広い一次元空間で、一点に集中していた物質が拡散する問題を考えます。
初期条件をデルタ関数で表すと、その解はガウス分布になります。
G（x，t）＝1／√（4πDt）exp｛－x²／（4Dt）｝
この関数が拡散方程式の基本解であり、グリーン関数として利用されます。
時刻が小さいと分布は原点付近に鋭く集中しています。
時間が経過すると幅が広がり、最大値は低下する一方、全空間で積分した総量は一定に保たれます。
分布の広がりを表す標準偏差は時間の平方根に比例します。
代表的な拡散距離は√（2Dt）程度です。
つまり、拡散距離を二倍にするには時間を四倍にする必要があります。
一定速度で移動する運動とは異なり、拡散の広がり方は時間に比例しない点が特徴です。
一般的な初期分布との畳み込みを行う
初期分布がデルタ関数ではなく任意の関数fで与えられる場合、それぞれの位置に微小な点源が存在すると考えます。
各点源から広がる基本解を重ね合わせれば、全体の解を構成できます。
u（x，t）＝∫マイナス無限大から無限大 G（x－ξ，t）f（ξ）dξ
この積分はグリーン関数と初期分布の畳み込みです。
ξは初期時刻における点源の位置を表します。
位置ξにあった量が時刻tまでにxへ到達する影響をG（x－ξ，t）が示します。
初期分布が複数のピークを持つ場合も、それぞれに対するガウス分布を足し合わせることで解を得られます。
畳み込みの考え方は、フーリエ変換による解法とも密接に関係しています。
フーリエ空間では畳み込みが積へ変わるため、計算を簡略化できる場合があります。
発生項を含む非同次方程式へ応用する
内部で物質や熱が発生する場合、拡散方程式の右辺に発生項が加わります。
∂u／∂t＝D∂²u／∂x²＋q（x，t）
qは単位時間および単位体積当たりの発生量です。
この場合は、初期分布による影響と、各時刻に発生した源による影響を分けて考えます。
u（x，t）＝∫G（x－ξ，t）f（ξ）dξ＋∫0からt∫G（x－ξ，t－τ）q（ξ，τ）dξdτ
二つ目の積分では、過去の時刻τに位置ξで生じた量が、残り時間t－τだけ拡散した影響を合計しています。
この考え方はデュアメルの原理とも関係します。
時間的に変化する加熱、連続的な薬剤供給、汚染物質の放出などを解析する際に有効です。
グリーン関数は、一点かつ一瞬の入力に対する応答を基本単位として、複雑な初期条件や発生項の影響を重ね合わせる方法です。
境界が存在する場合には、鏡像法を使ったり、境界条件に対応したグリーン関数を構成したりします。
有限領域では固有関数展開によってグリーン関数を表すことも可能です。
差分法などの数値解法で拡散方程式を計算する方法
続いては、解析解が得にくい場合に用いる数値解法を確認していきます。
陽解法では現在の値から次の時刻を計算する
差分法では、空間と時間を細かな格子へ分割し、微分を近くの格子点の値で近似します。
空間刻みをΔx、時間刻みをΔtとします。
時間微分を前進差分、空間二階微分を中心差分で近似すると、次の更新式が得られます。
u i 上付きn＋1＝u i 上付きn＋r（u i＋1 上付きn－2u i 上付きn＋u i－1 上付きn）
r＝DΔt／（Δx）²
現在時刻nの値が分かれば、次の時刻n＋1の値を直接計算できます。
連立方程式を解く必要がなく、実装が簡単な点が陽解法の利点です。
しかし、時間刻みを大きくしすぎると数値解が振動したり発散したりします。
一次元拡散方程式では、一般にrが二分の一以下となる条件が必要です。
陽解法ではDΔt／（Δx）²≦1／2という安定条件を必ず確認する必要があります。
空間格子を半分に細かくすると、安定性を保つための時間刻みはおよそ四分の一にしなければなりません。
高精度化のために格子を細かくすると計算回数が急増する点に注意が必要です。
陰解法では次の時刻の値を連立方程式で求める
陰解法では、空間二階微分を次の時刻n＋1の値で評価します。
－r u i－1 上付きn＋1＋（1＋2r）u i 上付きn＋1－r u i＋1 上付きn＋1＝u i 上付きn
次の時刻における複数の格子点が同じ式に含まれるため、連立一次方程式を解かなければなりません。
一次元問題では係数行列が三重対角行列になるため、トーマス法などを使えば効率よく計算できます。
陰解法は拡散方程式に対して無条件安定とされ、陽解法より大きな時間刻みを選びやすい方法です。
ただし、時間刻みを極端に大きくすると安定であっても精度が低下します。
安定性と正確性は別の条件であることを理解しておきましょう。
毎時刻で連立方程式を解く計算負荷と、大きな時間刻みを使用できる利点を比較して選択します。
クランク・ニコルソン法や有限要素法を使い分ける
クランク・ニコルソン法は、現在時刻と次の時刻における空間微分の平均を利用する方法です。
陽解法や単純な陰解法より時間方向の精度が高く、拡散問題で広く用いられています。
一方、複雑な形状や不均一な材料を扱う場合には有限要素法が適しています。
有限要素法では領域を小さな要素へ分割し、弱形式と呼ばれる積分形を使って近似解を構成します。
曲線境界や三次元構造にも対応しやすく、工学解析ソフトの多くで採用されています。
| 数値解法 | 主な長所 | 注意点 |
|---|---|---|
| 陽的差分法 | 式が単純で実装しやすい | 時間刻みに安定条件がある |
| 陰的差分法 | 大きな時間刻みでも安定しやすい | 連立方程式が必要 |
| クランク・ニコルソン法 | 時間方向の精度が高い | 急な変化で振動が出る場合がある |
| 有限要素法 | 複雑形状や物性分布に強い | 定式化とメッシュ作成が複雑 |
数値計算では、格子幅を変えたときに解が収束するかを確認する格子依存性の検証も重要です。
解析解が知られている簡単な問題でプログラムを検証した後、複雑な問題へ適用すると信頼性が高まります。
拡散方程式を解く際に注意したいポイント
続いては、計算過程で起こりやすい間違いと確認方法を見ていきます。
初期条件と境界条件の整合性を確認する
初期条件と境界条件が角の部分で一致しない問題も存在します。
例えば、初期時刻では区間全体の濃度がC0である一方、時刻が正になると両端をゼロにする条件です。
この場合、時刻ゼロかつ端点では条件同士が一致しません。
それでも弱い意味や極限として解を考えることはできますが、端点付近では初期直後に非常に急な変化が生じます。
フーリエ級数は不連続点で左右の極限値の平均へ収束するため、端点そのものの値に注意が必要です。
数値計算でも初期時刻直後に大きな勾配が生まれ、細かな格子が必要になることがあります。
次元解析で数式の誤りを見つける
拡散係数Dの単位は、長さの二乗を時間で割ったものです。
したがって、指数関数の中に現れるD（nπ／L）²tは無次元になります。
指数の中に単位が残る場合、Lの二乗や時間tを付け忘れている可能性があります。
代表的な拡散時間は次の尺度で見積もれます。
拡散時間はL²／D程度です。
長さを二倍にすると必要時間はおよそ四倍です。
計算結果がこの尺度と大きく異なる場合は、単位換算や係数を再確認するとよいでしょう。
センチメートルとメートル、秒と分を混在させる間違いも頻繁に起こります。
長時間挙動と保存則を使って答えを検算する
両端の値がゼロであれば、十分に時間が経過した後の解は通常ゼロへ近づきます。
一方、両端が断熱または無流束であれば、全量が保存され、最終的には初期分布の空間平均へ近づきます。
解析解や数値解がこれらの長時間挙動と一致するか確認すると、符号や境界条件の誤りを発見できます。
無流束境界で総量が時間とともに大きく減少する場合は、離散化や境界処理に問題があるかもしれません。
また、物質濃度や絶対温度を扱う問題では、本来負にならない量が負になっていないかも確認します。
数値解法によっては、格子が粗い場合や時間刻みが不適切な場合に負値や振動が生じる可能性があります。
まとめ
拡散方程式の基本的な解き方は、解を空間関数と時間関数の積に置く変数分離法から始まります。
空間側では境界条件を満たす固有値と固有関数を求め、時間側では各モードの指数関数的な減衰を計算します。
その後、初期分布をフーリエ級数で展開し、固有関数の係数を決定すれば有限区間の解を構成できます。
無限領域ではガウス型のグリーン関数が基本となり、初期分布や発生項との畳み込みによって一般的な解が得られます。
複雑な形状や変化する物性を含む場合は、陽解法、陰解法、クランク・ニコルソン法、有限要素法などの数値解法が有効です。
方程式、初期条件、境界条件、保存則、長時間挙動を一つずつ整理することが、拡散方程式を正確に解くための近道といえるでしょう。
