陈氏超混沌序列生成的Python实现出现NaN、无穷值异常问题

• 核心逻辑偏差：当前代码混淆了连续动力系统的导数值与系统实际状态值。陈氏超混沌系统是一组连续常微分方程，方程中定义的xdot/ydot/zdot/qdot是状态变量对时间的导数，并非下一时刻的状态值，直接将导数赋值为新状态迭代，属于完全错误的离散逻辑，会导致数值快速指数发散，仅10步左右就会溢出为无穷值或NaN。
• 缺失数值求解环节：连续常微分方程必须搭配离散化的数值积分方法（如欧拉法、四阶龙格-库塔法）与合理的固定时间步长才能求解，现有代码完全跳过积分环节，实际构造了一个数值稳定性极差的离散映射，完全复现不了原超混沌系统的动力学特性。
• 系统敏感性适配缺失：超混沌系统对离散误差、步长大小极为敏感，没有稳定的积分逻辑加持时，哪怕极小的初始偏差都会被快速放大，最终触发数值溢出。

实现步骤 #

• 保留原方程的导数计算逻辑，不将导数直接作为状态值迭代
• 选用稳定性更强的四阶龙格-库塔（RK4）积分方法，相比一阶欧拉法能大幅降低离散累积误差
• 选取适配系统动力学的时间步长，常规区间为0.001~0.01，步长过大依然会触发数值发散
• 每一步的状态更新，基于当前状态、导数值计算积分增量推导，而非直接替换为导数
• 增加异常值检测，一旦出现数值溢出即可定位发散位置，及时调整步长或初始值

修正后可稳定运行的代码 #

import math

# 系统参数与你提供的配置保持一致
a = 36
b = 3
c = 28
d = 16
k = 0.2
# 新增适配连续ODE求解的固定时间步长
dt = 0.005

def chen_deriv(state):
    """输入当前状态元组，返回四个方向的导数值"""
    x, y, z, q = state
    xdot = a * (y - x)
    ydot = (-x * z) + (d * x) + (c * y) - q
    zdot = (x * y) - (b * z)
    qdot = x + k
    return (xdot, ydot, zdot, qdot)

def chaotic_seq(x0, y0, z0, q0, length):
    sequence = []
    current_state = (x0, y0, z0, q0)
    sequence.append(current_state)
    
    for i in range(length - 1):
        # 四阶龙格-库塔积分计算过程
        k1 = chen_deriv(current_state)
        mid_state1 = tuple(current_state[j] + 0.5 * dt * k1[j] for j in range(4))
        k2 = chen_deriv(mid_state1)
        mid_state2 = tuple(current_state[j] + 0.5 * dt * k2[j] for j in range(4))
        k3 = chen_deriv(mid_state2)
        guess_end_state = tuple(current_state[j] + dt * k3[j] for j in range(4))
        k4 = chen_deriv(guess_end_state)
        
        # 计算得到更新后的下一时刻状态
        next_state = tuple(
            current_state[j] + (dt / 6) * (k1[j] + 2 * k2[j] + 2 * k3[j] + k4[j])
            for j in range(4)
        )
        # 异常值检测
        for val in next_state:
            if math.isnan(val) or math.isinf(val):
                print(f"数值发散位置：第{i}步")
                return sequence
        sequence.append(next_state)
        current_state = next_state
    return sequence

if __name__ == '__main__':
    x0 = 0.3
    y0 = -0.4
    z0 = 1.2
    q0 = 1
    length = 2048
    seq = chaotic_seq(x0=x0, y0=y0, z0=z0, q0=q0, length=length)
    print(f"最终生成有效序列长度：{len(seq)}")

关键注意事项 #

• 步长调试：如果出现数值发散，优先将时间步长调小（比如从0.005降到0.002），不要随意更改系统参数
• 初始值范围：初始状态只要在系统吸引域内（x/y/z值取在[-20,20]区间，q值取在[-5,5]区间），都能稳定生成超混沌序列，无需过度调整初始值
• 序列归一化：如果需要将超混沌序列用于伪随机数生成，可取序列中的单分量（如x分量）做线性映射，将数值范围缩放到[0,1]区间后使用

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