






























Finite Element Method; 
Varational Formulation; 
Numerical integration

An advection-diffusion equation (ADE) is
a mathematical model that has been used
to model the concentration of pollutants.
It gives the amount of pollutant concen-
tration fields a fter i nput o f t he velocity
data from the hydrodynamic model which
are derived from mass balances. Formally
the ADE equation is given by

ut + a∇u = ∇(D∇u) + f, (1)

where u is the concentration of the pol-
lutant, a is the velocity of the considered

particle, D the diffusion coefficient and f
defines the sources and sinks due to dif-
ferent processes.

East African Journal of Biophysical and Computational Sciences 

Journal homepage : https://journals.hu.edu.et/hu-journals/index.php/eajbcs 
College of Natural & Com

H
pu
aw

tat
as

ion
sa U

al 
n
Sc
ive

ie
r
nc
sit

e
y
s

Year 2021

Volume xx No xx

KEYWORDS: ABSTRACT 

* Corresponding author:
Email: kassahunm@hu.edu.et  +251918491415  https://dx.doi.org/10.4314/eajbcs.v4i1.5S

East Afr. J. Biophys. Comput. Sci. (2023), Vol. 4, Issue. 1, 52-74

For the vast majority of geome-
tries and problems, Eq. 1 cannot be
solved with analytical methods, and an
approximation of the equations can be
constructed with different types of dis-
cretizations. Many numerical schemes
have been implemented to approximately
solve the ADE (Lima et al., 2021; Mah-
mud, 2012; Pochai and Deepana, 2011;
Lian et al., 2016; Szymkiewicz and Gka-

Kassahun Getnet Mekonen*, Zerihun Kinfe Birhanu

Department of Mathematics, College of Natural and Computational Sciences, Hawassa University, 
Hawassa, Ethiopia. 

In this paper, we have implemented the finite element
method for the numerical solution of a boundary and initial
value problems, mainly on solving the one and two
dimensional advection-diffusion equation with constant
parameters. In doing so, the basic idea is to first rewrite the
problem as a variational equation, and then seek a solution
approximation from the space of continuous piece-wise
linear’s. This discretization procedure results in a linear
system that can be solved by using a numerical algorithm
for systems of these equations. The techniques are based on
the finite element approximations using Galerkin’s method in
space resulting system of the first order ODE’s and then
solving this first order ODE’s using backward Euler
descritization in time. For the two dimensional problems,
we use the ODE solver ODE15I to descritize time. The
validity of the numerical model is verified using different
test examples. The computed results showed that the use of
the current method is very applicable for the solution of the
advection-diffusion equation.

Numerical Solutions of Advection Diffusion Equations Using Finite Element Method 

INTRODUCTION

 
Research article



equation by a test function and integra-
tion by parts (Green-Gauss Theorem) to
reduce second order derivatives to first or-
der terms, i.e., weak formulation. Then
we represent the approximate solution by
the linear combination of basis functions,
by constructing a set of basis functions
based on the mesh of our domain. That
is, the solution u can be approximated by
a function uh using the linear combina-
tions of the basis functions ϕi according
to the following expressions:

u ≈ uh =
∑

uiϕi. (2)

and we solve the linear system of equa-
tions to obtain the approximate solution.

To derive the advection diffusion equa-
tions for the application of pollution mod-
els, consider an elementary water body. 
Water quality within this body depends 
on the polluting substance mass present 
there. The water quality models describe 
the change in the mass of a polluting 
substance within the water body. The 
change is calculated as the difference be-
tween mass-flows (mass fluxes) entering 
and leaving this water body, considering 
also the effects of internal sources and 
sinks of the substance, if any. The mecha-
nism of mass transfer into and out of this 
water body includes the following pro-
cesses:

� Mass is transported by the flow, a,
of the velocity vector. This pro-
cess is termed as the advective mass
transfer. The transfer of mass, that
is the mass flux can be calculated as
u × a, where u is the concentration
of the substance in the water.

� The dispersive mass transfer is usu-
ally expressed by the law of Fix
which states that the transport of

East Afr. J. Biophys. Comput. Sci. (2023), Vol. 4, No. 1, 52-74

52

Here,  ϕi denotes the basis functions
and ui denotes the coefficients of the func-
tions that approximate u with uh. After
this we get a system of linear equations

siorowski, 2021). In numerical method, a
discrete approximations for the solution
is computed by descritizing the given do-
main in to different sets of sub domains.
In this paper we were focus on finite ele-
ment method to solve the PDE given in
Eq. 1. We will implement the method for
one and two dimensional PDE’s. In this
method, first we develop a weak formula-
tion, from which we derive the discretiza-
tion by multiplying both sides of the ADE

The development of finite element
method has favored by the progress of
computer technology and numerical cal-
culus, and originally applied for mechan-
ical structures (Lima et al., 2021; Mah-
mud, 2012; Donea and Huerta, 2003) .
Several procedures have been tried to in-
terpret separately the advection and dif-
fusion pollutant transport. The FEM can
help to face more complex problems and
the privilege importance of the method is
that it can be adapted to complex geom-
etry domains, but the element wise inter-
vals can assume any form of size, and ob-
viously there is an expense of more bur-
densome calculations. We had to do the
descritization process of the finite element
method for the 1D and 2D Poisson equa-
tions in which it is an auxiliary step in
solving the advection diffusion equation
with the FEM. Usually the numerical so-
lutions of PDE’ including the equation
(Eq. 1) are done with the finite difference
method.

MATHEMATICALMODEL 
FORMULATION



the substance in the direction of a
space is proportional to the concen-
tration of this substance in that di-
rection and the proportionality fac-
tor being the coefficient of disper-
sion, D∇u.

transport by advection = audA ,

transport by diffusion = Dx
∂u
∂x

,

where dA is an elemental cross-sectional
area of the cubic element, and Dx is the
diffusion coefficient in the x-direction.

∂

∂t
u(x, t) +

∂

∂x
(a(x, t)u(x, t)) =

∂

∂x

(
d(x, t)

∂

∂x
u(x, t)

)
+ f(t, x). (3)

space interval Ω ⊂ R with time t ≥ 0.
An initial condition u(x, 0) will be given
and we also assume that suitable bound-
ary conditions are provided, and for our
work we consider the velocity field and the
diffusion term as constants.

pressions of the advection equation a∂u
∂x

which has a first-order derivative, and
the diffusion equation D ∂2u

∂x2 that has a
second-order derivative.

East Afr. J. Biophys. Comput. Sci. (2023), Vol. 4, No. 1, 52-74

52

Both advection and diffusion move a
pollutant material from one place to an-
other, but each accomplishes this differ-
ently. The essential difference of the ad-
vection and diffusion is that advection
moves the pollutants in one way (down-
stream) but diffusion goes in both ways
(regardless of a stream direction). This
is seen in the respective mathematical ex-

Questions are arise like, can we have
cases of fast advection and relatively weak
diffusion and other cases of negligible ad-
vection and fast diffusion? To answer this

Assuming that the two components
(advection and diffusion) may be super-
posed, the total amount of material trans-
ported parallel to any given direction is
obtained by summing the advective and
diffusive transports. Using the mass bal-
ance approach by equating the difference
between the mass of material entering a
volume element and that leaving the ele-
ment (i.e., net influx of mass) to the rate
of accumulation of mass inside the vol-
ume.

By considering a volume element of
porous mediums in three dimensional
cartesian coordinates the equations are
derived (Bajellan, 2015). Since we are
considering advection and diffusion as
the two modes of transport of a fluid
within the porous medium, we can rep-
resent these two transport modes in the
x-direction mathematically as:

The multidimensional advection-
diffusion equation is used for analyzing
mixing problems in rivers. One of the
practical difficulties is that the equation
requires some prior information about wa-
ter depths, velocities, and diffusion co-
efficients, which could not conveniently
be gathered in field experiments. In
some particular mixing problems, how-
ever, some of the terms in the multi-
dimensional advection-diffusion equation
are negligibly small, so that the problem
can be simplified by reducing the model to
one dimension (Lima et al., 2021; Hunds-
dorfer, 1996) and the one dimensional
advection diffusion equation is given in
Eq. 3.

The time-derivative term expresses
accumulation of mass at a point in space,
the advection term a∇u transport of mass
with the flow, and the diffusion term
D∇2u reflects transport of mass due to
molecular diffusion (Langtangen, 1999).
We shall consider the equation in the



question, we must compare the sizes of
the a∂u

∂x
and D ∂2u

∂x2 terms to each other,
and this is accomplishes by introducing
scales.

Using these scales, we can derive es-
timates of the sizes of the different terms.
Since the derivative ∂u

∂x
is expressing the

difference in concentration over a distance
of infinitesimal limit, we can estimate it
to be approximately U

X
, and the advection

term scales as:

a
∂u

∂x
∼ V

U

X
.

Similarly, the second derivative ∂2u
∂x2 rep-

resents the difference of a gradient over
a specified distance and is estimated at
( U
X
)

X
= U

X2 , and the diffusion term scales
as:

D
∂2u

∂x2
∼ D

U

X2
.

Equipped with these estimates, we can
then compare the two processes by form-
ing the ratio of their scales:

Advecton

Diffusion
=

V U
X

D U
X2

=
V X

D

This ratio is dimensionless and Tradition-
ally, it is called the Peclet number and is
denoted by Pe:

Pe =
V X

D
.

If Pe ≪ 1 (if Pe < 0.1): the advection
term will result significantly smaller than
the diffusion term. Physically, diffusion
dominates and advection is negligible. So,
spreading occurs symmetrically despite of
the flow of the directional bias. If we wish
to simplify the problem, we may drop the
a∂u
∂x

term, as if a were nil (no amount at
all). The relative error occurred in the so-
lution is expected to be on the order of the

Peclet number, and the smaller Pe leads
to the smaller error. The solutions estab-
lished with diffusion only were based on
such simplification and are thus valid as
long as Pe ≪ 1.

If Pe ≫ 1 (if Pe > 10): the advec-
tion term is significantly bigger than the
diffusion term. Physically, the diffusion
term is negligible and advection domi-
nates, and spreading is existent, with the
patch (small area) of pollutant being sim-
ply moved along by the flow. If we wish
to simplify the problem, we may drop the
D ∂2u

∂x2 term, as if D were zero. The rela-
tive error occurred in doing the solution is 
expected as the order of the Peclet num-
ber inverse (1/P e), and the larger P e will 
result the smaller error.

East Afr. J. Biophys. Comput. Sci. (2023), Vol. 4, No. 1, 52-74

55

If P e ∼ 1 (in practice, if 0.1 < P e < 10):
 the advection and diffusion terms are not 
significantly different which results for the 
non  dominance  of  the  two  in  the  pro-
cess. The full equation must be utilized as 
there will no approximation to the equa-
tion will be justified.

Variable scale choice of value
u U The concentration value such as initial, boundary, or average value
a V The maximum velocity value
x X Approximate length of the domain or size of release location

In this paper we use the finite element 
method (FEM) to approximate the solu-
tion of the advection diffusion equation. 
The method is examined as an emerging 
tool for the approximate solution of differ-
ential equations describing different phys-
ical processes (Yang et al., 2020). It is 
based on the basic finite e lement proce-
dures, those are: the variation form for-

NUMERICAL METHOD



mulation of the problem, the discretiza-
tion of the formulation in a finite element,
and the solution of the resulting finite el-
ement equations. FEM cuts a given do-
main into several elements (pieces of the
domain) and connected in a finite number
of nodal points.

Finite Element Implementa-
tion of the 1D Governing
Equation

Let us consider the one dimensional
Advection-Diffusion equation given by:

∂u

∂t
+ a

∂u

∂x
= D

∂2u

∂x2
+ f, u(0) = u(1) = 0,

(4)
where u is the concentration of the pol-
lutant, a is the velocity, f is the source
term, and D is the diffusion coefficient,
with all the three variables be constants.

Implementation of the 1D Advection
Equation

Let we first consider only an advection
equation, that is the diffusion term does
not exist;

∂u

∂t
+ a

∂u

∂x
= 0.

The first step is constructing a variational
or weak formulation, by multiplying both
sides of the differential equation by a
test function v(x) satisfying the bound-
ary conditions (BC) v(0) = 0, v(1) = 0
and v ∈ H1

0 (0, 1), where H1
0 (0, 1) is the

Sobolev space,

H1
0 (0, 1) = {v ∈ L2(0, 1); v′ ∈ L2(0, 1)},

and it is a function space where all the
functions are bounded. Let us now define
a sub-space of H where we can find our
solution u. We call this V and

V = {v ∈ H(X) : v|∂Ω = 0}, where Ω is our domain,

Then multiplying and integrating both
sides in the domain we have that:

∂u

∂t
.v + a

∂u

∂x
.v = 0,

∫ 1

0

(
∂u

∂t
.v + a

∂u

∂x
.v

)
= 0.

East Afr. J. Biophys. Comput. Sci. (2023), Vol. 4, No. 1, 52-74

56

The FEM is based on the integration of
the terms in the equation to be solved, in
form of point discretization schemes.
It utilizes the method of weighted residu-
als and integration by parts (Green-Gauss
Theorem) to reduce second order deriva-
tives to first order terms. The solu-
tion domain is discretized into individ-
ual elements and these elements are op-
erated upon individually and then solved
globally using matrix solution techniques.
Such a task could be done automatically
by a computer, but it necessitates an
amount of mathematical skill that to day
still requires human involvement, (Bren-
ner et al., 2008).

The theories of finite element meth-
ods provided the reasons why it worked
well for the class boundary/initial value

problems (Lima et al., 2021; Ahsan,
2012; Larson and Bengzon, 2010). Ex-
tension of the mathematical basis to non-
linear and non-structural problems was
achieved through the method of weighted
residuals (MWR), originally conceived by
Galerkin in the early 20th century. The
basics of the method requires multiplying
of the governing differential equation by a
set of predetermined weights and integrat-
ing the resulting product over a region.
Most of the finite element method uses
the Galerkin’s method to establish the ap-
proximations of the governing equations,
(Aragonés et al., 2019; Lima et al., 2021;

Ahsan, 2012; Yang et al., 2020; Bren-

ner et al., 2008). It allows us to convert a
continuous form of the problem, such as
the weak formulation for the partial dif-
ferential equation into a discrete problem
that may be solved numerically.



