Mathematics Asked by Gaganov Victor on November 1, 2021
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:
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…
Ok, after some digging I think I have it figured out, at least partially:
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)$
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:
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.
Answered by Gaganov Victor on November 1, 2021
Get help from others!
Recent Questions
Recent Answers
© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP