# Numerical Solutions of Singularly Perturbed Reaction Diffusion Equation with Sobolev Gradients

Authors: Nauman Raza; Asma Rashid Butt
Journal: Journal of Function Spaces and Applications
 (2013)

---

## Abstract

Critical points related to the singular perturbed reaction diffusion models are calculated using weighted Sobolev gradient method in finite element setting. Performance of different Sobolev gradients has been discussed for varying diffusion coefficient values. A comparison is shown between the weighted and unweighted Sobolev gradients in two and three dimensions. The superiority of the method is also demonstrated by showing comparison with Newton's method.


---

## Body

## 1. Introduction

Many problems in biology and chemistry can be represented in terms of partial differential equations (PDEs). One of the important models is reaction diffusion problems. Much attention has been devoted for the solutions of these problems. In the literature it is shown that numerical solutions of these problems can be computed provided the diffusion coefficients, reaction excitations, and initial and boundary data are given in a deterministic way. The solution of these PDEs is extremely challenging when they has singularly perturbed behavior. In this paper, we discuss the numerical solutions of(1)pt=ϵΔp+p(1-p),p∈Ω,
where ϵ is a small and strictly positive parameter called diffusion coefficient, Ω is some two- or three-dimensional region. The Dirichlet boundary conditions are used to solve the equation. Numerous numerical algorithms are designed to solve these kind of systems [1, 2]. We are also using some numerical techniques to solve these systems based on the Sobolev gradient methods. A weighted Sobolev gradient approach is presented, which provides an iterative method for nonlinear elliptic problems.Weighted Sobolev gradients [3] have been used for the solution of nonlinear singular differential equations. It is shown that [4] significant improvement can be achieved by careful considerations of the weighting. Numerous Sobolev norms can be used as a preconditioner strategies. In Sobolev gradient methods, linear operators are formed to improve the condition number in the steepest descent minimization process. The efficiency of Sobolev gradient methods has been shown in many situations, for example, in physics [4–11], image processing [12, 13], geometric modelling [14], material sciences [15–20], Differential Algebraic Equations, (DAEs) [21] and for the solution of integrodifferential equations [22].We refer [23] for motivation and background for Sobolev gradients. For some applications and open problems in this area, an interesting article is written by Renka and Neuberger [24]. For the computational comparison an Intel Pentium 1.4 GHZ core i3 machine with 1 GB RAM was used. All the programs were written in FreeFem++ [25] which is freely available to solve these kind of problems.The paper is organized as follows. In Section2, we introduce the basic Sobolev gradient approach. This section contains some important results for the existence and convergence of the solution. In Section 3, we find the expression for Sobolev and weighted Sobolev gradients. Section 4 is composed of numerical results. Summary and conclusion are discussed in Section 5.


## 2. Sobolev Gradient Approach

