Quadratic Bubble Upwinding Method

• Quadratic Bubble Upwinding Method is a Petrov–Galerkin finite element scheme that introduces quadratic bubble functions into the test space to stabilize convection–diffusion discretizations.
• It converts part of the convection term into effective artificial diffusion, bridging finite element upwinding and finite-difference schemes through explicit matrix identities.
• The method achieves optimal convergence and improved stability, extends to multidimensional problems, and effectively suppresses nonphysical oscillations in convection-dominated regimes.

Searching arXiv for papers on the quadratic bubble upwinding Petrov–Galerkin method and closely related convection–diffusion discretizations. The Quadratic Bubble Upwinding Method is a Petrov–Galerkin finite element discretization for convection–diffusion problems in which the standard continuous piecewise linear trial space is paired with a bubble-modified test space. In the treatment developed for the one-dimensional singularly perturbed model problem

$$-\varepsilon u''(x)+u'(x)=f(x), \qquad 0<x<1,\qquad u(0)=u(1)=0,$$

the method uses local quadratic bubbles of the form $B(x)=\beta x(h-x)$ or, equivalently under the paper’s scaling, $B^q(x)=\frac{4\beta}{h^2}x(h-x)$, and thereby converts part of the convection contribution into an effective artificial diffusion term (Bacuta et al., 2024, Bacuta et al., 2024, Bacuta, 4 Sep 2025). Across the 2024–2025 analysis, the method is presented not only as a stabilized finite element scheme, but also as an algebraic bridge between finite element upwinding and finite-difference artificial-diffusion schemes, with explicit matrix identities, discrete norm estimates, and a later convergence theory based on Green-matrix representations (Bacuta et al., 2024, Bacuta, 4 Sep 2025).

1. Model problem and variational framework

The foundational setting is a one-dimensional singularly perturbed convection–diffusion boundary value problem on $(0,1)$,

$$-\varepsilon u''(x)+k\,u'(x)=f(x),\qquad u(0)=u(1)=0,$$

with $\varepsilon>0$, $k>0$, and, without loss of generality, $k=1$, so that the model becomes

$-\varepsilon u''(x)+u'(x)=f(x),\qquad u(0)=u(1)=0.$

The analysis focuses on the convection-dominated regime $\varepsilon<1$, with $B(x)=\beta x(h-x)$0 assumed continuous on $B(x)=\beta x(h-x)$1, and in later estimates sometimes $B(x)=\beta x(h-x)$2, $B(x)=\beta x(h-x)$3, or $B(x)=\beta x(h-x)$4 depending on the result (Bacuta et al., 2024, Bacuta et al., 2024, Bacuta, 4 Sep 2025).

The standard weak form is posed on $B(x)=\beta x(h-x)$5: $B(x)=\beta x(h-x)$6 where

$B(x)=\beta x(h-x)$7

A central theme of the 2024 comparison paper is that different discretizations are best understood through the optimal trial norms induced by their test spaces; within that framework, the bubble upwinding Petrov–Galerkin method is singled out as the most performant discretization for the one-dimensional model (Bacuta et al., 2024).

On a uniform mesh $B(x)=\beta x(h-x)$8, $B(x)=\beta x(h-x)$9, $B^q(x)=\frac{4\beta}{h^2}x(h-x)$0, the trial space is the standard continuous piecewise linear finite element space

$B^q(x)=\frac{4\beta}{h^2}x(h-x)$1

or equivalently $B^q(x)=\frac{4\beta}{h^2}x(h-x)$2, depending on notation (Bacuta et al., 2024, Bacuta et al., 2024, Bacuta, 4 Sep 2025). The distinguishing feature of the quadratic bubble method lies entirely in the choice of the test space.

2. Bubble-modified test space and Petrov–Galerkin structure

The bubble construction begins from a generating local bubble $B^q(x)=\frac{4\beta}{h^2}x(h-x)$3 satisfying

$B^q(x)=\frac{4\beta}{h^2}x(h-x)$4

together with a mean-value identity of the form

$B^q(x)=\frac{4\beta}{h^2}x(h-x)$5

or, in the later notation,

$B^q(x)=\frac{4\beta}{h^2}x(h-x)$6

By translation, one obtains local bubbles $B^q(x)=\frac{4\beta}{h^2}x(h-x)$7 supported on $B^q(x)=\frac{4\beta}{h^2}x(h-x)$8 (Bacuta et al., 2024, Bacuta, 4 Sep 2025).

The bubble-modified test space is then defined by basis functions

$B^q(x)=\frac{4\beta}{h^2}x(h-x)$9

so that

$(0,1)$0

