File Exchange

image thumbnail

Trigonometric function errors at pi or pi/2 such as sin(pi) or cos(pi/2)

version 1.7.0.0 (3.85 KB) by Sung-Eun Jo
The constant pi in matlab is a floating-point approximation of 'pi'. There is a remedy for it!

1 Download

Updated 29 Sep 2014

View License

The pi in matlab is not a real 'pi',
but it is only a floating-point number close to 'pi'.
Trigonometric functions around pi may have errors close to machine precision.
See:
cos(pi/2) = 6.12323399573677e-17
sin(pi) = 1.22464679914735e-16
exp(i*pi/2) = 6.12323399573677e-17 + 1i
exp(i*pi) = -1 + 1.22464679914735e-16i
Note that sin(pi) returns '1.22464679914735e-16' and it is true somehow.
These unwanted results do not come from function evaluations. It is because pi value cannot be represented in floating-point number.
One approach to go around it is to use symbolic toolbox: sin(sym(pi)) returns 0. Note that sym(pi) is the exact pi.
See http://www.mathworks.co.kr/kr/help/symbolic/sin.html
There is another way without using symbolic toolbox.
You had better use degrees rather than radians
when your application needs to use correct values around such angles.
Note that sind(180) returns '0'.
The functions sind and cosd use degrees:
cosd(90)
sind(180)
Or, we can define functions with radian input:
cos_rad = @(theta) cosd(theta/pi*180)
sin_rad = @(theta) sind(theta/pi*180)
cos_rad(pi/2) = 0
sin_rad(pi) = 0
The exp function is not easy to handle with degree.
You can try the following definition:
exp_rad = @(theta) exp(real(theta)).*(cos_rad(imag(theta)) + i* sin_rad(imag(theta)))
exp_rad(i*pi/2)
exp_rad(i*pi)
The function exp_rad slightly abuses the definition pi.
This would be a trick for handling 'pi' as if it is a real 'pi'.
pi is not a rational number, but the 180 degrees are rational number.
Converting radians to degrees makes the exact result around pi.
The exp_rad also responses with matrix input:
A = [log(2) 0; i*pi/2 i*pi];
exp(A)
expm(A)
exp_rad(A)
[V,D]=eig(A);
V*diag(exp_rad(diag(D)))/V % expm_rad(A) example

Now we can provide 4 functions with abused pi definition:
cos_rad.m
sin_rad.m
exp_rad.m
expm_rad.m

Those show correct results around pi/2 and pi points if we want to take 'pi' as the real 'pi'.
These files are very short sample codes for understanding.
You can use them when your calculation needs to be exact zero with sin(N*pi) or cos((N+1/2)*pi).

For large angle > 6*pi radians, the *_rad functions happen to show non-zero at the intended points due to further loss of precision of large multiple of pi. See the comments on sin(pi) at http://www.mathworks.co.kr/company/newsletters/articles/the-tetragamma-function-and-numerical-craftsmanship.html

We can count the number of points where abs(sin(k*pi))<eps or abs(cos(k*pi+pi*0.5))<eps in order to see how accurate the trigonometric functions calculate for large k's.
N=5000;
tt=linspace(-N,N,8*N+1)*pi;
% count the points where the output is less than epsilon.
sum(abs(sin(tt))<eps(tt))
ans = 10001
sum(abs(cos(tt))<eps(tt))
ans = 10000
These points can be assumed to be 0 without loss of sin or cos natures.
Note that sum(abs(sin(tt))==0) and sum(abs(cos(tt))==0) return 1 and 0 because of floating-point errors.
Find more examples at test_pi_functions.m

Again we define *_eps functions that return zeros when their outputs are less than eps:
More functions are added:
cos_eps.m
sin_eps.m
exp_eps.m
expm_eps.m

Due to output comparison, _eps functions are slow a bit. If we know the phase angles are less 6*pi, we can still use _rad functions.

--jo

Cite As

Sung-Eun Jo (2020). Trigonometric function errors at pi or pi/2 such as sin(pi) or cos(pi/2) (https://www.mathworks.com/matlabcentral/fileexchange/47894-trigonometric-function-errors-at-pi-or-pi-2-such-as-sin-pi-or-cos-pi-2), MATLAB Central File Exchange. Retrieved .

Comments and Ratings (9)

Thank you for the submission. As of 2018b, 'cospi' and 'sinpi' is introduced to handle this mismatch. https://www.mathworks.com/help/matlab/ref/cospi.html

You can use sym in cos.
cos(sym(pi/2))

Varoujan

Thanks but the above explanation doesn't explain the following:
1000000000000000.0 * cosd(0.01) = 1.0000e+15
while
num2str(cosd(0.010),15) = 0.999999984769129

The later is correct. So, what would be the correct form to use sin/cos or sind/cosd in math near the 0, 90, 180, 270, 360 degree values?

Matthew

Thanks John for adding some sanity.

John D'Errico

If, in your case, A is a diagonal matrix, why do you think you need to compute an eigen-decomposition to compute the matrix exponential? That alone will cause unwanted numerical trash to accumulate.

Sung-Eun Jo

@John Thank you for reading. exp_rad and expm_rad are silly for naming and definition... But they have a very limited purpose such that...
Given A = diag([pi/2 pi]),
we want to calculate 'expm(i*A)' as a complex-phase rotating matrix.
expm(i*A) is theorentically diag([i -1]), but
it returns diag([6.1232e-17+1i -1+1.2246e-16i]).
These numerical errors are propagating unwantedly.
In my previous application, A = diag([anyangle1 any angle2 ... ]), angles are determined by problem and happen to be 'pi' easily.
As John's comment, such trivial cases can be done with matlab functions sind and cosd. However, we can introduce some silly mapping functions exp_rad or expm_rad, without changing radians into degrees.

John D'Errico

This is just a confusing set of tools, with very little purpose behind them.

Matrix exponentials, with degree input versus radians? Are you kidding? Perhaps you misunderstand what a matrix exponential is, or what it is used for. Anyway, there is already a function in MATLAB called EXPM, which does a very credible job for the purpose. And since it was Cleve who wrote the paper called "19 Dubious Ways to Compute the Matrix Exponential", maybe they know how to solve the problem in general.

As far as sin_rad goes, it is simply better to learn about floating point arithmetic and how to deal with it than to worry about whether sin(pi) returns a non-zero.

While the functions do have help in them, the tools themselves are just a silly idea.

Updates

1.7.0.0

update comments

1.6.0.0

typo

1.5.0.0

Thanks to others feedback, I added more functions available with large phase angles.

1.4.0.0

typo

1.3.0.0

update comments.

1.2.0.0

comments

1.1.0.0

typo

MATLAB Release Compatibility
Created with R2012b
Compatible with any release
Platform Compatibility
Windows macOS Linux