// TestExp5.cpp
//
// Using ODEs for optimisation. The very very simplest example to show how it works in the scalar
// case.
// DJD
// Systems? this is an execise.
//
//
// http://www.unige.ch/~hairer/preprints/gradientflow.pdf
//
// dx/dt = - nabla(U(x)) - lambda * nabla(c*x))
//
// Transform f(x) = 0 to a least-squares problem.
//
// (C) Datasim Education BV 2018-2022
//
#include <iostream>
#include <vector>
#include <cmath>
#include <complex>
#include <boost/numeric/odeint.hpp>
using value_type = double;
using state_type = std::vector<value_type>;
template <typename T>
using FunctionType = std::function<T(const T& arg)>;
using CFunctionType = FunctionType<std::complex<value_type>>;
value_type CSM(const CFunctionType& f, value_type x, value_type h)
{ // df/dx at x using tbe Complex step method
std::complex<value_type> z(x, h); // x + ih, i = sqrt(-1)
return std::imag(f(z)) / h;
}
// We choose min (f(x) - a)*(f(x)-a), f(x) = exp(x)
// Change the value of a to suit needs
const std::complex<value_type> a(std::exp(5.0), 0.0);
auto func = [](const std::complex<value_type>& z)
{
return (z - a) * (z - a);
};
value_type h = 1.0e-9;
// Free function to model RHS in dy/dt = RHS(t,y)
void Ode(const state_type &x, state_type &dxdt, const value_type t)
{
// dx/dt = - grad(f(x)
// We need to compute gradient
// 1. Exact
//
value_type a = std::exp(5.0);
dxdt[0] = 2.0 * (x[0] - a);
// 2. Complex Step method
//dxdt[0] = CSM(func, x[0], h);
// Transform semi-infinite time to (0,1)
// I use this technique in many places, dx/dtau = dx/dt (dt/dtau)
value_type tau = (1.0 - t) * (1.0 - t);
dxdt[0] = -dxdt[0]/tau;
}
int main()
{
namespace Bode = boost::numeric::odeint;
// Initial condition
value_type L = 0.0;
value_type T = 0.99;
value_type dt = 0.01001;
// Initial condition
state_type x{0.8};
// Middle-of-the-road solver
std::size_t steps = Bode::integrate(Ode, x, L, T, dt);
std::cout << "Number of steps Cash Karp 54: " << std::setprecision(16) << steps
<< "approximate: " << x[0] << '\n';
{
state_type x{0.8};
// Adaptive solver
Bode::bulirsch_stoer<state_type, value_type> bsStepper; // O(variable), controlled stepper
steps = Bode::integrate_const(bsStepper, Ode, x, L, T, dt);
std::cout << "Number of steps Bulirsch-Stoer: " << steps << ", approximate: " << x[0] << '\n';
}
return 0;
}
