LORENZ

Overview

The LORENZ function numerically solves the Lorenz system, a set of three coupled ordinary differential equations (ODEs) that serve as a foundational model in chaos theory. Developed by meteorologist Edward Lorenz in 1963, this system was originally designed as a simplified model for atmospheric convection and famously gave rise to the concept of the butterfly effect—the idea that small changes in initial conditions can lead to vastly different outcomes.

The system describes the evolution of three state variables over time:

\frac{dx}{dt} = \sigma(y - x)

\frac{dy}{dt} = x(\rho - z) - y

\frac{dz}{dt} = xy - \beta z

The three parameters control the dynamics: σ (sigma) is the Prandtl number representing the ratio of fluid viscosity to thermal conductivity, ρ (rho) is the Rayleigh number representing the temperature difference driving convection, and β (beta) is a geometric factor related to the aspect ratio of convection cells. With Lorenz’s classic parameter values (σ = 10, ρ = 28, β = 8/3), the system exhibits chaotic behavior and trajectories converge to the famous butterfly-shaped Lorenz attractor.

This implementation uses the SciPy library’s solve_ivp function, which provides robust numerical integration using adaptive Runge-Kutta methods. The default RK45 method (Dormand-Prince 4(5)) offers a good balance between accuracy and computational efficiency. For stiff parameter regimes, alternative methods like Radau or BDF are available. For more details on the underlying algorithm, see the SciPy integrate documentation.

The function returns time-series data for all three state variables, enabling visualization and analysis of chaotic dynamics directly within Excel.

This example function is provided as-is without any representation of accuracy.

Excel Usage

=LORENZ(x_initial, y_initial, z_initial, sigma, rho, beta, t_start, t_end, timesteps, solve_ivp_method)
  • x_initial (float, required): Initial X value for the system state.
  • y_initial (float, required): Initial Y value for the system state.
  • z_initial (float, required): Initial Z value for the system state.
  • sigma (float, required): Prandtl number (rate of fluid convection).
  • rho (float, required): Rayleigh number (temperature difference driving convection).
  • beta (float, required): Geometric factor (aspect ratio of convection cells).
  • t_start (float, required): Start time of integration.
  • t_end (float, required): End time of integration.
  • timesteps (int, optional, default: 10): Number of timesteps to evaluate.
  • solve_ivp_method (str, optional, default: “RK45”): Integration method to use.

Returns (list[list]): 2D list with header [t, X, Y, Z], or error message string.

Example 1: Classic Lorenz with required arguments only

Inputs:

x_initial y_initial z_initial sigma rho beta t_start t_end
1 1 1 10 28 2.6666666666666665 0 1

Excel formula:

=LORENZ(1, 1, 1, 10, 28, 2.6666666666666665, 0, 1)

Expected output:

t X Y Z
0 1 1 1
0.111111 2.41153 5.09514 1.20106
0.222222 8.3164 17.1331 6.41432
0.333333 19.3183 23.3774 39.0297
0.444444 8.4803 -7.11685 39.9857
0.555556 -2.97791 -8.36678 28.5779
0.666667 -6.45809 -8.3455 25.2686
0.777778 -8.27414 -9.78053 25.2989
0.888889 -9.58941 -10.2266 27.6946
1 -9.3814 -8.42619 29.3065
Example 2: Classic Lorenz with default RK45 method over 2 seconds

Inputs:

x_initial y_initial z_initial sigma rho beta t_start t_end timesteps
1 1 1 10 28 2.6666666666666665 0 2 10

Excel formula:

=LORENZ(1, 1, 1, 10, 28, 2.6666666666666665, 0, 2, 10)

Expected output:

t X Y Z
0 1 1 1
0.222222 8.3164 17.1331 6.41432
0.444444 8.4803 -7.11685 39.9857
0.666667 -6.45809 -8.3455 25.2686
0.888889 -9.58941 -10.2266 27.6946
1.11111 -7.9846 -6.74366 27.9798
1.33333 -7.66622 -8.67649 24.5397
1.55556 -9.83123 -9.72624 28.7982
1.77778 -7.42036 -6.5087 26.9516
2 -8.16794 -9.51754 24.6739
Example 3: Classic Lorenz with RK23 method over 2 seconds

Inputs:

x_initial y_initial z_initial sigma rho beta t_start t_end timesteps solve_ivp_method
1 1 1 10 28 2.6666666666666665 0 2 10 RK23

Excel formula:

=LORENZ(1, 1, 1, 10, 28, 2.6666666666666665, 0, 2, 10, "RK23")

Expected output:

