数值积分,使用 C 中的梯形法则

Numerical Integration, using the trapezium rule in C

我正在使用梯形法则来计算函数在 0 和无穷大之间的积分。我可以计算给定 N 值的积分值,现在我试图将 N 从 2 循环到给定值,但它不起作用。它不断计算 N 为 2 时的积分值并重复而不是 N 的新值。我认为问题出在 main() 的 for 循环中。

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

double f(double x) {
 double a;
  a =1/((1+x)*pow(x,0.5));
  return a;}

double tra(double upper, double lower, int N) {
 double sum, step, integral,lowest;
  step=(upper-lower)/(N-1);
  lower=lower+step;

    if(lower==0) {
       lowest=DBL_EPSILON;}
    else {
    lowest=lower;}

    while(lower<upper) {
         sum=sum+f(lower);
         lower=lower+step;}

 integral=step*(sum+(f(upper)/2)+(f(lowest)/2));
 sum=0;
 return integral;}


main() {
int N;
double upper=DBL_EPSILON*2, lower=0, total=0;
for(N=2;N<20000;N+=100) {   /*Here im trying to loop N so that the        integral is calculated for increasing values of N*/

   while(upper<FLT_MAX) {
        total=total+tra(upper, lower, N);
        lower=upper;
        upper=upper*2;}
        printf("Integral is %.10f\n", total);
}
}

我建议您将变量初始化移动到 for 循环内,如下所示:

int main(void) {
    int N;
    double upper, lower, total;
    for(N=2;N<20000;N+=100) {
        upper = DBL_EPSILON*2;
        lower = 0;
        total = 0;
        while(upper<FLT_MAX) {
            total=total+tra(upper, lower, N);
            lower=upper;
            upper=upper*2;
        }
        printf("Integral is %.10f\n", total);
    }
    return 0;
}