Python实现非线性弹簧力Leapfrog积分遇double_scalars溢出问题求助

我正在求解含非线性项α的简谐运动（SHM）方程，非线性项被加入弹簧力中。已编写采用Leapfrog积分法的Python代码求解该非线性振子方程的多个振荡周期，但运行时出现数值溢出和无效值警告，希望得到解决思路。

代码实现 #

import matplotlib.pyplot as plt
import numpy

k      = 1.0 # Spring constant
m      = 1.0 # Mass
cycles = 2.0 # No. of periods to integrate over
alpha  = 0.5 # Nonlinear Parameter
x0     = 1.0 # Initial displacement
v0     = 0.0 # Initial velocity

def leapfrog( steps ):
    """Solve the nonlinear oscillator equation for several 
    oscillation cycles"""
    omega = (k/m)**0.5
    delta = 2.0*cycles*numpy.pi/omega/steps
    x     = numpy.empty( steps+1 )
    v     = numpy.empty( steps+1 )
    t     = numpy.empty( steps+1 )
    t[0]  = 0.0
    x[0]  = x0
    v[0]  = v0 + 0.5*delta*x[0]*(omega**2*(1+alpha*x[0])*x[0])
    for i in range(steps):
        t[i+1] = t[i] + delta
        x[i+1] = x[i] + delta*v[i]
        v[i+1] = v[i] - delta*(omega**2*(1+alpha*x[i+1])*x[i+1])
    return t, x, v

def manufactured_solution(t, i):
    omega = numpy.sqrt(k/m)
    M = alpha*omega**2*x0**2*numpy.cos(omega*t[i])**2
    x_exact = x0*numpy.cos(omega*t[i]) + M/omega**2
    return x_exact

def l2_error_norm( t , x ):
    """Calculate the L2 relative error norm."""
    steps  = len( x ) - 1
    omega  = (k/m)**0.5
    l2_err = 0.0
    l2     = 0.0
    for i in range(steps):
        x_exact = manufactured_solution(t, i)
        l2_err += (x[i] - x_exact)**2.0
        l2     += x[i]**2.0
    return (l2_err/l2)**0.5

n = 14
steps = 8                                                    
delta    = numpy.empty( n )
l2_error = numpy.empty( n )

for i in range(0,n):
    t, x, v     = leapfrog( steps )
    delta[i]    = (k/m)**0.5*(t[1]-t[0])
    l2_error[i] = l2_error_norm( t , x )
    plt.plot( t , x )
    # plt.show()
    steps *= 2

# Switch to a new plotting window, and plot the L2 error norm,
# with guidelines for first, second, and third order accuracy.
plt.figure()
plt.loglog( delta , l2_error , 'o' )
plt.loglog( delta , l2_error[0]*(delta/delta[0])**1.0 )
plt.loglog( delta , l2_error[0]*(delta/delta[0])**2.0 )
plt.loglog( delta , l2_error[0]*(delta/delta[0])**3.0 )
plt.show()

运行警告信息 #

nhm.py:38: RuntimeWarning: overflow encountered in double_scalars
v[i+1] = v[i] - delta*(omega**2*(1+alpha*x[i+1])*x[i+1])
nhm.py:49: RuntimeWarning: overflow encountered in double_scalars
l2_err += (x[i] - x_exact)**2.0
nhm.py:50: RuntimeWarning: overflow encountered in double_scalars
l2 += x[i]**2.0
nhm.py:51: RuntimeWarning: invalid value encountered in double_scalars return (l2_err/l2)**0.5

解决思路 #

• 修正Leapfrog初始速度计算：原代码初始半步速度公式缺少负号，加速度公式应为 (a = -\omega^2(1+\alpha x)x)，正确的初始速度更新应为：

  v[0] = v0 + 0.5 * delta * (-omega**2 * (1 + alpha*x[0]) * x[0])
  

  符号错误会导致位移和速度不断累积增大，最终引发数值溢出。

• 验证人工制造解的合理性：当前manufactured_solution的形式是否满足非线性振子方程？非线性方程的解析解通常并非该形式，需代入原方程验证一致性，否则误差计算会出现无意义的偏差，甚至在积分发散后加剧溢出问题。

• 限制积分步数增长范围：循环中steps *= 2执行14次后，步数会达到65536。若小步长阶段已出现积分发散，后续计算会累积错误。可先测试前5次循环，确认小步长下积分稳定后再逐步增加次数。

• 添加数值稳定性检查：在Leapfrog循环中加入阈值判断，当位移或速度超过合理范围（如1e10）时终止循环并输出警告，避免溢出扩散：

  for i in range(steps):
      t[i+1] = t[i] + delta
      x[i+1] = x[i] + delta*v[i]
      if abs(x[i+1]) > 1e10:
          print(f"Overflow detected at step {i+1}, x = {x[i+1]}")
          break
      v[i+1] = v[i] - delta*(omega**2*(1+alpha*x[i+1])*x[i+1])
  

• 调整数值精度或步长：若初始问题修复后仍有溢出，可尝试使用numpy.float128等更高精度的浮点数类型，或调整初始步长，避免步长过大导致积分不稳定。

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