We study how well the Gaussian approximation is valid for computing the covariance matrices of the convergence power and bispectrum in weak gravitational lensing analyses . We focus on its impact on the cosmological parameter estimations by comparing the results with and without non-Gaussian error contribution in the covariance matrix . We numerically derive the covariance matrix as well as the cosmology dependence of the spectra from a large set of N -body simulations performed for various cosmologies and carry out Fisher matrix forecasts for tomographic weak lensing surveys with three source redshifts . After showing the consistency of the power and bispectra measured from our simulations with the state-of-the-art fitting formulas , we investigate the covariance matrix assuming a typical ongoing survey across 1500 deg ^ { 2 } with the mean source number density of 30 arcmin ^ { -2 } at the mean redshift z _ { s } = 1.0 . Although the shape noise contributes a significant fraction to the total error budget and it mitigates the impact of the non-Gaussian error for this source number density , we find that the non-Gaussian error degrades the cumulative signal-to-noise ratio up to the maximum multipole of 2000 by a factor of about 2 ( 3 ) in the power ( bi- ) spectrum analysis . Its impact on the final cosmological parameter forecast with 6 parameters can be as large as 15 \% in the size of the one-dimensional statistical error . This can be a problem in future wide and deep weak lensing surveys for precision cosmology . We also show how much the dark energy figure of merit is affected by the non-Gaussian error contribution and demonstrate an optimal survey design with a fixed observational time .