解决静电问题傅里叶变换中的零频率问题

Question

“我正在寻求利用离散傅里叶变换 (DFT) 来解决静电问题，但我遇到了零频率问题。为了说明这一点，我尝试了一个基本问题。如果您能指导我如何解决处理除法并获得正弦函数，这是该方程的正确解，它会减轻我肩上的负担。”

限制为零

当我做这样的事情时：

$$ \\frac{d\\phi(x)}{dx} = cos(x) $$
which solution is $\\phi(x) = sin(x)$ 

如果将限制设为 $k$ 接近零，则第一个频率将变为零。因此，您可以使用零作为第一个值，这样就可以解决所有问题。

import numpy as np
import matplotlib.pyplot as plt

# Define the parameters

num_points = 10000  # Number of points in the signal
delta_x = 0.01      # Spacing between points in the space domain
x = np.linspace(0, (num_points - 1) * delta_x, num_points)

# Define the cosine function of x

y = np.cos(x)

# Compute the Fourier transform of the cosine of x

Y = np.fft.rfft(y)

# Compute the corresponding frequencies

k_values = 2*np.pi * np.fft.rfftfreq(num_points, delta_x)

# Apply a mask to avoid the first value in the frequencies

k_values_masked = k_values[1:]  # Avoid the first value

# Compute the Fourier transform of the cosine of x divided by ik, handling the limit at k = 0

masked_Y_divided = np.ones_like(Y, dtype=np.complex128)
masked_Y_divided[1:] = Y\[1:] / (1j * k_values_masked)

# Perform the inverse Fourier transform to obtain the solution

reconstructed_y = np.fft.irfft(masked_Y_divided)

# Plot the result

plt.figure(figsize=(10, 6))

plt.subplot(2, 1, 1)
plt.plot(x, y, label='Cosine of x')
plt.title('Cosine of x')
plt.xlabel('x')
plt.ylabel('y')
plt.grid(True)
plt.legend()

plt.subplot(2, 1, 2)
plt.plot(x, reconstructed_y.real, label='Reconstructed sine of x')
plt.plot(x, np.sin(x), label='Expected result')
plt.title('Reconstructed sine of x')
plt.xlabel('x')
plt.ylabel('y')
plt.grid(True)
plt.legend()

plt.tight_layout()
plt.show()

我希望获得正弦函数，但我遇到了一般情况的问题，特别是当我不确定答案并且不知道为什么振幅的解也会受到影响时。此外，我获得的解决方案与应有的解决方案相比有所倾斜。

Answer 1

您的一些问题源于 FFT 长度选择不当。大胆猜测一下，这更接近您的期望：

import numpy as np
import matplotlib.pyplot as plt

# Number of points in the signal, multiple of the cosine's angular velocity
num_points = int(round(1e3 * np.pi))
delta_x = 0.01       # Spacing between points in the space domain
x = np.linspace(start=0, stop=(num_points - 1)*delta_x, num=num_points)
y = np.cos(x)
Y = np.fft.rfft(y)

# Actually the angular velocities
omega = 2*np.pi * np.fft.rfftfreq(n=num_points, d=delta_x)

# Avoid the first value in the frequencies
Y[0] = 0
omega[0] = 0.1

# Compute the Fourier transform of the cosine of x divided by ik, handling the limit at k = 0
masked_Y_divided = Y / (1j * omega)

# Perform the inverse Fourier transform to obtain the solution
reconstructed_y = np.fft.irfft(masked_Y_divided)

fig, (ax_t, ax_f) = plt.subplots(ncols=2)
ax_t.plot(x, y, label='Cosine of x')
ax_t.plot(x, reconstructed_y.real, label='Reconstructed sine of x')
ax_t.set_xlabel('x')
ax_t.set_ylabel('y')
ax_t.grid(True)
ax_t.legend()
ax_f.plot(omega[:40], np.abs(Y[:40]))
ax_f.set_xlabel('omega')
ax_f.set_ylabel('Y')

plt.show()