Numerical Solution of Heat Equation with Singular Robin Boundary Condition^†

ABSTRACT

In this work we study the numerical solution of one-dimensional heat diffusion equation subject to Robin boundary conditions multiplied with a small parameter epsilon greater than zero. The numerical evidences tell us that the numerical solution of the differential equation with Robin boundary condition are very close in certain sense of the analytic solution of the problem with homogeneous Dirichlet boundary conditions when ε tends to zero.

Keywords:
Eigenvalue Problems; Finite Difference Method; Robin Boundary Conditions; Numerical Solutions

RESUMO

Neste trabalho, obtemos soluções numéricas da equação diferencial de difusão do calor unidimensional com um parâmetro pequeno ε 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 ε tende a zero.

Palavras-chave:
Problema de Autovalores; Método de Diferenças Finitas; Condições de Contorno de Robin; Soluções Numéricas

1 INTRODUCTION

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 ^{1), (2}). However, in most cases, the geometry or boundary conditions make it impossible to apply analytic techniques to solve the heat diffusion equation.

In this work we use the Crank-Nicolson Finite Difference Method (FDM) (see ⁹) to solve the 1D heat diffusion equation in transient regime with Robin boundary conditions given by

$\{\begin{array}{c}{u}_t=u_{xx},-1<x<\mathrm{1,}t>\mathrm{0,}\\ -\epsilon u_x\left(-1\right)=u\left(-1\right),\begin{array}{cc}& \end{array}\\ \epsilon u_x\left(1\right)=u\left(1\right),\begin{array}{ccc}& & \end{array}\end{array}$ (1.1)

where ε ∈ (0,1].

If ε = 0 in (1.1) we have the classical problem with homogeneous Dirichlet boundary conditions for the heat equation which is well known.

There is great interest on heat problems and much work was done considering different boundary conditions. Nevertheless, the particular (1.1) problem with singular boundary conditions, 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 ε.

Many problems in the industry are modeled by the heat equation subject to specific initial and boundary conditions, and sometimes it is not possible to get the analytic solution. Actually many researchers use different numerical techniques to understand the behavior of the solution (for more details see ^{(5, 6, 10)} and in references therein).

This paper is organized as follows. In Section 2, we state that the equation (1.1) has unique solution for each ε > 0. Also we study asymptotic behaviour of the eigenvalues of (ℰ _ε ) when ε tends to zero. In Section 3 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 4, and finally the conclusions are presented in Section 5.

2 ON THE EXISTENCE OF SOLUTION

Let Ω = (-1, 1) and the Lebesgue space X = L ² (Ω). Let A _ε: D(A _ε ) ⊂ L ²(Ω) → L ²(Ω) be an unbounded linear operator defined through

$D\left(A_{\epsilon }\right)=\left\{u\in H^2\left(\Omega \right):-\epsilon u_x\left(-1\right)=u\left(-1\right),\epsilon u_x\left(1\right)=u\left(1\right)\right\},$ (2.1)

$A_{\epsilon }u=u".$ (2.2)

Thus we can write the equation (1.1) as an evolution equation in L ²(Ω) (see ⁷) as follows

$\{\begin{array}{c}\dot{u}=A_{\epsilon }u,t>\mathrm{0,}\begin{array}{cc}& \end{array}\\ u\left(0\right)=u_0\in H^1\left(\Omega \right).\end{array}$ (2.3)

Theorem 1.For each u₀ ∈ H ¹(Ω) there exists a unique solution u = u(t;u ₀) of(2.3)defined on its maximal interval of existence$\left[\mathrm{0,}{\tau }_{u_0}\right)$which mean that either${\tau }_{u_0}=+\infty$, or if${\tau }_{u_0}<\infty$them$\underset{t\to {\tau }_{u_0}}{\lim \sup }\Vert u\left(t,u_0\right)\Vert =+\infty$.

Proof. See ^1),(7. □

It is well known that for a fixed value ε > 0, the problem (2.3) generates a well-defined linear semigroup in H ¹(Ω), the solutions enter W ^1.p (Ω) for any p such that 1< p < ∞ and are classical for t positive (for more details see ^1),(7).

2.1 Equilibrium solution for (1.1)

The equilibrium solution of (1.1) satisfy the elliptic boundary value problem

$\{\begin{array}{c}{u}_{xx}=\mathrm{0,}t>\mathrm{0,}-1<x<\mathrm{1,}\\ -\epsilon u_x\left(-1\right)=u\left(-1\right),\begin{array}{cc}& \end{array}\\ \epsilon u_x\left(1\right)=u\left(1\right).\begin{array}{ccc}& & \end{array}\end{array}$ (2.4)

Theorem 2.For every ε > 0 the unique equilibrium solution of(2.4)is u _ε ≡ 0.

Proof. The solution of the problem (2.4) is given by

$u\left(x\right)=\left(\frac{u\left(1\right)-u\left(-1\right)}{2}\right)x+\frac{u\left(1\right)+u\left(-1\right)}{2},\begin{array}{c}\end{array}x\in \Omega .$ (2.5)

By the boundary conditions (2.4) in (2.5) we have that u(?1) and u(1) satisfy

$u\left(1\right)=\epsilon \frac{u\left(1\right)-u\left(-1\right)}{2}=-u\left(-1\right).$ (2.6)

Thus we have u(?1) + u(1) = 0, and (2.5) provides

$u\left(x\right)=\left(\frac{u\left(1\right)-u\left(-1\right)}{2}\right)x.$ (2.7)

Again, by using the boundary conditions (2.4) we have

$\frac{u\left(1\right)-u\left(-1\right)}{2}=\epsilon \frac{u\left(1\right)-u\left(-1\right)}{2},$ (2.8)

and for ε ≠ 1 we have

$u\left(1\right)-u\left(-1\right)=0.$ (2.9)

Finally, from (2.9) we conclude that u _ε ≡ 0. □

Remark 3. From the Theorem 2 we can say that u_ε ≡ 0 converges to the solution u ≡ 0 of the problem with homogeneous Dirichlet bounday conditions.

2.2 Eigenvalue problem

Consider the eigenvalue problem associated with the linear operator A _ε

$\{\begin{array}{c}{u}_{xx}={\lambda }^{\epsilon }u,\begin{array}{c}\end{array}\text{in}\begin{array}{c}\end{array}\Omega ,\\ -\epsilon u_x\left(-1\right)=u\left(-1\right),\\ \epsilon u_x\left(1\right)=u\left(1\right).\begin{array}{cc}& \end{array}\end{array}$

(ℰ ₍ )

For each ε > 0, the problem (ℰ _ε ) has a sequence of real eigenvalues ${\left\{{\lambda }_n^{\epsilon }\right\}}_{n=1}^{\infty }$, with an L ²(Ω)-orthogonal and complete associated system of eigenfunctions ${\left\{{\phi }_n^{\epsilon }\right\}}_{n=1}^{\infty }$. This conclusion is possible because the operator A_ε is densely defined closed and self-adjoint in L ²(Ω).

From (2.4) we conclude that zero is not an eigenvalue for (ℰ _ε ). The variational formulation for ${\lambda }_n^{\epsilon }$ is given by

${\lambda }_n^{\epsilon }=\underset{\psi \in Y_n}{\min }\left\{\frac{\epsilon \left(\psi \text{'}{\left(1\right)}^2+\psi \text{'}{\left(-1\right)}^2\right)-{\int }_{\Omega }{\left|\psi \text{'}\right|}^2}{\underset{\Omega }{\int }{\psi }^2}\right\},$ (2.10)

where

$Y_n=\left\{w\in C^2\left(\Omega \right):w\ne \mathrm{0,}w{|}_{-1}^1=\pm \epsilon w{|}_{-1}^1,{\displaystyle {\int }_{\Omega }w{\phi }_j^{\epsilon }=\mathrm{0,}\forall j=\mathrm{1,2,...,}n-1}\right\}.$

For more details about the variational formulation of boundary value problems see for example [³, Chapter 8] and [⁴, Chapter 5].

Let λ^ε = ω², ω ∈ ℝ, and φ ^ε (x) = cosh(ωx) the eigenfunction associated with λ⁽. From the boundary conditions for φ ^ε , we have

$\tanh \omega =\frac{1}{\epsilon \omega }.$ (2.11)

The solutions of (2.11) can be determined numerically. They can also be obtained approximately by sketching the graphs of Ψ ₁(ω) = tanh ω and ${\psi }_2\left(\omega \right)=\frac{1}{\epsilon \omega }$ for ε = 0.1, 0.09, 0.08, 0.07, and identifying the points of intersection of the curves (see in Fig. 1). Let ω ₁(ε) be the interception points of the curves Ψ ₁(ω) and Ψ ₂(ω). Thus ${\lambda }_1^{\epsilon }={\omega }_1{\left(\epsilon \right)}^2$.

Figure 1:
Graphical solution of $\tanh \left(\omega \right)=\frac{1}{\epsilon \omega }$ for ε = 0:1; 0:09; 0:08; 0:07.

Now, taking the eigenfunction φ ^ε (x) = sinh(ωx) and using the boundary condition, we have

$\tanh \omega =\epsilon \omega .$ (2.12)

In Fig. 2 we also have plotted the graphs of ϕ ₁(ω) = tanh ω and ϕ₂(ω) = εω for ε = 0.1, 0.09, 0.08, 0.07, as function of ω. Let ω ₂(ε) be the interception points of the curves ϕ₁(ω) and ϕ ₂(ω). Thus ${\lambda }_2^{\epsilon }={\omega }_2{\left(\epsilon \right)}^2$.

Figure 2:
Graphical solution of tanh(ω) = εω for ε = 0:1; 0:09; 0:08; 0:07.

From Figs. 1 and 2 we can observe that the eigenvalues ${\lambda }_1^{\epsilon }$ and ${\lambda }_2^{\epsilon }$ increase continually when ε tends to zero, respectively.

Lemma 4.Let ε > 0. Then${\lambda }_1^{\epsilon }>{\lambda }_2^{\epsilon }$. Moreover${\lambda }_1^{\epsilon },{\lambda }_2^{\epsilon }\to \infty$when ε → 0.

Proof. We know that the functions tanh ω and coth ω can be write as

$\tanh \omega =1-2e^{-2\omega }+2e^{-4\omega }-2e^{-6\omega }+\mathrm{...}$ (2.13)

$\coth \omega =1+2e^{-2\omega }+2e^{-4\omega }+2e^{-6\omega }+\mathrm{...}$ (2.14)

From the equation (2.11), we have

$\omega =\frac{1}{\epsilon }\frac{1}{\tanh \omega }.$ (2.15)

When ω is large enough, tanh ω approaches 1. Thus, from (2.15) we have

$\omega =\frac{1}{\epsilon }\coth \frac{1}{\epsilon }=\frac{1}{\epsilon }\left(1+2e^{-2/\epsilon }+2e^{-4/\epsilon }+2e^{-6/\epsilon }+\mathrm{...}\right).$ (2.16)

Since ${\lambda }_1^{\epsilon }={\omega }_1{\left(\epsilon \right)}^2$, we have

${\lambda }_1^{\epsilon }=\frac{1}{{\epsilon }^2}+\frac{4}{{\epsilon }^2}{e}^{-2/\epsilon }+\frac{8}{{\epsilon }^2}{e}^{4/\epsilon }+\cdots .$ (2.17)

From (2.12) we obtain

$\omega =\frac{1}{\epsilon }\tanh \omega .$ (2.18)

When ω is large enough, tanh ω approaches 1. Thus, from (2.18) we have

$\omega =\frac{1}{\epsilon }\tanh \frac{1}{\epsilon }=\frac{1}{\epsilon }\left(1-2e^{-2/\epsilon }+2e^{-4/\epsilon }-2e^{-6/\epsilon }+\cdots \right).$ (2.19)

Since ${\lambda }_2^{\epsilon }={\omega }_2{\left(\epsilon \right)}^2$, we have

${\lambda }_2^{\epsilon }=\frac{1}{{\epsilon }^2}-\frac{4}{{\epsilon }^2}{e}^{-2/\epsilon }+\frac{8}{{\epsilon }^2}{e}^{-4/\epsilon }-\cdots .$ (2.20)

From (2.17) and (2.20) follow the results. □

Also from (2.17) and (2.20) we have that the gap between ${\lambda }_1^{\epsilon }$ and ${\lambda }_2^{\epsilon }$ is given by

${\lambda }_1^{\epsilon }-{\lambda }_2^{\epsilon }\approx \frac{8}{{\epsilon }^2}{e}^{-2/\epsilon }.$ (2.21)

Remark 5.When ε tends to zero, the problem (ℰ _ε) becomes

$\{\begin{array}{c}{u}_{xx}={\lambda }^{0}u,\begin{array}{c}\end{array}-1<x<\mathrm{1,}\\ u\left(-1\right)=u\left(1\right)=0.\begin{array}{cc}& \end{array}\end{array}$

(ℰ ₀)

We known that the eigenvalues of the problem (ℰ ₀) are given by${\lambda }_n^0=-n^2{\pi }^2,n\in \mathbb{N}$.

When ε tends to infinity, the problem (ℰ _ε) becomes

$\{\begin{array}{c}{u}_{xx}={\lambda }^{0}u,\begin{array}{c}\end{array}-1<x<\mathrm{1,}\\ u_x\left(-1\right)=u\left(1\right)=0.\begin{array}{cc}& \end{array}\end{array}$

(ℰ _∞)

3 THE FDM APPROACH

In this section we present the numerical schemes to solve (1.1) by applying the finite difference method (FDM) combined with the classic and unconditionally stable Crank-Nicolson method (see⁸).

To solve (1.1), the spatial domain [-1; 1] is discretized with uniform grid of n divisions of size h, where each spatial nodal points are x _i = ih - 1. Similarly, the temporal domain [0, T] is divided in m parts of size k, where T > 0 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 (1.1) at x = - 1 and for ε > 0 we write

$u_x\left(-\mathrm{1,}t\right)=\frac{-u\left(-\mathrm{1,}t\right)}{\epsilon }.$

Assuming that at x = -1, the function u is twice differentiable in x, so that we can write

$u_{xx}\left(-\mathrm{1,}t\right)=\frac{-u_x\left(-\mathrm{1,}t\right)}{\epsilon }=\frac{u\left(-\mathrm{1,}t\right)}{{\epsilon }^2},$

and in the same way at x =1 we get

$u_{xx}\left(\mathrm{1,}t\right)=\frac{u_x\left(\mathrm{1,}t\right)}{\epsilon }=\frac{u\left(\mathrm{1,}t\right)}{{\epsilon }^2}.$

Also, by using the finite difference approach for u _xx (t) (see ⁹), the problema (1.1) can be write

$\frac{d{u}_i\left(t\right)}{dt}=\frac{u_{i+1}\left(t\right)-2u_i\left(t\right)+u_{i-1}\left(t\right)}{h^2}$ (3.1)

for i = 1, ..., n - 1 , and for i = 0 and i = n, we have

$\frac{d{u}_0\left(t\right)}{dt}=\frac{u_0\left(t\right)}{{\epsilon }^2},$ (3.2)

and

$\frac{d{u}_n\left(t\right)}{dt}=\frac{u_n\left(t\right)}{{\epsilon }^2},$ (3.3)

respectively. Thus by using (3.1)-(3.3), the problem (1.1) can be transformed into the first-order matrix differential equation given by

$\frac{\partial }{\partial t}\left[\begin{array}{c}{u}_0\\ u_1\\ u_2\\ ⋮\\ u_{n-2}\\ u_{n-1}\\ u_n\end{array}\right]=\frac{1}{h^2}\left[\begin{array}{ccccccc}\frac{h^2}{{\epsilon }^2}& 0& 0& \cdots & 0& 0& 0\\ 1& -2& 1& \ddots & & & 0\\ 0& 1& -2& \ddots & \ddots & & 0\\ ⋮& \ddots & \ddots & \ddots & \ddots & \ddots & ⋮\\ 0& & \ddots & \ddots & -2& 1& 0\\ 0& & & \ddots & 1& -2& 1\\ 0& 0& 0& \cdots & 0& 0& \frac{h^2}{{\epsilon }^2}\end{array}\right]\left[\begin{array}{c}{u}_0\\ u_1\\ u_2\\ ⋮\\ u_{n-2}\\ u_{n-1}\\ u_n\end{array}\right],$ (3.4)

or in compact way

$\frac{dU}{dt}=\mathbf{L}U,$ (3.5)

where, the matrix L represents the discrete counterpart of the corresponding differential operator given in (1.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 (3.5): the explicit and implicit finite difference schemes. The type of implicit scheme adopted in this work was the Crank-Nicolson algorithm (see ⁸).

4 NUMERICAL RESULTS

For the numerical examples we consider the initial condition given by

$u\left(x\mathrm{,0}\right)=\cos \left(\frac{\pi }{2}x\right),$ (4.1)

where for homogeneous Dirichlet boundary condition, which is obtained from (1.1) with ε = 0, the analytical solution is

$u\left(x,t\right)=\exp \left(-\frac{{\pi }^2}{4}t\right)\cos \left(\frac{\pi }{2}x\right).$ (4.2)

The examples were solved using n = 20 divisions in x axis and m = 100 divisions in temporal axis t. In Fig. 3 we display the spatial behavior of the solution for t = 1.5 considering several values of ε: ε = 0.08, ε = 0.05, ε = 0.025, ε = 0.0075. The absolute error and the maximum error norm between analytical and numerical solution, u _a and u _n , are calculated and shown in Tables 1 and 2, respectively. In Fig. 4 we show the temporal behavior of the solution for the same values of ε used in the Fig. 3 but at x = 0.8. From these simulations, shown in Figs. 3 and 4 and in Tables 1 and 2, we can see that the numerical solutions of (1.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 ε tends to zero.

Figure 3:
Numerical solutions of (1.1) when ε → 0 at t = 1.5.

Table 1:
Absolute error ||u _a - u _n || in function of x calculated from the data of Fig. 3.

	0:0075	0:025	0:05	0:08
0	0:0012	0:0040	0:0076	0:0115
0:1000	0:0012	0:0039	0:0076	0:0114
0:2000	0:0012	0:0038	0:0076	0:0111
0:3000	0:0011	0:0036	0:0070	0:0105
0:4000	0:0010	0:0034	0:0065	0:0098
0:5000	0:0010	0:0031	0:0059	0:0088
0:6000	0:0009	0:0027	0:0052	0:0077
0:7000	0:0007	0:0023	0:0044	0:0064
0:8000	0:0006	0:0019	0:0034	0:0050
0:9000	0:0005	0:0014	0:0025	0:0034
1:0000	0:0003	0:0008	0:0014	0:0018

Table 2:
Maximum error norm of the difference of solutions u _a and u _n .

Simulation	ε	||u a - u n ||
1	0.0075	0.00121
2	0.0250	0.00395
3	0.0500	0.00765
4	0.0800	0.01151

Figure 4:
Numerical solutions of (1.1) when ε → 0 at x = 0.8.

The Figs. 5 and 6 display the temporal evolution of the solutions with the initial condition given by (4.1), the numerical solution for ε = 0.05 in Fig. 5 and analytical solution in 6. It is worth noting that the numerical solution presented in Fig. 5 for ε = 0.05 agrees very well with analytic solution of the problem with homogeneous Dirichlet boundary conditions presented in Fig. 6.

Figure 5:
Temporal evolution of the numerical solution obtained with initial condition given by (4.1) and for ε = 0.05.

Figure 6:
Temporal evolution of the analytical solution obtained in (4.2), with initial condition given by (4.1).

5 CONCLUSIONS

In the first part of the paper, we prove several results on the well-posedness of the system (1.1) and the associated stationary problem. In the second part, the application of the Crank-Nicolson for the numerical solution of (1.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 (1.1) for several values of ε. Analytical approaches and numerical simulations have clearly illustrated the asymptotic behaviour of the solution of (1.1) when ε tends to zero. Continuous dependence of the solutions of (1.1) from the ε parameter has also been demonstrated.

ACKNOWLEDGEMENTS

G. Lozada-Cruz is partially supported by the São Paulo Research Foundation (FAPESP, grant: 2015/24095-6). C.E. Rubio-Mercedes is supported by The Mato Grosso do Sul Research Foundation (FUNDECT, grant: 219/2016). J. Rodrigues-Ribeiro is partially supported by Coordination for the Improvement of Higher Education Personnel (CAPES, grant: PICME-9725290/M).

The authors are grateful to the referee for valuable remarks improving the original version of the paper.

REFERENCES

• ¹
  J.M. Arrieta, A.N. Carvalho & A. Rodríguez-Bernal. Attractors of parabolic problems with nonlinear boundary conditions. Uniform bounds. Comunications in Partial Differential Equations, 25(1-2) (2000), 1-37.
• ²
  T.L. Bergman, A.S. Lavine, I.F. P. & D.P. Dewitt. “Fundamentals of Heat and Mass Transfer”. 7th ed., John Wiley and Sons, Inc. (2007).
• ³
  H. Brezis. “Functional analysis, Sobolev spaces and partial differential equations”. Springer, New York (2011).
• ⁴
  G. Buttazzo, M. Giaquinta & S. Hildebrandt. “One-dimensional variational problems. An introduction”, volume 15. The Clarendon Press, Oxford University Press, New York (1998).
• ⁵
  L.D. Chiwiacowsky & H. de Campos Velho. Different Approaches for the Solution of a Backward Heat Conduction Problem. Inverse Problems in Engineering, (3) (2010), 471-494.
• ⁶
  R. Das, S. Mishra, T. Kumar & R. Uppaluri. An Inverse Analysis for Parameter Estimation Applied to a Non-Fourier Conduction-Radiation Problem. Heat Transfer Engineering, 32(6) (2011), 455-466.
• ⁷
  D. Henry. “Geometric theory of semilinear parabolic equations”, volume 840. Springer-Verlag, Berlim (1981).
• ⁸
  H.E. Hernández-Figueroa & C.E. Rubio-Mercedes. Transparent boundary for the finite element simulation of temporal soliton propagation. IEEE Transaction on Magnetics, 34(5) (1998).
• ⁹
  J.D. Hoffman. “Numerical Methods for Engineers and Scientists”. Marcel Dekker, Inc., New York, 2nd edition (2001).
• ¹⁰
  L.B.L. Santos, L.D. Chiwiacowsky & H.F. Campos-Velho. Genetic algorithm and variational method to identify initial conditions: worked example in hyperbolic heat transfer. Tendências em Matemática Aplicada e Computacional, 14(2) (2013), 265-276.

• †
  Work presented at the XXXVI National Congress of Applied and Computational Mathematics, Gramado, RS, Brazil, 2016.

Publication Dates

• Publication in this collection
  May-Aug 2018

History

• Received
  17 Apr 2017
• Accepted
  23 Nov 2017