I have a faced a research problem where I would need to recover a frame field given its connection forms. More precisely, I begin with an orthonormal frame field (given by data) in $\Re^3$ written as $$ \mathbf F=\begin{pmatrix}\vec f_1\\\vec f_2\\\vec f_3\end{pmatrix} $$ where $\vec f_i:\Re^3\rightarrow\Re^3$ are vectors fields with $\vec f_i\cdot\vec f_j=\delta_{ij}$, and $\delta_{ij}$ is the Kronecker delta. I then obtain the connection forms $$ \omega_{ij}=\text{d}\vec f_i\cdot\vec f_j, $$ which yields $$ \text{d}\mathbf F=\Omega\mathbf F $$ where $\Omega=[\omega_{ij}]\in\Re^{3\times3}$ is the skew-symmetric matrix of connection forms.
In my application, I then proceed by computing the interior product of the 1-forms $\omega_{ij}$ onto the frame fields themselves (the Christoffel symbols), i.e., I compute $$ \omega_{ijk}\equiv\omega_{ij}\langle\vec f_k\rangle\in\Re $$ where $\langle\cdot\rangle$ denotes the standard interior product for forms. I thus obtain 9 different measurements at each point, i.e., $\omega_{121},\omega_{122},\omega_{123},\omega_{131},\omega_{132},\omega_{133},\omega_{231},\omega_{232},\omega_{232}$, each of which has a very precise meaning in the application at hand.
Now, I would like to do the converse, i.e., I would like to solve for $\mathbf F$ given the list of interior contractions $\omega_{121},\omega_{122},\omega_{123},\omega_{131},\omega_{132},\omega_{133},\omega_{231},\omega_{232},\omega_{232}$.
The Frobenius theorem states the unique existence of $\mathbf F$ in the neighborhood of $0$ if we set $\mathbf F(0)=I$ and if the following are satisfied: $$ \Omega=\text{d}(F)F^{-1} $$ $$ \text{d}\Omega-\Omega^2=0. $$ but I'm unsure about where to start. Would anybody have a suggestion on how to approach this problem? Is there a formal name to what I'm trying to do?
Perhaps a better way to formulate the problem is to enumerate the frame axis differentials directly: \begin{align} \text{d}\vec f_1&=\omega_{12}\vec f_2+\omega_{13}\vec f_3\\ \text{d}\vec f_2&=-\omega_{12}\vec f_1+\omega_{23}\vec f_3\\ \text{d}\vec f_3&=-\omega_{13}\vec f_1-\omega_{23}\vec f_2. \end{align}
Writing the direction of contraction as $$\vec v=\sum_i^3(\vec v\cdot\vec f_i)\vec f_i=\sum_i^3v_i\vec f_i$$ we obtain \begin{aligned} \text{d}\vec f_1\langle\vec v\rangle&=\left(\sum_i^3v_i\omega_{12i}\right)\vec f_2+\left(\sum_i^3v_i\omega_{13i}\right)\vec f_3\\ \text{d}\vec f_2\langle\vec v\rangle&=-\left(\sum_i^3v_i\omega_{12i}\right)\vec f_1+\left(\sum_i^3v_i\omega_{13i}\right)\vec f_3\\ \text{d}\vec f_3\langle\vec v\rangle&=-\left(\sum_i^3v_i\omega_{13i}\right)\vec f_1-\left(\sum_i^3v_i\omega_{23i}\right)\vec f_2. \end{aligned}
Which gives us three coupled differential equations to solve.
You are in fact asking how to integrate the following vector differential equation: $$ df = \Omega f, $$
where $f$ is, say, vector $f_1$ (to get $f_2$ and $f_3$ just use different initial conditions). If I understand correctly, you would like to find a matrix function $M$ such that: $$ f(x) = M(x)f(0). $$
On $\Bbb{R}^3$, the first equation is integrable if and only if both members are exact forms. So, if and only if: $$ 0=d(\Omega f) = d\Omega\, f - \omega\,df = d\Omega\,f - \Omega(\Omega\,f) = (d\Omega -\Omega^2)f. $$
This is what Frobenius's formula means. Now, $\Bbb{R}^3$ is flat, and it should be no trouble checking that this condition is indeed satisfied.
So how do you find $M$ from $\Omega$? What you do is called either product integral or ordered exponential (it is a good exercise to check that they are in fact equivalent).
Intuitively, consider a curve $t\mapsto \gamma(t)$ on $\Bbb{R}^3$. Since we have seen that our equation is integrable, only the endpoints matter. Suppose have our vectors of our frame field at $\gamma(0)$. We want to transport them, using our connection form matrix $\Omega$, to vectors at $\gamma(t)$.
Now after a small time $\Delta t$, our vector will have slightly rotated, so: $$ f(\Delta t) \approx \big( 1 + \Omega\,\Delta t \big)f(0). $$
Therefore, our vector at time $t$ will be the composition of "many" such small rotations: $$ f(t) \approx \big( 1 + \Omega(t_n)\,\Delta t \big)\big( 1 + \Omega(t_{n-1})\,\Delta t \big)...\big( 1 + \Omega(t_1)\,\Delta t \big)f(0). $$
If $\Delta t$ is small enough, we can write a matrix exponential: $$ \big( 1 + \Omega(t)\,\Delta t \big) \approx \exp\big( \Omega(t)\,\Delta t \big). $$
So we get: $$ f(t) \approx \exp\big( \Omega(t_n)\,\Delta t \big)...\exp\big(\Omega(t_1)\,\Delta t \big)f(0). $$
The matrix exponential has the rule $e^Ae^B=e^{A+B}$ only if $A$ and $B$ commute. So if the $\Omega_{t_i}$ all commute with each other, the product of exponentials becomes an exponential of a sum: $$ f(t) \approx \exp\big( \Omega(t_n)\,\Delta t +...+ \Omega(t_1)\,\Delta t \big)f(0), $$
so in the limit we get exactly an exponential of an integral: $$ f(t) \approx \exp\bigg( \int_0^t \Omega(t)\,dt \bigg)f(0). $$
This is, for example, what happens in a $2$-dimensional space. In 3 dimensions, in general, the $\Omega(t_i)$ do not commute with each other, and the product above does not turn into an ordinary integral. What you get is a special generalization of the exponential, called the ordered exponential. Click on the link above to see how to work with it.
So in the end: $$ M(x) = Pexp\bigg( \int_0^x \Omega \bigg)f(0). $$ In physics, it is used to reconstruct the global rotation after some time, when you have the angular velocity expressed as a function of $t$ (can you see why the problem is very similar to ours?).