PDE Matrix Generation

The PDE module provides tools for generating matrices from discretized partial differential equations in 1D, 2D, and 3D across multiple physical domains.

Overview

The module supports 17 different PDE types across various fields of physics and engineering:

Classical PDEs

• Poisson/Laplace equation

• Heat/Diffusion equation

• Wave equation

• Convection-Diffusion equation

• Helmholtz equation

• Biharmonic equation

Fluid Mechanics

• Stokes equations (incompressible viscous flow)

• Navier-Stokes equations (incompressible flow with convection)

• Burgers equation (1D & 2D nonlinear wave)

• Advection-Diffusion equation (thermal transport)

Solid Mechanics

• Linear Elasticity equations

Electromagnetics

• Maxwell equations (time-harmonic)

Reaction-Diffusion

• Linear Reaction-Diffusion

• Fisher-KPP equation (population dynamics)

• Gray-Scott model (pattern formation)

Quantum Mechanics

• Schrödinger equation

• Klein-Gordon equation (relativistic)

Features

• Dimensions: 1D, 2D, 3D problems

• Discretization: Finite differences (2nd, 4th, 6th order accuracy)

• Boundary Conditions: Dirichlet, Neumann, Periodic, Robin, Mixed

• Random Parameters: Generate batches with varying coefficients

• Backends: SciPy, NumPy, CuPy, JAX, PyTorch

• Formats: CSR, CSC, COO, Dense, etc.

• Time Integration: Various θ-methods (explicit, implicit, Crank-Nicolson)

Quick Start

Basic 1D Poisson Equation

from matrix_toolkit.pde import PDEMatrixGenerator, PDEConfig

# Configure
config = PDEConfig(
    dimension=1,
    mesh_size=100,
    backend='scipy',
    format='csr'
)

# Generate
gen = PDEMatrixGenerator('poisson', config)
A = gen.generate()

print(f"Shape: {A.shape}")
print(f"Sparsity: {1 - A.nnz/A.shape[0]**2:.2%}")

2D Navier-Stokes

config = PDEConfig(
    dimension=2,
    mesh_size=64,
    coefficients={
        'viscosity': 0.01,
        'dt': 0.01,
        'velocity_x0': 1.0,  # Background flow
        'velocity_y0': 0.0
    }
)

gen = PDEMatrixGenerator('navier_stokes', config)
A = gen.generate()  # Returns saddle-point system

Quantum Schrödinger

config = PDEConfig(
    dimension=1,
    mesh_size=512,
    coefficients={
        'potential': 'harmonic',
        'potential_strength': 1.0
    }
)

gen = PDEMatrixGenerator('schrodinger', config)
H = gen.generate()  # Hamiltonian matrix

Configuration

PDEConfig Class

Complete configuration reference:

from matrix_toolkit.pde import PDEConfig, BoundaryCondition

config = PDEConfig(
    # Domain
    dimension=2,                           # 1D, 2D, or 3D
    domain=(0.0, 1.0, 0.0, 2.0),          # Domain bounds
    mesh_size=(64, 128),                   # Grid points per dimension
    mesh_type='uniform',                   # uniform, stretched, random

    # Discretization
    discretization='finite_difference',
    stencil_order=2,                       # 2, 4, or 6

    # Boundary conditions
    boundary_condition=BoundaryCondition.DIRICHLET,

    # PDE coefficients
    coefficients={
        'viscosity': 0.01,
        'velocity': (1.0, 0.5),
        'dt': 0.01
    },

    # Random parameters
    random_params=False,
    random_seed=42,
    param_ranges={
        'viscosity': (0.001, 0.1)
    },

    # Output format
    backend='scipy',                       # scipy, numpy, cupy, jax, torch
    format='csr',                          # csr, csc, coo, dense
    dtype='float64',                       # float32, float64, complex64, complex128
)

Equation Types Reference

Classical PDEs

Poisson Equation

