Трапециевидная интеграция в C

Question

Я пытаюсь вычислить интеграл от функции f (x) = (1-x ^ 2) ^ (1/2) от x = 0 до x = 1.Ответ должен быть примерно пи / 4.В настоящее время я получаю 2.

Моя текущая реализация правила трапеции выглядит следующим образом:

double
def_integral(double *f, double *x, int n)
{
  double F;
  for (int i = 0 ; i < n ; i++) {
    F += 0.5 * ( x[i+1] - x[i] ) * ( f[i] + f[i+1] );
  }
  return F;
}

Я создаю N делений для аппроксимации области под кривой между x_1 = 0 и x_N = 1 путем циклического перехода от i к N с x_i = i / N.

int
main(int argc, char **argv)
{
  int N = 1000;
  double f_x[N];
  double x[N];

  for (int i = 0 ; i <= N ; i++) {
    double x = i * 1. / N;
    f_x[i] = sqrt(1. - pow(x, 2.));
    //printf("%.2f %.5f\n", x, f_x[i]); //uncomment if you wanna see function values
  }

  double F_x = def_integral(f_x, x, N);

  printf("The integral is %g\n", F_x);
}

Результат 2, который я получаю в настоящее время, должен зависеть от числа делений N, однако, независимо от того, получу я N = 10000 или N = 100, я все равно получу 2.

Есть предложения?

Answer 1

В этом цикле for вы также забыли массив обновлений x.

for (int i = 0 ; i <= N ; i++) {
    double x = i * 1. / N;
    f_x[i] = sqrt(1. - pow(x, 2.));
    //printf("%.2f %.5f\n", x, f_x[i]); //uncomment if you wanna see function values
}

Итак, петлю for следует заменить на

for (int i = 0 ; i <= N ; i++) {
    double xi = i * 1. / N;
    x[i] = xi;
    f_x[i] = sqrt(1. - pow(xi , 2.));
    //printf("%.2f %.5f\n", x, f_x[i]); //uncomment if you wanna see function values
}

Answer 2

В вашем основном коде вы вызываете def_integral с двойным (x), и в функции ожидается массив x (double * x).Возможно (это то, что я предполагаю), проблема заключается в том, что вам нужна формула x (i + 1) -x (i) , но вы используете постоянный шаг.Действительно, x (i + 1) -x (i) = step_x является константой, поэтому вам не нужно каждое x (i), а только значение: 1./N
Другое замечание, с постоянным шагом ваша формула может бытьупрощено до:
F_x = step_x * ( 0,5 * f_x (x0) + f_x (x1) + ... + f_x (xn-1) + 0,5 * f_x (xn) ).Это помогает упростить код и написать более эффективный.Все прокомментировано в коде выше.Я надеюсь, что это может помочь вам.
С уважением.

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

double
def_integral(double *f, double step_x, int n)
{
  double F;
  for (int i = 0 ; i < n ; i++) {
    F += 0.5 * ( step_x ) * ( f[i] + f[i+1] );
  }
  return F;
}

int main()
{
  int N = 1001; // 1001 abscissas means 1000 intervalls (see comment on array size and indices)
  double f_x[N]; // not needed for the simplified algorithm
  double step_x = 1. / N; // x(i+1)-x(i) is constant

  for (int i = 0 ; i < N ; i++) { // Note : i<N and not i<=N
    double xi = i * step_x; // abscissa calculation
    f_x[i] = sqrt((1. - xi )*(1. + xi )); // cf chux comment
  }

  double F_x = def_integral(f_x, step_x, N);
  printf("The integral is %.10g\n", F_x);

// simplified algorithm
// F_x=step_x*( 0.5*f_x(x0)+f_x(x1)+...+f_x(xn-1)+0.5f_x(xn) )
  double xi;
  xi=0; // x(0)
  F_x=0.5*sqrt((1. - xi )*(1. + xi ));
  for (int i=1 ; i<=N-1 ; i++) {
    xi=step_x*i;
    F_x+=sqrt((1. - xi )*(1. + xi ));
  }
  xi=step_x*N;
  F_x+=0.5*sqrt((1. - xi )*(1. + xi ));
  F_x=step_x*F_x;
  printf("The integral is %.10g\n", F_x);

}