[Python] カオスアトラクタのアルゴリズムを実装する方法

pip install numpy scipy matplotlib

import numpy as np
# 配列の作成
array = np.array([1, 2, 3, 4, 5])
print(array)
# 配列の要素の合計
total = np.sum(array)
print(total)

[1 2 3 4 5]
15

import matplotlib.pyplot as plt
# データの準備
x = [1, 2, 3, 4, 5]
y = [2, 3, 5, 7, 11]
# プロット
plt.plot(x, y)
plt.title("Sample Plot")
plt.xlabel("X-axis")
plt.ylabel("Y-axis")
plt.show()

from scipy.integrate import odeint
# 定義された方程式
def model(y, t):
    x, z = y
    dxdt = -x + z
    dzdt = -x - z
    return [dxdt, dzdt]
# 初期条件
y0 = [1, 0]
t = np.linspace(0, 10, 100)
# 数値解法
solution = odeint(model, y0, t)
print(solution)

def lorenz_equations(state, t, sigma, rho, beta):
    x, y, z = state
    dxdt = sigma * (y - x)
    dydt = x * (rho - z) - y
    dzdt = x * y - beta * z
    return [dxdt, dydt, dzdt]

def euler_method(lorenz_eq, initial_state, t, sigma, rho, beta):
    dt = t[1] - t[0]
    state = np.zeros((len(t), len(initial_state)))
    state[0] = initial_state
    for i in range(1, len(t)):
        state[i] = state[i - 1] + np.array(lorenz_eq(state[i - 1], t[i - 1], sigma, rho, beta)) * dt
    return state

def runge_kutta_method(lorenz_eq, initial_state, t, sigma, rho, beta):
    dt = t[1] - t[0]
    state = np.zeros((len(t), len(initial_state)))
    state[0] = initial_state
    for i in range(1, len(t)):
        k1 = np.array(lorenz_eq(state[i - 1], t[i - 1], sigma, rho, beta))
        k2 = np.array(lorenz_eq(state[i - 1] + 0.5 * k1 * dt, t[i - 1] + 0.5 * dt, sigma, rho, beta))
        k3 = np.array(lorenz_eq(state[i - 1] + 0.5 * k2 * dt, t[i - 1] + 0.5 * dt, sigma, rho, beta))
        k4 = np.array(lorenz_eq(state[i - 1] + k3 * dt, t[i - 1] + dt, sigma, rho, beta))
        state[i] = state[i - 1] + (k1 + 2 * k2 + 2 * k3 + k4) * dt / 6
    return state

import matplotlib.pyplot as plt
# パラメータの設定
sigma = 10.0
rho = 28.0
beta = 8/3
initial_state = [1.0, 1.0, 1.0]
t = np.linspace(0, 100, 10000)
# 数値解法の実行
state_euler = euler_method(lorenz_equations, initial_state, t, sigma, rho, beta)
state_rk = runge_kutta_method(lorenz_equations, initial_state, t, sigma, rho, beta)
# 軌道のプロット
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot(state_rk[:, 0], state_rk[:, 1], state_rk[:, 2], lw=0.5)
ax.set_title("Lorenz Attractor")
ax.set_xlabel("X-axis")
ax.set_ylabel("Y-axis")
ax.set_zlabel("Z-axis")
plt.show()

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
def lorenz_equations(state, t, sigma, rho, beta):
    x, y, z = state
    dxdt = sigma * (y - x)
    dydt = x * (rho - z) - y
    dzdt = x * y - beta * z
    return [dxdt, dydt, dzdt]
# パラメータの設定
sigma = 10.0
rho = 28.0
beta = 8/3
initial_state = [1.0, 1.0, 1.0]
t = np.linspace(0, 100, 10000)
# 数値解法の実行
state = odeint(lorenz_equations, initial_state, t, args=(sigma, rho, beta))
# 軌道のプロット
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot(state[:, 0], state[:, 1], state[:, 2], lw=0.5)
ax.set_title("Lorenz Attractor")
ax.set_xlabel("X-axis")
ax.set_ylabel("Y-axis")
ax.set_zlabel("Z-axis")
plt.show()

import numpy as np
import matplotlib.pyplot as plt
def henon_map(x, y, a=1.4, b=0.3):
    x_new = 1 - a * x**2 + y
    y_new = b * x
    return x_new, y_new
# 初期条件
x, y = 0, 0
iterations = 10000
x_values, y_values = [], []
for _ in range(iterations):
    x, y = henon_map(x, y)
    x_values.append(x)
    y_values.append(y)
# プロット
plt.figure(figsize=(8, 6))
plt.scatter(x_values, y_values, s=0.1, color='black')
plt.title("Henon Map")
plt.xlabel("X-axis")
plt.ylabel("Y-axis")
plt.axis('equal')
plt.show()

def logistic_map(r, x):
    return r * x * (1 - x)
# パラメータ設定
r_values = np.linspace(2.5, 4.0, 1000)
x = 0.5
iterations = 1000
last = 100
# プロットの準備
plt.figure(figsize=(10, 6))
for r in r_values:
    x = 0.5  # 初期条件
    for _ in range(iterations):
        x = logistic_map(r, x)
    # 最後の100回の値をプロット
    plt.plot([r] * last, [logistic_map(r, x) for _ in range(last)], ',k', alpha=0.25)
plt.title("Logistic Map Bifurcation Diagram")
plt.xlabel("r")
plt.ylabel("x")
plt.xlim(2.5, 4.0)
plt.ylim(0, 1)
plt.show()

def riccati_equation(t, y, a=1.0):
    return a + y**2
# 初期条件
y0 = 0.5
t = np.linspace(0, 10, 100)
# 数値解法
from scipy.integrate import solve_ivp
solution = solve_ivp(riccati_equation, [t[0], t[-1]], [y0], t_eval=t)
# プロット
plt.plot(solution.t, solution.y[0], label='Riccati Equation')
plt.title("Riccati Equation Solution")
plt.xlabel("Time")
plt.ylabel("y")
plt.legend()
plt.show()

def chua_circuit(state, t, params):
    x, y, z = state
    a, b, c = params
    dxdt = a * (y - x - (1/2) * (x**3 - x))
    dydt = x - y + z
    dzdt = -b * y
    return [dxdt, dydt, dzdt]
# パラメータ設定
params = [15.0, 28.0, 5.0]
initial_state = [0.7, 0.0, 0.0]
t = np.linspace(0, 100, 10000)
# 数値解法
solution = odeint(chua_circuit, initial_state, t, args=(params,))
# プロット
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot(solution[:, 0], solution[:, 1], solution[:, 2], lw=0.5)
ax.set_title("Chua Circuit")
ax.set_xlabel("X-axis")
ax.set_ylabel("Y-axis")
ax.set_zlabel("Z-axis")
plt.show()