This section is devoted to show the working of Sobolev gradient methods. A detailed analysis regarding the construction of Sobolev gradients can be seen in [23] and these lines are also taken from the same reference.Let us consider thatk is a positive integer and ψ is a real valued C1 function on Rk. The gradient ∇ψ is defined by
(2)limt→01tRk(ψ(x+th)-ψ(x))=ψ′(x)h=〈h,∇ψ(x)〉Rk,x,h∈Rk.
For ψ as in (2) but with 〈·,·〉S an inner product on Rk different from the usual inner product 〈·,·〉Rk, there is a function ∇sψ:Rk→Rk which satisfies
(3)ψ′(x)h=〈h,∇Sψ(x)〉S,x,h∈Rk.
The linear functional ψ′(x) can be represented using any inner product on Rk. Say that ∇Sψ is the gradient of ψ with respect to the inner product 〈·,·〉S and note that ∇Sψ has the same properties as ∇ψ.By applying a linear transformation:(4)A:Rk⟶Rk
we can relate these two inner products by
(5)〈x,y〉S=〈x,Ay〉Rk,x,y∈Rk,
and by applying a reflection we have
(6)(∇Sψ)(x)=A-1∇ψ(x),x∈Rk.
To each x∈Rk an inner product 〈·,·〉x is associated on Rk. Therefore, for x∈Rk, one can define ∇xψ:Rk→Rk such that
(7)ψ′(x)h=〈h,∇xψ(x)〉x,x,h∈Rk.
So corresponding to each metric, we can define a gradient for a function ψ and these gradients may have vastly different numerical properties. Therefore, the choice of metric plays a vital role in steepest descent minimization process. A gradient of a function ψ is said to be Sobolev gradient, when the underlying space is Sobolev. Readers who are unfamiliar with Sobolev spaces are referred to [26]. Steepest descent can be considered both in discrete as well as in continuous versions.Supposeψ is a real-valued C1 function defined on a Hilbert space H and ∇Sψ is its gradient with respect to the inner product 〈·,·〉S defined on H. Discrete steepest descent is a procedure of establishing a sequence {xk} such that for given x0 we have
(8)xj=xj-1-δj(∇ψ)(xj-1),j=1,2,…,
where for each j, δj is chosen in a way that it minimizes, if possible,
(9)ψ(xj-1-δj(∇ψ)(xj-1)).
In continuous steepest descent, a function z:[0,∞)→H is constructed so that
(10)dzdt=-∇ψ(z(t)),z(0)=zinitial.
Under suitable conditions on ψ, z(t)→z∞ where ψ(z∞) is the minimum value of ψ.So we conclude that discrete steepest descent (8) can be considered as a numerical method for approximating solutions to (10). Continuous steepest descent provides a theoretical starting point for proving convergence for discrete steepest descent.Using (8), one seeks p=limk→∞xj, so that
(11)ψ(p)=0or(∇Sψ)(p)=0.
Using (10), we seek p=limt→∞zt so that (11) holds.For the solution of PDEs, we can formulateψ by a variational principle, such that p satisfies the given PDE if and only if p is a critical point of ψ.To find a zero of the gradient ofψ we try to use steepest descent minimization process. The energy functional associated with (1) is given by
(12)ψ(p)=∫Ω‍δtp33+(1-δt)p22-fp+δtϵ2|∇p|2.
For the solution of PDEs, other functionals are also possible and one of the prime examples in this direction is the least square formulation. Such functions are shown in [11, 18]. The existence and convergence of z(t)→z(∞) for different linear and nonlinear forms of F is discussed in [23].Here we quote a result [23] for convergence of z(t) under a convexity assumption.Theorem 1 (see [23]).
LetH be a Hilbert space and let F be a nonnegative C2 function defined on H, for which there is ϵ>0 such that
(13)〈ψ′′(x)(∇ψ(x)),(∇ψ(x))〉≥ϵ∥∇ψ(x)∥H2,x∈H.
Suppose also that (10) holds and x∈H, ∇ψ(x)≠0. Then p=limt→∞zt exists and ∇ψ(p)=0.In this paper, only discrete steepest descents are used and for numerical computation, discretized versions of functionsψ are considered.


## 3. Gradients and Minimization

