We introduce a special class of functions for mathematical modeling of periodic signals of special shape with irregularly spaced arguments. This method was developed for determination of phenomenological characteristics of the light curves, which are necessary for registration in the "General Catalogue of Variable Stars" and other databases. For eclipsing binary stars with smooth light curves - of types EB and EW - it is recommended a trigonometric polynomial of optimal degree in a complete or symmetric form. For eclipsing binary systems with narrow minima (EA-type), statistically optimal is an approximation in a class of non-polynomial spline functions. It is used a combination of the second-order trigonometric polynomial (TP2, what describes effects of "reflection", "ellipsoidality" and "spotness") and localized contributions of minima (parametrized in depth and profile separately for primary and secondary minima). Such an approach is characterized by a statistical accuracy of the smoothing curve, which is up to ~1.5-2 times better than the trigonometric polynomial of statistically optimal degree, and the absence of false "waves" in the light curve associated with the effect of Gibbs. In addition to the minimum width, which can not be determined by a trigonometric polynomial approximation, the method allows to determine with better accuracy its depth, as well as to separate the effects of the eclipse and out-of-eclipse parts. For multi-color observations, improving the accuracy of the smoothing of the curve in each filter will allow to obtain with better accuracy the curves of the color index variations. Effectivity of the proposed method increases with decreasing eclipse depth. The method called NAV ("New Algol Variable"), was applied to eclipsing binary systems VSX J022427.8-104034=USNO-B1.0 0793-0023471 and BM UMa. For VSX0224, an alternative model of "double period" is discussed.