在python中查找勒讓德多項式的根

我正在編寫一個程序，它解決了勒讓德 - 高斯求積分的積分問題。 n階正交算法需要在一個點上找到n階Legendre多項式Pn（x）的根，並將它們分配給數組Absc（'橫座標'）。 Pn是在區間[-1,1]上具有n個獨立實根的n階多項式。我希望能夠計算根，而不是從某個庫中導入它們。我可以創建一個給出多項式係數的數組，我稱之爲PCoeff。爲了找到根我試圖在python中查找勒讓德多項式的根

Absc = numpy.roots(PCoeff)

這對於高至約N = 40，但除此之外，它開始出現故障，給複雜的根源時，它確實不應該。我也試着定義使用

P = numpy.poly1d(PCoeff) 
Absc = P.r

多項式，但是這給了同樣的問題，大概是因爲它使用相同的numpy的求根算法。

另一種看起來很有前途的方法是使用scipy.optimize.fsolve（Pn，x0），其中x0是我在根上猜測的一個n元素數組。問題在於，根據我的x0選擇，這種方法可能會給一個特定的根不止一次替代其他根。我試着在[-1,1]

x0 = numpy.zeros(n) 
step = 2./(n-1) 
for i in xrange(n): 
    x0[i] = -1. + i*step

填充X0爲等距點，但一旦我到n = 5，fsolve給出了一些根重複，而忽略其他。我也嘗試使用numpy.roots的結果作爲x0。然而，在問題區域，其中np.roots得到複雜的價值觀，這些引起fsolve錯誤

TypeError: array cannot be safely cast to required type

我在網上看到有一個scipy.optimize.roots（）函數可以工作，但它不是在我電腦上的scipy庫中。更新很麻煩，因爲我沒有權限在這臺電腦上下載東西。

我希望能夠以64的順序運行正交以獲得高精度，但是這個根找不到。有任何想法嗎？

來源

2012-08-03 johndmalcolm

您是否知道numpy.polynomial.legendre模塊？它包含找到根源的根源方法。你可以在這裏檢查： http://docs.scipy.org/doc/numpy/reference/routines.polynomials.legendre.html – Teudimundo 2012-08-03 12:16:48

這是一個自豪的能夠自己計算它，但這個模塊工作很好地爲我所需要的。比迄今爲止我所能做到的要好得多。謝謝一堆。 – johndmalcolm 2012-08-03 14:30:11

由於np.roots依賴於「找到伴隨矩陣的特徵值」，如文檔所述，您可能會遇到錯誤傳播問題，從而導致根上的非零虛部。也許你可以放棄使用np.real函數的虛部。

你可以嘗試用不同的方式來計算使用根泰勒近似根：

https://math.stackexchange.com/questions/12160/roots-of-legendre-polynomial

來源

2012-08-03 12:10:26

我考慮了你所做的第一點，但不幸的是，計算出的根的實際部分與實際值不匹配。我想知道模數是否...我應該檢查一下。無論如何，在上面的評論中指出的legendre模塊有一個計算我需要的根的例程。雖然使用它有點像作弊。通過您提供的鏈接，我找到了一篇論文，描述了計算根的方法，我計劃調查以確定我是否可以自己真正計算這些東西。謝謝！ – johndmalcolm 2012-08-03 14:34:09

我並不知道legendre模塊。實際上，legroots的代碼非常短：它執行3個步驟：查找伴隨矩陣，計算特徵值並對它們進行排序。不知道爲什麼np.roots無法執行。也許問題在於配對矩陣的計算。 – 2012-08-03 15:42:09

你可以找到its documentation一個簡單的例子，實現您的問題SymPy：

>>> for n in range(5): 
...  nprint(polyroots(taylor(lambda x: legendre(n, x), 0, n)[::-1])) 
... 
[] 
[0.0] 
[-0.57735, 0.57735] 
[-0.774597, 0.0, 0.774597] 
[-0.861136, -0.339981, 0.339981, 0.861136]

正如前面的答案所示，這個例子使用了多項式的泰勒展開式。

來源

2015-05-09 10:05:44

在python中查找勒讓德多項式的根

回答

相關問題