Equation: \(-\nabla^2 u = f\)

gen = PDEMatrixGenerator('poisson', config)

Applications: Electrostatics, steady heat conduction, potential flow

Coefficients:

• diffusion: Coefficient (default: 1.0)

Heat Equation

Equation: \(\frac{\partial u}{\partial t} = \alpha \nabla^2 u\)

config = PDEConfig(
    coefficients={
        'diffusion': 1.0,
        'dt': 0.01,
        'theta': 0.5  # Crank-Nicolson
    }
)
gen = PDEMatrixGenerator('heat', config)

Time Discretization:

• θ = 0: Forward Euler (explicit, conditionally stable)

• θ = 0.5: Crank-Nicolson (second-order, unconditionally stable)

• θ = 1: Backward Euler (implicit, unconditionally stable)

Wave Equation

Equation: \(\frac{\partial^2 u}{\partial t^2} = c^2 \nabla^2 u\)

config = PDEConfig(
    coefficients={
        'wave_speed': 1.0,
        'dt': 0.01
    }
)
gen = PDEMatrixGenerator('wave', config)

Helmholtz Equation

Equation: \(-\nabla^2 u - k^2 u = f\)

config = PDEConfig(
    coefficients={'wavenumber': 10.0}
)
gen = PDEMatrixGenerator('helmholtz', config)

Applications: Acoustics, wave propagation, electromagnetics

Biharmonic Equation

Equation: \(\nabla^4 u = f\) (or \(\nabla^2 \nabla^2 u = f\))

gen = PDEMatrixGenerator('biharmonic', config)

Applications: Plate bending, Stokes flow

Fluid Mechanics

Stokes Equations

Equations:

\[\begin{split}-\mu \nabla^2 \mathbf{u} + \nabla p = \mathbf{f} \\ \nabla \cdot \mathbf{u} = 0\end{split}\]

config = PDEConfig(
    dimension=2,  # or 3
    coefficients={'viscosity': 1.0}
)
gen = PDEMatrixGenerator('stokes', config)
A = gen.generate()  # Saddle-point system

Matrix Structure:

\[\begin{split}\begin{bmatrix} A & B^T \\ B & 0 \end{bmatrix} \begin{bmatrix} \mathbf{u} \\ p \end{bmatrix} = \begin{bmatrix} \mathbf{f} \\ 0 \end{bmatrix}\end{split}\]

where:

• \(A = \mu \nabla^2\) (velocity Laplacian)

• \(B = \nabla \cdot\) (divergence operator)

• \(\mathbf{u}\) = velocity (2 or 3 components)

• \(p\) = pressure

DOFs:

• 2D: \(2n_{vel} + n_{pres}\) where \(n_{vel} = n_x \times n_y\)

• 3D: \(3n_{vel} + n_{pres}\) where \(n_{vel} = n_x \times n_y \times n_z\)

Navier-Stokes Equations

Equations:

\[\begin{split}\frac{\partial \mathbf{u}}{\partial t} + (\mathbf{u} \cdot \nabla)\mathbf{u} = -\nabla p + \nu \nabla^2 \mathbf{u} + \mathbf{f} \\ \nabla \cdot \mathbf{u} = 0\end{split}\]

Linearized around background flow \(\mathbf{u}_0\):

\[\frac{\partial \mathbf{u}}{\partial t} + (\mathbf{u}_0 \cdot \nabla)\mathbf{u} = -\nabla p + \nu \nabla^2 \mathbf{u} + \mathbf{f}\]

config = PDEConfig(
    dimension=2,
    mesh_size=64,
    coefficients={
        'viscosity': 0.01,
        'dt': 0.01,
        'theta': 0.5,
        'velocity_x0': 1.0,  # Background flow
        'velocity_y0': 0.0
    }
)
gen = PDEMatrixGenerator('navier_stokes', config)

Reynolds Number: \(Re = \frac{UL}{\nu}\)

