超混沌 Lorenz 系统的四阶龙格-库塔法数值解

import numpy as np
import matplotlib.pyplot as plt

def lorenz(t, x, a, b, c, r):
    '超混沌 Lorenz 系统的微分方程'
    dx = a * (x[1] - x[0]) + x[3]
    dy = x[0] * c - x[1] - x[0] * x[2]
    dz = x[0] * x[1] - b * x[2]
    dw = -x[1] * x[2] + r * x[3]
    return np.array([dx, dy, dz, dw])

def runge_kutta_lorenz(a, b, c, r, x0, dt, num_steps):
    '使用四阶龙格-库塔法离散化 Lorenz 系统'
    t = np.zeros(num_steps + 1)
    x = np.zeros((num_steps + 1, 4))
    x[0] = x0

    for i in range(num_steps):
        t[i+1] = t[i] + dt

        k1 = lorenz(t[i], x[i], a, b, c, r)
        k2 = lorenz(t[i] + dt/2, x[i] + (dt/2) * k1, a, b, c, r)
        k3 = lorenz(t[i] + dt/2, x[i] + (dt/2) * k2, a, b, c, r)
        k4 = lorenz(t[i] + dt, x[i] + dt * k3, a, b, c, r)

        x[i+1] = x[i] + (dt/6) * (k1 + 2*k2 + 2*k3 + k4)

    return t, x

# 参数设置
a = 10
b = 8/3
c = 28
r = -1
x0 = np.array([0.1, 0, 0, 0])
dt = 0.01
num_steps = 10000

# 使用四阶龙格-库塔法离散化 Lorenz 系统
t, x = runge_kutta_lorenz(a, b, c, r, x0, dt, num_steps)

# 绘制轨迹
fig = plt.figure()
ax = fig.add_subplot(121, projection='3d')
ax.set_title('super Lorenz System - runge_kutta_lorenz', ha='center')
ax.plot(x[:, 0], x[:, 1], x[:, 2], lw=1)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')

ax = fig.add_subplot(122, projection='3d')
ax.set_title('super Lorenz System - runge_kutta_lorenz', ha='center')
ax.plot(x[:, 0], x[:, 1], x[:, 3], lw=1)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('w')

plt.show()

这段代码使用四阶龙格-库塔法数值求解超混沌 Lorenz 系统，并绘制其轨迹的三维图。代码中定义了 lorenz 函数用于计算 Lorenz 系统的微分方程，runge_kutta_lorenz 函数使用四阶龙格-库塔法求解微分方程。最终，代码使用 matplotlib 库绘制了系统轨迹在 (x, y, z) 和 (x, y, w) 空间中的三维图。

代码说明：

lorenz(t, x, a, b, c, r) 函数计算 Lorenz 系统的微分方程，其中 t 是时间，x 是状态向量，a, b, c, r 是系统参数。
runge_kutta_lorenz(a, b, c, r, x0, dt, num_steps) 函数使用四阶龙格-库塔法数值求解 Lorenz 系统，其中 x0 是初始状态向量，dt 是时间步长，num_steps 是时间步数。
代码最后使用 matplotlib 库绘制了系统轨迹的三维图。

如何运行代码：

确保已经安装了 numpy 和 matplotlib 库。
将代码保存为 .py 文件。
在命令行中运行 python your_file_name.py。

代码将生成包含两个三维图的窗口，分别显示系统轨迹在 (x, y, z) 和 (x, y, w) 空间中的变化情况。

代码优化：

使用了更具描述性的变量名。
在函数定义中添加了 docstring，以便于理解代码。
将代码组织得更加清晰易读。
使用 ha='center' 将图表的标题居中。

这段代码展示了如何使用 Python 求解超混沌 Lorenz 系统，并使用 matplotlib 库绘制其轨迹。您可以通过修改代码中的参数和初始条件来探索不同的 Lorenz 系统动力学。