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.1 2.13316 4.47121 1.11428
0.2 6.53431 13.7007 4.19251
0.3 16.7 27.2348 26.183
0.4 15.4555 1.27525 46.9318
0.5 1.23198 -8.91266 32.5575
0.6 -4.84354 -8.08004 26.7579
0.7 -7.04699 -8.71664 24.9578
0.8 -8.60266 -10.0354 25.6474
0.9 -9.65321 -10.1321 27.9586
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.2 6.53431 13.7007 4.19251
0.4 15.4555 1.27525 46.9318
0.6 -4.84354 -8.08004 26.7579
0.8 -8.60266 -10.0354 25.6474
1 -9.38131 -8.42656 29.3074
1.2 -7.21285 -6.80228 26.042
1.4 -8.47725 -9.80609 25.1204
1.6 -9.6297 -8.83409 29.3803
1.8 -7.2442 -6.57594 26.4335
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.2 6.53363 13.7146 4.16899
0.4 15.3295 1.16572 46.6518
0.6 -4.76281 -7.98779 26.5741
0.8 -8.67439 -10.1645 25.57
1 -9.37985 -8.30481 29.4198
1.2 -7.15133 -6.791 25.9084
1.4 -8.53679 -9.89862 25.1561
1.6 -9.57153 -8.69011 29.4074
1.8 -7.22035 -6.64725 26.2765
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.01 0.0950919 1.00259 1.94785
0.02 0.182556 1.02877 1.89799
0.03 0.265305 1.07675 1.85037
0.04 0.345687 1.14565 1.80504
0.05 0.425899 1.23498 1.76208
0.06 0.508028 1.34467 1.72166
0.07 0.593848 1.47532 1.68402
0.08 0.684947 1.62805 1.64949
0.09 0.78298 1.80431 1.61851
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.
    """
    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)

      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'."

      t_eval = np.linspace(t0, t1, ntp + 1)

      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]

      sol = solve_ivp(
        lorenz_ode,
        [t0, t1],
        [x0, y0, z0],
        method=solve_ivp_method,
        dense_output=False,
        t_eval=t_eval
      )
      if not sol.success:
        return f"Error: Integration failed: {sol.message}"

      result = [['t', 'X', 'Y', 'Z']]
      for i in range(len(sol.t)):
        vals = [float(sol.t[i]), float(sol.y[0][i]), float(sol.y[1][i]), float(sol.y[2][i])]
        if not all(np.isfinite(v) for v in vals):
          return "Error: Invalid output: nan or inf detected."
        result.append(vals)
      return result
    except Exception as e:
      return f"Error: {str(e)}"

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.