Getting Time-Stepping Right for Wave Equations
Transient wave problems are one of those things where every textbook looks clean until you actually code it and half your solutions blow up. The issue isn't usually the math itself. It's picking the right method, matching your time step to your spatial discretization, and then dealing with boundary conditions that don't reflect back into your domain like they should. I spent about three years trying to get a staggered-grid finite-difference scheme to stay stable in 3D before I stopped fighting the CFL condition and started respecting it properly.
Order Numerical Methods For Transient Wave Equations
This is really about two decisions: what order your spatial discretization is, and what order your time integrator is. For the wave equation, these two have to play nice together. A fourth-order spatial scheme paired with a first-order forward Euler time integrator will not save you. It just makes the instability prettier before it crashes. The standard approach that actually works in practice is a second-order central difference in time with a second or fourth-order spatial stencil. That gives you a method that is second-order accurate in time and whatever order you built into your space derivative. Explicit schemes are the default because the wave equation is hyperbolic and the physics propagate along characteristics. Implicit methods like backward Euler or Crank-Nicolson are possible but introduce artificial numerical damping that kills high-frequency content, which is exactly what you're trying to track in a transient wave simulation. Here is what that looks like on a basic 1D form of the problem. If you write the wave equation as u_tt = c^2 u_xx and discretize it, your update at time step n+1 depends on the current step n and the previous step n-1. The scheme is:
u_i^(n+1) = 2u_i^n - u_i^(n-1) + (c*dt/dx)^2 * (u_{i+1}^n - 2u_i^n + u_{i-1}^n) The ratio c*dt/dx is your Courant number. If that number is greater than 1, your solution will diverge regardless of anything else you do. This is not a suggestion. It is the von Neumann stability analysis telling you exactly where the line is.
Get the Full Details
Higher-Order Spatial Discretizations
Once you have the second-order scheme working, you will naturally want to improve accuracy without slashing your time step. You can swap the central difference for a fourth-order stencil, which uses two points on either side instead of one. The spatial update becomes: u_xx (u_{i+2} + 16u_{i+1} 30u_i + 16u_{i1} u_{i2}) / (12dx^2) This gives you much better dispersion properties. The second-order scheme has significant numerical dispersion even at CFL numbers well below 1. Waves of different frequencies travel at different speeds in the discretization, which means a sharp pulse spreads out artificially over time. The fourth-order stencil reduces that error substantially. You are trading a wider stencil, more memory, and slightly more complicated boundary handling for accuracy that actually matters when you run long simulations.
I ran into a specific problem with this on a 3D acoustic propagation model a few years ago. I had a domain with irregular geometry and I was using a fourth-order spatial scheme with a second-order leapfrog time integrator. Near the boundaries, where I had to truncate the stencil, I was getting reflections that contaminated the interior solution by about 15 percent. Standard absorbing boundary conditions weren't killing them. What actually worked was switching to a perfectly matched layer approach with a graded conductivity profile inside the PML region. I made the layer about ten cells thick and ramped the absorption coefficient quadratically from zero at the interior boundary to a maximum at the outer edge. That dropped spurious reflections to below 0.1 percent and the solution quality improved immediately.
Time Integration Options Beyond Leapfrog
Leapfrog is fine for pure wave propagation with no source terms or damping. The moment you add anything that breaks the symmetry of the problem, you need a different time integrator. Source terms, variable material properties, and nonlinear effects all show up in real simulations. A reasonably robust alternative is a third-order Runge-Kutta method, specifically the Shu-Osher variant. It handles source terms cleanly and does not require storing the previous time step like leapfrog does. The tradeoff is that it needs three sub-steps per time advance, so it is computationally more expensive per step. In practice though, you can sometimes take a larger time step with RK3 than you could with leapfrog if your CFL constraint is tight, and the overall wall time can come out comparable. For fourth-order spatial accuracy, people often pair it with a four-stage Runge-Kutta method. This keeps the temporal error from becoming the dominant source of inaccuracy. If your spatial error is O(dx^4) and your time error is O(dt^2), you are burning half your accuracy budget on time integration. That is a common mistake I see in student code and in some research implementations too.
Dispersion and Stability in Practice
The thing about wave equations that no one warns you about early enough is how quickly numerical dispersion ruins results. Even when your scheme is formally stable, the phase error accumulates. A wave that should arrive at a point at t=0.5 might arrive at t=0.52 with the fourth-order scheme and t=0.6 with the second-order scheme. Over many wavelengths of propagation, this difference becomes the dominant error. A useful rule of thumb is to resolve each wavelength with at least ten to twelve grid points for second-order schemes and five to six for fourth-order. That is significantly more than the minimum required for stability. Stability just means the solution does not explode. Accuracy is a separate question. I also learned the hard way that variable wave speed breaks the simple CFL analysis. When c varies spatially, you have to use the maximum value of c in your domain for the stability limit, not the average. A friend of mine was simulating seismic wave propagation through layered earth models and kept getting instabilities that he could not explain. He had computed his time step based on the average velocity. The actual limiting factor was a thin high-velocity layer that occupied less than two percent of his domain. Switching to the maximum velocity for his CFL calculation fixed it immediately.
Boundary Conditions and Their Problems
Boundary conditions are where most wave simulations go wrong. Dirichlet and Neumann conditions are straightforward at first but they create reflections at domain edges that are physically meaningless. If your domain is finite, waves hit the boundary and bounce back unless you do something about it. Beyond PML layers, there are first-order absorbing boundary conditions which approximate the behavior of outgoing waves. They work reasonably well for waves arriving near normal incidence but perform poorly for grazing angles. Second-order Engquist-Majda conditions are better across a wider range of angles but add complexity to your implementation. For most practical engineering work, PML is the standard and it is worth the extra effort to implement correctly. The catch with PML is that it only works properly if your coordinate stretching is smooth. Discontinuous transitions between the physical domain and the PML region generate reflections that are worse than what you would get from a simple truncated domain. I wasted about two weeks debugging reflections that turned out to be caused by a step discontinuity in the absorption profile at the PML interface. Making it continuous fixed everything in an afternoon.
What This Method Cannot Do Well
Explicit finite-difference methods for the wave equation struggle with very fine spatial features relative to your domain size. If you need to resolve sub-millimeter details in a meter-scale domain, your grid becomes unmanageable and your time step shrinks to something impractical. There is no workaround within the explicit framework except to refine, and refinement in all three dimensions scales poorly. Another limitation is that these methods assume you know the wave speed field in advance. If you are doing coupled problems where the medium properties change dynamically based on the wave field itself, you enter nonlinear territory where the standard approaches need modification and additional stability analysis. Some people try to use these methods for nonlinear problems without rechecking the stability conditions, and it rarely goes well. If you are dealing with complex geometries, explicit finite differences on structured grids become awkward. Body-fitted grids or unstructured meshes require more sophisticated treatments of the stencil near boundaries. In those cases, finite element methods with explicit time integration or discontinuous Galerkin methods are often better choices, though they come with their own implementation overhead.
A Practical Starting Point
If you are building a solver from scratch, start with the 1D second-order scheme. Get it stable, verify it against an analytic solution like a Gaussian pulse propagating in a lossless medium, and confirm that your convergence rate is actually second-order in both space and time. Then move to 2D. Then add the fourth-order spatial stencil. Then implement absorbing boundaries. Each step introduces new failure modes that you need to diagnose separately. The common error at each stage is assuming that what worked in one dimension carries over directly. Boundary conditions in 2D and 3D are not trivial extensions. The PML implementation in 3D requires splitting the wave equation into components and applying the absorption independently in each direction, which is straightforward but easy to get wrong if you are writing it for the first time. I found that writing a 2D version first and testing it against the known analytical solution for a point source in a homogeneous medium caught most of the bugs before I ever attempted 3D. There is no universal implementation you can just download and trust for this kind of problem. The stability properties depend on your specific grid, your wave speed distribution, and your boundary treatment. Verify everything against a case where you know the answer before you trust it on a case where you do not.
