ルジャンドル・ガウス求積法で積分を解くプログラムを書いています。n 次求積のアルゴリズムでは、ある時点で、n 次ルジャンドル多項式 Pn(x) の根を見つけ、それらを配列 Absc (「abscissa」) に割り当てる必要があります。Pn は区間 [-1,1] に n 個の独立した実根を持つ n 次多項式です。ライブラリからルートをインポートするだけでなく、ルートを計算できるようにしたいと考えています。PCoeff と呼ぶ多項式係数を与える配列を作成できます。私が試したルーツを見つけるために
Absc = numpy.roots(PCoeff)
これは約 n = 40 までは機能しますが、それを超えると失敗し始め、本来あるべきではない複雑な根を与えます。また、次を使用して多項式を定義しようとしました
P = numpy.poly1d(PCoeff)
Absc = P.r
しかし、これはおそらく同じnumpyルート検索アルゴリズムを使用しているため、同じ問題を引き起こします。
有望と思われる別の方法は、scipy.optimize.fsolve(Pn, x0) を使用することでした。ここで、x0 は根での私の推測の n 要素配列です。これに関する問題は、私の x0 の選択によっては、このメソッドが 1 つの特定のルートを他のルートの代わりに複数回与える可能性があることです。[-1,1] 上の等距離点として x0 を入力しようとしました
x0 = numpy.zeros(n)
step = 2./(n-1)
for i in xrange(n):
x0[i] = -1. + i*step
しかし、n = 5 になると、fsolve はいくつかの根を繰り返し与え、他の根を無視します。また、numpy.roots の結果を x0 として使用してみました。ただし、np.roots が複雑な値を与える問題領域では、これらは fsolve でエラーを引き起こします。
TypeError: array cannot be safely cast to required type
動作する可能性のある scipy.optimize.roots() ルーチンがあることをオンラインで見ましたが、私のコンピューターの scipy ライブラリにはありません。このコンピューターにダウンロードする権限がないため、更新が面倒です。
高精度のために 64 次で求積を実行できるようにしたいのですが、このルート検索は失敗を引き起こします。何か案は？