t X Y Z
0 1 1 1
0.222222 8.30806 17.1114 6.39721
0.444444 8.41259 -7.00184 39.7804
0.666667 -6.40742 -8.32007 25.0478
0.888889 -9.69009 -10.3138 27.8124
1.11111 -7.89666 -6.64393 27.9095
1.33333 -7.69375 -8.75884 24.4832
1.55556 -9.82329 -9.62556 28.9107
1.77778 -7.37348 -6.54972 26.7815
2 -8.29039 -9.63569 24.8495
Example 4: Short time span with different initial conditions

Inputs:

x_initial y_initial z_initial sigma rho beta t_start t_end timesteps
0 1 2 10 28 2.6666667 0 0.1 10

Excel formula:

=LORENZ(0, 1, 2, 10, 28, 2.6666667, 0, 0.1, 10)

Expected output:

t X Y Z
0 0 1 2
0.0111111 0.105129 1.00437 1.9422
0.0222222 0.201265 1.03759 1.88721
0.0333333 0.292251 1.09743 1.835
0.0444444 0.381231 1.18284 1.78564
0.0555556 0.471165 1.29339 1.7393
0.0666667 0.564734 1.42938 1.69624
0.0777778 0.664163 1.59213 1.65687
0.0888889 0.771692 1.78351 1.62176
0.1 0.889676 2.00588 1.59164

Python Code

Show Code 
from scipy.integrate import solve_ivp
import numpy as np

def lorenz(x_initial, y_initial, z_initial, sigma, rho, beta, t_start, t_end, timesteps=10, solve_ivp_method='RK45'):
    """
    Numerically solves the Lorenz system of ordinary differential equations for chaotic dynamics.

    See: https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html

    This example function is provided as-is without any representation of accuracy.

    Args:
        x_initial (float): Initial X value for the system state.
        y_initial (float): Initial Y value for the system state.
        z_initial (float): Initial Z value for the system state.
        sigma (float): Prandtl number (rate of fluid convection).
        rho (float): Rayleigh number (temperature difference driving convection).
        beta (float): Geometric factor (aspect ratio of convection cells).
        t_start (float): Start time of integration.
        t_end (float): End time of integration.
        timesteps (int, optional): Number of timesteps to evaluate. Default is 10.
        solve_ivp_method (str, optional): Integration method to use. Valid options: RK45, RK23, DOP853, Radau, BDF, LSODA. Default is 'RK45'.

    Returns:
        list[list]: 2D list with header [t, X, Y, Z], or error message string.
    """
    # Validate input types and ranges
    try:
        x0 = float(x_initial)
        y0 = float(y_initial)
        z0 = float(z_initial)
        sig = float(sigma)
        rh = float(rho)
        bet = float(beta)
        t0 = float(t_start)
        t1 = float(t_end)
        ntp = int(timesteps)
    except Exception:
        return "Error: Invalid input: All initial values, parameters, and timesteps must be numbers."
    if t1 <= t0:
        return "Error: Invalid input: t_end must be greater than t_start."
    if ntp <= 0:
        return "Error: Invalid input: timesteps must be a positive integer."
    if solve_ivp_method not in ['RK45', 'RK23', 'DOP853', 'Radau', 'BDF', 'LSODA']:
        return "Error: Invalid input: solve_ivp_method must be one of 'RK45', 'RK23', 'DOP853', 'Radau', 'BDF', 'LSODA'."
    # Create time array for evaluation
    t_eval = np.linspace(t0, t1, ntp)

    # Lorenz system definition
    def lorenz_ode(t, xyz):
        x, y, z = xyz
        dxdt = sig * (y - x)
        dydt = x * (rh - z) - y
        dzdt = x * y - bet * z
        return [dxdt, dydt, dzdt]
    # Integrate ODE
    try:
        sol = solve_ivp(
            lorenz_ode,
            [t0, t1],
            [x0, y0, z0],
            method=solve_ivp_method,
            dense_output=False,
            t_eval=t_eval
        )
    except Exception as e:
        return f"Error: Integration error: {e}"
    if not sol.success:
        return f"Error: Integration failed: {sol.message}"
    # Format output: header row then data rows
    result = [['t', 'X', 'Y', 'Z']]
    for i in range(len(sol.t)):
        t_val = float(sol.t[i])
        x_val = float(sol.y[0][i])
        y_val = float(sol.y[1][i])
        z_val = float(sol.y[2][i])
        # Disallow nan/inf using math.isfinite for clarity
        vals = [t_val, x_val, y_val, z_val]
        if not all(v == v and abs(v) != float('inf') for v in vals):
            return "Error: Invalid output: nan or inf detected."
        result.append(vals)
    return result

Online Calculator

Initial X value for the system state.
Initial Y value for the system state.
Initial Z value for the system state.
Prandtl number (rate of fluid convection).
Rayleigh number (temperature difference driving convection).
Geometric factor (aspect ratio of convection cells).
Start time of integration.
End time of integration.
Number of timesteps to evaluate.
Integration method to use.