Abstract

The aim of this research is to examine Courant-Friedrichs-Lewy condition for the analysis of convergence and stability of explicit forward time central space scheme for a three-dimensional wave equation. The wave equation, which models physical phenomena such as sound and electromagnetic wave propagation, is discretized using finite difference methods in both time and space. A central difference scheme is implemented to approximate the second-order derivatives across spatial and temporal domains. The CFL condition is derived as a criterion to ensure numerical stability and is shown to depend on the wave speed and spatial grid resolution. The explicit update scheme is constructed and analyzed under uniform grid spacing. Through von Neumann stability analysis, the amplification factor is expressed using Fourier modes and Euler’s identity. The characteristic equation for the scheme is derived, and its roots are examined to determine the conditions for numerical stability. The CFL number $\lambda =\frac{c\text{Δ}t}{h}$ is introduced and bounded to prevent error magnification. The analysis confirms that for stability, the time step must satisfy $\text{Δ}t\le \frac{h}{c\sqrt{3}}$ . Convergence is discussed in the context of satisfying both consistency and stability criteria. The initial and boundary conditions necessary for realistic modeling are incorporated. This work validates that adherence to the CFL condition is essential for reliable and accurate simulation of three-dimensional wave propagation using explicit finite difference methods.

Keywords

Courant-Friedrichs-Lewy Condition: Convergence, Stability, Explicit FTCS, Wave Equation

Share and Cite:

1. Introduction

Accurate numerical solutions of partial differential equations are critical for simulating wave propagation. The explicit Forward-Time Central-Space (FTCS) method offers simplicity but requires careful stability and convergence analysis. The Courant-Friedrichs-Lewy (CFL) condition establishes the maximum time step relative to spatial discretization to ensure stability. In three-dimensional wave equations, the multi-directional propagation increases computational challenges, making CFL analysis essential. This study investigates the CFL condition for three-dimensional FTCS schemes, providing guidelines for stable and convergent simulations.

Traditional finite-difference methods (FDMs) find it hard approximating acoustic wave propagation when the CFL number goes beyond 0.707 in 2D or 0.577 in 3D for equally spaced grids. This limits how large the time step can be. To address this, researchers have developed a variable-length temporal and spatial operator approach that allows wave modeling beyond these limits without losing accuracy. The idea is to make the temporal operators slightly longer than the spatial operators in high-velocity areas, which helps maintain stability at larger time steps. At the same time, both operator lengths are adjusted according to local velocity to keep the results accurate. With this method, the CFL number can be pushed up to 1.25 in 2D and 1.0 in 3D for high-velocity contrast cases. Tests in both simple and complex media show that the method works well [1].

[2] analyzed the transient diffusion equation in one dimension with diffusion coefficients that vary in both space and time. These types of transport equations, which can be derived from the Fokker-Planck equation, are essential for understanding diffusion mechanisms in general, such as those occurring in carbon nanotubes. Using the classical self-similar Ansatz, the authors obtained new, nontrivial analytical solutions, which they then reproduced using 16 explicit numerical time-integration methods-11 of them recent and unconditionally stable. The findings showed that certain algorithms, such as the leapfrog-hopscotch method, proved highly efficient and, in some cases, outperformed the standard FTCS method.

[3] examined the stability of three-dimensional numerical evolutions of the Einstein equations, comparing the standard ADM formulation with variations of a conformal-traceless (CT) formulation that separates the conformal and traceless parts of the system. The authors developed a CT implementation with improved stability for evolving both weak and strong gravitational fields, in vacuum and in spacetimes coupled to matter sources. Their tests included weak and strong gravitational wave packets, black holes, boson stars, and neutron stars. They identified the conditions under which the CT approach produced better results than ADM in 3D simulations. Overall, their CT implementation yielded more stable long-term evolutions in all cases studied, although it was less accurate in the short term for the range of resolutions used.

[4]-[6] addressed the often-overlooked effect of tidal movement on the shoreline by incorporating the moving boundary into a shallow water model through a coordinate transformation and a Lax-Friedrichs time-explicit scheme. Applied to the Ameland inlet system, the model derived from the 3D Navier-Stokes equations produced realistic results, capturing both steady conditions at the seaward side and nonlinear effects near the landward side. Sensitivity tests across wave amplitude, water depth, basin length, and resistance confirmed the model’s stability and accuracy.

With advances in distributed computing, interest has grown in explicit time-integration methods that offer greater stability, avoiding the parallelization challenges of implicit schemes [7]. This work introduced a weighted difference scheme combining delayed and conventional explicit methods which achieved a higher stability limit of 1.5 and eliminated the checkerboard instability of the delayed approach.

This work is an extension of [8] who examined the time-space domain explicit FDM that numerically solves the wave equation by approximating its spatial and temporal derivatives but often faces stability issues. The results showed that the maximum stable CFL number depends on the peak value of the spatial FD dispersion relation. While conventional methods determine spatial FD coefficients by matching the dispersion relation within a specific wavenumber range, indicating that outside this range, the dispersion and the CFL number are uncontrolled. their work involed a 2D wave propagation but we have considered a 3D case for better generalization and also set conditions for stability and convergence for FTCS scheme as these are interconnected.

It is important to clarify the terminology used in this work. The classical Forward-Time Central-Space (FTCS) scheme refers to a first-order accurate Euler forward discretization in time combined with second-order central differences in space. However, the wave equation ${\partial }^{2}u/\partial t^2=c^2{\nabla }^{2}u$ contains a second-order time derivative, which necessitates a different treatment.

2. Mathematical Formulation

The wave equation in three spatial dimensions models the propagation of waves in physical systems such as acoustics, electromagnetics, and elastic media. The classical 3D wave equation is given by

$\frac{{\partial }^{2}u}{\partial t^2}=c^2{\nabla }^{2}u,$ (1)

