VARIATIONAL STRUCTURE OF THE OPTIMAL ARTIFICIAL DIFFUSION
METHOD FOR THE ADVECTION-DIFFUSION EQUATION
K. B. NAKSHATRALA AND A. J. VALOCCHI

arXiv:0905.4771v3 [] 15 Dec 2010

Abstract. In this research note, we provide a variational basis for the optimal artificial diffusion method, which has been a cornerstone in developing many stabilized methods. The optimal
artificial diffusion method produces exact nodal solutions when applied to one-dimensional problems with constant coefficients and forcing function. We first present a variational principle for
a multi-dimensional advective-diffusive system, and then derive a new stable weak formulation.
When applied to one-dimensional problems with constant coefficients and forcing function, this
resulting weak formulation will be equivalent to the optimal artificial diffusion method. We present
representative numerical results to corroborate our theoretical findings.

1. INTRODUCTION
Many transport-related processes are modeled as advective-diffusive system. For example, transport of contaminants in subsurface flows is modeled as an advection-diffusion equation coupled with
Darcy flow. Except for very simple and limited problems, one cannot find analytical solutions and
hence must resort to numerical solutions. However, it is well-known that great care should be taken
in developing numerical formulations in order to avoid spurious oscillations due to the advective
term.
Many numerical formulations have been proposed, which fall under the realm of stabilized methods. One such method is the optimal artificial diffusion method, which has been the basis for
developing many stabilized methods and has also served as a benchmark for comparison [6, 2].
The optimal artificial diffusion method is derived by imposing the condition that it should produce
exact nodal solutions for one-dimensional problems with constant coefficients and forcing function.
This derivation, as one can see, does not have a variational basis.
Herein we outline a variational structure behind the optimal artificial diffusion method. We
start with a variational statement for the advection-diffusion equation (which is not a self-adjoint
operator), and derive a stable weak formulation. The resulting weak formulation when applied
to one-dimensional problems with constant coefficients and forcing function produces the same
difference equation as the optimal artificial diffusion method, which produces exact nodal solutions
Date: December 21, 2017.
Key words and phrases. variational principles, Euler-Lagrange equations, advection-diffusion equation, optimal
artificial diffusion.
1

when applied to such problems. This shows that the optimal artificial diffusion method has a firm
variational basis. This paper also highlights other possible routes in developing stable formulations
for non-self-adjoint operators.
An outline of this short paper is as follows. In Section 2, we present governing equations for an
advective-diffusive system, and also outline the standard weighted residual method. In Section 3, we
describe the optimal artificial diffusion method. In Section 4, we briefly discuss Vainberg’s theorem,
which provides a connection between the weighted residual statement and its corresponding scalar
functional (if it exists). In Section 5, we present a variational principle for the advection-diffusion
equation, and then derive a stable weak formulation. We then show that the resulting weak
formulation produces the same difference equation as the optimal artificial diffusion method for onedimensional problems with constant coefficients and forcing function. Finally, we draw conclusions
in Section 6.
Remark 1.1. It should be noted that there is a huge literature on developing stabilized (finite
element) formulations for an advective-diffusive system. For example, see [1, 3, 12, 4] and references
therein. A thorough discussion of these works is beyond the scope of this paper. In addition, these
discussions are not relevant to the subject matter of this paper as none of them discuss a variational
principle (that is, constructing a scalar functional) for an advection-diffusion system.
2. GOVERNING EQUATIONS: ADVECTIVE-DIFFUSIVE SYSTEM
Let Ω ⊂ Rd be a smooth and bounded domain, where “d” denotes the number of spatial dimensions, and ∂Ω denotes its smooth boundary. As usual, the boundary is divided into ΓD (the
part of the boundary on which Dirichlet boundary condition is prescribed) and ΓN (the part of
the boundary on which Neumann boundary condition is prescribed) such that ΓD ∪ ΓN = ∂Ω and
ΓD ∩ ΓN = ∅. Let x ∈ Ω denote the position vector, and the gradient and divergence operators
are denoted as “grad” and “div”, respectively. Let u : Ω → R denote the concentration, v(x) the
velocity vector field, and k(x) the symmetric and positive-definite diffusivity tensor. For further
discussion consider the following steady advective-diffusive system
(2.1)

