Mathematica関連 Menu

Mathematicaの落とし穴

プログラミング中、実際に遭遇した落とし穴の事例。
(注意:Ver.9以降では動作が異なる可能性もあります。)
Mathematicaの落とし穴の絵

CoefficientListによるリストの長さの判定

 任意の関数を純関数で指定する引数を持ち、その関数を Maclaurin 展開したときの係数と、その係数リストの長さを内部計算で必要とするようなコードを書いていたとする。ここでは簡単のため、0次~10次の係数 (つまり、11個の要素からなる長さ11のリスト) を仮定する。
 このとき、引数関数が余弦関数ならば、
  • CoefficientListを余弦関数に使用した場合

となり、意図した結果となる。しかし正弦関数ならば、
  • CoefficientListを正弦関数に使用した場合

となってしまう。もし、内部コードでリストの長さの判定条件にこれを使用していたら、エラーが発生する。よって Join[list, {0}] 等を用いて、リストの末尾に要素0を補う必要がある。
 CoefficientList は、元々多項式の係数を求める組み込み関数である。そこで、次のように SeriesCoefficient を使用するのが本来の方法であり、かつオールマイティーな方法であると思われる。
  • SeriesCoefficientを正弦関数に使用した場合

 なお、他の理由、例えば計算速度に関して CoefficientList と SeriesCoefficient のどちらを使うのが望ましいかは、個々の問題や条件によって異なってくると思われるので、実際には動作テストで確認した方がいい。

行列を指定するパターンマッチングでの混乱

 変数がリスト (ベクトル) のときは、パターンマッチングで次のように指定する。
  • リスト(ベクトル)のパターンマッチング

 一方、変数が行列のとき、パターンマッチングで次のように指定すると、うまくいかない。
 (パターンにマッチしないと判断されたため、入力した命令文がそのまま返されるだけとなる。)
  • (誤った)行列のパターンマッチング

 行列のときは、次のように指定すれば意図した動作をする。
 ( Mathematica の仕様によるものであり、全く理由はない。)
  • (正しい)行列のパターンマッチング

 これは特に、行列変数とベクトル変数が混在するようなコード内で混同しやすい

ループ型? or 関数型?

 冪級数の形の級数を While ループで記述しようとして、次の2通りのコードを書いた。
 (級数展開式の初項が1となるように定数倍した、重みk正規化 Eisenstein 級数E k (τ)の例。)
  • 2通りのループ型プログラミング例

 後者は意図した結果になる。しかし、前者は級数冪級数になってしまう。

 前者のコードを用いて Klein の楕円モジュラー関数J(τ)を求めようとしたが、誤った結果になっている。
  • 楕円モジュラー関数の(誤った)プログラミング結果

 後者のコードを用いれば、正しいJ(τ)のグラフになる。
  • 楕円モジュラー関数の(正しい)プログラミング結果

 勿論、局所変数q1は必須ではなく、もっと簡単な次のコードでも正しい結果になる。
  • 簡略化したループ型プログラミング例

 しかしこのコードは、前の例よりも計算が遅いのである!。
  • 楕円モジュラー関数の(正しいが遅い)プログラミング結果

 高速化の定石は、ループ型ではなく、関数型プログラミングを用いることである。しかし、後者にも思わぬ落とし穴がある。以下、その一例を示す。
 前にループ型でプログラミングしたコードes2が、虚部の大きさに応じて反復回数をどのように決めているかを確認する。
  • ループの反復回数を分析するコード

 実部が異なるいくつかの虚軸方向で描画した、反復回数のグラフ。概ね一次分数関数に近いように見える。
  • ループの反復回数のグラフ

 そこで、この曲線に概ね適合する具体的な一次分数関数を求める。
  • 反復回数曲線の適合関数を求める

 先の反復回数のグラフと重ねて描画してみると、よくフィットしていることが分かる。
  • 反復回数曲線の適合関数のグラフ

 この結果を用いて、関数型プログラミングを作成する。
  • 関数型プログラミング例1

 これまでで最も高速である。
  • 関数型プログラミング例1の結果

 もし、上記のような加算回数の見積を行わず、単に固定回数で関数型プログラミングを作成すると、どのような結果になるか確認する。
  • 関数型プログラミング例2

 実軸の近傍では加算回数不足のため不正確になっている。逆に虚部が大きい領域では回数が冗長なため余分に時間がかかり、全体として遅くなっている。関数型プログラミングは記述しやすく高速だが、被プログラミング関数の特性をよく踏まえないと、この事例のように本来の利点が発揮できない場合もある (つまり万能ではない) 。
  • 関数型プログラミング例2の結果

等高線は悩ましい