The discrete Petrov–Galerkin problem is: $(0,1)$1 Because $(0,1)$2, the method is genuinely Petrov–Galerkin; because the bubble correction is designed to bias the test functions in the convection direction, it is described as upwinding (Bacuta et al., 2024, Bacuta et al., 2024, Bacuta, 4 Sep 2025).

A generic test function $(0,1)$3 can be decomposed as

$(0,1)$4

with $(0,1)$5. The bubble identities then yield the key relations

$(0,1)$6

and, in the 2024 comparison paper’s presentation,

$(0,1)$7

In the 2025 formulation, the bilinear form reduces to

$(0,1)$8

In the 2024 formulation, the corresponding reformulation is expressed as

$(0,1)$9

The common conclusion is that bubble modification produces an effective artificial diffusion

$$-\varepsilon u''(x)+k\,u'(x)=f(x),\qquad u(0)=u(1)=0,$$0

which is the algebraic mechanism behind the upwinding effect (Bacuta et al., 2024, Bacuta et al., 2024, Bacuta, 4 Sep 2025).

3. Quadratic bubble specialization

In the quadratic-bubble specialization, the generating bubble on $$-\varepsilon u''(x)+k\,u'(x)=f(x),\qquad u(0)=u(1)=0,$$1 is chosen as

$$-\varepsilon u''(x)+k\,u'(x)=f(x),\qquad u(0)=u(1)=0,$$2

or, in the later scaled form,

$$-\varepsilon u''(x)+k\,u'(x)=f(x),\qquad u(0)=u(1)=0,$$3

The mean-value computation gives

$$-\varepsilon u''(x)+k\,u'(x)=f(x),\qquad u(0)=u(1)=0,$$4

so that the 2024 connection paper summarizes the resulting constant as

$$-\varepsilon u''(x)+k\,u'(x)=f(x),\qquad u(0)=u(1)=0,$$5

and therefore

$$-\varepsilon u''(x)+k\,u'(x)=f(x),\qquad u(0)=u(1)=0,$$6

In the 2025 convergence paper, the corresponding average is written as

$$-\varepsilon u''(x)+k\,u'(x)=f(x),\qquad u(0)=u(1)=0,$$7

hence again

$$-\varepsilon u''(x)+k\,u'(x)=f(x),\qquad u(0)=u(1)=0,$$8

Accordingly, the quadratic bubble test space is

$$-\varepsilon u''(x)+k\,u'(x)=f(x),\qquad u(0)=u(1)=0,$$9

or equivalently

$\varepsilon>0$0

This construction is precisely what the 2024 synthesis identifies as the Quadratic Bubble Upwinding Method (Bacuta et al., 2024), with the 2024 comparison and 2025 convergence papers embedding it into a broader Petrov–Galerkin framework (Bacuta et al., 2024, Bacuta, 4 Sep 2025).

One important parameter choice is $\varepsilon>0$1. In that case,

$\varepsilon>0$2

which coincides with the standard upwind artificial-diffusion correction in the finite-difference formulation (Bacuta et al., 2024). The 2024 comparison paper also presents a scaling in which the quadratic bubble constants satisfy $\varepsilon>0$3 and $\varepsilon>0$4, again identifying the same stabilized coefficient (Bacuta et al., 2024). This suggests that the specific normalization used in the papers differs, while the operative stabilized diffusion $\varepsilon>0$5 is the common practical reference point.

4. Algebraic equivalence with finite-difference upwinding

A major result of the 2024 connection paper is that the bubble-Petrov–Galerkin finite element method and a corresponding finite-difference upwinding scheme can be written so that they have the same stiffness matrix (Bacuta et al., 2024). On the finite-difference side, the standard upwind scheme for $\varepsilon>0$6 uses a backward difference because $\varepsilon>0$7. It can be expressed in an artificial-diffusion form with

$\varepsilon>0$8

and the standard upwind choice is

$\varepsilon>0$9

yielding

$k>0$0

For the quadratic bubble $k>0$1, the corresponding finite-difference artificial diffusion is matched by

$k>0$2

and when $k>0$3, the method coincides with standard upwinding (Bacuta et al., 2024).

On the finite element side, the bubble-PG stiffness matrix is

$k>0$4

under the paper’s scaling conventions, with

$k>0$5

The paper then chooses the bubble so that

$k>0$6

which implies

$k>0$7

The finite-difference and bubble-PG finite element methods thus share the same stiffness matrix; the distinction lies only in the right-hand side (Bacuta et al., 2024).

The 2025 convergence paper gives the same matrix mechanism in a more explicit finite element form. The one-dimensional PG system is

$k>0$8

with

