### Python 实现不依赖 SciPy 的龙格库塔法数值求解分数阶蔡氏混沌系统 #### 方法概述 分数阶蔡氏混沌系统的数值求解可以通过扩展经典的四阶龙格-库塔方法来完成。由于分数阶微分方程涉及历史记忆效应，通常采用 Grünwald-Letnikov 定义或其他离散化技术将其转化为差分形式[^1]。在此基础上，可以设计适用于分数阶系统的改进型 Runge-Kutta 方法。 以下是基于纯 NumPy 和 Matplotlib 的实现方案： --- #### 代码实现 ```python import numpy as np import matplotlib.pyplot as plt def fractional_runge_kutta(alpha, f, y0, t_span, h): """ 使用改进的分数阶龙格-库塔方法求数值解。 参数: alpha (float): 分数阶导数阶次 f (function): 微分方程右端项函数 y0 (array-like): 初始条件 t_span (tuple): 时间区间 [t_start, t_end] h (float): 步长 返回: ts (list): 时间序列 ys (list of lists): 解的时间演化 """ t_start, t_end = t_span n_steps = int((t_end - t_start) / h) ts = np.linspace(t_start, t_end, n_steps + 1) ys = np.zeros((n_steps + 1, len(y0))) ys[0] = y0 # 计算权重系数（Grünwald-Letnikov） weights = [] w_prev = 1.0 for k in range(n_steps + 1): weight = (-alpha * w_prev) / k if k != 0 else 1.0 weights.append(weight) w_prev = weight for i in range(1, n_steps + 1): sum_term = 0 for j in range(i): sum_term += weights[i-j-1] * f(ts[j], ys[j]) ys[i] = ys[i-1] + h**alpha * sum_term / np.math.gamma(2-alpha) return ts, ys # 蔡氏混沌系统的定义 def chua_system(t, state, a=15.6, b=28, c=-7/5): x, y, z = state g_x = -(c/7)*x + (c/2)*(abs(x+1)-abs(x-1)) dx_dt = a*(y-x-g_x) dy_dt = x-y+z dz_dt = -b*y-z return np.array([dx_dt, dy_dt, dz_dt]) # 参数设置 initial_state = [0.1, 0.1, 0.1] time_interval = (0, 100) step_size = 0.01 fractional_order = 0.95 # 求解 ts, states = fractional_runge_kutta(fractional_order, lambda t, s: chua_system(t, s), initial_state, time_interval, step_size) # 提取变量轨迹 xs, ys, zs = states[:, 0], states[:, 1], states[:, 2] # 绘制时间序列图 plt.figure(figsize=(12, 6)) plt.plot(ts, xs, label="x(t)", color='blue') plt.plot(ts, ys, label="y(t)", color='green', linestyle="--") plt.plot(ts, zs, label="z(t)", color='red', linestyle="-.") plt.title("Chua System Time Evolution", fontsize=16) plt.xlabel("Time", fontsize=14) plt.ylabel("State Variables", fontsize=14) plt.legend() plt.grid(True) plt.show() # 绘制相空间图 fig = plt.figure(figsize=(12, 8)) ax = fig.add_subplot(projection='3d') ax.plot(xs, ys, zs, lw=0.5, color='purple') ax.set_title("Phase Space Trajectory of Chua System", fontsize=16) ax.set_xlabel('X Axis', fontsize=14) ax.set_ylabel('Y Axis', fontsize=14) ax.set_zlabel('Z Axis', fontsize=14) plt.show() # 分叉图分析 param_values = np.linspace(10, 20, 500) # 可变参数范围 last_points = [] for param in param_values: def modified_chua_system(t, state): return chua_system(t, state, a=param) _, temp_states = fractional_runge_kutta(fractional_order, lambda t, s: modified_chua_system(t, s), initial_state, time_interval, step_size) last_point = temp_states[-1][0] # 获取最终状态的第一个坐标 last_points.append(last_point) plt.figure(figsize=(12, 6)) plt.scatter(param_values, last_points, s=1, color='black') plt.title("Bifurcation Diagram of Chua System", fontsize=16) plt.xlabel("Parameter 'a'", fontsize=14) plt.ylabel("Last Point of X", fontsize=14) plt.grid(True) plt.show() ``` --- #### 结果解释 上述代码实现了以下功能： 1. **分数阶龙格-库塔方法**：通过引入 Grünwald-Letnikov 权重计算方式，适应于分数阶微分方程的数值积分。 2. **蔡氏混沌系统建模**：利用非线性函数 `g_x` 描述电路特性，并构建三维动力学模型。 3. **几何行为分析**：绘制了时间序列图和三维相空间轨迹，直观展示系统动态特征。 4. **分叉图生成**：改变控制参数 `a` 并记录稳态响应，揭示系统随参数变化的行为模式。 --- #### 注意事项 - 上述实现假设初始条件已知且步长适当选取以保证精度。 - 如果需要更高的效率或更复杂的场景支持，则可能需进一步优化算法结构或并行化处理。 - 对于某些极端情况下的稳定性验证，建议结合理论推导补充说明。 ---