c - 从 0 到无穷大的数值积分

Question

我的目标是用 C 编程语言计算电子与氢原子核距离的概率分布函数 (PDF) 的数值积分。我已经编写了一个示例代码，但是由于我无法按我认为的必要增加限制，因此无法正确找到数值。我还包括了库，但我不能使用以下帖子中所述的值作为积分边界:min and max value of data type in C .这种情况下的补救措施是什么？也许应该切换到另一种编程语言？感谢任何帮助和建议，在此先感谢。

编辑:在一些值之后，我得到了错误段错误。我已经用 Wolframalpha 检查了积分的实际结果为 0.0372193。除此之外，如果我以较小的数量增加 k，我会因此得到零，这就是我定义 r[k]=k 的原因，我知道它应该更小以提高精度。

#include <stdio.h>
#include <math.h>
#include <limits.h>
#define a0 0.53 
int N = 200000;
// This value of N is the highest possible number in long double
// data format. Change its value  to adjust the precision of integration
// and computation time.
// The discrete integral may be defined as follows:
long double trapezoid(long double x[], long double f[]) {
  int i;
  long double dx = x[1]-x[0];
  long double sum = 0.5*(f[0]+f[N]);

    for (i = 1; i <  N; i++) 
        sum+=f[i];
  return sum*dx;  
}
main() {

    long double P[N], r[N], a;
    // Declare and initialize the loop variable
    int k = 0;
    for (k = 0; k <  N; k++)
    {
        r[k] = k ;
        P[k] = r[k] * r[k] * exp( -2*r[k] / a0);
        //printf("%.20Lf \n", r[k]);
        //printf("%.20Lf \n", P[k]);
    }
    a = trapezoid(r, P);    
    printf("%.20Lf \n", a);
}

最后的代码:

#include <stdio.h>
#include <math.h>
#include <limits.h>
#include <stdlib.h>
#define a0 0.53 
#define N LLONG_MAX
// This value of N is the highest possible number in long double
// data format. Change its value  to adjust the precision of integration
// and computation time.
// The discrete integral may be defined as follows:
long double trapezoid(long double x[],long double f[]) {
  int i;
  long double dx = x[1]-x[0];
  long double sum = 0.5*(f[0]+f[N]);

    for (i = 1; i <  N; i++) 
        sum+=f[i];
  return sum*dx;  
}
main() {
    printf("%Ld", LLONG_MAX);
    long double * P = malloc(N * sizeof(long double));
    long double * r = malloc(N * sizeof(long double));
    // Declare and initialize the loop variable
    int k = 0;
    long double integral;
    for (k = 1; k <  N; k++)
        {
        P[k] = r[k] * r[k] * expl( -2*r[k] / a0);
        }
    integral = trapezoid(r, P);
    printf("%Lf", integral);

}

编辑最后的代码工作:

#include <stdio.h>
#include <math.h>
#include <limits.h>
#include <stdlib.h>
#define a0 0.53 
#define N LONG_MAX/100
// This value of N is the highest possible number in long double
// data format. Change its value  to adjust the precision of integration
// and computation time.
// The discrete integral may be defined as follows:
long double trapezoid(long double x[],long double f[]) {
  int i;
  long double dx = x[1]-x[0];
  long double sum = 0.5*(f[0]+f[N]);

    for (i = 1; i <  N; i++) 
        sum+=f[i];
  return sum*dx;  
}
main() {
    printf("%Ld \n", LLONG_MAX);
    long double * P = malloc(N * sizeof(long double));
    long double * r = malloc(N * sizeof(long double));
    // Declare and initialize the loop variable
    int k = 0;
    long double integral;
    for (k = 1; k <  N; k++)
        {
        r[k] = k / 100000.0;
        P[k] = r[k] * r[k] * expl( -2*r[k] / a0);
        }
    integral = trapezoid(r, P);
    printf("%.15Lf \n", integral);
    free((void *)P);
    free((void *)r);
}

特别是我通过在除法运算中使用 float 来更改 r[k] 的定义，结果得到一个 long double 并且正如我在上一条评论中所述，我不能选择大于 Ns 的值LONG_MAX/100，我认为我应该进一步调查代码和 malloc 以解决问题。我找到了通过取极限分析获得的确切值；除了自己做之外，我已经用 TI-89 Titanium 和 Wolframalpha(在数值上和分析上)确认了结果。当间隔大小减小时，梯形法则效果很好。非常感谢这里所有的发帖者提出的想法。顺便说一句，LONG_MAX 的值为 2147483647 并不像我预期的那么大，极限不应该在 10 到 308 的幂左右吗？

Answer 1

数值观点

通常的梯形法不适用于不正确的积分。因此，高斯求积规则要好得多，因为它们不仅提供 2n-1 的精确性(即，对于 2n-1 次多项式，它们将返回正确的解)，而且还通过使用正确的权重函数来管理不正确的积分.

如果你的两边积分都不对，你应该试试Gauss-Hermite quadrature , 否则使用 Gauss-Laguerre quadrature .

“溢出”错误

long double P[N], r[N], a;

P大小约为 3MB，r 也是如此.内存太大了改为分配内存:

long double * P = malloc(N * sizeof(long double));
long double * r = malloc(N * sizeof(long double));

不要忘记包含 <stdlib.h>并使用 free在两个P和 r如果您不再需要它们。此外，您可能无法访问第 N 个条目，因此 f[N]是错误的。

使用高斯-拉盖尔积分

现在 Gauss-Laguerre 使用 exp(-x)作为权重函数。如果您不熟悉高斯求积:E(f) 的结果是 w * f 的积分, 其中w是权重函数。

你的 f看起来像这样，并且:

f x = x^2 * exp (-2 * x / a)

等一下。 f已经包含 exp(-term) , 所以我们可以用 t = x * a /2 代替 x并得到

f' x = (t * a/2)^2 * exp(-t) * a/2

自 exp(-t)已经是我们权重函数的一部分，您的函数现在完全适合 Gauss-Laguerre 正交。结果代码是

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

/* x[] and a[] taken from 
 * https://de.wikipedia.org/wiki/Gau%C3%9F-Quadratur#Gau.C3.9F-Laguerre-Integration
 * Calculating them by hand is a little bit cumbersome
*/
const int gauss_rule_length = 3;
const double gauss_x[] = {0.415774556783, 2.29428036028, 6.28994508294};
const double gauss_a[] = {0.711093009929, 0.278517733569, 0.0103892565016};

double f(double x){
    return x *.53/2 * x *.53/2 * .53/2;
}

int main(){
    int i;
    double sum = 0;
    for(i = 0; i < gauss_rule_length; ++i){
        sum += gauss_a[i] * f(gauss_x[i]);
    }
    printf("%.10lf\n",sum); /* 0.0372192500 */
    return 0;
}