Nownuri
Nownuri

Reputation: 739

precision of trapezoidal rule in Python and C

I wrote integration codes in Python and again in C using the trapezoidal rule. The task was to read data from this file and integrate the velocities over time. After copying the data file into my working directory, I compiled and ran the following programs. The C program gave me 8.218921, while the Python programs gave 8.218924000000003. These two values are different from the 6th decimal point.

Where does the difference come from?

This is my C program.

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

int
main (void)
{

    FILE *fptr;
    if((fptr = fopen("velocities.txt", "r"))==NULL){
        printf("cannot open the file\n");
        exit(1);
    }

    printf("file opened\n");

    int t[200] = {};
    double v[200] = {};
    int n = 0;

    while( (fscanf(fptr, "%d\t %lf\n", &t[n], &v[n])!=EOF) && (n<200) ){
        printf("%d\t %lf\n", t[n], v[n]);
        n++;
    }

    int h = 0;
    float area = 0.0;
    printf("calculating the travel distance\n");
    for(int i=1;i<n;i++){
        h = t[i]-t[i-1];
        area += (v[i]+v[i-1])*h/2;
    }
    printf("area : %lf\n", area);


    fclose(fptr);

    return 0;
}

Here are the Python codes.

import numpy as np

arr = np.loadtxt("./velocities.txt")
arr = arr.T

h, area = 0, 0

for i in range(1, 101):
    h = arr[0,i] - arr[0,i-1]
    area += (arr[1,i]+arr[1,i-1])/2

print(area)

Upvotes: 0

Views: 147

Answers (2)

ForceBru
ForceBru

Reputation: 44878

The difference is because even though double v[200] = {}; has double precision, float area = 0.0; is a floating-point number, so half the precision is lost.

Python will use double by default and thus keep the precision.


As you can see, floating-point numbers quickly lose precision:

#include <stdio.h>

int main()
{
    float pi_flt = ((float)22)/7;
    double pi_dbl = ((double)22)/7;

    printf("%.8f, %.8lf\n", pi_flt, pi_dbl);
    return 0;
}

Compile and run:

$ clang -Weverything test.c
test.c:9:29: warning: implicit conversion increases floating-point precision: 'float'
      to 'double' [-Wdouble-promotion]
    printf("%.8f, %.8lf\n", pi_flt, pi_dbl);
    ~~~~~~                  ^~~~~~
1 warning generated.
$ ./a.out
3.14285707, 3.14285714

Upvotes: 1

MarcMush
MarcMush

Reputation: 1488

the result is probably correct, this is due to the precision with floating numbers.

This happens in pretty much every language, your C program probably just hides it here.

Upvotes: 0

Related Questions