I'm trying to implement numerical integration using the trapezoidal approximation using this formula :
My problem is I don't get how to implement this correctly. To test I wrote a file with 22050
double values all equal to 2 like :
....................
value =2.0;
for ( index = 0 ; index < 22050;index++){
fwrite(&value,sizeof(double),1,inp2);
}
to keep the question simple, say I want to the integral value of each 100 samples:
X Area integral value
0-100 should be 200
100-200 should be 200
..... ...........
22000-22050 should be 100
to do that I 've wrote a program that should do that but the result that get is 4387950
for 100 samples
here is my code :
..............................
// opening the files
double* inputData= NULL;
unsigned int N = 100;
double h= 0.0;
unsigned int index= 0;
FILE* inputFile=NULL;
double value =0.0;
int i =0,j=0;
inputFile = fopen("sinusD","rb");
outputFile=fopen("Trapez","wb+");
if( inputFile==NULL || outputFile==NULL){
printf("Couldn't open the files \n");
return -1;
}
inputData = (double*) malloc(sizeof(double)*N);
h=22050/2;
while((i = fread(inputData,sizeof(double),N,inputFile))==N){
value += inputData[0] +inputData[N];
for(index=1;index<N;index++){
value+=(2*inputData[index]);
}
value *=h;
fprintf(outputFile,"%lf ",value);
value =0;
}
if(i!=0){
value = 0;
i=-i;
printf("i value %i\n", i);
fseek(inputFile,i*sizeof(double),SEEK_END);
fread(inputData,sizeof(double),i,inputFile);
for(index=0;index<-i;index++){
printf("index %d\n",index);
value += inputData[0] +inputData[i];
value+=(2*inputData[index]);
}
value *=h;
fprintf(outputFile,"%lf ",value);
value =0;
}
fclose(inputFile);
fclose(outputFile);
free(inputData);
return 0;}
any idea how to do that ?
UPDATE
while((i = fread(inputData,sizeof(double),N,inputFile))==N){
value = (inputData[0] + inputData[N])/2.0;
for(index=1;index<N;index++){
value+=inputData[index];
}
value *=h;
fprintf(outputFile,"%lf ",value);
printf(" value %lf\n",value);
value =0;
}
I get 199.000
as a result for each segment .
Why you didn't start with something simple. Let's say you have the following data {1,2,3,4,5,6,7,8,9,10}
and assume h = 1
. This is simple,
#include <stdio.h>
#define SIZE 10
int main()
{
double a[SIZE] = {1,2,3,4,5,6,7,8,9,10}, sum = 0.0, trapz;
int h = 1;
int i = 0;
for ( i; i < SIZE; ++i){
if ( i == 0 || i == SIZE-1 ) // for the first and last elements
sum += a[i]/2;
else
sum += a[i]; // the rest of data
}
trapz = sum*h; // the result
printf("Result: %f \n", trapz);
return 0;
}
This is the result
Result: 49.500000
Double check your work with Matlab:
Y = [1 2 3 4 5 6 7 8 9 10];
Q = trapz(Y)
Q =
49.5000
$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
Edit: For your question in the comment: This is the matlab code:
X = 0:pi/100:pi; % --> h = pi/100
Y = sin(X); % get values as many as the size of X
Q = trapz(X,Y);
Q =
1.9998
Now to fulfil same scenario in C, do the following
#include <stdio.h>
#include <math.h>
#define SIZE 101
#define PI 3.14159265358
int main()
{
double X[SIZE], Y[SIZE], incr = 0.0, h = PI/100.0, sum = 0.0, trapz;
int i = 0, k = 0, j = 0;
// Generate samples
for ( i; i < SIZE; ++i)
{
X[i] = incr;
incr += h;
}
// Generate the function Y = sin(X)
for ( k; k < SIZE; ++k)
{
Y[k] = sin(X[k]);
}
// Compute the integral of sin(X) using Trapezoidal numerical integration method
for ( j; j < SIZE; ++j){
if ( j == 0 || j == SIZE-1 ) // for the first and last elements
sum += Y[j]/2;
else
sum += Y[j]; // the rest of data
}
trapz = sum * h; // compute the integral
printf("Result: %f \n", trapz);
return 0;
}
The result is
Result: 1.999836