taylor series with error at most 10^-3

Question

I'm trying to calculate the the taylor series of cos(x) with error at most 10^-3 and for all x ∈ [-pi/4, pi/4], that means my error needs to be less than 0.001. I can modify the x +=in the for loop to have different result. I tried several numbers but it never turns to an error less than 0.001.

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

float cosine(float x, int j)
{
    float val = 1;
    for (int k = j - 1; k >= 0; --k)
        val = 1 - x*x/(2*k+2)/(2*k+1)*val;
    return val;
}

int main( void )
{
   for( double x = 0; x <= PI/4; x += 0.9999 )
   {

       if(cosine(x, 2) <= 0.001)
       {
           printf("cos(x) : %10g    %10g    %10g\n", x, cos(x), cosine(x, 2));
       }
       printf("cos(x) : %10g    %10g    %10g\n", x, cos(x), cosine(x, 2));
    }

    return 0;
}

I'm also doing this for e^x too. For this part, x must in [-2,2] .

float exponential(int n, float x)
{
    float sum = 1.0f; // initialize sum of series

    for (int i = n - 1; i > 0; --i )
        sum = 1 + x * sum / i;

    return sum;
}

int main( void )
{
    // change the number of x in for loop so you can have different range
    for( float x = -2.0f; x <= 2.0f; x += 1.587 )
    {
        // change the frist parameter to have different n value
        if(exponential(5, x) <= 0.001)
        {
            printf("e^x = %f\n", exponential(5, x));
        }
    printf("e^x = %f\n", exponential(5, x));
    }

    return 0;
}

But whenever I changed the number of terms in the for loop, it always have an error that is greater than 1. How am I suppose to change it to have errors less than 10^-3?

Thanks!

Answer 1

My understanding is that to increase precision, you would need to consider more terms in the Taylor series. For example, consider what happens when you attempt to calculate e(1) by a Taylor series.

$e(x) = \sum\limits_{n=0}^{\infty} frac{x^n}{n!}$

we can consider the first few terms in the expansion of e(1):

n             value of nth term           sum
0        x^0/0! = 1                   1
1        x^1/1! = 1                   2
2        x^2/2! = 0.5                 2.5
3        x^3/3! = 0.16667             2.66667
4        x^4/4! = 0.04167             2.70834

You should notice two things, first that as we add more terms we are getting closer to the exact value of e(1), also that the difference between consecutive sums are getting smaller.

So, an implementation of e(x) could be written as:

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

typedef float (*term)(int, int);
float evalSum(int, int, int, term);
float expTerm(int, int);
int fact(int);
int mypow(int, int);
bool sgn(float);

const int maxTerm = 10;         // number of terms to evaluate in series
const float epsilon = 0.001;    // the accepted error

int main(void)
{
    // change these values to modify the range and increment 
    float   start = -2;
    float   end = 2;
    float   inc = 1;

    for(int x = start; x <= end; x += inc)
    {
        float value = 0;
        float prev = 0;

        for(int ndx = 0; ndx < maxTerm; ndx++)
        {
            value = evalSum(0, ndx, x, expTerm);

            float diff = fabs(value-prev);
            if((sgn(value) && sgn(prev)) && (diff < epsilon))
                 break;
            else
                 prev = value;
        }

        printf("the approximate value of exp(%d) is %f\n", x, value);
    }

    return 0;
}

I've used as a guess that we will not need to use more then ten terms in the expansion to get to the desired precision, thus the inner for loop is where we loop over values of n in the range [0,10].

Also, we have several lines dedicated to checking if we reach the required precision. First I calculate the absolute value of the difference between the current evaluation and the previous evaluation, and take the absolute difference. Checking if the difference is less than our epsilon value (1E-3) is on of the criteria to exit the loop early. I also needed to check that the sign of of the current and the previous values were the same due to some fluctuation in calculating the value of e(-1), that is what the first clause in the conditional is doing.

float evalSum(int start, int end, int val, term fnct)
{
    float sum = 0;
    for(int n = start; n <= end; n++)
    {
        sum += fnct(n, val);
    }

   return sum;
}