v(x) · grad[u(x)] − div [k(x) grad[u(x)]] = f (x) in Ω

(2.2)

u(x) = up (x) on ΓD

(2.3)

n(x) · k(x)grad[u(x)] = tp (x)

on ΓN

where up (x) is the prescribed Dirichlet boundary condition, tp (x) is the prescribed Neumann
boundary condition, f (x) is the prescribed volumetric source, and n(x) denotes the unit outward
normal to the boundary. Note that the advective-diffusive operator is not self-adjoint. It is wellknown that, even under smooth functions for k(x), v(x) and f (x), the solution u(x) to equations
2

(2.1)–(2.3) may exhibit steep gradients close to the boundary especially for advection-dominated
problems [6] (and also see Figure 1, which will be described later in Section 3).
2.1. Notation and preliminaries. In the next section, we present a weak formulation for an
advective-diffusive system. To this end, let us define the following function spaces:

U := u(x) ∈ H 1 (Ω) | u(x) = up (x) on ΓD
(2.4a)

W := w(x) ∈ H 1 (Ω) | w(x) = 0 on ΓD
(2.4b)

where H 1 (Ω) is a standard Sobolev space defined on Ω. Note that the inner product for the above
vector spaces is the standard L2 inner product. That is,
Z
(2.5)
a(x) · b(x) dΩ
(a; b) :=
Ω

Similarly, one can define the weighted
(2.6)

L2

inner product
Z
µ(x)a(x) · b(x) dΩ
(a; b)µ :=
Ω

where µ : Ω →

R+

(where

R+

denotes the set of positive real numbers) is the (scalar) weight

function or measure density. Note that a weight function (which is used to define the weighted
inner product) should not be confused with weighting functions (which are sometimes referred to
as test functions). Using this weighted inner product one can define weighted Sobolev spaces as
n
o
(2.7a)
L2µ (Ω) := u(x) | (u; u)µ < +∞
n
o
(2.7b)
Hµ1 (Ω) := u(x) ∈ L2µ (Ω) | (grad[u]; grad[u])µ < +∞

Corresponding to the function spaces given in (2.4), we can define the following weighted function
spaces

Uµ := u(x) ∈ Hµ1 (Ω) | u(x) = up (x) on ΓD

Wµ := w(x) ∈ Hµ1 (Ω) | w(x) = 0 on ΓD

(2.8a)
(2.8b)

2.2. Standard weighted residual method and the Galerkin formulation. Let w(x) denote
the weighting function corresponding to u(x). A weak formulation based on the standard weighted
residual method for the advective-diffusive system given by equations (2.1)-(2.3) can be written as
Find u(x) ∈ U such that F(w; u) = 0 ∀w(x) ∈ W

(2.9)

where the bilinear functional F is defined as
Z
Z
grad [w(x)] · k(x)grad [u(x)] dΩ
w(x)v(x) · grad [u(x)] dΩ +
F(w; u) :=
Ω
Ω
Z
Z
(2.10)
w(x)tp (x) dΓ
w(x)f (x) dΩ −
−
Ω

ΓN

3

A finite element formulation corresponding to the standard weighted residual formulation can be
written as
(2.11)

Find uh (x) ∈ U h such that F(wh ; uh ) = 0 ∀wh (x) ∈ W h

