倭算数理研究所

科学・数学・学習関連の記事を、「倭マン日記」とは別に書いていくのだ!

spire で多項式 Polynomial を使ってみる (5) : 多項式で表される特殊関数

spire を使ってみるシリーズ(目次)。 今回は多項式で表される特殊関数に関連する spire.math.poly.SpecialPolynomials クラスのメソッドを見ていきます。

SpecialPolynomials クラスでは

  • ルジャンドル多項式legendre()
  • ラゲールの多項式laguerre()
  • チェビシェフの多項式
    • 第1種: chebyshevsFirstKind()
    • 第2種: chebyshevsSecondKind()
  • エルミートの多項式
    • 物理学で使うバージョン: physHermite()
    • 確率論で使うバージョン: probHermite()

の6つが使えます (spire 0.10)。 使い方はどれも同じで、上記のメソッドに Int 値を与えて、長さがその Int 値の Stream を返します。 この Stream の n 番目の要素は各特殊関数のパラメータの値が n となっています。 なぜ Stream を返すのに n 個で必ず切らないといけなくしているのかよく分かりませんが*1、まぁ使い方は別に難しくはないと思います。

各特殊関数の箇所に漸化式を書いているものは spire の実装で使われていることのメモでもあります (spire 0.10)。

以下のサンプルコードでは下記の import 文が書かれているものとします:

import spire.math._
import spire.math.poly.SpecialPolynomials

【参考】

ルジャンドル多項式

ルジャンドル (Legendre) の多項式  { P_n(x) }

  { \displaystyle\begin{align*}
  P_n(x) &= \frac{1}{2^n n!}\frac{d^n}{dx^n}(x^2 - 1)^n
\end{align*}}

で与えられ、以下を満たします:

  { \displaystyle\begin{align*}
  P_0(x) &= 1 \\
  P_1(x) &= x \\
  P_{n+1}(x) &= \frac{1}{n+1}\Big\{(2n+1) xP_n(x) - nP_{n-1}(x)\Big\}
\end{align*}}

spire では SpecialPolynomials.legendre() メソッドで取得できます:

  // Legendre Polynomials
  val ps: Stream[Polynomial[Rational]] = SpecialPolynomials.legendres(5)

  assert( ps(0) == Polynomial("1") )
  assert( ps(1) == Polynomial("x") )
  assert( ps(2) == Polynomial("3/2x^2 - 1/2") )
  assert( ps(3) == Polynomial("5/2x^3 - 3/2x") )
  assert( ps(4) == Polynomial("35/8x^4 - 15/4x^2 + 3/8") )

ラゲールの多項式

ラゲール (Laguerre) の多項式  { L_n(x) }

  { \displaystyle\begin{align*}
  L_n(x)
    &= \sum_{r=0}^n(-1)^r\binom{n}{r}\frac{x^r}{r!} \\[2mm]
    &= \frac{e^x}{n!}\frac{d^n}{dx^n}\left(e^{-x}x^n\right)
    = \frac{1}{n!}\left(\frac{d}{dx} - 1\right)^nx^n
\end{align*}}

で定義され、以下の漸化式を満たします:

  { \displaystyle\begin{align*}
  L_{n+1}(x) = \frac{1}{n+1}\Big\{(2n+1-x)L_n(x) - nL_{n-1}(x)\Big\}
\end{align*}}

spire では SpecialPolynomial.laguerre() メソッドで取得できます:

  // Laguerre Polynomials
  val ls: Stream[Polynomial[Rational]] = SpecialPolynomials.laguerres(5)

  assert( ls(0) == Polynomial("1") )
  assert( ls(1) == Polynomial("-x + 1") )
  assert( ls(2) == Polynomial("1/2x^2 - 2x + 1") )
  assert( ls(3) == Polynomial("-1/6x^3 + 3/2x^2 - 3x + 1") )
  assert( ls(4) == Polynomial("1/24x^4 - 2/3x^3 + 3x^2 - 4x + 1") )

チェビシェフの多項式

第1種チェビシェフの多項式  { T_n(x) }余弦 { n } 倍角の公式が

  { \displaystyle\begin{align*}
  \cos n\theta = T_n(\cos\theta)
\end{align*}}

と表されるように定義されます(余弦 { n } 倍角の公式は『三角関数の n 倍角の公式を導く ~余弦編~』参照)。  { T_n(x) } は以下を満たします:

  { \displaystyle\begin{align*}
  T_0(x) &= 1 \\
  T_1(x) &= x \\
  T_{n+1}(x) &= 2x T_n(x) - T_{n-1}(x)
\end{align*}}

