Jeremy Walton
The Numerical Algorithms Group, Ltd.
Wilkinson House, Jordan Hill Road
Oxford OX2 8DR, United Kingdom
目次
1 概要
2 偏微分方程式 － 技術的背景
3 例
4 他の手法
5 結論
6 謝辞
7 参考文献
1 概要
偏微分方程式（PDE）は複数の変数の関数と偏導関数を含む数学的関係式です。多変数関数を含む問題を公式化するために（また問題の解法に役立つために）偏微分方程式は使用されています。これらの偏微分方程式は様々な分野で現れます。例えば物理学では、音や熱の伝播、静電気学、電気力学、流体や弾力性を説明するのに使用されています。さらに金融においても金融オプションのプライシングのモデリングに使用されてきました。その結果、それらの性質や手法についての研究は非常に注目を集めてきました。
[1]
このレポートでは偏微分方程式の解法へのnAG ライブラリ[2] の応用について述べています。分類（セクション2.1）、付帯条件（セクション2.2）、解法（セクション2.3）を含む技術的な詳細を説明し、続いてセクション3では偏微分方程式の例とその解法をいくつか説明しています。これらは方程式の種類によって分類されます；解を得るために使用される適切な nAG ルーチンの詳細についても説明されています。例題の結果は nAG Toolbox for MATLAB® [3] のルーチンを使用して生成され、その環境のツールを用いてプロットされます。
セクション3で述べられている問題は有限差分数値法（例えば、nAGライブラリの偏微分方程式のチャプター[4]で実装されているような）を使用して処理できます。これは偏微分方程式が定義される領域の形状の特性によります。他の種類の問題（例えば、不規則な形状の領域をもつ問題）はいわゆる有限要素法の応用が必要です。セクション4では、この有限要素法の２つの構成要素、すなわち問題領域に渡るメッシュの生成（セクション4.1）と偏微分方程式の連立方程式への転換（セクション4.2）について説明しています。またnAGライブラリのルーチン[5] [6] [7]がそれらの実装[8]にいかに役立っているかを示しています。最後のセクション（セクション5）では結論が述べられています。
2 偏微分方程式 − 技術的背景
2.1 分類
偏微分方程式の例は２次元のラプラス方程式です：
 |
(1) |
ここで u(x,y) は２次元空間の位置関数です。この方程式は重力、静電気、流体力学といった領域のポテンシャル場の動きを表します。これは以下のようなより簡潔な表記で書くことができます:
| uxx + uyy = 0 | (2) |
あるいはもっと簡潔に∇2u=0 と書くことができます。ラプラスの方程式はいわゆる平衡偏微分方程式一つの例です。なぜなら時間に依存しない関数を説明しているためです。したがって不変の物理システムをモデル化するために使用できます。（例えば定常状態の任意の範囲の熱分布を決定します。）それに対し、一次元の波動方程式のような偏微分方程式は時間依存のプロセスをモデル化する動的偏微分方程式です：
| utt = ν2uxx | (3) |
これは波の速度に対応した定数 をもつ、時間 t と位置 x の関数である変位u(t,x)の動きを表しています。これは引き伸ばされた弦の平衡状態からの変位あるいは管の電界の大きさを説明するのに使用されます。最終的に偏微分方程式の次数はその式が含む最上位の導関数の次数です。(2) と (3) 両方とも二階微分の偏微分方程式です。科学的応用で現れる偏微分方程式の多くはこの種類です。線形二階偏微分方程式の一般式は以下のように表すことができます。
| αuxx + 2βuxy + γuyy + δux + εuy + ζu + η = 0 | (4) |
ここで係数α、β、γ、δ、ε、ζ、η は２つの独立変数 x と y のみの関数です。（もちろん、動的偏微分方程式ではそれらのうちの一つが t です。）もしも判別式αγ−β2 が正の場合は、偏微分方程式は「楕円型」と呼ばれます。負の場合は「双曲型」と呼ばれます。0 の場合は、「放物型」と呼ばれます。この定義によると、α(x,y)=γ(x,y)=1、β(x,y)=0 となるラプラス方程式(2)は楕円型であるとわかります。 α(t,x)=1、γ(t,x)=−1、β(t,x)=0 となる場合の波動方程式(3)は双曲型です。放物型偏微分方程式の例は以下の一次元の熱伝導方程式で表されます。
| ut = κuxx | (5) |
この式は任意の範囲の時間と位置の関数であるu(t,x) で表される発熱（あるいは同じく温度）を決定します。ここで は熱拡散率に関連する定数です。この偏微分方程式の変形が金融オプションプライシングのブラックショールズ偏微分方程式の解で現れる点にご注意ください。
2.2 付帯条件
一般に任意の偏微分方程式を満たす多くの関数があります。例えば(5)の場合、これらのうちの２つは以下の式です。
 |