Each inner product corresponds a gradient and a descent direction. The performance of steepest descent minimization process can be improved by defining gradients in a suitably chosen Sobolev space. It is still an open problem, how to choose the best inner product which is suitable for the problem.A Sobolev spaceH1(Ω) is defined as
(14)H1(Ω)={p∈L2(Ω):Dαp∈L2(Ω),0≤α≤1},
where Dα denotes the weak derivative of p of order α. H1(Ω) is a Hilbert space with the norm defined by
(15)∥p∥H12=∫Ω‍|∇p|2+|p|2.By inspiration from Mahavier's work on weighted Sobolev gradients, we define a new norm as(16)∥p∥H1,w2=∫Ω‍|w∇p|2+|p|2.
This norm takes care of w=ϵδt which is affecting the derivative term in (24). We call the Sobolev space with the norm (16) a weighted Sobolev space.Now we show that the weighted Sobolev space with the norm defined by (16) is a Hilbert space denoted by H1,w. For this first, we show it satisfies the properties of inner product space and then we show it is complete.Let us denote the inner product〈·,·〉w associated with the norm defined by (16):
(17)〈p,p〉w=〈p,p〉+〈w∇p,w∇p〉,
where 〈·,·〉 is the usual L2 inner product in the vector space RM defined by 〈p,q〉=∑i‍p(i)q(i). Consider the following:
(18)〈p,p〉w=∫Ω‍|w∇p|2+|p|2>0
if p=0, then 〈p,p〉w=0. 						(19)〈p,q〉w=〈p,q〉+〈w∇p,w∇q〉〈p,q〉w=〈q,p〉+〈w∇q,w∇p〉=〈q,p〉w,〈p,q+q^〉w=〈p,q+q^〉+〈w∇p,w∇(q+q^)〉〈p,q+q^〉w=〈p,q〉+〈p,q^〉〈p,q+q^〉w+〈w∇p,w∇q〉+〈w∇p,w∇q^〉〈p,q+q^〉w=〈p,v〉+〈w∇p,w∇q〉〈p,q+q^〉w+〈p,q^〉+〈w∇p,w∇q^〉〈p,q+q^〉w=〈p,q〉w+〈p,q^〉w,〈αp,q〉w=〈αp,q〉+〈wα∇p,w∇q〉〈αp,q〉w=α〈p,q〉+α〈w∇p,w∇q〉〈αp,q〉w=α〈p,q〉w.
Let p1,p2,… be a Cauchy sequence in H1,w. Since L2 is complete so p∈L2 such that ∥pn-p∥→0 also ∥∇pn-∇p∥→0.Therefore(20)∥pn-p∥2=∫Ω‍|w∇(pn-p)|2+|pn-p|2∥pn-p∥2=∫Ω‍|w|2|∇pn-∇p|2+|pn-p|2⟶0.
Hence, H1,w is complete by the norm defined by (16).To incorporate the Dirichlet boundary conditions, we define a perturbation subspaceL02(Ω) of functions such that
(21)L02(Ω)={p∈L02(Ω):p=0onΓ}.
Here Γ denotes the boundary of the domain Ω. Perturbation subspaces related to H1 and H1,w would be H01=L02⋂‍H1 and H01,w=L02⋂‍H1,w, respectively.We need to find the gradient∇ψ(p) of a functional ψ(p) associated with the original problem and to find zero of the gradient.Given an energy functional:(22)ψ(p)=∫Ω‍R(p,∇p),
the Fréchet derivative of ψ(p) is a bounded linear functional ψ′(p) defined by
(23)ψ′(p)h=limα→0ψ(p+αh)-ψ(p)αforh∈H01(Ω).The energy functional in our case is given by(24)ψ(p)=∫Ω‍δtp33+(1-δt)p22-fp+δtϵ2|∇p|2,
then according to (23) we have
(25)ψ′(p)h=limα→01α∫Ω‍{δt(p+αh)33+(1-δt)(p+αh)22mmmmm,,,,-f(p+αh)+δtϵ2|∇(p+αh)|2mmmmm,,,-δtp33-(1-δt)p22+fp-δtϵ2|∇p|2(p+αh)33}.
Simplifying and applying the limit we have
(26)ψ′(p)h=∫Ω‍(δtp2+(1-δt)p-f)h-∫Ω‍ϵδt∇p·∇h.Let∇ψ(p), ∇ψs(p), and ∇ψw(p) denote the gradients in L2, H1, and H1,w, respectively. By using (7) we can write
(27)ψ′(p)h=〈∇ψ(p),h〉L2=〈∇ψs(p),h〉H1=〈∇ψw(p),h〉H1,w.
Thus the gradient in L2 is
(28)∇ψ(p)=(δtp2+(1-δt)p-f)h+ϵδt∇2p.To satisfy the boundary conditions, we are looking for gradients that are zero on the boundary ofΩ. So instead of using ∇ψ(p), we use π∇ψ(p). π is a projection that sets values of test function h zero on the boundary of the system. For implementing the Sobolev gradient method, a software FreeFem++ [25] is used, which is designed to solve PDEs using the finite element method. Therefore, to find π∇ψ(p) we need to solve the following:
(29)π∫Ω‍(δtp2+(1-δt)p-f)h-∫Ω‍ϵδt∇p·∇h=π∫Ω‍∇ψ(p)h+∫Ω‍∇∇ψ(p)·∇h.
In other words, for p in L2(Ω) find ∇ψ(p)∈L02, such that
(30)〈∇ψ(p),h〉L2=〈A(p)-f,h〉L2,∀h∈L02(Ω),(31)〈A(p),h〉L2=∫Ω‍(δtp2+(1-δt)p)h-∫Ω‍ϵδt∇p·∇h.This gradient suffers from the CFL conditions, as something well known for numerical analysts. By using (15) and (27) we can relate the L2 gradient and unweighted Sobolev gradient in the weak from as
(32)〈(1-∇2)∇sψ(p),h〉L2=〈A(u)-f,h〉L2.
Similarly using (16) and (27) the weighted Sobolev gradient can be related with the L2 gradient:
(33)〈(1-w2∇2)∇wψ(p),h〉L2=〈A(u)-f,h〉L2.
In order to find gradients, we represent the above systems in weak formulation as
(34)π∫Ω‍(δtp2+(1-δt)p-f)h-∫Ω‍ϵδt∇p·∇h=π∫Ω‍∇sψ(p)h+∫Ω‍∇∇sψ(p)·∇h,(35)π∫Ω‍(δtp2+(1-δt)p-f)h-∫Ω‍ϵδt∇p·∇h=π∫Ω‍∇wψ(p)h+∫Ω‍∇∇wψ(p)·∇h.It is seen that when using the weighted Sobolev gradient, the step sizeλ does not have to be reduced as the numerical grid becomes finer and the number of minimization steps remains reasonable. At the same time, the conceptual simplicity and elegance of the steepest descent algorithm has been retained.


