orthopoly (Maxima 5.47post Manual)

71 orthopoly

Introduction to orthogonal polynomials
Functions and Variables for orthogonal polynomials

71.1 Introduction to orthogonal polynomials

orthopolyは、 Chebyshev, Laguerre, Hermite, Jacobi, Legendre, 超球(Gegenbauer) 多項式を含むいくつかの種類の直交多項式のシンボリックな評価と数値評価のためのパッケージです。さらに、orthopolyには球Bessel, 球Hankel, 球調和関数のサポートが含まれます。

ほとんどの部分に関して、 orthopolyは Abramowitz and StegunのHandbook of Mathematical Functions, Chapter 22 (10th printing, December 1972)の慣例に従います; 加えて、 Gradshteyn and RyzhikのTable of Integrals, Series, and Products (1980 corrected and enlarged edition)と Eugen MerzbacherのQuantum Mechanics (2nd edition, 1970)を使います。

University of Nebraska at Kearney (UNK)のBarton Willisが orthopolyパッケージとドキュメンテーションを書きました。パッケージはGNU General Public License (GPL)の下で公開されています。

Orthogonal polynomials · Share packages · Package orthopoly ·

Getting Started with orthopoly
Limitations
Floating point Evaluation
Graphics and orthopoly
Miscellaneous Functions
Algorithms

71.1.1 Getting Started with orthopoly

load ("orthopoly")はorthopolyパッケージをロードします。

3次のLegendre多項式を見つけるには、

(%i1) legendre_p (3, x);
                      3             2
             5 (1 - x)    15 (1 - x)
(%o1)      - ---------- + ----------- - 6 (1 - x) + 1
                 2             2

これをxの冪の和として表すには、 ratsimpかratを結果に適用してください。

(%i2) [ratsimp (%), rat (%)];
                        3           3
                     5 x  - 3 x  5 x  - 3 x
(%o2)/R/            [----------, ----------]
                         2           2

あるいはまた、 legendre_pの第二引数 (「主」変数)を正準有理形(CRE)にしてください。

(%i1) legendre_p (3, rat (x));
                              3
                           5 x  - 3 x
(%o1)/R/                   ----------
                               2

浮動小数点評価に関して、 orthopolyはランニング誤差解析を使って誤差の上限を推定します。例えば、

(%i1) jacobi_p (150, 2, 3, 0.2);
(%o1) interval(- 0.062017037936715, 1.533267919277521E-11)

区間(interval)は形式form interval (c, r)を取ります。ここで、cは中央値、rは区間の半径です。 Maximaは区間上の算術をサポートしていないので、グラフィックスなどいくつかの状況では誤差を抑制し、区間の中央値だけ出力したいでしょう。これをするには、オプション変数orthopoly_returns_intervalsを falseに設定してください。

(%i1) orthopoly_returns_intervals : false;
(%o1)                         false
(%i2) jacobi_p (150, 2, 3, 0.2);
(%o2)                  - 0.062017037936715

更に知るにはセクションsee Floating point Evaluationを参照してください。

orthopolyのほとんどの関数はgradefプロパティを持ちます; 例えば、

(%i1) diff (hermite (n, x), x);
(%o1)                     2 n H     (x)
                               n - 1
(%i2) diff (gen_laguerre (n, a, x), x);
              (a)               (a)
           n L   (x) - (n + a) L     (x) unit_step(n)
              n                 n - 1
(%o2)      ------------------------------------------
                               x

二番目の例の単位階段関数は、 nが0で評価することによって、そうでなければ生じる誤差を抑制します。

(%i3) ev (%, n = 0);
(%o3)                           0

gradefプロパティは「主」変数にのみ適用されます; 他の引数に関する導関数は、普通、エラーメッセージに帰着します; 例えば、

(%i1) diff (hermite (n, x), x);
(%o1)                     2 n H     (x)
                               n - 1
(%i2) diff (hermite (n, x), n);

Maxima doesn't know the derivative of hermite with respect the first
argument
 -- an error.  Quitting.  To debug this try debugmode(true);

一般に、orthopolyの関数はリストや行列上に写像します。写像を完全に評価するには、オプション変数doallmxopsとlistarithはともに true(デフォルト値)でなければいけません。行列上への写像を見るには、以下を考えてください。

