In the optical region of frequencies, the excitation spectrum of interacting polariton waves in molecular crystals has been studied by means of the polariton Hamiltonian, which includes both cubic and quartic anharmonic terms with respect to the polariton operators. A general expression for Dyson's equation is developed which describes the system of interacting polariton fields and the polarization operator is a function of two-, three-, and mixed-polariton Green's functions. The polarization operator and, therefore, the polariton Green's function is then evaluated in successive approximations. In the zero approximation, the quartic anharmonic terms in the polariton Hamiltonian renormalize the frequencies of each polariton mode. The polarization operator is calculated in the lowest order via a zeroth-order renormalized Hamiltonian; the derived expression contains terms arising from the cubic and quartic anharmonicity and is a function of the renormalized frequencies of the polariton modes. In this approximation, the spectral function for the polariton spectrum is found to have a shape of an asymmetric Lorentzian line when the frequency dependence of the energy shift and spectral width is neglected. Resonances occur (i) when either two renormalized polaritons are created or one is created and the other is absorbed (two-polariton process), and (ii) when three renormalized polaritons are created or two are created and one is absorbed (three-polariton process) by a single incident renormalized polariton and vice versa. Formulas are derived for the probability amplitudes corresponding to the third- and fourth-order nonlinear scattering processes in the transparent range of frequencies of the crystal; an expression for the energy of excitation is developed where the frequency of the polariton mode is dressed by the field of all the others in the first approximation and propagates in the crystal without damping. A discussion is given on how the derived results can be improved by calculating the polarization operator in the next order in such a self-consistent manner as to become a function of the frequencies of the polariton modes, which are correct in the first approximation.