where U h and W h are finite dimensional subspaces of U and W, respectively. One obtains the
Galerkin formulation by using the same function space for both U h and W h (except on ΓD ).
The Galerkin formulation produces node-to-node spurious oscillations for advection-dominated
problems [6]. The Galerkin method loses its best approximation property when the non-symmetric
convective operator dominates in the transport equation (for example, see the performance of the
Galerkin formulation in Figure 1). In other words, the Galerkin method is not optimal for solving
advection-dominated problems. In principle, it is possible to choose a small enough grid such that
the element Péclet number is less than one, and avoid spurious oscillations under the Galerkin
method. However, it may not always be practical to choose such a fine grid, and therefore, one
needs to employ a stabilized formulation to avoid unphysical oscillations and get meaningful results
on coarse grids. To understand this anomalous behavior several theoretical and numerical studies
have been performed, see [4, 6, 2] and references therein. One of those studies is a simple method
that gives nodally exact solutions for one-dimensional problems with constant coefficients, which is
commonly referred to as the optimal diffusion method.
3. A NODALLY EXACT FORMULATION AND OPTIMAL ARTIFICIAL DIFFUSION
Consider the following one-dimensional advection-diffusion equation with homogeneous Dirichlet
boundary conditions, and constant coefficients (that is, velocity, diffusivity and forcing function are
constants):
du
d2 u
−k 2 =f
dx
dx

∀x ∈ (0, 1)

(3.1)

v

(3.2)

u(x = 0) = 0 and u(x = 1) = 0

The analytical solution for the above problem is given by


1 − exp(vx/k)
f
(3.3)
x−
u(x) =
v
1 − exp(v/k)
which is plotted in Figure 1 for various value of v/k. As one can see, for large v/k we have
steep gradients near the outflow boundary. The above problem has been used as a benchmark for
developing many stabilized formulations [6].
For a numerical solution, let us divide the unit interval into N equal-sized elements (and hence
N +1 nodes), and define h := 1/N . Let the nodes be numbered as j = 0, · · · , N . Then, the position
4

vector of node j is xj = jh. The difference equation at an intermediate node (j = 1, · · · , N − 1)
arising from the Galerkin formulation for the above problem can be written as


v
2
Peh − 1
Peh + 1
uj−1 + h uj +
uj+1 = f
−
(3.4)
2h
Peh
Pe
Peh
where Peh = vh/(2k) is the element Péclet number, and uj is an approximate numerical solution at
node j. That is,
uj ≈ u(xj ) ∀j = 0, · · · , N

(3.5)

The difference equation (3.4) can be rearranged as
(3.6)

v

uj+1 − 2uj + uj−1
uj+1 − uj−1
−k
=f
2h
h2

The above equation basically reveals that the Galerkin formulation approximates the first- and
second-derivatives using a central difference approximation at intermediate nodes.
It is well-known that the Galerkin formulation is under-diffusive when applied to the advectiondiffusion equation, which is considered to be the reason why the formulation gives node-to-node
spurious oscillations. In Figure 1, we compare the numerical solution from the Galerkin formulation
with the analytical solution for the aforementioned one-dimensional problem. As one can see from
the figure, the Galerkin formulation produces spurious node-to-node oscillations for high Péclet
numbers.
To also understand this anomalous behavior of the Galerkin formulation, we now outline a nodally
exact formulation for the above problem, which is commonly referred to as the optimal artificial
diffusion formulation in the literature. To this end, we start with a difference equation at node j
of the form
(3.7)

β−1 uj−1 + β0 uj + β1 uj+1 = f

and the unknown coefficients (β−1 , β0 , and β1 ) are determined so that the above difference equation
gives nodally exact solution for the model problem given by equations (3.1)-(3.2) on the uniform
computational mesh described above. After simplification, the difference equation for the nodally
exact formulation can be written as: (see References [2, 6, 16]),

(3.8)

i
v h
−(1 + coth(Peh ))uj−1 + (2 coth(Peh ))uj + (1 − coth(Peh ))uj+1 = f
2h

where “coth” denotes the hyperbolic cotangent function. By rearranging the terms one can write
the above equation as follows
(3.9)

v

