Differential holonomic constraints Differential holonomic constraint is an integrable homogeneous first order differential equation: 
$$\sum_{i}\mathcal{E}_{i}(q)\frac{dq_{i}}{d\tau}=0;$$ 
in which $\sum_{i}\mathcal{E}_{i}(q)dq_{i}$ is a complete differential, that is the conditions  $$\partial\mathcal{E}_{i}(q)/\partial q_{j}= \partial\mathcal{E}_{j}(q)/\partial q_{i}$$ are met for all $i, j$.
In case when the potential function $E(q)$ for the corresponding form, such that: $$\mathcal{E}_{i}(q)=\partial E(q)/\partial q_{i};$$  is known, the general one-parametric solution of the above differential equation is obviously: $$E(q_{1},...q_{n})=Const;$$ 


*

*What can we say about the solution of this differential equation when the potential function is not known, can not be given in a closed explicit form (does not exist)?
Can we still expect in this case that the dimension of the configuration space is effectively reduced by this constraint, or maybe even there still exist in some implicit sense a general solution, as a one-parametric family of relations between variables $q_{i}$?

*If so, than can we find / build those solutions?

*How to deal with such a constraint? We can not just add it in its differential form to the Lagrangian as a Lagrangian multiplier, as we can do for both explicit holonomic and non-holonomic constraints.    
I appreciate your suggestions, specific references.
 A: The case that $E$ does not exist cannot happen on topologically trivial configuration spaces, and even then, it exists locally by the Poincaré lemma.
Your "complete conditions" say nothing but that the form $\mathcal{E}=\mathcal{E}_i\mathrm{d}q^i$ is closed, i.e. $\mathrm{d}\mathcal{E} = 0$. By the Poincaré lemma, every closed form is locally exact, i.e. there is some function $E$ such that $\mathcal{E} = \mathrm{d}E$.
So, on a topologically non-trivial configuration space, it can happen that a constraint given by a closed form is not exact, so it will not be holonomic.
This doesn't say much about your constraint, it says something about the configuration space - it has topological defects for which the system will change when it travels "around" them. If your configuration space is topologically non-trivial, many of the standard methods only work locally. However, finding local solutions to the equations of motion is enough, since you can glue them together. So, just restrict to a contractible neighbourhood around your initial conditions, and write the Lagrangian with the constraint added by a Lagrange multiplier as is standard. Solve those local equations of motion. When you near the boundary of that neighbourhood, take the endpoint of this local solution to the e.o.m., and choose another contractible neighbourhood around that endpoint to propagate the solution. Repeat.
