Subtraction Doubles in C++

Question

I'm doing a homework assignment where I need to use the McLaurin Series to evaluate for cosh(x) ans sin(x). Each time the program runs, the user inputs a number from 0-10 to specify how many terms to evaluate the series too. The user then inputs the value of x to evaluate to. Once the user hits enter, 10 increments of the series will be calculated and printed out. In addition, I have to find the exact error as well. To do this, I created several different functions for both sinh and cosh that performed the calculations of a specific term. For example: the code below evaluates the term up to 10!.

void cosFunction10(double power, double value)
{
     double cosh_x = 0;
     double math_cosh = 0;
     double exact_error;
     double x;

     if (value < 0)
         x = -0.1*0;

     for (int i = 0; i < 11; ++i)
     {
         cosh_x = (1.0 + ((x*x)/(2.0*1.0)) + ((x*x*x*x)/(4.0*3.0*2.0*1.0)) + ((x*x*x*x*x*x)/(6.0*5.0*4.0*3.0*2.0*1.0)) + ((x*x*x*x*x*x*x*x)/(8.0*7.0*6.0*5.0*4.0*3.0*2.0*1.0)) + ((x*x*x*x*x*x*x*x*x*x)/(10.0*9.0*8.0*7.0*6.0*5.0*4.0*3.0*2.0*1.0)));
         math_cosh = cosh(x);
         exact_error = (math_cosh - cosh_x);

         cout << setiosflags(ios::scientific) << setprecision(3) << x << "\t";
         cout << setiosflags(ios::scientific) << setprecision(5) << cosh_x << "\t";
         cout << math_cosh << "\t";
         cout << exact_error << endl;
         x = x + (0.1*value);
     }
 }

If I run the program and put in these values: 10 (the nth term) and -6 (value for x).

This is what I should get: (2nd increment)

 x                Series            Exact (using the cmath functions)    Exact % Error
 -6.000e-001      1.18547e+000      1.18547e+000                        0.00000e+000

The next 2 out puts would be the same for the exact error, until I get to the 4th increment

 x                Series            Exact (using the cmath functions)    Exact % Error
 -1.800e+100      3.10747e+000      3.10747e+000                        -7.67243e-005

When I run my code, I do not get the above results: (2nd increment)

 x                Series            Exact (using the cmath functions)    Exact % Error
-6.000e-001      1.18547e+000      1.18547e+000                         4.55325e-012

So it seems that something is going wrong because my exact error isn't even close to the same as to what I'm given. I understand that there may be slightly different values, which is expected but not to this extent.

Any suggestions would be much appreciated,

Thank you

Answer 1

There MUST be such a difference because cosh(x) is defined in a different way which is according to the formula given below:

Please see here for more details.

However, if you are calculating the amount based on computation of MacLaurin series, the problem is that you are not accumulating the values in the right way. The below implementation is perhaps what you actually needed:

void cosFunction(double x)
{
    double cosh_x = 1;
    double math_cosh;
    double exact_error;
    double factorial = 1;
    double x_pow = 1;

    for (int i = 0; i < 11; ++i) {
            cosh_x = x_pow/factorial;
            x_pow = x*x*x_pow; 
            factorial = factorial* (2 * i + 1) * (2 * i + 2);
    }

    math_cosh = cosh(x);
    exact_error = (math_cosh - cosh_x);
    cout << setiosflags(ios::scientific) << setprecision(5) << exact_error << endl;

}

Answer 2

As revealed in the comments, the error you get is actually extremely small, less than 10^-11. In addition, you must take into account that the way you compute cosh and the way it's done by the built-in function is not the exact same. Both methods are ultimately approximations, but they work in different ways, hence the rounding errors and approximation errors they produce will be slightly different.

You cannot expect the error to be "flat 0".

Answer 3

The correct answer may be, in a mathematical sense, zero. But in a practical computing sense, it won't be zero:

Imagine if you're doing arithmetic with fixed decimal precision, say 8 digits. You would represent 1.0/3.0 as .33333333. But that would make 3.0 * (1.0/3.0) turn into 3.0 * .33333333 which would be .999999999. So while in a mathematical sense, 1.0 - (3.0 * (1.0 / 3.0)) should be zero, with fixed decimal precision, you may well get something like .00000001.

The same rules apply to fixed binary precision, which is what you get when you use a double.

Answer 4