uj+1 − uj−1
uj+1 − 2uj + uj−1
− (k + k̄)
=f
2h
h2
5

where the artificial diffusion coefficient k̄ (which depends on the mesh size, and medium and flow
properties) is given by
vh
k̄ :=
2

(3.10)



coth(Peh )

1
− h
Pe



The difference equation (3.9) reveals that the nodally exact formulation (or the optimal artificial diffusion method) also employs a central difference approximation for the first- and second-derivatives
but solves a modified advective-diffusive system with an additional artificial diffusion given by the
coefficient k̄.
It should be noted that the above derivation for the optimal artificial diffusion method is not
based on a variational principle. In this paper we present a variational structure behind the
optimal artificial diffusion method, which, to the best of our knowledge, has not been reported in
the literature. We start by presenting a variational principle for an advective-diffusive system and
then derive a weak formulation. A finite element approximation of this weak form gives the optimal
artificial diffusion method.
We now present a variational principle for an advective-diffusive system. We then outline the
underlying variational structure behind the optimal artificial diffusion method.
4. VAINBERG’S THEOREM AND EXISTENCE OF A SCALAR FUNCTIONAL
Vainberg’s theorem [15] provides a connection between scalar functionals (also referred to as
“energy” functionals) and weighted residual statements. The theorem provides a criterion to establish when a scalar functional exists. The theorem also provides a formula to compute the scalar
functional (if it exists) from the weighted residual statement, which we will not invoke in this paper.
Vainberg’s theorem can be stated as follows. Let G(w; u) be a weighted residual functional,
which is linear with respect to w (but need not be linear with respect to u). There exists a scalar
functional E(u) such that



d
G(w; u) =
E(u + ǫw)
dǫ
ǫ=0
if and only if
(4.1)






d
d
G(w, u + ǫw̄)
=
G(w̄, u + ǫw)
dǫ
dǫ
ǫ=0
ǫ=0

For a simple proof of Vainberg’s theorem see Reference [9, Chapter 9].
Returning back to our discussion on advective-diffusive system, using Vainberg’s theorem one
can conclude that there exists no scalar functional E(u) whose directional derivative along w(x)
gives the (bilinear) functional in the standard weighted residual method (which is given by equation
6

(2.10)). That is, there is no scalar functional E(u) such that


d
F(w; u) =
E(u + ǫw)
dǫ
ǫ=0
However, we should note conclude that the advective-diffusive system does not have a variational
principle. All we concluded above (invoking Vainberg’s theorem) is that the standard weighted
residual statement (2.10) cannot be obtained as a directional derivative of a scalar functional. It
may be possible to construct a different weighted residual statement with different weight function
(or measure) that can be obtained as a directional derivative of a scalar functional. Stated differently, just because the operator is non-self-adjoint does not mean that there exists no variational
statement.
In a seminal paper, Tonti [14] has correctly highlighted the point that the symmetry of a bilinear
form in the sense of equation (4.1) (which according to Vainberg’s theorem is a necessary and
sufficient condition for the existence of a variational formulation) depends on the choice of the inner
product. (Note that Tonti has used the acronym “bilinear functional” to mean inner product.) One
of the ways to meet the symmetry requirement is by changing the underlying inner product. Tonti
[14] has also shown that changing the inner product is equivalent to transforming the given problem
into a different problem that has the same solution(s). Another related work is by Thangaraj and
Venkatarangan [13] who have presented dual variational principles with applications to magnetohydrodynamics. A related and important work is by Magri [11] who has developed a procedure to
select an appropriate inner product for any given linear operator to meet the symmetry requirement.
Following the discussion in References [14, 13], in the next section we present a variational principle
for an advection-diffusion system.
Remark 4.1. For self-adjoint operators (e.g., the Laplacian operator) one can show that the standard weighted residual method can be obtained as a directional derivative of a scalar functional.
5. A VARIATIONAL PRINCIPLE FOR THE ADVECTION-DIFFUSION EQUATION
Consider the following constrained extremization problem
Z
(5.1)
F (x, u(x), grad[u(x)]) dΩ
extremize
u(x)

