海洋表面波建模中Stokes波代码实现结果与维基百科不符的bug排查

代码核心问题 #

• 混淆了Stokes波色散关系修正项与波面升高展开项：你在基波余弦项前增加的(1-1/16*(k*a)**2)系数属于相速度的修正项，不属于波面升高的三阶展开内容，该额外系数会压低基波幅值，导致高次谐波占比异常，最终波形不符合预期。
• 小参数适用范围注意：Stokes波是基于ka<<1的小参数展开，当ka>0.4时展开会逐渐偏离真实波形，你测试用的a最大值为0.2，对应k=2π时ka≈1.256，已经远超小参数假设的适用范围，也会导致波形失真。

修正后代码 #

import numpy as np
import matplotlib.pyplot as plt

# 补充原代码缺失的通用绘图参数定义
phi = (1 + np.sqrt(5)) / 2
fig_width = 8

def stokes(pos, a=.3/(2*np.pi), k=2*np.pi):
    # ka为无量纲小参数，需满足ka << 1
    elevation = np.cos(k*pos)  # 基波无额外修正系数
    elevation += 1/2 * (k*a) * np.cos(2*k*pos)  # 二阶谐波项
    elevation += 3/8 * (k*a)**2 * np.cos(3*k*pos)  # 三阶谐波项
    return a * elevation

N_pos = 501
pos = np.linspace(0, 1, N_pos, endpoint=True)    
fig, ax = plt.subplots(figsize=(fig_width, fig_width/phi**3))
ax.plot(pos, stokes(pos))
ax.set_xlim(np.min(pos), np.max(pos))
plt.show()

fig, ax = plt.subplots(figsize=(fig_width, fig_width/phi**3))
# 调整a的范围，保证ka最大值不超过0.4，符合小参数假设
for a in np.linspace(.001, 0.4/(2*np.pi), 10, endpoint=True):
    elevation = stokes(pos, a=a)
    ax.plot(pos, elevation)
ax.set_xlim(np.min(pos), np.max(pos))
plt.show()

修改后生成的波形将符合Stokes波“波峰尖锐、波谷平缓”的典型特征，当小参数ka控制在合理范围内时，和标准Stokes波参考波形完全一致。

内容的提问来源于stack exchange，提问作者meduz