spire では SpecialPolynomial.chebyshevsFirstKind() メソッドで取得できます:

  // Chebyshev Polynomials of the First Kind
  val ts: Stream[Polynomial[Rational]] = SpecialPolynomials.chebyshevsFirstKind(5)

  assert( ts(0) == Polynomial("1") )
  assert( ts(1) == Polynomial("x") )
  assert( ts(2) == Polynomial("2x^2 - 1") )
  assert( ts(3) == Polynomial("4x^3 - 3x") )
  assert( ts(4) == Polynomial("8x^4 - 8x^2 + 1") )

第2種チェビシェフの多項式  { U_n(x) } は正弦の  { n+1 } 倍角の公式が

  { \displaystyle\begin{align*}
  \frac{\sin (n+1)\theta}{\sin\theta} = U_n(\cos \theta)
\end{align*}}

と表されるように定義されます。  { U_n(x) } は以下を満たします:

  { \displaystyle\begin{align*}
  U_0(x) &= 1 \\
  U_1(x) &= 2x \\
  U_{n+1}(x) &= 2x U_n(x) - U_{n-1}(x)
\end{align*}}

これらは  { n = 1 } の場合だけが  { T_n(x) = x } と異なるだけで、特に漸化式は同じです。

spire では SpecialPolynomial.chebyshevsSecondKind() メソッドで取得できます:

  // Chebyshev Polynomials of the Second Kind
  val us: Stream[Polynomial[Rational]] = SpecialPolynomials.chebyshevsSecondKind(5)

  assert( us(0) == Polynomial("1") )
  assert( us(1) == Polynomial("2x") )
  assert( us(2) == Polynomial("4x^2 - 1") )
  assert( us(3) == Polynomial("8x^3 - 4x") )
  assert( us(4) == Polynomial("16x^4 - 12x^2 + 1") )

エルミートの多項式

エルミート (Hermite) の多項式には分野によって2種類の定義があります。 物理学で使われるものは

  { \displaystyle\begin{align*}
  H_n(x) = (-1)^ne^{x^2}\frac{d^n}{dx^n}e^{-x^2} = \left(2x - \frac{d}{dx}\right)^n1
\end{align*}}

で定義されます。 漸化式の形で書けば以下のように書けます:

  { \displaystyle\begin{align*}
  H_0(x) &= 1 \\
  H_1(x) &= 2x \\
  H_{n+1}(x) &= \left(2x - \frac{d}{dx}\right)H_n(x)
\end{align*}}

エルミートの多項式がエルミートの微分方程式の解として求まることは『エルミートの微分方程式とその級数解』で示しました。

spire では SpecialPolynomials.physHermites() メソッドで取得できます:

  // Hermite Polynomials (physics)
  val hs: Stream[Polynomial[Rational]] = SpecialPolynomials.physHermites(5)

  assert( hs(0) == Polynomial("1") )
  assert( hs(1) == Polynomial("2x") )
  assert( hs(2) == Polynomial("4x^2 - 2") )
  assert( hs(3) == Polynomial("8x^3 - 12x") )
  assert( hs(4) == Polynomial("16x^4 - 48x^2 + 12") )

確率論で使われるものは

  { \displaystyle\begin{align*}
  He_n(x) = (-1)^ne^\frac{x^2}{2}\frac{d^n}{dx^n}e^\frac{-x^2}{2} = \left(x - \frac{d}{dx}\right)^n1
\end{align*}}

で定義され、漸化式の形で以下のようにも書けます:

  { \displaystyle\begin{align*}
  He_0(x) &= 1 \\
  He_1(x) &= x \\
  He_{n+1}(x) &= \left(x - \frac{d}{dx}\right)He_n(x)
\end{align*}}

spire では SpecialPolynomials.probHermites() メソッドで取得できます:

  // Hermite Polynomials (prob)
  val hes: Stream[Polynomial[Rational]] = SpecialPolynomials.probHermites(5)

  assert( hes(0) == Polynomial("1") )
  assert( hes(1) == Polynomial("x") )
  assert( hes(2) == Polynomial("x^2 - 1") )
  assert( hes(3) == Polynomial("x^3 - 3x") )
  assert( hes(4) == Polynomial("x^4 - 6x^2 + 3") )

これで多項式関連の API は大体見てきました。 spire.math.poly パッケージには他にもいくつかクラスなどが定義されていますが、直接使う必要はなさそうなもの(根を求めるなど)なので、このシリーズはこの辺で終了。

Scalaスケーラブルプログラミング第3版

Scalaスケーラブルプログラミング第3版

*1:ソースコードを見ると、漸化式で無限 Stream を作った後に単に take() を読んでるだけですが・・・