Why am I getting the wrong sign on my cos(x) approximation?

Question

I am writing a Matlab script that will approximate sin(x) and cos(x) using their Maclaurin polynomials.

When I input arg = (5*pi)/4 I expect to get the correct approximations for

sin((5*pi)/4) = -0.7071067811865474617

cos((5*pi)/4) = -0.7071067811865476838.

Instead I get the following when running the script:

Approximation of sin(3.92699)    >> -0.7071067811865474617
Actual sin(3.92699)               = -0.7071067811865474617
Error approximately               = 0.0000000000000000000 (0)
----------------------------------------------------------
Approximation of cos(3.92699)    >> 0.7071067811865474617
Actual cos(3.92699)               = -0.7071067811865476838
Error approximately               = 0.0000000000000001110 (1.1102e-16)

I am getting the correct answers for sin but incorrect for cosine when the argument (angle) is in quadrant 3 or 4. The problem is that I am getting the wrong sign on the cos(arg) value. Where have I messed up?

---

CalculatorForSineCosine.m

% Argument for sine/cosine in radians.
arg = (5*pi)/4;

% Move the argument x so it's within [0, pi/2].
newArg = moveArgumentV2(arg);

% Calculate what degree we need for our Taylorpolynomial.
TOL = 0; % If 0, assume we want Machine Epsilon.
r = findDegreeV2(TOL);

% Plot nth degree Taylorpolynomial around x = 0 for sine.
% and calculate approximation of sin(x).
[approximatedSin, errorSin] = sin_taylorV2(r, newArg);
eS = num2str(errorSin); % errorSin in string format

% Plot nth degree Taylorpolynomial around x = 0 for cosine.
% and calculate approximation of cos(x).
[approximatedCos, errorCos] = cos_taylorV2(r, newArg);
eC = num2str(errorCos); % errorCos in string format

% Print out the result.
fprintf('\nApproximation of sin(%.5f)\t >> %.19f\n', arg, approximatedSin);
fprintf('Actual sin(%.5f)\t\t\t\t  = %.19f\n', arg, sin(arg));
fprintf('Error approximately\t\t\t\t  = %.19f (%s)\n', errorSin, eS);
disp("----------------------------------------------------------")
fprintf('Approximation of cos(%.5f)\t >> %.19f\n', arg, approximatedCos);
fprintf('Actual cos(%.5f)\t\t\t\t  = %.19f\n', arg, cos(arg));
fprintf('Error approximately\t\t\t\t  = %.19f (%s)\n\n', errorCos, eC);

sin_taylorV2.m

function [approximatedSin, errorSin] = sin_taylorV2(r, x)
%% sss

    % Q_2n+1(x) where 2n+1 = degree of polynomial.
    n = (r - 1)/2;

    % Approximate sin(x) using its Taylorpolynomial.
    approximatedSin = 0;
    for k = 0:n
        approximatedSin = approximatedSin + (((-1).^k) .* (x.^(2.*k+1)))./(factorial(2.*k+1));
    end

    % Calculate the error.
    errorSin = abs(sin(x) - approximatedSin);

end

cos_taylorV2.m

function [approximatedCos, errorCos] = cos_taylorV2(r, x)
%% sss

    % Q_2n+1(x) where 2n+1 = degree of polynomial and n = # terms.
    n = (r - 1)/2;

    % Approximate cos(x) using its Taylorpolynomial.
    approximatedCos = 0;
    for k = 0:n
        approximatedCos = approximatedCos + (((-1).^k) .* (x.^(2.*k)))./(factorial(2.*k));
    end

    % Calculate the error.
    errorCos = abs(cos(x) - approximatedCos);

end

moveArgumentV2.m

function newArg = moveArgumentV2(arg)
%% Moves the argument x to the interval [0, pi/2].

    % Make use of sines periodocity and choose n as ceil( (x-pi)/2pi) )
    n = ceil((arg-pi)/(2*pi));
    x1 = arg - 2*pi*n; % New angle will be in [-pi, pi]

    x2 = abs(x1); % Angle will be in [0, pi]
    if (x2 < pi/2) && (x2 > 0)
        x3 = x2;
    else
        x3 = pi - x2;
    end
    newArg = x3*sign(x1); % Angle will be in [0, pi/2]