• Low Re: Stokes flow (viscous dominated)

• High Re: Turbulent flow (inertia dominated)

Burgers Equation

1D Burgers: \(\frac{\partial u}{\partial t} + u\frac{\partial u}{\partial x} = \nu \frac{\partial^2 u}{\partial x^2}\)

config = PDEConfig(
    dimension=1,
    mesh_size=256,
    coefficients={
        'viscosity': 0.01,
        'dt': 0.001,
        'u0': 0.0  # Linearization point
    }
)
gen = PDEMatrixGenerator('burgers', config)

2D Burgers:

\[\begin{split}\frac{\partial u}{\partial t} + u\frac{\partial u}{\partial x} + v\frac{\partial u}{\partial y} = \nu \nabla^2 u \\ \frac{\partial v}{\partial t} + u\frac{\partial v}{\partial x} + v\frac{\partial v}{\partial y} = \nu \nabla^2 v\end{split}\]

gen = PDEMatrixGenerator('burgers2d', config)

Applications: Shock wave formation, turbulence models

Advection-Diffusion

Equation: \(\frac{\partial T}{\partial t} + \mathbf{u} \cdot \nabla T = \alpha \nabla^2 T + Q\)

config = PDEConfig(
    dimension=2,
    coefficients={
        'diffusivity': 1.0,      # Thermal diffusivity
        'velocity': (1.0, 0.5),  # Flow velocity (given)
        'dt': 0.01
    }
)
gen = PDEMatrixGenerator('advection_diffusion', config)

Applications: Heat transport in flows, contaminant transport

Peclet Number: \(Pe = \frac{UL}{\alpha}\)

Solid Mechanics

Linear Elasticity

Equations: \(-\nabla \cdot \boldsymbol{\sigma}(\mathbf{u}) = \mathbf{f}\)

Stress tensor:

\[\boldsymbol{\sigma} = \lambda (\nabla \cdot \mathbf{u}) \mathbf{I} + \mu (\nabla \mathbf{u} + \nabla \mathbf{u}^T)\]

Lamé parameters:

\[\lambda = \frac{E\nu}{(1+\nu)(1-2\nu)}, \quad \mu = \frac{E}{2(1+\nu)}\]

config = PDEConfig(
    dimension=2,  # or 3
    mesh_size=32,
    coefficients={
        'youngs_modulus': 200e9,  # E (Pa) - e.g., steel
        'poisson_ratio': 0.3      # ν (dimensionless)
    }
)
gen = PDEMatrixGenerator('elasticity', config)

Material Properties:

Material	Young’s Modulus E (GPa)	Poisson’s Ratio ν
Steel	200	0.30
Aluminum	70	0.33
Copper	120	0.34
Rubber	0.01-0.1	0.48-0.50

DOFs:

• 2D: \(2n\) (ux, uy at each point)

• 3D: \(3n\) (ux, uy, uz at each point)

Electromagnetics

Maxwell Equations

Time-harmonic Maxwell equations:

\[\begin{split}\nabla \times \mathbf{E} - i\omega\mu \mathbf{H} = 0 \\ \nabla \times \mathbf{H} + i\omega\varepsilon \mathbf{E} = \mathbf{J}\end{split}\]

E-field formulation:

\[\nabla \times (\nabla \times \mathbf{E}) - \omega^2 \mu\varepsilon \mathbf{E} = -i\omega\mu\mathbf{J}\]

config = PDEConfig(
    dimension=3,  # Must be 3D
    mesh_size=16,
    coefficients={
        'frequency': 1e9,              # ω (rad/s)
        'permittivity': 8.854e-12,    # ε (F/m)
        'permeability': 1.257e-6      # μ (H/m)
    }
)
gen = PDEMatrixGenerator('maxwell', config)

Wave number: \(k = \omega\sqrt{\mu\varepsilon}\)

Wavelength: \(\lambda = \frac{2\pi}{k}\)

Reaction-Diffusion

