It is known that dual representation of problems (through the main function and its adjoint in the Lagrange sense) makes it possible to formulate an effective perturbation theory on which the successive approximation method in the inverse problem theory relies. Let us suppose that according to preliminary predictions, a solution to the inverse problem (for example, the structure of medium of interest) belongs to a certain set A. Then selecting a suitable (trial, reference) element a0 as an unperturbed one and applying the perturbation theory, one can approximately describe the behavior of a solution to the forward problem in this domain and find a subset A0 that matches the measurement data best. However, as the accuracy requirements increase, the domain of applicability of the first approximation A0 is rapidly narrowing, and its expansion via addition of higher terms of the expansion complicates the solving procedure. For this reason, a number of works have searched for unperturbed approaches, including the method of variational interpolation (VI method). In this method, not one but several reference problems a1, a2,...,an are selected, from which a linear superposition of the principal function and the adjoint one is constructed, followed by determination of coefficients from the condition of stationarity of the form of the desired functional representation. This paper demonstrates application of the VI method to solving inverse problems of cosmic rays astrophysics in the simplest statement.