That is;∫ 1

0

∂u

∂t
.v +

∫ 1

0

a
∂u

∂x
.v = 0, (5)

which is the weak formulation of the one
dimensional advection equation.

Advantages of weak form compared
to strong form

Equation 5 is the final weak formulation.
It is equivalent to the strong form, since
we can reverse all the steps, and get back
to the original equation. Firstly, if we
look at the strong form, we have two
separate partial derivatives of u, so the
strong form requires that u be continu-
ously differentiable until at least second
partial derivative. Our new formulations
has lowered this requirement to only first
partial derivatives by transforming one of
the partial derivatives onto the weight-
function v(x, y). This is the first big ad-
vantage of a weak formulation. The sub-
space V is not difficult to understand; it
is a subspace of H because our weak form
requires that the functions are in H; our

strong form requires that u be 0 along the
boundary, so V is the subspace of all func-
tion which are zero on the boundary.

ih, i = 0, 1, ..., n, where h = 1
n
, and

ϕi(xj) = δij, i, j = 1, ..., n− 1, (6)

where δij being the Kronecker delta. The
function ϕi is therefore piece wise linear
and are fix with one node (vertex) and

Figure 1: The basis (hat) function ϕi associated to node xj, in this figure φi on a mesh.
Also shown is the half hat ϕ0.

Its expression is given by

ϕi(x) =


x−xi−1

xi−xi−1
, if xi−1 ≤ x ≤ xi

xi+1−x
xi+1−xi

, if xi ≤ x ≤ xi+1

0, other wise

, for i = 1, 2, . . . , n− 1. (7)

That is,

ϕi(xj) = δij =

