🌟
直交多項式のゼロ点の数値計算（行列の固有値問題）
行列の固有値問題として求める
ルジャンドル多項式やラゲール多項式などの直交多項式は三項漸化式をを持つ。この性質を利用して、直交多項式のゼロ点を行列の固有値問題に帰着させて求める。
直交多項式
ここで、
である。これを行列の形に表すと、
となる。これを以降は、
と表記する。
となり、ゼロ点
したがって、直交多項式のゼロ点を求める問題は、行列の固有値問題に帰着された。
さらに、対称行列の固有値問題まで帰着させることもできる。
両辺に、適当な定数
となる。対称行列になるには、
を満たせば良い。すなわち、
よって、
となる。したがって，
の固有値が
Rust による実装
Rust で実装する。固有値問題を解く部分は、線形代数ライブラリnalgebraを使用する。
（2022/10/12 追記）
nalgebra ライブラリの SymmetricEigen構造体によって固有値と固有ベクトルを求めると，自由度 N が大きいときに誤差が大きくなるようだ．
nalgebra はnalgebra-lapackクレートで，LAPACK を用いた実装も提供している．今回はこちらを使ったほうが良さそうである．
まずは，任意の直交多項式のゼロ点を計算する trait Zeros を実装する．
use nalgebra::DMatrix;
use nalgebra_lapack::SymmetricEigen;
trait Zeros {
// 直交多項式の自由度 N を求める
fn degree(&self) -> usize;
// 三項漸化式の係数 a_n, b_n, c_n を求める
fn coefficients(&self) -> (Vec<f64>, Vec<f64>, Vec<f64>);
// 対称三重対角行列 S_N を求める
// 返り値の型は，nalgebra::DMatrix<f64>
fn symmetric_jacobi(&self) -> DMatrix<f64> {
let n = self.degree();
let (a, b, c) = self.coefficients();
let mut symmetric_coeff = Vec::with_capacity(n - 1);
for i in 0..n - 1 {
symmetric_coeff.push(libm::sqrt(c[i + 1] * a[i]));
}
DMatrix::<f64>::from_fn(n, n, |r, c| {
if r == c {
b[r]
// ex. (2,1), (3,2), etc..
} else if r == c + 1 {
symmetric_coeff[c]
// ex. (1,2), (2,3), etc..
} else if r == c - 1 {
symmetric_coeff[r]
} else {
0.
}
})
}
// ゼロ点を求める
// 返り値はVec<f64>
fn zeros(&self) -> Vec<f64> {
let sym_j = self.symmetric_jacobi();
let eigen = SymmetricEigen::new(sym_j);
eigen.eigenvalues.data.into()
}
}
例: ラゲール多項式
次に，具体的な直交多項式に対するゼロ点を計算する構造体を実装する．
今回は例として，（一般化された）ラゲール多項式のゼロ点を計算する．
(一般化された)ラゲール多項式
今回は，Laguerre 構造体を定義し，この Laguerre 構造体に対して上記で実装した Zeros trait を実装する．
// Laguerre 構造体
// n は自由度
// alpha は一般化されたラゲール多項式の係数 alpha > -1
pub struct Laguerre {
n: usize,
alpha: f64,
}
impl Laguerre {
// Laguerre 構造体を生成する関連関数
pub fn new(n: usize, alpha: f64) -> Self {
if alpha <= -1. {
panic!("can't define Laguerre polynomials for alpha <= -1")
}
Self { n, alpha }
}
}
// Laguerre 構造体に対するZeros traitの実装
impl Zeros for Laguerre {
// 自由度 N を取得する
fn degree(&self) -> usize {
self.n
}
fn coefficients(&self) -> (Vec<f64>, Vec<f64>, Vec<f64>) {
let mut a = Vec::with_capacity(self.n);
let mut b = Vec::with_capacity(self.n);
let mut c = Vec::with_capacity(self.n);
for i in 0..self.n {
// a_i = -(i + 1)
a.push(-(i as f64 + 1.));
// b_i = 2i + alpha + 1
b.push(2. * i as f64 + self.alpha + 1.);
// c_i = -(i + alpha)
c.push(-(i as f64 + self.alpha));
}
(a, b, c)
}
}
これで，ラゲール多項式のゼロ点を求めることができる．
#[cfg(test)]
mod tests {
use super::{Laguerre, Zeros};
#[test]
fn test_laguerre() {
let l = Laguerre::new(2, 0.);
let zeros = l.zeros();
let zeros_rounding: Vec<f64> = zeros
.iter()
.map(|z| (z * 10000.).floor() / 10000.)
.collect();
let want = vec![3.4142, 0.5857];
println!("{:?}", zeros);
assert_eq!(zeros_rounding, want);
}
running 1 test
[3.414213562373095, 0.5857864376269049]
test special::tests::test_laguerre ... ok
test result: ok. 1 passed; 0 failed; 0 ignored; 0 measured; 9 filtered out; finished in 0.00s
参考文献
- A. Gil, J. Segura, and, N.M. Temme, Numerical Methods for Special Functions. Society for Industrial and Applied Mathematics, 2007.
Discussion