\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}