Two Dimensional Wave FDM Simulation.¶
Mathematical Backgroud¶
If you haven’t already read the baby steps mathematical explanation used to create this animation, which focuses on how the Finite Difference Method is applied to solve the one-dimensional wave partial differential equation—rather than providing a general explanation of how the method itself works—please do so first. At the very least, make sure you understand the basic idea of how the Finite Difference Method is used in this context before continuing with this mathematical explanation.
Now that we did baby steps with one dimension, lets solve a more general version of the PDE mentioned above in a $\Omega \subseteq \mathbb{R}^2$ domain and the $(0,T]$ time period.
I am going to create a FDM scheme for the partial differential equation
\begin{align*} u_{tt} &= c^2 (u_{xx} + u_{yy} ) + f && (x,y) \in \Omega , t \in (0,T] \\ u(x,y,0) &= I(x,y) && \text{Initial Position in } (x,y) \in \Omega \\ u_t(x,y,0) &= V(x,y) && \text{Initial Velocity } (x,y)\in \Omega \\ u(x,y,t) &= 0 && \text{Dirichlet Boundary Conditions } (x,y) \in \partial \Omega, t \in (0,T] \end{align*}
For this problem we will create a mesh of the type
\begin{align*} x_i && \text{for } && i = 0, \dots ,N_x && \text{where } x_0 = 0 < x_1 < \dots < x_{N_x} = L \\ y_i && \text{for } && i = 0, \dots ,N_y && \text{where } y_0 = 0 < y_1 < \dots < y_{N_y} = L \\ t_n && \text{for } && n = 0, \dots ,N_t && \text{where } t_0 = 0 < t_1 < \dots < t_{N_t} = T \end{align*}
The solution $u(x,y,t)$ is sought at the meshpoints. For notation, we will denote $u_{i,j}^n$ the approximation at the point $(x_i,y_i,t_n)$. Given that the notation has been set, we will approximate $u_{tt},u_{xx}, u_{yy}$ as
\begin{align*} u_{tt} &\approx \frac{u_{i,j}^{n+1} - 2u_{i,j}^{n} + u_{i,j}^{n-1} }{\Delta t^2} \\ u_{xx} &\approx \frac{u_{i+1,j}^{n} - 2u_{i,j}^{n} + u_{i-1,j}^{n} }{\Delta x^2} \\ u_{yy} &\approx \frac{u_{i,j+1}^{n} - 2u_{i,j}^{n} + u_{i,j-1}^{n} }{\Delta x^2} \end{align*}
And we will denote
\begin{align*} D_x D_x = \frac{u_{i+1,j}^{n} - 2u_{i,j}^{n} + u_{i-1,j}^{n} }{\Delta x^2} \\ D_y D_y = \frac{u_{i,j+1}^{n} - 2u_{i,j}^{n} + u_{i,j-1}^{n} }{\Delta x^2} \end{align*}
Just like before we will use the centered finite difference to approximate $u_{t}(x_i,y_i,0)$ as
$$ u_{t}(x_i,y_i,0) \approx \frac{u_{i,j}^1 - u_{i,j}^{-1}}{2 \Delta t} $$
Similar to what we did in the 1D case, lets define the Courant numbers as
$C_x = \Delta t / \Delta x , C_y = \Delta t / \Delta y$ . Then, using the initial velocity condition we will find that
$$ u_{i,j}^1 = \frac{1}{2} \Delta t V_{i,j} + u_{i,j}^0 + \frac{1}{2} C_x^2 (u_{i+1,j}^0 - 2u_{i,j}^0 + u_{i-1,j}^{0}) + \frac{1}{2} C_y^2 (u_{i,j+1}^0 - 2u_{i,j}^0 + u_{i,j-1}^{0}) $$
And for the rest of the points the finite difference scheme will be
$$ u_{i,j}^{n+1} = -u_{i,j}^{n-1} + 2 u_{i,j}^{n} + c_x^2 (u_{i+1,j}^n - 2 u_{i,j}^n + u_{i-1,j}^n) + c_y^2 (u_{i,j+1}^n - 2 u_{i,j}^n + u_{i,j-1}^n) + \Delta t^2 f_{i,j}^n $$
And we produce the following python class
import numpy as np
from numpy import ndarray
from typing import Tuple, Callable
"""
Two-Dimensional Wave Equation Solver using Finite Difference Method (FDM)
This solver implements a numerical solution for the 2D wave equation with a source term:
u_tt = c^2 (u_xx + u_yy) + f(x,y,t)
where:
- (x,y) ∈ Ω (spatial domain)
- t ∈ (0, T] (time domain)
- c is the wave speed
- f(x,y,t) is a source/forcing function
Initial conditions:
u(x,y,0) = I(x,y) (initial displacement)
u_t(x,y,0) = V(x,y) (initial velocity)
Boundary conditions:
u = 0 on ∂Ω for all t ∈ (0, T] (Dirichlet boundary conditions)
The solver uses an explicit finite difference scheme with second-order accuracy
in both space and time.
"""
class TwoDimensionalWaveFDM:
_I: Callable[[float, float], float]
_V: Callable[[float, float], float]
_f: Callable[[float, float, float], float]
_Nx: int
_Ny: int
_Nt: int
_c: float
_Lx: float
_Ly: float
_T: float
_dx: float
_dy: float
_dt: float
_Cx: float
_Cy: float
_u: ndarray
_grid: Tuple[ndarray, ndarray, ndarray]
def __init__(self,I:Callable[[float,float],float],V:Callable[[float,float],float],f:Callable[[float,float,float], float],Nx:int, Ny:int, Nt:int,c:float=1.0,Lx=1.0,Ly=1.0,T=1.0) -> None:
self.I = I # Initial displacement function
self.V = V # Initial velocity function
self.f = f # Source/forcing function
self.Nx = Nx # Number of spatial grid points in x
self.Ny = Ny # Number of spatial grid points in y
self.Nt = Nt # Number of time steps
self.c = c # Wave speed
self.Lx = Lx # Length of the domain in x
self.Ly = Ly # Length of the domain in y
self.T = T # Total time
self.dx = Lx / (Nx - 1) # Spatial step size in x
self.dy = Ly / (Ny - 1) # Spatial step size in y
self.dt = T / (Nt - 1) # Time step size
self.Cx = (c * self.dt / self.dx) ** 2 # Courant number in x
self.Cy = (c * self.dt / self.dy) ** 2 # Courant number in y
self.u = np.zeros((Nx, Ny, Nt)) # Solution array
self.grid = self.create_grid() # Create spatial and temporal grid
self.apply_initial_conditions() # Apply initial conditions
def create_grid(self) -> Tuple[ndarray, ndarray, ndarray]:
x = np.linspace(0, self.Lx, self.Nx) # Spatial points in x
y = np.linspace(0, self.Ly, self.Ny) # Spatial points in y
t = np.linspace(0, self.T, self.Nt) # Time points
return np.meshgrid(x, y, t, indexing='ij')
def apply_boundary_conditions(self) -> None:
self.u[0,:,:] = 0 # u(0,y,t) = 0
self.u[-1,:,:] = 0 # u(Lx,y,t) = 0
self.u[:,0,:] = 0 # u(x,0,t) = 0
self.u[:,-1,:] = 0 # u(x,Ly,t) = 0
def apply_initial_displacement(self) -> None:
X,Y,_ = self.grid
xr = X[:,:,0]
yr = Y[:,:,0]
self.u[:,:,0] = self.I(xr, yr)
def apply_initial_velocity(self) -> None:
X,Y,_ = self.grid
for i in range(1, self.Nx - 1):
for j in range(1, self.Ny - 1):
xi = X[i,0,0]
yj = Y[0,j,0]
ft = 0.5*self.dt * self.V(xi, yj)+self.u[i,j,0]
p_x = 0.5*self.Cx * (self.u[i+1,j,0] - 2*self.u[i,j,0]+ self.u[i-1,j,0])
p_y = 0.5*self.Cy * (self.u[i,j+1,0] - 2*self.u[i,j,0]+ self.u[i,j-1,0])
lt = 0.5*self.dt**2 * self.f(xi, yj, 0)
self.u[i,j,1] = ft + p_x + p_y + lt
def solve(self) -> None:
X,Y,T = self.grid
for n in range(1,self.Nt-1):
for i in range(1,self.Nx-1):
for j in range(1,self.Ny-1):
xi = X[i,0,0]
yj = Y[0,j,0]
tn = T[0,0,n]
ft = -self.u[i,j,n-1] + 2*self.u[i,j,n]
p_x = self.Cx * (self.u[i+1,j,n] - 2*self.u[i,j,n] + self.u[i-1,j,n])
p_y = self.Cy * (self.u[i,j+1,n] - 2*self.u[i,j,n] + self.u[i,j-1,n])
lt = self.dt**2 * self.f(xi, yj, tn)
self.u[i,j,n+1] = ft + p_x + p_y + lt
def apply_initial_conditions(self) -> None:
self.apply_initial_displacement()
self.apply_boundary_conditions()
self.apply_initial_velocity()
# Getters and Setters
@property
def I(self) -> Callable[[float, float], float]:
"""Get the initial displacement function."""
return self._I
@I.setter
def I(self, I: Callable[[float, float], float]) -> None:
"""Set the initial displacement function."""
self._I = I
@property
def V(self) -> Callable[[float, float], float]:
"""Get the initial velocity function."""
return self._V
@V.setter
def V(self, V: Callable[[float, float], float]) -> None:
"""Set the initial velocity function."""
self._V = V
@property
def f(self) -> Callable[[float, float, float], float]:
"""Get the source/forcing function."""
return self._f
@f.setter
def f(self, f: Callable[[float, float, float], float]) -> None:
"""Set the source/forcing function."""
self._f = f
@property
def Nx(self) -> int:
"""Get the number of spatial grid points in x."""
return self._Nx
@Nx.setter
def Nx(self, Nx: int) -> None:
"""Set the number of spatial grid points in x."""
if Nx <= 0:
raise ValueError("Nx must be positive")
self._Nx = Nx
@property
def Ny(self) -> int:
"""Get the number of spatial grid points in y."""
return self._Ny
@Ny.setter
def Ny(self, Ny: int) -> None:
"""Set the number of spatial grid points in y."""
if Ny <= 0:
raise ValueError("Ny must be positive")
self._Ny = Ny
@property
def Nt(self) -> int:
"""Get the number of time steps."""
return self._Nt
@Nt.setter
def Nt(self, Nt: int) -> None:
"""Set the number of time steps."""
if Nt <= 0:
raise ValueError("Nt must be positive")
self._Nt = Nt
@property
def c(self) -> float:
"""Get the wave speed."""
return self._c
@c.setter
def c(self, c: float) -> None:
"""Set the wave speed."""
if c <= 0:
raise ValueError("Wave speed c must be positive")
self._c = c
@property
def Lx(self) -> float:
"""Get the length of the domain in x."""
return self._Lx
@Lx.setter
def Lx(self, Lx: float) -> None:
"""Set the length of the domain in x."""
if Lx <= 0:
raise ValueError("Lx must be positive")
self._Lx = Lx
@property
def Ly(self) -> float:
"""Get the length of the domain in y."""
return self._Ly
@Ly.setter
def Ly(self, Ly: float) -> None:
"""Set the length of the domain in y."""
if Ly <= 0:
raise ValueError("Ly must be positive")
self._Ly = Ly
@property
def T(self) -> float:
"""Get the total time."""
return self._T
@T.setter
def T(self, T: float) -> None:
"""Set the total time."""
if T <= 0:
raise ValueError("T must be positive")
self._T = T
@property
def dx(self) -> float:
"""Get the spatial step size in x."""
return self._dx
@dx.setter
def dx(self, dx: float) -> None:
"""Set the spatial step size in x."""
self._dx = dx
@property
def dy(self) -> float:
"""Get the spatial step size in y."""
return self._dy
@dy.setter
def dy(self, dy: float) -> None:
"""Set the spatial step size in y."""
self._dy = dy
@property
def dt(self) -> float:
"""Get the time step size."""
return self._dt
@dt.setter
def dt(self, dt: float) -> None:
"""Set the time step size."""
self._dt = dt
@property
def Cx(self) -> float:
"""Get the Courant number in x."""
return self._Cx
@Cx.setter
def Cx(self, Cx: float) -> None:
"""Set the Courant number in x."""
self._Cx = Cx
@property
def Cy(self) -> float:
"""Get the Courant number in y."""
return self._Cy
@Cy.setter
def Cy(self, Cy: float) -> None:
"""Set the Courant number in y."""
self._Cy = Cy
@property
def u(self) -> ndarray:
"""Get the solution array."""
return self._u
@u.setter
def u(self, u: ndarray) -> None:
"""Set the solution array."""
self._u = u
@property
def grid(self) -> Tuple[ndarray, ndarray, ndarray]:
"""Get the spatial and temporal grid."""
return self._grid
@grid.setter
def grid(self, grid: Tuple[ndarray, ndarray, ndarray]) -> None:
"""Set the spatial and temporal grid."""
self._grid = grid
Error Analysis.¶
For the error analysis, we will follow the same patterns as before, i.e.
- We will perform $m$ simulations
- The first simulation will be with a given $\Delta x_0,\Delta y_0, \Delta t_0$
- In each simulation $i=1,2,3,...$ we will decreased the size of $\Delta t_0$ by $\Delta t_i = 2^{-i} \Delta t_0$
- The value of the $u_{\text{approx}}(x_{N_x},y_{N_y},t_{N_x})$ approximation will be compared to the analytical value $u_{\text{exact}}(x_{N_x},y_{N_y},t_{N_x})$
- That will produce an array of errors $E_i$ which will be used with the $L^2$ norm to study the convergence error as $\Delta t$ decreases
But, this time, we will go an step further, we will estimate the convergence rates $r_{i}$ based on the consecutive experiments i.e.
- We will take the experiments $(\Delta t_{i-1},E_{i-1}),(\Delta t_{i},E_{i})$
- Based on the fact that we expect $E_i = c \left( \Delta t_i \right)^r$ and $E_{i-1} = c \left( \Delta t_{i-1} \right)^r$ where $r$ is the convergence rate
- The following sequence
$$ r_i = \frac{ \ln(E_{i-1}/E_i) }{\ln( \Delta t_{i-1}/t_i )} $$
Should converge to the convergence rate, as we are decreasing each $\Delta t_0$ by 2, then the sequence should converge to 2.
Lets do a code that implement this entire error analysis. In this case we will solve the PDE for the following conditions.
\begin{align*} I(x,y) &= 2 \sin(\pi x) \sin(\pi y) \\ V(x,y) &= 0 \\ f(x,y,t) &= (-1 + 2 \pi^2) \sin(\pi x) \sin(\pi y) \cos(t) \end{align*}
The domain region is $\Omega = [0,1] \times [0,1]$ and the time interval will be $T=(0,1]$
And the analytical solution will be
$$ u_e(x,y,t) = \sin(\pi x) \sin(\pi y) \cos(t) $$
import numpy as np
import matplotlib.pyplot as plt
import warnings
warnings.filterwarnings('ignore')
# --- Corrected Manufactured Solution ---
# u = sin(pi*x) * sin(pi*y) * cos(t)
# This ensures u_tt is not zero, allowing us to see O(dt^2)
def u_e(x, y, t):
return np.sin(np.pi * x) * np.sin(np.pi * y) * np.cos(t)
def I(x, y):
return np.sin(np.pi * x) * np.sin(np.pi * y)
def V(x, y):
# Velocity is du/dt = -sin(pi*x)*sin(pi*y)*sin(t). At t=0, this is 0.
return 0.0
def f(x, y, t):
# f = u_tt - c^2(u_xx + u_yy). Here c = 1.
# u_tt = -sin(pi*x) * sin(pi*y) * cos(t)
# u_xx = -pi^2 * sin(pi*x) * sin(pi*y) * cos(t)
# u_yy = -pi^2 * sin(pi*x) * sin(pi*y) * cos(t)
# f = [-1 - (1^2 * (-pi^2 - pi^2))] * u_e
return (-1.0 + 2.0 * np.pi**2) * np.sin(np.pi * x) * np.sin(np.pi * y) * np.cos(t)
# --- Simulation Parameters ---
m = 5 # Number of refinements
Nx0, Ny0, Nt0 = 10, 10, 40 # Start with a slightly finer base for stability
delta_t_array = []
l2_errors = []
for i in range(m):
# Halve steps by doubling intervals
Nx = (Nx0 - 1) * 2**i + 1
Ny = (Ny0 - 1) * 2**i + 1
Nt = (Nt0 - 1) * 2**i + 1
# Initialize and solve
simulator = TwoDimensionalWaveFDM(I, V, f, Nx, Ny, Nt, c=1.0, T=1.0)
simulator.solve()
X, Y, T = simulator.grid
dx, dy, dt = simulator.dx, simulator.dy, simulator.dt
# --- Corrected Error Analysis ---
# Compare the solution at the FINAL time step only
u_approx_final = simulator.u[:, :, -1]
u_exact_final = u_e(X[:, :, -1], Y[:, :, -1], T[:, :, -1])
# L2 Error at time T
error_final = u_exact_final - u_approx_final
l2_e = np.sqrt(dx * dy * np.sum(error_final**2))
l2_errors.append(l2_e)
delta_t_array.append(dt)
print(f"Level {i}: Grid {Nx}x{Ny}, dt={dt:.5f}, L2_Error={l2_e:.2e}")
# --- Calculate Rates ---
convergence_rates = [
np.log(l2_errors[i-1]/l2_errors[i]) / np.log(delta_t_array[i-1]/delta_t_array[i])
for i in range(1, len(l2_errors))
]
# --- Plotting Results ---
# Plot 1: Error vs Dt (Log-Log)
fig1, ax1 = plt.subplots(figsize=(10, 6))
ax1.loglog(delta_t_array, l2_errors, 'o-', label='L2 Error', linewidth=2.5,
markersize=8, color='#2E86AB')
ax1.set_xlabel('Time step size (Δt)', fontsize=12, fontweight='bold')
ax1.set_ylabel('L2 Error', fontsize=12, fontweight='bold')
ax1.set_title('Error Convergence: 2D Wave Equation FDM', fontsize=14, fontweight='bold', pad=15)
ax1.grid(True, alpha=0.3, linestyle='--', linewidth=0.7, which='both')
ax1.legend(loc='best', framealpha=0.9, fontsize=10, shadow=True)
ax1.set_facecolor('#F8F9FA')
fig1.patch.set_facecolor('white')
plt.tight_layout()
plt.show()
# Plot 2: Convergence Rates
fig2, ax2 = plt.subplots(figsize=(10, 6))
ax2.plot(range(1, len(l2_errors)), convergence_rates, 's-', color='#A23B72',
linewidth=2.5, markersize=8, label='Convergence Rate')
ax2.axhline(y=2.0, color='#2E86AB', linestyle='--', linewidth=1.5,
alpha=0.7, label='Expected O(h²)')
ax2.set_xlabel('Refinement Level', fontsize=12, fontweight='bold')
ax2.set_ylabel('Convergence Rate', fontsize=12, fontweight='bold')
ax2.set_title('Convergence Rate Analysis', fontsize=14, fontweight='bold', pad=15)
ax2.grid(True, alpha=0.3, linestyle='--', linewidth=0.7)
ax2.legend(loc='best', framealpha=0.9, fontsize=10, shadow=True)
ax2.set_facecolor('#F8F9FA')
fig2.patch.set_facecolor('white')
plt.tight_layout()
plt.show()
Level 0: Grid 10x10, dt=0.02564, L2_Error=4.45e-03 Level 1: Grid 19x19, dt=0.01282, L2_Error=1.09e-03 Level 2: Grid 37x37, dt=0.00641, L2_Error=2.70e-04 Level 3: Grid 73x73, dt=0.00321, L2_Error=6.74e-05 Level 4: Grid 145x145, dt=0.00160, L2_Error=1.68e-05
