% !TEX encoding = IsoLatin %% Antes de processar este arquivo LaTeX (LaTeX2e) deve %% verificar que o arquivo TEMA.cls esta na mesma %% pasta. O arquivo TEMA.cls pode ser obtido do %% endereco http://tema.sbmac.org.br. \documentclass{TEMA} %\usepackage[brazil]{babel} % para texto em Português \usepackage[english]{babel} % para texto em Inglês \usepackage[latin1]{inputenc} % para acentuação em Português IsoLatin1 %\usepackage[utf8]{inputenc} % para acentuação em Português UTF-8 \usepackage{graphicx} % insert here the call for the packages your document requires \usepackage{amsmath,amsfonts,amssymb} \usepackage{latexsym} \usepackage[T1]{fontenc} % para PDF \usepackage{hyperref} \usepackage{color} \usepackage{ae} \usepackage{longtable} \usepackage{dsfont} % para exibir melhor o conjuntos dos numeros \usepackage{mathrsfs} %%%%%%% Tipo de Letra \usepackage{epstopdf} %%%%%% transforma eps em pdf %%%%%%% \usepackage{hyperref} \hypersetup{ colorlinks=true, linkcolor=blue, %trocar para black depois (Junior) citecolor=black, filecolor=black, urlcolor=blue , } %%%%%% \newcommand{\prova}{\noindent\textbf{Proof. }} %começo de uma prova \newcommand{\fimprova}{$\hfill\square$} %fim de uma prova \newtheorem{theorem}{Theorem}[section] \newtheorem{remark}{Remark}[section] \newtheorem{lemma}{Lemma}[section] % Atenção: O presente arquivo usa a codificação IsoLatin1. Para % utilizar a codificação UTF-8, converter o arquivo usando seu % programa favorito, então troque o argumento do inputenc % adequadamente. %\usepackage[dvips]{graphics} %\usepackage{subfigure} \usepackage{graphicx} \usepackage{epsfig} \usepackage{hyperref} \usepackage{framed} %\usepackage{psfrag} %\usepackage{tikz} \newcommand{\B}{{\tt\symbol{92}}} \newcommand{\til}{{\tt\symbol{126}}} \newcommand{\chap}{{\tt\symbol{94}}} \newcommand{\agud}{{\tt\symbol{13}}} \newcommand{\crav}{{\tt\symbol{18}}} \begin{document} %******************************************************** \title{ Numerical solution of heat equation with singular \\ Robin boundary condition% \thanks{Paper presented at the XXXVI National Congress of Applied and Computational Mathematics, Gramado-RS, Brazil, 2016.} } \author{ G. LOZADA-CRUZ% \thanks{ Partially supported by the São Paulo Research Foundation (FAPESP, grant: 2015/24095-6), E-mail: german.lozada@sjrp.unesp.br}, São Paulo State University (Unesp), Institute of Biosciences, Humanities and Exact Sciences, Departament of Mathematics, 15054-000, São José do Rio Preto, SP, Brazil. \\ \\ C. E. RUBIO-MERCEDES% \thanks{Supported by The Mato Grosso do Sul Research Foundation (FUNDECT, grant: 219/2016), E-mail: cosme@uems.br}, Mathematics and Engineering Physics Programs, UEMS - Universidade Estadual do Mato Grosso do Sul, Dourados, MS, Brazil. \\ \\ J. RODRIGUES-RIBEIRO% \thanks{Partially supported by Coordinator for the Improvement of Higher Education Personnel (CAPES, grant: PICME-9725290/M), E-mail: jrodrib@usp.br} Institute of Mathematics Sciences and Computing, USP - Universidade de São Paulo, São Carlos, SP, Brazil } \criartitulo \runningheads{G.Lozada-Cruz, C.Rubio-Mercedes, J.Rodrigues-Ribeiro}% {Numerical solution of heat equation} \begin{abstract} {\bf Abstract}. In this work we study the numerical solution of one-dimensional heat diffusion equation with a small parameter $\epsilon>0$ subject to Robin boundary conditions. The simulations examples lead us to conclude that the numerical solutions of the differential equation with Robin boundary condition are very close of the analytic solution of the problem with homogeneous Dirichlet boundary conditions when $\epsilon $ tends to zero. {\bf Keywords}. Eigenvalue Problems, \and Finite Difference Method, \and Robin Boundary Conditions, \and Numerical Solutions. \end{abstract} %******************************************************** \section{Introduction} \label{intro} It is well known that the diffusion differential equation models the transient conduction phenomenon that occurs in numerous engineering applications and may be analyzed by using different analytic and numerical methods. Many transient problems involving geometry and simple boundary value conditions, their analytic solution are known explicitly, especially the one-dimensional (1D) case. Still for the two-dimensional (2D) and three-dimensional (3D) cases some of the analytic solutions are known (see \cite{JMAANCARB,ASMAR,INCROPERA,DIPRIMA}). However, in most cases, the geometry or boundary conditions make it impossible to apply analytic techniques to solve the heat diffusion equation (see \cite{ASMAR,DIPRIMA}). In this work we use the Crank-Nicolson Finite Difference Method (FDM) (see \cite{CUMINATO,HOFFMAN}) to solve the 1D heat diffusion equation in transient regime with Robin boundary conditions given by \begin{equation}\label{eq:1} \begin{cases} u_t =u_{xx}, \; -10,\\ -\epsilon u_x(-1) = u(-1), &\\ \epsilon u_x(1) =u(1), & \end{cases} \end{equation} where $\epsilon\in (0,1]$. If $\epsilon=0$ in \eqref{eq:1} we have the classical problem with homogeneous Dirichlet boundary conditions for the heat equation which is well known. The problem \eqref{eq:1} with singular boundary conditions at the boundary of the domain, depending on a positive parameter, has not been studied previously neither analytically nor numerically. Our little contribution with this kind of problems which depend of a small parameter is to show numerical solutions when we vary the values of $\epsilon$. This paper is organized as follows. In Section \ref{sec:1}, we state that the equation \eqref{eq:1} has unique solution for each $\epsilon>0$. Also we study asymptotic behaviour of the eigenvalues of \eqref{eq:10} when $\epsilon $ tends to zero. In Section \ref{sec:4} a brief description of the problem with Robin boundary condition in conjunction with the FDM approach is presented. The numerical experiments are discussed in Section \ref{sec:5}, and finally the conclusions are presented in Section \ref{sec:6}. \section{On the existence of solution}\label{sec:1} Let $\Omega=(-1, 1)$ and the Lebesgue space $X=L^2(\Omega)$. Let $A_\epsilon: D(A_\epsilon)\subset L^2(\Omega)\to L^2(\Omega)$ be an unbounded linear operator defined through \begin{align}\label{eq:dea} D(A_\epsilon) &=\big\{ u\in H^2(\Omega): -\epsilon u_x(-1) = u(-1), \, \epsilon u_x(1) =u(1)\big\},\\ A_\epsilon u &= u''. \end{align} Thus we can write the equation (\ref{eq:1}) as an evolution equation in $L^2(\Omega)$ (see \cite{HENRY}) as follows \begin{equation}\label{eq:abstract} \begin{cases} \dot{u} =A_\epsilon u, \; t>0,\\ u(0)= u_0\in H^1(\Omega). \end{cases} \end{equation} \begin{theorem} For each $u_0\in H^1(\Omega)$, there exists a unique solution $u=u(t, u_0)$ of (\ref{eq:abstract}) defined on its maximal interval of existence $[0, \tau_{u_0})$ which mean that either $\tau_{u_0}=+\infty$, or if $\tau_{u_0}<\infty$ them $\limsup\limits_{t\to \tau_{u_0}^{-}}\|u(t, u_0)\|=+\infty$. \end{theorem} \prova See \cite{JMAANCARB,HENRY}. \fimprova \medskip It is well known that for a fixed value $\epsilon>0$, the problem (\ref{eq:abstract}) generates a well-defined linear semigroup in $H^1(\Omega)$, the solutions enter $W^{1,p}(\Omega)$ for any $p$ such that $10, \; -10$ the unique equilibrium solution of (\ref{eq:2}) is $u_\epsilon\equiv 0$. \end{theorem} \prova The solution of the problem \eqref{eq:2} is given by \begin{equation}\label{eq:sol-eigenvalue} u(x)=\Big(\frac{u(1)-u(-1)}{2}\Big) x +\frac{u(1)+u(-1)}{2}, \quad x\in \Omega. \end{equation} By the boundary conditions (\ref{eq:2}) in (\ref{eq:sol-eigenvalue}) we have that $u(-1)$ and $u(1)$ satisfy \begin{equation}\label{eq:3} u(1) = \epsilon \frac{u(1)-u(-1)}{2}= -u(-1). \end{equation} Thus we have $u(-1)+ u(1)=0$, and (\ref{eq:sol-eigenvalue}) provides \begin{equation}\label{eq:6} u(x)= \Big(\frac{u(1)-u(-1)}{2}\Big) x. \end{equation} Again, by using the boundary conditions (\ref{eq:2}) we have \begin{equation}\label{eq:7} \frac{u(1)-u(-1)}{2} = \epsilon \frac{u(1)-u(-1)}{2}, \end{equation} and for $\epsilon \neq 1$ we have \begin{equation}\label{eq:9} u(1)-u(-1)=0. \end{equation} Finally, from \eqref{eq:9} we conclude that $u_\epsilon\equiv 0$. \fimprova \begin{remark} From the Theorem \ref{teo:2} we can say that $u_\epsilon\equiv 0$ converges to the solution $u\equiv 0$ of the problem with homogeneous Dirichlet bounday conditions. \end{remark} \subsection{Eigenvalue problem} \label{sec:3} Consider the eigenvalue problem associated with the linear operator $A_\epsilon$ \begin{equation}\label{eq:10} \begin{cases} u_{xx} = \lambda^\epsilon\, u,\;\; \hbox{in}\; \Omega,\\ -\epsilon u_x(-1) = u(-1),\\ \epsilon u_x(1) = u(1). \end{cases} \tag{$\mathcal{E}_\epsilon$} \end{equation} For each $\epsilon>0$, the problem (\ref{eq:10}) has a sequence of real eigenvalues $\{\lambda_n^\epsilon\}_{n=1}^\infty$, with an $L^2(\Omega)$-orthogonal and complete associated system of eigenfunctions $\{\varphi_n^\epsilon\}_{n=1}^\infty$. This conclusion is possible because the operator $A_\epsilon$ is densely defined closed and self-adjoint in $L^2(\Omega)$. From (\ref{eq:2}) we conclude that zero is not an eigenvalue for (\ref{eq:10}). The variational formulation for $\lambda_n^\epsilon$ is given by \begin{equation}\label{eq:lamdan} \lambda_n^\epsilon = \underset{\psi\in Y_n}{\min}\left\{ \frac{\epsilon (\psi'(1)^2 + \psi'(-1)^2) -\int_\Omega |\psi'| ^2}{\int\limits_\Omega \psi^2}\right\}, \end{equation} where $$Y_n=\Big\{w\in C^2(\Omega): w\neq 0, w\big|_{-1}^1= \pm \epsilon w\big|_{-1}^1, \int_\Omega w \varphi_j^\epsilon=0, \forall j=1, 2, \cdots, n-1\Big\}.$$ For more details about the variational formulation of boundary value problems see for example \cite[Chapter 8]{Brezis} and \cite[Chapter 5]{Butazzo}. \bigskip Let $\lambda^\epsilon= \omega^2$, $\omega \in \mathds{R}$, and $\varphi^\epsilon(x)=\cosh(\omega x)$ the eigenfunction associated with $\lambda^\epsilon$. From the boundary conditions for $\varphi^\epsilon$, we have \begin{equation}\label{eq:12} %%%%% \tanh \omega =\frac{1}{\epsilon \omega}. \end{equation} %%%%%%%%%% \begin{figure}[h!] \centering \includegraphics[width=0.8\textwidth]{Fig1.pdf} \caption[]{Graphical solution of $\tanh (\omega) = \frac{1}{\epsilon \omega}$ for $\epsilon=0.1, 0.09, 0.08, 0.07$.} \label{Fig1} \end{figure} The solutions of (\ref{eq:12}) can be determined numerically. They can also be obtained approximately by sketching the graphs of $\psi_1(\omega)=\tanh \omega$ and $\psi_2(\omega)=\frac{1}{\epsilon \omega}$ for $\epsilon=0.1, 0.09, 0.08, 0.07$, and identifying the points of intersection of the curves (see in Fig. \ref{Fig1}). Let $\omega_1(\epsilon)$ be the interception points of the curves $\psi_1(\omega)$ and $\psi_2(\omega)$. Thus $\lambda_1^\epsilon=\omega_1(\epsilon)^2$. Now, taking the eigenfunction $\varphi^\epsilon(x)=\sinh(\omega x)$ and using the boundary condition, we have \begin{equation}\label{eq:13} \tanh \omega =\epsilon \omega. \end{equation} In Fig. \ref{Fig2} we also have plotted the graphs of $\phi_1(\omega)=\tanh \omega$ and $\phi_2(\omega)=\epsilon \omega $ for $\epsilon=0.1, 0.09, 0.08, 0.07$, as function of $\omega$. Let $\omega_2(\epsilon)$ be the interception points of the curves $\phi_1(\omega)$ and $\phi_2(\omega)$. Thus $\lambda_2^\epsilon=\omega_2(\epsilon)^2$. %%%%%%%%%%% \begin{figure}[h!] \centering \includegraphics[width=0.8\textwidth]{Fig2.pdf} \caption[]{Graphical solution of $\tanh (\omega) =\epsilon \omega$ for $\epsilon=0.1, 0.09, 0.08, 0.07$.} \label{Fig2} \end{figure} \smallskip From Figs. \ref{Fig1} and \ref{Fig2} we can observe that the eigenvalues $\lambda_1^\epsilon$ and $\lambda_2^\epsilon$ increase continually when $\epsilon$ tends to zero, respectively. \begin{lemma} Let $\epsilon>0$. Then $\lambda_1^\epsilon >\lambda_2^\epsilon$. Moreover $\lambda_1^\epsilon,\lambda_2^\epsilon\to \infty$ when $\epsilon\to 0$. \end{lemma} \prova. We know that the functions $\tanh \omega$ and $ \coth \omega$ can be write as \begin{eqnarray} \tanh \omega &=& 1-2\, e^{-2 \omega} + 2 \, e^{-4 \omega} -2\, e^{-6\omega} + \cdots \\ \coth \omega &=& 1+2\, e^{-2 \omega} + 2 \, e^{-4 \omega} +2\, e^{-6\omega} + \cdots \end{eqnarray} From the equation (\ref{eq:12}), we have \begin{equation}\label{eq:14} \omega = \frac{1}{\epsilon} \frac{1}{\tanh \omega}. \end{equation} When $\omega$ is large enough, $\tanh \omega$ approaches $1$. Thus, from (\ref{eq:14}) we have \begin{equation}\label{eq:lamda1} \omega = \frac{1}{\epsilon} \coth \frac{1}{\epsilon} = \frac{1}{\epsilon} \left(1+ 2\, e^{-2/\epsilon} + 2 \, e^{-4/\epsilon} +2\, e^{-6/\epsilon} + \cdots \right). \end{equation} Since $\lambda_1^\epsilon= \omega_1(\epsilon)^2$, we have \begin{equation}\label{eq:15} \lambda_1^\varepsilon= \frac{1}{\epsilon^2} + \frac{4}{\epsilon^2} e^{-2/\epsilon} + \frac{8}{\epsilon^2} e^{-4/\epsilon} + \cdots. \end{equation} From (\ref{eq:13}) we obtain \begin{equation}\label{eq:16} \omega = \frac{1}{\epsilon} \tanh \omega. \end{equation} When $\omega$ is large enough, $\tanh \omega$ approaches $1$. Thus, from (\ref{eq:16}) we have \begin{align} \omega = \frac{1}{\epsilon} \tanh \frac{1}{\epsilon} = \frac{1}{\epsilon} \left(1 - 2\, e^{-2/\epsilon} + 2 \, e^{-4/\epsilon} -2\, e^{-6/\epsilon} + \cdots \right). \end{align} Since $\lambda_2^\epsilon= \omega_2(\epsilon)^2$, we have \begin{equation}\label{eq:17} \lambda_2^\varepsilon= \frac{1}{\epsilon^2} - \frac{4}{\epsilon^2} e^{-2/\epsilon} + \frac{8}{\epsilon^2} e^{-4/\epsilon} - \cdots. \end{equation} From (\ref{eq:15}) and (\ref{eq:17}) follow the results. \fimprova Also from (\ref{eq:15}) and (\ref{eq:17}) we have that the gap between $\lambda_1^\epsilon$ and $\lambda_2^\epsilon$ is given by \begin{equation}\label{eq:18} \lambda_1^\epsilon- \lambda_2^\epsilon \approx \frac{8}{\epsilon^2} e^{-2/\epsilon}. \end{equation} \begin{remark} When $\epsilon$ tends to zero, the problem (\ref{eq:10}) becomes \begin{equation}\label{eq:11} \begin{cases} u_{xx} = \lambda^0 \, u,\quad -10$ and the temporal nodal points are indexed by $t_j=jk$. With this indexes, we can use the following notation for the values of $u$: $u_{ij}=u(x_i,t_j)$ and $u_{i}=u(x_i,t)$. From the boundary condition given in \eqref{eq:1} at $x=-1$ and for $\epsilon >0$ we write \[u_x(-1,t)=\frac{-u(-1,t)}{\epsilon}.\] Assuming that at $x=-1$, the function $u$ is twice differentiable in $x$, so that we can write \[u_{xx}(-1,t)=\frac{-u_x(-1,t)}{\epsilon}=\frac{u(-1,t)}{\epsilon^2},\] and in the same way at $x=1$ we get \[u_{xx}(1,t)=\frac{u_x(1,t)}{\epsilon}=\frac{u(1,t)}{\epsilon^2}.\] Also, by using the finite difference approach for $u_{xx}(t)$ (see \cite{CUMINATO}), the problema \eqref{eq:1} can be write \begin{equation} \frac{d u_i(t)}{d t} = \frac{u_{i+1}(t) -2u_{i}(t)+u_{i-1}(t)}{h^2} \label{CEN} \end{equation} for $i=1,...,n-1$, and for $i=0$ and $i=n$, we have \begin{equation} \frac{d u_0(t)}{d t} =\frac{u_0(t)}{\epsilon^2}, \label{CIleft} \end{equation} and \begin{equation} \frac{d u_n(t)}{d t}=\frac{u_n(t)}{\epsilon^2},\label{CIrigth} \end{equation} respectively. Thus by using \eqref{CEN}-\eqref{CIrigth}, the problem \eqref{eq:1} can be transformed into the first-order matrix differential equation given by \begin{equation}\label{matriz4nova} \dfrac{\partial}{\partial t} \begin{bmatrix} u_{0} \\ u_{1} \\ u_{2} \\ \vdots \\ u_{n-2} \\ u_{n-1} \\ u_{n} \end{bmatrix} =\dfrac{1}{h^2} \begin{bmatrix} \frac{h^2}{\epsilon^2} & 0 & 0 & \cdots & 0 & 0 & 0 \\ 1 & -2 & 1 & \ddots & & & 0 \\ 0 & 1 & -2 & \ddots & \ddots & & 0 \\ \vdots & \ddots & \ddots & \ddots & \ddots & \ddots & \vdots \\ 0 & &\ddots & \ddots & -2 & 1 & 0 \\ 0 & & & \ddots & 1 & -2 & 1 \\ 0 & 0 & 0 & \cdots & 0 & 0 & \frac{h^2}{\epsilon^2} \end{bmatrix} \begin{bmatrix} u_{0} \\ u_{1} \\ u_{2} \\ \vdots \\ u_{n-2} \\ u_{n-1} \\ u_{n} \end{bmatrix}, \end{equation} or in compact way \begin{equation}\label{matriz4} \frac{d U}{d t}=\mathbf{L}U, \end{equation} where, the matrix $\mathbf{L}$ represents the discrete counterpart of the corresponding differential operator given in \eqref{eq:1}, which is a tridiagonal band matrix and $U$ denotes a vector with the unknown values of $u$ defined over the nodes of the spatial mesh. \\ There are two general methods to solve \eqref{matriz4}: the explicit and implicit finite difference schemes. The type of implicit scheme adopted in this work was the Crank-Nicolson algorithm, see \cite{CUMINATO,COSME}. \section{Numerical Results}\label{sec:5} For the numerical examples we consider the initial condition given by \begin{equation} u(x,0)=\cos\left(\frac{\pi}{2} x \right),\label{CI} \end{equation} where for homogeneous Dirichlet boundary condition, which is obtained from \eqref{eq:1} with $\epsilon=0$, the analytical solution is (see \cite{DIPRIMA}) \begin{equation} u(x,t)=\exp\left(-\frac{\pi ^2}{4} t \right) \cos\left(\frac{\pi}{2} x \right)\label{solucao}. \end{equation} %%%%%%%%%%%%% \begin{figure}[h!] \centering \includegraphics[width=0.8\textwidth]{Fig3.pdf} \caption{Numerical solutions of (\ref{eq:1}) when $\epsilon\rightarrow 0$ at $t=1.5$.} \label{Fig3} \end{figure} %%%%%%%%% \begin{figure}[h!] \centering \includegraphics[width=0.8\textwidth]{Fig4.pdf} \caption{Numerical solutions of (\ref{eq:1}) when $\epsilon\rightarrow 0$ at $x=0.8$.} \label{Fig4} \end{figure} The examples were solved using $n=20$ divisions in $ x$ axis and $m=100$ divisions in temporal axis $t$. In Fig. \ref{Fig3} we display the spatial behavior of the solution for $t=1.5$ considering several values of $\epsilon$: $\epsilon =0.08$, $\epsilon =0.05$, $\epsilon =0.025$, $\epsilon =0.0075$. In Fig. \ref{Fig4} we show the temporal behavior of the solution for the same values of $\epsilon$ used in the Fig. \ref{eq:1} but at $x=0.8$. From Figs. \ref{Fig3} and \ref{Fig4} we can see that the numerical solutions of (\ref{eq:1}), which is a boundary value problem with Robin boundary conditions, converges to the exact solution of the problem with homogeneous Dirichlet boundary conditions, when $\epsilon$ tends to zero. The Figs. \ref{Fig5} and \ref{Fig6} display the temporal evolution of the solutions with the initial condition given by \eqref{CI}, the numerical solution for $\epsilon =0.05$ in Fig. \ref{Fig5} and analytical solution in \ref{Fig6}. It is worth noting that the numerical solution presented in Fig. \ref{Fig5} for $\epsilon =0.05$ agrees very well with analytic solution of the problem with homogeneous Dirichlet boundary conditions presented in Fig. \ref{Fig6}. %%%%%%%%%%%%%%%%%%%%%%% \begin{figure}[!htb] \centering \includegraphics[width=0.8\textwidth]{Fig5.pdf} \caption{Temporal evolution of the numerical solution obtained with initial condition given by \eqref{CI} and for $\epsilon =0.05$.} \label{Fig5} \end{figure} %%%%%%%%%%%%%%%%% \section{Conclusions}\label{sec:6} %%%%%%%%%%%%%%%%%%%%%%%% \begin{figure}[!h] \includegraphics[width=0.8\textwidth]{Fig6.pdf} \caption{Temporal evolution of the analytical solution obtained in \eqref{solucao}, with initial condition given by \eqref{CI}.} \label{Fig6} \end{figure} In the first part of the paper, we prove several results on the well-posedness of the system \eqref{eq:1} and the associated stationary problem. In the second part, the application of the Crank-Nicolson for the numerical solution of \eqref{eq:1}, involving the Robin boundary condition, has been presented here. The usefulness of this numerical technique has also been demonstrated by means of examples involving the solution of \eqref{eq:1} for several values of $\epsilon$. Analytical approaches and numerical simulations have clearly illustrated the asymptotic behaviour of the solution of \eqref{eq:1} when $\epsilon$ tends to zero. Continuous dependence of the solutions of \eqref{eq:1} from the $\epsilon$ parameter has also been demonstrated. %%%%%%%%%%%%%% %\bibliographystyle{ieeetr} %\bibliography{Modelo-TEMA} %%%%%%%%%%%%%%%%%%% \medskip \noindent {\bf Resumo}. Neste trabalho, obtemos soluções numéricas da equação diferencial de difusão do calor unidimensional com um parâmetro pequeno $ \epsilon $ nas condições de contorno de Robin. Exemplos numéricos demonstram que as soluções numéricas do problema com condição de contorno de Robin convergem para a soluções analíticas do problema de Dirichlet homogêneo quando $ \epsilon $ tende a zero. %%%%%%%%%%%%%%%%%%%%%% %\section*{References} \begin{thebibliography}{99} \bibitem{JMAANCARB} J. M. Arrieta, A. N. Carvalho, and A. Rodriguez-Bernal, A., {\em Attractors for Parabolic Problems with Nonlinear Boundary Condition. Uniform boundeds}, Comunications in Partial Differential Equations {\bf 25}, 1-2, pp. 1-37, 2000. \bibitem{ASMAR} N. H. Asmar, {\em Partial Differential Equations with Fourier Series and Boundary Value Problems.} 2nd. ed., Pearson Education, Inc., Upper Saddle River, NJ, 2005. \bibitem{INCROPERA} T. L. Bergman, A. S. Lavine, F. P. Incropera, and D. P. Dewitt, {\em Fundamentals of Heat and Mass Transfer.} 7th. ed., John Wiley and Sons, Inc., Danvers, MA, 2007. \bibitem{DIPRIMA} W. E. Boyce and R. C. DiPrima, R. C. {\em Elementary Differential Equations and Boundary Value Problems.} 7th. ed., John Wiley and Sons, Inc., New York, 2001. \bibitem{Butazzo} G. Buttazzo, M. Giaquinta, S. Hildebrandt, {\em One-dimensional variational problems. An introduction}. Oxford Lecture Series in Mathematics and its Applications, 15. The Clarendon Press, Oxford University Press, New York, 1998. \bibitem{Brezis} H. Brezis, \textit{Functional analysis, Sobolev spaces and partial differential equations}. Universitext. Springer, New York, 2011. \bibitem{CUMINATO} J. A. Cuminato and Jr. M. Meneguette, {\em Discretização de Equaçõess Diferenciais Parciais: Técnicas de Diferenças Finitas.} 1. ed. Rio de Janeiro: SBM, 2013. \bibitem{HENRY} D. Henry, {\em Geometry theory of semilinear parabolic equations}, Lecture Notes in Mathematics {\bf 840}, Springer-Verlag, Berlim, 1981. \bibitem{COSME} H. E. Hernández-Figueroa and C. E. Rubio-Mercedes,{\em Transparent Boundary for the Finite-Element Simulation of Temporal Soliton Propagation}, IEEE Transaction on Magnetics, Vol. 34, No. 5, 1998. \bibitem{HOFFMAN} J. D. Hoffman {\em Numerical Methods for Engineers and Scientists.} 2nd. ed., Marcel Dekker, Inc., New York, 2001. \end{thebibliography} \end{document}