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.

Examples

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.1111	2.412	5.095	1.201
0.2222	8.316	17.13	6.414
0.3333	19.32	23.38	39.03
0.4444	8.48	-7.117	39.99
0.5556	-2.978	-8.367	28.58
0.6667	-6.458	-8.345	25.27
0.7778	-8.274	-9.781	25.3
0.8889	-9.589	-10.23	27.69
1	-9.381	-8.426	29.31

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.2222	8.316	17.13	6.414
0.4444	8.48	-7.117	39.99
0.6667	-6.458	-8.345	25.27
0.8889	-9.589	-10.23	27.69
1.111	-7.985	-6.744	27.98
1.333	-7.666	-8.676	24.54
1.556	-9.831	-9.726	28.8
1.778	-7.42	-6.509	26.95
2	-8.168	-9.518	24.67

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.2222	8.308	17.11	6.397
0.4444	8.413	-7.002	39.78
0.6667	-6.407	-8.32	25.05
0.8889	-9.69	-10.31	27.81
1.111	-7.897	-6.644	27.91
1.333	-7.694	-8.759	24.48
1.556	-9.823	-9.626	28.91
1.778	-7.373	-6.55	26.78
2	-8.29	-9.636	24.85

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.01111	0.1051	1.004	1.942
0.02222	0.2013	1.038	1.887
0.03333	0.2923	1.097	1.835
0.04444	0.3812	1.183	1.786
0.05556	0.4712	1.293	1.739
0.06667	0.5647	1.429	1.696
0.07778	0.6642	1.592	1.657
0.08889	0.7717	1.784	1.622
0.1	0.8897	2.006	1.592

Python 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 "Invalid input: All initial values, parameters, and timesteps must be numbers."
    if t1 <= t0:
        return "Invalid input: t_end must be greater than t_start."
    if ntp <= 0:
        return "Invalid input: timesteps must be a positive integer."
    if solve_ivp_method not in ['RK45', 'RK23', 'DOP853', 'Radau', 'BDF', 'LSODA']:
        return "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"Integration error: {e}"
    if not sol.success:
        return f"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 "Invalid output: nan or inf detected."
        result.append(vals)
    return result

Overview

Excel Usage

Examples

Python Code

Online Calculator