返回勒让德多项式系数的函数不够准确

Question

我正在尝试编写一个返回第 n 个勒让德多项式的系数的函数。

// store previous polynomials to make the recursive process faster
unordered_map<int, vector<double>> legendre_cache; 

// Recursive function to compute the coefficients of the nth Legendre polynomial
vector<double> legendre_polynomial(int n) { // uses Bonnet's recursion formula P_{n+1} = ( (2n+1) x P_n - n P_{n-1} )/(n+1)
    // Check if the result is already in the cache
    if (legendre_cache.find(n) != legendre_cache.end()) {
        return legendre_cache[n];
    }

    // Base cases
    if (n == 0) {
        return { 1.0 }; // P0(x) = 1
    }
    if (n == 1) {
        return { 1.0, 0.0 }; // P1(x) = x
    }

    // Get coefficients for P_(n-1) and P_(n-2)
    vector<double> Pn_minus_1 = legendre_polynomial(n - 1);
    vector<double> Pn_minus_2 = legendre_polynomial(n - 2);

        // Get the sizes of the vectors to n+1
    Pn_minus_1.push_back(0.0); // P_{n-1} gets multiplied by x
    Pn_minus_2.insert(Pn_minus_2.begin(), 0.0); // P_{n-2} needs two leading digits 
    Pn_minus_2.insert(Pn_minus_2.begin(), 0.0);

    // Calculate coefficients for P_n
    vector<double> Pn(n + 1, 0.0); // Initialize Pn with n+1 zeros

    // Using the recurrence relation to fill coefficients
    for (int k = 0; k <= n; k++) {
        Pn[k] += ((2 * n - 1) * Pn_minus_1[k] - (n - 1) * Pn_minus_2[k]) / n;
    }

    // Store the result in the cache before returning
    legendre_cache[n] = Pn;
    return Pn;
}

该算法有效，因为我将其结果与 n = 10 以内的其他来源进行了比较。

然而，这个函数很快就会变得不准确，而且它的递归性质并没有给它带来好处：

在这里，我评估我的多项式并将其与“boost/math”库中的类似函数进行比较：

int main() {
    cout << fixed << setprecision(15);
    // Example inputs
    int n = 48;
    double x = 0.9;

    cout << "Result of my function:" << endl;
    vector<double> coeffs = legendre_polynomial(n);
    double result = 0;
    int l = coeffs.size();
    for (int i = l - 1; i >= 0; i--) {
        result += coeffs[l - i - 1] * pow(x, i);
    }
    cout << result << endl;

    cout << "Result of boost library:" << endl;
    cout << boost::math::legendre_p(n, x) << endl;
    return 0;
}

n 越大，我的函数就越不准确，直到 n = 45 左右时完全崩溃：

测试代码的输出：

Result of my function:
0.151876441197771
Result of boost library:
-0.106892037065314

---

我从 boost 网站 知道 boost 函数使用与我相同的递归方法。

因此，假设他们的功能有效，那么也应该可以改进我的功能。

感谢您的帮助。

---

编辑1：

需要明确的是，我评估多项式只是为了证明我的系数是错误的；我主要需要多项式的系数（勒让德高斯求积）。

据我了解库中的函数，我认为那里从未提到过系数；递归公式仅应用于特定的 x 值。 坦率地说，我还没有达到重写代码来提取系数的水平。

---

编辑2：

我所说的“该算法是有效的，因为我将它的结果与其他来源的结果进行了比较，最多为n = 10”，我怀疑如果没有准确性错误，该函数将会起作用。我认为由于舍入、除法、浮点导致的错误会因为递归而像滚雪球一样不断增加，直到函数变得不可用。 （我的问题是如何处理这个问题）

Answer 1

我认为您无法像您一样计算勒让德多项式在给定点的值，即直接使用多项式的系数。事实上，您将面临数值错误，因为当多项式的阶数很大时，

double

无法处理

x

与“巨大”系数的乘法。

相反，您可以直接依赖这些多项式的三项递推关系来计算给定点的值：

#include <iostream>
#include <iomanip>
#include <boost/math/special_functions/legendre.hpp>

int main(int argc, char* argv[]) 
{
    std::cout << std::fixed << std::setprecision(15);

    int    n = argc>=2 ? atoi(argv[1]) : 48;
    double x = 0.9;
    
    double p0=1;
    double p1=x;
    for (size_t k=1; k<n; k++)  
    {
       double tmp = ( (2*k+1)*x*p1 - k*p0) / (k+1);
       p0=p1;
       p1=tmp;
   }

   std::cout << "Result of my function:" << std::endl;
   std::cout << p1 << std::endl;

   std::cout << "Result of boost library:"    << std::endl;
   std::cout << boost::math::legendre_p(n, x) << std::endl;
    
   return 0;
}

可能的输出：

Result of my function:
-0.106892037065314
Result of boost library:
-0.106892037065314