Ω

where F : Ω × R × Rd → R is a known scalar functional. It is well-known (for example, see [5, 10])
that the Euler-Lagrange equation for the above extremization problem is given by
(5.2)

Fu (x, u(x), grad[u(x)]) − div [Fp (x, u(x), grad[u(x)])] = 0 ∀x ∈ Ω

where Fu denotes the partial derivative with respect to u(x), and Fp denotes the vector of partial
derivatives with respect to the components of grad[u(x)].
7

We now construct a scalar functional for the advection-diffusion equation. To this end, let us
define the scalar function α(x) as
 Z
α(x) := exp −

(5.3)

k

−1



(x)v(x) · dx

Note that the integral inside the exponential operator in equation (5.3) is an indefinite line integral.
Since the tensor k(x) is positive-definite, its inverse exists, and hence the scalar function α(x) is
well-defined. It is important to note that α(x) > 0 ∀x ∈ Ω (as we have assumed the domain to be
bounded and smooth). Now define the functional for the advection-diffusion equation as


1 p
F (x, u(x), grad[u(x)]) := α(x)
(5.4)
∀x ∈ Ω
k k(x)grad [u(x)] k2 − u(x)f (x)
2
p
where k · k denotes the standard Euclidean norm, and k(x) denotes the square root of the
symmetric and positive-definite tensor k(x). (Note that the square root of a symmetric positive-

definite tensor is well-defined, and is also symmetric and positive-definite. See Halmos [8] and also
Gurtin [7, page 13].) A straightforward calculation shows that the Euler-Lagrange equation of the
functional (5.4) is
(5.5)

v(x) · grad[u(x)] − div [k(x) grad[u(x)]] = f (x)

which is the (multi-dimensional) advection-diffusion equation. For the advective-diffusive system
given by equations (2.1)-(2.3) (that is, including the boundary conditions), the variational statement
can be written as
I(u)

(5.6)

extremize

(5.7)

subject to u(x) = up (x) on ΓD

u(x)

where the functional I(u) is defined as
Z
Z
(5.8)
F (x, u(x), grad[u(x)]) dΩ −
I(u) :=
Ω

α(x)u(x)tp (x) dΓ
ΓN

and the scalar function α(x) is same as before, equation (5.3).
We now use this variational principle to derive a new weak formulation for the advective-diffusive
system. We later show that for constant coefficients and homogeneous Dirichlet boundary conditions, this formulation is same as the optimal artificial diffusion method.
Remark 5.1. It is well-known that, using the least-squares method one can construct a minimization problem for a given partial differential equation of the form Lu − f = 0 as follows:
Z
(5.9)
kLu − f k2 dΩ
minimize
u(x)

Ω

8

But it should be noted that the Euler-Lagrange equation of the above minimization problem need
not be the original differential equation Lu − f = 0. For example, if we consider the differential
equation to be d2 u/dx2 = 0 in Ω = (0, 1), the corresponding minimization problem based on the
least-squares method takes the following form:
(5.10)

minimize

Z

0

1  2 2
d u

dx2

dx

However, the Euler-Lagrange equation of the above minimization problem is d4 u/dx4 = 0 (which
is not same as the differential equation we considered). Also, another difference between the minimization problem based on the least-squares method and the variational problem is in the regularity
requirements. For example, in the minimization problem based on the least-squares method (5.10),
u should be twice differentiable (or, more precisely, u ∈ H 2 ((0, 1))). On the other hand, the (standard) variational principle for the differential equation d2 u/dx2 = 0 requires that u be differentiable
once (that is, u ∈ H 1 ((0, 1))).
5.1. A new stable formulation. A necessary condition that the extremum of the functional I(u)
has to satisfy is
(5.11)




d
I(u(x) + ǫw(x))
= 0 ∀w(x)
dǫ
ǫ=0

