A dual mutation differential evolution algorithm for singularly perturbed problems with two small parameters

Abstract

A novel high accuracy numerical method for a singularly perturbed convection-diffusion problem with two small parameters is presented. At first, the given problem is discretized by using a rational spectral collocation method in barycentric form with sinh transformation. By sinh transform, the original Chebyshev points are mapped into the transformed ones clustered near to the boundary layers of the domain. Them, a nonlinear unconstrained optimization problem is designed to determine the widths of boundary layers in sinh transform. Finally, a dual mutation differential evolution (DMDE) algorithm is proposed to solve the nonlinear unconstrained optimization problem, and a comparison of the DMDE algorithm with other algorithms has been made, which shows that DMDE algorithm can get fast convergence and find better results.

Keywords

Singularly perturbed convection-diffusion problems rational spectral collocation method optimization problem differential evolution algorithm

1. Introduction

Singularly perturbed problems arise in various branches of engineering and applied mathematics such as geophysical fluid dynamics, astrophysics, oceanography, meteorology, pollutant and sediment transport, and chemical reactor, etc. These problems described by differential equations contain one or more small positive parameters. Due to the presence of these parameters, the solutions of these problems possess boundary layers which are basically thin regions in the neighborhood of the boundary of the domain. The numerical methods of singularly perturbed problems with a small parameter have attracted much attention over the past decades, see the monographs [1, 2].