(6) |
及び
 |
(7) |
どの解法が適切であるかを判断するにはより多くの情報が必要です。それらには、u が定義される対象領域（領域の境界を示す曲線 C の仕様を含む）D の詳細と、u とその導関数の付帯条件が含まれます。後者は領域の境界上の解を特定する境界条件、及び／または開始時の解を表す初期条件となって現れます。
厳密に言えば、(2)のような平衡方程式の場合、付帯条件により C の全ての点の解についての情報が得られます。そのため解は以下の条件で求められます。
| u(x,y) = f(x,y) ∀(x,y)∈C | (8) |
ここでf(x,y) は任意の関数です；これはディリクレ境界条件として知られています。C の解の導関数の値を特定する条件は以下です。
| ∇u(x,y)•n(x,y) = g(x,y) ∀(x,y)∈C | (9) |
この条件はノイマン境界条件と呼ばれます。ここでg(x,y) は任意の関数で、∇は勾配ベクトル演算子であり、ドットは内積を表し、n は C の法線です。
この場合、唯一の解を決定するのに追加情報（特定の位置での解の値など）が必要となります。最終的には、ロビン境界条件は他の二つ境界条件を重みづけして組み合わせたものと言えます：
| ∇u•n + uf1 = f2 ∀(x,y)∈C | (10) |
ここでf1(x,y) とf2(x,y) は任意の関数です。
(3)のような動的方程式の場合、開始時点で解についての情報が付帯条件によって得られます。そのため、例えば(3) が のときに満たされる場合、解が以下の条件で求められます。
| u(0,x) = f(x), ut(0,x) = g(x) | (11) |
ここでf(x) とg(x) は任意の関数です。このように公式化された問題は初期値問題（コーシー問題として呼ばれる場合があります）として知られています。付帯条件のその他の例では空間制約と時間的制約の両方を指定している場合があります。(5)がt>0 と0<x<1 の場合に満たされていると仮定すると、以下を満たす解が求められます。
| u(0,x) = f(x), 0 < x < 1, | (12) |
及び
|
u(0,x) = 0 ut(t,1) = 1, t > 0 |
(13) |
これは混合型初期値／境界値問題の例です。
2.3 解法
偏微分方程式を解くために様々な手法が開発されました [1] [9]。
有限差分や有限要素は最も利便性のある数値法[11]の一部ですが、解析手法には[10] 変数分離法、特性曲線法、積分変換（例えばフーリエ分析）が含まれます。(2)のような線形平衡偏微分方程式の場合、その応用は線形システムソルバーを用いて処理される連立方程式をもたらします。動的偏微分方程式は（連続したままの）空間座標を離散化するためにたびたびこれらの手法を用いて対処されます。それらは適切な技術（セクション3.2、3.3及び3.4以降）で解決することができる常微分方程式（ODE）を生成します。特種な問題では別の解法を利用いただくことができます。
3 例
このセクションでは偏微分方程式の問題のいくつかの例とnAGライブラリのd03 チャプター [4] のルーチンを用いた解法について説明しています。
3.1 ３次元楕円偏微分方程式
３次元のラプラス方程式の解を求めます：
| uxx + uyy + uzz = 0 | (14) |
各方向に一様でない格子間隔をもつ長方形を形成し、以下のデリクレット条件をもちます。
 |
(15) |
ここで Y は y 方向のボックスの長さを示します。このような問題は有限差分法により処理されます。
つまり、これは D 上でのメッシュの生成をもたらします。またメッシュポイントの離散集合の値による関数の表現をもたらします；偏微分方程式は隣接点の値の間の差分について近似することができます。
これは偏微分方程式をnAGルーチンd03ec [12]で解くことができる連立代数方程式に変換します。これは方程式の解に陰的手法を積極的に採用しています。図１は等密度面（関数が特定の値をもつボックス内の全ての点と接続する３次元空間の曲線）として表示される解を表しています。また、曲面プロットとして表示される解を表しています。その曲面プロットはボックスから切りとられる２次元スライスによってどのように関数が変化するかを表します。図２は固定された z値でボックスの y 方向に沿って関数がどのように変化するか、また一連の x の値に対してどのように変化するかを表しています。
3.2 移流拡散方程式
一次元のいわゆる移流拡散方程式の最も簡易な式 は以下で表されます。
| ut = Cux + Duxx | (16) |
これは種濃度u(t,x)の動的進化が移流項と拡散項の総和によって与えられることを示しています［この偏微分方程式を拡散方程式(5)と比較します］。拡散が濃度勾配に従って種の広がりを示す一方、下層の流体の大量の運動によって移流が種の移送となっていることがあらためてわかります。
セクション2.3で述べたように、このような偏微分方程式は空間次元（時間ではなく）によって解かれます。それは偏微分方程式を硬いソルバーによって解が得られる常微分方程式に変換します。実際にnAGルーチンは特に以下のような、(16)の式よりも一般的な式を解きます。
 |
(17) |
ここで n は偏微分方程式の数とui(t,x) の解の数を表します。Diは（特に）∂ui ⁄∂x の関数です。そのため右辺の最初の項は∂2ui ⁄∂x2 に比例する項（例えば拡散項）を含んでいます。同様にFiはuiに依存し、そのため左辺の２番目の項は移流項となり、一方Si はソース項となります。例題は（単一成分からなる）無数の物質の拡散と移流のシンプルなモデルです：
| ut + φux = δuxx | (18) |
x∈[0,1] の場合、境界条件は以下です。
| u(t,0) = u(t,1) = 0 | (19) |
また初期条件は以下です。
 |
(20) |
この場合、a=0.2、b=0.4 であり、それ以外では u(0,x)=0 で、φとδは物理パラメータです。より一般的な式である(17) と(18)を比較することでφuを持つF とδuxを持つ D を特定します。
上記で述べましたように、このような問題は偏微分方程式を常微分方程式に変換することによって（特にこの例ではいわゆる線の方法を用いて）解くことができます。その方程式は例えば後退微分公式（BDF）を用 いて解くことができます。nAGルーチン d03ps[13] はこの技術を使用しています。また空間格子の自動の適応的再メッシュを組み込んでいます。図３は d03ps によって計算された(18)の解を示しています；時間が増加するにつれてx の値が大きい方向に濃度ピークが動いているのがはっきりと表れています。
3.3 放物型偏微分方程式
一空間次元の一般的な放物型偏微分方程式は以下のように記述されます。
 |
(21) |
ここで n は偏微分方程式の数であり、項Pij、Qi とRi は一般的に t、x、u とux の関数です。
パラメータ m は異なる座標系の処理、特にカルテシアン座標（m=0）、円柱座標（m=1）、球座標（m=2）の処理を可能にします。
このような方程式は空間微分係数を離散化することによって解くことができます。また硬いソルバー（例えば後退微分公式）を用いて結果として生じた常微分方程式を解くことにより解を求めることができます。離散化に使用できる様々な手法があります。例えば nAG ルーチンd03pd はチェビシェフ配置法を使用しますが、d03pc は有限差分を使用します。
例 [14]は 楕円放物型の偏微分方程式です：
 |
(22) |
t、r ∈ [0,1] に対して、境界条件は以下です。
 |
(23) |
また初期条件は以下です。
| u1(0,r) = 2ar, u2(0,r) = 1, 0 ≤ r ≤ 1 | (24) |
(22)と(21)とを比較すると以下がわかります。
|
P22 = (1−r2), それ以外は0 Q1 = 4a(u2 + r∂u2 ⁄ ∂r), Q2 = 0 R1 = r∂u1 ⁄ ∂r, R2 = ∂u2 ⁄ ∂r − u1u2 |
(25) |
図４はd03pcで計算されたこの偏微分方程式の解を表しています。 このルーチンはまた 以下の解で[16]使用されています：
 |
(26) |
ここで D、h、g は u （モデルシステムの密度）の関数です。３つの関数は拡散、移流、反応 に対応しています； (21)と比較すると以下がわかります。
| R = Dux − h, Q = −g | (27) |
3.4 ２次元の時間依存型偏微分方程式
いくつかの偏微分方程式は２次元空間をもつ時間依存型として分類することができます。一般的 にそれらは以下のように記述することができます。
| Fj(t,x,y,u,ut,ux,uy,uxx,uxy,uyy) = 0 j = 1,2,...,n | (28) |
ここで n は偏微分方程式の数であり、uj(t,x,y) の解の数です。私たちは矩形の空間領域 D 上で定義された方程式に関心があります：
| xmin ≤ x ≤ xmax; ymin ≤ y ≤ ymax | (29) |
このような問題の初期処理では D 上でメッシュをかけることによる有限差分法が（前述のように）使用されます。これは偏微分方程式を、可変のステップサイズで陰的後退微分公式を使用して時間積分される常微分方程式に変換します。 時間積分から生じる非線形方程式は修正ニュートン法を用いて解かれます。nAG ルーチンd03ra [17]はこれらの技術を使用し、解が急速に変化をしている領域で解の精度を改善するために局所的一様格子の改良を具体化します。最初の例題は空間の２次元領域[17]の２つの化学製品の混合物の単一の一段階反応のモデルです。
混合物の動的局部温度 u(t,x,y) に対する偏微分方程式は以下のように表されます。
 |