//Part of the problem is that in the line after the 'for' nothing changes (typo?.) The
//following works for x = 0.5`
//Note that the loop works properly for all terms including the first which is 1.0

#include <iostream>
#include <math.h>

using namespace std;

int factorial(int n)
{
    return (n == 1 || n == 0) ? 1 : factorial(n - 1) * n;
}

int main()
{
    // use taylor series to approximate cosh(x) - use 7 terms
    double x = 0.5;
    double top, bot;
    double s = 0.0;

    for (int i = 0; i < 7; i++)
    {
        top = pow(x,double(i)*2.0);
        bot = double(factorial(i*2));
        s  += (top / bot);
    }
    cout << "myCosh     " << s << endl;
    cout << "cosh       " << cosh(x)<< endl;
    cout << "difference " << fabs(s - cosh(x)) << endl;
    return 0;
}

Answer 5

I ran your function with the cout operation modified as follows:

cout << setiosflags(ios::fixed) << setprecision(12);
cout << "for x="<<x<<": "<<math_cosh<<" - "<<cosh_x<<" = "<<exact_error<<"\n";

I chose my own value for value, since you didn't indicate what you were using.

I get the following result, which demonstrates that math_cosh and cosh_x are not the same. Your problem seems to be that you didn't see enough digits to notice that they were different.

for x=0.003233200000: 1.000005226796 - 1.000005226791 = 0.000000000005
for x=0.006466400000: 1.000020907237 - 1.000020907164 = 0.000000000073
for x=0.009699600000: 1.000047041489 - 1.000047041120 = 0.000000000369
for x=0.012932800000: 1.000083629824 - 1.000083628658 = 0.000000001166
for x=0.016166000000: 1.000130672624 - 1.000130669778 = 0.000000002846
for x=0.019399200000: 1.000188170381 - 1.000188164480 = 0.000000005901
for x=0.022632400000: 1.000256123697 - 1.000256112765 = 0.000000010932
for x=0.025865600000: 1.000334533282 - 1.000334514632 = 0.000000018650
for x=0.029098800000: 1.000423399955 - 1.000423370081 = 0.000000029875
for x=0.032332000000: 1.000522724646 - 1.000522679112 = 0.000000045534
for x=0.035565200000: 1.000632508392 - 1.000632441726 = 0.000000066667

Regarding the homework assignment, to calculate cosh(x) using the Maclaurin series: Your function, as written, is not very precise since it only computes the first three terms of the series 1 + 0 + x*x/2. You are not computing x^n or n! beyond n=2. To generate a more precise result, you need to be accumulating more terms of the series -- as in the following...

int NUMBER_OF_TERMS=10; // the number of terms to compute (including zero terms)
for ( int i=0; i<NUMBER_OF_TERMS; i++ )
   {
   cosh_x += (i%2==0?1:0)*term; // even derivatives are 1, odd derivatives are 0
   term *= x/(i+1); // apply incremental power and factorial
   }

This still doesn't guarantee exactness. Due to rounding errors, you can't be sure to perfectly match the built-in function cosh(x) unless you duplicate its algorithm exactly. Both functions will only produce an approximation, since for most input values, the pure mathematical output would be a transcendental number -- and such a number cannot be exactly represented in a double.

You should also keep in mind that the Maclaurin series is a polynomial-based approximation which is most accurate near x=0. As is the nature of polynomials, an incomplete summation may eventually (as x moves away from 0) diverge very rapidly from the correct answer. This is why, in the example output above, the error values are smaller for smaller values of x.

---

Regarding scientific notation vs. fixed-point precision:

Scientific notation shows you the most significant non-zero digits of a number -- even if those digits start 12 places to the right of the decimal point. When reading scientific notation, it is important to be aware of the information after the e since this describes the true magnitude of the number. double-precision values can store approximately 20 decimal digits of real information. This is nice for doing many repeated calculations, but it usually isn't very meaningful when comparing two numbers for equality.

What you can do, although this may not be in the scope of your assignment, is limit the fixed-point precision of the difference. Here is an example of how you might do that:

const double PRECISION=1e-4; // indicates the number of fractional digits to keep
exact_error = round(exact_error/PRECISION)*PRECISION; // limit precision

...or it might be simpler to just check for the close-to-zero case ...

if ( fabs(exact_error)<PRECISION/2 ) exact_error = 0; // close enough to zero