---
title: Chorin's Projection Method
url: https://www.emergentmind.com/topics/chorin-s-projection-method
type: topic
---

# Chorin's Projection Method

Searching arXiv for the cited projection-method convergence papers to ground the article in current arXiv records.
Chorin's projection method is a fractional-step procedure for the incompressible Navier–Stokes equations in which the momentum update and the incompressibility constraint are separated into distinct subproblems. In the formulation emphasized by recent convergence analyses, one first computes a tentative velocity that is generally not divergence-free, then enforces incompressibility by a projection onto a divergence-free subspace, with the pressure appearing as the potential part of the correction. The method was first studied by Chorin in the framework of a finite difference method and by Temam in the framework of a finite element method; contemporary arXiv analyses establish convergence, up to subsequences, to Leray–Hopf weak solutions under low regularity assumptions and make precise the role of discrete energy estimates, compactness, and mesh–time scaling in both finite-difference and conforming finite-element settings [1809.04383], [2508.13416], [2009.05229].

## 1. Historical development and continuous problem

The incompressible Navier–Stokes problem considered in the modern convergence literature is posed on a bounded Lipschitz domain $\Omega \subset \mathbb{R}^3$ or, in the finite-element treatment, $\Omega \subset \mathbb{R}^d$ with $d=2,3$, over a time interval $(0,T)$. In the notation of Kuroki–Soga, one seeks $(u,p)$ satisfying
\[
\partial_t u + (u\cdot \nabla)u - \nu \Delta u + \nabla p = f,
\qquad
\nabla \cdot u = 0,
\qquad
u|_{\partial \Omega}=0,
\]
together with divergence-free initial data $u_0 \in L^2(\Omega)^3$ and external force $f \in L^2$ [1809.04383]. Weber formulates the same structure with viscosity parameter $\mu>0$ and the standard no-slip boundary condition [2508.13416].

The weak-solution concept used in these analyses is the Leray–Hopf notion. In the finite-difference treatment, a Leray–Hopf weak solution is a vector field
\[
u \in L^\infty(0,T;L^2(\Omega)^3)\cap L^2(0,T;H^1_{0,\sigma}(\Omega)^3)
\]
with $\nabla\cdot u=0$ and satisfying the standard weak form against divergence-free test functions [1809.04383]. In the finite-element setting, the definition is sharpened by including
\[
\partial_t u \in L^{4/3}(0,T;(H^1_{div})^*)
\]
and the energy inequality
\[
\frac12 \|u(t)\|_{L^2}^2 + \mu \int_0^t \|\nabla u(s)\|_{L^2}^2\,ds
\le
\frac12 \|u_0\|_{L^2}^2 + \int_0^t (f,u)\,ds
\]
[2508.13416].

Historically, Chorin's original work addressed periodic boundary conditions and convergence under smoothness assumptions, while Temam gave an abstract route to Leray–Hopf weak solutions on bounded domains with no-slip boundary conditions. The later finite-difference analysis extends Chorin's result with full details to arbitrary bounded Lipschitz domains in $\mathbb{R}^3$, no-slip boundary conditions, and external forcing [1809.04383].

## 2. Fractional-step structure of the method

The defining feature of Chorin's method is the decomposition of one incompressible update into a prediction stage and a projection stage. In the finite-difference scheme of Kuroki–Soga, given a discrete divergence-free $u^n$ vanishing on the boundary, one computes an intermediate velocity $u^{n+\frac12}$ from
\[
\frac{u^{n+\frac12}-u^n}{\tau}
-\nu \Delta_h u^{n+\frac12}
+\text{convective term}^n
=
f^{n+1}
\]
on interior grid points, with $u^{n+\frac12}=0$ on the discrete boundary [1809.04383]. The pressure is then obtained from the discrete Poisson problem
\[
\Delta_h p^{n+1}
=
\frac1{\tau}\nabla_h\cdot u^{n+\frac12},
\qquad
p^{n+1}=0 \text{ on } \partial \Omega_h,
\]
and the divergence-free update is
\[
u^{n+1}=u^{n+\frac12}-\tau \nabla_h p^{n+1}.
\]
Equivalently, the method applies the discrete Helmholtz–Hodge operator $P_h$ and writes $u^{n+1}=P_h u^{n+\frac12}$ [1809.04383].

