Lorenz System

Overview

This document focuses on demonstrating the Lorenz system, originally introduced by Edward N. Lorenz in his seminal (1963) paper. The Lorenz system comprises three coupled, nonlinear ordinary differential equations that model atmospheric convection, effectively illustrating the chaotic nature of weather patterns.

The dynamics of the model are represented by the following set of first-order, nonlinear differential equations:

\[ \begin{aligned} \frac{dx}{dt} &= \sigma(y - x), \\ \frac{dy}{dt} &= x(\rho - z) - y \\ \frac{dz}{dt} &= xy - \beta z \end{aligned} \] In these equations:

• \(x\) represents the rate of convection;
• \(y\) denotes the horizontal temperature variation;
• \(z\) indicates the vertical temperature variation;
• \(\sigma\), \(\rho\), and \(\beta\) are system parameters corresponding to the Prandtl number, Rayleigh number, and specific physical dimensions of the fluid layer.

To learn more about the Lorenz system, see Lorenz (2008).

Setting up the environment

library(checkmate, quietly = TRUE)
library(deSolve, quietly = TRUE)
library(dplyr, quietly = TRUE)
library(ggplot2, quietly = TRUE)
library(gg3D, quietly = TRUE) # remotes::install_github("AckerDWM/gg3D")
library(latex2exp, quietly = TRUE)
library(magrittr, quietly = TRUE)

Numerical solution of the equations

lorenz_system <- function(
    x = 1, 
    y = 1, 
    z = 1, 
    sigma = 10, 
    rho = 28, 
    beta = 8 / 3, 
    from = 0, 
    to = 100,
    by = 0.01
  ) {
  checkmate::assert_number(x)
  checkmate::assert_number(y)
  checkmate::assert_number(z)
  checkmate::assert_number(sigma)
  checkmate::assert_number(rho)
  checkmate::assert_number(beta)
  checkmate::assert_number(from, lower = 0)
  checkmate::assert_number(to, lower = from)
  checkmate::assert_number(by, lower = 0)
  
  fun <- function (t, y, parms) {
    list2env(as.list(y), envir = environment())
    list2env(as.list(parms), envir = environment())
    
    list(
      c(
        dx = sigma * (y - x),
        dy = x * (rho - z) - y,
        dz = (x * y) - (beta * z)
      )
    )
  }
  
  initial_values <- c(x = x, y = y, z = z)
  parameters <- list(sigma = sigma, rho = rho, beta = beta)
  time <- seq(from = from, to = to, by = by)
  
  data <- 
    deSolve::ode(
      y = initial_values,
      times = time, 
      func = fun,
      parms = parameters
    ) |>
    dplyr::as_tibble() |>
    dplyr::mutate(dplyr::across(dplyr::everything(), ~ as.numeric(.x)))
  
  list(
    data = data,
    initial_values = as.list(initial_values),
    parameters = as.list(parameters)
  ) |>
  invisible()
}

lorenz_system() |> magrittr::extract2("data")
#> # A tibble: 10,001 × 4
#>    time     x     y     z
#>   <dbl> <dbl> <dbl> <dbl>
#> 1  0     1     1    1    
#> 2  0.01  1.01  1.26 0.985
#> 3  0.02  1.05  1.52 0.973
#> 4  0.03  1.11  1.80 0.965
#> 5  0.04  1.19  2.09 0.962
#> 6  0.05  1.29  2.40 0.964
#> # ℹ 9,995 more rows

Plotting system dynamics

plot_system_dynamics <- function(
    x = 1, 
    y = 1, 
    z = 1, 
    sigma = 10, 
    rho = 28, 
    beta = 8 / 3, 
    from = 0, 
    to = 100,
    by = 0.01
  ) {
  checkmate::assert_number(x)
  checkmate::assert_number(y)
  checkmate::assert_number(z)
  checkmate::assert_number(sigma)
  checkmate::assert_number(rho)
  checkmate::assert_number(beta)
  checkmate::assert_number(from, lower = 0)
  checkmate::assert_number(to, lower = from)
  checkmate::assert_number(by, lower = 0)
  
  lorenz_system(x, y, z, sigma, rho, beta, from, to, by) |> 
    list2env(envir = environment())
  
  plot <- 
    data |>
    ggplot2::ggplot(ggplot2::aes(x = time)) +
    ggplot2::geom_line(
      ggplot2::aes(y = y, color = "Horizontal temperature variation (Y)"),
      linewidth = 0.5,
      alpha = 0.75
    ) +
    ggplot2::geom_line(
      ggplot2::aes(y = z, color = "Vertical temperature variation (Z)"),
      linewidth = 0.5,
      alpha = 0.75
    ) +
    ggplot2::geom_line(
      ggplot2::aes(y = x, color = "Rate of convection (X)"),
      linewidth = 0.5,
      alpha = 0.75
    ) +
    ggplot2::labs(
      title = "Lorenz System Dynamics",
      subtitle = latex2exp::TeX(
        paste0(
          "$X_0$ = ", x, " | ",
          "$Y_0$ = ", y, " | ",
          "$Z_0$ = ", z, " | ",
          "$\\sigma$ = ", round(sigma, 2), " | ",
          "$\\rho$ = ", round(rho, 2), " | ",
          "$\\beta$ = ", round(beta, 2)
        ),
      ),
      x = "Time", 
      y = "Values",
      color = ggplot2::element_blank()
    ) +
    ggplot2::scale_color_manual(
      breaks = c(
        "Rate of convection (X)", 
        "Horizontal temperature variation (Y)", 
        "Vertical temperature variation (Z)"
      ),
      values = c("blue", "red", "black")
    )
  
  print(plot)
  
  invisible()
}

plot_system_dynamics()

Phase space visualization

plot_phase_space <- function(
    x = 1, 
    y = 1, 
    z = 1, 
    sigma = 10, 
    rho = 28, 
    beta = 8 / 3, 
    from = 0, 
    to = 100,
    by = 0.001,
    theta = 180,
    phi = 0
  ) {
  checkmate::assert_number(x)
  checkmate::assert_number(y)
  checkmate::assert_number(z)
  checkmate::assert_number(sigma)
  checkmate::assert_number(rho)
  checkmate::assert_number(beta)
  checkmate::assert_number(from, lower = 0)
  checkmate::assert_number(to, lower = from)
  checkmate::assert_number(by, lower = 0)
  
  lorenz_system(x, y, z, sigma, rho, beta, from, to, by) |> 
    list2env(envir = environment())

  plot <-
    data |>
    ggplot2::ggplot(ggplot2::aes(x = x, y = y, z = z, colour = time)) +
    gg3D::stat_3D(theta = theta, phi = phi, geom = "path") +
    theme_void() +
    theme(legend.position = "none")
  
  print(plot)
  
  invisible()
}

plot_phase_space()

References

Lorenz, E. N. (1963). Deterministic nonperiodic flow. Journal of the Atmospheric Sciences, 20(2), 130–141. https://doi.org/10.1175/1520-0469(1963)020<0130:DNF>2.0.CO;2

Lorenz, E. N. (2008). The essence of chaos. UCL Press.