{
1, if i = j

0, if i ̸= j
. (8)

East Afr. J. Biophys. Comput. Sci. (2023), Vol. 4, No. 1, 52-74

57

The next step is to generate a mesh,
let be a uniform Cartesian mesh xi =

associate the value one to this node and
zero at the remaining nodes of the par-
tition (see Fig. 1, (Larson and Bengzon,
2010)).

we define the intervals as [xi−1, xi], i =
1, 2, ..., n.
After generating a mesh we construct a
set of basis functions based on the mesh
for each intervals, such as the piece wise
linear functions for i = 1, 2, ..., n− 1. The
characteristic basis functions are charac-
terized by the following property, (Quar-
teroni and Quarteroni, 2009)



uh(x) =
n−1∑
j=1

cjϕi(x), (9)

where the coefficients cj are the unknowns
to be determined. since the hat (basis)
functions are piece wise linear, uh(x) is
also a piece wise linear function, although
this is not usually the case for the true

solution u(x), and here we have,

uh(xj) =
n−1∑
i=1

cjϕi(xj) = cj.

We then derive a linear system of equa-
tions for the coefficients by substituting
the approximate solution uh(x) for the ex-
act solution u(x) in the weak form:

u(t, x) =
N−1∑
i=1

uiϕi(x), (10)

v(t, x) =
N−1∑
j=1

vjϕj(x). (11)

Now substituting Eq.10 and Eq.11
in the weak formulation of the equation
Eq.5, we have:

∫ 1

0

∂

∂t

N−1∑
i=1

uiϕi.

N−1∑
j=1

vjϕj +

∫ 1

0

a
∂

∂x

(
N−1∑
i=1

uiϕi

)
.

N−1∑
j=1

vjϕj = 0.

which then implies,

N−1∑
j=1

vj

(
∂

∂t

N−1∑
i=1

ui

∫ 1

0

ϕi.ϕj + a
N−1∑
i=1

ui

∫ 1

0

ϕ′
i.ϕj

)
= 0.

That is
∂

∂t

N−1∑
i=1

ui

∫ 1

0

ϕi.ϕj + a

N−1∑
i=1

ui

∫ 1

0

ϕ′
i.ϕj = 0.

In a matrix form it can be written as:

MU̇ + aBU = 0, (12)

where,

East Afr. J. Biophys. Comput. Sci. (2023), Vol. 4, No. 1, 52-74

58

We use this hat (basis) functions
through out the 1D space of this research
using equal spaced step size (xi+1−xi = h,
for all i). We say that the functions are
basis for the following reasons. If we want
to approximate our continuous function u
with a piece wise continuous linear func-

Then let the approximate solution for u
be given by a linear combination of ba-
sis functions ϕi = δij, as given in Eq. 10,
and also for v as in Eq. 11. Now we find
a finite element solution of the discrete
problem by using the hat functions ϕi(x)
defined in Eq. 7. For the given basis func-
tion the approximation of u and v can be
written as :

tion u′, these functions are what we need.
These functions are linearly independent
of each other; it is not possible to make
one out of a combination of others. For
example, only one of these functions, ϕi,
is non-zero (equal to 1) at node i. The
next step in approximating a PDE with
FEM is represent the approximate (FE)
solution by the linear combination of such
basis functions, (Quarteroni and Quar-
teroni, 2009) as



� M is the mass matrix with entries:

Mi,j =

∫ 1

0

ϕi(x)ϕj(x)dx.

� B is a matrix with entries:

Bij =

∫ 1

0

ϕ′
i(x)ϕj(x)dx =


0, if i = j
−1
2
, if i− j = 1

1
2
, if j − i = 1

0, other wise

. (13)

M

(
Un+1 − Un

∆t

)
+ aBUn+1 = 0.

Here Un denotes U at time t = tn = ∆tn,
and ∆t is the time step. Rearranging the
terms we obtain the system:(
M

∆t
+ aB

)
Un+1 =

1

∆t
MUn, n = 0, 1, 2, . . . ,

to be solved for Un+1 by using initial con-
dition for U0 = U(x, t = 0).

Assembly of the mass matrix
M in 1D

∫
Ω
ϕiϕjdx ex-

actly. Moreover, since the hats ϕi and ϕj

Figure 2: Illustration of the hat functions ϕi−1 and ϕi, in this figure φ, and their support.

East Afr. J. Biophys. Comput. Sci. (2023), Vol. 4, No. 1, 52-74

59

The next step is to descritize the sys-tem
Eq.12 in time. Here we were con-
sider finite difference approximations spe-
cially the implicit Euler (Back ward Eu-
ler) method. By using the back ward Eu-
ler scheme, the system Eq.12 results the
following system of algebraic equations:

Euler Back ward represents an im-
plicit scheme which is stable for all choices
of ∆t (?). Since the scheme is implicit, we
have to solve a system of algebraic equa-
tions at each time step.

lack common support for |i− j| > 1 only
Mi,i, Mi,i+1 , and Mi+1,i need to be calcu-
lated. All other matrix entries are zero by
default. This is clearly seen from Figure 2,
(Larson and Bengzon, 2010) showing two
neighboring hat functions and their sup-
port. As a consequence, the mass matrix
M is tridiagonal.

Let us now go through the details of how
to assemble the mass matrix M . We be-
gin by calculating the entries Mi,j of the
mass matrix, which involve products of
hat functions given in Eq. 7. Since each
hat is a linear polynomial, the product
of two hats is a quadratic polynomials.
Thus, Simpson’s formula (Eq. 18) can
be used to integrate Mi,j =

Here M and A are tridiagonal matrices and Eq.12 is a simple system of ODE.



a

f(x)dx =
f(a) + 4f(a+b

2
) + f(b)

6
(b−a).

(14)
Let we start on the diagonal entries of M ,
Mi,i using Simpson’s formula:

Mi,i =

∫
Ω

ϕ2
i dx,

=

∫ xi

xi−1

ϕ2
i dx+

∫ xi+1

xi

ϕ2
i dx,

=
0 + 4(1

2
)2 + 1

6
hi +

1 + 4(1
2
)2 + 0

6
hi+1,

=
hi

3
+

hi+1

3
, for i = 1, 2, . . . , n− 1,

where hi+1 = xi+1−xi and hi = xi−xi−1,
but in our case we use a uniform mesh
length h = xi+1 − xi = xi − xi−1 and
Mii = 2h

3
. The first and last diagonal

entry are M00 = h1

3
and Mnn = hn

3
, re-

spectively, since the hat functions ϕ0 and
ϕn are only half.

formula we have

Mi+1,i =

∫
Ω

ϕi+1ϕidx,

=

∫ xi+1

xi

ϕi+1ϕidx,

=
1.0 + 4(1

2
)2 + 0.1

6
hi+1

=
hi+1

6
, for i = 0, 1, 2, . . . , n− 1.

By using a similar calculation the super
diagonal entries Mi,i+1 have the values

Mi,i+1 =
hi+1

6
, for i = 0, 1, 2, . . . , n− 1.

Hence, the mass matrix m takes the form

M =



h1

3
h1

6
h1

6
h1

3
+ h2

3
h2

6
h2

6
h2

3
+ h3

3
h3

6
. . . . . . . . .

hn−1

6
hn−1

3
+ hn

3
hn

6
hn

6
hn

3


. (15)

The global mass matrix M can be written as a sum of n simpler elemental matrices as:

M =


h1

3
h1

6
h1

6
h1

3

+


h2

3
h2

6
h2

6
h2

3

+ . . .+

 hn

3
hn

6
hn

6
hn

3

 .

i.e., M = MΩ1 +MΩ2 + . . .+MΩn .

e

M e = h
6

[
2 1
1 2

]
, where h is the length

of e. We refer to M e as the local element
mass matrix. The summation of the ele-
ment mass matrices into the global mass
matrix is called assembling.

East Afr. J. Biophys. Comput. Sci. (2023), Vol. 4, No. 1, 52-74

60

On the interval I = (a, b), Simpson’s
formula is of the form, (Larson and Beng-
zon, 2010)∫ b We know continue with the sub di-

agonal entries Mi+1,i still using Simpson’s

, i = 1, 2, . . . , n, is
obtained by restricting the integration to
one sub interval or element Ωe and is
therefore called an element mass matrix.
From the sum we see that on each ele-
ment e this small block takes the form:

Each matrix MΩ



Implementation of the 1D Diffusion
Equation

Let we now consider the diffusion equa-
tion

∂u

∂t
= D

∂2u

∂x2
.

We multiply the equation by a test func-
tion v(x), which satisfies the boundary
conditions v(0) = v(1) = 0 and then inte-

grating by parts we have that:

∂u

∂t
.v = D

∂2u

∂x2
.v

∫ 1

0

(
∂u

∂t
.v

)
=

∫ 1

0

(
D
∂2u

∂x2
.v

)
Hence by using integration by parts, we
obtain∫ 1

0

∂u

∂t
.v = −D

∫ 1

0

∂u

∂x
.
∂v

∂x
, (16)

which is the weak formulation of the one
dimensional Diffusion equation.

∫ 1

0

∂

∂t

N−1∑
i=1

uiϕi.
N−1∑
j=1

vjϕj = −D

∫ 1

0

∂

∂x

N−1∑
i=1

uiϕi.
∂

∂x

N−1∑
j=1

vjϕj.

which then implies,

N−1∑
j=1

vj

(
∂

∂t

N−1∑
i=1

ui

∫ 1

0

ϕi.ϕj

)
=

N−1∑
j=1

vj

(
−D

N−1∑
i=1

ui

∫ 1

0

ϕ′
i.ϕ

′
j

)
.

That is
∂

∂t

N−1∑
i=1

ui

∫ 1

0

ϕi.ϕj = −D

N−1∑
i=1

ui

∫ 1

0

ϕ′
i.ϕ

′
j.

In a matrix form it can be written as:

MU̇ +DAU = 0, (17)

Assembly of the stiffness matrix in 1D

The stiffness matrix A is symmetric for this simple problem, which makes the computa-
tion of the matrix faster since we don’t have to compute all of the elements, symmetric
matrices are also much faster to invert. Here ϕi’s are the hat functions given in Eq. 7,
the entries of each element of the stiffness matrix A is given by

East Afr. J. Biophys. Comput. Sci. (2023), Vol. 4, No. 1, 52-74

61

Now the system of equation, that is
Eq.12 is a simple system of ordinary dif-
ferential equations. To solve this system
of ode’s, we have to use a back ward Euler
method and a Matlab soft ware to solve
the system at each time steps using an
initial condition.

where M is the mass matrix with entries given in Eq. 15 and A is the stiffness matrix
with entries given in Eq. 19.

Here the matrix A is often referred to as the stiffness matrix, a name coming from
corresponding matrices in the context of structural problems.

Now substituting Eq.10 and Eq.11
in the weak formulation of the equation
Eq.16, we have:



Ai,j =

∫ 1

0

ϕ′
iϕ

′
jdx,

=
N−1∑
e=1

∫
Ωe

ϕ′
iϕ

′
jdx,

=
N−1∑
e=1

Ae
ij.

similarly the load vector

Fi =

∫ 1

0

fϕidx =
N−1∑
e=1

∫
Ωe

fϕidx =
N−1∑
e=1

F e
i .

a

f(x)dx =
f(a) + 4f(a+b

2
) + f(b)

6
(b− a). (18)

Then we can be illustrate by the hat functions and Simpsons formula as follows:

ϕi(x) =


x−xi−1

hi
, if xi−1 ≤ x ≤ xi

xi+1−x
hi+1

, if xi ≤ x ≤ xi+1

0, other wise

, ϕi−1(x) =


x−xi−2

hi−1
, if xi−2 ≤ x ≤ xi−1

xi−x
hi

, if xi−1 ≤ x ≤ xi

0, other wise

,

and ϕi+1(x) =


x−xi

hi+1
, if xi ≤ x ≤ xi+1

xi+2−x
hi+2

, if xi+1 ≤ x ≤ xi+2

0, other wise

.

Now

Ai,i = a(ϕi, ϕi)
e =

∫
Ωe

ϕ′
iϕ

′
idx,

=

∫ xi

xi−1

ϕ′
iϕ

′
idx+

∫ xi+1

xi

ϕ′
iϕ

′
idx,

=

∫ xi

xi−1

1

hi

1

hi

dx+

∫ xi+1

xi

−1

hi+1

−1

hi+1

dx,

=
hi

6

(
1

h2
i

+
4

h2
i

+
1

h2
i

)
+

hi+1

6

(
1

h2
i+1

+
4

h2
i+1

+
1

h2
i+1

)
,

=
1

hi

+
1

hi+1

.

,

Ai−1,i = a(ϕi−1, ϕi)
e =

∫
Ωe

ϕ′
i−1ϕ

′
idx,

=

∫ xi

xi−1

ϕ′
i−1ϕ

′
idx+

∫ xi+1

xi

ϕ′
i−1ϕ

′
idx,

=

∫ xi

xi−1

−1

hi

1

hi

dx,

=
hi

6

(
−1

h2
i

+
−4

h2
i

+
−1

h2
i

)
=

−1

hi

.

East Afr. J. Biophys. Comput. Sci. (2023), Vol. 4, No. 1, 52-74

62

On the interval I = (a, b) Simpson’s formula is of the form, (Larson and Bengzon, 2010)∫ b



and

Ai,i+1 = a(ϕi+1, ϕi)
e =

∫
Ωe

ϕ′
i+1ϕ

′
idx,

=

∫ xi

xi−1

ϕ′
i+1ϕ

′
idx+

∫ xi+1

xi

ϕ′
i+1ϕ

′
idx,

=

∫ xi+1

xi

−1

hi+1

1

hi+1

dx,

=
hi+1

6

(
−1

h2
i+1

+
−4

h2
i+1

+
−1

h2
i+1

)
=

−1

hi+1

.

Each generic interior element contributes to the stiffness matrix of a 2× 2 sub matrix.

A =

∫ 1

0

ϕ′
iϕ

′
jdx =

N−1∑
e=1

Ae =



1
h1

−1
h1−1

h1

1
h1

+ 1
h2

−1
h2

−−1
h2

1
h2

+ 1
h3

−1
h3

. . . . . . . . .
−1

hn−1

1
hn−1

+ 1
hn

−1
hn

−1
hn

1
hn


.

(19)
The global stiffness matrix A can be written as a sum of n simpler elemental matrices as:

A = 1
h1


1 −1
−1 1

+ 1
h2


1 −1
−1 1

+ . . .+ 1
hn

 1 −1
−1 1

 .

i.e., A = AΩ1 + AΩ2 + . . .+ AΩn .
Each matrix AΩe , e = 1, 2, . . . , n, is ob-
tained by restricting the integration to
one sub interval or element Ωe and is
therefore called an element stiffness ma-
trix. From the sum we see that on each
element e this small block takes the form:

Ae = 1
h

[
1 −1
−1 1

]
, where h is the length

of e. We refer to Ae as the local element
stiffness matrix.

Implementation of One Dimensional
Advection Diffusion equation

Let us now solve the 1D governing (ad-
vection diffusion) equation

∂u

∂t
+ a

∂u

∂x
= D

∂2u

∂x2
+ f, u(0) = u(1) = 0,

by using FEM. Let we find the weak for-
mulation of the equation by multiplying
the equation with a test function v(x),
which satisfies the boundary conditions
v(0) = v(1) = 0 and then integrating by
parts as the same procedure above. we
have that:

∂u

∂t
.v + a

∂u

∂x
.v = D

∂2u

∂x2
.v + f.v,∫ 1

0

(
∂u

∂t
.v + a

∂u

∂x
.v

)
=

∫ 1

0

(
D
∂2u

∂x2
.v + f.v

)
.

Using integration by parts, we have then

East Afr. J. Biophys. Comput. Sci. (2023), Vol. 4, No. 1, 52-74

63

Now the system of equation, that is
Eq.17 also is a simple system of ordinary
differential equations and we can solve
this system of Ode’s to get the solution
of the original PDE.



∫ 1

0

∂u

∂t
.v +

∫ 1

0

a
∂u

∂x
.v = −D

∫ 1

0

∂u

∂x
.
∂v

∂x
+

∫ 1

0

f.v, (20)

∫ 1

0

∂

∂t

N−1∑
i=1

uiϕi.
N−1∑
j=1

vjϕj +

∫ 1

0

a
∂

∂x

(
N−1∑
i=1

uiϕi

)
.
N−1∑
j=1

vjϕj

= −D

∫ 1

0

∂

∂x

N−1∑
i=1

uiϕi.
∂

∂x

N−1∑
j=1

vjϕj +

∫ 1

0

f.

N−1∑
j=1

vjϕj.

which then implies,

N−1∑
j=1

vj

(
∂

∂t

N−1∑
i=1

ui

∫ 1

0

ϕi.ϕj + a
N−1∑
i=1

ui

∫ 1

0

ϕ′
i.ϕj

)
=

N−1∑
j=1

vj

(
−D

N−1∑
i=1

ui

∫ 1

0

ϕ′
i.ϕ

′
j +

∫ 1

0

f.ϕj

)
.

(21)
That is

∂

∂t

N−1∑
i=1

ui

∫ 1

0

ϕi.ϕj + a
N−1∑
i=1

ui

∫ 1

0

ϕ′
i.ϕj = −D

N−1∑
i=1

ui

∫ 1

0

ϕ′
i.ϕ

′
j +

∫ 1

0

f.ϕj.

In a matrix form it can be written as:

MU̇ + aBU +DAU = F, (22)

Assembly of the Load Vector
in 1D

the interval [a, b] with h = b−a
N

, xi = a+ih
for each i = 0, 1, . . . , N :

∫ b

a

f(x)dx =
h

3

f(x0) + 2

N
2
−1∑

j=1

f(x2j) + 4

N
2∑

j=1

f(x2j−1) + f(xN)

 , j = 1, 2, ..., (
N

2
)− 1.

(23)

East Afr. J. Biophys. Comput. Sci. (2023), Vol. 4, No. 1, 52-74

64

which is the weak formulation of
the one dimensional advection-Diffusion
equation. Substituting Eq.10 and Eq.11

where M is the mass matrix with entries
given in Eq. 15, B is a matrix with entries
given in Eq. 13, A is the stiffness matrix
given in Eq. 19 and F is a load vector
given in Eq. 25.

The right-hand-side, load vector of Eq. 20
contains an integral over a function f(x).

In general, exactly computing this inte-
gral is very difficult, so another numerical
approximation is required. We can use
a well known integration rule composite
simpson rule to approximate these inte-
gration whose formula (for more informa-
tion you can see, (?)) is given in Eq. 23,
by selecting a set of distinct N nodes in

in the weak formulation of the equation
Eq.20, we have:



∫ b

a

f(x)dx =
f(a) + f(b)

2
(b− a), (24)

we have

fi =

∫
I

fϕidx,

=

∫ xi+1

xi−1

fϕidx,

=

∫ xi

xi−1

fϕidx+

∫ xi+1

xi

fϕidx,

≈ f(xi−1)ϕi(xi−1) + f(xi)ϕi(xi)

2
hi +

f(xi+1)ϕi(xi+1) + f(xi)ϕi(xi)

2
hi+1,

=
0 + f(xi)

2
hi +

f(xi) + 0

2
hi+1,

= f(xi)

(
hi

2
+

hi+1

2

)
.

Now using this trapezoidal method,
the approximate load vector takes the
form

F =



f(x0)
h1

2

f(x1)
(
h1+h2

2

)
f(x2)

(
h2+h3

2

)
...

f(xn−1)
(

hn−1+hn

2

)
f(xn)

hn

2


. (25)

Splitting F into a sum over the elements
yields the n global element load vectors
FΩe :

F =


f(x0)
f(x1)


h1

2
+


f(x1)
f(x2)

 h2

2
+

 f(x2)
f(x3)

 h3

2
+ . . .+

 f(xn−1)
f(xn)


hn

2

i.e., F = FΩ1 + FΩ2 + . . . + FΩn . Each

vector FΩe , e = 1, 2, . . . , n, is formally de-
rived by restricting the integration to ele-
ment Ωe.

MU̇ + (aB +DA)U = F, (26)

which is a simple system of ordinary dif-
ferential equations. For solving this sys-
tem of Ode’s, we have to use a Matlab soft
ware in which it has a number of tools for
numerically solving ordinary differential
equations. We would focus on the back
ward Euler method to descritize time.

FEM implementation of the
Two Dimensional AD equa-
tion

The 2D advection diffusion equation with
the same and constant velocity and diffu-
sion term is given by

ut+a(ux+uy) = D(uxx+uyy)+f, 
(x, y)  ∈ Ω = [0, 1],

East Afr. J. Biophys. Comput. Sci. (2023), Vol. 4, No. 1, 52-74

65

Using the another quadrature rule for simplicity, for instance, using the Trapezoidal rule,
(Larson and Bengzon, 2010)

Now the system of equation, that is
Eq.22 can be written in the form:



with homogeneous boundary conditions.
That is

ut + a∇u = D∇2u+ f.

Now to find a weak formulation for this
2D equation, we multiply both sides of
the equation with a test function v =
v(x, y) ∈ V which satisfies the boundary
conditions.

ut.v + a∇u.v = D∇2u.v + fv.

Here integrating this over the domain Ω
yields the following:

∫
Ω

(ut.v + a∇u.v) =

∫
Ω

(
D∇2u.v + fv

)
.

We see from the 2D FEM of Poisson equa-
tion (using Gauss theorem and the trans-
formation of a surface integral to a line in-
tegral) that

∫
Ω
v∇2u = −

∫
Ω
∇v∇u, and

hence we get:

∫
Ω

(ut.v + a∇u.v) = −D

∫
Ω

(∇v∇u) +

∫
Ω

fv.∫
Ω

ut.v = −
∫
Ω

(D∇v∇u+ a∇u.v) +

∫
Ω

fv. (27)

(ut, v) = l(u, v) + (f, v) ∀v ∈ V, (28)

where (ut, v) =
∫
Ω
ut.v and l(u, v) = −

∫
Ω
(D∇v∇u+ a∇u.v).

Given a FE space V , with ϕi(x, y), i = 1, 2, . . . , N denoting a set of basis functions
for V , we seek the FE solution of form

uh(x, y, t) =
N∑
j=1

uj(t)ϕj(x, y). (29)

And taking the test function v(x, y) as a linear combination of basis functions

v(x, y) =
N∑
j=1

vjϕj, with vj are constants. (30)

∑
j=1

u′
j(t)ϕj(x, y),

N∑
j=1

vjϕj

)
= l

