Ask Your Question
1

Numerical approximation error involving cosh and sinh

asked 2021-05-20 15:54:31 +0200

Mathark gravatar image

updated 2021-05-21 20:08:08 +0200

slelievre gravatar image

Hi. I'm new to Sage and to computational mathematics in general. I have been trying to code something, but I am struggling to understand some numerical approximation errors. My problem seems to boil down to the following phenomenon:

sage: x = exp(12*pi*i*i)
sage: x.n()
0.000000000000000

Which is as expected, $x=e^{-12\pi }$ is "very close" to 0. However, I do not understand why Sage gives me

sage: (5+x).n()
4.00000000000000

If I give it more bits of precision, then it works:

sage: (5+x).n(200)
5.0000000000000000424115118301607754401746440550887456940681

But I can't understand why it is making such grotesque error. Sage knows $x$ is very small (it returned 0.000000000000000). So why did it approximate $5+x$ to 4?

What is also weird, is that if I repeat the same computation with 4+x or 6+x instead, then it works just fine:

sage: (6+x).n()
6.00000000000000

My apologies if this is something trivial in computational maths, but can someone clarify what is going on here?

In a situation like this, is there a way for me to know if I have given Sage enough bits of precision to get at least the first few digits correct?

(Remark: I wrote exp(12*pi*i*i) instead of exp(-12*pi) to make Sage use sinh and cosh.)

edit retag flag offensive close merge delete

Comments

Welcome to Ask Sage! Thank you for your question.

slelievre gravatar imageslelievre ( 2021-05-20 17:10:03 +0200 )edit

Congratulations! You found a bug!

slelievre gravatar imageslelievre ( 2021-05-20 17:10:13 +0200 )edit

Thanks for the report. This reminds me of

although it seems quite distinct.

slelievre gravatar imageslelievre ( 2021-05-20 17:13:29 +0200 )edit

So this is a bug? I was trying to do some computations with modular forms, but was getting some weird errors when evaluating them at some specific points of the upper halfplane. After a long time trying to understand what was going on, it boiled down to this issue here. So, should I open a report ticket? Or does this post already counts as one? Many thanks!!

Mathark gravatar imageMathark ( 2021-05-20 18:58:02 +0200 )edit

Oops. I got that wrong. Not a bug!

slelievre gravatar imageslelievre ( 2021-05-21 20:10:10 +0200 )edit

2 Answers

Sort by ยป oldest newest most voted
3

answered 2021-05-21 04:26:32 +0200

Juanjo gravatar image

updated 2021-05-21 04:31:52 +0200

A Sage real number $a$ can be written in the form $a=s\cdot m\cdot 2^e$, where $s$, $m$ and $e$ are its sign, mantissa and exponent, given by the method .sign_mantissa_exponent. For example,

sage: a = -12.345                                                               
sage: s, m, e = a.sign_mantissa_exponent(); s, m, e                             
(-1, 6949617174986097, -49)
sage: s*m*2^e.n()                                                               
-12.3450000000000

Let us see what is $x=\exp(-12\pi)=\exp(12\pi\cdot i\cdot i)$:

sage: x = exp(12*pi*i*i); x                                                     
cosh(12*pi) - sinh(12*pi)

We next compute the mantissa and exponent of $5+\cosh(12\pi)$ and $\sinh(12\pi)$:

sage: (5 + cosh(12*pi)).n().sign_mantissa_exponent()                            
(1, 5894625992139550, 1)
sage: sinh(12*pi).n().sign_mantissa_exponent()                                  
(1, 5894625992139548, 1)

So SageMath gets $5+x$ as follows:

$$\begin{aligned}5+x&=5+\cosh(12\pi)-\sinh(12\pi) \\ &\approx 5894625992139550\cdot 2^1 - 5894625992139548\cdot 2^1=2\cdot 2=4.\end{aligned}$$

This justifies the result you get. You can repeat computations to check what happens with $6+x$:

sage: (6 + cosh(12*pi)).n().sign_mantissa_exponent()
(1, 5894625992139551, 1)

Hence: $$\begin{aligned}6+x&=6+\cosh(12\pi)-\sinh(12\pi) \\ &\approx 5894625992139551\cdot 2^1 - 5894625992139548\cdot 2^1=3\cdot 2=6.\end{aligned}$$

The problem here is that you are reaching the limits of the double precision floating point arithmetic. Just use more digits:

sage: x = exp(12*pi*i*i)                                                        
sage: R = RealField(100)                                                        
sage: R(5 + x)                                                                  
5.0000000000000000000000000000
edit flag offensive delete link more

Comments

Thanks you so much!! This answer was incredibly helpful. So, it was not a bug!! Though, is there a way for me to know if I am using enough bits of precision, so that the first few digits of a given calculation are correct? I guess I could be using .sign_mantissa_exponent to check if the exponents in my calculations are small, right?

For example, having exponent 1 means the computer can only see things with precision $2^1$. So 4 and 5 are equal in this case. However, if I give Sage enough bits to make all the exponents in my calculations lower than $-10$, then I would have precision of $2^{-10}$. Is my reasoning correct?

Mathark gravatar imageMathark ( 2021-05-21 15:48:18 +0200 )edit

Now, with this in mind, I think I could code some kind of precision handling into many parts of my code. But I wonder, doesn't Python/Sage have some sort of command to deal with this kind of issue automatically?

Say, something that I could initiate at the start of my program. And then, if one of my numerical calculations gets exponents too big, it gives me a warning massage. Or something similar, to make sure I have used enough precision. Once again, many thanks!!

Mathark gravatar imageMathark ( 2021-05-21 15:59:49 +0200 )edit
2

answered 2021-05-21 20:34:29 +0200

mwageringel gravatar image

Have a look at real or complex ball fields, which keep track of approximation error.

sage: x = exp(12*pi*i*i)
sage: RBF(x)
[+/- 67.1]
sage: 5 + RBF(x)  # not accurate enough
[+/- 72.1]
sage: 5 + RealBallField(80)(x)  # more precision
[5.00000 +/- 1.75e-6]
edit flag offensive delete link more

Comments

That is what I needed. Thanks!

Mathark gravatar imageMathark ( 2021-05-22 15:06:15 +0200 )edit

Your Answer

Please start posting anonymously - your entry will be published after you log in or create a new account.

Add Answer

Question Tools

Stats

Asked: 2021-05-20 15:54:31 +0200

Seen: 73 times

Last updated: May 21