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, (Aragoné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)