A high-order, pressure-robust, and decoupled finite difference method for the Stokes problem

In this paper, we consider the Stokes problem with Dirichlet boundary conditions and the constant kinematic viscosity ν in an axis-aligned domain Ω. We decouple the velocity u and pressure p by deriving a novel biharmonic equation in Ω and third-order boundary conditions on ∂Ω. In contrast to the fourth-order streamfunction approach, our formulation does not require Ω to be simply connected. For smooth velocity fields u in two dimensions, we explicitly construct a finite difference method (FDM) with sixth-order consistency to approximate u at all relevant grid points: interior points, boundary side points, and boundary corner points. The resulting scheme yields two linear systems A₁u_h⁽¹⁾=b₁ and A₂u_h⁽²⁾=b₂, where A₁,A₂ are constant matrices, and b₁,b₂ are independent of the pressure p and the kinematic viscosity ν. Thus, the proposed method is pressure- and viscosity-robust. To accommodate velocity fields with less regularity, we modify the FDM by removing singular terms in the right-hand side vectors. Once the discrete velocity is computed, we apply a sixth-order finite difference operator to first approximate the pressure gradient locally, and then calculate the pressure itself locally with sixth-order accuracy, both without solving any additional linear systems. In our numerical experiments, we test both smooth and non-smooth solutions (u,p) in a square domain, a triply connected domain, and an L-shaped domain in two dimensions. The results confirm sixth-order convergence of the velocity, pressure gradient, and pressure in the ℓ_∞-norm for smooth solutions. For non-smooth velocity fields, our method achieves the expected lower-order convergence. Moreover, the observed velocity error ‖uh−u‖_∞ is independent of the pressure p and viscosity ν.