A closely related fully discrete formulation is analyzed in the 2020 continuation paper. There, on a uniform Cartesian grid, one solves implicitly for $u^{n+1}$ with a centered convective discretization, then projects by
\[
u^{n+1}\leftarrow P_h u^{n+1}.
\]
The pressure is not advanced as an explicit independent variable; rather, it is hidden in the discrete Helmholtz–Hodge decomposition $u=P_hu + D\phi$ [2009.05229].

In the conforming finite-element setting, the same conceptual split is expressed in variational form. The prediction step computes an intermediate velocity field $\widetilde u_h^{m+1}$ that is generally not divergence-free:
\[
\frac{\widetilde u_h^{m+1}-u_h^m}{\Delta t}
+
b(u_h^m,\widetilde u_h^{m+1},v)
+
\mu(\nabla \widetilde u_h^{m+1},\nabla v)
=
(P_h f^m,v)
\]
for all $v \in U_h$. The projection step then finds $(u_h^{m+1},p_h^{m+1})\in U_h\times P_h$ such that
\[
\frac{u_h^{m+1}-\widetilde u_h^{m+1}}{\Delta t}
=
\nabla p_h^{m+1}
\quad\text{in weak form,}
\qquad
\operatorname{div} u_h^{m+1}=0
\]
[2508.13416]. Weber explicitly distinguishes this original Chorin form from incremental variants in which $\nabla p^m$ is added in the prediction and the pressure increment is solved for in the projection [2508.13416].

## 3. Discrete realizations and operators

The recent literature treats two principal discretization frameworks.

| Aspect | Finite difference | Conforming finite element |
|---|---|---|
| Spatial setting | $\Omega_h \subset h\mathbb{Z}^3$ on a uniform Cartesian grid | Shape-regular simplicial meshes $\{\mathcal{T}_h\}$ |
| Incompressibility enforcement | Discrete Helmholtz–Hodge projector $P_h$ | Mixed projection step in $U_h\times P_h$ |
| Pressure representation | Explicit Poisson solve or potential part in $u=P_hu+D\phi$ | Pressure variable $p_h^{m+1}$ in the Darcy-form projection |

In the finite-difference analysis of Kuroki–Soga, the mesh is
\[
\Omega_h=\{x\in h\mathbb{Z}^3:\text{ the closed cube of side } h \text{ centered at } x \text{ lies in } \Omega\},
\]
with interior points $\Omega_h\setminus \partial \Omega_h$. The discrete first derivatives are
\[
D_i^+\phi(x)=\frac{\phi(x+h e_i)-\phi(x)}{h},
\qquad
D_i^-\phi(x)=\frac{\phi(x)-\phi(x-he_i)}{h},
\]
the discrete gradient is $\nabla_h\phi=(D_1^+\phi,D_2^+\phi,D_3^+\phi)$, the discrete divergence is $\nabla_h\cdot w=D_1^-w_1+D_2^-w_2+D_3^-w_3$, and the discrete Laplacian is $\Delta_h=\sum_i D_i^+D_i^-$ [1809.04383]. The 2020 paper uses a closely related uniform-grid discretization, with centered discrete derivatives $D_i=\tfrac12(D_i^+ + D_i^-)$ and no staggering: each velocity component and pressure live on the same grid [2009.05229].

The discrete nonlinear term in the 2018 implicit scheme is a Ladyzhenskaya-type form,
\[
\text{convective term}_i^n(x)
=
\sum_{j=1}^3
\bigl[
u_j^n(x-he_j)D_j^+u_i^n(x-he_j)
+
u_j^n(x)D_j^-u_i^{n+\frac12}(x)
\bigr],
\]
chosen so that the discrete energy argument closes under the divergence-free constraint [1809.04383].

