生成随机速度场时功率谱E(k)呈现k^(-(n-1))而非k^(-n)的原因及修正方案咨询

Hey there! I see exactly what's going on here—this is a super common pitfall when working with 2D power spectra, tied to mixing up total shell energy with the per-k energy density. Let's break this down clearly, then fix your code.

Why You're Seeing k^-(n-1) Instead of k^-n #

First, let's clarify two critical definitions that are getting mixed up:

1. Per-wavevector energy: You defined k_energy as the energy for each individual (kx, ky) pair in Fourier space, scaling with k_mag^(-n). That part is correct for your target spectrum.
2. Energy spectrum E(k): This is the energy per unit k interval (formally, dE/dk), so the total energy of the system is the integral of E(k) over all k from 0 to ∞.

The key issue is 2D k-space geometry: In 2D, the number of wavevectors (or the "volume element" of k-space) in a thin shell between k and k+dk scales linearly with k. Think of a ring around the origin—its circumference (and thus the number of discrete wavevectors it contains) is 2πk, so the count grows with k.

Right now, your code calculates E_k[i] as the total energy summed across all wavevectors in the shell. That total energy is:

Total shell energy ≈ (number of wavevectors in shell) × (energy per wavevector)
Total shell energy ∝ k × k^(-n) = k^(-n+1)

That's exactly the k^-(n-1) scaling you're seeing! You're plotting total shell energy instead of the normalized per-k energy density.

How to Fix the Code #

We need two key adjustments:

1. Correct FFT normalization: Numpy's FFT/IFT has implicit scaling that we need to account for to get accurate energy magnitudes.
2. Normalize by k-shell measure: Divide the total shell energy by the size of the shell (proportional to k × dk for 2D) to get the proper per-k energy spectrum E(k).

Here's the corrected version of your code, with comments highlighting changes:

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

# 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 with correct FFT normalization
# Numpy's ifft2 divides by n_grid², so we reverse that to get back to original ux_k
ux_k_recon = np.fft.fft2(ux) / n_grid**2
uy_k_recon = np.fft.fft2(uy) / n_grid**2
# Energy density per wavevector (matches our initial k_energy definition)
energy_density = 0.5 * (np.abs(ux_k_recon)**2 + np.abs(uy_k_recon)**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))
dk = k_bins[1] - k_bins[0]  # Width of each k-bin

for i in range(len(k_bin_centers)):
    mask = (k_mag >= k_bins[i]) & (k_mag < k_bins[i+1])
    total_shell_energy = np.sum(energy_density[mask])
    # Normalize by 2D k-shell measure (circumference × bin width) to get E(k) = dE/dk
    # We can omit the 2π constant since it won't affect the power-law scaling
    E_k[i] = total_shell_energy / (k_bin_centers[i] * dk)

# 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"Original result: k^(-{n-1})")
plt.xlabel("k")
plt.ylabel("E(k)")
plt.legend()
plt.show()

What This Does #

• The FFT normalization fixes the amplitude of the spectrum to match your initial k_energy definition.
• Dividing by k_bin_centers[i] * dk accounts for the 2D k-shell geometry, converting total shell energy into the per-k energy density E(k) that scales with k^-n as you expected.

When you run this corrected code, the computed E(k) should align perfectly with your target k^-n power law!

备注：内容来源于stack exchange，提问作者Mr. Who