Linear Reaction-Diffusion

Equation: \(\frac{\partial u}{\partial t} = D\nabla^2 u + \alpha u\)

config = PDEConfig(
    coefficients={
        'diffusion': 1.0,
        'reaction': -0.1,  # α < 0: decay
        'dt': 0.01
    }
)
gen = PDEMatrixGenerator('reaction_diffusion', config)

Fisher-KPP Equation

Equation: \(\frac{\partial u}{\partial t} = D\nabla^2 u + ru(1-u)\)

config = PDEConfig(
    coefficients={
        'diffusion': 1.0,
        'growth_rate': 1.0  # r
    }
)
gen = PDEMatrixGenerator('fisher_kpp', config)

Applications: Population dynamics, gene propagation, epidemic spreading

Gray-Scott Model

Two-species system:

\[\begin{split}\frac{\partial u}{\partial t} = D_u \nabla^2 u - uv^2 + F(1-u) \\ \frac{\partial v}{\partial t} = D_v \nabla^2 v + uv^2 - (F+k)v\end{split}\]

config = PDEConfig(
    dimension=2,
    mesh_size=128,
    coefficients={
        'diffusion_u': 2e-5,
        'diffusion_v': 1e-5,
        'feed_rate': 0.055,    # F
        'kill_rate': 0.062     # k
    }
)
gen = PDEMatrixGenerator('gray_scott', config)

Pattern Types (vary F and k):

• Spots: F=0.055, k=0.062

• Stripes: F=0.035, k=0.065

• Waves: F=0.014, k=0.054

Quantum Mechanics

Schrödinger Equation

Time-dependent: \(i\hbar \frac{\partial \psi}{\partial t} = \hat{H}\psi\)

Hamiltonian: \(\hat{H} = -\frac{\hbar^2}{2m}\nabla^2 + V(\mathbf{x})\)

config = PDEConfig(
    dimension=1,
    mesh_size=512,
    domain=(-10.0, 10.0),
    coefficients={
        'hbar': 1.0,
        'mass': 1.0,
        'potential': 'harmonic',      # or 'free', 'barrier', 'well'
        'potential_strength': 1.0
    }
)
gen = PDEMatrixGenerator('schrodinger', config)
H = gen.generate()  # Hamiltonian matrix

Potential Types:

• 'free': V(x) = 0

• 'harmonic': V(x) = kx²/2 (quantum harmonic oscillator)

• 'barrier': Rectangular potential barrier

• 'well': Potential well (infinite or finite)

Time Evolution: \(\psi(t+\Delta t) = e^{-i\hat{H}\Delta t/\hbar}\psi(t)\)

Using Crank-Nicolson:

\[\left(I + \frac{i\hat{H}\Delta t}{2\hbar}\right)\psi^{n+1} = \left(I - \frac{i\hat{H}\Delta t}{2\hbar}\right)\psi^n\]

Klein-Gordon Equation

Relativistic wave equation:

\[\frac{\partial^2 \phi}{\partial t^2} - c^2\nabla^2\phi + \frac{m^2c^4}{\hbar^2}\phi = 0\]

config = PDEConfig(
    dimension=2,
    coefficients={
        'mass': 1.0,
        'wave_speed': 1.0,
        'hbar': 1.0,
        'dt': 0.01
    }
)
gen = PDEMatrixGenerator('klein_gordon', config)

Applications: Relativistic quantum mechanics, particle physics

Advanced Features

Batch Generation with Random Parameters

Generate multiple matrices with varying parameters:

from matrix_toolkit.pde import PDEMatrixGenerator, PDEConfig

config = PDEConfig(
    dimension=2,
    mesh_size=64,
    random_params=True,
    random_seed=42,
    param_ranges={
        'viscosity': (0.001, 0.1),
        'velocity_x': (-2.0, 2.0),
        'velocity_y': (-2.0, 2.0)
    }
)

