变分数阶Lorenz系统分岔图绘制Python代码验证与实现求助

问题背景 #

我编写了用于绘制变分数阶Lorenz系统分岔图的Python代码，采用分数阶Adams-Bashforth-Moulton(ABM)数值方法求解，初始条件为x0=[-2.5, 6.6, -15.0]，genOrds(start, end, step)用于生成阶数数组。目前已得到结果图，但不确定代码正确性，希望验证或获取正确实现。

现有代码问题分析 #

现有代码存在几个关键问题：

• 未定义Lorenz系统核心参数：sigma、rho、beta是Lorenz系统的必要参数，原代码未赋值，通常取标准值sigma=10、rho=28、beta=8/3。
• 分岔图点选取错误：直接绘制所有时间步的数值，未剔除暂态阶段的数据，会干扰分岔特征的观察。
• 依赖库未显式导入：代码使用了NumPy、SciPy、Matplotlib的功能，但未导入相关库。
• ABM算法细节偏差：部分系数计算的求和逻辑需调整，以匹配标准ABM预测-校正步骤。

修正后的完整代码 #

1. 导入依赖与定义基础参数 #

import numpy as np
from scipy.special import gamma
import matplotlib.pyplot as plt

# Lorenz系统标准参数
sigma = 10.0
rho = 28.0
beta = 8.0 / 3.0
x0 = [-2.5, 6.6, -15.0]

def genOrds(start, end, step):
    return np.arange(start, end + step, step)

2. 修正后的ABM数值求解实现 #

def x_fn(x, y, z):
    return sigma * (y - x)

def y_fn(x, y, z):
    return x * (rho - z) - y

def z_fn(x, y, z):
    return x * y - beta * z

def ABM(alpha, N=10000, h=0.004, s=x0):
    x = np.zeros(N + 1)
    y = np.zeros(N + 1)
    z = np.zeros(N + 1)
    x[0], y[0], z[0] = s

    # 预计算ABM算法系数
    b = np.zeros(N + 1)
    a = np.zeros(N + 1)
    for k in range(1, N + 1):
        b[k] = k**alpha - (k - 1)**alpha
        a[k] = (k + 1)**(alpha + 1) - 2 * k**(alpha + 1) + (k - 1)**(alpha + 1)

    factor_p = h**alpha / gamma(alpha + 1)
    factor_f = h**alpha / gamma(alpha + 2)

    for j in range(1, N + 1):
        # 预测步（Adams-Bashforth）
        sum_x = np.sum(b[1:j+1] * np.array([x_fn(x[k], y[k], z[k]) for k in range(j)]))
        sum_y = np.sum(b[1:j+1] * np.array([y_fn(x[k], y[k], z[k]) for k in range(j)]))
        sum_z = np.sum(b[1:j+1] * np.array([z_fn(x[k], y[k], z[k]) for k in range(j)]))
        
        x_p = x[0] + factor_p * sum_x
        y_p = y[0] + factor_p * sum_y
        z_p = z[0] + factor_p * sum_z

        # 校正步（Adams-Moulton）
        sum_a_x = np.sum(a[1:j] * np.array([x_fn(x[k], y[k], z[k]) for k in range(1, j)])) if j > 1 else 0
        sum_a_y = np.sum(a[1:j] * np.array([y_fn(x[k], y[k], z[k]) for k in range(1, j)])) if j > 1 else 0
        sum_a_z = np.sum(a[1:j] * np.array([z_fn(x[k], y[k], z[k]) for k in range(1, j)])) if j > 1 else 0
        
        term_x = ((j - 1)**(alpha + 1) - (j - 1 - alpha) * j**alpha) * x_fn(x[0], y[0], z[0])
        term_y = ((j - 1)**(alpha + 1) - (j - 1 - alpha) * j**alpha) * y_fn(x[0], y[0], z[0])
        term_z = ((j - 1)**(alpha + 1) - (j - 1 - alpha) * j**alpha) * z_fn(x[0], y[0], z[0])

        x[j] = x[0] + factor_f * (x_fn(x_p, y_p, z_p) + term_x + sum_a_x)
        y[j] = y[0] + factor_f * (y_fn(x_p, y_p, z_p) + term_y + sum_a_y)
        z[j] = z[0] + factor_f * (z_fn(x_p, y_p, z_p) + term_z + sum_a_z)

    return x, y, z

3. 分岔图数据生成与绘制 #

def get_bifurcation_data():
    values = []
    ords = genOrds(0.8, 1.0, 0.001)
    ord_list = []
    current_state = x0.copy()
    
    # 暂态步数：剔除前N_trans个点，只保留稳态数据
    N_trans = 2000
    N_total = 3000

    for alpha in ords:
        x, y, z = ABM(alpha, N=N_total, h=0.004, s=current_state)
        # 提取稳态阶段的Z值（可替换为X/Y值）
        steady_vals = z[N_trans:]
        values.extend(steady_vals)
        ord_list.extend([alpha] * len(steady_vals))
        # 更新初始状态为当前阶数的最终状态，加速收敛
        current_state = [x[-1], y[-1], z[-1]]
    
    return np.array(ord_list), np.array(values)

def plot_bifurcation(ords, values, title="变分数阶Lorenz系统分岔图", xlabel="导数阶数", ylabel="Z值"):
    plt.figure(figsize=(10, 6))
    plt.scatter(ords, values, color='darkblue', s=0.3, alpha=0.4)
    plt.title(title)
    plt.xlabel(xlabel)
    plt.ylabel(ylabel)
    plt.grid(alpha=0.3)
    plt.show()

# 生成并绘制分岔图
ords, z_values = get_bifurcation_data()
plot_bifurcation(ords, z_values)

关键修正说明 #

1. 补充核心参数：添加Lorenz系统标准参数，确保系统行为符合经典定义。
2. 剔除暂态数据：设置N_trans=2000，只保留系统达到稳态后的点，避免暂态数据干扰分岔特征。
3. 优化ABM算法：调整求和逻辑，严格匹配ABM预测-校正步骤，提升数值计算精度。
4. 状态延续：将前一个阶数的最终状态作为下一个阶数的初始值，加速系统收敛到稳态。

预期结果 #

修正后的代码生成的分岔图会清晰展示：当导数阶数从0.8上升到1.0时，系统从周期态逐渐过渡到混沌态，阶数接近1时呈现经典Lorenz系统的混沌特征，与理论预期一致。

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