Python: Is there any way to get the n-th order antiderivative of a periodic 3D signal/field without padding data?

Question

As stated in the title, I want to obtain the n-th (e.g. 4-th) order antiderivative of a 3D field (e.g. array with shape (1024,1024,1024) ) with period L on each side. If I need the antiderivative values beyond period L, I have to pad the data periodically ("wrap") before solving the antiderivative. However, this approach is highly memory consumption, particularly for the 1024^3 array. My Python code snippet is shown below:

import numpy as np
from scipy.interpolate import make_interp_spline
data = np.load("xxx.npy") # shape: (1024,1024,1024)
data_pad = np.pad(data, ( (1024,1024), (1024,1024), (1024,1024) ), "wrap")

# 4th antiderivative of the padded data along x
interp_x = make_interp_spline(x_pad, data_pad, k=3, axis=0)
interp_x4 = interp_x.antiderivative(4)

Is there other way to achieve the same goal without padding data?

Answer 1

As I mentioned in the comments for periodic signals you can get antiderivatives from the Fourier transfrom.

If your function is (n+1) times differentiable, you can get n derivatives, after that you will have problems related to the Gibbs phenomenon

It is sad that we don't have equations support here

but using the translation identities you can get e.g. (f(x+a) - f(x-a)) / (2*a), as D(omega) * F(omega), and if D(omega) is non-zero everywhere, you can invert the derivation.

Let's start with a 1D example

import numpy as np

y = np.linspace(-4, 4, 256)
Y = np.fft.fft(y)
dt = np.exp(-np.fft.fftfreq(len(Y)) * 2j * np.pi)
dyOp = (dt**-0.5 - dt**0.5) * len(y) / (4 - (-4))
y_derivative = np.fft.ifft(dyOp * Y).real 
y_a1 = np.fft.ifft(np.divide(Y, dyOp, where=dyOp!=0)).real
y_a2 = np.fft.ifft(np.divide(Y, dyOp**2, where=dyOp!=0)).real
y_a3 = np.fft.ifft(np.divide(Y, dyOp**3, where=dyOp!=0)).real

plt.plot(y, y)
plt.plot(y, y_a1)
plt.plot(y, y_a2)
plt.plot(y, y_a3)
plt.ylim([-8, 8])
plt.legend(['y=t/4', '$d^{-1}ydt$', '$d^{-2}y dt^2$', '$d^{-3}y dt^3$'], ncol=2)

If your data has multiple dimensions you can apply this for for each axis individually

Assuming that your dimensions are 0 <= x,y,z < 1 with period 1

then you could write

period_x = 1
period_y = 1
period_z = 1

# used 1j pi, to avoid taking the square root down there
dx = np.exp(- 1j * np.pi * np.fft.fftfreq(data.shape[0])).reshape(-1, 1, 1)
dy = np.exp(- 1j * np.pi * np.fft.fftfreq(data.shape[1])).reshape(1, -1, 1)
dz = np.exp(- 1j * np.pi * np.fft.fftfreq(data.shape[2])).reshape(1, 1, -1)

data_dx = (dx.conj() - dx) * data.shape[0] / period_x
data_dy = (dx.conj() - dx) * data.shape[1] / period_y
data_dz = (dx.conj() - dx) * data.shape[2] / period_z

# maybe you will want to set some tolerance arround zero
data_dx[data_dx == 0] = 1
data_dy[data_dy == 0] = 1
data_dz[data_dz == 0] = 1

Suppose you want to solve the Poison's equation, you can do it by

D = np.fft.nfftn(data)
divergent_op = (data_dx**2 + data_dy**2 + data_dz**2)
data_divergent = np.fft.ifftn(D / divergent_op).real