Home Menu

いくつかの数値計算法

関数計算に用いた数値計算法のうち、比較的珍しい方法や独自の方法を中心に説明します※1。
なお、詳細やコード実装例は、別頁 Mathematica関連 にある個別のファイルを参照して下さい。
※1:このほか Euler - Maclaurin 総和法などを使用していますが、これらはより詳しい説明が他にも多数あるので省略します。なお、独自の方法に対して当サイトは種々の知的財産権を保有しています。注意事項 (Notes)!の内容に留意願います。

Padé 近似法

 H. Padé が1890年頃に開発・発展させた関数近似法であるが、その発想は合理的な級数近似法を模索した F. G. Frobenius の研究にまで遡れる。
 冪級数近似では、級数展開の中心点から最も近い所にある極などの特異点を超えた範囲では収束しない (有限の収束半径が存在する)。しかし、多項式ではなく有理関数で近似し、しかも両者が共に収束する有界な範囲で結果が一致していれば、有理関数近似は特異点を超えて収束すると期待できる。
 Padé 近似は、多項式と多項式の有理関数で上記の効果が実現されるような近似法であって、その定係数を具体的に求める際に正方行列を用いる。行列の形にはいくつかの種類がある。
 当サイトでは、極を持つ関数のみならず、級数の収束が遅いだけの関数等にも使用している。

計算法の説明

 Padé 近似を適用しようとする関数f(z)が、冪級数
冪級数展開式
に展開されるものとする。以降、これの2n項まで
冪級数展開式の有限項
の係数を使用する。二つの行列
  • 冪級数展開式の係数から作られる行列A,B
を作り、①線形代数方程式を解いて{b1, b2, …, bn}を求め、②それを用いて{a1, a2, …, an}を求める。
  • {b1, b2, …, bn}等を求める線形代数方程式
これらの結果を用いて、関数f(z)
  • Pade近似
で近似する。

簡単な計算比較の例

 正接関数 (タンジェント) を、Mathematica実装関数、冪級数近似 (50項)、Padé 近似 (50項/50項) で計算する。特異点よりも絶対値が小さい範囲内ならば、3種類とも同じ結果になる。
  • 数値計算結果の比較

 しかし、冪級数近似の収束範囲はAbs(z)<π/2を超えることができない。
  • 冪級数近似の場合の計算結果

 Padé 近似は、(冪級数と同じ50項使用でも) はるかに広い範囲で収束する。(しかも、約4秒しかかからない。これは、前のような非収束域の描画に時間を要しないため。)
  • Pade近似の場合の計算結果

適用例 ~ 第1種 Painlevé 超越関数

 近似しようとする関数が分岐点を持たず、収束半径が0でない冪級数に展開できて、その係数を具体的に求められる場合であれば、Padé 近似を適用できる。そのような例は非常に多いため、極めて有用な数値計算法といえる※1。
 例えば、z=μに2位の極がある第1種 Painlevé 超越関数は、その特異点を中心とする Laurent 級数
  • 第1種Painleve超越関数のLaurent級数展開
に展開できる。よって、有理関数項を除いた (冪級数の) 部分にPadé 近似を適用できる。
 次のグラフは、この方法によって第1種Painleve超越関数を計算した結果である。(行列の要素数は1500×1500。この事例では、絶対値が約15以上になると収束していない。)

  • 第1種Painleve超越関数のグラフ

Mathematica Code 正接関数および Painlevé 超越関数に Padé 近似を使用したコードの例を、Mathematica Code の頁に掲載しています。)

【註記】
※1 ただし欠点もある。この方法や後述の Fair - Luke 法,および三重対角行列法は、いずれも要素数の多い行列計算が伴うが、大抵の数値・数式処理システムでは負荷の掛かる計算になる。
 また、上記の描画事例に見られるような収束の限界をもっと遠方に広げる場合は、行列の要素数を大幅に増やす必要がある。(本来は漸近級数の併用が理想的だが、無数の特異点が散在する上記のような事例では、複素数値計算に適した漸近級数展開が確立されていないことも多い。)
* * *

Fair - Luke 法

 ある特定の形に限られるが、代数的非線形微分方程式の解に対してPadé 近似が適用できるよう、W. Fair と Y. L. Luke が開発した方法。近似に使う有理関数は、Padé 近似と同じく多項式の有理関数が選択できる他、連分数も選択できる。当サイトでは専ら後者を用いているので、それを説明する。
 また、私はこの方法が適用できる代数的非線形微分方程式の種類を若干広げた、「拡張された Fair - Luke 法」を開発したので、併せてこれも説明する。
 当サイトでは Painlevé 超越関数、および3階の Painlevé 超越関数に使用している。

「Fair - Luke 法」の説明

 2階非線形常微分方程式、
  • Fair-Luke法が使用可能な微分方程式
の解wが、連分数によって
  • Fair-Luke法の連分数
と表わされるとき、係数αnは漸化式
  • Fair-Luke法の漸化式
によって求められる係数を用いて
Fair-Luke法の連分数の係数
と定義される。なお、漸化式[3]は、微分方程式[1]において
Fair-Luke法を導出する変換式
なる変換を行うことで容易に確かめられる。
 ゆえに、[2]によって解wの数値計算が可能となる。これを「Fair - Luke 法」という。ここに条件
  • Fair-Luke法の連分数係数に課せられる条件
が課せられる。W. Fair と Y. L. Luke は論文※1において、第1種、第2種 Painlevé 方程式
  • 第1種・第2種Painleve方程式
は、この方法が適用できる方程式であることを指摘している。また V. Novokshenov は、実際にこのことを数値計算で検証し、複素領域における第1種、第2種 Painlevé 超越関数の動く特異点の詳細を明らかにした※2。具体的には、[7],[8]に
  • Fair-Luke法へ代入する関数形
を代入し、それぞれ
  • Fair-Luke法によって得られた初期値
に変形したのち、 Fair - Luke 法を適用する。明らかに、ω0, ω1
  • 第1種・第2種Painleve方程式の初期値
で定義される。

【註記】
※1 W. Fair and Y. L. Luke,
「Rational Approximations to the Solution of the Second Order Riccati Equation」
Mathematics of Computation, 20 (1968), pp.602-605.

※2 V. Yu. Novokshenov,
「Padé approximations for Painlevé I and II transcendents」
Theoretical and Mathematical Physics, 159 (No.3) (2009), pp.853-862.

「拡張された Fair - Luke 法」の説明

(2013年9月17日 掲載記事)
 微分方程式[1]を拡張したもので、その解が、連分数[2]によって表わされるような、 Fair - Luke 法の拡張を考える。
 始めに、第3種、および第4種 Painlevé 方程式
  • 第3種・第4種Painleve方程式
が取り扱えるような微分方程式の形を決定する。この場合、非線形な項で最高次のものはw^4で十分となり、かつ[5]の変換を施したときに生じる項が、変換前の微分方程式にあった項の種類内で収まるようなものとして
  • 拡張型Fair-Luke法が適用可能な微分方程式の候補
が候補となることが分かる。さらに検証によって、実際に条件を満たすためにはM=-2Cでなければならないことが判明する。よって
  • 拡張型Fair-Luke法が適用可能な微分方程式
が得られた。微分方程式[16]に変換[5]を施して得られる漸化式は
  • 拡張型Fair-Luke法の漸化式
となり、連分数[2]の係数はこの場合
  • 拡張型Fair-Luke法の漸化式に課せられる条件
と定義される。これを「拡張された Fair - Luke 法」と呼ぶこととする。
 第3種 Painlevé 方程式[13]の場合は、原点に動かない分岐点を持つため若干テクニックが必要となる。よってこれは後述し、先に第4種 Painlevé 方程式[14]の場合を述べる。[14]も[9]を代入すれば、拡張された Fair - Luke 法における漸化式の初期値として
  • 拡張型Fair-Luke法の漸化式の初期値
を採用できることが分かる。
 第4種の場合は、第1種、第2種の場合に比べて漸化式とその初期値が複雑であるため、計算可能な複素領域は狭くなる。
 第3種 Painlevé 方程式[13]の場合は、原点が動かない分岐点なので[12]のような指定値は使えない。また、分枝切断線が生じるが、そのままの Fair - Luke 法ではこれが 「極が密集したもの」 となるため、その近傍で値が不正確になる※1。さらに、[13]は分数形の項があるため、これも両辺にかけることで解消する必要がある。
 そこで[13]の解をu(z)とすれば、u(z)=sqrt(z)*y(sqrt(z)-1)となるようなy(z)が満たす微分方程式
  • 最適化された第3種Painleve方程式