In the conforming finite-element analysis, the spaces satisfy
\[
U_h\subset H_0^1(\Omega)^n,
\qquad
P_h\subset L_0^2(\Omega),
\]
together with the usual approximation properties of the $L^2$-projection and the discrete inf–sup condition
\[
\inf_{q\in P_h}\sup_{v\in U_h}
\frac{(\operatorname{div}v,q)}{\|v\|_{H^1}\|q\|_{L^2}}
\ge \beta^* >0.
\]
The convection term is written using the skew-symmetric trilinear form
\[
b(w,u,v)=\int_\Omega (w\cdot \nabla)u\cdot v
+\frac12 (\operatorname{div}w)(u\cdot v)\,dx
\]
[2508.13416]. The implementation-oriented discussion explicitly lists conforming LBB-stable pairs such as P2–P1 Taylor–Hood, MINI, and $P_k$–$P_{k-1}$ for $k\ge 2$ [2508.13416].

## 4. Solvability, stability, and a priori estimates

A central result of the finite-difference convergence theory is unconditional solvability of the implicit scheme. At each time step, the combined system defined by the tentative velocity equation, the pressure Poisson problem, and the projection is equivalent to a linear system $Ay=b$, and a discrete energy argument shows that $A$ is invertible for every $h,\tau>0$; no restriction on the time step enters the solvability argument [1809.04383]. The 2020 paper states the same conclusion for its fully discrete formulation: for each fixed $h,\tau>0$ and each divergence-free $u^n$ vanishing on the boundary, the implicit step for $u^{n+1}$ has a unique solution [2009.05229].

The associated discrete energy inequalities furnish uniform bounds. In the 2018 analysis one obtains
\[
\|u^{n+1}\|_{\ell_h^2}^2
+
2\nu \tau \sum_{j=1}^3 \|D_j^+u^{n+\frac12}\|_{\ell_h^2}^2
\le
\|u^n\|_{\ell_h^2}^2
+
\tau \|f^{n+1}\|_{\ell_h^2}^2,
\]
and summation yields uniform bounds in $\ell_t^\infty \ell_x^2$ and $\ell_t^2 H_h^1$ [1809.04383]. The 2020 continuation derives related estimates,
\[
\|u^{n+1}\|_{\ell_h^2}\le \|u^n\|_{\ell_h^2}+\tau \|f^n\|_{\ell_h^2},
\]
together with a discrete $H^1$-type bound obtained by summation by parts [2009.05229].

For the conforming finite-element discretization, testing the prediction by $v=\widetilde u_h^{m+1}$ and exploiting
\[
b(u_h^m,\widetilde u_h^{m+1},\widetilde u_h^{m+1})=0
\]
gives the discrete energy identity
\[
\frac12\|u_h^M\|^2
+\frac12\sum_{m=0}^{M-1}\|u_h^{m+1}-u_h^m\|^2
+\frac12\sum_{m=0}^{M-1}\|u_h^{m+1}-\widetilde u_h^{m+1}\|^2
+\mu \Delta t\sum_{m=0}^{M-1}\|\nabla \widetilde u_h^{m+1}\|^2
=
\frac12\|u_h^0\|^2
+\Delta t \sum_{m=0}^{M-1}(f^m,\widetilde u_h^{m+1}),
\]
from which one obtains uniform bounds
\[
u_h,\widetilde u_h \in L^\infty(0,T;L^2),
\qquad
\widetilde u_h \in L^2(0,T;H_0^1),
\]
together with a bound on the discrete time derivative in a dual space [2508.13416].

A recurring point of clarification is that unconditional solvability and basic energy stability are stronger statements than unconditional strong convergence. In the finite-difference paper, the authors explicitly impose no CFL-type upper bound on $\tau$, and the lower bound
\[
h^{3-\alpha}\le \tau,\qquad \alpha\in(0,2],
\]
is used only in the strong-convergence proof, not in solvability or basic energy estimates [1809.04383].

## 5. Compactness mechanisms and convergence theory

