递归计算Legendre多项式出错（段错误）[C++]求助

问题描述 #

使用C编写了基于递归公式计算N阶勒让德多项式及其导数的函数，运行时返回零值、错误数值或出现“segmentation fault”（段错误）。编译器为g，已尝试拆分多项式与导数的计算逻辑，但问题仍未解决。

原代码如下：

#include <iostream>
#include <cmath>
#include <iomanip>
#include <vector>
#include <functional>

using namespace std;

struct LegendreValues
{
  double function_value;
  double derivative_value;
};

LegendreValues calculateLegendrePolynomAndDerivative(const int&, const double&);
double calculateLegendrePolynom(const int&, const double&);
double calculateLegendreDerivative(const int&, const double&);

int main()
{
  int N;
  LegendreValues values;
  
  cout << "Unesite red Legendreovog polinoma: " << endl;
  cin >> N;
     
  cout << calculateLegendrePolynom(N, 0.5) << "\n" << calculateLegendreDerivative(N, 0.5) << endl;
}

LegendreValues calculateLegendrePolynomAndDerivative(const int& N, const double& x)
{
  LegendreValues results;
  if (N==0)
  {
    results.function_value = 1;
    results.derivative_value = 0;

  }
  
  else if (N==1)
  {
    results.function_value = x;
    results.derivative_value = 1;

  }

  else if (N>1)
  {
    results.function_value = (1/N) * ((2*N-1) * x * calculateLegendrePolynomAndDerivative(N-1, x).function_value - (N-1) * calculateLegendrePolynomAndDerivative(N-2, x).function_value);
    
    results.derivative_value = N/(1-x*x) * ( calculateLegendrePolynomAndDerivative(N-1, x).function_value - x * calculateLegendrePolynomAndDerivative(N, x).function_value);
    
  }

  return results;
}

double calculateLegendrePolynom(const int& N, const double& x)
{
  double fval;
  if (N==0)
  {
    fval =  1;
  }
  
  else if (N==1)
  {
    fval =  x;
  }
  
  else if (N>1)
  {
    fval = (1/N) * ((2*N-1) * x * calculateLegendrePolynom(N-1, x) - (N-1) * calculateLegendrePolynom(N-2, x));
  }
  return fval;
}

double calculateLegendreDerivative(const int& N, const double& x)
{
  double dfval;
  if (N==0)
  {
    dfval =  0;
  }
  
  else if (N==1)
  {
    dfval =  1;
  }
  
  else if (N>1)
  {
    dfval = N/(1-x*x) * ( calculateLegendreDerivative(N-1, x) - x * calculateLegendreDerivative(N, x));
  }
  return dfval;
}

使用的递归公式：
勒让德多项式递推公式：$nP_n(x) = (2n-1)xP_{n-1}(x) - (n-1)P_{n-2}(x)$
导数递推公式：$(1-x^2)P'n(x) = n[P{n-1}(x) - xP_n(x)]$

错误分析 #

1. 整数除法导致计算结果为0
   代码中(1/N)是整数除法，当N≥1时，1除以整数N的结果为0，直接导致多项式的计算结果全部错误。

2. 导数函数的无限递归
   在calculateLegendreDerivative函数中，当N>1时的公式调用了calculateLegendreDerivative(N, x)，这会触发无限递归，最终导致栈溢出，出现“segmentation fault”错误。同样，calculateLegendrePolynomAndDerivative中的导数计算也调用了自身的N阶计算，引发相同问题。

3. 未处理x=±1的除零情况
   当x=1或x=-1时，1-x*x的值为0，会触发除零错误，导致程序异常。

修正方案与代码 #

关键修正点 #

• 将整数除法改为浮点数除法（如1.0/N）。
• 修正导数计算逻辑，避免无限递归：利用已计算的多项式值来计算导数，而非递归调用导数函数。
• 添加x=±1的特殊处理，避免除零错误。
• 优化递归逻辑，避免重复计算（可选：改用迭代递推，提升效率并减少栈开销）。

修正后的代码 #

#include <iostream>
#include <cmath>
#include <iomanip>

using namespace std;

struct LegendreValues
{
    double function_value;
    double derivative_value;
};

// 计算多项式值
double calculateLegendrePolynom(const int& N, const double& x)
{
    if (N == 0)
        return 1.0;
    if (N == 1)
        return x;
    // 改用浮点数除法
    return ((2.0 * N - 1.0) * x * calculateLegendrePolynom(N - 1, x) - (N - 1.0) * calculateLegendrePolynom(N - 2, x)) / N;
}

// 计算导数值
double calculateLegendreDerivative(const int& N, const double& x)
{
    const double eps = 1e-10;
    // 处理x=1的情况
    if (fabs(x - 1.0) < eps)
        return N * (N + 1.0) / 2.0;
    // 处理x=-1的情况
    if (fabs(x + 1.0) < eps)
        return pow(-1.0, N+1) * N * (N + 1.0) / 2.0;
    
    if (N == 0)
        return 0.0;
    if (N == 1)
        return 1.0;
    // 利用多项式值计算导数，避免递归调用导数函数
    double Pn = calculateLegendrePolynom(N, x);
    double Pn_1 = calculateLegendrePolynom(N-1, x);
    return N * (Pn_1 - x * Pn) / (1.0 - x*x);
}

// 同时计算多项式和导数（优化：避免重复计算Pn和Pn-1）
LegendreValues calculateLegendrePolynomAndDerivative(const int& N, const double& x)
{
    LegendreValues results;
    const double eps = 1e-10;
    
    if (N == 0)
    {
        results.function_value = 1.0;
        results.derivative_value = 0.0;
        return results;
    }
    if (N == 1)
    {
        results.function_value = x;
        results.derivative_value = 1.0;
        return results;
    }
    
    // 先计算低阶多项式，避免重复计算
    double Pn_2 = calculateLegendrePolynom(N-2, x);
    double Pn_1 = calculateLegendrePolynom(N-1, x);
    results.function_value = ((2.0 * N - 1.0) * x * Pn_1 - (N - 1.0) * Pn_2) / N;
    
    // 处理x=±1的情况
    if (fabs(x*x - 1.0) < eps)
    {
        if (x > 0)
            results.derivative_value = N * (N + 1.0) / 2.0;
        else
            results.derivative_value = pow(-1.0, N+1) * N * (N + 1.0) / 2.0;
    }
    else
    {
        results.derivative_value = N * (Pn_1 - x * results.function_value) / (1.0 - x*x);
    }
    
    return results;
}

int main()
{
    int N;
    cout << "请输入勒让德多项式的阶数: " << endl;
    cin >> N;
    if (N < 0)
    {
        cout << "阶数不能为负数" << endl;
        return 1;
    }
    
    double x = 0.5;
    // 使用优化后的同时计算函数，减少重复计算
    LegendreValues values = calculateLegendrePolynomAndDerivative(N, x);
    cout << fixed << setprecision(8);
    cout << "多项式值: " << values.function_value << endl;
    cout << "导数值: " << values.derivative_value << endl;
    
    return 0;
}

额外优化建议 #

递归计算勒让德多项式会存在大量重复计算（比如计算Pn时会重复计算Pn-1、Pn-2多次），对于高阶N（如N>20），建议改用迭代递推的方式，从P0、P1开始逐步计算到Pn，同时记录中间值，既能提升计算效率，也能避免递归栈溢出的问题。

内容的提问来源于stack exchange，提问作者Josip Gregorić