【その1:Throw::nocatch エラー】
 MeshFunctions は、等高線を引く位置を決める関数を指定するオプションである。その位置関数は 「純関数」 で表現することになっており、純関数は Function または Slot (#&) の2通りの方法で記述できる。
 次のコードは、Function で記述した場合であるが、PlotRange の上限値が0に近いグラフでは、絶対値の等高線が乱れる。(以降、絶対値の等高線は黒、偏角の等高線は白で描画する。)
  • 等高線をFunctionで記述した場合

 純関数を Slot で記述した場合は、エラー Throw::nocatch が発生して計算が強制終了する。
  • 等高線をSlotで記述した場合

 ParabolicCylinderDHermiteH で記述できる。そうすれば Throw::nocatch は発生しないので、このエラーの原因は ParabolicCylinderD にあることが分かる。(しかし、その理由は不明。)
 因みに、純関数を Slot で記述した場合、絶対値の等高線は乱れない。
  • ParabolicCylinderDをHermiteHで記述した場合

 後述するとおり、ParabolicCylinderD と同様のエラーが発生する Mathematica 組込の複素関数は、いくつか存在する。
 どの組込関数であってもエラーが発生せず、しかも常に絶対値の等高線が乱れない方法として、次の (やや奇妙な) 折衷策が有効となる。この方法は、MathematicalFunctionPlot.m にも採用している。
  • 等高線をSlotとFunctionで記述した場合

 偏角の等高線にはなっていないが、純関数 Slot で#1および#2を用いること自体に問題はないことが、次のコードから分かる。
  • #1,#2自体に問題はない事の確認描画

 恐らく、#1および#2を ParabolicCylinderD 内で使用することに問題がある。(つまり、局所記号fが原因ではない。)
  • ParabolicCylinderDに問題がある事の確認描画

 因みに ColorFunction も純関数で指定するオプションである。これも同様に、純関数を Function で記述した場合はエラーにならないが、Slot の場合はエラーが発生する。(MathematicalFunctionPlot.m は Function を使用している。)
  • ColorFunctionをFunctionで記述した場合
  • ColorFunctionをSlotで記述した場合

 そこで、Mathematica 組込の複素関数のほとんど全て (ただし、Ver.8 時点) について、ParabolicCylinderD と同様の問題が発生するか調査した。
 次は、属性が 「NumericFunction」 である組込関数名をシステムから抽出する方法である。
  • システムから組込関数名を抽出する

 上記のうち、複素変数グラフを描画する可能性のある組込関数を手動で選び★、その各々について次の2種類の描画テストを実行し、「Test-1」 は問題なく描画できる事を、「Test-2」 は問題が発生するかどうかを確認する。
 (f[x_,y_]:=の右辺、およびprMaxは、その都度 書き換えて実行する。★で取り除いた関数の例として、変数が自然数のみで定義される AiryAiZero、複素関数を考える事ができても Mathematica では引数が実数に限る InverseErf, BSplineBasis 等がある。)
  • 描画Test-1
  • 描画Test-2

(調査結果)
 2種類の描画テストのうち、Test-1 は以外の全ての組込関数で描画可能、Test-2 は以外の組込関数ならば描画可能であることが判明した。
次の2個の組込関数は、Throw::nocatch が発生する。
 ParabolicCylinderD, SiegelTheta

次の9個の組込関数は、Throw::nocatch は発生しないものの計算が終わらず、カーネルを強制終了するしかなくなる (「評価を放棄」 が効かなくなる) 問題が発生する。その発生タイミングから、恐らくと同様に MeshFunctions における偏角等高線の純関数に原因があると思われる。
 StruveH, StruveL, AngerJ, WeberE, MarcumQ, OwenT, HypergeometricPFQ, HypergeometricPFQRegularized, MeijerG

次の3個の組込関数は、描画テストで使用したプロット方法では、複素関数自体がうまく描画できないか、制御が難しい組込関数である。この問題は、当頁で関心のあるとは全く別物である。
 EllipticExp, EllipticExpPrime, EllipticLog
 (については、描画範囲, 描画精度, 関数や引数の定数倍, 等高線の本数または間隔等を、様々 (整数の他にも、無理数を使用する等) に取り替えてテストしましたが、問題の発生状況に違いは無いようです。このことから、ある特定の値が Throw::nocatch を発生させている可能性も低いです。なお、Ver.9以降での動作、他のPC環境での動作がどうなるか大変気になりますが、これについては確認していません。)


【その2:偏角等高線の彩色】
 偏角の等高線を色相環に基づいて彩色する。しかし、明らかに ColorFunction の色相と一致していない。
  • 色相が一致していない描画例

 等高線の本数は、この問題と関係ない。
  • 色相が一致していない描画例(等高線本数の増)

 MeshFunctions の純関数を Slot (#&) ではなく Function で記述しても、色相の不一致が起きる。
  • 色相が一致していない描画例(純関数の変更)

 Mesh の指定から、j=0となるインデックスのみを取り除くと、なぜか色相が一致する。(理由は不明。)
  • インデックスからj=0のみを除去した描画例

 よって、次のように記述すれば良い。(ただし、偏角が±πの等高線は太くなる事が、別の問題として残る。)
  • 色相が一致する描画例

 正弦関数の場合。
  • 色相が一致する描画例(正弦関数の場合)

 色相の不一致は、Plot3D でも起きる。
  • 色相が一致していない描画例(Plot3Dの場合)

 DensityPlot での方法を応用すれば、同様に色相は一致する。(しかし、依然として理由は分からない…。一方、偏角が±πの等高線のみ太くなる理由ならばこの図から分かる。すなわち、複数の等高線が重なっていることによる。)
  • 色相が一致する描画例(Plot3Dの場合)


*******


  • Mathematicaでの誤った入力の例
何となく悪い予感が…

  • Mathematicaでの誤った入力の例2
なぜ意図した結果にならないの? orz

Mathematica関連 Menu