1
$\begingroup$

I have a simple differential equation of form $\frac{df(x)}{dx}=u(x)$ and u(x) is given to me in a sampled form so I need to solve this equation numerically. I decided to do this via DFT and I came up with following plan:

Since Fourier Transform of $\frac{df(x)}{dx}$ is given by $2\pi i \xi F(\xi)$ I have to:

  1. Compute the DFT of my sampled derivative
  2. Divide each entry of the DFT by $2\pi i \xi$ with corresponding frequency (except for the zero frequency component that kind of correpondes to the initial condition of my equation)
  3. Compute the inverse DFT of the resulting signal and this should be the solution

I have coded this up with Python and numpy and implemented a simple unit test based on a function with known derivative (I used $sinc(x)$). The resulting solution differs alot from $sinc(x)$ and I can't figure out what am I doing wrong here. Here is the code of all this with visualization using matplotlib that illustrates the difference between target function and the result and it's not even close...

import numpy as np
from numpy.fft import fft, ifft
import matplotlib.pyplot as plt
import cmath as cm
import math

def integrate(f, fdx, a, b):
    N = 1024
    X = np.linspace(a, b, N)
    step = (b - a) / (N - 1)
    D = [step * fdx(x) for x in X]
    S = [f(x) for x in X]

    ft = fft(D)
    w = [(2 * math.pi * 1j * i / N) for i in range(0, N)]
    w[0] = 1    # dc component, ignore it
    ft = ft / w
    res = np.real(ifft(ft))

    plt.figure()
    plt.plot(X, S, 'g--', linewidth = 2.0)
    plt.plot(X, res, 'r-', linewidth = 1.0)
    plt.show()

f = lambda x : math.sin(x) / x if math.fabs(x) > 1e-10 else 1
fdx = lambda x : (x * math.cos(x) - math.sin(x)) / (x ** 2) if math.fabs(x) > 1e-10 else 0
integrate(f, fdx, 0, 25)

Do you have any idea on what's wrong here?

UPDATE: Here is the plot of sinc function and numeric solution obtained by this method. There are two subplots displaying the same thing with different ylimit

Fig.1: Visualization of described experiment

Another interesting thing to observe here is the amplitude of the signal spectrum before and after the division by $2\pi i \xi$ is applied. Here is the corresponding plot of $log$ amplitudes of the spectrum:

Fig.2: The spectrum amplitude before and after the procedure

