数値計算でn次多項式の実根を求める方法を書く前に,まずは多項式の値をプログラムで計算するための方法についてまとめます.
n次多項式の値の計算
n次多項式の係数が次数が低い方から$$a_{0}, a_{1}, a_{2}, \ldots , a_{n},$$($$a_{0}$$が定数項)とすると, n次多項式$$f(x)$$は以下のように表せる.
\(f(x)=\sum_{i=0}^{n}a_{i}x^{i}=a_{0}+a_{1}x+a_{2}x^{2}+\ldots+a_{n}x^{n} \ldots (1)\)
これを素直にプログラムに書き下すと以下のようになります.
float a[n+1] = { /* coefficients */ }; float x = /* some value */; float fx = a[0]; for(int i = 1; i <= n; ++i) fx += a[i] * pow(x, i);
このやり方を採用すると,xの冪乗の計算量が無駄になります. (例えば$$x^4$$は$$x^3\cdot x$$として計算できるのに,上のコードだとxの冪乗を個別に計算しています.「ある値の冪乗」に相当する機械語はありませんから,powを呼ぶたびにループが発生していると考えられます.).
ところで式(1)は以下のように変形できます.
\(f(x)=a_{0}+a_{1}x+a_{2}x^{2}+\ldots+a_{n}x^{n}\)
\(f(x)=a_{0}+x(a_{1}+a_{2}x+\ldots+a_{n}x^{n-1})\)
\(f(x)=a_{0}+x(a_{1}+x(a_{2}+\ldots+a_{n}x^{n-2}))\)
\(f(x)=a_{0}+x(a_{1}+x(a_{2}+x(a_{3}\cdots+x(a_{n-1}+x\cdot a_{n})\cdots)))\)
たとえば4次多項式なら
\(f(x)=a_{0}+a_{1}x+a_{2}x^{2}+a_{3}x^{3}+a_{4}x^{4}\)
\(=a_{0}+x(a_{1}+a_{2}x+a_{3}x^{2}+a_{4}x^{3})\)
\(=a_{0}+x(a_{1}+x(a_{2}+a_{3}x+a_{4}x^{2}))\)
\(=a_{0}+x(a_{1}+x(a_{2}+x(a_{3}+a_{4}x)))\)
というような形になります.
このような変形によって計算を楽にする方法をホーナー法というようです. 以下のような級数を考えると簡単にプログラムに書き下すことが出来ます.
\(
f(x)=P_{0} \\
P_{0}=a_{0}+xP_{1} \\
P_{1}=a_{1}+xP_{2} \\
\vdots \\
P_{i}=a_{i}+xP_{i+1} \\
\vdots \\
P_{n}=a_{n} \\
\)
float a[n+1] = { /* coefficients */ }; float x = /* some value */; float fx = a[n]; for(int i = n-1; i >= 0; --i) fx = a[i] + x * fx;
n次多項式の1階導関数の値の計算
(書きかけ)
\(f'(x)=P’_{0}\)
\(P’_{0}=P_{0}+xP’_{1}\)
\(P’_{1}=P_{1}+xP’_{2}\)
\(\vdots\)
\(P’_{i}=P_{i}+xP’_{i+1}\)
\(\vdots\)
\(P’_{n-1}=P_{n-1}\)
float a[n+1] = { /* coefficients */ }; float x = /* some value */; float fx = a[n]; float dfx = 0; for(int i = n-1; i >= 0; --i) { dfx = fx + x * dfx; fx = a[i] + x * fx; }
n次多項式のp階導関数の値の計算
(書きかけ)
\( \\
f^{(p)}(x)=p!\cdot P^{(p)}_{0} \\
P^{(p)}_{0}=P^{(p-1)}_{0}+xP^{(p)}_{1} \\
P^{(p)}_{1}=P^{(p-1)}_{1}+xP^{(p)}_{2} \\
\vdots \\
P^{(p)}_{i}=P^{(p-1)}_{i}+xP^{(p)}_{i+1} \\
\vdots \\
P^{(p)}_{n-p}=P^{(p-1)}_{n-p} \\
\)
おまけ:n次多項式のp階導関数の一般形
\(f^{(p)}(x)=\sum^{n}_{i=p}(n-i)^{\underline{n}}a_{i}x^{i-p}\)
下降階乗
\(x^{\underline{n}}=\prod^{n-1}_{i=0}(x-i)=\frac{x!}{(x-n)!}\)