(
N∑
j=1

uj(t)ϕj(x, y),
N∑
j=1

vjϕj

)
+ (f,

N∑
j=1

vjϕj). (31)

Then we get the linear system of ordinary differential equations in the uj(t) as:(
u′
1(t)ϕ1(x, y) +

N∑
j=2

u′
j(t)ϕj(x, y),

N∑
j=1

vjϕj

)
= l

(
u1(t)ϕ1(x, y) +

N∑
j=1

uj(t)ϕj(x, y),
N∑
j=1

vjϕj(x, y)

)
+

(f,
N∑
j=1

vjϕj(x, y)),

East Afr. J. Biophys. Comput. Sci. (2023), Vol. 4, No. 1, 52-74

66

Substituting this expression (29 and 30) into eq.28, we obtain(
N

Here Eq. 27 is the weak formulation and it can be simplified as



(
u′
1(t)ϕ1,

N∑
j=1

vjϕj

)
+

(
N∑
j=2

u′
j(t)ϕj,

N∑
j=1

vjϕj

)
= l

(
u1(t)ϕ1,

N∑
j=1

vjϕj

)
+

l

(
N∑
j=1

uj(t)ϕj,

N∑
j=1

vjϕj

)
+ (f,

N∑
j=1

vjϕj),

N∑
j=1

(
u′
j(t)ϕj,

N∑
j=1

vjϕj

)
=

N∑
j=1

l

(
uj(t)ϕj,

N∑
j=1

vjϕj

)
+ (f,

N∑
j=1

vjϕj),

N∑
j=1

N∑
j=1

(
u′
j(t)ϕj, vjϕj

)
=

N∑
j=1

N∑
j=1

l (uj(t)ϕj, vjϕj) +
N∑
j=1

(f, vjϕj).

Since vj’s are constants we have also that:

N∑
j=1

vj

N∑
j=1

(ϕj, ϕj)u
′
j(t) =

N∑
j=1

vj

N∑
j=1

l (ϕj, ϕj)uj(t) +
N∑
j=1

vj(f, ϕj).

The corresponding problem can therefore be expressed as

V TMU̇ = V TAU + V TF.

That is
MU̇ = AU + F. (32)

Where, M =


(ϕ1, ϕ1) (ϕ1, ϕ2) . . . (ϕ1, ϕN)
(ϕ2, ϕ1) (ϕ2, ϕ2) . . . (ϕ2, ϕN)

...
...

...
(ϕN , ϕ1) (ϕN , ϕ2) . . . (ϕN , ϕN)

,

A =


l(ϕ1, ϕ1) l(ϕ1, ϕ2) . . . l(ϕ1, ϕN)
l(ϕ2, ϕ1) l(ϕ2, ϕ2) . . . l(ϕ2, ϕN)

...
...

...
l(ϕN , ϕ1) l(ϕN , ϕ2) . . . l(ϕN , ϕN)

, and F =


(f, ϕ1)
(f, ϕ2)

...
(f, ϕN)

 .

In this section, we are compared for the 
advection-diffusion equations with an ex-
act solution for the given finite element 
methods and then we solve the equa-
tion with out knowing the exact solution. 
The comparison is carried out by means 
of computed solutions for a wide range

East Afr. J. Biophys. Comput. Sci. (2023), Vol. 4, No. 1, 52-74

67

There are many methods to solve the
above problem involving the system of
first order ODE. We can use FD meth-
ods that will descritize in time by us-
ing Explicit Euler method, Implicit Euler
method or the Crank-Nicolson method,
by considering an appropriate initial con-
dition, (Johnson, 2012). But for this pa-
per in the two dimensional case, we can
use the ODE Suite in Matlab which is the
Matlab build in system of ODE solver,
ODE15I.

RESULTS AND DISCUSSION



of characteristic parameters. Linear el-
ements are employed at the discretiza-
tion in case of one-dimensional problem
and bilinear elements in case of two-
dimensional problems. Dirichlet and gen-
eral boundary conditions are considered
with different initial conditions and dif-
ferent size of computational domain.

x2) + t cos(πx)(1 − 2x) and the veloc-
ity parameter a = 3, with homogeneous