$k>0$9

hence

$k=1$0

For the quadratic bubble,

$k=1$1

A notable structural point, stated explicitly in that paper, is that the matrix depends only on $k=1$2, $k=1$3, and the mean value $k=1$4 of the bubble (Bacuta, 4 Sep 2025).

The right-hand sides remain different. The finite element side uses

$k=1$5

whereas the finite-difference side uses pointwise values $k=1$6 or quadrature approximations thereof. The 2024 connection paper interprets the FD method as the FE bubble-PG method with the dual vector approximated by a quadrature rule, typically the composite trapezoidal rule, and notes that better quadratures can improve the FD variant, including a Cavalieri–Simpson FD method (Bacuta et al., 2024).

5. Stability, optimal norms, and approximation properties

The 2024 comparison paper develops the analysis in the language of continuous and discrete optimal trial norms (Bacuta et al., 2024). For the continuous problem,

$k=1$7

where $k=1$8 is defined by

$k=1$9

In one dimension with $-\varepsilon u''(x)+u'(x)=f(x),\qquad u(0)=u(1)=0.$0, the paper derives

$-\varepsilon u''(x)+u'(x)=f(x),\qquad u(0)=u(1)=0.$1

hence

$-\varepsilon u''(x)+u'(x)=f(x),\qquad u(0)=u(1)=0.$2

so that

$-\varepsilon u''(x)+u'(x)=f(x),\qquad u(0)=u(1)=0.$3

The stated purpose of this framework is to compare how different test spaces alter the norm in which the discrete method is quasi-optimal (Bacuta et al., 2024).

For the discrete upwinding Petrov–Galerkin setting, the same paper gives a discrete optimal norm on $-\varepsilon u''(x)+u'(x)=f(x),\qquad u(0)=u(1)=0.$4,

$-\varepsilon u''(x)+u'(x)=f(x),\qquad u(0)=u(1)=0.$5

with $-\varepsilon u''(x)+u'(x)=f(x),\qquad u(0)=u(1)=0.$6 determined by bubble moment identities. The 2024 connection paper expresses a related norm as

$-\varepsilon u''(x)+u'(x)=f(x),\qquad u(0)=u(1)=0.$7

In both presentations, the essential point is that bubble enrichment changes the discrete norm so that the convection-dominated term is controlled more favorably (Bacuta et al., 2024, Bacuta et al., 2024).

The resulting quasi-optimal estimate for the general UPG method is

$-\varepsilon u''(x)+u'(x)=f(x),\qquad u(0)=u(1)=0.$8

under the mesh condition

$-\varepsilon u''(x)+u'(x)=f(x),\qquad u(0)=u(1)=0.$9

For the quadratic bubble case, the comparison paper reports the explicit estimate

$\varepsilon<1$0

along with a rescaled mixed $\varepsilon<1$1/energy version (Bacuta et al., 2024).

The same paper contrasts this with the standard $\varepsilon<1$2 Galerkin method, which is stable only in a weaker optimal norm and may exhibit nonphysical oscillations when $\varepsilon<1$3, and with the SPLS discretization, which is well posed and optimally accurate in its norm but can still oscillate near the boundary for data with nonzero mean forcing (Bacuta et al., 2024). Within that comparison, the bubble upwinding Petrov–Galerkin method is concluded to be the most performant discretization for the one-dimensional model (Bacuta et al., 2024).

6. Exponential bubbles, special quadratic scaling, and convergence results

Although the topic is specifically quadratic bubbles, the modern understanding of the method in these papers depends heavily on comparison with the exponential bubble discretization. The exponential bubble is defined by

$\varepsilon<1$4

In the 2024 connection paper, it produces

$\varepsilon<1$5

leading to the classical Il’in–Allen–Southwell / Scharfetter–Gummel discretization (Bacuta et al., 2024). The same paper proves the exact recovery theorem: if $\varepsilon<1$6 is the FE solution of the bubble-PG method with exponential test space, then

$\varepsilon<1$7

equivalently $\varepsilon<1$8, the nodal interpolant of the exact solution (Bacuta et al., 2024). The 2025 convergence paper restates the exact nodal property as

$\varepsilon<1$9

and further proves that the inverse matrix of the exponential-bubble discretization is exactly the Green matrix: $B(x)=\beta x(h-x)$00 That exact inverse becomes the algebraic backbone for the quadratic-bubble convergence proof (Bacuta, 4 Sep 2025).

The 2025 innovation is a special choice of the quadratic parameter $B(x)=\beta x(h-x)$01 so that the quadratic bubble has the same average as the exponential bubble: $B(x)=\beta x(h-x)$02 that is,

