我应该使用多少个正交点？

计算科学 正交 一体化 gsl
2021-12-26 23:14:57

我正在尝试计算以下积分

∫0∞e−yya/2Lcb(y)Led(y)dy
使用 GNU 科学图书馆中的广义 Gauss-Laguerre 正交例程。这里L是广义拉盖尔函数，并且a>−1,b>0,d>0. 我不确定如何选择正交点的数量，不同的数字给出不同的值。

有什么建议？

（编辑：MWE）

#include <stdio.h>
#include <stdlib.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_sf_laguerre.h>

struct data { double b; int c; double d; int e; };

double f(double y, void* userdata) {
    struct data *d = (struct data *) userdata;
    return gsl_sf_laguerre_n(d->c, d->b, y) * gsl_sf_laguerre_n(d->e, d->d, y);
}

int main() {
    const double a = -0.5;
    const double b = 0.5;
    const int c = 0;
    const double d = 0.5;
    const int e = 0;
    const size_t num_nodes = (c + e + 1) / 2;

    const gsl_integration_fixed_type *T = gsl_integration_fixed_laguerre;
    gsl_integration_fixed_workspace *w
        = gsl_integration_fixed_alloc(T, num_nodes, 0.0, 1.0, 0.5 * a, 0.0);

    struct data params = { b, c, d, e };
    gsl_function F;
    F.function = &f;
    F.params = &params;

    double result;
    gsl_integration_fixed(&F, &result, w);
    printf("%12.e\n", result);

    gsl_integration_fixed_free(w);

}
2个回答

以下程序会将您的积分与 GSL 集成。我无法重现您在 Mathematica 中引用的结果。请注意，当c=e=0，拉盖尔多项式总是1对于任何值x. 下面的代码做了一个稍微有趣的例子c=2,e=3. 输出是：

result = 1.539250228072e+00

Mathematica 的输出是1.539250228072043. 什么时候c=e=0正如您最初发布的那样，GSL 和 Mathematica 都会产生结果1.225416702465178.

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#include <gsl/gsl_math.h>
#include <gsl/gsl_specfunc.h>
#include <gsl/gsl_integration.h>

struct data
{
  double b;
  int c;
  double d;
  int e;
};

double
f(double x, void * params)
{
  struct data *d = (struct data *) params;
  double Lcb = gsl_sf_laguerre_n(d->c, d->b, x);
  double Led = gsl_sf_laguerre_n(d->e, d->d, x);
  return Lcb * Led;
}

int
main()
{
  const gsl_integration_fixed_type * T = gsl_integration_fixed_laguerre;
  const double a = -0.5;
  const double b = 0.5;
  const int c = 2;
  const double d = 0.5;
  const int e = 3;
  const size_t n = (c + e + 1) / 2;
  gsl_integration_fixed_workspace * w =
    gsl_integration_fixed_alloc(T, n, 0.0, 1.0, 0.5 * a, 0.0);
  gsl_function F;
  double result;
  struct data data_params = { b, c, d, e };

  F.function = &f;
  F.params = &data_params;

  gsl_integration_fixed(&F, &result, w);
  fprintf(stderr, "result = %.12e\n", result);

  gsl_integration_fixed_free(w);

  return 0;
}

printf 格式的工作方式，像这样的调用

printf("%12.e\n", result);

指示 printf 打印一个包含 12 位数字/字符/空格（宽度为12）的数字，小数点后没有数字，就像在精度说明符中一样"12."。

将其更改为

printf("%.12e\n", result);

导致打印1.225416702465e+00，这可能是预期的，并且与 vibe 的答案相匹配。

请参阅https://en.cppreference.com/w/cpp/io/c/fprintf：

（可选的） 。后跟整数或 *，或者两者都不指定转换的精度。在使用 * 的情况下，精度由 int 类型的附加参数指定。如果此参数的值为负，则将其忽略。如果既没有使用数字也没有使用 *，则精度为零。有关精度的确切影响，请参见下表。[强调补充]

其它你可能感兴趣的问题
上一篇SciPy 3d 绘图积分∫X是的dX∫xydx为了是的y在[ - 4 , 4 ][−4,4] 下一篇使用牛顿法和二次插值求根 - 这本教科书有错误吗？