An explicit method can suffer from numerical instability in case when Courant number is sufficiently large. The Courant number is inversely proportional to square of grid spacing. This case shows the occurrence of numerical instability that starts at spots where the grid spacing is exceptionally small. Finally, the code gives "floating point overflow" numerical error.
Equation to solve
+ 10∆Ψ = 0.
Cylinder. The radius of the cylinder is R=1 and the length L=4.
A mesh for a cylinder constructed with use of blockMesh utility.
The cylinder is divided in four prism-like blocks and a central block, which is a rectangular parallelepiped with a 16×16×32 grading. Four side blocks are curvilinear hexahedrons with a similar grading.
The function at zero time is Gaussian distribution,
Ψ(x,y,z) = exp(−16 [(x−0.4)2−y2−z2]).
Blue color corresponds to Ψ = 0, and red to Ψ = 1. The cylinder is sliced in order to show the interior.
Neumann (zero gradient normal to the boundary).
The time step is ∆t=1.8·10−5. Calculation time is around 1 minute.
Numerical divergence caused by large Courant number. The sub-figures show calculated absolute value of the wave function. Numerical instability starts to develop at t=0.002 in a region where the grid spacing is smallest. Then, it spreads to an entire region.