In recent years, there has been tremendous interest in developing numerical methods for singularly perturbed convection-diffusion problems with two small parameters, which is given by ${\begin{array}{l} - ε u^{″} (x) - μ a (x) u^{'} (x) + b (x) u (x) = f (x), \\ u (0) = A, u (1) = B, \end{array}$ (1) where x ∈ (0, 1), ɛ, μ are two small parameters, and A, B are two given constants. We assume that the functions a (x), b (x), f (x) are all in C¹ [0, 1], and a (x) and b (x) satisfying $\begin{matrix} a (x) \geq σ > 0, b (x) \geq γ > 0, \\ \forall x \in [0, 1], ρ = min_{0 \leq x \leq 1} \frac{b (x)}{a (x)} . \end{matrix}$

Under these assumptions, the problem (1) has a unique solution, which exhibits boundary layers at both ends of the interval [0, 1]. As a consequence, standard finite difference or finite element methods with uniform meshes applied to (1) may yield inaccurate results. One effective technique to obtain accurate numerical results is to use some suitable nonuniform meshes that are very fine where layers appear in the solution. Such nonuniform meshes can be divided into two classes. One is the use of highly non-uniform layer-adapted meshes (called Shishkin mesh and Bakhvalov mesh), see, e.g., [3 –7]. Another is the use of adaptive mesh generated by equidistributing a monitor function over the domain of the problem, see, e.g., [8, 9]. It should be pointed out that these methods used to solve (1) only obtain first-order or second order rate of convergence, independent of perturbation parameters. Thus, it is very necessary to design a high accuracy numerical method to solve problem (1).

Compared with finite difference and finite element methods, spectral collocation methods are more suitable to solve singularly perturbed boundary layer problems due to the fact that spectral collocation points are naturally well-adapted to the representation of boundary layers. Take the example of Chebyshev-Gauss-Lobatto points x_k = cos(kπ/N), we have $\begin{matrix} \begin{matrix} | x_{1} - x_{0} | & = | x_{N} - x_{N - 1} | = | cos (π / N) - 1 | \\ \approx \frac{1}{2} {(\frac{π}{N})}^{2} \approx \frac{5}{N^{2}} . \end{matrix} \end{matrix}$ It means that the distance between the points near the boundaries is of order O (N^-2). However, when the parameter ɛ tends to 0, still large N is required to get reliable numerical results. Therefore, some improved spectral collocation methods were proposed in [10, 11]. In these methods, the authors constructed a certain transformation x ↦ x^t = t (x), which made the collocation points more clustered in the boundary layer region. However, as pointed out in [10], the derivatives in the underlying differential equation have to be transformed into new coordinates after such a transformation, which would greatly increase the computational complexity. Recently, the authors in [12, 13] used rational spectral method in barycentric form with sinh transform to solve some singularly perturbed problems, where, ∀x ∈ [-1, 1], the sinh transform was defined as follows: $\begin{matrix} g (x) = λ + δ sinh [({sinh}^{- 1} (\frac{1 - λ}{δ}) \\ + {sinh}^{- 1} (\frac{1 + λ}{δ})) \frac{x - 1}{2} + {sinh}^{- 1} (\frac{1 - λ}{δ})] . \end{matrix}$ (2) Here, the parameters λ and δ represent the location and width of the boundary layer, respectively. In general, the location of the boundary layer λ can be determined by using the prior information of solution. For the width of boundary layer, Chen, et.al., [12, 13] chosen δ = βɛ, where β is a positive constant. Recently, Liu, et.al., [14] also used the rational spectral collocation method to solve the singularly perturbed problems with an interior layer, and used differential evolution to obtain the optimal width of the interior layer. It should be pointed out that the sinh transform given in (2) is only suitable to solve some singularly perturbed problems with one boundary layer. Thus, the scholars in [15] presented a novel collocation method for a coupled system of singularly perturbed boundary value problems with two boundary layers, where the new sinh transform was given by ${\tilde{g}}_{δ} (x) = {\begin{matrix} 1 / 2 [g_{- 1, δ} (2 x + 1) - 1], x \in [- 1, 0), \\ 1 / 2 [g_{1, δ} (2 x - 1) + 1], x \in [0, 1] . \end{matrix}$ For the width of boundary layer, they chosen $δ = - 2 \frac{ln \sqrt{ɛ}}{α} \sqrt{ɛ}$ , where ɛ is the perturbation parameter and α is a positive constant. It is noted that the method presented in [15] is unsuited to solve the singularly perturbed problems with different widths of boundary layers.

The purpose of this paper is to study a new rational spectral collocation method in barycentric form for the above singularly perturbed convection-diffusion problems with two small parameters (1). Obviously, for different parameters ɛ and μ, the solution of problem (1) have two boundary layers near the boundary of the domain. Furthermore, the widths of these boundary layers are different. Therefore, in order to solve the problem (1) by using rational spectral collocation method, we shall define a modified sinh transform as follows: ${\tilde{g}}_{δ_{1}, δ_{2}} (x) = {\begin{matrix} 1 / 2 [g_{- 1, δ_{1}} (2 x + 1) - 1], x \in [- 1, 0), \\ 1 / 2 [g_{1, δ_{2}} (2 x - 1) + 1], x \in [0, 1] . \end{matrix}$ (3) where δ_i, i = 1, 2 represent the width of two boundary layers, respectively. Here, in order to determine the widths of the two boundary layer, we choose δ₁ = θ₁ɛ, δ₂ = θ₂μ, where θ_i (i = 1, 2) are two constants. Obviously, it is very important to obtain the optimal parameters θ_i. For this reason, this paper will design a new swarm intelligence algorithm to calculate the optimal parameters θ_i.

As far as we know, research community has implemented a variety of evolutionary intelligence algorithms for finding accurate and reliable solution of singularly perturbed problems and two point boundary value problems. For instance, singularly perturbed problems with Robin boundary conditions [16], two-point boundary value problems [17], singularly perturbed convection-diffusion problem with two small parameters [18].

It is well known that DE algorithm is a fast and simple yet powerful population-based stochastic search technique, which is an efficient and effective global optimizer in the continuous search domain. It is a relatively new nonlinear search and optimization method, which is suitable to solve complicate optimization problems. Recently, optimization technique of DE has been used in a number of optimization problems in various fields such as mechanical engineering [19], signal processing [20], pattern recognition [21], parameter estimation problems [22, 23].

In this paper, a dual mutation differential evolution algorithm is proposed to solve the singularly perturbed convection-diffusion problems with two small parameters. Firstly, the given problem will be discretized by using rational spectral collocation method in barycentric form with sinh transform (3). Then, the original problem (1) is posed as an optimization problem, where the goal is to minimize the fitness function about absolute error of the numerical solution, which is solved by using the dual mutation differential evolution. Finally, we can obtain the optimal parameters θ_i, i = 1, 2, and the corresponding numerical solution of problem (1).

2. Rational spectral collocation method in barycentric form

2.1 The barycentric form of rational interpolation

Assume that v (x) is a sufficiently smooth in the interval [-1, 1]. Let p_N (x) be a rational function in barycentric form which interpolates v (x) at N + 1 distinct points ${x_{k}}_{k = 0}^{N}$ . Then, we have $p_{N} (x) = \sum_{j = 0}^{N} L_{j} (x) v (x_{j}),$ (4) where $\begin{matrix} L_{j} (x) = \frac{\frac{ω_{j}}{x - x_{j}}}{\sum_{k = 0}^{N} \frac{ω_{k}}{x - x_{k}}} \end{matrix}$ are the interpolation basis functions and ${ω_{k}}_{k = 0}^{N}$ are nonzero numbers called barycentric weights. Here, we choose the Chebyshev-Gauss-Lobatto points x_k = cos(kπ/N), k = 0, 1, ⋯ , N as the collocation points for our presented rational spectral collocation method, and the corresponding barycentric weights are chosen as follows: $ω_{0} = \frac{1}{2}, ω_{N} = \frac{(- 1)^{N}}{2}, ω_{k} = (- 1)^{k},$ (5) where k = 1, 2, ⋯ , N - 1.

By simple calculation, the nth order derivatives of p_N (x) at the point x_i can be given as $\begin{matrix} p_{N}^{(n)} (x_{i}) = \sum_{j = 0}^{N} L_{j}^{(n)} (x_{i}) v (x_{j}) = \sum_{j = 0}^{N} D_{ij}^{(n)} v (x_{j}), \end{matrix}$ where n = 1, 2, ⋯, ${D_{ij}^{(n)}}_{i, j = 0}^{N}$ are the nth order differentiation matrices. The advantage of the barycentric form is the simplicity of the formula for the entries of first and second order differential matrices which are given by $D_{jk}^{(1)} = {\begin{matrix} \frac{ω_{k}}{ω_{j} (x_{j} - x_{k})}, j \neq k, \\ - \sum_{i \neq j} D_{ji}^{(1)}, j = k, \end{matrix}$ (6) $D_{jk}^{(2)} = {\begin{matrix} 2 D_{jk}^{(1)} (D_{jj}^{(1)} - \frac{1}{x_{j} - x_{k}}), j \neq k, \\ - \sum_{i \neq j} D_{ji}^{(2)}, j = k . \end{matrix}$ (7)

2.2 Numerical discretization

First, by introducing the transformation $x = \frac{1}{2} (y + 1)$ and defining $\tilde{u} (y) = u (x) = u (\frac{1}{2} (y + 1))$ , then $2 {\tilde{u}}^{'} (y) = u^{'} (x), 4 {\tilde{u}}^{″} (y) = u^{″} (x)$ , we can rewrite the problem (1) as ${\begin{array}{l} - 4 ε {\tilde{u}}^{″} (y) - 2 μ \tilde{a} (y) {\tilde{u}}^{'} (y) + \tilde{b} (y) \tilde{u} (y) = \tilde{f} (y), \\ \tilde{u} (- 1) = A, \tilde{u} (1) = B, \end{array}$ (8) where $\begin{matrix} \tilde{a} (y) = a (\frac{1}{2} (y + 1)), \\ \tilde{b} (y) = b (\frac{1}{2} (y + 1)), \\ \tilde{f} (y) = f (\frac{1}{2} (y + 1)) . \end{matrix}$

Obviously, it follows from Lemma 1 of literature [9] that the problem (1) have two boundary layers at x = 0 and x = 1. Thus, the solution of above Equation (8) also have two boundary layers at y = -1 and y = 1.

As is stated in [12, 13], the rational spectral collocation method based on Chebyshev-Gauss-Lobatto points y_k = cos(kπ/N) can’t obtain satisfactory numerical results for very small ɛ and μ. Therefore, we utilize the sinh transform given in Equation (3) to obtain the following transformed Chebyshev points ${\hat{y}}_{k} = {\tilde{g}}_{δ_{1}, δ_{2}} (cos (k π / N) \in [- 1, 1] .$ (9)

Next, let S_N (y) be a rational function in barycentric form which interpolates exact solution $\tilde{u} (y)$ of problem (8) at the above transformed Chebyshev points ${\hat{y}}_{k}$ , then we have $S_{N} (y) = \sum_{j = 0}^{N} L_{j} (y) \tilde{u} ({\hat{y}}_{j}) .$ (10)

Substituting Equation (10) into Equation (8), we get $\begin{matrix} - 4 ɛ \sum_{j = 0}^{N} {L_{j}}^{″} (y) {\tilde{U}}_{j} - 2 μ \tilde{a} (y) \sum_{j = 0}^{N} {L_{j}}^{'} (y) {\tilde{U}}_{j} \\ + \tilde{b} (y) \sum_{j = 0}^{N} L_{j} (y) {\tilde{U}}_{j} = \tilde{f} (y), \end{matrix}$ (11)

where ${\tilde{U}}_{j}$ is the approximation solution of $\tilde{u} (y_{j})$ .

Let $y = {\hat{y}}_{k}$ , k = 0, 1, ⋯ , N, then Equation (11) can be written as follows: $[- 4 ɛ D^{(2)} - 2 μ {PD}^{(1)} + Q] U = F,$ (12)

where $\begin{matrix} P = diag (\tilde{a} (y_{0}), \dots, \tilde{a} (y_{N})), \\ Q = diag (\tilde{b} (y_{0}), \dots, \tilde{b} (y_{N})), \\ F = {(\tilde{f} (y)_{0}, \tilde{f} (y)_{1}, \dots, \tilde{f} (y)_{N})}^{T}, \\ U = {({\tilde{U}}_{0}, {\tilde{U}}_{1}, \dots, {\tilde{U}}_{N})}^{T} . \end{matrix}$

By defining K = -4ɛD⁽²⁾ - 2ωPD⁽¹⁾ + Q, we have $K U = F .$ (13)

Given two row vectors with length N + 1: e₁ = (1, 0, ⋯ , 0) , e_N+1 = (0, 0, ⋯ , 1) and replaced the first and last rows of matrix K with e₁ and e_N+1, respectively, we can obtained the new matrix $\tilde{K}$ . Similarly, the ${\tilde{f}}_{0}$ and ${\tilde{f}}_{N}$ should be replaced by A and B to get the new vector $\tilde{F}$ . Finally, we obtain the following linear system $\tilde{K} U = \tilde{F} .$ (14)

2.3 Fitness function

Here, for given parameters θ_i, i = 1, 2, by using conjugate gradient method to solve the above linear system (14), we can obtain the numerical results of problem (8). Obviously, for different parameters θ_i, we can get different numerical solution of problem (8). Therefore, in order to obtain the optimal parameters θ_i, we construct a fitness function about parameters θ_i as follows:

To give the absolute error formula, we first define two grids ${\bar{Ω}}^{N}$ and ${\bar{Ω}}^{2 N}$ as follows: ${\bar{Ω}}^{N} = {y_{k}}_{k = 0}^{N} = {g (cos (k π / N))}_{k = 0}^{N},$ (15) ${\bar{Ω}}^{2 N} = {z_{j}}_{j = 0}^{2 N} = {g (cos (j π / 2 N))}_{j = 0}^{2 N} .$ (16) It follows from (15)-(16) that the mesh ${\bar{Ω}}^{2 N}$ can be obtained by bisecting ${\bar{Ω}}^{N}$ . Thus, for 0 ≤ i ≤ N, we have y_i = z_2i. Next, let ${\tilde{U}}^{N} (y_{k})$ and ${\tilde{U}}^{2 N} (z_{j})$ be the solutions to the discrete problem (14) computed on ${\bar{Ω}}^{N}$ and ${\bar{Ω}}^{2 N}$ , respectively. Then, we use the following formula to estimate the maximum point-wise error $E_{ɛ, μ}^{N} = max_{0 \leq j \leq N} | {\tilde{U}}^{N} (y_{j}) - {\tilde{U}}^{2 N} (z_{2 j}) | .$ (17) Obviously, for each ɛ, μ and N, the maximum point-wise error $E_{ɛ, μ}^{N}$ is a function about variables θ_i, i = 1, 2. In practical computation, we hope that the value of $E_{ɛ, μ}^{N}$ as small as possible. So, we may construct the following fitness function $E_{ɛ, μ}^{N} (θ_{1}, θ_{2}) = max_{0 \leq j \leq N} | U^{N} (y_{j}) - U^{2 N} (z_{2 j}) | .$ (18) In conclusion, the original numerical solution of singularly perturbed problem (8) is now posed as an optimization problem, where the goal is to minimize the above fitness function (18). It should be pointed out that the above fitness function (18) is an implicit function about variables θ_i, i = 1, 2. Obviously, it is hard to derive the partial derivative of function $E_{ɛ}^{N} (θ_{1}, θ_{2})$ . So, it is difficulty to solve it by using some traditional optimization methods. In the next section, we will give a dual mutation differential evolution algorithm to solve the above problem (18).

3. A dual mutation differential evolution (DMDE) algorithm

DE is good at exploring the search space and locating the region of global minimum. However, it is slow exploiting of the solution [24]. In order to balance the exploration and the exploitation of DE, we propose a dual mutation differential evolution (DMDE) algorithm, which combines two different mutation operators in one algorithm. Our proposed DMDE algorithm is described as follows:

A. Mutation

For each target vector X_i,g, two kinds of mutation operations are generated. The first mutation is focused on exploration. The mutant vector V_i,g is generated according to $V_{i, g + 1} = X_{r_{1}, g} + F * (X_{i, g} - X_{r_{2}, g}),$ (19) where i = 1, ⋯ , NP, r₁, r₂ ∈ {1, ⋯ , NP} are randomly selected and satisfy: r₁ ≠ r₂ ≠ i, F ∈ [0, 2] is a scaling factor, X_i,g is known as the target vector, and V_i,g is known as the donor vector in the population at generation g.

The second mutation is focused on exploration. The mutant vector $V_{i, g}^{'}$ is generated according to $V_{i, g + 1}^{'} = X_{i, g} + γ L (λ) * (X_{best, g} - X_{i, g}) .$ (20) Here, γ is a scaling factor to control the step size, L (λ) is the Lévy-flights-based step size, the definition refers to [25], X_best,g is the best individual vector with the best fitness in the population at generation g.

B. Crossover

The target vector is mixed with the mutated vector, using the following scheme, to yield the trial vector $U_{i, g + 1} = (u_{i 1, g + 1}, u_{i 2, g + 1}, \dots, u_{iD, g + 1}),$ (21) where $u_{ij, g + 1} = {\begin{matrix} v_{ij, g + 1}, if r (j) \leq CR or j = rn (i), \\ v_{ij, g + 1}^{'}, if r (j) > CR and j \neq rn (i), \end{matrix}$ for j = 1, 2, ⋯ , D . r (j) ∈ [0, 1] is the j-th evaluation of a uniform random generator number. CR ∈ [0, 1] is the crossover constant, which has to be determined by the user. rn (i) ∈ (1, 2, ⋯ , D) is a randomly chosen index which ensures that U_i,g+1 gets at least one element form V_i,g+1.

C. Selection

If, and only if, the trial vector U_i,g+1 yields a better cost function value than X_i,g, then X_i,g+1 is set to U_i,g+1; otherwise, the old value X_i,g is retained.

According to Storn et al. [26], DE algorithm is very sensitive to the choice of F and CR. The crossover rate plays a very important role on the performance of DMDE algorithm and using it dynamically makes the algorithm more reliable in exploring and exploiting the search space. The general equation for dynamic switch probability CR used here is given by: $CR = a - \frac{g}{g_{\max}} * b,$ (22) where a and b are two constants, which satisfy 1 ≥ a ≥ b > 0. g is the current iteration and g_max is the total number of iterations.

Besides, to keep the diversity of the population, a strategy is used here, which is the same as one in the artificial bee colony (ABC) algorithm [28]. In DMDE algorithm, if the fitness value of a individual keeps unchanged for at least limit times, this individual is considered to be abandoned, and a new individual is generated to replace it as follows. $x_{i}^{j} = x_{\min}^{j} + rand (0, 1) (x_{\max}^{j} - x_{\min}^{j}),$ (23) where $[x_{\min}^{j}, x_{\max}^{j}]$ is the boundary constraint for the jth dimensional variable, and rand (0, 1) is a uniformly distributed random number within the range [0, 1].

The pseudocode of DMDE algorithm is described as follows:

Begin
Initialize the population X_i, i = 1, ⋯ , NP;
Calculate the objective function value f (X_i) for
all X_i; g = 1;
Repeat
Calculate the CR by (22).
fori = 1 to NPdo
Generate 2 random integers,r₁, r₂∈
1, 2, ⋯ , NP, r₁ ≠ r₂ ≠ i;
Generate a random integer j_rand∈
1, 2, ⋯ , D;
forj = 1 to Ddo
ifrand (0, 1) < CR or j_rand = = jthen
Generated the mutant variable v_ij by
(19) and the trial variable u_ij by (21)
else
Generated the mutant variable $v_{ij}^{'}$ by
(20) and the trial variable u_ij by (21)
end if
Check whether the trial variable u_ij is
within range
end for
Calculate the objective function value for
vector U_i
Replace X_i with the trial individual U_i, if U_i is better
Generated new individuals by (23) for those whose limit to
be improved has been reached
end for
g = g + 1
Until g=Maxgeneration

4. Numerical experiments and discussion

4.1 Test example and error evaluation approach

In this section, we consider the following test example to demonstrate the effectiveness of the proposed method $- ɛ u^{″} (x) - μ u^{'} (x) + u (x) = exp (x - 1),$ (24) $u (0) = 1, u (1) = 0,$ (25) where -1 < x < 1.

Let U^N and U^2N be the optimal numerical solutions to (24)-(25) calculated on meshes ${\bar{Ω}}^{N}$ and ${\bar{Ω}}^{2 N}$ , respectively. Then, we use the following formula to estimate the maximum error $E_{ɛ, μ}^{N} = U^{N} - U^{2 N}_{\infty} .$ (26)

4.2 Experimental environment and parameter settings

All algorithms are implemented in Matlab2014a to program a m-file for implementing the algorithms on a PC with a 32-bit windows 8 operating system, a 8GB RAM, and a 3.10 GHz-core (TM) i5-based processor.

The DMDE algorithm is validated by comparing it with five optimization algorithms, namely, flower pollination algorithm(FPA) [25], differential evolution (DE) [26], particle swarm optimization (PSO) [27], artificial bee colony (ABC) [28], covariance matrix adaption evolution strategy (CMA-ES) algorithm [29]. For the sake of fairness, we use a fixed population size for all algorithm: n = 20 individuals; all algorithms are executed in 30 independent runs; the number of iterations in each run for each algorithm equals 100 iterations; the range of parameters θ_i are [1, 100]. The details about parameters settings are listed as follows:

DMDE: F = 0.5, γ = 0.5,limit=10;

DE: F = 0.5, CR = 0.8;

FPA: switch probability P = 0.8, γ = 0.5, λ = 1.5;

PSO: learning factors c₁ = c₂ = 2, inertia factor dropped from 0.9 to 0.4;

ABC: limit=10;

CMA-ES: standard deviation sigma=0.5.

4.3 Experimental results and analyses

In this section, to illustrate the DMDE algorithm’s effectiveness, comparison experiments are done using FPA,DE,PSO,ABC,CMA-ES and DMDE algorithms. At first, when ɛ = 10^-10, μ = 10^-8 and N = 8, 16, 32, 64, 128, the statistical datas about the maximum, minimum, average and variances of maximum errors

E_{ɛ, μ}^{N}

are given in Tables 1–5, respectively. Obviously, when N = 8, 16, 32, 64, the variances of maximum errors E_ɛ,μ calculated by DMDE algorithm is smallest than that of other algorithms. When N = 128, the numerical results of DMDE, FPA, and PSO algorithms are almost the same. Meanwhile, the variances of maximum errors E_ɛ,μ obtained by DMDE algorithm is much lower than that of ABC,CMA-ES and DE algorithms. In all, DMDE algorithm works very well in the numerical results for problem (24)-(25), and DMDE algorithm has higher precision and stronger robustness by comparing with the FPA,DE,PSO,CMA-ES and ABC algorithms.

Table 1
Numerical results by using different algorithms (N = 8)

Method	Maximum error	Minimum error	Mean value	Variance
PSO	1.2569e-02	4.9412e-03	6.3392e-03	1.6946e-03
FPA	5.0878e-03	4.9374e-03	4.9676e-03	3.6388e-05
ABC	4.9382e-03	4.9371e-03	4.9372e-03	2.0772e-07
CMA-ES	5.2319e-03	4.9371e-03	5.0452e-03	1.4447e-04
DE	5.2319e-03	4.9371e-03	5.0021e-03	1.2092e-04
DMDE	4.9371e-03	4.9371e-03	4.9371e-03	2.2592e-12

Table 2

Numerical results by using different algorithms(N = 16).

Method	Maximum error	Minimum error	Mean value	Variance
PSO	1.0990e-02	8.3919e-03	9.2936e-03	5.8662e-04
FPA	8.9689e-03	7.6169e-03	8.5631e-03	4.4938e-04
ABC	8.9168e-03	7.9076e-03	8.8154e-03	2.7783e-04
CMA-ES	9.2964e-03	8.9168e-03	9.0052e-03	1.6309e-04
DE	8.9168e-03	7.4759e-03	8.8687e-03	2.6307e-04
DMDE	8.9168e-03	7.4759e-03	8.3884e-03	4.9874e-07

Table 3

Numerical results by using different algorithms(N = 32).

Method	Maximum error	Minimum error	Mean value	Variance
PSO	2.8571e-03	2.3057e-03	2.5006e-03	1.7809e-04
FPA	2.4788e-03	2.3043e-03	2.3541e-03	5.3971e-05
ABC	2.3191e-03	2.3012e-03	2.3022e-03	3.2786e-06
CMA-ES	3.9895e-03	2.3012e-03	2.6042e-03	3.5566e-04
DE	2.6965e-03	2.3012e-03	2.4277e-03	1.6231e-04
DMDE	2.3012e-03	2.3012e-03	2.3012e-03	9.8535e-12

Table 4

Numerical results by using different algorithms(N = 64).

Method	Maximum error	Minimum error	Mean value	Variance
PSO	6.7558e-05	6.3524e-05	6.3889e-05	7.8329e-07
FPA	6.3639e-05	6.3518e-05	6.3542e-05	2.5643e-08
ABC	6.3518e-05	6.3517e-05	6.3518e-05	1.5783e-10
CMA-ES	2.5946e-04	6.3517e-05	9.9498e-05	5.0959e-05
DE	6.3517e-05	6.3517e-05	6.3517e-05	5.3278e-14
DMDE	6.3517e-05	6.3517e-05	6.3517e-05	4.0599e-15

Table 5

Numerical results by using different algorithms(N = 128).

Method	Maximum error	Minimum error	Mean value	Variance
PSO	1.0184e-07	1.0184e-07	1.0184e-07	0
FPA	1.0184e-07	1.0184e-07	1.0184e-07	0
ABC	1.0184e-07	1.0184e-07	1.0184e-07	5.082e-15
CMA-ES	8.1893e-07	1.0184e-07	2.0606e-07	2.0334e-07
DE	1.0192e-07	1.0184e-07	1.0184e-07	1.1048e-11
DMDE	1.0184e-07	1.0184e-07	1.0184e-07	0

In addition, to test the convergence speed of the presented DMDE algorithm, the convergence profiles are given in Figs. 1 and 2. These Figures show the superior performance of DMDE algorithm as compared to other algorithms.

Fig. 1

Convergence profiles of various algorithms with ɛ = 10^-10, μ = 10^-8 and N = 16.

Fig. 2

Convergence profiles of various algorithms with ɛ = 10^-10, μ = 10^-8 and N = 32.

5. Conclusions

In this paper, a new high accuracy numerical method based on the rational spectral collocation method and a dual mutation differential evolution (DMDE) algorithm have been developed for solving singularly perturbed convection-diffusion problems with two small parameters. The key to the success of this approach is that the DMDE algorithm is presented to determine the widths of boundary layers in the sinh transform of the rational spectral collocation method in barycentric form. For this reason, we first construct a nonlinear unconstrained optimization problem based on minimum absolute error. Then, the DMDE algorithm is proposed to solve the presented optimization problem. The details of the new presented method is easy to implement on a computer. It is noted that the method presented in this paper can be extend to other type of singularly perturbed problems.

Footnotes

Acknowledgments

This work is supported by National Science Foundation of China (11761015), the National Natural Science Foundation Mathematics Tianyuan Foundation of China (11826211, 11826212), the Natural Science Foundation of Guangxi (2017GXNSFBA198183), the key project of Guangxi Natural Science Foundation (2017GXNSFDA198014), the key pro-ject of Anhui natural science research (KJ2016A514), the Natural Science Foundation of Anhui (1708085QA13), and the key Natural Science Projects of Chizhou University (CZ2018ZRZ06).

References

Linß , Layer-adapted meshes for reaction-convectiondiffusion problem, volume 1985 of lecture Notes in Mathematics. Springer-Verlag, Berlin, 2010.

H.-G.

Roos ,

Stynes and

Tobiska , Numerical methods for singularly perturbed differential eqations, 2nd edn., Springer-Verlag, Berlin, 2008.

J.L.

Gracica ,

O'Riordan and

M.L.

Pickett , A parameter robust second order numerical method for a singularly perturbed two-parameter problem, Appl Numer Math 56 (2006), 962–980.

Teofanov and

Zarin , Superconvergence analysis of a finite element method for a two-parameters singularly perturbed problem, BIT Numer Math 49 (2009), 743-765.

Linß and

Roos , Analysis of a finite-difference scheme for a singularly perturbed problem with two small parameters, J Math Anal Appl 289 (2004), 355–366.

H.-G.

Roos and

Uzelac , The SDFEM for a convection-diffusion problem with two small parameters, Comput Meth Appl Math 3(3) (2003), 443–458.

M.K.

Kadalbajoo and

A.S.

Yadaw , B-spline collocation mehod for a two-parameter singularly pertubed convection-diffusion boundary value problems, Appl Math Comput 201 (2008), 504–513.

Wu ,

Zhang and

Yuan , A robust adaptive method for singularly perturbed convection-diffusion problem with two small parameters, Comput Math Appl 66 (2013), 996–1009.

Linß , A posteriori error estimation for a singularly perturbed problem with two small parameters, Int J Numer Anal Model 7 (2010), 491–506.

10.

Tang and

M.R.

Trummer , Boundary layer resolving pseudospectral methods for singular perturbation problems, SIAM J Sci Comput 17 (1996), 430–438.

11.

Liu and

Tang , Error analysis for a Galerkin-spectral method with coordinate transformation for solving singularly perturbed problems, Appl Numer Math 38 (2001), 315–345.

12.

Chen ,

Wang , A rational spectral collocation method for third-order singularly perturbed problems, J Comput Appl Math 307 (2016), 93–105.

13.

Wang ,

Chen and

Wu , A rational spectral collocation method for solving a class of parameterized singular perturbation problems, J Comput Appl Math 233 (2010), 2652–2660.

14.

L.-B.

Liu ,

Long ,

Huang , et al., Rational spectral collocation and differential evolution algorithms for singularly perturbed problems with an interior layer, J Comput Appl Math 335 (2018), 312–322.

15.

Chen ,

Wang ,

Wu , Rational spectral collocation method for a coupled system of singularly perturbed boundary value problems, J Comput Math 29 (2011), 458–473.

16.

L.-B.

Liu ,

Long ,

Ouyang , et al., Numerical solution of a singularly perturbed problem with Robin boundary conditions using particle swarm optimization algorithm, J Intell Fuzzy Syst 33 (2017), 1785–1795.

17.

M.F.

Fatech ,

Zameer ,

N.M.

Mirze , et al., Biologically inspired computing framework for solving two-point boundary value problems using differential evolution, Neural Comput Appl 28 (2016), 1–15.

18.

Luo ,

L.-B.

Liu ,

Ouyang , et al., B-spline collocation and self-adapting differential evolution (jDE) algorithm for a singularly perturbed convection-diffusion problem, Soft Comput (2017). https://doi.org/10.1007/s00500-017-2523-9.

19.

Abderazek ,

Ferhat and

Atanasovska , A differential evolution algorithm for tooth profile optimization with respect to balancing specific sliding coefficients of involute cylindrical spur and helical gears, Adv Mech Eng 7 (2015), 1–11.

20.

Liu and

Lampinen , A fuzzy adaptive differential evolution algorithm, Soft Computing 9 (2005), 448–462.

21.

Das and

Konar , Automatic image pixel clustering with an improve differential evolution, Appl Soft Comput 9 (2009), 226–236.

22.

J.J.

Guedes ,

M.F.

Castoldi ,

Goedtel , et al., Parameters estimation of three-phase induction motors using differential evolution, Electric Power Systems Research 154 (2018), 204–212.

23.

Lai ,

Y.D.

Ji ,

S.C.

Zhong , et al., Sequential parameter identification of fractional-order duffing system based on differential evolution algorithm, Mathematical Problems in Engineering 2017 (2017), 1–13.

24.

Noman and

Iba , Accelerating differential evolution using an adaptive local search, IEEE T Evolut Comput 12 (2008), 107–125.

25.

X. S.

Yang , Flower pollination algorithm for global optimization. in: Unconventional Computation and Natural Computation, Lecture Notes in Computer Science 7445 (2012), 240–249.

26.

Storn and

Price , Differential evolution: A simple and efficient adaptive scheme for global optimization over continuous spaces, J Global Optim 11 (1997), 341–359.

27.

Kennedy and

R.C.

Eberhart , Particle swarm optimization, in: Proceeding of the IEEE international conference on neural networks (Perth, Australia), IEEE Service center, Piscataway, NJ, 1995, pp. 1942–1948.

28.

Karaboga ,

Basturk , A powerful and efficient algorithm for numerical function optimization: Artificial bee colony (ABC) algorithm. J Global Optim 39 (2007), 459–471.

29.

Hansen ,

S.D.

Muller and

Koumoutsakos , Reducing the time complexity of the derandomized evolution strategy with covariance matrix adaptation (CMA-ES), Evolutionary Computation 11 (2014), 1–18.

A dual mutation differential evolution algorithm for singularly perturbed problems with two small parameters

Abstract

Keywords

1. Introduction

2.1 The barycentric form of rational interpolation

4.1 Test example and error evaluation approach

4.3 Experimental results and analyses

Table 1 Numerical results by using different algorithms (N = 8)

Footnotes

Acknowledgments

References

Table 1
Numerical results by using different algorithms (N = 8)