where $u\left(\zeta ,t\right)$ is the wave displacement, $\zeta =\left(x,y,z\right)\in \Omega \subset {ℝ}^3$ is the spatial coordinate vector, $t\ge 0$ is time, and $c$ is the constant wave speed [5] [9]. As discussed in [10]-[12], the Laplacian operator in three dimensions is

${\nabla }^{2}u=\frac{{\partial }^{2}u}{\partial x^2}+\frac{{\partial }^{2}u}{\partial y^2}+\frac{{\partial }^{2}u}{\partial z^2}.$ (2)

$\frac{{\partial }^{2}u}{\partial t^2}=c^2\left(\frac{{\partial }^{2}u}{\partial x^2}+\frac{{\partial }^{2}u}{\partial y^2}+\frac{{\partial }^{2}u}{\partial z^2}\right)$ (3)

For one dimensional wave equation describing vibrations on a string given by

$\frac{{\partial }^{2}u}{\partial t^2}=c^2\frac{{\partial }^{2}u}{\partial {\xi }^2}$ (4)

has a general solution consists of two traveling waves:

$u\left(\xi ,t\right)=f\left(\xi -ct\right)+g\left(\xi +ct\right)$ (5)

where:

$f\left(\xi -ct\right)$ is a right-moving wave, $g\left(\xi +ct\right)$ is a left-moving wave, $c$ is the wave propagation speed, $f$ and $g$ are the functions describing displacement and velocity of a wave respectively.

However, for a 3D wave function, the propagation is in all three spatial directions, therefore, the general solution is given by:

$u\left(x,y,z,t\right)=f\left(k_{x}x+k_{y}y+k_{z}z-\omega t\right),$ (6)

and

$c=\frac{\omega }{\left|k\right|},\text{where}\left|k\right|=\sqrt{k_x^2+k_y^2+k_z^2}.$ (7)

which defines the magnitude of the velocity at which the wavefronts move through space. where $k_x,k_y,k_z$ are spatial components of the wave vector $k=\left(k_x,k_y,k_z\right)$ , which determine the direction of wave propagation in 3D space, $\omega$ is the angular frequency, which determines how fast the wave oscillates in time, $k_{x}x+k_{y}y+k_{z}z-\omega t$ is the wave phase, which moves with constant speed.

2.1. Discretization Scheme

The numerical method employed in this study discretizes both temporal and spatial derivatives using second-order accurate central difference approximations. Unlike the classical Forward-Time Central-Space (FTCS) scheme—which uses a first-order forward Euler approximation for $\partial u/\partial t$ —the wave equation requires approximation of the second-order time derivative ${\partial }^{2}u/\partial t^2$ . This is accomplished by three-level central difference:

${\frac{{\partial }^{2}u}{\partial t^2}|}_{t^n}\approx \frac{u^{n+1}-2u^n+u^{n-1}}{\text{Δ}{t}^2},$

which is second-order accurate in time. This scheme is more precisely classified as a central-time central-space (CTCS) or leapfrog scheme. The update formula requires solution values from two previous time levels ($n$ and $n-1$ ), making it a two-step explicit method. We discretize the cubic spatial domain $\Omega =\left[a_x,b_x\right]\times \left[a_y,b_y\right]\times \left[a_z,b_z\right]$ uniformly:

$x_i=a_x+i\Delta x,y_j=a_y+j\Delta y,z_k=a_z+k\Delta z,$

for integers

$i=0,1,2,\cdots ,N_x$ ,

$j=0,1,2,\cdots ,N_y$ ,

and

$k=0,1,2,\cdots ,N_z$ ,

where

$\Delta x=\frac{b_x-a_x}{N_x},\Delta y=\frac{b_y-a_y}{N_y},\Delta z=\frac{b_z-a_z}{N_z}.$

Time is discretized as

$t^n=n\text{Δ}t,n=0,1,2,\cdots ,N_t$

with time step $\text{Δ}t$ .

The numerical solution approximating $u\left(x_i,y_j,z_k,t^n\right)$ is denoted by $u_{i,j,k}^n$ .

We define a uniform grid over the domain:

$x_i=i\Delta x,y_j=j\Delta y,z_k=k\Delta z,t^n=n\Delta t,$

and let $u_{i,j,k}^n\approx u\left(x_i,y_j,z_k,t^n\right)$ . The following are the central differences in both time and space according to [9] [12]-[14]

$\frac{{\partial }^{2}u}{\partial t^2}\approx \frac{u_{i,j,k}^{n+1}-2u_{i,j,k}^n+u_{i,j,k}^{n-1}}{\text{Δ}{t}^2},$ (8)

$\frac{{\partial }^{2}u}{\partial x^2}\approx \frac{u_{i+1,j,k}^n-2u_{i,j,k}^n+u_{i-1,j,k}^n}{\text{Δ}{x}^2},$ (9)

$\frac{{\partial }^{2}u}{\partial y^2}\approx \frac{u_{i,j+1,k}^n-2u_{i,j,k}^n+u_{i,j-1,k}^n}{\text{Δ}{y}^2},$ (10)

$\frac{{\partial }^{2}u}{\partial z^2}\approx \frac{u_{i,j,k+1}^n-2u_{i,j,k}^n+u_{i,j,k-1}^n}{\text{Δ}{z}^2}.$ (11)

Substituting these into Equation (3), the explicit finite difference updated scheme is:

$\frac{u_{i,j,k}^{n+1}-2u_{i,j,k}^n+u_{i,j,k}^{n-1}}{{\left(\Delta t\right)}^2}$ (12)

