Explicit time integration

Predicting the evolution of dynamic systems—from the collision of galaxies to the fracture of a bone—is a cornerstone of modern science and engineering. The fundamental challenge lies in translating the laws of physics into a step-by-step computational recipe that can march a system forward in time. Explicit time integration offers a powerful and intuitive solution to this problem, forming the backbone of countless simulations of transient phenomena. This approach is built on a simple premise: if we know everything about a system right now, we can directly calculate what it will look like a moment later.

This article explores the world of explicit time integration, demystifying its core concepts and showcasing its vast reach. We will begin by examining the foundational "Principles and Mechanisms," exploring simple algorithms like the Forward Euler and central difference methods. You will learn how information propagates through a computational grid and understand the universal "speed limit" imposed by the Courant–Friedrichs–Lewy (CFL) condition, a critical concept for ensuring a simulation's stability. We will also uncover the practical challenges these methods face, most notably the crippling problem of "stiffness."

Following this, the section on "Applications and Interdisciplinary Connections" will take you on a tour of the diverse fields where these methods are indispensable. From capturing shock waves in aerospace engineering to simulating galactic evolution in astrophysics, we will see how the same principles apply. This journey will reveal how scientists and engineers overcome the inherent limitations of explicit methods through clever adaptations, highlighting the ongoing innovation at the frontiers of computation, from reduced-order modeling to the integration of machine learning.

At the heart of predicting the future, whether it's the trajectory of a planet or the ripple of a pond, lies a simple and profound idea: if we know everything about a system right now—its state—and we know the laws that govern its change, we can figure out what it will look like a moment later. Explicit time integration is the computational embodiment of this philosophy, a method for stepping through time, one small leap at a time.

Imagine you are tracking a particle. Its state at any time $t$ can be described by its position $\boldsymbol{u}(t)$ and its velocity $\boldsymbol{\dot{u}}(t)$. The fundamental question is, given the state at time $t$, what is the state at a slightly later time, $t + \Delta t$? The most straightforward guess is to assume the velocity is constant over this small interval. This gives us the famous ​​Forward Euler​​ method:

But where does the rate of change, $\boldsymbol{\dot{u}}(t)$, come from? For a mechanical system, this involves acceleration, which is determined by the forces acting on the particle at that instant, according to Newton's second law, $\boldsymbol{F} = m\boldsymbol{a}$. If the forces themselves depend only on the current positions and velocities, we can compute the acceleration, update the velocity, and then update the position. We have everything we need at time $t$ to find the state at $t+\Delta t$. This is the very definition of an ​​explicit​​ method: the future is calculated solely from the present.

For a system with many degrees of freedom, we can write this process in a more powerful matrix form. If we let $\mathbf{u}^k$ be a vector representing the entire state of our system at time step $k$, the explicit update can be seen as the action of a special matrix, the ​​amplification matrix​​ $G$, that propels the system forward in time:

Simulating the system for many steps is then equivalent to repeatedly applying this matrix: $\mathbf{u}^{k+N} = G^N \mathbf{u}^{k}$. The entire fate of our numerical world is encoded in the properties of this single matrix, $G$. For the simulation to be stable and not explode into infinity, the eigenvalues of $G$ must be carefully controlled, a point we will return to with great consequence.

While Forward Euler is intuitive, it has certain numerical quirks. A more popular and robust choice for simulating physical phenomena like waves and vibrations is the ​​explicit central difference​​ method. It gets its name from how it approximates derivatives, but a more intuitive picture is that of a game of leapfrog.

To find the position at the next time step, $\boldsymbol{u}^{n+1}$, the scheme looks at the current position, $\boldsymbol{u}^{n}$, and also takes a glance backward at the previous position, $\boldsymbol{u}^{n-1}$. The core of the method is a beautifully symmetric update rule derived directly from a Taylor expansion of time:

Here, $\boldsymbol{a}^n$ is the acceleration at the current time step $n$. Since the acceleration is computed from the forces at time $n$, which in turn depend only on the state at time $n$ (and perhaps previous steps), this scheme is fully explicit. There is no need to solve complex equations involving the unknown future state $\boldsymbol{u}^{n+1}$. This scheme is not only second-order accurate, a step up from Forward Euler, but it also has excellent energy conservation properties over long simulations, making it a favorite among physicists and engineers for modeling dynamic events.

Physical reality is a continuum, but our computers can only handle a finite list of numbers. So, we discretize space, replacing the continuum with a grid of points, often called a mesh. The governing partial differential equations (PDEs), which describe how properties like pressure or displacement vary continuously in space and time, become a large, coupled set of ordinary differential equations (ODEs)—one for each point on the grid.

Consider a simple elastic bar, modeled as a line of nodes connected by springs. The force on an interior node $j$ depends on its displacement relative to its neighbors, nodes $j-1$ and $j+1$. When we compute the acceleration at node $j$, we only need to look at this immediate neighborhood. This local recipe for computing forces is called a ​​stencil​​.

This locality has a profound consequence. In a single explicit time step $\Delta t$, information from node $j$ can only influence its immediate neighbors defined by the stencil. Information propagates across the grid one layer of nodes at a time per step, like a line of dominoes. The region of the grid at time $t^n$ that influences the solution at a single point $(x_j, t^{n+1})$ is called the ​​numerical domain of dependence​​. After one time step, this domain is just the stencil; after $N$ steps, it's a region roughly $N$ grid points wide.

Now we arrive at one of the most fundamental principles in numerical simulation. The physical world has its own rules about how fast information can travel. For the acoustic wave equation, the speed limit is the speed of sound, $c$. For an elastic bar, it's the elastic wave speed $c = \sqrt{E/\rho}$. The true solution of the PDE at a point $(x, t+\Delta t)$ is determined by the data within a specific region at the earlier time $t$. This region is the ​​physical domain of dependence​​.

The ​​Courant–Friedrichs–Lewy (CFL) condition​​ is a powerful statement of causality: for a simulation to be valid, the numerical domain of dependence must be large enough to contain the physical domain of dependence. In other words, the numerical algorithm must have access to all the physical information required to determine the answer.

If a physical wave can travel from node $j-1$ to node $j+1$ in a single time step $\Delta t$, but our numerical stencil only connects immediate neighbors, the scheme cannot possibly capture the correct physics. The physical signal has outrun the numerical dominoes. This leads to a stability constraint on the size of the time step. For a simple 1D wave on a grid with spacing $h$, the condition is:

The dimensionless quantity $\mathrm{CFL} = c \Delta t / h$ is the Courant number. The stability limit depends on the details of the scheme and the dimension $d$. For a standard finite-difference scheme for the wave equation, this limit is $\Delta t \le h / (c\sqrt{d})$. For the heat diffusion equation, the constraint is even more severe, scaling with the square of the grid spacing: $\Delta t \le \Delta x^2 / (2\alpha)$. This means that if you refine your mesh by a factor of 10 to get more spatial detail, you must take 10 times more steps for a wave problem, but a staggering 100 times more steps for a diffusion problem!

The CFL condition reveals a trade-off: higher spatial resolution demands smaller time steps. Despite this, the supreme advantage of explicit methods is their computational efficiency per step. Because the update for each node's state is a local calculation based on known values, there is no need to solve a massive, interconnected system of linear equations. This makes the workload "embarrassingly parallel"—we can assign different regions of the mesh to different processors, and they can all compute their updates simultaneously, only needing to communicate boundary data between steps.

This efficiency can be threatened by the formalism of some discretization methods, like the Finite Element Method (FEM). In FEM, the equations of motion often take the form $\mathbf{M} \ddot{\mathbf{d}} = \mathbf{F}$, where $\mathbf{M}$ is the ​​mass matrix​​. To find the accelerations $\ddot{\mathbf{d}}$, we must compute $\mathbf{M}^{-1} \mathbf{F}$. If $\mathbf{M}$ is a large, fully-coupled matrix (a so-called ​​consistent mass matrix​​), this would require solving a linear system at every step, destroying the method's efficiency.