に対して、拡張された Fair - Luke 法を適用する。これによりy(z)を経由してu(z)が求められる。そのため、指定値は[12]の代わりに
  • 最適化された第3種Painleve方程式の初期値
を採用することとなる。
 以後、漸化式の初期値は、前述の要領で求められるので省略するが、さらに複雑な結果になる。第5種、第6種 Painlevé 方程式も動かない分岐点を持つので、同様のテクニックが必要になる。
 第5種、第6種 Painlevé 方程式、および特定の3階非線形常微分方程式に適用できる[16]相当の微分方程式型はそれぞれ
  • 拡張型Fair-Luke法が適用可能な他の微分方程式
となるが、漸化式や初期値は極めて長大な式となるため、ここに表記できない。別頁 Mathematica関連 にあるコード「Painleve.m」および「Chazy.m」をダウンロードのうえ、ご確認下さい。

【註記】
※1 参考までに、この不都合が生じた状況を再現すると次のようになる。
  • 分枝切断線の近傍の不具合状況

前述のテクニックを用いれば、この不都合を回避できる。
  • 分枝切断線の近傍の不具合が回避された状況
* * *

無限次元三重対角行列の固有値

 元々は、線形代数における三重対角行列法 (または L. Thomas 法) を基礎とする数値計算法である。三重対角行列法とは、三項間の漸化式
  • 三重対角行列法の漸化式
dnXnを求めるために、線形代数方程式
  • 三重対角行列法の線形代数方程式
を解く方法をいう。これはさらに、Gauss の消去法が元になっている。(→ Wikipedia「Tridiagonal matrix algorithm」を参照。)
 この三重対角行列法において行列を無限次元にし、2階線形常微分方程式の固有値を求める場合に適用したものが、標題の数値計算法である。

計算法の説明

 三重対角行列法において、漸化式の項が正負両方向に無限個で、dn=Λ*Xnならば、[1],[2]は次のようになる。
  • 無限次元三重対角行列法(-∞~∞)
この漸化式[3]と線形代数方程式[4]は、2階線形常微分方程式
  • 固有値Λを持つ2階線形常微分方程式
の固有値Λを求める際によく現れる。特に、[5]が楕円体関数の微分方程式のときに効果的で、そのような例を扱った論文が数多くある。当サイトのコードも、その内容を参考にしている※1。
 なお、漸化式が[1]と同様に、
  • 無限次元三重対角行列法(1~∞)の漸化式
のように途中から開始されるときの線形代数方程式は、[2]においてN→∞とした場合を応用した
  • 無限次元三重対角行列法(1~∞)の線形代数方程式
となる。整数次の Mathieu 固有値関数を求める場合などが、この具体例となっている。

【註記】
※1:例えば、
宮崎佳典、浅井信吉、蔡東生、池辺八洲彦
「Mathieu微分方程式の逆固有値問題」
日本応用数理学会論文誌 Vol.8, No.2, (1998), p199-222

P. E. Falloon
「Theory and Computation of Spheroidal Harmonics with General Arguments」
Australia: University of Western Australia, (2001)

P. Abbott
「Tricks of the Trade (Legendre Functions and Spheroial Harmonics)」
The Mathematica Journal Vol.7, Issue1, (1997), p18-22

計算比較(Mathematicaと異なる分枝切断線)

 無限次元三重対角行列による固有値計算法を用いた場合の Mathieu 固有値関数は、Mathematica の実装関数とは異なる分枝切断線となる。(どちらが望ましいということではなく、単に好みの問題である。本来、分枝切断線の形状は任意に決めることができる。)
 因みに Mathematica は、「Newton 法の一般化形式」 を使用しているとのことらしい。

Mathieu 固有値関数Mathieu固有値関数の記号(無限次元三重対角行列による固有値計算法)
  • 無限次元三重対角行列法の場合のグラフ

Mathieu 固有値関数Mathieu固有値関数の記号(Mathematica)
  • Mathematicaの場合のグラフ

Home Menu