Meshless solutions of PDE inverse problems on irregular geometries

TL;DR

The paper addresses the challenge of solving inverse PDE problems on irregular geometries by introducing a mesh-free, time-inclusive spectral representation on a hyper-rectangular envelope and optimizing the tensor of basis coefficients through a Physics-Informed loss that enforces the PDE, boundary conditions, and data terms. Time is treated as just another spectral coordinate, enabling simultaneous forward, inverse, and control tasks across a wide range of PDEs, including stiff time-dependent problems and geophysical data assimilation, with empirical exponential convergence and substantially fewer parameters than neural-network baselines. The approach is demonstrated across Laplace, Allen–Cahn, nonlinear Schrödinger, wave, diffusion on curved surfaces, ice-sheet viscosity inversion, heat control, and optimal transport problems, highlighting flexibility in irregular domains via Fourier extension embeddings. While nonconvex optimization introduces multiple local minima and potential convergence variability, the method offers a scalable, mesh-free framework that substantially improves accessibility for inverse design and data assimilation in complex PDE settings.

Abstract

Solving inverse and optimization problems over solutions of nonlinear partial differential equations (PDEs) on complex spatial domains is a long-standing challenge. Here we introduce a method that parameterizes the solution using spectral bases on arbitrary spatiotemporal domains, whereby the basis is defined on a hyperrectangle containing the true domain. We find the coefficients of the basis expansion by solving an optimization problem whereby both the equations, the boundary conditions and any optimization targets are enforced by a loss function, building on a key idea from Physics-Informed Neural Networks (PINNs). Since the representation of the function natively has exponential convergence, so does the solution of the optimization problem, as long as it can be solved efficiently. We find empirically that the optimization protocols developed for machine learning find solutions with exponential convergence on a wide range of equations. The method naturally allows for the incorporation of data assimilation by including additional terms in the loss function, and for the efficient solution of optimization problems over the PDE solutions.

Figures (24)

• Figure 1: Schematic view of basis fitting scheme. We assume that any function can be written as a tensor product of basis terms with some set of coefficients, which are not initially known. We define a loss function over some geometry by computing the residuals of some physical equation and boundary conditions at a series of collocation points. The loss function can also include terms that require the solution to match data. Find a set of optimized coefficients $c_{ij}$ that minimize this loss function using a differentiable code and gradient-based nonlinear optimizer.
• Figure 2: (a) Solution to the Laplace equation in an irregular domain with a hole as found by fitting polynomial coefficients to the equation residuals with $N=60$ coefficients in both $x$ and $y$. (b) The same system solved using COMSOL with an extremely high mesh density. (c) Convergence of the basis fitting approach to the COMSOL solution as the number of basis modes is increased. (d) Solution of Laplace's equation when the boundary condition is not provided but 379 data samples from the interior of the domain (shown as the black dots) included in the loss with $N=45$ coefficients. (e) Boundary condition for the Laplace equation on a very irregular domain described by the perimeter of Lake Taal in the Philippines. (f) Solution of the Laplace equation on the surface of Lake Taal with $N=100$ basis modes in both $x$ and $y$.
• Figure 3: Solutions to stiff PDEs shown versus reference solutions. (a) shows a polynomial fit to the Allen-Cahn equations, with 200 basis modes in $x$ and 10 modes in $t$. (b) shows a reference solution from WangS2023, while (c) shows the convergence of the solution to the reference as the number of basis modes in both $x$ and $t$ are increased. (d) similarly shows a polynomial fit to the nonlinear Schrödinger equations with $N=45$ in both $x$ and $t$ compared to a reference solution (e) from Raissi2019. (f) shows the convergence of the polynomial fit to the reference as a function of the number of basis modes.
• Figure 4: Solution to the wave equation on an irregular domain at several time steps evaluated using a polynomial fit (top row) with $N=12$ polynomials in each dimension of $x$, $y$, and $t$, and COMSOL (bottom row) solved using a finite element method with extremely fine mesh density.
• Figure 5: Diffusion equation solved on the surface of a sphere. The surface was embedded in 3D space and the surface diffusion equation was solved with $N=10$ Chebyshev modes in each of $x$, $y$, $z$, and $t$. Simulations were run to $t=2.0$ with $D = 0.5$. (a) and (b) show the concentration field on the sphere from the non-axisymmetric initial condition at $t=0$ to the increasingly axisymmetric profile at $t=0.5$. (c) compares the our simulation results in the open circles with the exact analytic result (solid line) expressed in spherical harmonic functions with 20 modes. Results are shown at several time points in the simulation.