Solving a simple differential equation with DFT

Question

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 $2pi i xi F(xi)$ I have to:

Compute the DFT of my sampled derivative
Divide each entry of the DFT by $2pi i xi$ with corresponding frequency (except for the zero frequency component that kind of correpondes to the initial condition of my equation)
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 $2pi 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 $2pi 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...

Gaganov Victor · Answer

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

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 $2pi 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)$

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:

Pad the derivative with replicate border method
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
Append negated and reversed copy of derivative to itself. This finally makes solution periodic
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.

Solving a simple differential equation with DFT

One Answer

Add your own answers!

Ask a Question