(%i1) hermite (2, x);
                                     2
(%o1)                    - 2 (1 - 2 x )
(%i2) m : matrix ([0, x], [y, 0]);
                            [ 0  x ]
(%o2)                       [      ]
                            [ y  0 ]
(%i3) hermite (2, m);
               [                             2  ]
               [      - 2        - 2 (1 - 2 x ) ]
(%o3)          [                                ]
               [             2                  ]
               [ - 2 (1 - 2 y )       - 2       ]

二番目の例では、値のi, j要素はhermite (2, m[i,j])です; 次の例で見るように、これは、計算-2 + 4 m . mと同じではありません。

(%i4) -2 * matrix ([1, 0], [0, 1]) + 4 * m . m;
                    [ 4 x y - 2      0     ]
(%o4)               [                      ]
                    [     0      4 x y - 2 ]

定義域外の点で関数を評価すると、一般に、orthopolyは未評価関数を返します。例えば、

(%i1) legendre_p (2/3, x);
(%o1)                        P   (x)
                              2/3

orthopolyはTexへの翻訳をサポートします; 端末上での2次元出力も行います。

(%i1) spherical_harmonic (l, m, theta, phi);
                          m
(%o1)                    Y (theta, phi)
                          l
(%i2) tex (%);
$$Y_{l}^{m}\left(\vartheta,\varphi\right)$$
(%o2)                         false
(%i3) jacobi_p (n, a, a - b, x/2);
                          (a, a - b) x
(%o3)                    P          (-)
                          n          2
(%i4) tex (%);
$$P_{n}^{\left(a,a-b\right)}\left({{x}\over{2}}\right)$$
(%o4)                         false

71.1.2 Limitations

式がいくつかの直交多項式を記号順で含む時、式が実際に0になる可能性がありますが、まだMaximaはそれを零に整理することができません。もしそんな量で割るなら、困ったことになるでしょう。例えば、以下の式は 1より大きな整数 nで零になりますが、まだMaximaはそれを零に整理することができません。

(%i1) (2*n - 1) * legendre_p (n - 1, x) * x - n * legendre_p (n, x)
      + (1 - n) * legendre_p (n - 2, x);
(%o1)  (2 n - 1) P     (x) x - n P (x) + (1 - n) P     (x)
                  n - 1           n               n - 2

特定の nでは式を零に換算できます。

(%i2) ev (% ,n = 10, ratsimp);
(%o2)                           0

一般に、直交多項式の多項式形は浮動小数点評価に関しては不適当です。以下は例です。

(%i1) p : jacobi_p (100, 2, 3, x)$

(%i2) subst (0.2, x, p);
(%o2)                3.4442767023833592E+35
(%i3) jacobi_p (100, 2, 3, 0.2);
(%o3)  interval(0.18413609135169, 6.8990300925815987E-12)
(%i4) float(jacobi_p (100, 2, 3, 2/10));
(%o4)                   0.18413609135169

真値は約0.184です; この計算は極端な減算消去誤差に苦しみます。多項式を展開し評価すると、よりよい結果を与えます。

(%i5) p : expand(p)$
(%i6) subst (0.2, x, p);
(%o6) 0.18413609766122982

これは一般的な規則ではありません; 多項式を展開することはいつも、数値評価により適した式を生じるわけではありません。数値評価する最もよい方法は断然、1つ以上の関数引数を浮動小数点数にすることです。それをすることで、特別な浮動小数点アルゴリズムが評価に使われます。

Maximaの float関数は幾分でたらめです; もし floatを記号次数や順序パラメータを持つ直交多項式を含む式に適用するなら、これらのパラメータは浮動小数点に変換されるかもしれません; その後、式は完全には評価されません。以下を考えてください。

(%i1) assoc_legendre_p (n, 1, x);
                               1
(%o1)                         P (x)
                               n
(%i2) float (%);
                              1.0
(%o2)                        P   (x)
                              n
(%i3) ev (%, n=2, x=0.9);
                             1.0
(%o3)                       P   (0.9)
                             2

(%o3)の式は浮動小数点に評価されません; orthopolyは整数を要求するところで浮動小数点値を認識しません。同様に、数値評価 recognize floating point values where it requires an integer. Similarly, numerical evaluation of the pochhammer_max_indexを越える位数の pochhammer関数は迷惑かもしれません; 以下を考えてください。