The principal convergence result of the 2018 finite-difference analysis is Theorem 7.1: under the scaling condition
\[
h^{3-\alpha}\le \tau,
\qquad
\alpha\in(0,2],
\]
the piecewise-constant-in-time-and-space reconstructions of the discrete velocity converge, up to a subsequence, strongly in $L^2(0,T;L^2(\Omega)^3)$ to a limit $v$, and that limit is a Leray–Hopf weak solution of the incompressible Navier–Stokes equations with no-slip boundary condition and external force [1809.04383]. The argument proceeds by unconditional solvability, discrete energy bounds, weak compactness, a new interpolation-type inequality for step functions, and passage to the limit in a weak formulation of the fully discrete scheme [1809.04383].

The crucial analytical issue is the convective term. Weak convergence is insufficient for passing to the nonlinear term in the standard way, so strong convergence in $L^2$ is required. In Kuroki–Soga, this is obtained by a compactness method based on a new interpolation inequality for step functions, described as a discrete analog of Aubin–Lions and tied to the structure of the projection scheme and the lower-bound scaling relation [1809.04383].

The 2025 finite-element analysis isolates an analogous but structurally different difficulty: the intermediate velocity $\widetilde u_h$ has good $H^1$ bounds but no time continuity, whereas the divergence-free velocity $u_h$ has time continuity only in a very weak dual norm and only $L^2$ bounds in space. Weber therefore proves a modified Aubin–Lions lemma. If sequences $\{w_n\}$ and $\{\widehat w_n\}$ satisfy
\[
w_n-\widehat w_n \to 0 \text{ in } L^2((0,T)\times \Omega),
\]
with
\[
w_n\in L^2(0,T;L^2_{div}),
\quad
\partial_t w_n\in L^2(0,T;(H^2\cap H^1_{div})^*),
\quad
\widehat w_n\in L^2(0,T;H_0^1),
\]
then both are precompact in $L^2((0,T)\times \Omega)$ [2508.13416]. Combined with the discrete energy inequality, this yields
\[
u_h,\bar u_h,\widetilde u_h \to u \text{ strongly in } L^2((0,T)\times \Omega),
\]
while $u_h \rightharpoonup u$ in $L^2(0,T;H^1)$ and $\partial_t u_h \rightharpoonup \partial_t u$ in $L^{4/3}(0,T;(H^1)^*)$, from which one recovers the Leray–Hopf weak formulation and the continuous energy inequality [2508.13416].

The 2020 continuation paper extends the finite-difference program in two directions. First, on fixed time intervals it obtains weak convergence in $L^2([0,T];H^1)$ and strong convergence in $L^2([0,T];L^2)$ under the condition
\[
h^{4-\alpha}\le \tau,
\qquad
0<\alpha\le 2,
\]
with the limit again a Leray–Hopf weak solution [2009.05229]. Second, under the uniform forcing bound
\[
\|f\|_{L^2([n-1,n];L^2(\Omega))}\le C \quad \text{uniformly in } n,
\]
the paper proves time-global solvability and obtains a time-global Leray–Hopf solution in the limit by a diagonal subsequence argument [2009.05229].

## 6. Error estimates, periodic solutions, and practical variants

The rough-solution convergence theory is primarily existential and compactness-based. Kuroki–Soga explicitly state that they do not derive quantitative error bounds or rates in terms of $h,\tau$ for approximation to a weak solution; the focus is convergence to a Leray–Hopf solution without any further regularity assumptions [1809.04383].

Quantitative error estimates appear in the 2020 continuation for smooth exact solutions. If the exact solution $(v,p)$ belongs to
\[
C^3([0,T]\times \overline{\Omega})\times C^2
\]
and
\[
\tau = O(h^{3/4}),
\]
then
\[
\max_{0\le n\le T/\tau}\|u^n-v(t^n)\|_{\ell_h^2}
+
\max_{0\le n\le T/\tau}\|\widetilde u^n-v(t^n)\|_{\ell_h^2}
=
O(h^{1/4}).
\]
The paper attributes the reduced rate to the boundary layer error in the projection step,
\[
\|v^n-P_h v^n\|_{\ell_h^2}=O(h),
\]
which forces the sup-error to be only $O(h^{1/4})$ under the stated time-step scaling [2009.05229]. Under periodic boundary conditions and higher regularity, the same paper reports recovery of the classical $O(h^2)$ error and a uniform $O(\sqrt h)$ sup-norm estimate under $\tau=O(h^2)$ [2009.05229].