$B(x)=\beta x(h-x)$03

With this “special scaling parameter,”

$B(x)=\beta x(h-x)$04

The paper motivates this construction by noting that the exponential bubble is computationally inconvenient and can be numerically fragile when $B(x)=\beta x(h-x)$05 underflows or becomes indistinguishable from zero in floating-point arithmetic; the quadratic bubble avoids exponential functions while preserving the same discrete operator (Bacuta, 4 Sep 2025).

Under the assumptions $B(x)=\beta x(h-x)$06 and

$B(x)=\beta x(h-x)$07

the paper proves the discrete infinity-norm estimate

$B(x)=\beta x(h-x)$08

A corollary states that if $B(x)=\beta x(h-x)$09, then

$B(x)=\beta x(h-x)$10

Using the decomposition $B(x)=\beta x(h-x)$11, together with the inverse-type inequality stated in the paper, the authors then derive

$B(x)=\beta x(h-x)$12

and

$B(x)=\beta x(h-x)$13

Provided the interpolant $B(x)=\beta x(h-x)$14 has standard approximation properties on the region of interest, the paper concludes that the method yields $B(x)=\beta x(h-x)$15 in $B(x)=\beta x(h-x)$16 and $B(x)=\beta x(h-x)$17 in $B(x)=\beta x(h-x)$18. An important caveat, stated explicitly, is that on domains containing boundary layers, any loss of order on the full domain comes from the interpolation of the exact layer solution on a uniform mesh, not from the UPG discretization itself (Bacuta, 4 Sep 2025).

7. Oscillation behavior, multidimensional extension, and significance

One of the principal motivations for bubble upwinding is suppression of the oscillatory behavior associated with convection-dominated discretizations. The 2024 comparison paper states that the standard linear discretization may exhibit nonphysical oscillations, and that SPLS can still oscillate near the boundary when $B(x)=\beta x(h-x)$19 (Bacuta et al., 2024). By contrast, the same paper reports that the bubble UPG method eliminates oscillations in the one-dimensional model (Bacuta et al., 2024).

For a special quadratic bubble parameter choice, the paper states that the matrix becomes bidiagonal or lower triangular, and that the discrete solution can be written explicitly so as to approximate

$B(x)=\beta x(h-x)$20

at the nodes, with a pointwise estimate of the form

$B(x)=\beta x(h-x)$21

which is used to conclude that the nodal values are close to the monotone exact profile and hence do not oscillate (Bacuta et al., 2024). This suggests that the method’s stabilization is not merely norm-theoretic but also visible in the qualitative shape of the discrete solution.

The multidimensional extension appears explicitly in the 2025 paper for

$B(x)=\beta x(h-x)$22

with convection direction $B(x)=\beta x(h-x)$23. The trial space is tensor-product $B(x)=\beta x(h-x)$24-type,

$B(x)=\beta x(h-x)$25

while the test space is bubble-modified only in the streamline direction: $B(x)=\beta x(h-x)$26 The corresponding 2D PG method is

$B(x)=\beta x(h-x)$27

For quadratic bubbles with the same special scaling as in 1D, the linear system is

$B(x)=\beta x(h-x)$28

with

$B(x)=\beta x(h-x)$29

where

$B(x)=\beta x(h-x)$30

and $B(x)=\beta x(h-x)$31 is the 1D upwind matrix (Bacuta, 4 Sep 2025).

The 2025 paper does not claim a fully sharp two-dimensional analogue of the one-dimensional theorem, but it states that the streamline stabilization mimics the 1D behavior and that optimal orders are observed numerically away from boundary layers. On subdomains excluding the boundary layer region, the reported observations are

$B(x)=\beta x(h-x)$32

and the discrete nodal error remains $B(x)=\beta x(h-x)$33 when $B(x)=\beta x(h-x)$34 (Bacuta, 4 Sep 2025). The paper summarizes the broader construction principle as using an efficient upwinding Petrov–Galerkin discretization along each streamline direction in combination with a standard discretization for the orthogonal direction or directions (Bacuta, 4 Sep 2025).

Taken together, the 2024–2025 papers define the Quadratic Bubble Upwinding Method as a structured stabilization strategy in which artificial diffusion emerges from test-space design rather than from an ad hoc modification of the discrete operator (Bacuta et al., 2024). The method preserves a tridiagonal stencil in one dimension, admits exact matrix matching with finite-difference upwinding, supports quasi-optimal estimates in discrete optimal norms, avoids the implementation issues associated with exponential bubbles, and extends naturally to streamline-aligned multidimensional formulations (Bacuta et al., 2024, Bacuta, 4 Sep 2025).