This is a utility function that I wrote to evaluate the first n-terms of a series. start is the starting value (which is this code always 0), and end is the ending value. The final parameter is a pointer to a function that represents how to calculate a given term. In this code, fnct can be a pointer to any function that takes to integer parameters and returns a float.

float expTerm(int n, int x)
{
    return (float)mypow(x,n)/(float)fact(n);
}

Buried down in this one-line function is where most of the work happens. This function represents the closed form of a Taylor expansion for e(n). Looking carefully at the above, you should be able to see that we are calculating $\fract{x^n}{n!}$ for a given value of x and n. As a hint, for doing the cosine part you would need to create a function to evaluate the closed for a term in the Taylor expansion of cos. This is given by $(-1)^n\fact{x^{2n}}{(2n)!}$.

int fact(int n)
{
    if(0 == n)
        return 1;             // by defination
    else if(1 == n)
        return 1;
    else
        return n*fact(n-1);
}

This is just a standard implementation of the factorial function. Nothing special to see here.

int mypow(int base, int exp)
{
    int result = 1;

    while(exp)
    {
        if(exp&1)              // b&1 quick check for odd power
        {
            result *= base;
        }

        exp >>=1;              // exp >>= 1 quick division by 2
        base *= base;
    }

    return result;
}

A custom function for doing exponentiation. We certainly could have used the version from <math.h>, but because I knew we would only be doing integer powers we could write an optimized version. Hint: in doing cosine you probably will need to use the version from <math.h> to work with floating point bases.

bool sgn(float x)
{
    if(x < 0) return false;
    else return true;
}

An incredibly simple function to determine the sign of a floating point value, returning true is positive and false otherwise.

This code was compiled on my Ubuntu-14.04 using gcc version 4.8.4:

******@crossbow:~/personal/projects$ gcc -std=c99 -pedantic -Wall series.c -o series
******@crossbow:~/personal/projects$ ./series
the approximate value of exp(-2) is 0.135097
the approximate value of exp(-1) is 0.367857
the approximate value of exp(0) is 1.000000
the approximate value of exp(1) is 2.718254
the approximate value of exp(2) is 7.388713

The expected values, as given by using bc are:

******@crossbow:~$ bc -l
bc 1.06.95
Copyright 1991-1994, 1997, 1998, 2000, 2004, 2006 Free Software Foundation, Inc.
This is free software with ABSOLUTELY NO WARRANTY.
For details type `warranty'. 
e(-2)
.13533528323661269189
e(-1)
.36787944117144232159
e(0)
1.00000000000000000000
e(1)
2.71828182845904523536
e(2)
7.38905609893065022723

As you can see, the values are well within the tolerances that you requests. I leave it as an exercise to do the cosine part.

Hope this helps,
-T

Answer 2

exp and cos have power series that converge everywhere on the real line. For any bounded interval, e.g. [-pi/4, pi/4] or [-2, 2], the power series converge not just pointwise, but uniformly to exp and cos.

Pointwise convergence means that for any x in the region, and any epsilon > 0, you can pick a large enough N so that the approximation you get from the first N terms of the taylor series is within epsilon of the true value. However, with pointwise convergence, the N may be small for some x's and large for others, and since there are infinitely many x's there may be no finite N that accommodates them all. For some functions that really is what happens sometimes.

Uniform convergence means that for any epsilon > 0, you can pick a large enough N so that the approximation is within epsilon for EVERY x in the region. That's the kind of approximation that you are looking for, and you are guaranteed that that's the kind of convergence that you have.

In principle you could look at one of the proofs that exp, cos are uniformly convergent on any finite domain, sit down and say "what if we take epsilon = .001, and the regions to be ...", and compute some finite bound on N using a pen and paper. However most of these proofs will use at some steps some estimates that aren't sharp, so the value of N that you compute will be larger than necessary -- maybe a lot larger. It would be simpler to just implement it for N being a variable, then check the values using a for-loop like you did in your code, and see how large you have to make it so that the error is less than .001 everywhere.

So, I can't tell what the right value of N you need to pick is, but the math guarantees that if you keep trying larger values eventually you will find one that works.