For convenience let us define the bilinear functional B(w; u) as


Z
d
α(x)grad [w(x)] · (k(x)grad [u(x)]) dΩ
I(u(x) + ǫw(x))
=
B(w; u) :=
dǫ
Ω
ǫ=0
Z
Z
(5.12)
α(x)w(x)tp (x) dΓ
− α(x)w(x)f (x) dΩ −
ΓN

Ω

A new weak formulation for the advective-diffusive system can be written as
(5.13)

Find u(x) ∈ Uα such that B(w; u) = 0 ∀w(x) ∈ Wα

where Uα and Wα are weighted function spaces defined in equation (2.8) with weight function
µ(x) = α(x). Recall that the scalar function α(x) is defined in equation (5.3). A corresponding
finite element formulation can be written as
(5.14)

Find uh (x) ∈ Uαh such that B(wh ; uh ) = 0 ∀wh (x) ∈ Wαh

where Uαh and Wαh are finite dimensional subspaces of Uα and Wα , respectively. It is important to
note the formulation (5.13) (and hence the formulation (5.14)) is valid even for spatial dimensions
d = 1, 2, 3.
9

Remark 5.2. Some notable differences between the bilinear functions F(w; u) and B(w; u) are as
R
follows. A non-symmetric term similar to Ω w(x)v(x) · grad[u(x)] dΩ is not present in the bilinear

functional B(w; u). The bilinear functional B(w; u) is symmetric in the sense that





d
d
B(w; u + ǫw̄)
=
B(w̄; u + ǫw)
dǫ
dǫ
ǫ=0
ǫ=0

which is not the case for the bilinear functional F(w; u). The symmetry of B(w; u) (which, according
to Vainberg’s theorem, is necessary and sufficient for the existence of a scalar functional) should
not be surprising as we defined the bilinear functional B(w; u) (5.12) as a directional derivative of
a scalar functional. Another notable difference is that the bilinear functional F(w; u) is based on
the standard L2 inner product, which has the weight function (or measure density) to be 1. The
bilinear functional B(w; u) is based on the weighted L2 inner product with the weight function equal
to α(x) > 0 (which is defined in equation (5.3)).
5.2. Application to the one-dimensional problem. We now apply the new formulation to the
one-dimensional problem defined in Section 3, and compare the difference equation produced by
this new stable formulation with the one produced by the optimal artificial diffusion method given
by equation (3.9) for an intermediate node. For the case at hand, the coefficients v and k are
constants, and hence α(x) = exp[−vx/k]. The variable u(x) is interpolated as
(5.15)

u(x) = uj−1 N−1 (x) + uj N0 (x) + uj+1 N+1 (x) xj−1 ≤ x ≤ xj+1

where the shape functions are defined as

(5.16a)

(5.16b)

(5.16c)

(x−xj )
h

N−1 (x) =