## 4. Numerical Results

Let us consider a two-dimensional domainΩ in the form of a circle with centre at the origin and radius 12. An elliptic region is removed that has border x(t)=6cos(t), y(t)=2sin(t) with t∈[0,2π]. To start the minimization process, we specify the value of p=2.0 as initial state and the boundary condition was that p=1 on the outer boundary and p=-1 on the inner boundary on the boundaries. We also let ϵ=0.01 and time step δt=0.95.For the implementation of FreeFem++, one needs to specify the number of nodes on each border. After that a mesh is formed by triangulating the region.FreeFem++ then solves the equations of the type, as given by (35) which can be used to compute gradients in H1 and H1,w. We performed numerical experiments by specifying different numbers of nodes on each border. For each time step δt the functional defined by (24) was minimized using steepest descent steps with both H1 and H1,w until the infinity norm of the gradient was less than 10-6. The system was evolved over 15 time steps. The results are reported in Table 3.Table1 shows that by refining the mesh, the weighted Sobolev gradient becomes more and more efficient as compared to unweighted Sobolev gradients. It is also observed that by decreasing ϵ, the performance of steepest descent in H1,w is better than in H1.Table 1
Numerical results of steepest descent inH1, H1,w using δt=0.95, ϵ=0.01 for fifteen time steps in the two-dimensional case.
λ
Iterations
CPU
M
Triangles
H
1
H
1
,
w
H
1
H
1
,
w
H
1
H
1
,
w
—
—
0.8
0.035
948
793
6.68
2.52
10
24
0.8
0.035
1782
1137
45.29
14.87
20
152
0.8
0.035
2627
1674
126.05
40.7
30
322
0.8
0.035
2309
1621
215.16
155.84
40
348Let us consider now a three-dimensional domainΩ in the form of a cube with center at the origin and radius 8. To start the minimization process, we specify the value of p=0.0 as initial state and the boundary condition was set as p=1 on the top and bottom faces and p=-1 on the front back, left, and right faces of the cube. We also let ϵ=0.1 and time step δt=0.95. Once again, functional was minimized using steepest descent with both H1 and H1,w until the infinity norm of the gradient was less than 10-6. Once again, the same software is used. Numerical experiments are performed by varying number of nodes specified on each border. The results are recorded in Table 2.Table 2
Numerical results of steepest descent inH1, H1,w using δt=0.95 for fifteen time steps in the three-dimensional case.
λ
Iterations
CPU
M
Triangles
H
1
H
1
,
w
H
1
H
1
,
w
H
1
H
1
,
w
—
—
1.2
0.5
39.62
12.45
2293
462
5
50
1.2
0.5
4529.1
1487.2
3302
474
10
200
1.2
0.5
25214.3
8140.3
3986
499
20
800Table 3
Comparison of Newton's method and steepest descent inH1,w for different values of ϵ.
ϵ
=
1.0
ϵ
=
0.1
ϵ
=
0.01
—
Newton
H
1
,
w
Newton
H
1
,
w
Newton
H
1
,
w
Error
10
13
30
10
NC
13
1
0
-
5
11
20
38
16
NC
22
1
0
-
7
11
28
NC
22
NC
34
1
0
-
9By refining the mesh with increasingM, the efficiency of weighted Sobolev gradient in terms of minimization steps taken for convergence and CPU time also increased. By reducing the value of ϵ, the performance of steepest descent in H1,w is more pronounced over than over H1.