To circumvent this, practitioners use a wonderfully pragmatic trick called ​​mass lumping​​. The idea is to approximate the consistent mass matrix with a diagonal one. Inverting a diagonal matrix is trivial: it is simply the element-wise division of the force vector by the nodal masses. This single trick can change the computational complexity of each time step from $\mathcal{O}(n^{1.5})$ or worse for solving a sparse system to a blissful $\mathcal{O}(n)$ for simple division. For a problem with a million nodes, this can mean a speedup of a thousand times or more per step.

Interestingly, some advanced methods like the ​​Discontinuous Galerkin (DG)​​ method are naturally suited for explicit integration. By design, DG methods produce a ​​block-diagonal mass matrix​​, where each block is a small matrix corresponding to a single element. Inverting these small, independent blocks is computationally cheap and perfectly parallel, preserving the key advantages of the explicit approach while offering high-order accuracy.

With their speed and simplicity, explicit methods seem almost too good to be true. And indeed, they have an Achilles' heel. The first hint of trouble is the realization that the CFL condition, while necessary, is not always sufficient for stability. It is possible to devise a scheme—like the forward-time, centered-space method—that respects the domain of dependence principle yet is unconditionally unstable. Its amplification factor has a magnitude greater than one for any non-zero time step, meaning errors are always amplified, leading to an inevitable explosion. The precise details of the discretization matter.

The truly formidable barrier, however, is a property known as ​​stiffness​​. A system is stiff when it involves processes occurring on vastly different time scales. Consider a chemical reaction where some species react in microseconds, while the overall mixture evolves over minutes, or a nuclear reactor where some isotopes decay in microseconds while others persist for thousands of years.

The stability of an explicit method is always dictated by the fastest phenomenon in the system. The time step $\Delta t$ must be small enough to resolve the microsecond-scale event. However, to capture the slow evolution of the entire system, the total simulation must run for minutes, or days, or years. The total number of steps required becomes the ratio of the slow physical timescale to the fast numerical timescale dictated by stability. This is the ​​stiffness ratio​​, $S = \tau_{slow} / \tau_{fast}$.

For a realistic nuclear reactor simulation, this ratio can be of the order of $10^{14}$. This means you would need to take on the order of $10^{10}$ tiny, microsecond-long time steps just to simulate one day of operation. The computational cost is astronomical. The explicit method is held hostage by the fastest, often fleeting, transient, forced to crawl at a snail's pace even long after that transient has vanished. This is the wall where the beautiful simplicity of explicit methods breaks down. To overcome stiffness, we must relinquish the comfort of calculating the future from the present and venture into the world of ​​implicit methods​​, where the future and present are solved for together.

Now that we have the basic machine in hand—this beautifully simple idea of taking small, explicit steps forward in time—let's take it for a spin. Where does this notion of marching cautiously through a simulation actually take us? The answer, it turns out, is almost everywhere. From the catastrophic failure of a bone under impact to the majestic dance of galaxies, from the flow of a river to the inner workings of an artificial intelligence, the same fundamental principles of explicit time integration appear again and again. The journey is not just a tour of different sciences, but a deeper look into the very nature of our numerical world, a world governed by its own universal speed limit.

Explicit methods are at their most natural and powerful when we are trying to understand phenomena that are, by their very nature, fast, violent, and brief. Think of a shock wave, an explosion, or a high-speed collision. The physics itself is a cascade of events happening on very short timescales, and our numerical method is perfectly suited to follow along, step-by-tiny-step.