gen = PDEMatrixGenerator('navier_stokes', config)

# Generate 100 matrices with random parameters
matrices = gen.generate(n_matrices=100, randomize=True)

Latin Hypercube Sampling

For better parameter space coverage:

from matrix_toolkit.pde.random_params import LatinHypercubeSampler

sampler = LatinHypercubeSampler()
sampler.add_parameter('youngs_modulus', low=1e9, high=300e9)
sampler.add_parameter('poisson_ratio', low=0.2, high=0.45)

# Generate 50 samples
samples = sampler.sample(n=50, seed=42)

# Create matrices for each sample
for params in samples:
    config = PDEConfig(
        dimension=2,
        mesh_size=32,
        coefficients=params
    )
    gen = PDEMatrixGenerator('elasticity', config)
    A = gen.generate()

Multi-Physics Coupling

Combine multiple PDEs:

# Thermo-elasticity: heat + elasticity

# 1. Temperature field (heat equation)
config_heat = PDEConfig(
    dimension=2,
    mesh_size=64,
    coefficients={'diffusion': 1.0, 'dt': 0.01}
)
gen_heat = PDEMatrixGenerator('heat', config_heat)
M_heat = gen_heat.generate()

# 2. Thermal stress (elasticity with thermal expansion)
config_elastic = PDEConfig(
    dimension=2,
    mesh_size=64,
    coefficients={
        'youngs_modulus': 200e9,
        'poisson_ratio': 0.3
    }
)
gen_elastic = PDEMatrixGenerator('elasticity', config_elastic)
K_elastic = gen_elastic.generate()

Backend Conversion

Automatic conversion to different computational backends:

# SciPy sparse (CPU)
config_scipy = PDEConfig(
    dimension=2,
    mesh_size=128,
    backend='scipy',
    format='csr'
)
gen = PDEMatrixGenerator('poisson', config_scipy)
A_scipy = gen.generate()

# CuPy sparse (GPU)
config_cupy = PDEConfig(
    dimension=2,
    mesh_size=128,
    backend='cupy',
    format='csr'
)
gen = PDEMatrixGenerator('poisson', config_cupy)
A_cupy = gen.generate()  # Matrix on GPU

# JAX (for automatic differentiation)
config_jax = PDEConfig(
    dimension=2,
    mesh_size=128,
    backend='jax',
    format='dense'
)
gen = PDEMatrixGenerator('poisson', config_jax)
A_jax = gen.generate()

Performance Considerations

Matrix Sizes

Typical Matrix Sizes

Problem	Mesh	Matrix Size	Non-zeros	Memory (CSR)
1D Poisson	10,000	10⁴ × 10⁴	~30K	240 KB
2D Poisson	100×100	10⁴ × 10⁴	~50K	400 KB
2D Stokes	64×64	12K × 12K	~60K	1 MB
3D Poisson	50×50×50	125K × 125K	~875K	7 MB
3D Elasticity	32×32×32	98K × 98K	~650K	10 MB

Solver Recommendations

For small to medium problems (n < 10,000):

• Use direct solvers: scipy.sparse.linalg.spsolve

For large problems (n > 10,000):

• Use iterative solvers:

  • Symmetric positive definite: CG (scipy.sparse.linalg.cg)

  • Symmetric indefinite: MINRES

  • Non-symmetric: GMRES, BiCGSTAB

• Use preconditioners: ILU, AMG

For saddle-point systems (Stokes, Navier-Stokes):

• Block preconditioners

• Uzawa iteration

• Schur complement methods

GPU Acceleration

For very large problems:

import cupy as cp
from cupyx.scipy.sparse.linalg import cg

# Generate on GPU
config = PDEConfig(
    dimension=3,
    mesh_size=64,
    backend='cupy',
    format='csr'
)
gen = PDEMatrixGenerator('poisson', config)
A = gen.generate()

# Solve on GPU
b = cp.ones(A.shape[0])
x, info = cg(A, b)

