A set of integro-differential equations in the Lagrangian renormalized approximation (Kaneda 1981) is rederived by applying a perturbation method developed by Kraichnan (1959), which is based upon an extraction of direct interactions among Fourier modes of a velocity field and was applied to the Eulerian velocity correlation and response functions, to the Lagrangian ones for homogeneous isotropic turbulence. The resultant set of integro-differential equations for these functions has no adjustable free parameters. The shape of the energy spectrum function is determined numerically in the universal range for stationary turbulence, and in the whole wavenumber range in a similarly evolving form for the freely decaying case. The energy spectrum in the universal range takes the same shape in both cases, which also agrees excellently with many measurements of various kinds of real turbulence as well as numerical results obtained by Gotoh et al. (1988) for a decaying case as an initial value problem. The skewness factor of the longitudinal velocity derivative is calculated to be −0.66 for stationary turbulence. The wavenumber dependence of the eddy viscosity is also determined.