Numerical simulations are performed to simulate atmospheric perturbations observed at ionospheric heights just after the 2011 off the Pacific coast of Tohoku Earthquake. A time-dependent, two-dimensional, nonlinear, non-hydrostatic, compressible and neutral, numerical model is developed to reproduce the atmospheric perturbations. An impulsive upward surface motion is assumed as the source of the perturbations. Simulated atmospheric perturbations at 300-km altitude show remarkable agreement with oscillations observed in the ionospheric total electron content (TEC) when the source width is about 250 km. In the vicinity of the source, the acoustic resonance modes between the ground surface and the lower thermosphere are dominant. They have three dominant frequencies for the interval between 20 and 60 min after the impulsive input. The perturbation with the maximum amplitude has a frequency of 4.4 mHz. The other dominant modes have frequencies of 3.6 and 5.1 mHz. The beats between the dominant modes are also seen. In the distance, the gravity modes are dominant. The horizontal phase velocities are about 220 to 300 m/s, and the horizontal wavelengths are about 200 to 400 km. The good agreement between the simulation and the observations indicates that ionospheric oscillations generated by the earthquake are mainly due to the motion of the neutral atmosphere.