The 2020 paper also analyzes time-periodic forcing. If $f(t+1)=f(t)$, then the discrete time-$1$ map leaves a large $\ell_h^2$-ball invariant, and Brouwer's fixed-point theorem yields a time-periodic discrete solution of period $T_1=1/\tau$; as the discretization parameters vanish, any such discrete periodic solution converges, up to a subsequence, to a time-periodic Leray–Hopf weak solution of period $1$ [2009.05229]. Under a smallness condition,
\[
\|\widetilde u^n\|_\infty \le \beta_0 < A^{-1}/4,
\]
a discrete periodic solution is unique in that class and every other solution converges exponentially fast to it in $\ell_h^2$; the limiting argument yields exponential stability for small Leray–Hopf weak solutions as well [2009.05229].

From an implementation viewpoint, the finite-difference algorithm consists of initialization by cell averages and projection, followed at each time step by a tentative solve,
\[
[I-\nu \tau \Delta_h]u^{n+\frac12} + \tau N_h(u^n,u^{n+\frac12})
=
u^n+\tau f^{n+1},
\]
a pressure Poisson solve, and a projection update [1809.04383]. Each step requires inversion of a sparse, symmetric positive definite matrix for the operators $[I-\nu\tau \Delta_h]$ and $\Delta_h$, along with gradient and divergence evaluations [1809.04383]. In the finite-element setting, the prediction step assembles the mass matrix, the skew-symmetric convection term, and the stiffness matrix, while the projection step may be written in Darcy form and reduced by elimination to a symmetric positive definite pressure-Poisson system, followed by a cheap local update for the velocity [2508.13416]. Weber also notes that one may replace the Darcy-form projection by the classical Poisson-projection form
\[
\Delta \phi = \operatorname{div}u^*,
\qquad
u^{m+1}=u^*-\Delta t \nabla \phi,
\]
with only minor changes in the convergence proof [2508.13416].

## 7. Interpretive issues in the modern literature

Several points recur across recent analyses. First, the phrase “unconditional” does not denote a fully scale-free strong-convergence theory in every formulation. In the finite-difference work of 2018, unconditional solvability holds for every $h,\tau>0$, while strong $L^2$ convergence requires the compactness-forcing relation $h^{3-\alpha}\le \tau$ [1809.04383]. In the finite-element proof, the method is described as unconditionally energy-stable, but the convergence proof assumes
\[
h \lesssim \sqrt{\Delta t}
\]
so that pressure-error terms from the projection remain $o(1)$ [2508.13416].

Second, the target of convergence is a Leray–Hopf weak solution, generally only up to a subsequence. This is a deliberately low-regularity framework: the convergence proofs do not assume more than square-integrable initial data in the finite-element analysis, and the 2018 finite-difference paper explicitly avoids further regularity assumptions [2508.13416], [1809.04383]. A plausible implication is that these results prioritize compactness and consistency over uniqueness, which is unavailable at the Leray–Hopf level in three dimensions.

Third, the role of pressure depends on the formulation. In projection methods written as Poisson correction schemes, pressure appears as the scalar potential generating the divergence-removing update. In the discrete Helmholtz–Hodge formulation of the 2020 finite-difference paper, pressure is not carried as a primary unknown in the time-marching formula but is hidden in the decomposition
\[
u=P_hu + D\phi,
\]
with $\phi$ determined by a discrete Poisson problem and a zero-mean constraint on parity classes [2009.05229].

Taken together, these results place Chorin's projection method at the intersection of operator splitting, incompressibility enforcement, and compactness-based weak convergence theory. The method's numerical appeal lies in the decoupling of momentum and incompressibility, while its modern mathematical analysis turns on discrete energy identities, projection structure, and nonstandard compactness arguments capable of handling the mismatch between regularity of tentative and projected velocities [1809.04383], [2508.13416], [2009.05229].

Source: https://www.emergentmind.com/topics/chorin-s-projection-method