为什么我的速度场的功率谱 E(k) 遵循 𝑘 ^(−(n−1)) 而不是 𝑘^(−n)？

Question

我正在生成一个随机速度场。通常，傅立叶空间中的能量定义为

E(k) = (vx(k)^2+vy(k)^2)/2

。在某些情况下，能量遵循幂律谱， 𝐸(𝑘)∼𝑘^(−𝑛) 在 2D 中，其中

k=np.sqrt(kx^2+ky^2)

。然而，当我根据生成的速度场计算功率谱时，我始终观察到 𝐸(𝑘)∼𝑘^(−(𝑛-1))。我怀疑这可能与维度或我计算功率谱的方式有关，但我不确定。

这是我的代码的简化版本：

import numpy as np
import matplotlib.pyplot as plt
import os

# Parameters
n_particles = 5000
n_grid = 1024
domain_size = 1.0
dt = 0.01
n = 1.4  # Target power-law exponent for energy spectrum, but you can put any value here

# Generate wave vectors and power-law energy spectrum
k = np.fft.fftfreq(n_grid, d=domain_size / n_grid)
kx, ky = np.meshgrid(k, k)
k_mag = np.sqrt(kx**2 + ky**2)
k_mag[0, 0] = 1.0  # Avoid division by zero
k_energy = k_mag**(-n)
k_energy[0, 0] = 0.0  # Zero out the mean component

# Generate velocity field in Fourier space with random phases
np.random.seed(0)
random_phase_x = np.exp(2j * np.pi * np.random.rand(n_grid, n_grid))
random_phase_y = np.exp(2j * np.pi * np.random.rand(n_grid, n_grid))
ux_k = random_phase_x * np.sqrt(k_energy)
uy_k = random_phase_y * np.sqrt(k_energy)

# Transform velocity field to real space
ux = np.fft.ifft2(ux_k).real
uy = np.fft.ifft2(uy_k).real

# Compute the power spectrum
energy_density = 0.5 * (np.abs(np.fft.fft2(ux))**2 + np.abs(np.fft.fft2(uy))**2)
kx = np.fft.fftfreq(n_grid).reshape(-1, 1)
ky = np.fft.fftfreq(n_grid).reshape(1, -1)
k_mag = np.sqrt(kx**2 + ky**2)

k_bins = np.arange(0, k_mag.max(), 0.01)
k_bin_centers = 0.5 * (k_bins[:-1] + k_bins[1:])
E_k = np.zeros(len(k_bin_centers))

for i in range(len(k_bin_centers)):
    mask = (k_mag >= k_bins[i]) & (k_mag < k_bins[i+1])
    E_k[i] = np.sum(energy_density[mask])

# Plot the computed spectrum
plt.loglog(k_bin_centers, E_k, label="Computed E(k)")
plt.loglog(k_bin_centers, k_bin_centers**(-(n)), '--', label=f"Expected: k^(-{n})")
plt.loglog(k_bin_centers, k_bin_centers**(-(n-1)), 'r.-', label=f"what I get: k^(-{n-1})")
plt.xlabel("k")
plt.ylabel("E(k)")
plt.legend()
plt.show()

为什么计算出的功率谱尺度为 𝐸(𝑘)∼𝑘^(−(𝑛-1)) 而不是 𝐸(𝑘)∼𝑘^(−𝑛)？如何修改我的代码或解释以确保速度场具有所需的幂律能谱？

Answer 1

在二维中，极坐标中的“面积”将是 k.dk.d（角度），而笛卡尔坐标中的“面积”是 dkx.dky。您实际上是在对角度进行积分，因此“面积”与周长乘以厚度成正比，即 k.dk。同样，你的面积与 k 成比例增加。因此，计算 E/k 作为“面积”能量密度。

如果您使用后者

E_k[i] = np.sum(energy_density[mask]) / ( 2 * k_bin_centers[i] )

然后你的图表就变成了