From the data on Cherenkov light from EASs collected over the period of more than 15 years, the spectrum of cosmic rays was obtained in the energy range 1.0E15-1.0E18 eV. This spectrum has two features: at ~3.0E15 eV (first knee) and at ~1.0E17 eV (second knee). The first knee is characterized by the indexes \gamma_1=2.7+/-0.03 and \gamma_2 = 3.12+/-0.03 and the second knee - by \gamma_1 = 2.92+/-0.03 and \gamma_2=3.24+/-0.04. In the second case the difference amounts to \Delta\gamma_23=0.32{\pm}0.03{\pm}0.05 which is less than the case of the first knee \Delta\gamma_12=0.42+/-0.03+/-0.05 (here dual errors represent statistics and systematics correspondingly). The lesser difference \Delta\gamma_23 can be explained with the influx of cosmic rays from Meta-galaxy and, hence, with some increase of the cosmic rays intensity in the energy range 5.0E16-1.0E18 eV, which compensates the escape of heavier nuclei from the Galaxy. The presence of the second knee could be confirmed by abrupt change in mass composition from from lnA~3 at ~1.0E17 eV to lnA~1.5 at 1.0E18eV found from the analysis of longitudinal development of EAS.

It is proposed a new approach for estimating the composition of cosmic rays. It is found that the zenith angle distributions and muon components of EAS for energies E>10^19 eV and E>4.10^19 eV differ from each other. It is shown that the cosmic rays above E>4.10^19 eV is heavier than the cosmic rays at energy E~ 10^19 eV. According to our estimation the SUGAR array detected 8 showers above 10^20 eV. It is concluded that no sign of Greisen-Zatsepin-Kuzmin (GZK) cut off in the spectrum of cosmic rays and all cosmic rays are galactic.

The ratio of the muon flux density to charged particle flux density at distances of 300 and 600 m from the shower axis ($\rhom(300)/\rhos(300)$ and $\rhom(600)/\rhos(600)$) is measured. In addition, the energy dependence of $\rhom(1000)$ is analysed for showers with energies above $10^{18}$ eV. A comparison between the experimental data and calculations performed with the QGSJET model is given for the cases of primary proton, iron nucleus and gamma- ray. We conclude that the showers with $\E\ge3\times10^{18}$ eV can be formed by light nuclei with a pronounced fraction of protons and helium nuclei. It is not excluded however that a small part of showers with energies above $10^{19}$ eV could be initiated by primary gamma-rays.

At the Yakutsk EAS array E_em is determined by using measurements of EAS Cherenkov light flux and charged particle flux. It is known from calculations that these characteristics depend on a sort of primary particle and, therefore, the estimation of E_em depends on a primary particle mass. In the work the dependence of the E_em/E_0 ratio on the energy is given and experimental data are compared with calculations by the QGSJET model. The calculations have been carried out for the primary proton and iron nucleus. The average calculated meaning of the value of E_em/E_0 ratio (between the proton and iron nucleus) within experimental errors is in agreement with experimental data that doesnt contradict to the mixed mass composition of primary cosmic radiation.

The experimental distributions of $\Xmax$ obtained with the Yakutsk EAS array at fixed energies of $5 \times10^{17}$, $1\times10^{18}$ and $5\times10^{18}$ eV are analysed. A recent version of the QGSJET model is used as a tool of our analysis. In the framework of this model, the most adequate mass composition of primary particles satisfying the experimental data on $\Xmax$ is selected.

The energy spectrum of primary cosmic rays with ultra high energies based on the Yakutsk EAS Array data is presented. For the largest events values of S600 and axis coordinates have been obtained using revised lateral distribution function. The affect of the arrival time distribution at several axis distance on estimated density for Yakutsk and AGASA is considered.

There are some discrepancies in the results on energy spectrum from Yakutsk, AGASA, and HiRes experiments. In this work differential energy spectrum of primary cosmic rays based on the Yakutsk EAS Array data is presented. For the largest events values of $S_{600}$ and axes coordinates have been obtained using revised lateral distribution function. Simulation of converter's response at large distances showed no considerable underestimation of particle density. Complex shape of spectrum in region of $\E > 10^{17}$ eV is confirmed. After adjustment of parameters and additional exposition at the Yakutsk array there are three events with energy $\E > 10^{20}$ eV.

In the work the data of the Yakutsk complex EAS array and their comparison with calculation in the case of primary nuclei of different chemical elements are presented. The calculation by QGSJET model have been used interpreting experimental data.

The energy spectrum of UHECRs is discussed on the basis of the Yakutsk array database analysis. In the region E=0.1 to 30 EeV the showers are detected under trigger-500, while at energies above 30 EeV the whole acceptance area for trigger-1000 is used in order to utilize all the data available in the region of GZK cutoff.

Arrival directions of extensive air showers by using world data are considered. It is shown that distributions of showers in zenith angle at E>10^{19}eV and E>4.10^{19}eV differ from each other. On this basis the conclusions are made: the estimation of the energy of inclined showers at E>4.10^{19}eV at Yakutsk and AGASA arrays is incorrect; the cosmic ray spectrum is not break up and most likely continued up to E=10^{20}eV and higher. It is found that ultrahigh energy cosmic rays are, probably, superheavy nuclei.