Skip to content

Building PDE Operators

The @operator macro lets you write PDE operators in mathematical notation. It translates symbolic expressions into composable operator objects, producing a single weight matrix that applies the full PDE operator in one matrix-vector multiply.

This is the recommended way to build operators the built-ins don't cover. If your operator can't be expressed here, see Custom Operators for the closure escape hatch.

julia
using RadialBasisFunctions
using StaticArrays

x = rand(SVector{2,Float64}, 100)
f(p) = sin(p[1]) * cos(p[2])
u = f.(x)

Recognized symbols

SymbolMeaning
∇², ΔLaplacian
∂(dim)First partial derivative in dimension dim
∂²(dim)Second partial derivative in dimension dim
∇ ⋅ (κ * ∇)Diffusion operator (scalar or vector κ)
c ⋅ ∇Advection operator (vector c)
f, IIdentity operator
Everything elseScalar coefficient

Standard arithmetic (+, -, *) and unary negation work as expected. Scalars can be literals, variables, or expressions like k^2 or c[1].

The macro produces rank-0 operators — see Understanding Rank (N) for what that means, and Custom Operators if you need rank 1.

Helmholtz Operator

The Helmholtz equation    appears in acoustics, electromagnetics, and quantum mechanics. The operator combines a Laplacian with a scaled identity.

julia
= 4.0

op = @operator ∇² +* f
helm_op = op(x)

# Verify against separate built-in operators
expected = laplacian(x)(u) .+.* u
maximum(abs, helm_op(u) .- expected)
2.0976553827267708e-12

Diffusion — Textbook Notation

The diffusion operator   appears in heat conduction, mass transfer, and many other physical models. The @operator macro recognizes the textbook form directly:

julia
κ = [2.0, 0.5]

op = @operator* ∇)
diff_op = op(x)

# Verify against separate built-in operators
expected = κ[1] .* partial(x, 2, 1)(u) .+ κ[2] .* partial(x, 2, 2)(u)
maximum(abs, diff_op(u) .- expected)
1.2363443602225743e-12

Scalar produces an isotropic operator (scaled Laplacian):

julia
op = @operator (3.0 * ∇)
diff_iso = op(x)  # equivalent to 3∇²f
expected = 3.0 .* laplacian(x)(u)
maximum(abs, diff_iso(u) .- expected)
4.035882739117369e-12

Anisotropic Diffusion — Explicit Partials

The same anisotropic diffusion can also be written with explicit per-dimension coefficients:

julia
κ_x = 2.0
κ_y = 0.5

op = @operator κ_x * ∂²(1) + κ_y * ∂²(2)
aniso_op = op(x)

# Verify against separate built-in operators
expected = κ_x .* partial(x, 2, 1)(u) .+ κ_y .* partial(x, 2, 2)(u)
maximum(abs, aniso_op(u) .- expected)
1.2363443602225743e-12

When  , this reduces to a scaled Laplacian.

Advection-Diffusion

The steady advection-diffusion equation     balances viscous diffusion against transport by a velocity field. It appears in fluid dynamics, pollutant transport, and thermal convection.

julia
ν = 0.01
c = SVector(1.0, 0.5)

op = @operator ν * ∇² - c 
advdiff_op = op(x)

# Verify against separate built-in operators
expected = ν .* laplacian(x)(u) .- c[1] .* partial(x, 1, 1)(u) .- c[2] .* partial(x, 1, 2)(u)
maximum(abs, advdiff_op(u) .- expected)
4.6629367034256575e-14

Sharing Stencils

When multiple operators act on the same point set, precompute the neighbor list once and pass it to avoid redundant nearest-neighbor searches:

julia
adjl = find_neighbors(x, 30)

helm_op  = (@operator ∇² +* f)(x; adjl=adjl)
aniso_op = (@operator κ_x * ∂²(1) + κ_y * ∂²(2))(x; adjl=adjl)
RadialBasisOperator
├─Operator: 2.0 × ∂ⁿf/∂xᵢ (n = 2, i = 1) + 0.5 × ∂ⁿf/∂xᵢ (n = 2, i = 2)
├─Data type: StaticArraysCore.SVector{2, Float64}
├─Number of points: 100
├─Stencil size: 30
└─Basis: Polyharmonic spline (r³) with degree 2 polynomial augmentation

Boundary Conditions on PDE Operators

Operators built with @operator accept the hermite keyword just like the built-ins — see Boundary Conditions for how to enforce Dirichlet, Neumann, or Robin conditions.