Numerical instabilities

Beutler, T. C., Mark, A. E., van Schaik, R. C., Gerber, P. R., van Gunsteren, W. F. Avoiding singularities and numerical instabilities in free energy calculations based on molecular simulations. Chem. Phys. Letters 222 (1994) 529-539... 

Buetler T C, A E Mark, R C van Schaik, P R Gerber and W F van Gunsteren 1994. Avoiding Singularities and Numerical Instabilities in Free Energy Calculations Based on Molecular Simulations. Chemical Physics Letters 222 529-539. 

Errors in advection may completely overshadow diffusion. The amplification of random errors with each succeeding step causes numerical instability (or distortion). Higher-order differencing techniques are used to avoid this instability, but they may result in sharp gradients, which may cause negative concentrations to appear in the computations. Many of the numerical instability (distortion) problems can be overcome with a second-moment scheme (9) which advects the moments of the distributions instead of the pollutants alone. Six numerical techniques were investigated (10), including the second-moment scheme three were found that limited numerical distortion the second-moment, the cubic spline, and the chapeau function. 

Fig. 6.7. The predicted, one-dimensional, mean-bulk temperatures versus location at various times are shown for a typical powder compact subjected to the same loading as in Fig. 6.5. It should be observed that the early, low pressure causes the largest increase in temperature due to the crush-up of the powder to densities approaching solid density. The "spike in the temperature shown on the profiles at the interfaces of the powder and copper is an artifact due to numerical instabilities (after Graham [87G03]).

Highly skewed cells should be avoided. The angles between the grid lines of a hexahedral mesh should be 90°. Angles < 40° or > 140° often imply a reduced accuracy or numerical instabilities. 

Especially for the electrons, the fluid model has the advantage of a lower computational effort than the PIC/MC method. Their low mass (high values of the transport coefficients) and consequent high velocities give rise to small time steps in the numerical simulation (uAf < Aa) if a so-called explicit method is used. This restriction is easily eliminated within the fluid model by use of an implicit method. Also, the electron density is strongly coupled with the electric field, which results in numerical Instabilities. This requires a simultaneous implicit solution of the Poisson equation for the electric field and the transport equation for the electron density. This solution can be deployed within the fluid model and gives a considerable reduction of computational effort as compared to a nonsi-multaneous solution procedure [179]. Within the PIC method, only fully explicit methods can be applied. 

For models described by a set of ordinary differential equations there are a few modifications we may consider implementing that enhance the performance (robustness) of the Gauss-Newton method. The issues that one needs to address more carefully are (i) numerical instability during the integration of the state and sensitivity equations, (ii) ways to enlarge the region of convergence. 

It should be emphasized here that it is unnecessary to integrate the state equations for the entire data length for each value of u. Once the objective function becomes greater than S(k0)). a smaller value for p can be chosen. By this procedure, besides the savings in computation time, numerical instability is also avoided since the objective function often becomes large very quickly and integration is stopped well before computer overflow is threatened. 

The transition from WA( ) = 1 to WA( ) = 0 as the distance from the nucleus A increases needs to be smooth enough such that numerical instabilities are avoided but at the same time also as abrupt as possible such that density peaks from nearby the nuclei are extinguished. The implementation of this concept in the three-dimensional space involves a special choice of coordinates - see Becke, 1988c, for details - but actually leads to a smoothened step function as schematically sketched in Figure 7-1 for the one-dimensional case. 

The surface-tension force Fa in Eq. (7) is smoothed and distributed into the thickness of the interface. In order to circumvent numerical instability, the fluid properties such as density and viscosity in the interface region are determined with a continuous transition ... 

Fig. 3a indicates that the bubble-rise velocity measured based on the displacement of the top surface of the bubble ( C/bt) quickly increases and approaches the terminal bubble rise velocity in 0.02 s. The small fluctuation of Ubt is caused by numerical instability. The bubble-rise velocity measured based on the displacement of the bottom surface of the bubble (Ubb) fluctuates significantly with time initially and converges to Ubt after 0.25 s. The overshooting of Ubb can reach 45-50 cm/s in Fig. 3a. The fluctuation of Ubb reflects the unsteady oscillation of the bubble due to the wake flow and shedding at the base of the bubble. Although the relative deviation between the simulation results of the 40 X 40 x 80 mesh and 100 x 100 x 200 mesh is notable, the deviation is insignificant between the results of the 80 x 80 x 160 mesh and those of the 100 X 100 x 200 mesh. The agreement with experiments at all resolutions is generally reasonable, although the simulated terminal bubble rise velocities ( 20 cm/s) are slightly lower than the experimental results (21 25 cm/s). A lower bubble-rise velocity obtained from the simulation is expected due to the no-slip condition imposed at the gas-liquid interface, and the finite thickness for the gas-liquid interface employed in the computational scheme.

Avoiding singularities and numerical instabilities in free energy calculations based on molecular simulations, Chem. Phys. Lett. 222 529 (1994). 

Each of these steps introduces numerical errors that must be carefully controlled (Jenny et al. 2001). In particular, taking spatial gradients in step (2) using noisy estimates from step (1) can lead to numerical instability. We will look at this question more closely when considering particle-field estimation below. 

The conclusions from this rather elementary survey of the symmetry constraint problem all point in the same general direction. The imposition of symmetry constraints (other than the Pauli principle) on a variationally-based model is either unnecessary or harmful. Far from being necessary to ensure the physical reality of the wave function, these constraints often lead to absurd results or numerical instabilities in the implementation. The spin eigenfunction constraint is only realistic when the electrons are in close proximity and in such cases comes out of the UHF calculation automatically. The imposition of molecular spatial symmetry on the AO basis is not necessary if that basis has been chosen carefully — i.e. is near optimum. Further, any breakdowns in the spatial symmetry of the AO basis are a useful indication that the basis has been chosen badly or is redundant. 

In real calculations, it is advisable to use normalized basis functions in order to avoid problems with numerical instabilities. Therefore, the overlap matrix elements, Su, are defined using normalized basis functions. After simplification, we obtain... 

Equation (54) can suffer from numerical instabilities. Therefore we adopt the immittance matrix approach. First we define the column vector of mode amplitudes p(p) in the same way as in Eq. (41). Then the immittance matrix is introduced by the relation... 

---