1

As a hobby project I'm trying to implement something with a Fourier series for an arbitrary function. I have a library that can integrate real functions numerically. I'm using the complex notation for a Fourier series, so I had to rewrite the expression for $c_n$ to: $$ c_n = \int_0^1 r(t) \cos{\lambda_n(t)} dt + i \int_0^1 r(t) \sin{\lambda_n(t)} dt $$ where $\lambda_n(t) \equiv \phi(t) - 2 \pi n t$ and some original function $f(t) = a(t) + i b(t) = r(t) \exp{i \phi(t)}$.

Then, I calculate some number $N$ coefficients and determine: $$ g(t) = \sum_{n=-N}^{+N} c_n \exp{2 \pi i n t} $$

Plotting this function $g(t)$ in white and the original function $f(t)$ in red, this results in: plot

As can be seen, the approximation has some weird curvature. I suspect this is caused by the boundary conditions or the fact that it is periodic and "jumps" from start to end. Any idea what is wrong here? I would like to achieve a similar result as in the last 3B1B video.

Edit: Additional discovery: for $N >> 1$ this effect blows up. See this for $N=75$: blowup 1

For $N=85$: blowup 2

kwantuM
  • 151
  • 2
    You should take a look at Gibbs phenomenon, that will tell you (unfortunately) that there's a limit to what you can achieve – postmortes Jul 02 '19 at 11:14
  • @postmortes https://en.wikipedia.org/wiki/Gibbs_phenomenon#/media/File:Gibbs_phenomenon_10.svg here the oscillations do not curve – kwantuM Jul 02 '19 at 11:17
  • 1
    @postmortes I think the OP is talking about how the white dots curve when moving through the discontinuities. – Matthew Cassell Jul 02 '19 at 11:18
  • 3
    What you see is a combination of the Gibbs phenomenon (causing the reconstructed signal to be oscillatory) and of aliasing (causing the sampling of the oscillatory signal to generate spurious patterns by accidental alignment). On top of that, the numerical errors explain the asymmetry. –  Jul 02 '19 at 13:00
  • @YvesDaoust I'm using 64 bit floats and integrating to less than 1e-9 error per coefficient. I'm having trouble believing that 64 bit precision causes gigantic issues like this. – kwantuM Jul 02 '19 at 13:20
  • @kwantuM: there are three effects combined, and the errors accumulate. –  Jul 02 '19 at 13:51
  • @YvesDaoust so how did 3b1b do this in python with N=350? – kwantuM Jul 02 '19 at 14:16
  • 1
    There is an error somewhere, in your code or in the numerical integrator (you don't say whether that's your own library). One possibility is that you're not actually calculating $c_n$ with as small an error as you think. – David C. Ullrich Jul 02 '19 at 14:50
  • 1
    Start with a simpler $f$, where the graph is just one line segment. So you can calculate the Fourier coefficients exactly, then compare with the values you get nummerically.. – David C. Ullrich Jul 02 '19 at 14:52
  • @DavidC.Ullrich I'm using https://docs.rs/quadrature/0.1.2/quadrature/double_exponential/fn.integrate.html. I could try to swap it out for another integrator. – kwantuM Jul 02 '19 at 14:52
  • But are you using it correctly? – David C. Ullrich Jul 02 '19 at 14:57
  • I suppose so? It just takes $f(t)$, boundaries and a target error, which I set to $1e-12$. – kwantuM Jul 02 '19 at 14:58
  • Btw, assuming that "it takes $f(t)$" means you pass the name of or a pointer to some function defined in whatever language, it seems clear to me that it's impossible for a numerical integration routine to guarantee a certain error, given that it just samples the function at finitely many points. (An integration routine could guarantee an error if you passed it information on the modulus of continuity of $f$...) – David C. Ullrich Jul 02 '19 at 15:12
  • If you want to find where the error is one thing you should do is start with a simpler $f$, where the graph is just one line segment, so you can calculate the Fourier coefficients exactly, then compare with the values you get numerically. You don't have to do that if you don't feel like it, but there's no point to this discussion until you do. (I missed the point to your question "What am I ignoring exactly?" - what you appear to be ignoring is my suggestion that you do this...) – David C. Ullrich Jul 02 '19 at 15:15
  • 1
    The code for the integrator comments "Assume only that error is decreasing". For an oscillatory integral there's no reason to think that the error is decreasing. – David C. Ullrich Jul 02 '19 at 15:22
  • 1
    What does it give for $\int_0^{2\pi}\cos(1000000t),dt$ with a target error of $0.1$? – David C. Ullrich Jul 02 '19 at 15:31
  • Same question for $\int_0^{2\pi} \cos(3628800t),dt$. – David C. Ullrich Jul 02 '19 at 15:34
  • Or maybe better, $\int_0^1\cos(3628800\pi t)$, with target_error=$0.1$. – David C. Ullrich Jul 02 '19 at 16:00

1 Answers1

1

This issue was caused by using a wrong integration algorithm. Just using https://en.wikipedia.org/wiki/Trapezoidal_rule with small steps (1e-4) MASSIVELY improves the results.

kwantuM
  • 151