An FFT-based solver with general boundary conditions for stationary diffusion problems based on Chebyshev collocation

Abstract: An efficient and robust FFT-based solver is proposed for diffusion-type problems with general Neumann and Dirichlet boundary conditions, based on a Chebyshev collocation framework. The method combines Chebyshev polynomial approximations with FFT-based operators to provide a matrix-free implementation of the discrete differential operator at the Chebyshev-Gauss-Lobatto points. The linear system of equations resulting from the Chebyshev discretization is solved using LGMRES. To overcome convergence problems on fine grids, a hierarchical refinement strategy based on modal prolongation is proposed, enabling the solution of very large 3D problems. The methodology is applicable to homogeneous and heterogeneous domains, as well as to linear and nonlinear constitutive equations. The accuracy of the proposed method is analyzed by solving the Poisson equation in homogeneous 1D and 3D domains with general boundary conditions, using manufactured analytical solutions as references. Convergence to the analytical solution is achieved in a few iterations, with smaller errors than those obtained using DCT/DST approaches. Discretizations of up to $256^3$ are achieved thanks to the hierarchical refinement strategy. In the case of heterogeneous domains, the accuracy and efficiency obtained are similar to those of a standard periodic FFT approach. It is found that the computational complexity of the method preserves the FFT scaling, of order $n\log n$, in all the cases studied.

• The paper presents an FFT-based matrix-free Chebyshev collocation solver that achieves exponential convergence for stationary diffusion PDEs.
• It employs a modal hierarchical refinement strategy to overcome ill-conditioning and enables efficient resolution of large-scale 3D problems.
• The method seamlessly enforces general, spatially-varying Dirichlet and Neumann boundary conditions while maintaining O(N log N) computational complexity.

FFT-Based Chebyshev Collocation Solver for Stationary Diffusion with General Boundary Conditions

Introduction

The paper presents an advanced Fourier-based spectral collocation method utilizing Chebyshev polynomial discretization for efficiently solving stationary diffusion-type PDEs subject to arbitrary Dirichlet and Neumann boundary conditions. This work directly addresses the constraints of existing FFT solvers, which are fundamentally tied to periodic boundary conditions due to their trigonometric basis, limiting their use for general boundary value problems (BVPs) on bounded, non-periodic domains. The proposed framework combines the computational efficiency of FFT algorithms with the flexibility and high-order accuracy of Chebyshev polynomial collocation, introducing a matrix-free paradigm that extends FFT methodologies to BVPs with non-periodic and mixed boundary conditions.

A key contribution is an efficient modal hierarchical refinement strategy, which leverages prolongation in the Chebyshev basis to overcome ill-conditioning at fine discretization levels—a traditional bottleneck in spectral collocation—thus enabling solution of very large-scale 3D problems. The methodology is implemented for both linear and nonlinear diffusion equations and accommodates material heterogeneity without explicit matrix assembly.

Numerical Formulation

Chebyshev Collocation and FFT-Enabled Operators

The solver employs Chebyshev polynomials as the global basis due to their superior approximation properties for smooth non-periodic solutions. Collocation is performed at Chebyshev-Gauss-Lobatto (CGL) nodes, ensuring spectral convergence for analytic fields. The discrete differentiation operators required for the PDE are evaluated in a matrix-free fashion using (I)DCT and (I)DST transforms, enabling $O(N \log N)$ complexity akin to FFT-based approaches.

In 1D, the Chebyshev representation reformulates the variational problem into a collocation system, with boundary conditions imposed directly via additional constraints on the degrees of freedom. Multidimensional extension is achieved using tensor-product collocation grids, with multidirectional differential operators efficiently handled as separable 1D operations along each axis.

Matrix-Free Linear and Nonlinear Solvers

The linear systems arising from discretization are solved iteratively with LGMRES, avoiding explicit formation or storage of dense differentiation matrices and thus achieving scalability to high resolutions. For nonlinear equations, outer Picard or Newton iterations are applied, with Krylov solvers addressing the linearized subproblems at each step. For nonlinear constitutive laws (as in temperature-dependent conductivity), the Chebyshev discretization is employed consistently through the Jacobian evaluation required by Newton's method.

Enforcement of General Boundary Conditions

The method supports fully general, spatially varying Dirichlet and Neumann boundary conditions, imposing them through a systematic reordering and elimination scheme at the algebraic level—without recasting the problem to an auxiliary homogeneous setting or requiring problem-specific function bases. Eliminations and prolongations are performed implicitly, maintaining the matrix-free property except for small-size (boundary-localized) manipulations.

Modal Hierarchical Refinement Strategy

Ill-conditioning of high-order spectral collocation operators, which typically undermines convergence of iterative solvers, is addressed by a hierarchical spectral prolongation strategy. Solutions computed at coarser resolutions are prolonged in modal (Chebyshev) space via zero-padding, providing accurate initial guesses for finer-level solves. Since the high-frequency Chebyshev coefficients decay rapidly for smooth solutions, this results in minimal corrections at each refinement step. The refinement process self-terminates based on the convergence of the normed difference between neighboring discretization levels.

Practically, this refinement mechanism not only accelerates convergence and avoids preconditioners but also makes high-resolution (over $20$ million DoF) solutions feasible, with computational scaling preserved at $O(N^d \log N)$ in $d$ dimensions.

Benchmark Results and Analysis

1D and 3D Linear Diffusion Problems

The solver's convergence and accuracy are validated for the Poisson equation with analytic solutions. Comparative analysis with established DCT/DST-based methods reveals that the Chebyshev collocation approach achieves exponential convergence, with errors rapidly reaching machine precision (order $10^{-14}$–$10^{-16}$) using moderate $N$ ($N=64$). By contrast, sine/cosine transform solvers demonstrate only algebraic decay unless the solution and boundary data are compatible with periodic (or appropriate even/odd) extensions. Therefore, for truly non-periodic smooth fields, Chebyshev collocation is strongly superior in convergence rate.

A 3D Poisson problem with mixed boundary conditions demonstrates the method's scalability and robustness, with errors below $10^{-10}$ attained at $N=256$. The computational cost matches theoretical FFT scaling ($20$0), and the number of solver iterations remains almost constant across discretization levels, a direct consequence of the hierarchical refinement.

In the absence of hierarchical refinement, convergence stagnates due to operator ill-conditioning at high $20$1, making the refinement essential for practical large-scale applications.

Heterogeneous Diffusion: Homogenization Benchmarks

In domains with strong material heterogeneity (e.g., spherical inclusion problems), the matrix-free Chebyshev solver delivers solution accuracy and homogenized effective properties matching those of periodic FFT-based reference solvers. Accuracy with respect to the Maxwell model for effective conductivity is maintained across contrasts varying from $20$2 to $20$3, with relative errors remaining below $20$4 for most contrasts and robust LGMRES convergence. The non-uniform distribution of CGL nodes is shown to incur only minor discrepancies in interface representation for coarsest grids; as $20$5 increases, agreement with regular-grid FFT solvers becomes exact.

Nonlinear Diffusion Problems

The method is extended to fully nonlinear diffusion (nonlinear thermal equilibrium) governed by temperature-dependent conductivity. Manufactured analytic solutions demonstrate that the method retains exponential convergence, achieving errors below $20$6 for $20$7. Computational cost scaling remains $20$8, confirming no loss in efficiency for nonlinear operator handling. Hierarchical refinement is found to reduce both the number of Newton iterations and total Krylov steps for large $20$9.

Implications and Directions for Future Research

The methodology generalizes FFT-based PDE solvers beyond the long-standing restriction to periodic or semi-periodic settings. Key implications include:

• Practical applicability to arbitrary BCs: The method natively handles fully general, spatially inhomogeneous Dirichlet/Neumann BCs without domain tricks, auxiliary extensions, or BC-adapted expansions.
• Compatibility with nonlinear constitutive physics: Matrix-free Chebyshev collocation readily incorporates strong nonlinearity, supporting advanced material behavior in computational micromechanics, heat transfer, and more.
• Scalability: Hierarchical refinement systematically overcomes spectral ill-conditioning, permitting resolution of high-fidelity 3D fields at practical computational cost.

Potential future extensions include preconditioned or multilevel approaches for further acceleration, hybridization with domain decomposition for complex geometries, or adaptation to time-dependent and coupled PDE systems within this matrix-free Chebyshev-FFT paradigm. Direct integration with high-performance computing architectures appears straightforward due to the reliance on established FFT libraries.

Conclusion

This work establishes a robust and efficient matrix-free Chebyshev collocation solver, leveraging FFT-based operators, for stationary diffusion-type problems under arbitrary boundary conditions. The method demonstrates exponential convergence for smooth, non-periodic solutions and achieves computational performance comparable to classic FFT-based periodic solvers. The introduction of hierarchical modal refinement resolves inherent conditioning limitations and is essential for large-scale applications. The framework is rigorously validated for linear and nonlinear diffusion as well as heterogeneous property fields and can serve as a foundation for extending FFT-based spectral methods across a wide class of BVPs prevalent in scientific and engineering computational modeling.