end

Answer 1

I would like to notice two things in your code.

First, you don't need the moveArgumentV2(arg) function, as, if you remember, the radius of convergence for the Maclaurin/Taylor series of the sin(x)/cos(x) is the set of all real numbers. That means the series should converge for any real x, disregarding the round-off errors inherently to every arithmetic operations done in a computer.

As a matter of fact, following your code, we can write a function that approximates the cos as:

function y = mycos(x,n)
    y = 0;

    for k=0:n
        term = (-1)^k*x.^(2*k)/factorial(2*k);
        y = y + term;
    end

end

Notice this function works for values outside the range [-pi,pi]:

x = -10*pi:0.1:10*pi;
ye = cos(x) % exact value
ya = mycos(x,100) % approximated value
plot(x,ye,x,ya,'o')

The values returned by the mycos function are close to the exact value given by the cos built-in function. This happens because I calculated the approximation with the first 100 terms. The error, however, for higher values of x, is extremely large if we use just a few terms.

ya = mycos(x,10) % approximated value with 10 terms only
plot(x,ye-ya); title('error')

The problem now is that we can't just increase the number of terms without running in another problem. If we increase the number of points, the mycos function crumbles due to round-off errors, because of the factorial function that overflows. A good idea is to try to change your code in order to avoid the use of the factorial function. Notice the recurrence between sucessive terms in the Maclaurin expansion of the cos function, and you can create another function without the use of the factorial:

function y = mycos2(x,n)
    term = 1;
    y = 1;

    for k=1:n
        term = -term.*x.^2/(2*k-1)/(2*k);
        y = y + term;
    end
end

Here, we calculate each term in the series expansion from the previous calculated term. We avoid the calculation of the factorial and make use of what we already have. This speeds the code and avoids overflow. As a matter of fact, if we now calculate the cos approximation with 500 terms, we get:

x = -10*pi:0.5:10*pi;
ye = cos(x); % exact value
ya = mycos(x,500); % approximated value
ya2 = mycos2(x,500); % approximated value
plot(x,ye,x,ya,'x',x,ya2,'s')
legend('ye','ya','ya2')

Notice in this figure the x marks are the calculations done with the mycos function, while the o marks are done without using the factorial function. The first function crumbles for values outside the range [-2,2], but the second one runs just fine. It works even when I use 1e5 terms. Increasing the number of terms reduces the errors, so you can estimate how much terms you will use on an approximation, given a desired tolerance. If this number is greater than 170, the first function will not work properly.

factorial(170) returns 7.2574e+306, but factorial(171) returns Inf, so any value that should be calculated with more than 170 terms will have problems in the first function. Avoid the calculation of factorial at all costs.

Answer 2

This is what I tried:

x = -3*pi:0.01:3*pi;
y = x;
for ii=1:numel(y)
   y(ii) = moveArgumentV2(y(ii)); % not vectorized
end
plot(sin(x))
hold on
plot(sin(y))

Both sin(x) and sin(y) produce the same plot. But:

plot(cos(x))
hold on
plot(cos(y))

Now we see that cos(x) and cos(y) are not the same! This is because moveArgumentV2 changes the angle to be in the first and fourth quadrant (in the range [-pi/2, pi/2]), which is what you need for the sin function, but is not adequate for the cos function.

I would modify sin_taylorV2 and cos_taylorV2 to call moveArgumentV2, so you don't rely on the caller to know what the valid input range is. In cos_taylorV2 you would need to call it this way:

   x = moveArgumentV2(x+pi/2) - pi/2;

and in sin_taylorV2 you'd call it the same way you do now.

Or, better, write cos_taylorV2 in terms of sin_taylorV2, which we know to be correct. This avoids code duplication.