## 5. Comparison with Newton's Method

In this section, a comparison is given between Newton's method and the Sobolev gradient method is with weighting. Newton's methods considered as one of the optimal methods to solve these kinds of problems. For convergence of Newton's method a good initial guess is required. For the comparison between the two methods, we take a good initial guess for implementation of Newton's method so that we fairly can compare the two methods. In variational form, the given nonlinear problem can be written as(36)〈ψ(p),h〉=∫Ω‍δtp33h+(1-δt)p22h-fph+δtϵ2|∇p|2h.
To apply Newton's method, we need to find Gâteaux derivative which is defined by
(37)〈F′(un)cn,v〉=〈F(p),v〉,∀v∈H01(Ω).
A linear solver is required to solve (37). Thus the Newton's iteration scheme looks like
(38)un+1=un-cn.An example is solved in two-dimensional case. We letΩ be in the form of a circle with centre at the origin and radius 12. An elliptic region is removed that has border x(t)=6cos(t), y(t)=2sin(t) with t∈[0,2π]. The initial state was p=2.0 and the boundary condition was that p=1 on the outer boundary and p=-1 on the inner boundary on the boundaries. The value of time step was set as δt=0.95. The system was evolved over 1 time step. The minimization was performed until the infinity norm of the gradient was less than some set tolerance. We specify 20 nodes on each border to obtain results.Different values ofϵ results are recorded in Table 3. The term NC denotes no convergence. From the table, one can observe that Newton's method perform better than the steepest descent in H1,w. For strict stopping criterion, Newton's method does not converge. Newton's method was superior when the minimization process start but not able to handle low frequency errors. On the other hand, steepest descent in H1,w takes more iterations but it continues to converge even for very tight stopping criteria. By decreasing the value of diffusion coefficient ϵ, the steepest descent in H1,w always manage to converge whereas this is not in the case of Newton's method.


## 6. Summary and Conclusions

In this draft, for the solution of concerned problem, a minimization scheme based on the Sobolev gradient methods [23] is developed. The performance of descent in H1,w is better than descent in H1, as the spacing of the numerical grid is refined. In this paper, we signify a systematic way to choose the underlying space and the importance of weighting to develop efficient codes.The convergence of Newton's method depend on the good initial guess which is sufficiently close to local minima. To implement Newton's method, we need to evaluate the inverse of Hessian matrix which is expensive at some times. A failure of the Newton's method has been shown in [11, 18] by taking large size of jump discontinuity and by refining the grid. But this is not in the case of steepest descent method. So a broader range of problems can be addressed by using steepest descent in appropriate Sobolev space.One of the advantages of thge Sobolev gradient methods is that it still manages to converge even if we take rough initial guess. The weighted Sobolev gradient offers robust alternative method for problems which are highly nonlinear and have large discontinuities. An interesting project can be done by comparing the performance of Sobolev gradient methods with some multigrid technique in some finite element setting.



---

*Source: 542897-2013-11-06.xml*