solver

solver.py — Spectral PDE solver for 1D/2D equations with pseudo‑differential operators

Overview

The solver module provides a high‑performance spectral solver for partial differential equations in one and two spatial dimensions. It is designed to handle a wide class of problems, including:

Linear and nonlinear PDEs with periodic or Dirichlet boundary conditions.
Time‑dependent evolution (first‑ or second‑order in time) and stationary problems.
Equations involving pseudo‑differential operators (via psiOp), such as fractional Laplacians or non‑local terms.
Source terms and arbitrary spatial derivatives.

The solver is built on Fourier spectral methods for spatial discretisation, ensuring exponential convergence for smooth solutions. Time stepping is performed with either a simple exponential integrator, a fourth‑order Exponential Time Differencing Runge–Kutta (ETD‑RK4) scheme, or explicit Runge–Kutta methods for nonlinear problems. For stationary pseudo‑differential equations, an asymptotic inversion of the operator is available.

Key features:

Symbolic parsing of the PDE using SymPy – the user provides an equation in standard mathematical form.
Automatic extraction of linear, nonlinear, source, and pseudo‑differential terms.
Computation of the Fourier symbol of the linear operator (dispersion relation) from the equation.
Dealiasing via sharp spectral cut‑off to control aliasing errors from nonlinear terms.
Efficient FFT‑based application of linear and pseudo‑differential operators.
Support for both 1D and 2D problems with uniform grids.
Built‑in energy monitoring, CFL checks, and spectral symbol analysis.
Rich visualisation: animated solutions, error plots, energy evolution, and 3D surface/2D image displays.
Testing framework for comparison with exact solutions.

Mathematical background

Spectral discretisation. The solution u(x,t) is approximated by its truncated Fourier series. Spatial derivatives are replaced by multiplication with the corresponding wavenumber in Fourier space:

∂/∂x → i k , ∂²/∂x² → -k² .

Nonlinear terms are evaluated in physical space and then transformed back (pseudo‑spectral approach). Dealiasing is applied by zeroing out the highest one‑third of the Fourier modes (the 2/3‑rule).

Linear operator. After parsing the PDE, the linear part is converted into a Fourier multiplier L(k) (or L(kx,ky)). For a first‑order‑in‑time equation

∂ₜ u = L u + N(u) + f(x,t),

the variation‑of‑constants formula motivates exponential integrators. The solver implements:

Exponential Euler (first order): uⁿ⁺¹ = e^{LΔt} uⁿ + Δt φ₁(LΔt) (N(uⁿ)+fⁿ) with φ₁(z) = (e^z-1)/z.
ETD‑RK4 (fourth order): a Runge–Kutta scheme that uses exponentials of the linear part, achieving high accuracy for stiff problems.

For second‑order equations

∂ₜ² u = L u + N(u) + f(x,t),

the system is reduced to first order in time, and a leap‑frog style integrator or a second‑order ETD scheme is used.

Pseudo‑differential operators. When the equation contains psiOp(p(x,ξ), u), the solver builds a PseudoDifferentialOperator object (from the companion psiop module) and applies it using Kohn–Nirenberg quantisation (fast Fourier for constant‑coefficient symbols, direct quadrature for spatially varying symbols). For stationary problems, an asymptotic right inverse is constructed via symbolic series expansion.

Stability. The solver checks the CFL condition based on the maximum group velocity derived from the dispersion relation. It also verifies that the symbol satisfies the necessary conditions for well‑posedness (Re L(k) ≤ 0, sufficient high‑frequency dissipation, moderate growth).

References

class solver.PDESolver(equation, time_scheme='default', dealiasing_ratio=0.6666666666666666)[source]

Bases: object

A partial differential equation (PDE) solver based on spectral methods using Fourier transforms.

This solver supports symbolic specification of PDEs via SymPy and numerical solution using high-order spectral techniques. It is designed for both linear and nonlinear time-dependent PDEs, as well as stationary pseudo-differential problems.

Key Features:

Symbolic PDE parsing using SymPy expressions
1D and 2D spatial domains with periodic boundary conditions
Fourier-based spectral discretization with dealiasing
Temporal integration schemes:
- Default exponential time stepping
- ETD-RK4 (Exponential Time Differencing Runge-Kutta of 4th order)
Nonlinear terms handled through pseudo-spectral evaluation
Built-in tools for:
- Visualization of solutions and error surfaces
- Symbol analysis of linear and pseudo-differential operators
- Microlocal analysis (e.g., Hamiltonian flows)
- CFL condition checking and numerical stability diagnostics

Supported Operators:

Linear differential and pseudo-differential operators
Nonlinear terms up to second order in derivatives
Symbolic operator composition and adjoints
Asymptotic inversion of elliptic operators for stationary problems

Example Usage:

