求勒让德多项式的递归程序

Recursive program to obtain Legendre Polynomials

我正在尝试在 C++ 上为 Legendre 多项式的构造函数编写一个函数,该函数打印出次数为 m 的多项式的系数。多项式遵循 simple recursion relation.

现在,我正在尝试实现该关系,并且它对每个 n 都完美工作,最多 8 个,从 9 开始,在第 7 次迭代周围,它会拾取一个垃圾值甚至不在任何先前的系数向量中。我想知道我该如何解决这个问题。我向您展示我的代码:

#include <vector>
#include <cmath>
std::vector<double> set_coeffs(int m){
std::vector<double> coeffs;

if (m == 0) //Casos iniciales para empezar la recursión
{               
    coeffs.push_back(1);
} else if (m == 1)
{
    coeffs.push_back(1);
    coeffs.push_back(0);
} else if (m == 2) //Puse también el caso 2 porque de otro modo ocurre el mismo problema pero desde n=5
{
    coeffs.push_back(1.5);
    coeffs.push_back(0);
    coeffs.push_back(-0.5);
} 
 else
{
    std::vector<double> v = set_coeffs(m-1);
    std::vector<double> u = set_coeffs(m-2);
    std::cout << "inicia cicle" << std::endl;

    double a = (2* ((double)m) -1)/((double)m);
    double b = (((double)m)-1)/((double)m);

    coeffs.push_back(a*v[0]);
    coeffs.push_back(a*v[1]);

    for (int i = 0; i < m-1; i++)
    {
        double c = a*v[i+2] - b*u[i];
        std::cout << m << " " << v[i+2] << " " << u[i] << " " << c <<std::endl;
        coeffs.push_back(c);
    }

    std::cout << "termina cicle" << std::endl;
} 
return coeffs;
}

正在调用 set_coeffs(n) returns 具有 n+1 个元素的向量。

因此,std::vector<double> v = set_coeffs(m-1); 有 m 个元素。

在您的 for 循环中,i0 运行到 m-2,然后您访问 v[i+2]。在最后一次迭代中,这将访问超出范围的 v[m]

interjay's 中所述,发布的代码没有正确考虑向量大小的差异。

为了使索引匹配更容易(至少恕我直言,但它也应该有助于涉及这些系数的进一步计算),我将以相反的顺序存储系数。我也会利用它们的特定模式(只有奇数 偶数不为零)。

可能值得一提的是,当应用于大向量时,这种递归方法效率低下。

也就是说,这是一个可能的实现:

std::vector<double> Legendre_coefficients(int m)
{
    if (m == 0) 
    {               
        return {1};
    }
    if (m == 1)
    {
        return {0, 1};
    }

    // Initialize with zero, only at most (half + 1) of the terms will be changed later
    std::vector<double> coeffs(m + 1);

    // Consider some form of memoization instead of this recursion
    std::vector<double> v = Legendre_coefficients(m - 1);
    std::vector<double> u = Legendre_coefficients(m - 2);

    // using literals of floating point type, 'm' is promoted by the compiler
    double a = (2.0 * m - 1.0) / m;
    double b = (m - 1.0) / m;

    int first = 1;
    // If 'm' is even, so is (m - 2) and its first element is zero. It can be skipped.
    // It also avoids annoying '-0' in the output 
    if ( m % 2 == 0 )
    {
        coeffs[0] = -b * u[0];
        first = 2;
    }
    for (int i = first; i < m - 1; i += 2)
    {
        coeffs[i] = (a * v[i - 1] - b * u[i]);
    }
    coeffs[m] = a * v[m - 1];

    return coeffs;
}

可测试here