Python模拟微分Jerk系统及多涡卷混沌吸引器空图问题求助

核心问题分析与修复 #

• 初始状态选择错误：混沌系统对初始条件极度敏感，全零初始状态会让系统陷入不动点（所有状态变量保持为0），无法激发混沌行为。需要给初始状态添加微小扰动，比如：

  initial_state = [0.1, 0.0, 0.0, 0.0, 0.0, 0.0]  # 给x1一个微小非零值
  

• 缺失绘图代码：你只完成了ODE求解，没有实现绘图逻辑。添加以下代码生成相图（以x₁ vs x₃为例，这是混沌吸引器常用的可视化维度）：

  import matplotlib.pyplot as plt
  
  # 提取状态变量
  x1 = solution[:, 0]
  x3 = solution[:, 2]
  
  # 绘制相图
  plt.figure(figsize=(8,6))
  plt.plot(x1, x3, linewidth=0.5)
  plt.xlabel('x₁')
  plt.ylabel('x₃')
  plt.title('Multi-scroll Chaotic Attractor')
  plt.show()
  

• 非线性函数f(x)验证：检查函数计算是否符合论文定义。可以单独测试f(x)的输出，比如：

  print(f(0))  # 验证是否符合论文预期值
  

  若参数与论文3.2(i)不符，需调整M、e、q的取值。

完整修复后的代码示例 #

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt

def system(state, t, alpha):
    x1, x2, x3, x4, x5, x6 = state
    dxdt = [
        x2,
        x3,
        x4,
        x5,
        x6,
        alpha[0] * (f(x1) - x1) - alpha[1] * x2 - alpha[2] * x3 - alpha[3] * x4 - alpha[4] * x5 - alpha[5] * x6
    ]
    return dxdt

def f(x):
    M = 1
    e = 1
    q = 0.01
    return (e / (2 * q)) * sum([np.abs(x - 2 * m * e + q) - np.abs(x - 2 * m * e - q) for m in range(-M, M + 1)])

# 参数设置
alpha_values = [7.5, 10, 10, 40, 10, 10]  
initial_state = [0.1, 0.0, 0.0, 0.0, 0.0, 0.0]  # 添加初始扰动
time = np.linspace(0, 100, 10000) 

# 求解ODE
solution = odeint(system, initial_state, time, args=(alpha_values,))

# 绘制相图
x1 = solution[:, 0]
x3 = solution[:, 2]
plt.figure(figsize=(8,6))
plt.plot(x1, x3, linewidth=0.5, color='darkblue')
plt.xlabel('x₁')
plt.ylabel('x₃')
plt.title('Multi-scroll Chaotic Attractor')
plt.show()

Jerk系统是三阶微分方程，通用形式为：

$x''' + a x'' + b x' + f(x) = 0$

模拟时需将其转换为一阶微分方程组，步骤如下：

1. 变量替换：令

   • $x_1 = x$
   • $x_2 = x'$
   • $x_3 = x''$

   转换后的一阶方程组为：
   $$
   \begin{cases}
   \dot{x_1} = x_2 \
   \dot{x_2} = x_3 \
   \dot{x_3} = -a x_3 - b x_2 - f(x_1)
   \end{cases}
   $$

2. 示例实现（经典分段线性Jerk系统）

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt

def jerk_system(state, t, a, b):
    x1, x2, x3 = state
    # 分段线性非线性函数f(x)
    def f(x):
        if x > 1:
            return 1
        elif x < -1:
            return -1
        else:
            return x
    dxdt = [
        x2,
        x3,
        -a*x3 - b*x2 - f(x1)
    ]
    return dxdt

# 参数设置（经典Jerk系统参数）
a = 0.5
b = 0.7
initial_state = [0.1, 0.0, 0.0]  # 同样需要非零初始扰动
time = np.linspace(0, 200, 20000)

# 求解ODE
solution = odeint(jerk_system, initial_state, time, args=(a, b))

# 绘制相图（x1 vs x2）
x1 = solution[:, 0]
x2 = solution[:, 1]
plt.figure(figsize=(8,6))
plt.plot(x1, x2, linewidth=0.5, color='darkred')
plt.xlabel('x')
plt.ylabel('x\'')
plt.title('Jerk System Chaotic Attractor')
plt.show()

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