(%i1) x :  pochhammer (1, 10), pochhammer_max_index : 5;
(%o1)                         (1)
                                 10

floatを適用することは xを浮動小数点に評価しません。

(%i2) float (x);
(%o2)                       (1.0)
                                 10.0

xを浮動小数点に評価するには、 pochhammer_max_indexを11以上にバインドして、 floatを xに適用する必要があります。

(%i3) float (x), pochhammer_max_index : 11;
(%o3)                       3628800.0

pochhammer_max_indexのデフォルト値は100です; orthopolyをロードした後、値を変えてください。

最後に、参考書は直交多項式の定義を変えることを承知してください; 一般的にAbramowitz and Stegunの慣例を使っています。

orhtopolyのバグを疑う前に、いくつかの特殊なケースをチェックして、あなたの定義がorthopolyが使っているものと一致しているかを明らかにしてください。

定義はしばしば規格化について異なります; 時々、著者は (-1, 1)以外の区間上で直交な族を作る関数の「シフト」版を使います。例えば、 (0, 1)上で直交するLegendre多項式を定義するのに、以下を定義します。

(%i1) shifted_legendre_p (n, x) := legendre_p (n, 2*x - 1)$

(%i2) shifted_legendre_p (2, rat (x));
                            2
(%o2)/R/                 6 x  - 6 x + 1
(%i3) legendre_p (2, rat (x));
                               2
                            3 x  - 1
(%o3)/R/                    --------
                               2

71.1.3 Floating point Evaluation

orthopolyの関数のほとんどは浮動小数点評価中の誤差を見積もるのに、ランニング誤差解析を使います; 例外は球Bessel関数と第二種Legendreの陪多項式です。数値評価のため、球Bessel関数はSLATEC関数をコールします。第二種Legendreの陪多項式の数値評価のために特別な方法は使われません。

ランニング誤差解析は (丸め単位としても知られている)計算機イプシロンの二次か高次の誤差を無視します。 2,3の他の誤差も無視します。 (ありそうにありませんが、)実際の誤差は推定を越える可能性があります。

区間は形式 interval (c, r)を持ちます。ここで、 cは区間の中心で、 rは半径です。区間の中心は複素数であり得ますし、半径はいつも正の実数です。

以下は例です。

(%i1) fpprec : 50$

(%i2) y0 : jacobi_p (100, 2, 3, 0.2);
(%o2) interval(0.1841360913516871, 6.8990300925815987E-12)
(%i3) y1 : bfloat (jacobi_p (100, 2, 3, 1/5));
(%o3) 1.8413609135168563091370224958913493690868904463668b-1

実際の誤差が誤差推定よりも小さいことをテストしましょう。

(%i4) is (abs (part (y0, 1) - y1) < part (y0, 2));
(%o4)                         true

なるほど、この例では、誤差推定は真の誤差の上限です。

Maximaは区間の算術をサポートしていません。

(%i1) legendre_p (7, 0.1) + legendre_p (8, 0.1);
(%o1) interval(0.18032072148437508, 3.1477135311021797E-15)
        + interval(- 0.19949294375000004, 3.3769353084291579E-15)

ユーザーは区間算数を行う計算をする演算子を定義できます。区間の足し算を定義するには、以下を定義できます。

(%i1) infix ("@+")$

(%i2) "@+"(x,y) := interval (part (x, 1) + part (y, 1), part (x, 2)
      + part (y, 2))$

(%i3) legendre_p (7, 0.1) @+ legendre_p (8, 0.1);
(%o3) interval(- 0.019172222265624955, 6.5246488395313372E-15)

引数が複素数の時、特殊な浮動小数点ルーチンがコールされます。例えば、

(%i1) legendre_p (10, 2 + 3.0*%i);
(%o1) interval(- 3.876378825E+7 %i - 6.0787748E+7, 
                                           1.2089173052721777E-6)

これを真値と比較しましょう。

(%i1) float (expand (legendre_p (10, 2 + 3*%i)));
(%o1)          - 3.876378825E+7 %i - 6.0787748E+7

更に、引数が多倍長浮動小数点の時、特殊な浮動小数点ルーチンがコールされます; しかしながら、多倍長浮動小数点は倍精度浮動小数点に変換され、最終結果は倍精度です。

(%i1) ultraspherical (150, 0.5b0, 0.9b0);
(%o1) interval(- 0.043009481257265, 3.3750051301228864E-14)

71.1.4 Graphics and `orthopoly`

直交多項式を含む式をプロットするには、 2つのことをしなければいけません:

オプション変数 orthopoly_returns_intervalsを falseに設定する。
orthopoly関数のすべてのコールをクォートする。

もし関数コールがクォートされていないなら、 Maximaはプロットする前にそれらを多項式に評価します; 結果として、特殊な浮動小数点コードはコールされません。以下は、Legendre多項式を含む式をどうやってプロットするかの例です。

(%i1) plot2d ('(legendre_p (5, x)), [x, 0, 1]),
                        orthopoly_returns_intervals : false;
(%o1)

式 legendre_p (5, x)全体をクォートします; これは 'legendre_p (5, x)を使って関数名をクォートするだけとは違います。

Plotting ·

71.1.5 Miscellaneous Functions

orthopolyパッケージは Pochhammerシンボルと単位階段函数を定義します。 orthopolyはgradef文の中で Kroneckerのデルタ函数と単位階段函数を使います。

Pochhammerシンボルをガンマ函数の商に変換するには、 makegammaを使ってください。

(%i1) makegamma (pochhammer (x, n));
                          gamma(x + n)
(%o1)                     ------------
                            gamma(x)
(%i2) makegamma (pochhammer (1/2, 1/2));
                                1
(%o2)                       ---------
                            sqrt(%pi)

Pochhammerシンボルの導函数は psi函数を使って与えられます。

(%i1) diff (pochhammer (x, n), x);
(%o1)             (x)  (psi (x + n) - psi (x))
                     n     0             0
(%i2) diff (pochhammer (x, n), n);
(%o2)                   (x)  psi (x + n)
                           n    0

(%o1)の式に注意する必要があります; psi函数の差分は x = -1, -2, .., -nの時多項式です。これらの多項式は nが正の整数の時、導函数をn - 1次多項式にするように、 pochhammer (x, n)の因子を相殺します。

Pochhammerシンボルはガンマ函数の商としての表現を通して負の位数で定義されます。以下を考えてください。

(%i1) q : makegamma (pochhammer (x, n));
                          gamma(x + n)
(%o1)                     ------------
                            gamma(x)
(%i2) sublis ([x=11/3, n= -6], q);
                               729
(%o2)                        - ----
                               2240

代わりに、この結果を直接得ることができます。

(%i1) pochhammer (11/3, -6);
                               729
(%o1)                        - ----
                               2240

単位階段函数は左連続です; 従って、

(%i1) [unit_step (-1/10), unit_step (0), unit_step (1/10)];
(%o1)                       [0, 0, 1]

もし零で左連続でも右連続でもない単位階段函数が必要なら、 signumを使って自分のものを定義してください; 例えば、

(%i1) xunit_step (x) := (1 + signum (x))/2$

(%i2) [xunit_step (-1/10), xunit_step (0), xunit_step (1/10)];
                                1
(%o2)                       [0, -, 1]
                                2

unit_step自身を再定義しないでください; orthopolyの中のあるコードは単位階段函数が左連続であることを要求します。

71.1.6 Algorithms

一般的に、 orthopolyは直交多項式の超幾何表現を使うことで記号評価をします超幾何函数は (ドキュメント化されていない)関数 hypergeo11と hypergeo21を使って評価されます。例外は半整数Bessel函数と第二種Legendreの陪函数です。半整数Bessel函数は明示的な表現を使って評価されます。第二種Legendreの陪函数は再帰を使って評価されます。

浮動小数点評価のために、函数のほとんどを超幾何形式に再び変換します; 順方向再帰を使って超幾何函数を評価します。ここでも、例外は半整数Bessel函数と第二種Legendreの陪函数です。数値的に、半整数Bessel函数はSLATECコードを使って評価されます。

71.2 Functions and Variables for orthogonal polynomials

関数: assoc_legendre_p (n, m, x) ¶

次数 nと位数 mの第一種Legendre陪函数。

参考文献: Abramowitz and Stegun, equations 22.5.37, page 779, 8.6.6 (second equation), page 334, and 8.2.5, page 333.