how to numerically integrate a variable that is being calculate in the program as a pointer (using e.g. trapezoidal rule) in c language

I have a code, that was not made by me. In this complex code, many rules are being applied to calculate a quantity, d(x). in the code is being used a pointer to calculate it.

I want to calculate an integral over this, like: W= Int_0 ^L d(x) dx ?

I am doing this:

#define DX 0.003
void WORK(double *d, double *W)
{
    double INTE5=0.0;
    int N_X_POINTS=333;
    double h=((d[N_X_POINTS]-d[0])/N_X_POINTS); 
    W[0]=W[0]+((h/2)*(d[1]+2.0*d[0]+d[N_X_POINTS-1])); /*BC*/
    for (i=1;i<N_X_POINTS-1;i++) 
    {
        W[i]=W[i]+((h/2)*(d[0]+2*d[i]+d[N_X_POINTS]))*DX;
        INTE5+=W[i];
    }
    W[N_X_POINTS-1]=W[N_X_POINTS-1]+((h/2)*(d[0]+2.0*d[N_X_POINTS-1]+d[N_X_POINTS-2])); /*BC*/
}

And I am getting "Segmentation fault". I was wondering to know if, I am doing 开发者_运维百科right in calculate W as a pointer, or should declare it as a simple double? I guess the Segmentation fault is coming for this.

Other point, am I using correctly the trapezoidal rule?

Any help/tip, will very much appreciate.

Luiz

I don't know where that code come from, but it is a lot ugly and has some limits hard-encoded (333 points and increment by 0.003). To use it you need to "sample" properly your function and generate pairs (x, f(x))...

A possible clearer solution to your problem is here.

Let us consider you function and let us suppose it works (I believe it does't, it's a really obscure code...; e.g. when you integrate a function, you expect a number as result; where's this number? Maybe INTE5? It is not given back... and if it is so, why the final update of the W array? It's useless, or maybe we have something meaningful into W?). How does would you use it?

The prototype

  void WORK(double *d, double *W);

means the WORK wants two pointers. What these pointers must be depends on the code; a look at it suggests that indeed you need two arrays, with N_X_POINTS elements each. The code reads from and writes into array W, and reads only from d. The N_X_POINTS int is 333, so you need to pass to the function arrays of at least 333 doubles:

   double d[333];
   double W[333];

Then you have to fill them properly. I thought you need to fill them with (x, f(x)), sampling the function with a proper step. But of course this makes no too much sense. Already said that the code is obscure (now I don't want to try to reverse engineering the intention of the coder...).

Anyway, if you call it with WORK(d, W), you won't get seg fault, since the arrays are big enough. The result will be wrong, but this is harder to track (again, sorry, no "reverse engineering" for it).

Final note (from comments too): if you have double a[N], then a has type double *.

A segmentation fault error often happens in C when you try to access some part of memory that you shouldn't be accessing. I suspect that the expression d[N_X_POINTS] is the culprit (because arrays in C are zero-indexed), but without seeing the definition of d I can't be sure.

Try putting informative printf debugging statements before/after each line of code in your function so you can narrow down the possible sources of the problem.

Here's a simple program that integrates $f(x) = x^2$ over the range [0..10]. It should send you in the right direction.

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

double int_trapezium(double f[], double dX, int n)
{
   int i;
   double sum;

   sum = (f[0] + f[n-1])/2.0;
   for (i = 1; i < n-1; i++)
       sum += f[i];
   return dX*sum;
}

#define N 1000
int main()
{
    int i;
    double x;
    double from = 0.0;
    double to = 10.0;
    double dX = (to-from)/(N-1);
    double *f = malloc(N*sizeof(*f));

    for (i=0; i<N; i++)
    {
        x = from + i*dX*(to-from);
        f[i] = x*x;
    }
    printf("%f\n", int_trapezium(f, dX, N));
    free(f);
    return 0;
}