GSL 使用勒让德多项式的问题

Problems with GSL using Legendre polynomials

我正在尝试更新旧代码,该代码使用的是带有已弃用函数的 GSL 版本,但我在寻找如何使用新版本的归一化勒让德多项式函数时遇到了麻烦。这是总结问题的片段:

#include <iostream>
#include <gsl/gsl_sf_legendre.h>
#include <cmath>

#define GSL_NEW

using namespace std;

int main() {

  int order = 17;
  int ntheta = 36;
  double theta_step = M_PI / ntheta;
  double c, theta;
  double legendre[ntheta][order+1];

  for( int m = 0; m <= order; m += 2) {
    for(int l = m; l <= ntheta; l += 2 ) {
      for( int t = 0; t < ntheta; t++ ) {
        theta = ( ntheta + 0.5 ) * theta_step;
        c = cos(theta);

        if( l == m ) {
#ifdef GSL_NEW
          gsl_sf_legendre_array( GSL_SF_LEGENDRE_SPHARM, order, c, &legendre[t][l] );
          cout << legendre[t][l] << endl;
#else
          gsl_sf_legendre_sphPlm_array(order, m, c, &legendre[t][l] );
          cout << legendre[t][l] << endl;
#endif
        }
      }
    }
  }
}

当我使用 GSL 1.9 编译时,我使用了已弃用的函数 gsl_sf_legendre_sphPlm_array,而当我使用 GSL 2.5 进行计算时,我使用了新函数 gsl_sf_legendre_array,它明确要求规范化关键字(GSL_SF_LEGENDRE_SPHARM) 并且不要求参数 m。旧版本给了我一致的结果,而新版本 returns 在 25 t 循环后出现分段错误。你们中的任何人都可以帮我更正代码并解释我做错了什么吗?

我尝试使用调试符号进行编译(下面命令中的 -g 标志假设您的程序名为 "main.cpp")

usr@cmptr $ g++ -g main.cpp -lgsl -lgslcblas -lm -o main

并且 运行 使用调试器的程序,例如gdb(GNU 调试器):

usr@cmptr $ gdb main
(gdb) run main
Starting program: /home/usr/Desktop/main 
0, 0, 0.282095
1, 0, 0.282095
2, 0, 0.282095
3, 0, 0.282095
4, 0, 0.282095
5, 0, 0.282095
6, 0, 0.282095
7, 0, 0.282095
8, 0, 0.282095
9, 0, 0.282095
10, 0, 0.282095
11, 0, 0.282095
12, 0, 0.282095
13, 0, 0.282095
14, 0, 0.282095
15, 0, 0.282095
16, 0, 0.282095
17, 0, 0.282095
18, 0, 0.282095
19, 0, 0.282095
20, 0, 0.282095
21, 0, 0.282095
22, 0, 0.282095
23, 0, 0.282095
24, 0, 0.282095

Program received signal SIGSEGV, Segmentation fault.
0x0000555555554ded in main () at main.cpp:26
26            cout << t << ", " << l << ", " << legendre[t][l] << endl;
(gdb) quit

错误在您写入屏幕的第 26 行。原因可能是您试图越界访问数组。也许看看如何在 manual for gsl Legendre polynomials 中分配内存 "properly"。我特别想到函数 gsl_sf_legendre_array_ngsl_sf_legendre_array_index

Note: I have not used this part of GSL myself, but I usually find that when using GSL it is useful to use temporary variables and arrays that you pack/unpack before/after calling functions. Not pretty but less error prone, and since you are using the "plusspluss" version of C you can always wrap the implementation in a class. Hopefully it is helpful.


编辑:

我尝试增加 legendre 数组的大小,当大小设置为:

时程序完成
double legendre[ntheta + 10][order + 4];

也许知道函数如何工作的人可以回答为什么...