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

Problems with GSL using Legendre polynomials

本文关键字：GSL 问题 多项式      更新时间：2023-10-16

我正在尝试更新使用具有已弃用函数的 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
}
}
}
}
}

当我使用 GSL1.9 编译时，我使用已弃用的函数gsl_sf_legendre_sphPlm_array，而当我使用 GSL 2.5 计算时，我使用新函数gsl_sf_legendre_array，它明确需要一个关键字进行规范化 (GSL_SF_LEGENDRE_SPHARM)，并且不要求参数m。旧版本给了我一致的结果，而新版本在 25 吨循环后返回了分割错误。你们中的任何人都可以帮我纠正代码并解释我做错了什么吗？

我尝试使用调试符号(以下命令中的-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 行。原因可能是您尝试越界访问数组。也许看看如何在 gsl 勒让德多项式手册中"正确"分配内存。我正在特别考虑gsl_sf_legendre_array_n和gsl_sf_legendre_array_index的功能。

注意：我自己没有使用过GSL的这一部分，但我通常发现在使用GSL时，使用在调用函数之前/之后打包/解压缩的临时变量和数组很有用。不漂亮，但不太容易出错，并且由于您使用的是 C 的"plusspluss"版本，因此您始终可以将实现包装在class中。 希望它对您有所帮助。

编辑：

我尝试增加legendre数组的大小，当大小设置为以下值时，程序完成：

double legendre[ntheta + 10][order + 4];

也许知道这些功能如何工作的人可以回答为什么......

相关文章：