3
$\begingroup$

Here's the sun's declination for a 10 year period (roughly 2000-2010):

sd = Table[AstronomicalData["Sun", {"Declination", DateList[t]}],
 {t, 3155716800, 3155716800+86400*365.2425*10, 86400}];

As expected, it follows a sinusoidal pattern:

ListPlot[sd, PlotRange->All]

enter image description here

When I take fast Fourier transform, however, the max amplitude is:

Max[Abs[Fourier[sd]]] // InputForm

InputForm= 0``-3.9244986797945733

(this also prints out as 0. * 10^4 (with the 4 superscripted))

The min amplitude is 0``-2.324244273207087 or 0. * 10^3, so I strongly suspect a loss of precision somewhere, but I'm not seeing where or how.

I've read through similar questions on this site, but I believe none of those cases apply to my situation.

Any thoughts?

Improve this question
$\endgroup$
5
  • 3
    $\begingroup$ The input is comprised of low precision "bignums", with precision around 6 digits. It appears to be a precision loss, and my guess is Fourier is not using fixed precision arithmetic (I do not know if that is good or bad). You might try converting to machine doubles: In[61]:= Max[Abs[Fourier[N[sd[[All, 1]]]]]] // InputForm Out[61]//InputForm= 703.066438973999 $\endgroup$
    – Daniel Lichtblau
    Jan 13, 2016 at 19:31
  • $\begingroup$ @DanielLichtblau Thanks. I ended up using Rationalize[sd,0], but I'm still not happy about why this happened. This should be well within $MachinePrecision, shouldn't it? $\endgroup$
    – user1722
    Jan 13, 2016 at 20:49
  • 1
    $\begingroup$ Apparently not. Some experimentation indicates that the FT loses around 3.5 digits in precision, likely from standard significance arithmetic during the computation. The input has very low precision (some values have only 3 or so digits), hence the catastrophic loss. $\endgroup$
    – Daniel Lichtblau
    Jan 13, 2016 at 21:54
  • $\begingroup$ @DanielLichtblau Write that up as an answer and I'll approve it. $\endgroup$
    – user1722
    Jan 15, 2016 at 16:56
  • $\begingroup$ I posted a response, eventually realizing as I wrote it that I had once addressed this in past (see end note). So this might actually end up getting closed as a duplicate. $\endgroup$
    – Daniel Lichtblau
    Jan 15, 2016 at 17:29

1 Answer 1

Reset to default
5
$\begingroup$

We'll do this at a few different precisions. First we show the result from the original input and also from machine doubles created from that input.

sd = Table[
   AstronomicalData["Sun", {"Declination", DateList[t]}], {t, 
    3155716800, 3155716800 + 86400*365.2425*10, 86400}];
sdvals = sd[[All, 1]];
Max[Abs[Fourier[sdvals]]]
Max[Abs[Fourier[N[sdvals]]]]

(* Out[266]= 0.*10^3

Out[267]= 703.066438974 *)

The result based directly on the input is showing full cancellation, that is, no significant digits (notice it has the right accuracy, that is, three orders of magnitude). We can understand what is happening by gauging "precision loss" during the course of operations. The initial precision is indicated below.

Precision[sdvals]

(* Out[271]= 3.31426802501 *)

We bump to 12 digits precision for input, and see a result with but two digits correct, so the loss is 10 digits.

ft12 = Fourier[SetPrecision[sdvals, 12]];
Max[Abs[ft12]]
Precision[ft12]

(* Out[277]= 703.06644

Out[278]= 2.15175705976 *)

(I had claimed 3.5 digits loss in a comment. That was only correct for the max value, others lost more.)

To understand the precision loss, well, it turns out this is really a duplicate. I provide some explanation in a response here.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.