\documentclass[reqno]{amsart}
\usepackage{hyperref}
\usepackage{graphicx}
\AtBeginDocument{{\noindent\small
\emph{Electronic Journal of Differential Equations},
Vol. 2017 (2017), No. 43, pp. 1--11.\newline
ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu}
\thanks{\copyright 2017 Texas State University.}
\vspace{8mm}}
\begin{document}
\title[\hfilneg EJDE-2017/43\hfil A semi-analytic spectral method]
{A semi-analytic spectral method for elliptic partial differential equations}
\author[I. Ali, M. T. Saleem \hfil EJDE-2017/43\hfilneg]
{Ishtiaq Ali, Maliha Tahseen Saleem}
\address{Ishtiaq Ali \newline
Department of Mathematics,
COMSATS Institute of Information Technology,
Park Road, Chak Shahzad, Islamabad 44000, Pakistan}
\email{ishtiaqali@comsats.edu.pk}
\address{Maliha Tahseen Saleem \newline
Department of Mathematics,
COMSATS Institute of Information Technology,
Park Road, Chak Shahzad, Islamabad 44000, Pakistan}
\email{malihasaleem@yahoo.com}
\dedicatory{Communicated by Giovanni Molica Bisci}
\thanks{Submitted October 26, 2016. Published February 10, 2017.}
\subjclass[2010]{35J25, 65N35}
\keywords{Semi-analytical technique; Chebyshev-spectral method;
\hfill\break\indent exponential matrix}
\begin{abstract}
In this article we present a semi-analytic method for solving
elliptic partial differential equations. The technique is based on
using a spectral method approximation for the second-order derivative
in one of the spatial directions followed by solving the resulting
system of second-order differential equations by an analytic method.
That is, the system of second-order two-point boundary-value problems
are solved analytically by casting them in first-order form and
solving the resulting set of first-order equations by using the matrix
exponential. An important aspect of our technique is that the solution
obtained is semi-analytic, e.i., using analytic method in $y$ and
spectral method in $x$. The new method can be used for both linear and
non-linear boundary conditions as well as for nonlinear elliptic
partial differential equations.
\end{abstract}
\maketitle
\numberwithin{equation}{section}
\allowdisplaybreaks
\section{Introduction}\label{sec:INTRODUCTION}
Elliptic partial differential equations (PDEs) are frequently used to model
a variety of engineering phenomena, such as steady-state heat conduction
in a solid, or reaction-diffusion type problems. However, computing a
solution of these type PDEs can sometimes be difficult or inefficient
using standard solvers. Techniques have been developed, including the
method of lines, which can solve parabolic PDEs using well developed
numerical solvers, but are not directly applicable to elliptic PDEs.
Laplace equation is the well known example of elliptic PDEs.
There are many numerical techniques to solve the Laplace equations.
Some of them are implicit alternative direction method, over relaxation
method \cite{c1}. Schiesser \cite{s1,s2} introduce the false transients
method in which a time derivative of dependent variable gets add to
the Laplace equation and then spatial derivatives are approximated
by the the finite differences, after that obtained system of equations
are solved by the method of lines \cite{c3,c5,d1,r1,s1,s2,s3,s6,s7,t1}.
Laplace equations have solutions that are well known as harmonic functions.
They are analytic in the domain where the equation is satisfied and if
there exists two functions that are solutions of a Laplace equation then
there sum is also a solution of that Laplace equation.
In this article, a semi-analytical method is used to solve the Laplace equation.
Spectral method is used to approximates the second-order derivatives in one
of the spacial direction and the resultant system is then solved by using
any analytical method suitable for a system of second order differential
equations. That is, a system of two points second order boundary-value
problems can analytically be solved by converting them into first-order by
using the exponential matrix. Spectral method are global method, generally
a faster method and is more accurate for sufficiently regular geometries
than other two methods \cite{p1}. In dealing with PDEs let say time dependent,
the solution with spectral method is obtained by writing it in the summation
of the basis function. There are several types of spectral method such
as Galerikn spectral method deals with the PDEs that has coefficient of
this expression. A spectral collocation method that deals with the direct
usage of grid points hence it considered as similar to finite difference
method \cite{l1,v1}. More details about spectral methods can be find in
\cite{s4,s5}.
The collocation method is the most useful method of them all in a sense that
it deals with the non linear terms more easily then any other methods.
Our semi-analytical technique is applicable for both non linear and linear
problems having boundary conditions at $y=0$ and $y=1$.
\section{Semi-analytical spectral method for elliptic PDEs}
\label{sec:Semianalytical method of lines for linear elliptic PDEs:}
Let there be a heat transfer in rectangle of height $H$ and length $L$.
For dimensionless temperature distribution the governing equations could be
\cite{c2,s1}.
\begin{equation} \label{e2.1}
\epsilon^2\frac{\partial^2 u}{\partial x^2}+\frac{\partial^2 u}{\partial y^2}=0 ,
\end{equation}
with boundary conditions
\begin{equation} \label{e2.2}
\begin{gathered}
u(0,y)=0 ,\quad\text{for }0\leq y\leq 1,\\
u(1,y)=0, \quad\text{for }0\leq y\leq 1,\\
u(x,0)=0, \quad\text{for }0\leq x\leq 1,\\
u(x,1)=\sinh(\epsilon \pi) \sin(\pi x) . \quad\text{for }0\leq x\leq 1.
\end{gathered}
\end{equation}
Here $\epsilon$ is the aspect ratio and its considered value is
$\epsilon =H/L$. The analytical solution of this problem is \cite{s1}.
\begin{equation} \label{e2.3}
u=\sinh(\epsilon \pi y) \sin(\pi x).
\end{equation}
First of all, discritizing the spatial coordinate in $x$ direction
i.e. $\partial^2 u/ \partial x^2$ by a Chebyshev approximation to gives.
\begin{equation}
\frac{\partial^2 u}{\partial y^2}
=- \epsilon^2 \Big[\sum^{N}_{l=0} \left(D^2_N\right)_{jl} \ u_l(y)\Big] ,
\end{equation}
here $N$ is known as interior node points and $j$ is from $1$ to $N-1$.
The boundary conditions are transformed into
\begin{equation}
\begin{gathered}
u_0=0,\\
u_{N+1}=0,\\
u(y=0)=0,\\
u(y=1)=\sinh(\epsilon \pi) \sin(\pi i x)\quad\text{for }i=1,\dots,N.
\end{gathered}
\end{equation}
Just for our convenience let $\xi=y\epsilon$ which results in the system
of equation
\begin{equation}
\frac{d^2 u_i}{d \xi^2} =\sum^{N}_{l=0} \left(D^2_N\right)_{jl} u_l(y),
\end{equation}
with
\begin{equation}
\begin{gathered}
u_0=0,\\
u_{N+1}=0,\\
u(\xi=0)=0,\\
u(\xi =\epsilon)=\sinh(\epsilon \pi) \sin(\pi i x).
\end{gathered}
\end{equation}
It is difficult to handle the $N$ second order equations as above so one
can transform them to $2N$ first order equations as follows \cite{r1,s8}
\begin{gather}
\frac{du_i}{d\xi} = u_{N+1+i}, \\
\frac{du_{N+1+i}}{d \xi} =\sum^{N}_{l=0} \left(D^2_N\right)_{jl} \ u_l(y),
\end{gather}
similarly the boundary equations should be transformed for $\xi= 0,\dots,2N$ as
\begin{equation}
\begin{gathered}
u_0=0,\\
u_{N+1}=0,\\
u(\xi=0)=0,\\
u_{N+1+i}(\xi=0)=c_i,\\
u(\xi =\epsilon)=\sinh(\epsilon \pi) \sin(\pi i x).
\end{gathered}
\end{equation}
Here the constants $c_i$ are found by integrating and by using the
boundary condition at $y=1$ in it. The required form of linear first order
$2N$ system of equations are obtained in the matrix form
\begin{equation}
\frac{dY}{d\xi}=AY+b(\xi),
\end{equation}
considering
\begin{equation}
Y=[u_1,u_2,\dots,u_N,u_{N+2},u_{N+3},\dots,u_{2N+1}]^T.
\end{equation}
where $u_{N+1}$ corresponds the boundary condition at $(x=1)=0$.
Also $A$ is a coefficient matrix that is symmetric and of order $2N$
having matrices $0$ (zeros matrix), $I$(identity matrix) and $a$ which
are of $N\times N$ order respectively.
\begin{equation}
A=\begin{bmatrix}
0 & I\\
a & 0
\end{bmatrix},
\end{equation}
$b(\xi)$ is $2N\times1$ order column matrix and zero vector in this
case because of boundary value problems i-e $u_0=u_{N+1}=0$.
$Y$ can be found by integrating using exponential matrix \cite{a2,s6,t1,v1}.
\begin{equation}
Y=\exp(A\xi)Y_0+\int^{\xi}_{0}\exp[A(\xi-\lambda)]b(\lambda)d\lambda.
\end{equation}
Where $\lambda$ is a dummy variable of integration.
$Y_0$is a vector of initial conditions.
Like in our case $Y_0$ is chosen as
\begin{equation}
\begin{aligned}
Y&=[u_1,u_2,\dots,u_N,u_{N+2},u_{N+3},\dots,u_{2N+1}]^{T}_{\xi=0}\\
&=[0,0,\dots,0,c_1,c_2,\dots,c_N]^T.
\end{aligned}
\end{equation}
Now, let us solve the example for the simplest case $N=1$,
this makes the dependent variable because of the boundary condition $x=0$ and
$1$, as $u_0=0$ and $u_2=0$. Considering that, for the chosen example
matrix $b(\xi)$ is zero and so, one gets the solution on interior node
point because $N=1$.
\begin{equation}
Y=\begin{bmatrix}
u_1\\
u_3
\end{bmatrix}=\begin{bmatrix}
u_1\\
\frac{du_1}{d\xi}
\end{bmatrix}=\begin{bmatrix}
\cosh(\sqrt{2}\xi) & \frac{1}{\sqrt{2}}\sinh(\sqrt{2}\xi)\\
\sqrt{2} \sinh(\sqrt{2}\xi) & \cosh(\sqrt{2}\xi)
\end{bmatrix}\begin{bmatrix}
0\\
c_1
\end{bmatrix},
\end{equation}
on more simplification gives
\begin{equation}
Y=\begin{bmatrix}
\frac{c_1}{\sqrt{2}}\sinh(\sqrt{2}\xi)\\
c_1 \sinh(\sqrt{2}\xi)
\end{bmatrix}.
\end{equation}
which is considered as analytical solution because the above solution
is analytical in $y$. Where $c_1$ is still not found yet.
It can be calculated by using the boundary condition at $y=1$,
where $\epsilon$ is considered as $1$.
\begin{equation}
c_1=\sqrt{2} \frac{\sinh(\pi)}{\sinh(\sqrt{2})}.
\end{equation}
As we are using $N=1$ so there is only one constant exists like $c_1$
but if $N$ is more then $1$ then we would get more than just one constants.
\section{Exponential matrix calculation}
In our problem, Maple coding is used for obtaining the exponential matrix.
By increasing $N$, calculation time of matrix also increases and after
some certain value it would take hours to calculate the matrix.
Finding Eigenvalues and Eigenvector this particular issue can be resolved \cite{v1}.
Changing the conditions on $y=0$ or $1$ can not effect the coefficient matrix,
but changing in $x$ boundary condition would change the coefficient matrix.
So one has to be very careful while dealing with the different problems.
For the calculation at large scale, one can convert the coefficient matrix
into the canonical form as \cite{a2,v1}.
\begin{equation}
A=PBP^{-1} ,
\end{equation}
For elliptic PDEs it can be considered that all the eigenvalues are
real and distinct, where $B$is the diagonal matrix of $2N\times2N$ order and
$P$ is the eigen matrix and is given by
\begin{equation}
P=[P_1,P_2,\dots,P_{2N}].
\end{equation}
Now taking exponent of equation, it results in
\begin{equation}
\exp(A\xi)=P \exp(B\xi)P^{-1},
\end{equation}
where
\begin{equation}
\exp(B\xi)=\begin{bmatrix}
e^(\lambda_1)^\xi & 0 & \cdots &0 &0\\
0 & e^(\lambda_2)^\xi& \cdots &0 &0\\
\vdots&\vdots&\vdots&\vdots&\vdots\\
0 & 0 &\cdots & e^(\lambda_{2N-1})^\xi &0\\
0 & 0 &\cdots & 0 & e^(\lambda_{2N})^\xi
\end{bmatrix} .
\end{equation}
Here comes the discussion that Maple takes very few time as in seconds
to calculate the eigenvalues, but still there is a long required time
for the calculation of eigenvectors. Hence lets find some particular
eigenvector $P_k$ instead of the whole vector as \cite{v1}.
\begin{equation}
(A-\lambda_k U)P_k=0 ,
\end{equation}
taking $U$ is an identity matrix and
$P_k=[\beta_1,\beta_2,\cdots,\beta_{2N}]^T$ and putting these values in above
equation by letting $N\geq 3$ gives
\begin{equation}
\begin{bmatrix}
-\lambda_k\beta_1+\beta_{N+1}\\
\dots\\
-\lambda_k\beta_N+\beta_{2N}\\
-\lambda_k\beta_{N+1}-\beta_2+2\beta_1\\
-\lambda_k\beta_{N+2}+2\beta_2-\beta_1-\beta_3\\
\dots\\
-\lambda_k\beta_{2N-1}+2\beta_{N-1}-\beta_{N-2}-\beta_N\\
-\lambda_k\beta_{2N}+2\beta_N-\beta_{N-1}
\end{bmatrix}=0 ,
\end{equation}
so the above equation can be solved by letting initially $\beta_1=1$,
which start to give the values of remaining expressions as.
\begin{equation}
\beta_2=(2-\lambda^2_k)\beta_1,
\beta_i=\beta_{i-2}+(2-\lambda^2_k)\beta_{i-1},
\beta_{N+i}=\lambda_k\beta_i .
\end{equation}
Here one only needs to find the exponential matrix only for once and it
would become automatically valid for any boundary condition that is at
$y=0$ or $1$ and at $x=1$ but if at $x=0$ conditions changes it would
change $\beta_2$ value.
Dimensionless temperature profile for the solved example is given in
figures with varying number of $N$. As $N$ increases, results become more
accurate but time taken for the computation increases drastically.
Few of the profiles obtained for different N are as under.
\begin{figure}[ht]
\begin{center}
\includegraphics[width=8cm,height=5cm]{fig1.jpg}
\end{center}
\caption{2D dimensionless temperature profile for linear elliptic
partial differential equation taking $N=2$} \label{fig1}
\end{figure}
\begin{figure}[ht]
\begin{center}
\includegraphics[width=8cm,height=5cm]{fig2.jpg}
\end{center}
\caption{3D dimensionless temperature profile for linear elliptic partial
differential equation taking $N=2$} \label{fig2}
\end{figure}
\begin{figure}[ht]
\begin{center}
\includegraphics[width=8cm,height=5cm]{fig3.jpg}
\end{center}
\caption{2D dimensionless temperature profile for linear elliptic partial
differential equation taking $N=3$} \label{fig3}
\end{figure}
\begin{figure}[ht]
\begin{center}
\includegraphics[width=8cm,height=5cm]{fig4.jpg}
\end{center}
\caption{3D dimensionless temperature profile for linear elliptic partial
differential equation taking $N=3$} \label{fig4}
\end{figure}
\begin{figure}[ht]
\begin{center}
\includegraphics[width=8cm,height=5cm]{fig5.jpg}
\end{center}
\caption{2D dimensionless temperature profile for linear elliptic
partial differential equation taking $N=8$} \label{fig5}
\end{figure}
\begin{figure}[ht]
\begin{center}
\includegraphics[width=8cm,height=5cm]{fig6.jpg}
\end{center}
\caption{3D dimensionless temperature profile for linear elliptic
partial differential equation taking $N=8$} \label{fig6}
\end{figure}
\subsection{Error estimates}
For error analysis, we work out the average flux at $x=0$ (along $y$)
using $\epsilon =1,$ in equation \eqref{e2.3}, to get \cite{s7}.
\begin{equation}
u=\cosh(\pi)-1.
\end{equation}
The associated error with our proposed method is obtained by using the
following formula
\begin{equation} \label{e3.9}
{\rm Error} (\%) = \frac{|{\rm flux}_{\rm semi-analytical}
-{\rm flux}_{\rm exact}|}{{\rm flux}_{\rm exact}}.
\end{equation}
The error obtained from equation \eqref{e3.9}
is decreasing as we increase the number of interior collocations points $N$
as shown in Figure \ref{fig7}. For $N=10$, the error of our proposed method
is less then 0.5\% for average flux as compare to the other methods,
which is clear evidence that our proposed method is most suitable of
these type of differential equations.
\begin{figure}[ht]
\begin{center}
\includegraphics[width=8cm,height=5cm]{fig7.jpg}
\end{center}
\caption{Error(\%) using false transient method (FTM), semi-analytical
Finite Difference Method (FDM) and semi-analytical spectral methods
(SPM) for different interior collocation points $N$} \label{fig7}
\end{figure}
\subsection{Nonlinear boundary conditions}
Let there be a boundary value problem having the nonlinear boundary condition as
\begin{equation}
\frac{\partial^2 u}{\partial x^2}+\frac{\partial^2 u}{\partial y^2}=0,
\end{equation}
with boundary conditions
\begin{equation}
\begin{gathered}
u(0,y)=0 ,\quad\text{for }0\leq y\leq 1,\\
\frac{\partial u}{\partial x}(1,y)=0 \quad\text{for }0\leq y\leq 1,\\
\frac{\partial u}{\partial y}(x,0)=u^4(x,0) \quad\text{for }0\leq x\leq 1,\\
u(x,1)=1\quad\text{for }0\leq x\leq 1.
\end{gathered}
\end{equation}
Exponential matrix for semi-analytical method do not get any change in
this case. The only change we are facing while solving the nonlinear
boundary condition is in calculating the constants i-e $c_i$'s.
The profile is obtained at $N=5$ and is shown in figure \ref{fig8}.
\begin{figure}[ht]
\begin{center}
\includegraphics[width=8cm,height=5cm]{fig8.jpg}
\end{center}
\caption{Dimensionless temperature distribution inside a rectangle for
Elliptic partial differential equation having non linear boundary condition}
\label{fig8}
\end{figure}
\subsection{Nonlinear elliptic partial differential equation}
The method discussed in this paper can also be useful for solving non-linear
elliptic PDEs. to illustrate this, let us have equations as follows:
\begin{equation}
\frac{\partial^2 u}{\partial x^2}+\frac{\partial^2 u}{\partial y^2}=u^2,
\end{equation}
with boundary conditions
\begin{equation}
\begin{gathered}
\frac{\partial u}{\partial x}(0,y)=0 \quad\text{for }0\leq y\leq 1,\\
u(1,y)=0 \quad\text{for }0\leq y\leq 1,\\
u(x,0)=0 \quad\text{for }0\leq x\leq 1,\\
\frac{\partial u}{\partial y}(x,1)=0\quad\text{for }0\leq x\leq 1.
\end{gathered}
\end{equation}
As one can see that because of term $u^2$ the presenting problem is acting
as nonlinear, so the only change we have in our program is by adding
the $u^2(x,y)$ as
\begin{equation}
\frac{du_{N+1+i}}{d \xi}
=\Big[\sum^{N}_{l=0} \left(D^2_N\right)_{jl} \ u_l(y)\Big]+u^2_i(\xi),
\end{equation}
that would cast change in matrix $b$ not $A$. Instead of zero matrix $b$
would be now
\begin{equation}
b=[0,0,\cdots,0,u^2_1(\xi),u^2_2(\xi),\cdots,u^2_N(\xi)]^T
\end{equation}
Further, after implying the required changes in maple coding, At $N=5$
one would gives the result as figure \ref{fig42}
\begin{figure}[ht]
\begin{center}
\includegraphics[width=8cm,height=5cm]{fig9.jpg}
\end{center}
\caption{Dimensionless temperature distribution inside a rectangle
for non linear Elliptic partial differential equation} \label{fig42}
\end{figure}
\subsection*{Conclusion}
A semi-analytical solution based on matrix exponential method is used for the
solution of elliptic PDEs. The method is then extended to non-linear elliptic
PDEs and elliptic PDEs with non-linear boundary conditions. A Chebyshev-spectral
method is used to discretize the given elliptic PDEs in $x$-direction combining
with any other analytical technique in $y$-direction. The numerical simulations
results was obtained for different numbers of $N$ to observe the accuracy of our
our proposed semi-analytical method. It was found that our numerical simulations
results have a very good agreement with the other available semi-analytical
techniques. Even though the technique has been developed for a single elliptic
linear and non-linear and with non-linear boundary conditions, the same concept
could be extended to couple linear and non-linear elliptic PDEs.
\begin{thebibliography}{99}
\bibitem{a1} Amen, S.; Bilokon, P.; Codd, A. B.; Fofaria, M.; Shah, T.;
\emph{Numerical solutions of Differential Equations}, publisher?? 2004.
\bibitem{a2} Amundson, N.R.;
\emph{Mathematical Methods in Chemical Engineering. In: Matrices and Their
Application}, Vol. 1. Englewood Cliffs, New Jersey, 1966.
\bibitem{b1} Bokil, Vrushali A.; Gibson, Nathan L.;
\emph{Finite Difference, Finite Element and Finite Volume Methods for
the Numerical Solution of PDEs}. Oregon State University Corvallis, 2007.
\bibitem{c1} Carnahan, B.; Luther, H. A.; Wilkes, J. O.;
\emph{Applied Numerical Methods}. John Wiley and sons, New York, 1969.
\bibitem{c2} Carslaw, H.S.; Jaeger, J. C.;
\emph{Conduction of Heat in Solids}. Oxford University Press, London, 1973.
\bibitem{c3} Constantinides, A.; Mostou, N.;
\emph{Numerical Methods for Chemical Engineers with Matlab Applications}.
Prentice-Hall, New Jersey, 1999.
\bibitem{c4} Costa, Bruno;
\emph{Spectral methods for partial differential equations}. Cubo-Revista de
Matematica,2004. pp. 1-32.
\bibitem{c5} Cutlip, M. B.; Shacham, M.;
\emph{The numerical method of lines for partial differential equations}.
CACHE News 47, 1998. pp. 18-21.
\bibitem{d1} Davis, M. E.;
\emph{Numerical Methods and Modeling for Chemical Engineers}.
John Wiley and sons, New York, 1984.
\bibitem{l1} von Luxburg, Ulrike;
\emph{A tutorial on spectral clustering}, Statistics and computing, 2007.
pp. 395-416.
\bibitem{m1} Moler, C.; van Loan, C.;
\emph{Nineteen dubious ways to compute the exponential of a matrix}.
SIAM Review 20 (5), 1978. pp. 801-836.
\bibitem{p1} Pfeiffer, H. P.; Kidder, L. E.; Scheel, M. A.; Teukolsky, S. A.;
\emph{ A multidomain spectral method for solving elliptic equations}.
Computer physics communications, 152(3) (2003), 253-273.
\bibitem{r1} Rice, R. G.; Do, D. D.;
\emph{Applied mathematics and modelling for chemical
Engineers}. John Wiley \& Sons Inc., New York, 1995.
\bibitem{s1} Schiesser, W. E.;
\emph{The Numerical Method of Lines}. Academic Press, New York, 1991.
\bibitem{s2} Schiesser, W. E.;
\emph{Computational mathematics in Engineering and Applied Science: ODEs,
DAEs and PDEs}. CRC Press, 1994.
\bibitem{s3} Schiesser, W. E.; Silebi, C. A.;
\emph{Computational Transport Phenomena}.
Cambridge University Press, New York, 1997.
\bibitem{s4} J. Shen, T. Tang;
\emph{Spectral and High-Order Methods with Applications}, Science Press, 2006.
\bibitem{s5} J. Shen, T. Tang; L.-L. Wang;
\emph{Spectral Methods, Algorithms, Analysis and Applications},
Springer Series in Computational Mathematics, vol. 41, Springer,
Berlin, Heidelberg and New York, 2011.
\bibitem{s6} Subramanian, V. R.; White, R. E.;
\emph{Solving differential equations with maple}.
Chemical Engineering Education 34, 2000a. pp. 328-336
\bibitem{s7} Subramanian, V. R.; White, R. E.;
\emph{Semi-analytical method of lines for solving elliptic partial
differential equations}, Chem. Eng. Sci., 59 (2004), pp. 781-88
\bibitem{s8} Subramanian, V. R.; White, R. E.;
\emph{A semianalytical method for predicting primary and secondary current
density distributions: linear and non-linear boundary conditions}.
Journal of Electrochemical Society 147 (5), 2000. pp. 1636-1644.
\bibitem{t1} Taylor, R.;
\emph{Engineering computing with maple: solution of PDEs via method of lines}.
CACHE News 49, 1999. pp. 5-8.
\bibitem{t2} Taylor, R.; Krishna, R.;
\emph{Multi Component Mass Transfer}. Wiley, New York, 1993.
\bibitem{v1} Varma, A.; Morbidelli, M.;
\emph{Mathematical Methods in Chemical Engineering}. Oxford University Press,
New York, 1997.
\bibitem{w1} Weisstein, E. W.;
\emph{CRC concise encyclopedia of mathematics}. CRC press, 2002. pp. 1451-1553.
\end{thebibliography}
\end{document}