$=c^2(\frac{u_{i+1,j,k}^n-2u_{i,j,k}^n+u_{i-1,j,k}^n}{{\left(\Delta x\right)}^2}+\frac{u_{i,j+1,k}^n-2u_{i,j,k}^n+u_{i,j-1,k}^n}{{\left(\Delta y\right)}^2}$ (13)

$+\frac{u_{i,j,k+1}^n-2u_{i,j,k}^n+u_{i,j,k-1}^n}{{\left(\Delta z\right)}^2}).$ (14)

$\begin{array}{c}{u}_{i,j,k}^{n+1}=2u_{i,j,k}^n-u_{i,j,k}^{n-1}+c^2\Delta t^2(\frac{u_{i+1,j,k}^n-2u_{i,j,k}^n+u_{i-1,j,k}^n}{\Delta x^2}\\ +\frac{u_{i,j+1,k}^n-2u_{i,j,k}^n+u_{i,j-1,k}^n}{\Delta y^2}+\frac{u_{i,j,k+1}^n-2u_{i,j,k}^n+u_{i,j,k-1}^n}{\Delta z^2}).\end{array}$ (15)

CFL condition for the stability of explicit finite difference schemes applied to the three-dimensional wave equation and uniform grid, the time step $\text{Δ}t$ must satisfy:

$c^2\text{Δ}{t}^2\left(\frac{1}{\text{Δ}{x}^2}+\frac{1}{\text{Δ}{y}^2}+\frac{1}{\text{Δ}{z}^2}\right)\le 1,$ (16)

$c\text{Δ}t\le \frac{1}{\sqrt{\frac{1}{\text{Δ}{x}^2}+\frac{1}{\text{Δ}{y}^2}+\frac{1}{\text{Δ}{z}^2}}}$

where $c$ = wave propagation speed, $\text{Δ}t$ = time step size and $\text{Δ}x,\text{Δ}y,\text{Δ}z$ = spatial grid spacing in $x,y,z$ respectively. Using a uniform grid spacing

$\text{Δ}x=\text{Δ}y=\text{Δ}z=h,$

the CFL condition reduces to

$\text{Δ}t\le \frac{h}{\sqrt{3}c}.$

Let the CFL number be:

$\lambda =\frac{c\text{Δ}t}{h}$

According to [15], the finite difference formula become:

$\frac{{\partial }^{2}u}{\partial t^2}\approx \frac{u_{i,j,k}^{n+1}-2u_{i,j,k}^n+u_{i,j,k}^{n-1}}{h^2},$ (17)

$\frac{{\partial }^{2}u}{\partial x^2}\approx \frac{u_{i+1,j,k}^n-2u_{i,j,k}^n+u_{i-1,j,k}^n}{h^2},$ (18)

$\frac{{\partial }^{2}u}{\partial y^2}\approx \frac{u_{i,j+1,k}^n-2u_{i,j,k}^n+u_{i,j-1,k}^n}{h^2},$ (19)

$\frac{{\partial }^{2}u}{\partial z^2}\approx \frac{u_{i,j,k+1}^n-2u_{i,j,k}^n+u_{i,j,k-1}^n}{h^2}.$ (20)

Substituting in (13)

$\begin{array}{c}{u}_{i,j,k}^{n+1}=2u_{i,j,k}^n-u_{i,j,k}^{n-1}+{\lambda }^2(u_{i+1,j,k}^n-2u_{i,j,k}^n+u_{i-1,j,k}^n+u_{i,j+1,k}^n\\ -2u_{i,j,k}^n+u_{i,j-1,k}^n+u_{i,j,k+1}^n-2u_{i,j,k}^n+u_{i,j,k-1}^n).\end{array}$ (21)

$\begin{array}{c}{u}_{i,j,k}^{n+1}=2u_{i,j,k}^n-u_{i,j,k}^{n-1}+{\lambda }^2(u_{i+1,j,k}^n+u_{i-1,j,k}^n+u_{i,j+1,k}^n\\ +u_{i,j-1,k}^n+u_{i,j,k+1}^n+u_{i,j,k-1}^n-6u_{i,j,k}^n)\end{array}$ (22)

To ensure stability, the CFL condition must be satisfied:

$\lambda \le \frac{1}{\sqrt{3}}\text{or}\text{equivalently}\text{Δ}t\le \frac{h}{c\sqrt{3}}$

2.2. Boundary Conditions