boundary conditions. The FM solution
using Matlab with a back ward Euler des-
critization in time is given in Fig. 3(a),
with N = 100 nodes. To saw our er-
ror the exact solution for this equation is
u(x, t) = t cos(πx)(x−x2) and its graph is

of nodes from N = 100 to N = 1000, then
our numerical solution becomes more ac-
curate and we saw that the error is an or-
der of 10−7 by modifying h = 0.001 from
the algorithm.

(a) FEM solution (b) Exact solution

Figure 3: The Matlab implementation of the advection dominated equation with f(x, t) =
cos(πx)− atπ sin(πx)(x− x2) + t cos(πx)(1− 2x).

tively large diffusion coefficient, the dif-
fusion process is faster. Hence our finite
element method is reasonable and accu-
rate.

East Afr. J. Biophys. Comput. Sci. (2023), Vol. 4, No. 1, 52-74

68

Fig. 3(b). The error of this equation is an
order of 10−4. If we increase the number

In (Bergara, 2011), there is a diffusion
(ut = Duxx) example solved with finite
difference methods, and let we solve that
equation with the finite element method.
He solves the diffusion (heat) equation by
using initial condition u(x, 0) = sin(πx)

First, consider the advection equation
ut + aux = f with sources and sinks func-
tion, f(x, t) = cos(πx) − atπ sin(πx)(x −

When we saw the results of those fig-
ures (Fig. 4 to Fig. 6) we can observe
the following. If we use a small amount
of diffusion coefficient, then it has a mall
diffusion process and if we consider rela-

and with homogeneous Dirichlet bound-
ary conditions in the interval 0 ≤ x ≤ 1.
If we consider different values for final
time and diffusion coefficient we get the
following simulations (Fig. 4 to Fig. 6)
with a similar descritization of space (x)
in to 50 nodes.



(a) FEM sol. with D = 0.1 and tf = 0.5 (b) FEM sol. with D = 1 and tf = 0.5

Figure 4: Solution of the diffusion equation using D = 0.1 and D = 1 for a constant final
time of tf = 0.5.

(a) FEM sol. with D = 10 and tf = 0.5 (b) FEM sol. with D = 0.1 and tf = 3

Figure 5: Solution of the diffusion equation using D = 10 with tf = 0.5 and D = 0.1 for
final time of tf = 3.

sin(πx) (1 +Dπ2t) + aπt cos(πx),
u(x, 0) = 0, u(0, t) = u(1, t) = 0. By
using the forward Euler descritization in

Pe = al
D

= 0.003 ≪ 1 which is diffusion
dominated for this choice of a and D.

East Afr. J. Biophys. Comput. Sci. (2023), Vol. 4, No. 1, 52-74

69

time the Matlab implementation with
a = 0.03 and D = 10 for this problem
is given in Fig. 7. The peclet’s number

Let we know consider the ADE
ut + aux − Duxx = f with f(x, t) =

Since we implements zero flux con-
centrations of the pollutants in the bound-
aries, the total mass of the pollution
should remain constant. These bound-
aries physically correspond to a system
where the species is enclosed inside a mesh
that it cannot penetrate, however, the
mesh allows the fluids to diffuse through

the flow field. We see that all masses even-
tually concentrate along the domain as we
have seen in Fig. 7. This steady state cor-
responds to diffusive and advective fluxes
balancing each other. The flow carries ad-
ditional mass towards the domain, but the
density gradient limits how much more
mass can be deposited.



(a) FEM sol. with D = 1 and tf = 3 (b) FEM sol. with D = 10 and tf = 3

Figure 6: Solution of the diffusion equation using D = 1 and D = 10 for a constant final
time of tf = 0.5.

(a) FEM solution (b) Bar plot of the solution

Figure 7: Implementation of 1D AD equation with the source functionf(x, t) =
sin(πx) (1 +Dπ2t) + aπt cos(πx) and homogeneous Dirichlet boundary conditions and
a = 0.03, D = 10, in which the peclet’s number Pe = al

D
= 0.003 ≪ 1 which is diffusion

dominated.

ary condition Dux(0, t) = k0(u(0, t) + g0
and −Dux(1, t) = k1(u(1, t) + g1 and by
considering different robin boundary pa-

East Afr. J. Biophys. Comput. Sci. (2023), Vol. 4, No. 1, 52-74

70

We see know on an arbitrary in-
put concentration by using Neumann and
Robin boundary conditions. Here, the
problem is undetermined because of an
unknown input concentration and hence
unknown exit concentration and depends
on the parameters. Authors of previous
works on problems of this type have done
on a known exit concentration by assum-
ing a continuous and constant concen-
tration at the flow boundary of Dirich-
let type. This yields we to consider a

problem by forcing the flow boundaries
using Neumann and Robin boundary con-
ditions. This conditions yields the flow of
the concentration to move freely. Here the
velocity term and the diffusion coefficient
highly affects the flow of the concentra-
tion. We saw it by giving a source func-
tion f(x, t) = a(2x + 1) − 2D − 1, 0 ≤
x ≤ 1, 0 ≤ t ≤ tf with a robin bound-



rameters.

2

(a) FEM solution with tf = 1 (b) FEM solution with tf = 3

Figure 8: Numerical solution using FEM for the advection diffusion equation with
f(x, t) = a(2x + 1) − 2D − 1 and robin boundary conditions, with a diffusion coeffi-
cient D = 0.02 and velocity term a = 1 with time tf = 1 for (a) and tf = 3 for (b).

u(t, x, y) = e−t sin(πx) sin(πy) be the
given exact solution for the diffusion equa-
tion ut − D∆u = f , then our source

function f becomes f(t, x, y) = (2Dπ2 −
1)e−t sin(πx) sin(πy). Now by using the

East Afr. J. Biophys. Comput. Sci. (2023), Vol. 4, No. 1, 52-74

71

Let we change the velocity term and
the diffusion coefficient for the given
source term using D = 0.02 and a = 1,
(g0, g1) = (0.02, 0.06). The solution for

sion equation with a source function given
as f(x, y) = Dπ2(sin(πx)+sin(πy)) to im-

Let we consider the 2D steady diffu-

In order to consider the finite el-
ement method of the time dependent
two dimensional equations, we use the
unsteady 2D diffusion equation. Let
we test it with a known exact function
and use the source terms and bound-
ary conditions from that function. Let

plement the 2D problem. The FEM nu-
merical simulation of the PDE is given in
Fig. 9.

ODE15I, ODE solver to integrate for time
in the final finite element method descriti-
zation the solutions for the equation is
given in Fig. 10. The surface plots of
those figures are the exact solutions, the
finite element method solutions and the
color bar of the finite element solution to
view its properties.

this conditions is given in Fig. 8(a) us-
ing a total time of tf = 1 and 8(b) using
a total time of tf = 3. This is an ad-
vection dominated example in which the
movement of the concentration is faster in
advection terms relative to diffusion.



(a) FEM solution (b) Bar plot of the solution

Figure 9: The numerical simulation using FEM for the steady 2D diffusion equation with
f(x, y) = Dπ2(sin(πx) + sin(πy)), with a diffusion coefficient D = 0.5.

(a) FEM & exact sol (b) Colour bar of the FEM solution

Figure 10: The surface of time dependent two dimensional diffusion equation with D =
100 with its color bar plot for its source term f , f(t, x, y) = (2Dπ2−1)e−t sin(πx) sin(πy).

method to solve any differential equations
numerically in any type of geometries.

East Afr. J. Biophys. Comput. Sci. (2023), Vol. 4, No. 1, 52-74

72

If we consider another linear test ex-
ample with exact solution u(x, y, t) = x+
y + t, then its source function becomes
f(x, y, t) = 1. The finite element solution

The numerical results in all the above
mentioned discussions using the finite el-
ement method are almost closed to the
analytical solution in the case of the test
examples. This shows us that the finite el-
ement method is one of the best numerical

for this 2D unsteady diffusion equation
with its exact solution and color bar is
given in Fig. 11 by ODE solver ODE15I.



(a) FEM & exact sol (b) Colour bar of the FEM solution

Figure 11: The surface of time dependent two dimensional diffusion equation with D = 1
with its color bar plot for its source term f , f(t, x, y) = 1.

sional and two dimensional cases with
the weighted residual (Galerkin) method
of finite elements with constant veloc-
ity term and diffusion coefficient. These
statements are supported by our test nu-
merical investigations for the one and two
dimensional Poisson equation, the advec-
tion equation, diffusion equation and the
advection-diffusion equation. We have
also seen the simulations of the equa-
tions by using the language of technical
computing called Matlab and for time
dependent equation using backward (im-
plicit) Euler finite difference method for
one dimensional problem and the built-
in function ode15i to solve the system of
ode’s for two dimensional equations.

East Afr. J. Biophys. Comput. Sci. (2023), Vol. 4, No. 1, 52-74

73

The generalization of the proposed
Galerkin method to the three-dimensional
advection diffusion is obtained within the
current framework and an interested body
will be done on this approach by including
the reaction term. To obtain the solutions
to this problems, the method can be ex-
tended to the least square finite element
method, and another interested body will
also be done on this approach. One can
also extend this method using variable
velocities and diffusion coefficients whichWe have solved the advection-

diffusion e quation f or b oth o ne dimen-

In this paper, we study the finite ele-
ment solution of the advection-diffusion 
equation. In the finite e lement analy-
sis, we approximate a function defined 
in a domain, Ω, with a set of orthogo-
nal basis functions with coefficients corre-
sponding to the functional values at some 
node points. We deal with the numeri-
cal simulation of the advection-diffusion 
equation using the finite e lement method 
scheme in space and the backward Euler 
method in time. The main focus of the 
paper has been on the variational formu-
lation techniques for the solution of the 
discrete Galerkin method and the other 
hand on computational analysis of differ-
ent one and two dimensional PDE’s. We 
have done numerical simulations for the 
above mentioned equations by consider-
ing the technique and using different test 
examples. The solution for the values at 
the nodes for the partial differential equa-
tions can be obtained by solving a linear 
system of equations involving the inver-
sion of the sparse matrices.

CONCLUSIONS 



East Afr. J. Biophys. Comput. Sci. (2023), Vol. 4, No. 1, 52-74

Hundsdorfer W.1996. Numerical solution of advection-
diffusion-reaction equations: lecture notes for Ph.
D. course, 1996, Thomas Stieltjes Institute.
Department of Numerical Mathematics [NM],
(9603).

Johnson C. 2012. Numerical solution of partial
differential equations by the finite element
method. Courier Corporation.

Langtangen H.P. 2003. Computational partial
differential equations: numerical methods and
diffpack programming (Vol. 2). Berlin: Springer.

Larson M.G. and Bengzon F. 2010. The finite element
method: theory, implementation, and practice.
Texts in Computational Science and Engineering,
10, 23-44.

Lian Y., Ying Y., Tang S., Lin S., Wagner G.J. and Liu
W.K. 2016. A Petrov–Galerkin finite element
method for the fractional advection– diffusion
equation. Comput. Methods Appl. Mech. Eng.,
309: 388-410.

Lima S.A., Kamrujjaman M. and Islam M.S. 2021.
Numerical solution of convection–diffusion–
reaction equations by a finite element method with
error correlation. AIP Advances 11(8): 085225.

Pochai N. and Deepana R. 2011. A numerical
computation of water quality measurement in a
uniform channel using a finite difference method.
Procedia Eng. 8: 85-88.

Quarteroni A. and Quarteroni S. 2009. Numerical
models for differential problems (Vol. 2). Milan:
Springer.

Szymkiewicz R. and Gąsiorowski D. 2021. Adaptive
method for the solution of 1D and 2D advection–
diffusion equations used in environmental
engineering. J. Hydroinformatics 23(6): 1290-
1311.

Yang W.Y., Cao W., Chung T.S. and Morris J.
2005. Applied Numerical Methods Using
MATLAB. A John Wiley & Sons. Hoboken, NJ.

vary with the time and space of the given
dimension. The finite element method de-
scretizes only in space and the time de-
scritization is  based on the finite differ-
ence methods. Hence the finite element
method needs human power still  now for
the descritization process of time. So this
is a very interesting area to do researches
in the future to include the time descriti-
zation in the method.

Ahsan M. 2012. Numerical solution of the advection-
diffusion equation using Laplace transform finite
analytical method. Int. J. River Basin Manag.
10(2): 177- 188.

Aragonés L., Pagán J.I., López I., Navarro- González
F.J. and Villacampa Y. 2019. Galerkin's
formulation of the finite elements method to
obtain the depth of closure. Sci. Total Environ.
660: 1256-1263.

Bajellan A.A.F. 2015. Computation of the convection
diffusion equation by the fourth order compact
finite difference method. Doctoral dissertation,
Izmir Institute of Technology, Turkey.

Bergara A. 2011. Finite Difference Numerical
Methods of Partial Differential Equations in
Finance with MATLAB. Master and Banca. New
York, Springer.

Brenner S.C., Scott L.R., & Scott L.R. 2008. The
mathematical theory of finite element methods.
Vol , pp 263-291. Springer, New York.

Donea J. and Huerta A. 2003. Finite element methods
for flow problems. John Wiley & Sons.

74

The  core of this work was completed dur-
ing the master study at Hawassa Uni-
versity and presented to   former School
of Mathematics and Statistical Science
(currently Department of Mathematics)
in partial fulfillment of the requirement
for the Degree of Master. The authers
would like to thank the Department of
Mathematics at the University  of Hawassa
University. This work of  has been par-
tially supported by Hawassa University.

References

Acknowledgments


	Juliet and Florence _ Beef handling practices at Abattoir and Bucher shops_DOI.pdf (p.1-17)
	Amanuel_ Study of Reaction Mechanisms in α + 69Ga reaction at ≈10 – 50 MeV_DOI.pdf (p.18-27)
	Buchale and Atnafu_Population dynamics and yield estimation of Common Carp.pdf (p.28-41)
	Haftom et al_Synthesis Characterization and Antibacterial activity of Benzimidazole.pdf (p.42-51)
	Kassahun and Zerihun_ Numerical Solutions of Advection Diffusion Equations.pdf (p.52-74)

