pip install numpy scipy matplotlib

# lorenz_system.py

def lorenz_deriv(t, coords, sigma, rho, beta):
    """
    Defines the Lorenz system of differential equations.

    Args:
        t (float): Current time (required by solve_ivp, though not used in the Lorenz equations themselves).
        coords (list): A list or array containing the current [x, y, z] coordinates.
        sigma (float): The Prandtl number.
        rho (float): The Rayleigh number.
        beta (float): The geometric factor.

    Returns:
        list: The derivatives [dx/dt, dy/dt, dz/dt] at the current state.
    """
    x, y, z = coords
    dx_dt = sigma * (y - x)
    dy_dt = x * (rho - z) - y
    dz_dt = x * y - beta * z
    return [dx_dt, dy_dt, dz_dt]

if __name__ == '__main__':
    # Example of how the function works (not solving yet)
    current_coords = [0.0, 1.0, 1.05]
    sigma_val = 10
    rho_val = 28
    beta_val = 8/3

    derivatives = lorenz_deriv(0, current_coords, sigma_val, rho_val, beta_val)
    print(f"Current coordinates: {current_coords}")
    print(f"Derivatives [dx/dt, dy/dt, dz/dt]: {derivatives}")

python lorenz_system.py

# solve_lorenz.py
import numpy as np
from scipy.integrate import solve_ivp
from lorenz_system import lorenz_deriv # Import our derivative function

# --- Parameters for the Lorenz system ---
sigma_val = 10.0
rho_val = 28.0
beta_val = 8/3.0 # Ensure this is a float

# --- Initial conditions ---
# Small changes here will lead to vastly different results later!
initial_conditions = [0.0, 1.0, 1.05]

# --- Time span for the simulation ---
# (start_time, end_time)
t_span = (0, 40) # Simulate from t=0 to t=40

# --- Time points where we want to evaluate the solution ---
# We'll generate 10,000 points evenly spaced within our time span
t_eval = np.linspace(t_span[0], t_span[1], 10000)

# --- Solve the ODEs ---
# The args tuple passes our sigma, rho, beta to lorenz_deriv
print(f"Solving Lorenz system with initial conditions: {initial_conditions} and parameters: sigma={sigma_val}, rho={rho_val}, beta={beta_val}")
solution = solve_ivp(
    fun=lorenz_deriv,
    t_span=t_span,
    y0=initial_conditions,
    args=(sigma_val, rho_val, beta_val),
    t_eval=t_eval # Provide specific points to evaluate the solution
)

# --- Check if the solver was successful and extract results ---
if solution.success:
    print("\nSolver completed successfully!")
    # solution.y contains the computed states [x, y, z] for each time point
    # It's shaped as (num_dimensions, num_time_points)
    x_coords, y_coords, z_coords = solution.y

    print(f"Shape of solution.y: {solution.y.shape}")
    print(f"First 5 x values: {x_coords[:5]}")
    print(f"Last 5 x values: {x_coords[-5:]}")
else:
    print(f"\nSolver failed: {solution.message}")

python solve_lorenz.py

# plot_lorenz.py
import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D # Import for 3D plotting

# Import the derivative function from our lorenz_system.py file
from lorenz_system import lorenz_deriv

# --- Parameters for the Lorenz system ---
sigma_val = 10.0
rho_val = 28.0
beta_val = 8/3.0

# --- Initial conditions ---
initial_conditions = [0.0, 1.0, 1.05]

# --- Time span for the simulation ---
t_span = (0, 40)
t_eval = np.linspace(t_span[0], t_span[1], 10000) # 10,000 points for smoothness

# --- Solve the ODEs ---
solution = solve_ivp(
    fun=lorenz_deriv,
    t_span=t_span,
    y0=initial_conditions,
    args=(sigma_val, rho_val, beta_val),
    t_eval=t_eval
)

# --- Plotting the results ---
if solution.success:
    x_coords, y_coords, z_coords = solution.y

    fig = plt.figure(figsize=(10, 8))
    # Add a 3D subplot to the figure
    ax = fig.add_subplot(111, projection='3d')

    # Plot the trajectory
    ax.plot(x_coords, y_coords, z_coords, lw=0.5, color='blue', alpha=0.9)

    # Set labels and title
    ax.set_xlabel("X-axis")
    ax.set_ylabel("Y-axis")
    ax.set_zlabel("Z-axis")
    ax.set_title("The Lorenz Attractor")
    ax.grid(True)

    # You can rotate the view programmatically if needed
    # ax.view_init(elev=20, azim=120)

    plt.show()
else:
    print(f"Solver failed: {solution.message}")

python plot_lorenz.py

# butterfly_effect.py
import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

from lorenz_system import lorenz_deriv

# --- Parameters for the Lorenz system ---
sigma_val = 10.0
rho_val = 28.0
beta_val = 8/3.0

# --- Initial conditions ---
initial_conditions_1 = [0.0, 1.0, 1.05]
# Slightly perturb the Z coordinate by a tiny amount (e.g., 1 part in 100,000)
initial_conditions_2 = [0.0, 1.0, 1.05 + 1e-5]

# --- Time span for the simulation ---
t_span = (0, 40)
t_eval = np.linspace(t_span[0], t_span[1], 10000)

print(f"Simulating trajectory 1 with initial conditions: {initial_conditions_1}")
print(f"Simulating trajectory 2 with initial conditions: {initial_conditions_2}")

# --- Solve for the first set of initial conditions ---
sol1 = solve_ivp(lorenz_deriv, t_span, initial_conditions_1,
                 args=(sigma_val, rho_val, beta_val), t_eval=t_eval)

# --- Solve for the second set of initial conditions ---
sol2 = solve_ivp(lorenz_deriv, t_span, initial_conditions_2,
                 args=(sigma_val, rho_val, beta_val), t_eval=t_eval)

if sol1.success and sol2.success:
    x1, y1, z1 = sol1.y
    x2, y2, z2 = sol2.y

    # --- Plotting the two trajectories ---
    fig = plt.figure(figsize=(12, 10))
    ax = fig.add_subplot(111, projection='3d')

    ax.plot(x1, y1, z1, lw=0.7, color='blue', label='Trajectory 1: [0, 1, 1.05]')
    ax.plot(x2, y2, z2, lw=0.7, color='red', linestyle='--', label='Trajectory 2: [0, 1, 1.05 + 1e-5]')

    ax.set_xlabel("X-axis")
    ax.set_ylabel("Y-axis")
    ax.set_zlabel("Z-axis")
    ax.set_title("Lorenz Attractor: Sensitive Dependence on Initial Conditions")
    ax.legend()
    ax.grid(True)
    plt.show()

    # --- Plotting the Euclidean distance between the two trajectories over time ---
    distance = np.sqrt((x1 - x2)**2 + (y1 - y2)**2 + (z1 - z2)**2)
    
    plt.figure(figsize=(10, 6))
    plt.plot(t_eval, distance, color='purple')
    plt.xlabel("Time")
    plt.ylabel("Euclidean Distance Between Trajectories")
    plt.title("Divergence of Two Lorenz Trajectories Over Time")
    plt.grid(True)
    plt.yscale('log') # Use a logarithmic scale to show exponential divergence clearly
    plt.show()

else:
    print("One or both solvers failed.")
    if not sol1.success: print(f"Solver 1 failed: {sol1.message}")
    if not sol2.success: print(f"Solver 2 failed: {sol2.message}")

python butterfly_effect.py