>>> from PDESolver import *
>>> u = Function('u')
>>> t, x = symbols('t x')
>>> eq = Eq(diff(u(t, x), t), diff(u(t, x), x, 2) + u(t, x)**2)
>>> def _initial(x): return np.sin(x)
>>> solver = PDESolver(eq)
>>> solver.setup(Lx=2*np.pi, Nx=128, Lt=1.0, Nt=1000, initial_condition=initial)
>>> solver.solve()
>>> ani = solver.animate()
>>> HTML(ani.to_jshtml())  # Display animation in Jupyter notebook

animate(component='abs', overlay='contour', mode='surface')[source]

Create an animated plot of the solution evolution over time.

This method generates a dynamic visualization of the stored solution frames self.frames. It supports:

1D line animation (unchanged),

2D surface animation (original behavior, ‘surface’),

2D image animation using imshow (new, ‘imshow’) which is faster and often clearer for large grids.

Parameters:

component (str, optional, one of {'real', 'imag', 'abs', 'angle'}) –
Which component of the complex field to visualize:
- ’real’ : Re(u)
- ’imag’ : Im(u)
- ’abs’ : |u|
- ’angle’ : arg(u)
Default is ‘abs’.
overlay (str or None, optional, one of {'contour', 'front', None}) –
For 2D modes only. If None, no overlay is drawn.
- ’contour’ : draw contour lines on top (or beneath for 3D surface)
- ’front’ : detect and mark wavefronts using gradient maxima
Default is ‘contour’.
mode (str, optional, one of {'surface', 'imshow'}) – 2D rendering mode. ‘surface’ keeps the original 3D surface plot. ‘imshow’ draws a 2D raster (faster, often more readable). Default is ‘surface’ for backward compatibility.

Returns:

A Matplotlib FuncAnimation instance (you can display it in a notebook or save it to file).

Return type:

FuncAnimation

Notes

The method uses the same time-mapping logic as before (linear sampling of stored frames to animation frames).
For ‘angle’ the color scale is fixed between -π and π.
For other components, color scaling is by default dynamically adapted per frame in ‘imshow’ mode (this avoids extreme clipping if amplitudes vary).
Overlays are updated cleanly: previous contour/scatter artists are removed before drawing the next frame to avoid memory/visual accumulation.
Animation interval is 50 ms per frame (unchanged).

plot_energy(log=False)[source]

Plot the time evolution of the total energy for wave equations. Visualizes the energy computed during simulation for both 1D and 2D cases. Requires temporal_order=2 and prior execution of compute_energy() during solve().

Parameters:: log – bool If True, displays energy on a logarithmic scale to highlight exponential decay/growth.

Notes

Energy is defined as E(t) = 1/2 ∫ [ (∂ₜu)² + |L¹⸍²u|² ] dx
Only available if energy monitoring was activated in solve()
Automatically skips plotting if no energy data is available

Displays:

Time vs. Total Energy plot with grid and legend
Appropriate axis labels and dimensional context (1D/2D)
Logarithmic or linear scaling based on input parameter

setup(Lx, Ly=None, Nx=None, Ny=None, Lt=1.0, Nt=100, boundary_condition='periodic', initial_condition=None, initial_velocity=None, n_frames=100, plot=True)[source]

Configure the spatial/temporal grid and initialize the solution field.

This method sets up the computational domain, initializes spatial and temporal grids, applies boundary conditions, and prepares symbolic and numerical operators. It also performs essential analyses such as:

CFL condition verification (for stability)

Symbol analysis (e.g., dispersion relation, regularity)

Wave propagation analysis for second-order equations

If pseudo-differential operators (ψOp) are present, symbolic analysis is skipped in favor of interactive exploration via interactive_symbol_analysis.

Parameters:

Lx (float) – Size of the spatial domain along x-axis.
Ly (float, optional) – Size of the spatial domain along y-axis (for 2D problems).
Nx (int) – Number of spatial points along x-axis.
Ny (int, optional) – Number of spatial points along y-axis (for 2D problems).
Lt (float, default=1.0) – Total simulation time.
Nt (int, default=100) – Number of time steps.
initial_condition (callable) – Function returning the initial state u(x, 0) or u(x, y, 0).
initial_velocity (callable, optional) – Function returning the initial time derivative ∂ₜu(x, 0) or ∂ₜu(x, y, 0), required for second-order equations.
n_frames (int, default=100) – Number of time frames to store during simulation for visualization or output.

Raises:

ValueError – If mandatory parameters are missing (e.g., Nx not given in 1D, Ly/Ny not given in 2D).

Notes

The spatial discretization assumes periodic boundary conditions by default.
Fourier transforms are computed using real-to-complex FFTs (scipy.fft.fft, fft2).
Frequency arrays (KX, KY) are defined following standard spectral conventions.
Dealiasing is applied using a sharp cutoff filter at a fraction of the maximum frequency.
For second-order equations, initial acceleration is derived from the governing operator.
Symbolic analysis includes plotting of the symbol’s real/imaginary/absolute values and dispersion relation.