Numerical Wave Simulation Is Mostly About Not Being Wrong
If you want to simulate wave phenomena using mathematical methods, you need to understand partial differential equations first, then learn how to discretize them without introducing errors that make your results meaningless. This is not a glamorous topic. It is also not difficult once you stop trying to find shortcuts. The core tools—finite difference, finite element, and spectral methods—are well understood. The hard part is choosing the right one and knowing when your mesh is too coarse or your time step is too large. The wave equation in its basic form is second order in both space and time. In one dimension it looks like u_tt = c^2 * u_xx. The constant c is the wave speed. That is the starting point. Everything after that is about approximation. I learned this early on when I was working with acoustic simulations for a research group in 2019, and I still see people miss the basics years later. I recommend starting with finite differences because they are the most transparent. You replace derivatives with difference quotients on a grid. The standard approach uses a central difference in both space and time. For the spatial derivative you write (u[i+1] - 2u[i] + u[i-1]) / dx^2. For the temporal derivative you write (u[t+1] - 2u[t] + u[t-1]) / dt^2. You then solve for u[t+1] using values from the two previous time steps.
This gives you an explicit scheme. The advantage is simplicity. The disadvantage is the Courant-Friedrichs-Lewy condition, which says your time step must satisfy dt
= dx / c in one dimension. If you violate this, your simulation blows up almost immediately. I have seen students ignore this constraint for hours before giving up. The fix is to reduce dt or coarsen your mesh, but coarsening the mesh introduces dispersion errors that you will notice as waves traveling at the wrong speed.
Spectral Methods Are More Accurate But More Expensive
When you need higher accuracy and your domain is simple, spectral methods are worth considering. Instead of representing the solution on a grid of point values, you represent it as a sum of global basis functions, usually trigonometric or polynomial. The discrete Fourier transform becomes your main tool. A pseudo-spectral method evaluates nonlinear terms in physical space and linear terms in spectral space, which is the standard trick. The accuracy gain over finite differences is real. With the same number of degrees of freedom, a spectral method can give you orders of magnitude better accuracy for smooth problems. The catch is that spectral methods assume periodic or very smooth boundary conditions. If your solution has discontinuities or sharp gradients, Gibbs phenomena will appear and those oscillations can destabilize everything. I ran into this exact issue when simulating shock waves in a duct. The oscillations grew until the simulation crashed. My workaround was to add a small artificial viscosity term, something like epsilon * u_xx with epsilon around 1e-5, which dampened the high-frequency noise without affecting the main solution.
Finite Element Methods Handle Complex Geometry
When your domain has irregular boundaries, finite element methods are the standard choice. You divide the domain into elements, usually triangles or quadrilaterals in 2D, and represent the solution as a piecewise polynomial on each element. The weak form of the wave equation is where this approach lives. You multiply by a test function, integrate by parts, and end up with a system that involves mass and stiffness matrices. The mass matrix can be either consistent or lumped. A lumped mass matrix is diagonal and allows explicit time stepping without solving a linear system at each step. The downside is reduced accuracy, especially for higher modes. I spent two weeks debugging a simulation where the eigenfrequencies were off by about eight percent compared to an analytical solution. The problem traced back to using a consistent mass matrix with a time integrator designed for lumped matrices. Switching to a lumped mass matrix fixed it immediately, but you lose some precision in the process. That is a real tradeoff you need to be aware of.
Time Integration Schemes Matter More Than You Think
People focus heavily on spatial discretization and underweight time integration. For the wave equation, the leapfrog scheme is popular and second order accurate. It is also symplectic, which means it preserves energy reasonably well over long simulations. The Newmark-beta method is another option, and it gives you control over numerical dissipation by adjusting the beta parameter. Setting beta to one quarter gives the average acceleration method, which is unconditionally stable for linear problems but less accurate than leapfrog. Explicit methods are fast per step but require small time steps for stability. Implicit methods allow larger time steps but require solving a linear system at each step. For large-scale wave problems, explicit methods are usually preferred unless your CFL condition is prohibitively restrictive. There are hybrid approaches like spectral-element methods that combine spectral accuracy within elements with the flexibility of finite elements for meshing.
Boundary Conditions Are Where Most Simulations Break
This is one of those things that seems simple until it is not. A wave hitting a boundary needs to either reflect, transmit, or be absorbed. If you want absorption, you need a perfectly matched layer or a sponge layer. A simple absorbing boundary condition might look like applying a damping term near the boundary, but these are approximate. The reflection coefficient depends on frequency, so a condition that works well for one frequency range will reflect waves at another frequency. I had a simulation where low-frequency waves were reflecting off the boundary and creating interference patterns that completely distorted the solution. The analytical reflection coefficient for my absorbing boundary was wrong because I had derived it for a different geometry. I ended up implementing a simpler sponge layer instead, where I gradually increased a damping coefficient from zero in the interior to a large value near the boundary over a distance of about ten grid points. It was not elegant, but it worked reliably across the frequency range I cared about.
Validation Should Not Be Optional
Before you trust any numerical wave simulation, validate it against an analytical solution if one exists. The 1D wave equation on a finite domain with fixed endpoints has a well-known eigenfunction expansion solution. Compare your numerical eigenfrequencies and mode shapes against the analytical ones. If your method does not converge to the right frequencies as you refine the mesh, something is wrong. I have seen papers publish results from simulations that failed this basic check, and the errors went undetected for months. If you are just starting out, use a library rather than writing everything from scratch. FEniCS handles finite element wave simulations well. deal.II is another robust option, though it has a steeper learning curve. For spectral methods, Dedalus or PDElab are reasonable choices. MATLAB's pdepe is fine for simple 1D problems but will not scale to anything complex. Python with NumPy and SciPy is sufficient for learning and small-scale problems. A typical workflow looks like this: define your domain and boundary conditions, choose a spatial discretization method, set your time stepping scheme, run a convergence study with at least three mesh refinements, compare against an analytical solution, and only then proceed to more complex scenarios. The convergence study takes time, usually a few hours on a modern laptop for a 2D problem, but skipping it is the fastest way to produce unreliable results.
Known Limitations
No single method works for all wave problems. Finite difference methods fail on complex geometries without significant mesh generation effort. Spectral methods struggle with discontinuities and non-periodic boundaries. Finite element methods become computationally expensive for high-frequency problems because you need many elements per wavelength to maintain accuracy, which means very large systems. A rule of thumb is that you need at least ten elements per wavelength for finite elements and about twenty grid points per wavelength for finite differences to keep dispersion errors below one percent. For very high-frequency wave propagation in large domains, neither traditional FEM nor spectral methods are efficient. Boundary element methods or high-order discontinuous Galerkin methods may be more appropriate, but they add significant complexity to the implementation. If your problem involves wave scattering from complex objects at high frequencies, consider whether a ray-based or asymptotic method might be more practical, even though it sacrifices some accuracy.
Common Pitfalls
The most common mistake is using a time step that is close to the stability limit without testing it. This leads to slow growth of high-frequency errors that are hard to detect until the simulation produces obviously wrong results. Another common mistake is ignoring the difference between phase velocity and group velocity in dispersive numerics. On a coarse grid, numerical dispersion causes different frequency components to travel at different speeds, which smears out wave packets over time. If you need to preserve the shape of a wave packet, you will need a finer mesh or a higher-order method. People also tend to underestimate the importance of initial conditions. A discontinuous initial condition in a finite difference method creates high-frequency content that the grid may not resolve. Adding a small amount of smoothing to the initial condition, like a Gaussian convolution, can prevent issues without significantly affecting the physics you care about.