(30) |
ここで d、R、α、δは反応メカニズムと関連した物理パラメータです。領域は以下のように定義されます。
| xmin = ymin = 0; xmax = ymax = 1 | (31) |
また初期条件は領域全体を通して u(0,x,y)=1 です。境界条件は以下のとおりです。
|
ux(t,xmin,y) = 0,
u(t,xmax,y) = 1;
ymin ≤ y ≤ ymax uy(t,x,ymin) = 0, u(t,x,ymax) = 1; xmin ≤ x ≤ xmax |
(32) |
図５は d03ra で計算された (30)の解を表しています。t=0.4 付近で発火が起こる前に、気 温が突然 2.0 に跳ね上がり、反応面が生じ外側に伝播するときに、原点のまわりの円領域（全 領域の四分円をモデリングしていると仮定）で気温が徐々に上昇していることがわかります。 これは以下の式の解にも適用されています [16]。
| ut = d(uxx + uyy) + f(u) | (33) |
ここで f は開区間(0,1)上で f(0)=f(1)=0 かつ f>0 となる、区間 [1,0] 上の凹関数です。
4 他の手法
上記で説明のあった d03 ルーチンは有限差分の技術を使用しており、規則的な形状をもつ領域 で解が求められる、限定された数の特定の偏微分方程式に適用されます。 不規則な形状を扱う 手法の一つにいわゆる有限要素法（FEM）があります。この手法は空間次元で方程式を最初に離散化することによって偏微分方程式（あるいは他の変分問題）を解く技術[18] です。これには２つの構成要素があります：一つ目は空間領域のメッシュへの分割であり、２つ目は偏微分方程式の連立方程式への変換です。ここでそれぞれ別々に考えてみます。
4.1 メッシュの生成
偏微分方程式が解かれる領域の分割は有限要素法の重要な要件です。多くの場合、解の精度と妥当性は潜在するメッシュの特性と強く結びついています。 領域に渡ってメッシュを生成する際には（通常）二段階のプロセス[5]があります。一つ目は、領域の境界をメッシュすることです。二つ目は、境界のメッシュ点を領域内のメッシュの生成の開始点として使用することです。
nAG ルーチン d06ba [19]は境界メッシュを生成するために使用することができます；境界は任意の形状であり、任意の領域（領域は内部境界によって特定される穴を含む可能性があります）に対して二つ以上の境界がある点に注意が必要です。様々な異なる手法が領域内のメッシュの生成に利用できます。これらのうちの３つが、増分解法（d06aa [20]で使用）、Delaunay-Voronoi 法 (d06ab [21] ) と前進フロント法 (d06ac [22] )です。
最後に、d03ma ルーチン[23]は１回の呼び出しで２次元の領域に渡り三角形メッシュを生成する点にご注意ください。
4.2 偏微分方程式の変換
領域をメッシュ化する場合、有限要素法の２番目のステップでは偏微分方程式を連立方程式へ変換する際にメッシュ（通常三角形）の要素を使用します。これは各要素内の方程式を表すために最初に一連の基底関数を選択することによって行われます。この基底関数は要素の各ノードの入力を同じ点の出力と関連づける行列方程式をもたらします。領域全体に対する連立方程式全体の係数の数はメッシュの大きさ－通常ほぼ何千程度－に比例します。しかし、対角要素と少数の非対角要素だけが非ゼロとなるため、係数行列は帯行列となります。
nAG ライブラリには、方程式の構造を活用するこの種の方程式を解くためのルーチンが含まれています。（f07 チャプター [6])
一方で、行列はスパースの場合もあるため、このような方程式を解くための特別な手法（例えば、非ゼロ要素を格納するだけでスパース性を利用する手法）を使用することができます。これは私たちがここで行っている手法であり、ライブラリのf11 チャプター [7] のルーチンを使用します。例えば f11jc [24] は、共役勾配法あるいはランチョス法を用いて実スパース対称一次方程式を解くための５つの他の f11 ルーチンを呼び出すいわゆるブラックボックスルーチンです。ルーチン f11ja [25] は行列の不完全コレスキー分解を計算するためにf11jc を使用する前に呼び出される必要があります（この前処理ステップはソルバーの速度と効率を増やします）；さらに、f11ja の前にf11zb [26]を呼び出すことは行列の非ゼロ要素の並べ替えによってソルバーの効率性をさらに増加させます。
d06 と f11 ルーチンは、任意の形状のユーザ指定の２次元領域に渡って解くための有限要素法(1)（そして非ゼロ右辺をもつ変数、すなわちポアソン方程式）を使用するデモ用偏微分方程式ソルバー[8]で使用されています。ソルバーは nAG Toolbox for MATLAB の一部として配布されています。
最後に、偏微分方程式を解くために nAG のメッシュ生成とスパースソルバールーチンを使用する技術的背景については[27]を参照下さい。
5 結論
このテクニカルレポートでは、偏微分方程式の使用、分類、付帯条件や解くためのいくつかの手法を含め、偏微分方程式の特徴について述べてきました。この中でnAGライブラリのルーチンが数値解法でどのように使用されるかを示してきました。これらのルーチンはライブラリの偏微分方程式チャプター (d03)だけでなく、メッシュ生成(d06) や帯状大規模連立一次方程式の解(f07)もしくはスパース大規模連立一次方程式の解(f11)を扱うチャプターのルーチンです；これらの３つのチャプターは有限要素法の実装に適用することができます。それらは偏微分方程式が解かれる領域の形状の複雑さが有限要素法の応用を妨げる場合に使用されます（例えば、d03 チャプター）。
最後に、ここで述べられた内容に加えてnAGライブラリ[2]が多種多様な数値問題や統計問題（例えば最適化、曲線及び曲面フィッティング、求積法、相関、回帰分析や乱数生成など）を扱うユーザ呼び出し可能なルーチンを含んでいることが注目されるでしょう。
6 謝辞
私たちは、いくつかの例をご提供いただいた Judith Perez-Velazquez 教授 (Institute of Biomathematics and Biometry, Helmholtz Zentrum, Munich) 、このレポートを作成するまでに有益な議論をいただいたLawrence Mulholland 教授とDavid Sayers 教授(nAG) 、さらに洞察力のあるコメントをいただいた David Cassell と Mike Hooper (nAG) に感謝申し上げます。
7 参考文献
| [1] | Fox, L, [ed.]. Numerical Solution of Ordinary and Partial Differential Equations. Oxford : Pergamon Press, 1962. |
| [2] | Numerical Algorithms Group. nAG 数値計算ライブラリコンポ―ネント。 [オンライン] http://www.nag.com/numeric/numerical_libraries.asp. |
| [3] | -. The nAG Toolbox for MATLAB. [オンライン] http://www.nag.com/numeric/MB/start.asp. |
| [4] | -. D03 - 偏微分方程式。nAG チャプター・イントロダクション。 [オンライン] http://www.nag.com/numeric/fl/nagdoc_fl23/pdf/D03/d03intro.pdf. |
| [5] | -. D06 - メッシュ生成。nAG チャプター・イントロダクション。 [オンライン] http://www.nag.com/numeric/fl/nagdoc_fl23/pdf/D06/d06intro.pdf. |
| [6] | -. F07 - 一次方程式。nAG チャプター・イントロダクション。 [オンライン] http://www.nag.com/numeric/fl/nagdoc_fl23/pdf/F07/f07intro.pdf. |
| [7] | -. F11 - 大規模線形方程式。 nAG チャプター・イントロダクション。 [オンライン] http://www.nag.com/numeric/fl/nagdoc_fl23/pdf/F11/f11intro.pdf. |
| [8] | Esteves, Nicolas, Fenton, Nathaniel and Walton, Jeremy. Using the nAG
Toolbox for MATLAB - Part 4. [オンライン] November 2009. http://www.nag.com/IndustryArticles/usingtoolboxmatlabpart4.asp. |
| [9] | Press, William H., et al. Numerical Recipes: The Art of Scientific Computing. Cambridge : Cambridge University Press, 1986. |
| [10] | Wikipedia. Analytical methods to solve PDEs. [オンライン] http://en.wikipedia.org/wiki/Partial_differential_equation#Analytical_methods_to_solve_PDEs. |
| [11] | -. Numerical methods to solve PDEs. [オンライン] http://en.wikipedia.org/wiki/Partial_differential_equation#Numerical_methods_to_solve_PDEs. |
| [12] | Numerical Algorithms Group. D03ECF nAG ルーチンドキュメント。 [オンライン] http://www.nag.com/numeric/fl/nagdoc_fl23/pdf/D03/d03ecf.pdf. |
| [13] | -. D03PSF nAG ルーチンドキュメント。 [オンライン] http://www.nag.com/numeric/fl/nagdoc_fl23/pdf/D03/d03psf.pdf. |
| [14] | -. D03PCF nAG ルーチンドキュメント。[オンライン] http://www.nag.com/numeric/fl/nagdoc_fl23/pdf/D03/d03pcf.pdf. |
| [15] | -. D03PDF nAG ルーチンドキュメント。 [オンライン] http://www.nag.com/numeric/fl/nagdoc_fl23/pdf/D03/d03pdf.pdf. |
| [16] | Perez-Velazquez, Judith. 2011. Personal communication. |
| [17] | Numerical Algorithms Group. D03RAF nAG ルーチンドキュメント。 [オンライン] http://www.nag.com/numeric/fl/nagdoc_fl23/pdf/D03/d03raf.pdf. |
| [18] | Strang, Glibert and Fix, George J. An Analysis Of The Finite Element Method. Englewood Cliffs, NJ : Prentice-Hall, 1973. |
| [19] | Numerical Algorithms Group. D06BAF nAG ルーチンドキュメント。 [オンライン] http://www.nag.com/numeric/fl/nagdoc_fl23/pdf/D06/d06baf.pdf. |
| [20] | -. D06AAF nAG ルーチンドキュメント。 [オンライン] http://www.nag.com/numeric/fl/nagdoc_fl23/pdf/D06/d06aaf.pdf. |
| [21] | -. D06ABF nAG ルーチンドキュメント。 [オンライン] http://www.nag.com/numeric/fl/nagdoc_fl23/pdf/D06/d06abf.pdf. |
| [22] | -. D06ACF nAG ルーチンドキュメント。 [オンライン] http://www.nag.com/numeric/fl/nagdoc_fl23/pdf/D06/d06acf.pdf. |
| [23] | -. D03MAF nAG ルーチンドキュメント。 [オンライン] http://www.nag.com/numeric/fl/nagdoc_fl23/pdf/D03/d03maf.pdf. |
| [24] | -. F11JCF nAG ルーチンドキュメント。[オンライン] http://www.nag.com/numeric/fl/nagdoc_fl23/pdf/F11/f11jcf.pdf. |
| [25] | -. F11JAF nAG ルーチンドキュメント。 [オンライン] http://www.nag.com/numeric/fl/nagdoc_fl23/pdf/F11/f11jaf.pdf. |
| [26] | -. F11ZBF nAG ルーチンドキュメント。 [オンライン] http://www.nag.com/numeric/fl/nagdoc_fl23/pdf/F11/f11zbf.pdf. |
| [27] | Bouhamou, Nadir. The use of nAG mesh generation and sparse solver routines
for solving partial differential equations. [オンライン] 2001. http://www.nag.com/doc/techrep/pdf/Tr1_01.pdf. |