1

我的目标是在 C 编程语言中计算电子到氢原子核的距离的概率分布函数 (PDF) 的数值积分。我已经编写了一个示例代码，但是由于我无法根据需要增加限制，因此无法正确找到数值。我还包含了该库，但我不能将以下帖子中所述的值用作整数边界：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);
}

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

4

1 回答 1

7

数值观点

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

如果你的积分两边都不合适，你应该尝试Gauss-Hermite 求积，否则使用Gauss-Laguerre 求积。

“溢出”错误

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>使用它们。此外，您可能无法访问第 N 个条目，所以是错误的。freePrf[N]

使用 Gauss-Laguerre 求积

现在 Gauss-Laguerreexp(-x)用作权重函数。如果您不熟悉高斯求积： 的结果E(f)是 的积分w * f，其中w是权重函数。

你f看起来像这样，并且：

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

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

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;
}
于 2013-11-01T10:11:38.173 回答