1

I have a question about adding the number 1 to very small numbers. Right now, I am trying to plot a circular arc in the complex plane centered around the real number 1. My code looks like:

arc = 1 + rho .* exp(1i.*theta);

The value rho is a very small number, and theta runs from 0 to pi, so whenever 1 is added to the real part of arc, MATLAB seems to just round it to 1, so when I type in plot(real(arc),imag(arc)), all I see is a spike instead of a semicircle around 1. Does anyone know how to remedy this so that MATLAB will not round 1 + real(arc) to 1, and instead conserve the precision?

Thanks

aka.nice
  • 9,100
  • 1
  • 28
  • 40
user1824904
  • 19
  • 1
  • Can you provide the exact code that demonstrates this problem? Not really any way to give a definitive solution unless we can reproduce the problem. – Dan Becker Nov 14 '12 at 21:53
  • Bitwise, sorry if you were offended by my downvote! But I don't think that I was the only person who didn't understand the relevance of your answer. I downvoted (and added a comment explaining why), you updated with a fuller explanation, and now everyone reading this understands your point better. Seems like an example of the community working to me. – Dan Becker Nov 15 '12 at 00:10
  • @DanBecker No problem and no hard feelings - my fault for not clarifying how the answer is relevant. – Bitwise Nov 15 '12 at 00:45

3 Answers3

3

rho=1e-6; theta=0:pi/100:pi; arc=1+rho*exp(1i.*theta); plot(arc); figure(); plot(arc-1);

Shows, that the problem is in plot, not in loss of precision. After rho<1e-13 there will be expected trouble with precision.

The two other possible misconceptions:
- doubles have finite precision. 16 decimal digits or 1+2^-52 is the limit with doubles.
- format short vs. format long -- matlab shows by default only 6 or 7 digits

It also happens to be that 6-7 digits is the limit of a 32-bit float, which could explain also that perhaps the plot function in Octave 3.4.3 is also implemented with floats.

Left: 1+1e-6*exp, Right: (1+1e-6*exp)-1

Aki Suihkonen
  • 19,144
  • 1
  • 36
  • 57
  • The left is for 1e-6, I assume? In freemat I need get a similar picture with 1e-15, but not earlier. – arne.b Nov 14 '12 at 21:10
  • Aki, is the code in your post the code that you used to generate the two plots you show? When I run that code, both plots look fine to me (R2012a on Windows 8). – Dan Becker Nov 14 '12 at 21:53
  • Ah, I see that you are using Octave. So it appears that you have demonstrated a problem in Octave, but matlab appears to not have this problem. – Dan Becker Nov 14 '12 at 21:55
  • @DanBecker Not quite. With MATLAB, I get no picture like the one on the left, but the half-circle starts becoming to narrow for `rho` between `1e-11` and `1e-12` (where it already looks quite "spiky"), so MATLAB's plot are somewhat unusual as well. (`plot(arc-1)` is still normal, as are both plots for these rho values in freemat, so I think Aki is on to something here.) – arne.b Nov 14 '12 at 22:11
2

There is a builtin solution for exactly this probem:

exp1m()

log1p()

explicitly:

log(arc)=log1p(rho*exp(1i*theta))

to get what you need.

Of course you need to work in log space to represent this precision, but this is the typical way this is done.

Bitwise
  • 7,577
  • 6
  • 33
  • 50
  • 1
    Useful functions, but not relevant for this question (the argument of `exp` is not a small number here). – Dan Becker Nov 14 '12 at 21:49
  • @DanBecker define eps=rho*exp(1i*theta) and then log(arc)=log1p(eps). I added it to the answer so it would be written explicitly. – Bitwise Nov 14 '12 at 23:14
1

In double precision floating point representations, the smallest number strictly greater than 1 that can be represented is 1 + 2^-52.

This is a limitation imposed by the way non-integer numbers are represented on most machines that can be avoided in software, but not easily. See this question about approaches for MATLAB.

Community
  • 1
  • 1
arne.b
  • 4,212
  • 2
  • 25
  • 44
  • @AkiSuihkonen Not in my octave... did you mean 1e-16 (which would be <`2^-53`)? – arne.b Nov 14 '12 at 20:58
  • I mean rho=1e-6 as in my answer. The plot (arc) is staggered, while plot(arc-1) is smooth. Only after rho<1e-14 1+exp(i*theta) can't hold the precision as is expected. My Octave is 3.4.3. – Aki Suihkonen Nov 14 '12 at 21:02
  • @AkiSuihkonen Ok, I just checked the real value of arc(1) and arc(2) manually, never used octave's plot. With freemat I get great half-circles up to `rho=1e-13`. – arne.b Nov 14 '12 at 21:07
  • While Aki's post is interesting, I strongly suspect that this is the real problem, and there is no real solution other than to decide to plot the arc around a more convenient point (like 0). – Dan Becker Nov 14 '12 at 21:59