Complete Examples

Example 1: Solving 2D Poisson

import numpy as np
import matplotlib.pyplot as plt
from scipy.sparse.linalg import spsolve
from matrix_toolkit.pde import PDEMatrixGenerator, PDEConfig

# Configure
config = PDEConfig(
    dimension=2,
    mesh_size=64,
    domain=(0.0, 1.0, 0.0, 1.0)
)

# Generate matrix
gen = PDEMatrixGenerator('poisson', config)
A = gen.generate()

# Create RHS: f(x,y) = 2π²sin(πx)sin(πy)
n = 64
h = 1.0 / 65
x = np.linspace(h, 1-h, n)
y = np.linspace(h, 1-h, n)
X, Y = np.meshgrid(x, y, indexing='ij')

f = 2 * np.pi**2 * np.sin(np.pi * X) * np.sin(np.pi * Y)
f = f.flatten() * h**2

# Solve
u = spsolve(A, f)
u_2d = u.reshape(n, n)

# Exact solution
u_exact = np.sin(np.pi * X) * np.sin(np.pi * Y)

# Plot
fig, axes = plt.subplots(1, 3, figsize=(15, 4))

im1 = axes[0].contourf(X, Y, u_2d, levels=20)
axes[0].set_title('Numerical Solution')
plt.colorbar(im1, ax=axes[0])

im2 = axes[1].contourf(X, Y, u_exact, levels=20)
axes[1].set_title('Exact Solution')
plt.colorbar(im2, ax=axes[1])

error = np.abs(u_2d - u_exact)
im3 = axes[2].contourf(X, Y, error, levels=20)
axes[2].set_title(f'Error (max={np.max(error):.2e})')
plt.colorbar(im3, ax=axes[2])

plt.tight_layout()
plt.show()

Example 2: Time-dependent Heat Equation

# Heat equation simulation
config = PDEConfig(
    dimension=2,
    mesh_size=64,
    coefficients={
        'diffusion': 0.1,
        'dt': 0.01,
        'theta': 0.5
    }
)

gen = PDEMatrixGenerator('heat', config)
M = gen.generate()  # System matrix

# Initial condition: Gaussian
n = 64
x = np.linspace(0, 1, n)
y = np.linspace(0, 1, n)
X, Y = np.meshgrid(x, y, indexing='ij')

u0 = np.exp(-50*((X-0.5)**2 + (Y-0.5)**2))
u = u0.flatten()

# Time stepping
n_steps = 100
for step in range(n_steps):
    # RHS for Crank-Nicolson
    # (I - θ dt D L) u^{n+1} = (I + (1-θ) dt D L) u^n
    # Simplification: just solve M u^{n+1} = u^n
    u = spsolve(M, u)

    if step % 20 == 0:
        u_2d = u.reshape(n, n)
        plt.figure()
        plt.contourf(X, Y, u_2d, levels=20)
        plt.colorbar()
        plt.title(f'Time step {step}')
        plt.show()

Example 3: Fluid Flow (Stokes)

# Lid-driven cavity flow
config = PDEConfig(
    dimension=2,
    mesh_size=32,
    coefficients={'viscosity': 0.01}
)

gen = PDEMatrixGenerator('stokes', config)
A = gen.generate()

# RHS: no body force
n_vel = 32 * 32 * 2  # 2 velocity components
n_pres = 32 * 32
b = np.zeros(n_vel + n_pres)

# Boundary conditions would be applied here
# (lid moving with unit velocity at top)

# Solve
sol = spsolve(A, b)

# Extract velocity and pressure
u = sol[:n_vel//2]
v = sol[n_vel//2:n_vel]
p = sol[n_vel:]

See Also

• PDE Module API Reference - Complete API reference

• ../examples/pde_examples - Basic examples

• ../examples/pde_advanced_examples - Advanced examples

• ../tutorials/pde_tutorial - Step-by-step tutorial