Introducción
Última actualización del artículo: 9 de noviembre de 2021.
El
sistema de Lorenz es un sistema dinámico determinista 3-dimensional descrito por las siguientes ecuaciones diferenciales ordinarias (EDOs) no-lineales [1]:
dx/dt = σ(y - x)
dy/dt = x(ρ - z) - y
dz/dt = xy - βz
donde σ, ρ y β son parámetros reales positivos. Dada una posición inicial p0 = (x0, y0, z0) ∈ ℝ3, el sistema seguirá una órbita en el espacio de fases parametrizada por el tiempo t ∈ [0, +∞).
En particular, los valores σ = 10, ρ = 28 y β = 8/3 conducen a un comportamiento caótico del sistema en el que la mayoría de las trayectorias tienden a un estructura denominada atractor de Lorenz, de aspecto similar a una mariposa (véase la animación superior). Una mínima modificación de las condiciones iniciales producirá órbitas irregulares completamente distintas.
En este artículo emplearemos técnicas de programación modernas propias de C++17/20/23 con el fin de proporcionar una aplicación de consola que represente en una ventana de perspectiva tridimensional diversas trayectorias en el espacio fásico, según los parámetros y las condiciones iniciales escogidas por el usuario. Para ello, utilizaremos el compilador GCC 11.2 bajo el flag -std=c++23, así como las bibliotecas no-estándar Boost.Numeric.Odeint [1], {fmt} [2], Range-v3 [3] y Dlib [4]. El lector puede encontrar en una serie de artículos anterior una guía de instalación de un entorno de trabajo para proyectos en C++ bajo MS Windows (64-bit) basado en MSYS2, el compilador mingw-w64 y CMake, así como los comandos de instalación de las tres primeras bibliotecas mencionadas a través del gestor de paquetes pacman.
La biblioteca Boost.Numeric.Odeint proporciona una solución flexible y de alto rendimiento para la resolución de problemas de valores iniciales en sistemas de 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. Sus funcionalidades gráficas nos permitirán obtener una representación tridimensional del atractor. Su instalación a través del gestor de paquetes pacman proporcionado por el entorno MSYS2 resulta inmediata. Basta abrir el shell MSYS2 MSYS y ejecutar el comando pacman -S mingw-w64-x86_64-dlib.
Biblioteca terminal.hpp
A lo largo del artículo haremos uso de una biblioteca para C++20 de adquisición de datos por la terminal,
terminal.hpp, analizada en un
post anterior. Ésta incluye dos funciones genéricas que gestionan de forma automática los posibles errores de formato en los datos proporcionados por el usuario:
- terminal::prompt: Imprime un mensaje en la terminal y aguarda a que el usuario introduzca el dato solicitado. Si el input es válido, éste es almacenado en una variable externa capturada por referencia. En caso contrario, vuelve a solicitarse el dato. Opcionalmente, es posible proporcionar una condición que deba ser satisfecha por el input.
- terminal::prompt_init: Una función auxiliar similar a la anterior que inicializa ella misma por defecto una variable del tipo a leer y procede a su lectura segura por la terminal. La función retorna el valor, que puede declararse constante de ser conveniente.
A modo de ejemplo:
// lectura de un double positivo previamente inicializado:
auto d = 0.0;
terminal::prompt("Insert a positive double: ", d, [](auto a){ return a > 0.0; });
// inicialización y lectura de un string no vacío:
auto const s = terminal::prompt_init<std::string>("Insert a non-empty string: ",
[](auto const& w){ return !w.empty(); });
Aplicación de representación tridimensional
Comenzaremos nuestra tarea de codificación incluyendo las cabeceras, tanto estándares como no-estándares, en las que se apoyará nuestra aplicación. En el caso de terminal.hpp, la definición previa de la macro TERMINAL_HPP_USE_FMTLIB garantiza que la emisión de mensajes por la terminal se realice a través de la función fmt::print() de la biblioteca {fmt}lib:
#include <algorithm>
#include <array>
#include <cmath>
#include <ctime>
#include <iterator>
#include <string>
#include <type_traits>
#include <vector>
#include <boost/algorithm/string.hpp>
#include <boost/numeric/odeint.hpp>
#include <boost/operators.hpp>
#include <dlib/gui_widgets.h>
#include <dlib/image_transforms.h>
#include <range/v3/view/enumerate.hpp>
#define TERMINAL_HPP_USE_FMTLIB
#include "terminal.hpp"
Como aspecto clave de nuestro código, definiremos un tipo R3_point cuyos objetos representen vectores de posición en el espacio euclídeo 3-dimensional. Las coordenadas serán accesibles a través de tres datos miembro públicos de tipo double y nombres x, y, z. Gracias a la sobrecarga de operadores, garantizaremos que dichos objetos satisfagan el álgebra habitual de adición de vectores y producto por escalares reales (+, +=, *, *=). R3_point derivará de clases proporcionadas en Boost.Operators con el fin de facilitar la definición de tales operaciones.
Las instancias de R3_point se emplearán en conjunción con algoritmos de integración numérica de EDOs con corrección adaptativa de longitud de pasos, por lo que de acuerdo con la documentación de la biblioteca Boost.Numeric.Odeint deberemos proporcionar una función abs() que retorne, para un vector dado, un nuevo vector cuyas entradas se calculen como los valores absolutos de las coordenadas del original. Asimismo, deberemos sobrecargar el operador operator/ con el fin de realizar el cociente, coordenada a coordenada (entry-wise), de dos vectores. Finalmente, habremos de definir una operación que proporcione la norma infinito de un vector, ‖∙‖∞, definida como el valor máximo del conjunto formado por los valores absolutos de sus coordenadas. Dichas operaciones serán utilizadas por el algoritmo de control de paso [5]. La implementación de la clase R3_point y sus funciones asociadas es directa y puede consultarse en las referencias [6, 7]:
struct R3_point : boost::additive1<R3_point,
boost::additive2<R3_point, double,
boost::multiplicative2<R3_point, double>>> {
double x = 0.0,
y = 0.0,
z = 0.0;
R3_point() = default;
R3_point(double a) noexcept : x{a}, y{a}, z{a} { }
R3_point(double a, double b, double c) noexcept : x{a}, y{b}, z{c} { }
auto& operator+=(R3_point const& p) noexcept
{
x += p.x; y += p.y; z += p.z;
return *this;
}
auto& operator-=(R3_point const& p) noexcept
{
x -= p.x; y -= p.y; z -= p.z;
return *this;
}
auto& operator*=(double c) noexcept
{
x *= c; y *= c; z *= c;
return *this;
}
auto& operator/=(double c) noexcept
{
x /= c; y /= c; z /= c;
return *this;
}
};
auto abs(R3_point const& p) noexcept
{
return R3_point{std::abs(p.x), std::abs(p.y), std::abs(p.z)};
}
auto operator/(R3_point const& p1, R3_point const& p2) noexcept
{
return R3_point{p1.x/p2.x, p1.y/p2.y, p1.z/p2.z};
}
namespace boost::numeric::odeint {
template<> struct vector_space_norm_inf<R3_point> {
using result_type = double;
auto operator()(R3_point const& p) const noexcept -> result_type
{
using std::max;
using std::abs;
return max(max(abs(p.x), abs(p.y)), abs(p.z));
}
};
}
Señalemos, en este punto, que la biblioteca Dlib proporciona una plantilla de clase dlib::vector<> para representar vectores 2- y 3-dimensionales con el álgebra usual, si bien ésta no cumple las condiciones anteriores impuestas por la biblioteca Boost.Numeric.Odeint.
Iniciaremos la implementación de la función principal main() con un mensaje de cabecera:
auto main() -> int
{
namespace ode = boost::numeric::odeint;
fmt::print("{:_^70}\n", "Plotting the Lorenz System");
La siguiente expresión lambda auxiliar, de nombre affirmative_answer_to, nos permitirá 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. La introducción de un input incorrecto producirá la emisión de un mensaje de error y una nueva solicitud de respuesta:
auto affirmative_answer_to = [](std::string_view question) -> bool {
auto typed_yes = bool{};
auto yes_or_no = [&typed_yes](std::string const& answer){
auto const lwr = boost::to_lower_copy(answer);
typed_yes = lwr == "y" or lwr == "yes";
return typed_yes or lwr == "n" or lwr == "no";
};
[[maybe_unused]] auto _ = terminal::prompt_init<std::string>(question, yes_or_no);
return typed_yes;
};
Sirviéndonos de la función anterior, proporcionaremos al usuario la opción de emplear los valores típicos de los parámetros σ = 10, ρ = 28 y β = 8/3 para los que el sistema exhibe comportamiento caótico (atractor de Lorenz), o bien establecer unos distintos. Éste y próximos bloques de código emplearán la técnica IILE (véase este post) a través de la cual una variable es inicializada mediante la invocación inmediata de una expresión lambda. En nuestro caso, la optimización NRVO debiera eliminar los posibles sobrecostes de tales operaciones:
auto gtr_than_zero = []<typename T> requires std::is_arithmetic_v<T> (T d){
return d > T{};
};
auto const [sigma, rho, beta] = /* IILE: parámetros del sistema */ [&]{
auto res = std::array<double, 3>{};
if (affirmative_answer_to(
"Will you use the default parameters σ = 10, ρ = 28, β = 8/3? [y/n]: ")
) {
res = {10.0, 28.0, 8.0/3.0};
}
else {
terminal::prompt("σ [>0.0]: ", res[0], gtr_than_zero);
terminal::prompt("ρ [>0.0]: ", res[1], gtr_than_zero);
terminal::prompt("β [>0.0]: ", res[2], gtr_than_zero);
}
return res;
}();
Asimismo, el programa preguntará al usuario si desea asignar una posición inicial p0 arbitraria alrededor del origen (único punto crítico asintóticamente estable cuando ρ < 1, no-único e inestable si ρ > 1) o si, por el contrario, prefiere establecer dichas coordenadas de forma manual:
auto const p_0 = /* IILE: punto de partida */ [&]{
auto res = R3_point{};
if (affirmative_answer_to(
"And a random initial position near (0,0,0)? [y/n]: ")
) {
auto rnd = [g = dlib::rand{std::time(nullptr)}] mutable {
return g.get_double_in_range(-5.0, 5.0);
};
res = {rnd(), rnd(), rnd()};
}
else {
terminal::prompt("x_0 in (-∞,+∞): ", res.x);
terminal::prompt("y_0 in (-∞,+∞): ", res.y);
terminal::prompt("z_0 in (-∞,+∞): ", res.z);
}
return res;
}();
Finalmente, solicitaremos al usuario que escoja un valor de tiempo orbital total Δt. La órbita será aproximada como una curva poligonal de segmentos rectos [e0, e1]. Éstos serán almacenados en un vector de nombre polyline cuyas entradas, de tipo dlib::perspective_window::overlay_line, podrán ser embebidas en un entorno de representación 3-dimensional. El algoritmo de integración numérica de EDOs que adoptaremos para obtener los vértices de los segmentos será el de Runge-Kutta-Dormand-Prince (RKDP, quinto orden), proporcionado por la biblioteca Boost.Numeric.Odeint. Este método realiza un control adaptativo del tamaño del paso dt con el fin de disminuir el error cometido [8]. En este caso, fijaremos las tolerancias de error absoluto y relativo a 1.0e-10 y estableceremos un tamaño máximo de paso de 5.0e-3 [5]:
auto const polyline =
/* IILE: segmentos de la órbita */ [&, q = p_0]
mutable {
auto res = std::vector<dlib::perspective_window::overlay_line>{};
// solicitamos el tiempo orbital total y reservamos memoria para un número
// mínimo de segmentos:
auto const time = terminal::prompt_init<
double>(
"Δt [>0.0]: ", gtr_than_zero);
auto const max_step_sz =
5.0e-3;
res.reserve(
static_cast<std::size_t>(time/max_step_sz));
// derivadas temporales de las coordenadas (miembros derecho de las EDOs):
auto lorenz = [&](R3_point
const& p, R3_point& ode_rhs,
[[maybe_unused]] double t){
auto const& [x, y, z] = p;
ode_rhs = {sigma*(y - x), x*(rho - z) - y, x*y - beta*z};
};
// inicializamos un 'stepper' con control adaptativo de paso que emplee RKDP:
using Rkdp = ode::runge_kutta_dopri5<R3_point,
double, R3_point,
double,
ode::vector_space_algebra>;
auto const rkdp = ode::make_controlled<Rkdp>(
1.0e-10,
1.0e-10, max_step_sz);
// calculamos y almacenamos en el vector los extremos de los segmentos [e_0, e_1],
// de tamaño diverso, partiendo de p_0 y cronometrando desde t = 0 hasta Δt:
auto [fst, lst] = ode::make_adaptive_range(rkdp, lorenz, q,
0.0, time, max_step_sz);
std::for_each(std::next(fst), lst, [&, e_0 = q](R3_point
const& e_1)
mutable {
res.emplace_back(dlib::vector<
double>{e_0.x, e_0.y, e_0.z},
dlib::vector<
double>{e_1.x, e_1.y, e_1.z});
e_0 = e_1;
});
// * insertar aquí el bloque de código que dota de color a los segmentos (ver abajo)
return res;
}();
|
La iteración inicial std::next en el bucle std::for_each evita la generación de un segmento innecesario de extremos coincidentes p_0 y longitud nula. El color adoptado por los segmentos producidos en dicho bucle es blanco por defecto, a la espera de conocer el número total de pasos realizados y poder dotar así a cada tramo de una tonalidad acorde al nivel de avance a lo largo de la curva. La animación proporcionada en la introducción de este artículo y la figura izquierda, por ejemplo, muestran el atractor de Lorenz correspondiente a los parámetros por defecto empleando el mapa de color dlib::colormap_heat. Sin embargo, con el fin de poder visualizar claramente todos los tramos de la órbita, es recomendable utilizar el mapa de color dlib::colormap_jet u otro alternativo a gusto del programador: |
// * bloque de código de color que insertar arriba:
for (
auto const num_segments = res.size();
auto [idx, segment] : res | ranges::views::enumerate
) {
segment.color = dlib::colormap_jet(idx,
0, num_segments);
}
En el código anterior hemos empleado la vista enumerate de la biblioteca Range-v3 para asignar un índice entero a cada segmento y, en base a su valor, adoptar un color específico.
Llegados a este punto, procederemos a abrir una ventana de perspectiva 3-dimensional dlib::perspective_window donde poder ubicar la órbita del sistema recién calculada. Entre otras acciones, ello nos permitirá rotar la representación para obtener una mejor perspectiva de la misma y/o realizar zoom sobre una región concreta de puntos. Con el fin de facilitar la visualización, incluiremos un sistema de referencia cartesiano cuyos ejes puedan distinguirse por sus diferentes colores:
fmt::print("Starting from ({:.2f}, {:.2f}, {:.2f})...", p_0.x, p_0.y, p_0.z);
// abrimos una ventana de perspectiva 3D:
auto wndw = dlib::perspective_window{};
wndw.set_title("Lorenz system - x axis (blue), y axis (green), z axis (white)");
wndw.set_size(640, 640);
// superponemos los ejes cartesianos:
auto const origin = dlib::vector<double>{}; // (0,0,0)
wndw.add_overlay(origin, {10.0, 0.0, 0.0}, dlib::rgb_pixel{0, 0, 255}); // eje x
wndw.add_overlay(origin, {0.0, 10.0, 0.0}, dlib::rgb_pixel{0, 255, 0}); // eje y
wndw.add_overlay(origin, {0.0, 0.0, 10.0}, dlib::rgb_pixel{255, 255, 255}); // eje z
// superponemos la curva poligonal:
wndw.add_overlay(polyline);
// aguardamos a que el usuario cierre la ventana:
wndw.wait_until_closed();
} // final de main()
A modo de ejemplo, la siguiente ejecución del programa genera una órbita estable para σ = 10, ρ = 14 y β = 8/3:
Fichero de construcción CMakeLists.txt
El siguiente fichero CMakeLists.txt puede ser empleado para construir el proyecto mediante CMake, bajo el supuesto de que el código de la sección anterior se encuentre en un archivo fuente de nombre lorenz.cpp. La biblioteca terminal.hpp y el fichero de configuración CMakeLists.txt deben ubicarse en el mismo directorio que lorenz.cpp.
cmake_minimum_required(VERSION 3.20)
project(lorenz_system
VERSION 0.1.0
DESCRIPTION "Plotting the Lorenz System - dgvergel.blogspot.com"
LANGUAGES CXX
)
add_executable(${PROJECT_NAME} lorenz.cpp terminal.hpp)
target_compile_features(${PROJECT_NAME} PRIVATE cxx_std_23)
set_target_properties(${PROJECT_NAME} PROPERTIES
CXX_EXTENSIONS OFF
RUNTIME_OUTPUT_DIRECTORY_DEBUG ${CMAKE_BINARY_DIR}/debug
RUNTIME_OUTPUT_DIRECTORY_RELEASE ${CMAKE_BINARY_DIR}/release
)
#-------------------------------------------------------------------
# integración de bibliotecas de terceros:
find_package(Boost 1.77.0 REQUIRED)
find_package(dlib 19.20 REQUIRED)
find_package(fmt 8.0 REQUIRED)
find_package(range-v3 0.11.0 REQUIRED)
target_link_libraries(${PROJECT_NAME} PRIVATE
Boost::headers
dlib::dlib
fmt::fmt
range-v3
)
#-------------------------------------------------------------------
# política de avisos específica de GCC:
target_compile_options(${PROJECT_NAME} PRIVATE -Wall -Wextra -pedantic -Werror)
#-------------------------------------------------------------------
# mensajes a emitir en construcción:
message("-- CMake Generator: ${CMAKE_GENERATOR}")
get_target_property(CFEATURES ${PROJECT_NAME} COMPILE_FEATURES)
message("-- Target compile features: ${CFEATURES}")
get_target_property(COPTIONS ${PROJECT_NAME} COMPILE_OPTIONS)
message("-- Target compile options: ${COPTIONS}")
Referencias bibliográficas
- Boost.Numeric.Odeint - https://www.boost.org/doc/libs/1_75_0/libs/numeric/odeint/doc/html/index.html
- {fmt} - https://fmt.dev/latest/index.html
- Range-v3 - https://ericniebler.github.io/range-v3/
- Dlib C++ library - http://dlib.net
- Boost.Numeric.Odeint - Controlled steppers - https://www.boost.org/doc/libs/1_75_0/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 - https://www.boost.org/doc/libs/1_75_0/libs/numeric/odeint/doc/html/boost_numeric_odeint/odeint_in_detail/state_types__algebras_and_operations.html
- Boost.Operators - Example - https://www.boost.org/doc/libs/1_75_0/libs/utility/operators.htm#example
- Boost.Numeric.Odeint - Stepper overview - https://www.boost.org/doc/libs/1_75_0/libs/numeric/odeint/doc/html/boost_numeric_odeint/odeint_in_detail/steppers.html#boost_numeric_odeint.odeint_in_detail.steppers.stepper_overview
No hay comentarios:
Publicar un comentario