Consider the propagation of a shock wave through the air, a problem central to aerospace engineering. The governing physics is captured by the Euler equations, a set of laws for the conservation of mass, momentum, and energy. When we simulate this using an explicit scheme, we find that the stability of our calculation is dictated by the famous Courant-Friedrichs-Lewy (CFL) condition. This condition is not some arbitrary numerical rule; it is the embodiment of a deep physical principle. It tells us that our numerical time step $\Delta t$ must be small enough that information—in this case, traveling at the speed of sound $c$ relative to the moving fluid $u$—does not jump over a whole computational cell in a single step. In a very real sense, the simulation is not allowed to get ahead of the physics. The speed of sound becomes the speed limit for our numerical universe.

This same principle applies with equal force to solids. Imagine a forensic investigation trying to understand how a skull fractures under a blunt-force impact. By modeling the bone using the Finite Element Method, engineers can simulate the event. Here, the "information" is a stress wave traveling through the material. Its speed, determined by the bone's stiffness ($E$) and density ($\rho$) as $c = \sqrt{E/\rho}$, now sets the speed limit. To capture the fracture process, which happens over milliseconds, the simulation must take incredibly small time steps, often on the order of nanoseconds, to respect the time it takes for the stress wave to cross the smallest element in their computer model. Sometimes, to make a simulation feasible, an analyst might artificially increase the density of the material in the model. This "mass scaling" slows down the computed stress waves, allowing for a larger, more economical time step, but at the cost of slightly distorting the true timing of events—a necessary, but carefully considered, compromise.

The real world, of course, rarely presents us with just a fluid or just a solid. What happens when a fluid shock wave slams into a solid structure? This is a classic fluid-structure interaction problem. Here, we see the first signs of trouble with our simple explicit approach. In a simulation of a shock wave hitting a plate, the fluid side might demand a time step of nanoseconds to resolve the shock, while the plate itself might vibrate with a characteristic period of microseconds. In a simple "partitioned" approach where both simulations must march in lockstep, the entire calculation is held hostage by the fluid's tiny time step. Simulating a few milliseconds of the plate's vibration could take hundreds of thousands of steps, a computationally Herculean task. This mismatch in timescales is a new kind of beast, a phenomenon called "stiffness," which will become a recurring theme in our journey.

While explicit methods are beautiful in their simplicity, they live under a harsh dictatorship. The stability of the entire simulation, no matter how vast, is dictated by the most restrictive condition anywhere in the domain. This can lead to what we might call the "tyranny of the smallest and the fastest."

One form of this tyranny is geometric. Imagine we are simulating the propagation of a fracture through rock using an advanced technique like the Extended Finite Element Method (XFEM). This method allows a crack to cut through the elements of our computational mesh. In doing so, it can create tiny, sliver-like sub-regions for calculation. While these "slivers" may be geometrically insignificant to the overall picture, they are not insignificant to our explicit integrator. The stability condition, remember, depends on the time it takes a wave to cross the smallest characteristic length in the model. A single, tiny sliver, perhaps a fraction of a millimeter across, can force the entire simulation of a kilometer-scale geological fault to take pitifully small time steps, grinding the calculation to a halt. The smallest part of the model holds the entire simulation hostage.

An even more profound form of tyranny is physical. Many systems in nature involve processes that occur at vastly different speeds. Consider a simple model of nutrients in the ocean. The chemical conversion of ammonium to nitrate (nitrification) can be extremely fast, with a characteristic timescale of less than a day. At the same time, the slow physical process of ventilation—mixing with the deep ocean—occurs on a timescale of decades. If we want to simulate how the ocean's nutrient balance evolves over a century, we are interested in the slow process. But if we use an explicit method, we find that our time step is limited by the stability of the fast chemical reaction. We are forced to take steps of minutes or hours, even though the changes we care about happen over years. This is the classic definition of a "stiff" system. The stability is governed by the fastest, often uninteresting, timescale, while the accuracy is determined by the slow timescale of interest. The computational cost becomes astronomical.

This challenge is not unique to biochemistry. It is everywhere.