$\{\begin{array}{l}u\left(0,y,z,t\right)=u\left(L_x,y,z,t\right)=0,\\ u\left(x,0,z,t\right)=u\left(x,L_y,z,t\right)=0,\\ u\left(x,y,0,t\right)=u\left(x,y,L_z,t\right)=0,\end{array}$ (23)

$\{\begin{array}{l}{u}_x\left(0,y,z,t\right)=u_x\left(L_x,y,z,t\right)=0,\\ u_y\left(x,0,z,t\right)=u_y\left(x,L_y,z,t\right)=0,\\ u_z\left(x,y,0,t\right)=u_z\left(x,y,L_z,t\right)=0.\end{array}$ (24)

Equations (23) and (24) impose both Dirichlet and Neumann conditions on identical boundaries. In the present work, this combination is employed to model a clamped boundary condition, as encountered in elastic membrane or plate problems where the boundary is rigidly fixed (zero displacement) and cannot rotate (zero slope). This scenario is physically realizable for fourth-order problems but is extended here to the second-order wave equation under the assumption of compatible initial data satisfying $f=0$ and $\partial f/\partial n=0$ on $\partial \Omega$ . The inclusion of both conditions also facilitates testing of the numerical scheme’s robustness under mixed boundary data.

3. Stability and Convergence Analysis

The paper now analyzes the stability of the explicit central-time central-space (leapfrog) scheme. It is worth emphasizing that the classical FTCS nomenclature can be misleading in this context. The scheme under investigation is a three-level method requiring two initial conditions ($u^0$ and $u^1$ ), whereas the standard FTCS scheme for parabolic problems is a two-level method. The von Neumann stability analysis presented below accounts for the quadratic characteristic equation arising from the three-level time discretization, which yields two amplification modes a feature absent in one-step FTCS methods. Consider the three-dimensional wave Equation (22) in finite difference form. Applying the von Neumann stability analysis by assuming the error solution of the form:

$\begin{array}{c}{u}_{i,j,k}^n=G^n{\text{e}}^{\text{i}\left(k_x{x}_i+k_y{y}_j+k_z{z}_k\right)}\\ =G^n{\text{e}}^{\text{i}\left(k_x\left(a_x+i\Delta x\right)+k_y\left(a_y+j\Delta y\right)+k_z\left(a_z+k\Delta z\right)\right)}\end{array}$ (25)

Substituting the Fourier mode into each term:

$u_{i,j,k}^{n+1}=G^{n+1}{\text{e}}^{\text{i}\left(k_x\left(a_x+i\Delta x\right)+k_y\left(a_y+j\Delta y\right)+k_z\left(a_z+k\Delta z\right)\right)}$ (26)

$u_{i,j,k}^{n-1}=G^{n-1}{\text{e}}^{\text{i}\left(k_x\left(a_x+i\Delta x\right)+k_y\left(a_y+j\Delta y\right)+k_z\left(a_z+k\Delta z\right)\right)}$ (27)

$\begin{array}{c}{u}_{i+1,j,k}^n=G^n{\text{e}}^{\text{i}\left[k_x\left(a_x+\left(i+1\right)\Delta x\right)+k_y\left(a_y+j\Delta y\right)+k_z\left(a_z+k\Delta z\right)\right]}\\ =G^n{\text{e}}^{\text{i}{k}_x\Delta x}{\text{e}}^{\text{i}\left[k_x\left(a_x+i\Delta x\right)+k_y\left(a_y+j\Delta y\right)+k_z\left(a_z+k\Delta z\right)\right]}\end{array}$ (28)

$\begin{array}{c}{u}_{i,j+1,k}^n=G^n{\text{e}}^{\text{i}\left[k_x\left(a_x+i\Delta x\right)+k_y\left(a_y+\left(j+1\right)\Delta y\right)+k_z\left(a_z+k\Delta z\right)\right]}\\ =G^n{\text{e}}^{\text{i}{k}_y\Delta y}{\text{e}}^{\text{i}\left[k_x\left(a_x+i\Delta x\right)+k_y\left(a_y+j\Delta y\right)+k_z\left(a_z+k\Delta z\right)\right]}\end{array}$ (29)

$\begin{array}{c}{u}_{i,j,k+1}^n=G^n{\text{e}}^{\text{i}\left[k_x\left(a_x+i\Delta x\right)+k_y\left(a_y+j\Delta y\right)+k_z\left(a_z+\left(k+1\right)\Delta z\right)\right]}\\ =G^n{\text{e}}^{\text{i}{k}_z\Delta z}{\text{e}}^{\text{i}\left[k_x\left(a_x+i\Delta x\right)+k_y\left(a_y+j\Delta y\right)+k_z\left(a_z+k\Delta z\right)\right]}\end{array}$ (30)

$u_{i-1,j,k}^n=G^n{\text{e}}^{\text{i}\left(k_x\left(a_x+\left(i-1\right)\Delta x\right)+k_y\left(a_y+j\Delta y\right)+k_z\left(a_z+k\Delta z\right)\right)}$ (31)

$=G^n{\text{e}}^{-\text{i}{k}_x\Delta x}{\text{e}}^{\text{i}\left(k_x\left(a_x+i\right)\Delta x+k_y\left(a_y+j\Delta y\right)+k_z\left(a_z+k\Delta z\right)\right)}$ (32)

$u_{i,j-1,k}^n=G^n{\text{e}}^{\text{i}\left(k_x\left(a_x+i\Delta x\right)+k_y\left(a_y+\left(j+1\right)\Delta y\right)+k_z\left(a_z+k\Delta z\right)\right)}$ (33)

$=G^n{\text{e}}^{-\text{i}{k}_y\Delta y}{\text{e}}^{\text{i}\left(\left(k_x\left(a_x+i\right)\Delta x\right)+k_y\left(a_y+j\Delta y\right)+k_z\left(a_z+k\Delta z\right)\right)}$ (34)

$u_{i,j,k-1}^n=G^n{\text{e}}^{\text{i}\left(k_x\left(a_x+i\Delta x\right)+k_y\left(a_y+j\Delta y\right)+k_z\left(a_z+\left(k+1\right)\Delta z\right)\right)}$ (35)

$=G^n{\text{e}}^{-\text{i}{k}_z\Delta z}{\text{e}}^{\text{i}\left(\left(k_x\left(a_x+i\right)\Delta x\right)+k_y\left(a_y+j\Delta y\right)+k_z\left(a_z+k\Delta z\right)\right)}$ (36)

where ${\text{i}}^2=-1$ , $G$ is the amplification factor and $k_x,k_y,k_z$ are wave numbers in each spatial direction.

Substituting into the scheme (21) leads to the relation

$\begin{array}{l}{G}^{n+1}{\text{e}}^{\text{i}\left(k_x\left(a_x+i\Delta x\right)+k_y\left(a_y+j\Delta y\right)+k_z\left(a_z+k\Delta z\right)\right)}\\ =2G^n{\text{e}}^{\text{i}\left(k_x\left(a_x+i\Delta x\right)+k_y\left(a_y+j\Delta y\right)+k_z\left(a_z+k\Delta z\right)\right)}-G^{n-1}{\text{e}}^{\text{i}\left(k_x\left(a_x+i\Delta x\right)+k_y\left(a_y+j\Delta y\right)+k_z\left(a_z+k\Delta z\right)\right)}\\ +{\lambda }^2(G^n{\text{e}}^{\text{i}{k}_x\Delta x}{\text{e}}^{\text{i}\left(k_x\left(a_x+i\Delta x\right)+k_y\left(a_y+j\Delta y\right)+k_z\left(a_z+k\Delta z\right)\right)}+G^n{\text{e}}^{-\text{i}{k}_x\Delta x}{\text{e}}^{\text{i}\left(k_x\left(a_x+i\Delta x\right)+k_y\left(a_y+j\Delta y\right)+k_z\left(a_z+k\Delta z\right)\right)}\\ +G^n{\text{e}}^{\text{i}{k}_y\Delta y}{\text{e}}^{\text{i}\left(k_x\left(a_x+i\Delta x\right)+k_y\left(a_y+j\Delta y\right)+k_z\left(a_z+k\Delta z\right)\right)}+G^n{\text{e}}^{-\text{i}{k}_y\Delta y}{\text{e}}^{\text{i}\left(k_x\left(a_x+i\Delta x\right)+k_y\left(a_y+j\Delta y\right)+k_z\left(a_z+k\Delta z\right)\right)}\\ +G^n{\text{e}}^{\text{i}{k}_z\Delta z}{\text{e}}^{\text{i}\left(k_x\left(a_x+i\Delta x\right)+k_y\left(a_y+j\Delta y\right)+k_z\left(a_z+k\Delta z\right)\right)}+G^n{\text{e}}^{-\text{i}{k}_z\Delta z}{\text{e}}^{\text{i}\left(k_x\left(a_x+i\Delta x\right)+k_y\left(a_y+j\Delta y\right)+k_z\left(a_z+k\Delta z\right)\right)}\\ -6G^n{\text{e}}^{\text{i}\left(k_x\left(a_x+i\Delta x\right)+k_y\left(a_y+j\Delta y\right)+k_z\left(a_z+k\Delta z\right)\right)})\end{array}$ (37)

Factor out the common term

${\text{e}}^{\text{i}\left(k_x\left(a_x+i\Delta x\right)+k_y\left(a_y+j\Delta y\right)+k_z\left(a_z+k\Delta z\right)\right)}$

from both sides of Equation (37) to get

$G^{n+1}=2G^n-G^{n-1}+{\lambda }^2{G}^n\left({\text{e}}^{\text{i}{k}_x\Delta x}+{\text{e}}^{-\text{i}{k}_x\Delta x}+{\text{e}}^{\text{i}{k}_y\Delta y}+{\text{e}}^{-\text{i}{k}_y\Delta y}+{\text{e}}^{\text{i}{k}_z\Delta z}+{\text{e}}^{-\text{i}{k}_z\Delta z}-6\right)$ (38)

Applying the Euler’s identity

${\text{e}}^{\text{i}\varphi }+{\text{e}}^{-\text{i}\varphi }=2\cos \left(\varphi \right),$

Equation (38) becomes

$G^{n+1}=2G^n-G^{n-1}+2{\lambda }^2{G}^n\left(\text{cos}\left(k_x\text{Δ}x\right)+\text{cos}\left(k_y\text{Δ}y\right)+\text{cos}\left(k_z\text{Δ}z\right)-3\right)$ (39)

Using a uniform grid spacing

$\text{Δ}x=\text{Δ}y=\text{Δ}z=h$

Equation (39) becomes

$G^{n+1}=2G^n-G^{n-1}+2{\lambda }^2{G}^n\left(\text{cos}\left(k_{x}h\right)+\text{cos}\left(k_{y}h\right)+\text{cos}\left(k_{z}h\right)-3\right)$ (40)

Suppose a parameter

$\mu =2{\lambda }^2\left(\text{cos}\left(k_{x}h\right)+\text{cos}\left(k_{y}h\right)+\text{cos}\left(k_{z}h\right)-3\right)$ (41)

Then Equation (40) becomes

$G^{n+1}=\left(2+\mu \right)G^n-G^{n-1}$

Assume a solution of the form $G^n={\xi }^n$ . Substituting gives:

${\xi }^{n+1}=\left(2+\mu \right){\xi }^n-{\xi }^{n-1}\Rightarrow {\xi }^2-\left(2+\mu \right)\xi +1=0$

This characteristic equation has solutions:

$\xi =\frac{\left(2+\mu \right)\pm \sqrt{{\left(2+\mu \right)}^2-4}}{2}$

For the numerical scheme to be stable, the amplification factor must satisfy:

$\left|\xi \right|\le 1\Rightarrow {\left(2+\mu \right)}^2-4\le 0\Rightarrow -4\le \mu \le 0$

Substituting the expression for $\mu$ :

$-4\le 2{\lambda }^2\left(\text{cos}\left(k_{x}h\right)+\text{cos}\left(k_{y}h\right)+\text{cos}\left(k_{z}h\right)-3\right)\le 0$ (42)

Since $\left|\text{cos}\left(m\right)\right|\le 1$ , then

${\lambda }^2\le \frac{1}{3}\Rightarrow \lambda \le \frac{1}{\sqrt{3}}\text{or}\lambda \ge \frac{-1}{\sqrt{3}}$ (43)

or equivalently:

$\left|\lambda \right|\le \frac{1}{\sqrt{3}}\Rightarrow \text{Δ}t\le \frac{h}{c\sqrt{3}}\left|\xi \right|\le 1$ (44)

The stability of the three-dimensional FTCS scheme is governed by the (CFL) condition (44). This fundamental stability criterion prevents numerical instabilities that could lead to oscillatory solutions or divergence in the computational results.

The operator family $C\left(\text{Δ}t\right)$ is said to be a convergent approximation to the true solution operator $E\left(t\right)$ if, for any initial value $u_0\left(x,y,z,t\right)\in ℬ$ , and for any sequence of time steps ${\text{Δ}}_{j}t$ and integers $n_j$ satisfying:

${\text{Δ}}_{j}t\to 0\text{and}{n}_j{\text{Δ}}_{j}t\to t\in \left[0,T\right],$

then the iterated approximation converges to the true solution:

$\Vert {\left(C\left({\text{Δ}}_{j}t\right)\right)}^{n_j}{u}_0\left(x,y,z,t\right)-E\left(t\right)u_0\left(x,y,z,t\right)\Vert \to 0,\text{for}\text{all}t\in \left[0,T\right].$

where $E\left(t\right)$ is the exact evolution operator, $u_0\left(x,y,z,t\right)\in ℬ$ is the initial condition, $C\left(\text{Δ}t\right)$ is a numerical approximation of the evolution over a time step $\text{Δ}t$ and

$u_{i,j,k}^n\left(x,y,z,t\right)=C{\left(\text{Δ}t\right)}^n{u}_0\left(x,y,z,t\right).$

after iterating the approximate operator $C\left(\text{Δ}t\right)$ $n$ times which gives the numerical approximation at time $t=n\text{Δ}t$ . The hope is that as $\text{Δ}t\to 0$ , this approximation converges to the true solution.

A numerical scheme such as FTCS is said to be stable if, under successive refinements of the time step ${\text{Δ}}_{j}t\to 0$ , the computed solution remains bounded over the time interval of interest.

For such a scheme, the operators used during the computation are drawn from the set

$\left\{C{\left({\text{Δ}}_{j}t\right)}^n\right\},j=1,2,3,\cdots ,$

where $n$ satisfies $0\le n{\text{Δ}}_{j}t<T$ . These operators act on the initial condition $u_0$ to produce approximations at discrete time levels.

Stability refers to the property that no component of the initial data is allowed to grow without bound due to the numerical procedure.

The approximation $C\left(\text{Δ}t\right)$ is said to be stable if the set of operators $\left\{C{\left({\text{Δ}}_{j}t\right)}^n\right\}$ is uniformly bounded. This means there exists a constant $K>0$ such that

$\Vert C{\left({\text{Δ}}_{j}t\right)}^n\Vert \le K,\text{for}\text{all}j\text{and}\text{all}n\text{such}\text{that}0\le n{\text{Δ}}_{j}t<T.$

The operator norm $\Vert C{\left(\text{Δ}t\right)}^n\Vert$ depends continuously on $\text{Δ}t$ for very small values. Therefore, $C\left(\text{Δ}t\right)$ is stable if there exists $\tau >0$ such that for all $0<\text{Δ}t\le \tau$ and all $n$ with $n\text{Δ}t\le T$ , the operators

$\left\{C{\left(\text{Δ}t\right)}^n\right\}$

remain uniformly bounded. Specifically, there exists a constant $K>0$ such that

$\Vert C{\left(\text{Δ}t\right)}^n\Vert \le K,\text{for}\text{all}\text{Δ}t\in \left(0,\tau \right]\text{and}\text{all}n\text{with}n\text{Δ}t\le T.$

Theorem 1. Lax Equivalence: For a uniformly solvable linear finite difference scheme that approximates a well-posed linear evolution problem, stability constitutes both a necessary and sufficient condition for convergence, provided the scheme is consistent

Proof. Proof of Sufficiency of the Lax Equivalence Theorem Let $u\left(x,y,z,t\right)\in \mathbb{M}$ be a sufficiently smooth solution of the 3D wave equation in finite differences, then,

$P_1\left(U^{m+1}-u^{m+1}\right)=P_0\left(U^m-u^m\right)-T^m,$

$U^{m+1}-u^{m+1}=\left(P_1^{-1}{P}_0\right)\left(U^m-u^m\right)-P_1^{-1}{T}^m.$

by truncation error where $\mathbb{M}$ is a Banach space, $P_0$ and $P_1$ are difference operators which are not functions of on $m$ and

$T^m=P_1{u}^{m+1}-\left[P_0{u}^m+F^m\right],$

is the truncation error of the finite difference scheme. Recursively, and by assuming $U^0=u^0$ , obtaining

$U^m-u^m=-{\displaystyle \sum }_{l=0}^{m-1}{\left(P_1^{-1}{P}_0\right)}^l{P}_1^{-1}{T}^{m-l-1}.$

Therefore, by the uniform solvability

$\Vert P_1^{-1}\Vert \le K\tau$

and the stability

$\Vert {\left(P_1^{-1}{P}_0\right)}^l\Vert \le K_1$ ,

this gives

$\Vert U^m-u^m\Vert \le K{K}_1\tau {\displaystyle \sum }_{l=0}^{m-1}\Vert T^l\Vert ,\forall m>0,K>0$

Assuming consistency, then

$\underset{\tau \left(h\right)\to 0}{\text{lim}}\Vert U^m-u^m\Vert =0,0\le m\tau \le t_{\text{max}}.$

For a general solution $u\left(x,y,z,t\right)$ , let ${\zeta }_{\alpha }\left(x,y,z,t\right)$ be the smooth solution sequence satisfying

$\underset{\alpha \to \infty }{\text{lim}}\Vert {\zeta }_{\alpha }^0-u^0\Vert \to 0.$

$\forall \epsilon >0$ , $\exists A>0$ , such that $\Vert {\zeta }_{\alpha }^0-u^0\Vert <\epsilon$ for all $\alpha >A$ .

For fixed $\beta >A$ , let ${\zeta }_{\beta }^m\left(x,y,z,t\right)$ be the solution of the difference.

For fixed $\beta >A$ , let ${\zeta }_{\beta }^m$ be the solution of the difference scheme with ${\zeta }_{\beta }^0={\zeta }_{1_{\beta }}^0$ .

$\forall \epsilon >0$ , $\exists h\left(\epsilon \right)>0$ , s.t. $\Vert {\zeta }_{\beta }^m-{\zeta }_{1_{\beta }}^m\Vert <\epsilon$ , for all $h<h\left(\epsilon \right)$ .

Thus, by the stability and the uniform invertibility of the scheme and the well-posedness of the problem that, if $h<h\left(\epsilon \right)$ , then

$\begin{array}{c}\Vert U^m-u^m\Vert \le \Vert U^m-{\zeta }_{\beta }^m\Vert +\Vert {\zeta }_{\beta }^m-{\zeta }_{1_{\beta }}^m\Vert +\Vert {\zeta }_{1_{\beta }}^m-u^m\Vert \\ \le \left(K_1+1+C\right)\epsilon .\end{array}$

then

$\underset{\tau \left(h\right)\to 0}{\text{lim}}\Vert U^m-u^m\Vert =0,0\le m\tau \le t_{\text{max}}.$

since $\epsilon$ is picked arbitrarily.

To complete the stability analysis and establish convergence of the explicit central-time central-space (leapfrog) scheme for the three-dimensional wave equation, we now examine its consistency. Consistency quantifies how accurately the discrete difference operators approximate the continuous partial differential equation at each grid point.

Local Truncation Error

Let $u\left(x,y,z,t\right)$ be a sufficiently smooth solution of the three-dimensional wave equation

$\frac{{\partial }^{2}u}{\partial t^2}-c^2\left(\frac{{\partial }^{2}u}{\partial x^2}+\frac{{\partial }^{2}u}{\partial y^2}+\frac{{\partial }^{2}u}{\partial z^2}\right)=0.$

local truncation error $T_{i,j,k}^n$ as the residual obtained when the exact solution is substituted into the finite difference scheme (21). Assuming uniform grid spacing $\text{Δ}x=\text{Δ}y=\text{Δ}z=h$ and time step $\text{Δ}t$ , the scheme is:

$\begin{array}{l}\frac{u^{n+1}-2u^n+u^{n-1}}{\text{Δ}{t}^2}-c^2(\frac{u_{i+1,j,k}^n-2u_{i,j,k}^n+u_{i-1,j,k}^n}{h^2}+\frac{u_{i,j+1,k}^n-2u_{i,j,k}^n+u_{i,j-1,k}^n}{h^2}\\ +\frac{u_{i,j,k+1}^n-2u_{i,j,k}^n+u_{i,j,k-1}^n}{h^2})=T_{i,j,k}^n.\end{array}$

Expanding each term about the point $\left(x_i,y_j,z_k,t^n\right)$ using Taylor series with remainder:

$u^{n+1}=u+\text{Δ}t u_t+\frac{\text{Δ}{t}^2}{2}{u}_{tt}+\frac{\text{Δ}{t}^3}{6}{u}_{ttt}+\frac{\text{Δ}{t}^4}{24}{u}_{tttt}+\mathcal{O}\left(\text{Δ}{t}^5\right),$

$u^{n-1}=u-\text{Δ}t u_t+\frac{\text{Δ}{t}^2}{2}{u}_{tt}-\frac{\text{Δ}{t}^3}{6}{u}_{ttt}+\frac{\text{Δ}{t}^4}{24}{u}_{tttt}+\mathcal{O}\left(\text{Δ}{t}^5\right),$

$u_{i\pm 1,j,k}=u\pm h u_x+\frac{h^2}{2}{u}_{xx}\pm \frac{h^3}{6}{u}_{xxx}+\frac{h^4}{24}{u}_{xxxx}+\mathcal{O}\left(h^5\right),$

$u_{i,j\pm 1,k}=u\pm h u_y+\frac{h^2}{2}{u}_{yy}\pm \frac{h^3}{6}{u}_{yyy}+\frac{h^4}{24}{u}_{yyyy}+\mathcal{O}\left(h^5\right),$

$u_{i,j,k\pm 1}=u\pm h u_z+\frac{h^2}{2}{u}_{zz}\pm \frac{h^3}{6}{u}_{zzz}+\frac{h^4}{24}{u}_{zzzz}+\mathcal{O}\left(h^5\right).$

Substituting these expansions into the finite difference stencil yields:

$\frac{u^{n+1}-2u^n+u^{n-1}}{\text{Δ}{t}^2}=u_{tt}+\frac{\text{Δ}{t}^2}{12}{u}_{tttt}+\mathcal{O}\left(\text{Δ}{t}^4\right),$

$\frac{u_{i+1,j,k}^n-2u_{i,j,k}^n+u_{i-1,j,k}^n}{h^2}=u_{xx}+\frac{h^2}{12}{u}_{xxxx}+\mathcal{O}\left(h^4\right),$

with analogous expressions for the $y$ and $z$ directions.

The local truncation error therefore evaluates to:

$\begin{array}{c}{T}_{i,j,k}^n=\left(u_{tt}-c^2\left(u_{xx}+u_{yy}+u_{zz}\right)\right)+\frac{\Delta t^2}{12}{u}_{tttt}-\frac{c^2{h}^2}{12}\left(u_{xxxx}+u_{yyyy}+u_{zzzz}\right)\\ +\mathcal{O}\left(\Delta t^4,h^4,\Delta t^2{h}^2\right).\end{array}$

The leading-order term in parentheses vanishes identically because $u$ satisfies the wave equation. Hence, the local truncation error reduces to:

$T_{i,j,k}^n=\frac{\text{Δ}{t}^2}{12}{u}_{tttt}-\frac{c^2{h}^2}{12}\left(u_{xxxx}+u_{yyyy}+u_{zzzz}\right)+\mathcal{O}\left(\text{Δ}{t}^4,h^4,\text{Δ}{t}^2{h}^2\right).$ (45)

4. Results and Discussion

To quantify stability, we monitor the discrete total energy. For the continuous wave equation, the total energy is conserved and given by:

$E\left(t\right)=\frac{1}{2}{\displaystyle {\int }_{\Omega }\left[u_t^2+c^2{\left|\nabla u\right|}^2\right]\text{d}\Omega }.$

The discrete counterpart is computed at each time level $n$ as:

$\begin{array}{c}{E}^n=\frac{1}{2}{\displaystyle \sum }_{i,j,k}[{\left(\frac{u_{i,j,k}^{n+1}-u_{i,j,k}^{n-1}}{2\Delta t}\right)}^2+c^2({\left(\frac{u_{i+1,j,k}^n-u_{i-1,j,k}^n}{2\Delta x}\right)}^2\\ +{\left(\frac{u_{i,j+1,k}^n-u_{i,j-1,k}^n}{2\Delta y}\right)}^2+{\left(\frac{u_{i,j,k+1}^n-u_{i,j,k-1}^n}{2\Delta z}\right)}^2)]\Delta x\Delta y\Delta z.\end{array}$

The first term approximates kinetic energy via central difference velocity; the remaining terms approximate potential energy via spatial gradients. Conservation of $E^n$ (boundedness with no secular growth) indicates numerical stability, while growth signals CFL violation. To demonstrate how the FTCS scheme performs under different stability regions, numerical simulations of the three-dimensional wave equation problem have been conducted, presented and analyzed. These sim-

ulations employ varying values of the combined stability parameter $\text{Δ}t\le \frac{h}{\sqrt{3}c}$ ,

showcasing the scheme’s behavior across different stability conditions fixing $\text{Δ}x=\text{Δ}y=\text{Δ}z=h=0.001$ and $\text{Δ}t=0.00001$ . The numerical investigation presents the evolution of total energy over time for a three-dimensional wave equation discretized using an explicit FDM [14]. The primary objective was to assess the stability and convergence behavior of the scheme under varying CFL numbers.

Numerical simulations of wave propagation exhibit strong dependence on the CFL number, as demonstrated by these results. Figures 1-4 show effects of varying CFL numbers (0.200, 0.300 and 0.500) less than the maximum CFL number needed for stability of the explicit numerical scheme, the results demonstrate consistent, oscillatory behavior in the total energy profiles [9] [16] [17]. The total energy remains bounded and exhibits no secular growth, which is characteristic of a numerically stable scheme. These results confirm that the discretization conserves energy to a satisfactory degree, thereby supporting both the stability and convergence of the method under sub-critical CFL conditions.

Figure 1. Energy distribution over time for CFL = 0.200.

Figure 2. Energy conservation for CFL = 0.200.

Figure 3. Energy conservation for CFL = 0.300.

Figure 4. Energy conservation for CFL = 0.500.

Figure 5. Energy conservation for CFL = 0.600.

Figure 6. Energy conservation for CFL = 0.800.

However, Figure 5 and Figure 6 show the effects of varying extreme CFL values (0.600 and 0.800) on the stability of the numerical method use and it can be observed that the numerical solution becomes unstable. This is evident in the monotonic growth of total energy over time. The presence of unbounded energy amplification is an indicator of numerical instability and a breakdown of the stability criterion.

4.1. Definition of Total Energy

To quantify stability, we monitor the discrete total energy. For the continuous wave equation, the total energy is conserved and given by:

$E\left(t\right)=\frac{1}{2}{\displaystyle {\int }_{\Omega }\left[u_t^2+c^2{\left|\nabla u\right|}^2\right]\text{d}\Omega }.$

The discrete counterpart is computed at each time level $n$ as:

$\begin{array}{c}{E}^n=\frac{1}{2}{\displaystyle \sum }_{i,j,k}[{\left(\frac{u_{i,j,k}^{n+1}-u_{i,j,k}^{n-1}}{2\text{Δ}t}\right)}^2+c^2({\left(\frac{u_{i+1,j,k}^n-u_{i-1,j,k}^n}{2\text{Δ}x}\right)}^2\\ +{\left(\frac{u_{i,j+1,k}^n-u_{i,j-1,k}^n}{2\text{Δ}y}\right)}^2+{\left(\frac{u_{i,j,k+1}^n-u_{i,j,k-1}^n}{2\text{Δ}z}\right)}^2)]\text{Δ}x\text{Δ}y\text{Δ}z.\end{array}$

4.2. Numerical Simulation Parameters

To ensure reproducibility of the numerical experiments presented in this section, the article documents below the complete set of simulation parameters, initial and boundary conditions, and discretization details used in all runs.

4.2.1. Computational Domain and Discretization

Domain extents: $\Omega =\left[0,L_x\right]\times \left[0,L_y\right]\times \left[0,L_z\right]$ with $L_x=L_y=L_z=1.0$ , Grid spacing: Uniform grid with $\text{Δ}x=\text{Δ}y=\text{Δ}z=h=0.001$ and Grid resolution: $N_x=N_y=N_z=1000$ grid points in each direction, yielding 10⁹ degrees of freedom per time level. Time step: Fixed at $\text{Δ}t=1\times {10}^{-5}$ s for all simulations. Wave speed: $c=1.0$ m/s (normalized). CFL number: Defined as $\lambda =c\text{Δ}t/h$ . For the above parameters, $\lambda =0.01$ is the baseline. To test stability limits, the reported CFL values (0.200, 0.300, 0.500, 0.600, 0.800) are achieved by artificially scaling the wave speed $c$ while holding $\text{Δ}t$ and $h$ fixed, or equivalently by reporting the $\lambda$ value directly. This approach isolates the effect of the CFL parameter without altering the spatial/temporal resolution. Final time: $T=0.01\text{s}$ (1000 time steps for baseline $\text{Δ}t$ ; proportionally fewer steps for larger $\text{Δ}t$ when $c$ is scaled).