Utilities

#include <audi/invert_map.hpp>

template<typename Cf, typename Monomial> inline taylor_map<Cf, Monomial> audi::invert_map(const taylor_map<Cf, Monomial> &map_in, bool verbose = false)

Inversion of Taylor maps.

Consider a system of non linear equations in \(n\) variables:

\[\begin{split}\left\{ \begin{array}{l} f_1(x_1, x_2, .., x_n) = 0 \\ ... \\ f_n(x_1, x_2, .., x_n) = 0 \\ \end{array} \right. \end{split}\]

and represent each of it via a Taylor polynomial truncated at order \(m\) (i.e. \(\in \mathcal P_m\)) around some point \(\overline {\mathbf x}\):

\[\begin{split}\left\{ \begin{array}{l} T_{f_1}(dx_1, dx_2, .., dx_n) = 0 \\ ... \\ T_{f_n}(dx_1, dx_2, .., dx_n) = 0 \\ \end{array} \right. \end{split}\]

The above system of, now, polynomial equations can be written as \(C + M + N = C + \mathcal M = 0\), where we have explicitly indicated the constant term \(C\), the linear term \(M\) and the rest \(N\). The symbol \(\mathcal M\) indicates what we call a Taylor map and its inverse is found by this function by performing some Picard iterations (fixed point) over the following simple formal equalities:

\[\mathcal M \circ \mathcal M^{-1} = \mathcal I \rightarrow (M + N) \circ \mathcal M^{-1} = \mathcal I \rightarrow M \circ \mathcal M^{-1} + N \circ \mathcal M^{-1} = \mathcal I \]

and hence:

\[ \mathcal M^{-1} = M^{-1} \circ \left(\mathcal I - N \circ \mathcal M^{-1}\right) \]

which can be used to make a fixed point iterative scheme once the inverse of the linear part \(M^{-1}\) is found using any linear algebraic tool. Note that the inverse is guaranteed to exist if the linear part is invertible.

Example usage:

gdual<double> x(0., "x1", 9);
auto sinx = sin(x);
auto asinx = asin(x);
auto inv_sinx = invert_map({sinx});
std::cout << "sin x:\t\t" << sinx << std::endl;
std::cout << "asin x:\t\t" << asinx << std::endl;
std::cout << "inv sinx:\t" << inv_sinx[0] << std::endl;

The code above will produce the output:

sin x:        dx1+2.75573e-06*dx1**9-0.166667*dx1**3+0.00833333*dx1**5-0.000198413*dx1**7
asin x:       dx1+0.0303819*dx1**9+0.166667*dx1**3+0.075*dx1**5+0.0446429*dx1**7
inv sinx:     dp0+0.0303819*dp0**9+0.166667*dp0**3+0.075*dp0**5+0.0446429*dp0**7

Parameters:

map_in – The input map represented as an std::vector of gduals . They can have a constant coefficient which will be neglected
verbose – when true some output is shown during the Picard iterations

Throws:

std::invalid_argument – if map_in is empty, if its linear part is not invertable, if the symbol set of the map components are not all equal, if the map is not square (i.e. it is of size n with n symbols) or if the order of the gduals are not all the same.

Returns:

The inverse map \(\mathcal M^{-1}\) represented as an std::vector of gduals with symbol set [p0, p1, .. pn]