• In a viscoelastic material, the time step might be limited by the minimum of two constraints: one from the speed of elastic waves ($h/c$) and another from the material's internal relaxation time ($\tau$).
• In a thermomechanical problem, the limit could be the minimum of a thermal diffusion timescale ($h^2/\alpha$) and a mechanical viscosity timescale ($\eta/E$).
• In a complex environmental simulation of a flood using Smoothed Particle Hydrodynamics (SPH), the time step is a three-way battle between the wave speed (CFL condition), the water's viscosity, and the acceleration due to gravity.

In all these cases, the overall time step is given by:

The fastest process wins, and our computational budget loses. When faced with this, and if we cannot switch to a different (implicit) method, we sometimes resort to other tricks. For instance, in impact simulations, we might introduce a small amount of non-physical "numerical damping" to kill off spurious, high-frequency oscillations that are excited by the model but are too expensive to resolve properly. This is a delicate balancing act: we add just enough damping to stabilize the calculation without corrupting the slower, physically important response we are trying to capture.

The story doesn't end with tyranny and compromise. Scientists and engineers are a clever bunch, and they have developed ingenious ways to adapt and overcome the limitations of explicit integration.

A wonderfully intuitive idea is that not every part of a simulation needs to run at the same speed. Consider simulating the evolution of a galaxy, a collection of a million stars held together by gravity. In the sparse outer regions of the galaxy, stars move slowly and their gravitational environment changes gradually. They can be simulated with large time steps. But in the dense central core, stars are whipping around each other, and the forces are changing violently from moment to moment. Here, we need tiny time steps. An adaptive timestep scheme gives each particle its own individual clock. A common criterion sets a particle's time step $\Delta t_i$ to be proportional to the ratio of its acceleration's magnitude to the rate of change of its acceleration (the "jerk"): $\Delta t_i \propto |\mathbf{a}_i|/|\dot{\mathbf{a}}_i|$. This is a measure of the characteristic time over which the force field is changing. It's a beautifully simple idea that allows the computer to focus its effort where the action is, making such simulations possible.

Another powerful modern approach is to not simulate the full, complex problem at all. Reduced-Order Models (ROMs) are based on the idea that the behavior of many complex systems, even if described by thousands or millions of variables, might actually evolve within a much smaller, simpler subspace. By projecting the governing equations onto this small subspace, we create a much smaller model that is far cheaper to solve. What's remarkable is that if this projection is done in a way that respects the underlying structure of the physics (a so-called Galerkin projection), the resulting ROM is guaranteed to be no less stable than the original full model. In fact, since the ROM inherently filters out the high-frequency components of the original system that were so restrictive, its own stability limit is often much more generous, allowing for significantly larger time steps. We can also see the elegance of this approach in its warnings: if one uses a method that does not preserve the physical structure, the resulting ROM can become unstable and produce complete nonsense.

Perhaps the most exciting frontier is the intersection of simulation with machine learning. What happens when we don't have a perfect physical law for a material, but instead have a model learned from experimental data by a neural network? Suppose a neural network, $\mathcal{N}$, learns to describe part of the evolution of a material's internal state. Can we still use our explicit integrator? Yes! And remarkably, the stability analysis follows the exact same logic. The stability of the time step is now found to depend on the properties of the neural network itself—specifically, on its Lipschitz constant $L$, which is a measure of the maximum "steepness" of the learned function. The final stability limit becomes a function of both the known physics (a relaxation time $\tau$) and the property of the AI model, yielding a limit like: $\Delta t_{\max} = \frac{2\tau}{1 + \tau L}$ This stunning result shows that the fundamental principles of numerical stability are universal. They apply just as well to the differential equations we derive from physical principles as they do to the functions learned by our most advanced artificial intelligence models.

From the speed of sound in air to the Lipschitz constant of a neural network, the simple requirement of stability for an explicit time integrator has taken us on a grand tour of science. It has revealed itself not as a mere numerical constraint, but as a deep principle connecting computation, physics, and the very flow of information. It forces us to think deeply about the multiple timescales of nature and rewards us with powerful tools to explore worlds that would otherwise remain forever beyond our sight.