Artículos de la serie:
- Representando el atractor de Lorenz con C++20, Boost.Numeric.Odeint y Dlib
- Análisis dimensional con mp-units, regresiones no-lineales con Dlib.Optimization y visualización de datos con Matplot++
- Rectas de mejor ajuste por mínimos cuadrados (mp-units, Matplot++)
- Ceros de una función
- Comprobando la ley de Zipf
- Instalación de Matplot++ y mp-units con Mingw-w64
En este post analizaremos la obtención del cero de una función F(x) = 0 empleando el conjunto de herramientas matemáticas Math Toolkit de la familia de bibliotecas Boost [1]. En primera aproximación, descartaremos la presencia de múltiples ceros. A modo de ejemplo, consideremos la denominada ecuación de Colebrook-White, que relaciona la pérdida de presión por rozamiento a lo largo de una tubería en términos de la velocidad media del flujo del fluido:
Aquí, 𝑓 denota el factor de fricción –una medida adimensional de la caída de presión—, 𝜖 es la rugosidad de la tubería, Dh su diámetro hidráulico y Re el número de Reynolds que caracteriza el movimiento del fluido. En particular, determinaremos el factor de fricción 𝑓 para 𝜖 = 0.00015 m, Dh = 0.10 m y Re = 5000 (régimen turbulento). Para ello, realizaremos el cambio de variable 𝑓 = 1/x2 con el fin de reducir el número de iteraciones y mejorar la estabilidad del algoritmo.
Nos proponemos realizar la búsqueda de ceros, en primer lugar, sin necesidad de recurrir al cálculo de derivadas. El programa resulta inmediato de codificar:
#include <array>
#include <cmath>
#include <cstdlib>
#include <limits>
#include <utility>
#include <boost/math/tools/roots.hpp>
#include <fmt/format.h>
auto main() -> int
try {
namespace bmt = boost::math::tools;
auto const epsilon = 0.00015, // en metros
dh = 0.10, // en metros
re = 5'000.0;
auto colebrook_white_eq = [&](double x){
return x + 2.0*std::log10(epsilon/(3.71*dh) + (2.51/re)*x);
};
auto const [x1, x2] = /* IILE */ [&]{
auto const [guess, factor] = std::array{1.0, 2.0};
auto const rising = true;
auto const tol = bmt::eps_tolerance<double>{std::numeric_limits<double>::digits - 3u};
auto const maxit = boost::uintmax_t{100}; // número máximo de iteraciones
auto it = maxit;
auto const res = bmt::bracket_and_solve_root(colebrook_white_eq,
guess, factor, rising, tol, it);
if (it >= maxit)
throw std::runtime_error{fmt::format("{} iterations exceeded", maxit)};
return res;
}();
auto const x = x1 + (x2 - x1)/2.0;
fmt::print("Numerical solution: x = {:.6f} f = {:.6f}\n", x, 1/(x*x));
}
catch (std::exception const& e) {
fmt::print(stderr, "\nexception occurred: {}", e.what());
return EXIT_FAILURE;
}
El método bracket_and_solve_root [2] aplica el algoritmo TOMS 748 para encontrar el cero de la función colebrook_white_eq. Se cumplen aquí los prerrequisitos exigidos por el método, a saber:
- colebrook_white_eq es una función unaria necesariamente monotónica en la mitad del eje real que contiene a la primera estimación (guess) proporcionada al algoritmo. En este caso, la función es creciente (en todo su dominio), punto éste que comunicamos al algoritmo estableciendo como true el booleano rising. Esta información es utilizada, junto al signo de colebrook_white_eq(guess), para determinar si el cero de la función es mayor o menor que guess.
- La estimación inicial guess debe poseer el mismo signo que la raíz que se desea encontrar, en este caso positiva.
La constante factor establece un factor de escala utilizado para situar el cero de la función dentro de un intervalo. Así, el valor estimado se multiplicará (o se dividirá, según corresponda) por dicho factor hasta encontrar los extremos del intervalo óptimo.
Nuestro programa emitirá una excepción de llegar a superarse el número máximo maxit de iteraciones permitidas. En cualquier otro caso, obtendremos un intervalo [x1,x2] conteniendo el cero de la función, verificándose:
colebrook_white_eq(x1)*colebrook_white_eq(x2) <= 0.0
Se cumple también que tol(x1,x2) == true, siendo tol una lambda binaria que determina la condición de terminación para la búsqueda del cero de la función. Éste es finalmente estimado como punto medio del intervalo [x1,x2]. La ejecución del programa conduce, así, a x = 5.061739 (𝑓 = 0.039030).
Como segundo análisis del problema planteado, consideremos ahora la búsqueda del cero F(x) = 0 descansando en el cálculo de la derivada primera F' de la función. Podemos realizar el cálculo empleando el método iterativo de Newton-Raphson [3] en la forma siguiente:
auto colebrook_white_eq_and_derivative = [&](double x){
return std::pair{
x + 2.0*std::log10(epsilon/(3.71*dh) + (2.51/re)*x), // función
1.0 + 2.0*3.71*2.51*dh/(std::log(10)*(2.51*3.71*dh*x + epsilon*re)) // derivada
};
};
auto const xs = /* IILE */ [&]{
auto const [guess, min, max] = std::array{1.0, 0.0, 15.0};
auto const digits = static_cast<int>(0.6*std::numeric_limits<double>::digits);
auto const maxit = boost::uintmax_t{100};
auto it = maxit;
auto const res = bmt::newton_raphson_iterate(colebrook_white_eq_and_derivative,
guess, min, max, digits, it);
if (it >= maxit)
throw std::runtime_error{fmt::format("{} iterations exceeded", maxit)};
return res;
}();
fmt::print("Numerical solution: x = {:.6f} f = {:.6f}", xs, 1/(xs*xs));
En este caso, observemos que la lambda colebrook_white_eq_and_derivative retorna una pareja de valores: el valor de la función evaluada en su argumento junto a su derivada primera. Las constantes min y max establecen los valores mínimo y máximo esperados para la raíz, respectivamente. Por su parte, digits indica el número deseado de dígitos de precisión. Los resultados numéricos alcanzados en este caso son coincidentes con los mostrados en la primera sección del post.
Referencias bibliográficas
- Boost - Root Finding & Minimization Algorithms - https://www.boost.org/doc/libs/1_76_0/libs/math/doc/html/root_finding.html
- Boost - Bracket and Solve Root - https://www.boost.org/doc/libs/1_76_0/libs/math/doc/html/math_toolkit/roots_noderiv/bracket_solve.html
- Boost - Root Finding With Derivatives: Newton-Raphson, Halley & Schröder - https://www.boost.org/doc/libs/1_76_0/libs/math/doc/html/math_toolkit/roots_deriv.html
No hay comentarios:
Publicar un comentario