(

−

N0 (x) =

(

(x−xj )+h
h
h−(x−xj )
h

N+1 (x) =

(

0

−h ≤ x − xj ≤ 0

(x−xj )
h

0 < x − xj ≤ +h

−h ≤ x − xj ≤ 0
0 < x − xj ≤ +h

0

−h ≤ x − xj ≤ 0
0 < x − xj ≤ +h

The weighting function is also similarly interpolated. The difference difference equation produced
by the new stable formulation at an intermediate node can be written as
(5.17)

γ−1 uj−1 + γ0 uj + γ+1 uj+1 = f
10

The coefficients γ−1 , γ0 and γ+1 can be written as
R +h
′
′
−h α(x)k(x)N0 (x)N−1 (x) dx
γ−1 =
(5.18a)
R +h
−h α(x)N0 (x) dx
 ′ 2
R +h
α(x)k(x)
N0 (x) dx
−h
(5.18b)
γ0 =
R +h
−h α(x)N0 (x) dx
R +h
′
′
−h α(x)k(x)N0 (x)N+1 (x) dx
(5.18c)
γ+1 =
R +h
−h α(x)N0 (x) dx

where a superscript prime denotes a derivative with respect to x. After simplification, the coeffi-

cients can be written as
(5.19)

γ−1 = −



v 
v
v 
1 + coth(Peh ) , γ0 = coth(Peh ), γ+1 =
1 − coth(Peh
2h
h
2h

As one can see, the coefficients γ−1 , γ0 and γ+1 are, respectively, the same as the coefficients
β−1 , β0 and β+1 (see equations (3.7) and (3.8)), which are obtained using the optimal artificial
diffusion method. In Figure 2, we compare the numerical solutions obtained using the new stable
formulation with the analytical solutions. As predicted by the theory, the new stable method
produces nodally exact solution for all Péclet numbers. Both the theoretical and numerical studies
have shown that the stable weak formulation proposed in the previous section is equivalent to the
optimal artificial diffusion method. Hence, we have provided a variational basis for the artificial
diffusion method, which has academic importance, and also has been the basis in developing and
testing many stabilized formulations.
Although the new stable formulation is valid even in two- and three-dimensions, preliminary
numerical studies have indicated that this formulation may not be practically feasible for largescale problems due to the following reasons: evaluating the exponential function is expensive, one
may require more Gauss points per element in order to evaluate integrals in the weak form, and
also the resulting stiffness matrix will be ill-conditioned without an appropriate preconditioner.
However, the new stable formulation, as shown in this paper, does have theoretical significance. It
also highlights a possible route of developing robust stable formulations. This paper also highlights
the relevance of the theoretical studies by Tonti on variational principles (for example, Reference
[14]) to develop new stabilized formulations.
6. CONCLUSIONS
It is well-known that the classical Galerkin formulation (which is based on the standard weighted
residual method) produces spurious node-to-node oscillations for the advection-diffusion equation
for advection-dominated problems. For self-adjoint operators (e.g., the Laplacian equation), the
11

standard weighted residual formulation can be obtained as a directional derivative of a scalar
functional. However, for the advection-diffusion equation (using Vainberg’s theorem) it can be
shown that the standard weighted residual method cannot be obtained as a directional derivative
of a scalar functional.
In this paper, we presented a variational principle for an advective-diffusive system, and derived a stable weak formulation. This resulting weak formulation when applied to one-dimensional
problems gives rise to the same difference equation as the optimal artificial diffusion method, which
produces exact nodal solutions when applied to one-dimensional problems with constant coefficients
and forcing function. Hence, we have provided a variational basis for the optimal artificial diffusion
method, which has been the cornerstone in developing many stabilized methods. We presented
representative numerical results to corroborate our theoretical findings
ACKNOWLEDGMENTS
The first author (K. B. Nakshatrala) acknowledges the financial support given by the Texas
Engineering Experiment Station (TEES). The second author (A. J. Valocchi) was supported by
the Department of Energy through a SciDAC-2 project (Grant No. DOE DE-FCO207ER64323).
The opinions expressed in this paper are those of the authors and do not necessarily reflect that of
the sponsors.
References
[1] R. Codina. Comparison of some finite element methods for solving the diffusion-convection-reaction equation.
Computer Methods in Applied Mechanics and Engineering, 156:185–210, 1998.
[2] J. Donea and A. Huerta. Finite Element Methods for Flow Problems. John Wiley & Sons, Inc., Chichester, UK,
2003.
[3] L. P. Franca, S. L. Frey, and T. J. R. Hughes. Stabilized finite element methods: I. Application to the advectivediffusive model. Computer Methods in Applied Mechanics and Engineering, 95:253–276, 1992.
[4] L. P. Franca, G. Hauke, and A. Masud. Revisiting stabilized finite element methods for the advective-diffusive
equation. Computer Methods in Applied Mechanics and Engineering, 195:1560–1572, 2006.
[5] I. M. Gelfand and S. V. Fomin. Calculus of Variations. Dover Publications, New York, USA, 2000.
[6] P. M. Gresho and R. L. Sani. Incompressible Flow and the Finite Element Method: Advection-Diffusion, volume 1.
John Wiley & Sons, Inc., Chichester, UK, 2000.
[7] M. E. Gurtin. An Introduction to Continuum Mechanics. Academic Press, San Diego, USA, 1981.
[8] P. R. Halmos. Finite-Dimensional Vector Spaces. Springer-Verlag, New York, USA, 1993.
[9] K. D. Hjelmstad. Fundamentals of Structural Mechanics. Springer Science+Business Media, Inc., New York,
USA, second edition, 2005.
[10] J. Jost and X. Li-Jost. Calculus of Variations. Cambridge University Press, Cambridge, UK, 1998.
[11] F. Magri. Variational formulation for every linear equation. International Journal of Engineering Science, 12:537–
549, 1974.
12

[12] A. Masud and R. A. Khurram. A multiscale/stabilized finite element method for the advection-diffusion equation.
Computer Methods in Applied Mechanics and Engineering, 193:1997–2018, 2004.
[13] D. Thangaraj and S. N. Venkatarangan. Dual variational principles for a class of non-linear partial differential
equations. IMA Journal of Applied Mathematics, 30:21–26, 1983.
[14] E. Tonti. Variational formulation for every nonlinear problem. International Journal of Engineering Science,
22:1343–1371, 1984.
[15] M. M. Vainberg. Variational Methods for the Study of Nonlinear Operators. Holden-Day, Inc., San Francisco,
USA, 1964.
[16] O. C. Zienkiewicz and R. L. Taylor. The Finite Element Method: Fluid Dynamics. John Wiley & Sons, Inc.,
New York, USA, fifth edition, 2000.

13

1.6
1.4
1.2

u(x)

1

h
e

P = 5.0
0.8
0.6
0.4

h
e

h
e

P = 0.5

P = 0.1
0.2

h

0
0

0.1

0.2

0.3

0.4

P = 0.05
e
0.5
0.6
0.7
x

0.8

0.9

1

Figure 1. In this figure, we compare the Galerkin formulation (which is denoted
using solid squares and dotted lines) with the analytical solution (which is denoted
using solid continuous lines) for one-dimensional advection-diffusion equation with
homogeneous Dirichlet boundary conditions for various v/k ratios. We have taken
the forcing function to be unity (i.e., f = 1), and h = 0.1. As expected, the Galerkin
formulation produced spurious node-to-node oscillations for high Péclet numbers.
Correspondence to: Kalyana Babu Nakshatrala, Department of Mechanical Engineering, 216 Engineering/Physics Building, Texas A&M University, College Station, Texas-77843, USA. TEL: +1-979845-1292
E-mail address: knakshatrala@tamu.edu
Albert Valocchi, Department of Civil and Environmental Engineering, University of Illinois at
Urbana-Champaign, Urbana, Illinois-61801, USA.
E-mail address: valocchi@uiuc.edu

14

1
h
e

P = 5.0
0.8

h
e

0.6
u(x)

P = 0.5

0.4
h
e

P = 0.1
0.2

0
0

h
P
e

0.1

0.2

0.3

0.4

0.5
x

= 0.05

0.6

0.7

0.8

0.9

1

Figure 2. In this figure, we compare the numerical solution using the new stable formulation (which is denoted using solid squares) with the analytical solution (which is denoted using solid continuous lines) for one-dimensional advectiondiffusion equation with homogeneous Dirichlet boundary conditions for various v/k
ratios. We have taken the forcing function to be unity (i.e., f = 1), and h = 0.1. As
predicted by the theory, the new stable formulation produces nodally exact solutions
for the chosen one-dimensional problem for all Péclet numbers. The continuous line
denotes the analytical solution, and solid squares denote the numerical solution from
the new stable formulation.

15

