1

I have a code that tries to solve an integral of a function in a given interval numerically, using the method of Trapezoidal Rule (see the formula in Trapezoid method ), now, for the function sin(x) in the interval [-pi/2.0,pi/2.0], the integral is waited to be zero.

In this case, I take the number of partitions 'n' equal to 4. The problem is that when I have pi with 20 decimal places it is zero, with 14 decimal places it is 8.72e^(-17), then with 11 decimal places, it is zero, with 8 decimal places it is 8.72e^(-17), with 3 decimal places it is zero. I mean, the integral is zero or a number near zero for different approximations of pi, but it doesn't have a clear trend.

I would appreciate your help in understanding why this happens. (I did run it in Dev-C++).

#include <iostream>
#include <math.h>
using namespace std;
#define pi 3.14159265358979323846
//Pi: 3.14159265358979323846
double func(double x){
    return sin(x);
    }

int main() {
    double x0 = -pi/2.0, xf = pi/2.0;
    int n = 4;
    double delta_x = (xf-x0)/(n*1.0);
    double sum = (func(x0)+func(xf))/2.0;
    
    double integral;
    
    for (int k = 1; k<n; k++){      
       // cout<<"func: "<<func(x0+(k*delta_x))<<" "<<"last sum: "<<sum<<endl;
        sum = sum + func(x0+(k*delta_x));     
       // cout<<"func + last sum= "<<sum<<endl;
    }
    integral = delta_x*sum;
    cout<<"The value for the integral is: "<<integral<<endl; 
    return 0;
}
Spherk
  • 59
  • 8
  • 1
    [`M_PI`](https://stackoverflow.com/questions/1727881/how-to-use-the-pi-constant-in-c) is a thing. – tadman Mar 22 '21 at 18:46
  • 2
    @tadman or the more modern [`std::numbers::pi_v`](https://en.cppreference.com/w/cpp/numeric/constants) if available. – Brian61354270 Mar 22 '21 at 18:47
  • 1
    Are you asking why using different values of Pi results in different conclusions? Yeah, it will likely do that. – tadman Mar 22 '21 at 18:47
  • 2
    @Brian Given that it uses MinGW, might be ten years before that happens. – tadman Mar 22 '21 at 18:48
  • What I say is that one can expect to get a more exact integral (in this case nearer to zero) when using a more exact value of pi, but this is not the case, because, for example, I get zero with pi= 3.141, and something different for pi defined with 14 decimal places. – Spherk Mar 22 '21 at 18:54
  • The fact it uses `double` or any sort of floating number to store the values, it will be inaccurate. – Ranoiaetep Mar 22 '21 at 19:01
  • 1. The fact the you get something that close to the exact result with only 4 partitions is because you are "lucky" that rather large errors exactly cancel each other in this case so we are down to limitations inherent to the `double` representation. 2. The `double` type has about 17 significant digits hence almost any operation will give a rounding error of relative magnitude around 1e-17. The cumulative effect of these rounding errors gives the difference to the exact result. It is expected that there is no trend. – nielsen Mar 22 '21 at 19:39
  • @nielsen You are right about the cumulative rounding effect. But when a change the interval to, x0= pi/2.0 and xf=(3pi)/2.0, where sin(x) is symmetric with x-axes, the more exact the value of pi, the more near to zero is the result, so there is a trend. – Spherk Mar 22 '21 at 21:42

1 Answers1

1

OP is integrating y=sin(x) from -a to +a. The various tests use different values of a, all near pi/2.

The approach uses a linear summation of values near -1.0, down to 0 and then up to near 1.0.

This summation is sensitive to calculation error with the last terms as the final math sum is expected to be 0.0. Since the start/end a varies, the error varies.

A more stable result would be had adding the extreme f = sin(f(k)) values first. e.g. sum += sin(f(k=1)), then sum += sin(f(k=3)), then sum += sin(f(k=2)) rather than k=1,2,3. In particular the formation of term x=f(k=3) is likely a bit off from the negative of its x=f(k=1) earlier term, further compounding the issue.

Welcome to the world or numerical analysis.

Problem exists if code used all float or all long double, just different degrees.

Problem is not due to using an inexact value of pi (Exact value is impossible with FP as pi is irrational and all finite FP are rational).

Much is due to the formation of x. Could try the below to form the x symmetrically about 0.0. Compare exactly x generated this way to x the original way.

x = (x0-x1)/2 + ((k - n/2)*delta_x)

Print out the exact values computed for deeper understanding.

printf("x:%a y:%a\n", x0+(k*delta_x), func(x0+(k*delta_x)));
chux - Reinstate Monica
  • 143,097
  • 13
  • 135
  • 256
  • I see, thanks for your answer, the only thing that is not clear is why is more stable to sum in order k=1,3,2. I don't see the difference. – Spherk Mar 22 '21 at 21:47
  • In the case `n==4` and `x=f(k=2) --> 0`, no difference, but consider `n = 100`. Code is adding lots of values up only to negated then later . To maintain precision, the `sum` would need extra bits. How about `x`? Was the `x0+(1*delta_x)` the exact same as `x0+(3*delta_x)`? – chux - Reinstate Monica Mar 22 '21 at 21:57
  • From what I see, the value of `x0+(1*delta_x)` is the same as that of `x0+(3*delta_x)`, and that would be the case also for `n==100`. Then the issue would be as you mentioned, in the formation of `x=f(k)`, could that be related to the negative sign? which would occupy some memory "changing " the actual value for the negative component? – Spherk Mar 22 '21 at 23:26
  • "what I see, the value of x0+(1*delta_x) is the same as that of x0+(3*delta_x)" --> by analysis or did you make code that confirms that exact equivalence? – chux - Reinstate Monica Mar 23 '21 at 01:24
  • @trsk `*1`, `*2`, `/4` do not introduce round-off error. I suspect `3*delta_x` introduces an unexpected round-off error. – chux - Reinstate Monica Mar 23 '21 at 15:14
  • By analysis, if the sign does not change the number in some way, as I mentioned before, then `x0+(k*delta_x)= -pi/2.0+(k*(pi/n)) = pi*(-1/2.0 + k/n)` in the case `n==4.0`, and `k==1,2,3` we have that; `x0+(1*delta_x)= -pi/4.0`, `x0+(2*delta_x) = 0` (No matter the value of pi), and `x0+(3*delta_x) = pi/4.0`, for different test of pi the term `x0+(k*delta_x)` for `k==1` and `k==3` is equally spaced on left and right from the origin. But, clearly f(k=1) and f(k=3) doesn't add cero for certain values of pi. Anyway, I'll try to check if that is true with a program, thanks a lot. – Spherk Mar 23 '21 at 15:14
  • 1
    You are right. 3*delta_x introduces that round-off error. with `#include` and `cout< – Spherk Mar 26 '21 at 01:09