Utilities

#include <audi/invert_map.hpp>

std::vector<gdual_d> invert_map(const std::vector<gdual_d> &map_in, bool verbose = false)

Inversion of Taylor series.

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 equation via a Taylor polynomial truncated at order \(m\) (i.e. \(\in \mathcal P_{n,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 is called 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. Note that the inverse is guaranteed to exist if the linear part is invertible.

Parameters

• map_in – The input map represented as an std::vector of gdual<double>. They can have a constant coefficient which will be neglected

• verbose – when true some output is shown during the iterations

Returns

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

Exception

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.

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