It can be seen, that before the procedure spectrum is symmetic (and that's cool cause the signal that I start with is real valued), however after the procedure is applied, the spectrum turns non symmetric and the inverse DFT of the signal isn't even real anymore, it's complex. So there is clearly some flaw in the procedure but I can't see where it is. The only idea that I have in mind is that the continuous formula $2\pi i \xi F(\xi)$ for the FT of first order derivative doesn't hold in the discrete case, but I can't see why and what are the alternatives...

$\endgroup$
8
  • $\begingroup$ Can you attach the figure as well? $\endgroup$
    – user619894
    Commented Jul 21, 2020 at 8:49
  • $\begingroup$ @user619894 Sure, I have attached the resulting figure to the original post. I hope this would help! $\endgroup$ Commented Jul 21, 2020 at 9:36
  • $\begingroup$ shouldn't you multiply by $w$, and not divide? $\endgroup$
    – user619894
    Commented Jul 21, 2020 at 9:55
  • $\begingroup$ @user619894 Multiplying by $w$ doesn't seem to be the right thing. After I apply the Fourier Transform to both sides of the equation I get $2\pi i \xi F(\xi)= \hat u(\xi)$ so $F(\xi)=\frac{\hat u(\xi)}{2\pi i \xi}$. However I tried to do what you propose and it doesn't give the correct result as well $\endgroup$ Commented Jul 21, 2020 at 10:12
  • 1
    $\begingroup$ I think that the numpy fft has the positive freqs up to $N/2$ and then negative freqs, so the $w$ vector should reflect that. $\endgroup$
    – user619894
    Commented Jul 21, 2020 at 12:53

1 Answer 1

1
$\begingroup$

Ok, after some digging I think I have it figured out, at least partially:

  1. There was a bug in the implementation that was kindly noted by user619894 and that is numpy FFT frequencies are not $[0,...,N - 1]$ but rather [0,..,N/2,-N/2,..,-1]. The most reliable way to implement division by $2\pi i \xi$ in numpy turns out to be using it's builtin numpy.fftfreq(N) function, this func returns an array of frequencies, corresponding to entries of FFT and with it the division procedure can be implemented in one line by:
ft = ft / np.append([1], (2 * math.pi * 1j * fftfreq(N)[1:]))

However this fix alone doesn't make the thing work for $sinc(x)$

  1. The second but probably more shocking thing turns out to be the pereodic nature DTF based integration procedure. After staring at the result of DFT integration of my derivative I noticed that the resulting solution kind of resambles the $sinc(x)$ in it's shape but is bent so that it converges to the same value that it starts with. After that I have realised that DFT actually represents an endless periodic signal and the last derivative in my sampled array of derivatives actually correspondes to the difference between samples $N - 1$ and $0$ since the thing is periodic. To veryfy this I replaced the last derivative by the corresponding difference and the results became much more like the $sinc(x)$ after that, but with minor high frequency distortions near the begining and the end of the signal.

After that I have conducted multiple experiments with different signals weighted by various windowing fucntions that zero out the signal near the end and the beginning and the procedure works perfect for such signals and produces nice results even in presence of noise.

So it seems that it's the periodizing nature of descrete version of the transform that doesn't allow to solve for such a function like $sinc(x)$ in the described setup. I don't see any simple way around it. Does anybody know how to work around this or the DFT is simply not apropriate for such a non-periodic problem?

UPDATE: In the end I used the following trick to do the integration with FFT here, so it turns out to be doable dispite periodicity constraint implicitly imposed by DFT:

  1. Pad the derivative with replicate border method
  2. Multiply your padding with some windowing function on both sides of the signal. This makes solution to augmented equation flatten at the beginning and at the end. But the solution is still not periodic
  3. Append negated and reversed copy of derivative to itself. This finally makes solution periodic
  4. Solve with FFT, then extract the portion of augmented signal that correspondes to the original part of the signal

Works perfect, although the code is kinda ugly and not as cool as an initial one-liner:

    half = int(N / 2) # padding size
    tails = np.hanning(2 * half) * np.array(([D[0]] * half) + ([D[-1]] * half))
    SP = np.hstack((tails[:half], D, tails[half:])) # flatten near boundaries
    SPP = np.append(SP, -np.flip(SP)); L = SPP.shape[0] # periodize
    res = np.real(ifft(fft(SPP) / np.append([1], (2 * math.pi * 1j * fftfreq(L)[1:]))))
    res = res[half:half + N] # extract original signal from augmented signal

Also this quadruples the size of FFT. Maybe a smaller padding can be used, but anyway this solution at least doubles the size of FFT.


$\endgroup$
3
  • $\begingroup$ Great catch. To solve the periodic issue, either increase the range, or multiply by a window function en.wikipedia.org/wiki/Window_function. $\endgroup$
    – user619894
    Commented Jul 22, 2020 at 5:09
  • 1
    $\begingroup$ An alternative is to detrend: decompose the signal into a linear term that is just a line between the two end points, and the original signal with this line subtracted. The subtracted signal is now periodic. You can FT the subtracted signal, divide by $w$ and iFT, then adding the integral of the linear signal as well. $\endgroup$
    – user619894
    Commented Jul 22, 2020 at 7:25
  • $\begingroup$ @user619894 Nice ideas, thanks! I padded + windowed the derivative to avoid high freq jitter and then mirrored and negated it to make the thing periodic. Works like a charm. $\endgroup$ Commented Jul 25, 2020 at 12:36

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .