El
sistema de Lorenz es un sistema dinámico determinista 3-dimensional descrito por el siguiente sistema de ecuaciones diferenciales ordinarias (EDOs) no-lineales [1]:
dx/dt = σ(y - x)
dy/dt = x(ρ - z) - y
dz/dt = xy - βz
donde x(t), y(t) y z(t) representan las variables de estado del sistema, mientras que σ, ρ y β son parámetros reales positivos. Dado un estado inicial p0 = (x0, y0, z0) ∈ ℝ3, las ecuaciones anteriores determinan de forma unívoca una trayectoria p(t) en el espacio de fases para t ∈ [0, +∞).
En particular, los valores σ = 10, ρ = 28 y β = 8/3 conducen a un régimen caótico en el que las trayectorias permanecen confinadas en una región acotada del espacio de fases conocida como atractor de Lorenz, cuya forma característica recuerda a una mariposa (véase la animación superior). Aunque el sistema es completamente determinista, una perturbación arbitrariamente pequeña de las condiciones iniciales puede dar lugar a trayectorias muy diferentes con el paso del tiempo.
En este artículo emplearemos técnicas modernas de programación en C++ con el fin de proporcionar una aplicación de consola capaz de representar distintas trayectorias del sistema de Lorenz en una ventana de perspectiva tridimensional, a partir de los parámetros y las condiciones iniciales especificadas por el usuario. Para ello utilizaremos el compilador GCC 16.1 con soporte para C++26 (-std=c++26), junto con las bibliotecas Boost.Numeric.Odeint [2] y Dlib [3].
La biblioteca Boost.Numeric.Odeint proporciona un marco flexible y de alto rendimiento para la resolución numérica de problemas de valores iniciales asociados a sistemas de ecuaciones diferenciales ordinarias (EDOs).
El kit Dlib, por su parte, contiene múltiples algoritmos y herramientas de aprendizaje automático en C++ con aplicaciones en robótica, dispositivos integrados, telefonía móvil y entornos de alto rendimiento, entre otras áreas de interés. En nuestro caso utilizaremos sus funcionalidades gráficas para obtener una representación tridimensional de las trayectorias del sistema de Lorenz.
El lector puede encontrar en
una serie de artículos anterior
una guía detallada para la instalación de un entorno de desarrollo para C++ en MS Windows (64 bits) basado en MSYS2, el compilador GCC y CMake. La instalación de Dlib mediante el gestor de paquetes pacman, incluido en el entorno MSYS2, es inmediata. Basta con abrir la consola MSYS2 UCRT64 y ejecutar el siguiente comando:
pacman -S mingw-w64-ucrt-x86_64-dlib
2. Biblioteca terminal.hpp
A lo largo del artículo utilizaremos una biblioteca sencilla de adquisición de datos por terminal para C++26,
terminal.hpp, analizada en un
post anterior. Su objetivo es simplificar la lectura segura de datos desde la entrada estándar, gestionando automáticamente los posibles errores de formato y solicitando de nuevo el dato de producirse una entrada inválida. La biblioteca proporciona dos funciones genéricas particularmente convenientes para aplicaciones de consola:
- terminal::prompt: Imprime un mensaje en la terminal y aguarda a que el usuario introduzca un valor. De ser válido, el resultado se almacena en una variable capturada por referencia. En caso contrario, se vuelve a solicitar el dato. Opcionalmente, puede especificarse un predicado adicional que el valor introducido deba satisfacer.
- terminal::prompt_init: Función auxiliar similar a la anterior que inicializa internamente un objeto del tipo solicitado antes de proceder a su lectura segura desde la terminal. En lugar de almacenar el resultado en una variable externa, la función devuelve directamente el valor leído una vez validado, lo que permite inicializar objetos constantes de forma directa.
A modo de ejemplo:
// lectura de un double positivo previamente inicializado:
auto num = 0.0;
terminal::prompt(
"Enter a positive number: ",
num,
[](double a){ return a > 0.0; }
);
// inicialización y lectura de un string no-vacío:
auto const str = terminal::prompt_init<std::string>(
"Enter a non-empty string: ",
[](std::string_view w){ return not w.empty(); }
);
3. Aplicación de representación tridimensional
Comenzaremos la implementación incluyendo las cabeceras, tanto estándar como no estándar, en las que se apoyará nuestra aplicación:
#include <algorithm>
#include <cmath>
#include <contracts>
#include <cstdio>
#include <cstdlib>
#include <ctime>
#include <print>
#include <ranges>
#include <string>
#include <string_view>
#include <utility>
#include <vector>
#include <boost/algorithm/string.hpp>
#include <boost/numeric/odeint.hpp>
namespace ode = boost::numeric::odeint;
#include <dlib/gui_widgets.h>
#include <dlib/image_transforms.h>
#include "terminal.hpp"
Como uno de los elementos centrales de nuestra implementación, definiremos un tipo R3_vec para representar el vector de estado del sistema. Cada instancia almacenará las tres coordenadas del espacio euclídeo 3-dimensional x, y, z como datos miembro públicos de tipo double. Con el fin de que R3_vec pueda participar en las operaciones algebraicas requeridas por Boost.Numeric.Odeint, implementaremos la adición de vectores y el producto por escalares reales mediante sobrecarga de operadores.
Las instancias de R3_vec se emplearán en combinación con algoritmos de integración numérica de EDOs con corrección adaptativa de la longitud de paso, por lo que de acuerdo con la documentación de Boost.Numeric.Odeint deberemos:
- Proporcionar una función abs() que, dado un vector, retorne otro vector cuyas coordenadas sean los valores absolutos de las correspondientes coordenadas originales.
- Sobrecargar el operador operator/ para realizar el cociente componente a componente (entry-wise) de dos vectores.
- Especializar una operación que proporcione la norma infinito de un vector, ‖v‖∞ = máx{∣x∣,∣y∣,∣z∣}, definida como el mayor de los valores absolutos de sus coordenadas. Esta norma será empleada por el mecanismo de control adaptativo del paso de integración para estimar el error numérico cometido en cada iteración [4].
La implementación de la clase R3_vec y sus funciones asociadas resulta directa [5]:
struct R3_vec {
double x{};
double y{};
double z{};
constexpr R3_vec() noexcept = default;
// conversión implícita requerida por Boost.Numeric.Odeint (vector_space_algebra):
// un escalar 'a' se interpreta como el vector '(a,a,a)'
constexpr R3_vec(double a) noexcept
: x{a}, y{a}, z{a}
{ }
constexpr R3_vec(double x_, double y_, double z_) noexcept
: x{x_}, y{y_}, z{z_}
{ }
constexpr auto& operator+=(R3_vec const& v) noexcept
{
x += v.x;
y += v.y;
z += v.z;
return *this;
}
constexpr auto& operator-=(R3_vec const& v) noexcept
{
x -= v.x;
y -= v.y;
z -= v.z;
return *this;
}
constexpr auto& operator*=(double a) noexcept
{
x *= a;
y *= a;
z *= a;
return *this;
}
constexpr auto& operator/=(double a) noexcept
pre(a != 0.0)
{
x /= a;
y /= a;
z /= a;
return *this;
}
};
[[nodiscard]]
constexpr auto operator+(R3_vec lhs, R3_vec const& rhs) noexcept -> R3_vec
{
lhs += rhs;
return lhs;
}
[[nodiscard]]
constexpr auto operator-(R3_vec lhs, R3_vec const& rhs) noexcept -> R3_vec
{
lhs -= rhs;
return lhs;
}
[[nodiscard]]
constexpr auto operator*(double a, R3_vec v) noexcept -> R3_vec
{
v *= a;
return v;
}
[[nodiscard]]
constexpr auto operator*(R3_vec v, double a) noexcept -> R3_vec
{
return a*v;
}
[[nodiscard]]
constexpr auto operator/(R3_vec v, double a) noexcept -> R3_vec
pre(a != 0.0)
{
v /= a;
return v;
}
[[nodiscard]]
constexpr auto abs(R3_vec const& v) noexcept -> R3_vec
{
using std::abs;
return {abs(v.x), abs(v.y), abs(v.z)};
}
[[nodiscard]]
constexpr auto operator/(R3_vec const& a, R3_vec const& b) noexcept -> R3_vec
pre(b.x != 0.0 and b.y != 0.0 and b.z != 0.0)
{
return {a.x / b.x, a.y / b.y, a.z / b.z};
}
namespace boost::numeric::odeint {
template<>
struct vector_space_norm_inf<R3_vec> {
using result_type = double;
[[nodiscard]]
constexpr auto operator()(R3_vec const& v) const noexcept -> result_type
{
using std::abs;
using std::max;
return max({abs(v.x), abs(v.y), abs(v.z)});
}
};
}
Conviene señalar que la biblioteca Dlib proporciona la plantilla de clase dlib::vector<> para representar vectores 2- y 3-dimensionales con las operaciones habituales del álgebra vectorial. Sin embargo, dicha clase no cumple directamente los requisitos que Boost.Numeric.Odeint impone sobre los tipos utilizados como vectores de estado en algoritmos de integración adaptativa, lo que justifica la definición de nuestro propio tipo R3_vec. Ello nos permite, asimismo, personalizar su interfaz, incorporar contratos a determinadas operaciones y utilizar de forma natural características del lenguaje como los structured bindings sobre sus coordenadas.
Una vez definido el tipo R3_vec que actuará como vector de estado, procederemos a modelar las ecuaciones diferenciales del sistema de Lorenz. Para ello, comenzaremos agrupando los parámetros σ, ρ y β en un agregado de datos Lorenz_parameters, cuyos valores por defecto coincidirán con los utilizados habitualmente en la literatura para obtener el atractor de Lorenz:
struct Lorenz_parameters {
double sigma = 10.0;
double rho = 28.0;
double beta = 8.0/3.0;
};
A continuación definiremos la estructura Lorenz_system encargada de representar el sistema de ecuaciones diferenciales. De acuerdo con las convenciones de Boost.Numeric.Odeint, éste se implementa como un function object (functor) cuyo operador de llamada recibe el estado actual del sistema, un vector donde almacenar sus derivadas temporales y el instante de tiempo correspondiente:
struct Lorenz_system {
Lorenz_parameters const params;
// derivadas temporales de las coordenadas x, y, z:
void operator()(
R3_vec const& state,
R3_vec& drdt,
[[maybe_unused]] double t
) const noexcept
{
auto const& [sigma, rho, beta] = params;
auto const& [x, y, z] = state;
drdt = {sigma*(y - x), x*(rho - z) - y, x*y - beta*z};
}
};
El parámetro state contiene las coordenadas actuales (x, y, z), mientras que drdt (capturado por referencia) recibe los valores (dx/dt, dy/dt, dz/dt) calculados a partir de las ecuaciones del sistema.
Observemos que el operador anterior no realiza aún ninguna integración numérica; se limita a describir la dinámica del sistema proporcionando las derivadas temporales.
En este artículo utilizaremos el método de Runge–Kutta–Dormand–Prince de orden 5(4) (RKDP), implementado en Boost.Numeric.Odeint mediante la plantilla runge_kutta_dopri5 [6]. Este método proporciona simultáneamente una aproximación de la solución y una estimación del error local cometido en cada paso, lo que permite ajustar automáticamente el tamaño del paso temporal Δt para satisfacer unas tolerancias de error prefijadas [7]. En nuestro caso, fijaremos las tolerancias de error absoluto y relativo en 1.0e-7 y estableceremos un tamaño máximo de paso de 5.0e-3.
La siguiente función lorenz_integration() construye un rango adaptativo mediante make_adaptive_range(), capaz de generar sucesivamente los estados del sistema entre t = 0 y el tiempo final de integración integration_time especificado por el usuario:
using Rkdp_stepper = ode::runge_kutta_dopri5<
R3_vec,
double,
R3_vec,
double,
ode::vector_space_algebra
>;
using Adaptive_iter = ode::adaptive_iterator<
ode::controlled_runge_kutta<Rkdp_stepper>,
Lorenz_system,
R3_vec
>;
[[nodiscard]]
auto lorenz_integration(
Lorenz_parameters
const& params,
R3_vec
& state,
// estado I/O (requerido por Boost.Numeric.Odeint)
double integration_time,
double max_step_sz
) -> std::pair<Adaptive_iter, Adaptive_iter>
{
// inicializamos un 'stepper' con control adaptativo de paso
// que emplee RKDP:
constexpr auto abs_tol =
1.e-7;
constexpr auto rel_tol =
1.e-7;
auto const stepper = ode::make_controlled<Rkdp_stepper>(
abs_tol,
rel_tol,
max_step_sz
);
// obtenemos los estados que actuarán como extremos de los segmentos
// [e_0, e_1], entre t = 0 e integration_time:
auto const initial_step_sz = max_step_sz /
10.0;
return ode::make_adaptive_range(
stepper,
Lorenz_system{params},
state,
// estado modificado in situ por el integrador
0.0,
integration_time,
initial_step_sz
);
}
El resultado de la función es un par de iteradores Adaptive_iter que delimitan el rango adaptativo de estados generado por el integrador. Estos iteradores modelan el concepto single-pass iterator [8], por lo que no permiten realizar múltiples recorridos sobre la secuencia de estados.
A continuación, implementaremos la función get_colored_polyline(), que actuará como puente entre la integración numérica del sistema y su representación gráfica. A partir de los parámetros del sistema (σ, ρ y β), un estado inicial y un tiempo total de integración, construiremos una curva poligonal (polilínea) tridimensional formada por segmentos rectos [e_0, e_1], cuyos vértices corresponderán a los sucesivos estados producidos por el integrador adaptativo. Dichos segmentos serán almacenados en un contenedor std::vector cuyas entradas, de tipo dlib::perspective_window::overlay_line, podrán ser embebidas posteriormente en un entorno de representación 3-dimensional:
using Trajectory = std::vector<dlib::perspective_window::overlay_line>;
[[nodiscard]]
auto get_colored_polyline(
Lorenz_parameters
const& params,
R3_vec init,
double integration_time
) -> Trajectory
{
auto polyline = Trajectory{};
constexpr auto max_step_sz =
5.e-3;
polyline.reserve(
static_cast<std::size_t>(integration_time/max_step_sz));
// preservamos el estado inicial:
auto e_0 = init;
// obtenemos el rango adaptativo:
auto const [first, last] = lorenz_integration(
params,
init,
integration_time,
max_step_sz
);
// construimos la curva poligonal a partir del rango adaptativo:
std::for_each(
std::next(first),
last,
[&polyline, &e_0](R3_vec
const& e_1) {
using Vec = dlib::vector<
double>;
polyline.emplace_back(
Vec{e_0.x, e_0.y, e_0.z},
Vec{e_1.x, e_1.y, e_1.z}
);
e_0 = e_1;
}
);
// coloreamos los segmentos de la trayectoria:
for (
auto const max_value =
static_cast<
double>(polyline.size()) -
1.0;
auto&& [idx, segment] : polyline | std::views::enumerate
) {
segment.color = dlib::colormap_jet(idx,
0.0, max_value);
}
return polyline;
}
La llamada a std::next(first) en el bucle std::for_each de la función anterior evita la generación de un segmento degenerado de extremos coincidentes e_0 y longitud nula.
La función anterior asigna un color a cada segmento empleando una paleta de colores. La animación mostrada en la introducción y la figura inferior representan el atractor de Lorenz correspondiente a los parámetros por defecto utilizando el mapa de color dlib::colormap_heat. No obstante, para distinguir con mayor claridad las distintas regiones de la trayectoria, puede resultar preferible emplear dlib::colormap_jet u otro mapa de color alternativo, según las preferencias del programador.
La siguiente función plot_system() creará una ventana de perspectiva tridimensional mediante la clase dlib::perspective_window, superponiendo sobre ella los distintos elementos gráficos que conforman la representación del sistema:
- La polilínea coloreada tridimensional generada por get_colored_polyline().
- Los ejes cartesianos como segmentos coloreados: el eje x en azul, el eje y en verde y el eje z en blanco.
Entre otras posibilidades, podremos rotar la representación para contemplarla desde distintas perspectivas, así como ampliar regiones concretas de la trayectoria con el fin de examinar su estructura con mayor detalle:
void plot_system(
R3_vec const& init,
Trajectory const& polyline
){
auto const& [x_0, y_0, z_0] = init;
std::print("Plotting from ({:.2f}, {:.2f}, {:.2f})...", x_0, y_0, z_0);
// abrimos una ventana de perspectiva 3D:
auto win = dlib::perspective_window{};
win.set_title("Lorenz system - x axis (blue), y axis (green), z axis (white)");
win.set_size(640, 640);
// superponemos los ejes cartesianos:
auto const origin = dlib::vector<double>{}; // (0,0,0)
win.add_overlay(origin, {10.0, 0.0, 0.0}, dlib::rgb_pixel{0, 0, 255}); // eje x
win.add_overlay(origin, {0.0, 10.0, 0.0}, dlib::rgb_pixel{0, 255, 0}); // eje y
win.add_overlay(origin, {0.0, 0.0, 10.0}, dlib::rgb_pixel{255, 255, 255}); // eje z
// superponemos la curva poligonal:
win.add_overlay(polyline);
// aguardamos a que el usuario cierre la ventana:
win.wait_until_closed();
}
Llegados a este punto en nuestra implementación, proporcionaremos al usuario la opción de emplear los valores clásicos σ = 10, ρ = 28 y β = 8/3 —para los que el sistema exhibe comportamiento caótico y da lugar al conocido atractor de Lorenz—, o bien especificar valores alternativos para dichos parámetros. Asimismo, el programa permitirá seleccionar un estado inicial p0 = (x0, y0, z0) aleatorio en las proximidades del origen o introducir manualmente dichas coordenadas. Finalmente, se solicitará al usuario un tiempo total de integración:
Conviene recordar que el origen constituye un punto de equilibrio del sistema, cuya estabilidad depende del valor del parámetro ρ. En particular, el origen es asintóticamente estable cuando ρ < 1, mientras que para ρ > 1 pierde dicha estabilidad y aparecen otros dos puntos de equilibrio distintos del origen.
[[nodiscard]]
auto affirmative_answer_to(std::string_view question) -> bool
{
auto is_yes = false;
[[maybe_unused]] auto _ = terminal::prompt_init<std::string>(
question,
[&is_yes](std::string_view ans) {
is_yes = boost::iequals(ans, "y")
or boost::iequals(ans, "yes");
return is_yes or boost::iequals(ans, "n")
or boost::iequals(ans, "no");
}
);
return is_yes;
}
auto main() -> int
try {
std::print("{:_^70}\n", "Plotting the Lorenz System");
auto greater_than_zero = [](double v){ return v > 0.0; };
auto const params = /* IILE: parámetros σ, ρ, β */ [&] -> Lorenz_parameters {
auto res = Lorenz_parameters{};
if (
not affirmative_answer_to(
"Use default Lorenz parameters (σ=10, ρ=28, β=8/3)? [y/n]: ")
) {
terminal::prompt(" · σ [> 0.0]: ", res.sigma, greater_than_zero);
terminal::prompt(" · ρ [> 0.0]: ", res.rho, greater_than_zero);
terminal::prompt(" · β [> 0.0]: ", res.beta, greater_than_zero);
}
return res;
}();
auto const p_0 = /* IILE: estado inicial */ [&]{
auto res = R3_vec{};
if (
affirmative_answer_to(
"Use random initial position near (0,0,0)? [y/n]: ")
) {
auto rnd = [rng = dlib::rand{std::time(nullptr)}] mutable {
return rng.get_double_in_range(-5.0, 5.0);
};
res = {rnd(), rnd(), rnd()};
}
else {
terminal::prompt(" · x_0 ∈ (-∞,+∞): ", res.x);
terminal::prompt(" · y_0 ∈ (-∞,+∞): ", res.y);
terminal::prompt(" · z_0 ∈ (-∞,+∞): ", res.z);
}
return res;
}();
// tiempo de integración:
auto const integration_time = terminal::prompt_init<double>(
"Integration time [> 0.0]: ",
greater_than_zero
);
// obtenemos la curva poligonal resultante:
auto const polyline = get_colored_polyline(params, p_0, integration_time);
// mostramos la trayectoria en una ventana de representación 3D:
plot_system(p_0, polyline);
return EXIT_SUCCESS;
}
catch (std::exception const& e) {
std::print(stderr, "exception: {}", e.what());
return EXIT_FAILURE;
}
La función auxiliar affirmative_answer_to permite plantear al usuario una pregunta dicotómica de tipo sí/no, retornando true en caso de que su respuesta sea afirmativa y false de ser negativa.
A modo de ejemplo, la siguiente ejecución del programa genera una órbita estable para σ = 10, ρ = 14 y β = 8/3:
4. Fichero de configuración CMakeLists.txt
El siguiente fichero CMakeLists.txt permite configurar la construcción del proyecto mediante CMake, suponiendo que el código desarrollado en la sección anterior se encuentre en un archivo fuente denominado lorenz.cpp. Tanto la cabecera terminal.hpp como el propio fichero de configuración CMakeLists.txt deberán ubicarse en el mismo directorio que lorenz.cpp. Se asume implícitamente el empleo de GCC 16 o superior:
cmake_minimum_required(VERSION 4.2)
project(
lorenz_system
VERSION 0.1.0
DESCRIPTION "Plotting the Lorenz System - dgvergel.blogspot.com"
LANGUAGES CXX
)
add_executable(${PROJECT_NAME})
target_sources(${PROJECT_NAME} PRIVATE
lorenz.cpp
terminal.hpp
)
target_compile_features(${PROJECT_NAME}
PRIVATE
cxx_std_26
)
set_target_properties(${PROJECT_NAME} PROPERTIES
CXX_EXTENSIONS OFF
RUNTIME_OUTPUT_DIRECTORY "${CMAKE_BINARY_DIR}/bin"
)
find_package(Boost REQUIRED)
find_package(dlib REQUIRED)
target_link_libraries(${PROJECT_NAME}
PRIVATE
Boost::boost
dlib::dlib
-lstdc++exp
)
target_compile_options(${PROJECT_NAME} PRIVATE
-fcontracts
-Wall
-Wextra
-Wpedantic
-Werror
)
message(STATUS "Generator: ${CMAKE_GENERATOR}")
message(STATUS "Compiler: ${CMAKE_CXX_COMPILER_ID} ${CMAKE_CXX_COMPILER_VERSION}")
message(STATUS "Boost version: ${Boost_VERSION}")
Referencias bibliográficas
- Wikipedia - Lorenz System - https://en.wikipedia.org/wiki/Lorenz_system
- Boost.Numeric.Odeint - https://www.boost.org/doc/libs/latest/libs/numeric/odeint/doc/html/index.html
- Dlib C++ library - http://dlib.net
- Boost.Numeric.Odeint - Controlled steppers - https://www.boost.org/doc/libs/latest/libs/numeric/odeint/doc/html/boost_numeric_odeint/odeint_in_detail/steppers.html#boost_numeric_odeint.odeint_in_detail.steppers.controlled_steppers
- Boost.Numeric.Odeint - State types, algebras, and operations - Point type - https://www.boost.org/doc/libs/latest/libs/numeric/odeint/doc/html/boost_numeric_odeint/odeint_in_detail/state_types__algebras_and_operations.html#boost_numeric_odeint.odeint_in_detail.state_types__algebras_and_operations.algebras_and_operations.point_type
- Wikipedia - Dormand–Prince method - https://en.wikipedia.org/wiki/Dormand%E2%80%93Prince_method
- Boost.Numeric.Odeint - Class template runge_kutta_dopri5 - https://www.boost.org/doc/libs/latest/libs/numeric/odeint/doc/html/doxygen/odeint_reference/classboost_1_1numeric_1_1odeint_1_1runge__kutta__dopri5.html#id51036
- Boost.Numeric.Odeint - Class template adaptive_iterator - https://www.boost.org/doc/libs/latest/libs/numeric/odeint/doc/html/doxygen/odeint_reference/classboost_1_1numeric_1_1odeint_1_1adaptive__iterator.html
Comentarios
Publicar un comentario