1. Introduction

Recently, Merkin et al. in ^[7] considered the following reaction-diffusion traveling waves system in region $I$ as follows: for quadratic autocatalytic reaction

$A+B \rightarrow 2 B (rate~k_{1} a b),$

(1.1)

together with a linear decay step

$B \rightarrow C (rate~k_{2} b),$

(1.2)

for cubic autocatalytic reaction

$A+2B \rightarrow 3 B (rate~k_{3} a b^{2})$

(1.3)

together with a linear decay step

$B \rightarrow C (rate~k_{4} b),$

(1.4)

where $a$ and $b$ are concentrations of reactant $A$ and auto-catalyst $B$, $k_{i}(i = 1, 4)$ are the rate constants and $C$ is some inert product of reaction. On the region II we assume that only the (1.1) and (1.3) are taking place for quadratic autocatalytic reaction and cubic autocatalytic reaction respectively. Here, we consider the following system for the dimensionless concentrations $(\alpha_{1}, \beta_{1})$ and $(\alpha_{2}, \beta_{2})$ in region $I$ and $II$ of species $A$ and $B$, respectively with $x>0$ and $t>0$:

$\frac{\partial\alpha_{1}}{\partial t} = \frac{\partial^{2}\alpha_{1}}{\partial x^{2}}-\alpha_{1} \beta_{1}^{2},$

(1.5)

$\frac{\partial\beta_{1}}{\partial t} = \frac{\partial^{2}\beta_{1}}{\partial x^{2}}+\alpha_{1} \beta_{1}^{2}-k \beta_{1}+ \gamma(\beta_{2}-\beta_{1}),$

(1.6)

$\frac{\partial\alpha_{2}}{\partial t} = \frac{\partial^{2}\alpha_{2}}{\partial x^{2}}-\alpha_{2} \beta_{2}^{2},$

(1.7)

$\frac{\partial\beta_{2}}{\partial t} = \frac{\partial^{2}\beta_{2}}{\partial x^{2}}+\alpha_{2} \beta_{2}^{2}+ \gamma(\beta_{1}-\beta_{2}),$

(1.8)

with the boundary conditions

$\alpha_{i}(0, t) = \alpha_{i}(L, t) = 1, ~~\beta_{i}(0, t) = \beta_{i}(L, t) = 0,$

(1.9)

where $k$ and $\gamma$ are the strength of the auto-catalyst decay and the coupling between the two regions respectively.

The present paper is organized as follows: In section-2, we described the idea of the standard HAM. Section-3, is devoted to the application of HAM to CIACS and Section-4, devoted to the numerical results. In the last section, we summarized the result in the conclusion.

---

2. Basic idea of HAM

In recent years, many authors presented homotopy analysis method and its application for differential equations in many ways (see, for example, ^{[6,13,14,26,27,28,29]} and also see for recent results ^{[2,12,20,21,22,23,24,25]}). After motivation with above mentioned works here we consider the following nonlinear differential equation:

$\mathcal{N}[y(t)] = 0, \, \, \, \, t \geq 0,$

(2.10)

where $\mathcal{N}$ is nonlinear differential operator and $y(t)$ is an unknown function. Liao ^[14] constructed the so-called zeroth-order deformation equation :

$(1-q)\mathcal{L} [\phi (t;q)-y_0(t)] = qhH(t)\mathcal{N}[\phi (t;q)],$

(2.11)

where in the following, $q \in[0, 1]$, $h \neq 0$, $H(t)\neq 0$, $\mathcal{L}$, $\phi (t;q)$ be the embedding parameter, auxiliary parameter, auxiliary function, auxiliary linear operator and, respectively, and $y_0(t)$ be an initial guess for $y(t)$ which satisfies the initial conditions. Clearly, when $q = 0$ and $q = 1$, the following relations hold respectively

$\phi (t;0) = y_{0}(t), ~~~~~~ \phi (t;1) = y(t).$

Expanding $\phi (t;q)$ in Taylor series with respect to $q$, one has

$\phi (t;q) = y_0(t)+\sum\limits_{m = 1}^{\infty}y_{m}(t)q^m, \label{h5}$

(2.12)

where

$y_m(t) = {\frac{1}{m!}\frac{\partial^{m}\phi (t;q)}{\partial q^m}}\biggr\rvert_{q = 0}. \label{h4}$

Let us assume that the $h$, $H(t)$, $y_{0}(t)$ and $\mathcal{L}$ are selected such that the series (2.12) converges at $q = 1$, and one has

$y(t) = y_0(t)+\sum\limits_{m = 1}^{\infty}y_{m}(t).$

(2.13)

We can deduce the governing equation from the zero order deformation equation by defining the vector

$\overrightarrow{y_n} = \{y_0(t), y_1(t), y_2(t), \ldots, y_n(t)\}.$

Differentiating (2.11), $m$-times with respect to $q$, then by choosing $q = 0$ and dividing by $m!$, we get the so-called $m$th-order deformation equation

$\mathcal{L}[y_m(t)-\chi_{m}y_{m-1}(t)] = \hbar H(t)R_m(\overrightarrow{y}_{m-1}(t)), \, \, \, \, m = 1, 2, 3, ...., n,$

(2.14)

where

$R_m(\overrightarrow{y}_{m-1}) = \frac{1}{(m-1)!}\frac{\partial^{m-1}\mathcal{N}[\phi(t;q)]}{\partial q^{m-1}}|_{q = 0},$

(2.15)

and

$\begin{equation*} \chi_m = \left\{ \begin{array}{rcl} 0, \, \, \, \,&m\leq 1 \\ 1 , \, \, \, \,&m>1 \end{array}\right. \label{h10} \end{equation*}$

More detailed analysis of HAM and the modified version of it together with various applications could be found in ^{[4,8,9,10,11,17,18,19]}.

---

3. HAM solution of CIACS

Here, we apply the HAM on CIACS. We take the initial conditions to satisfy the boundary conditions, namely

$\alpha_{i}(x, 0) = 1-\sum\limits_{n = 1}^{\infty} a_{ni} \cos(0.5 (L - 2 x) \lambda_{n}) \sin(\lambda_{n} L/2), (i = 1, 2),$

(3.16)

$\beta_{i}(x, 0) = \sum\limits_{n = 1}^{\infty} b_{ni} \cos(0.5 (L - 2 x) \lambda_{n}) \sin(\lambda_{n} L/2), (i = 1, 2),$

(3.17)

where $\lambda_{n} = \frac{n \pi}{L}.$ The HAM is based on a kind of continuous mapping

$\alpha_{i}(x, t)\rightarrow \phi_{i}(x, t;q), \quad \beta_{i}(x, t)\rightarrow \psi_{i}(x, t;q)$

such that, as the embedding parameter q increases from 0 to 1, $\phi_{i}(x, t;q), \psi_{i}(x, t;q)$ and $i = 1, 2$ varies from the initial approximation to the exact solution.

We define the nonlinear operators

$\begin{eqnarray*} \label{M10} \mathcal{N}_{i}(\phi_{i}(x, t;q)) & = & \phi_{i, t}(x, t;q)-\phi_{i, xx}(x, t;q)+ \phi_{i}(x, t;q) \psi_{i}^{2}(x, t;q), \\ \mathcal{M}_{i}(\psi_{i}(x, t;q)) & = & \psi_{i, t}(x, t;q)-\psi_{i, xx}(x, t;q)+(-2(i-1)k+ik) \psi_{i}(x, t;q)\\ && + (-1)^{i}\gamma (\psi_{1}(x, t;q)-\psi_{2}(x, t;q))-\phi_{i}(x, t;q) \psi_{i}^{2}(x, t;q), \end{eqnarray*}$

Now, we construct a set of equations, using the embedding parameter $q$

$\begin{equation*}\label{M13} (1-q) \mathcal{L}_{i} (\phi_{i}(x, t;q)-\alpha_{i0}(x, t)) = q h H(x, t) \mathcal{N}_{i}(\phi_{i}(x, t;q)), \\(1-q) \mathcal{L}_{i} (\psi_{i}(x, t;q)-\beta_{i0}(x, t)) = q h H(x, t) \mathcal{M}_{i}(\psi_{i}(x, t;q)), \end{equation*}$

with the initial conditions

$\phi_{i}(x, 0;q) = \alpha_{i0}(x, 0), \quad \psi_{i}(x, 0;q) = \beta_{i0}(x, 0), (i = 1, 2)$

Where $h\neq 0$ and $H(x, t) \neq 0$ are the auxiliary parameter and function, respectively. We expand $\phi_{i}(x, t;q)$ and $\psi_{i}(x, t;q)$ in a Taylor series with respect to $q$, and get

$\phi_{i}(x, t;q) = \alpha_{i0}(x, t)+\sum\limits_{m = 1}^{\infty} \alpha_{im}(x, t) q^{m},$

(3.18)

$\psi_{i}(x, t;q) = \beta_{i0}(x, t)+\sum\limits_{m = 1}^{\infty} \beta_{im}(x, t) q^{m},$

(3.19)

where

$\begin{equation*}\label{M18} \alpha_{im}(x, t) = \frac{1}{m!}\frac{\partial^{m}\phi_{i}(x, t;q)}{\partial q^{m}}|_{q = 0}, \\ \beta_{im}(x, t) = \frac{1}{m!}\frac{\partial^{m}\psi_{i}(x, t;q)}{\partial q^{m}}|_{q = 0}. \end{equation*}$

Let $q = 1$ into (3.18)--(3.19), the series become

$\begin{equation*}\label{M20} \alpha_{i}(x, t) = \alpha_{i0}(x, t)+\sum\limits_{m = 1}^{\infty} \alpha_{im}(x, t), \\ \beta_{i}(x, t) = \beta_{i0}(x, t)+\sum\limits_{m = 1}^{\infty} \beta_{im}(x, t). \end{equation*}$

Now, we construct the mth-order deformation equation from (2.14)--(2.15) as follows:

$\begin{equation*}\label{M22} \mathcal{L}_{i} (\alpha_{im}(x, t)-\mathcal{X}_{m} \alpha_{i(m-1)}(x, t)) = h H(x, t) R_{1}((\vec{\alpha}_{i(m-1)}, \vec{\beta}_{i(m-1)})), \\ \mathcal{L}_{i} (\beta_{im}(x, t)-\mathcal{X}_{m} \beta_{i(m-1)}(x, t)) = h H(x, t) R_{2}((\vec{\alpha}_{i(m-1)}, \vec{\beta}_{i(m-1)})), \end{equation*}$

with the initial conditions $\alpha_{im}(x, 0) = 0, \quad \beta_{im}(x, 0) = 0, m>1$ where

$R_{1}((\vec{\alpha}_{i(m-1)}, \vec{\beta}_{i(m-1)})) = \frac{\partial \alpha_{i(m-1)}}{\partial t}-\frac{\partial^{2} \alpha_{i(m-1)}}{\partial x^{2}}+\sum\limits_{r = 0}^{m-1}\sum\limits_{j = 0}^{r}\alpha_{i(m-1-r)} \beta_{i(j)} \beta_{i(r-j)},$

$\begin{eqnarray*} \label{M25} \begin{split} R_{2}((\vec{\alpha}_{i(m-1)}, \vec{\beta}_{i(m-1)})) = &\frac{\partial \beta_{i(m-1)}}{\partial t}-\frac{\partial^{2} \beta_{i(m-1)}}{\partial x^{2}}+(-2(i-1)k+ik)\beta_{i(m-1)} \\ & +(-1)^{i}\gamma (\beta_{1(m-1)}-\beta_{2(m-1)})-\sum\limits_{r = 0}^{m-1}\sum\limits_{j = 0}^{r}\alpha_{i(m-1-r)} \beta_{i(j)} \beta_{i(r-j)}. \end{split} \end{eqnarray*}$

If we take $\mathcal{L}_{i} = \frac{d}{dt}, (i = 1, 2)$ then the right inverse of $\mathcal{L}_{i} = \frac{d}{dt}$ will be $\int_{0}^{t}(.) d\tau$

$\alpha_{im} = \mathcal{X}_{m} \alpha_{i(m-1)}+h \int_{0}^{t} \left(\frac{\partial\alpha_{i(m-1)}}{\partial\tau}-\frac{\partial^{2} \alpha_{i(m-1)}}{\partial x^{2}}+\sum\limits_{r = 0}^{m-1}\sum\limits_{j = 0}^{r}\alpha_{i(m-1-r)} \beta_{i(j)} \beta_{i(r-j)} \right)d\tau,$

(3.20)

$\begin{eqnarray} \label{M27} \begin{split} \beta_{im} = &\mathcal{X}_{m}\beta_{i(m-1)}+h \int_{0}^{t}\left( \frac{\partial\beta_{i(m-1)}}{\partial\tau}-\frac{\partial^{2} \beta_{i(m-1)}}{\partial x^{2}}+(-2(i-1)k+ik)\beta_{i(m-1)}\right)d\tau \\ & +h\int_{0}^{t}\left((-1)^{i}\gamma (\beta_{1(m-1)}-\beta_{2(m-1)})-\sum\limits_{r = 0}^{m-1}\sum\limits_{j = 0}^{r}\alpha_{i(m-1-r)} \beta_{i(j)} \beta_{i(r-j)}\right)d\tau. \end{split} \end{eqnarray}$

(3.21)

Let the initial approximation

$\alpha_{i0}(x, t) = \alpha_{i0}(x, 0), \quad \beta_{i0}(x, t) = \beta_{i0}(x, 0).$

(3.22)

For $m = 1$, we obtain the first approximation as following:

$\alpha_{i1} = h \int_{0}^{t}\left( \frac{\partial\alpha_{i0}}{\partial\tau}-\frac{\partial^{2} \alpha_{i0}}{\partial x^{2}}+\alpha_{i0} \beta^{2}_{0i} \right) d\tau,$

(3.23)

$\beta_{i1} = h \int_{0}^{t} \left(\frac{\partial\beta_{i0}}{\partial\tau}-\frac{\partial^{2} \beta_{i0}}{\partial x^{2}}+(2-i)k\beta_{i0}+(-1)^{i} \gamma (\beta_{10}-\beta_{20})-\alpha_{i0} \beta^{2}_{0i}\right) d\tau.$

(3.24)

---

4. Numerical results

Here, we compute the average residual error and the residual error and investigate the intervals of convergence by the $h$-curves. Finally, we checked the accuracy of the HAM solutions by comparing with another numerical method. The first approximation of $\alpha_{i1}(x, t)$ and $\beta_{i1}(x, t)$ are

$\alpha_{i, 1}(x, t) = \sum\limits_{n = 1}^{\infty} \lambda_{n}^{2} a_{ni} \cos(\delta_{n}) \sin(\lambda_{n} L/2)h t+ \alpha_{i0}(x, t) \beta_{i0}(x, t)^{2}h t,$

(4.25)

$\begin{eqnarray}\label{n1} \beta_{i, 1}(x, t) & = & \sum\limits_{n = 1}^{\infty} \lambda_{n}^{2} b_{ni} \cos(\delta_{n}) \sin(\lambda_{n} L/2) h t+(2-i) k \sum\limits_{n = 1}^{\infty} b_{ni} \cos(\delta_{n}) \sin(\lambda_{n} L/2) h t \nonumber \\ && +(-1)^{i} \gamma \left(\sum\limits_{n = 1}^{\infty} b_{ni} \cos(\delta_{n}) \sin(\lambda_{n} L/2) -\sum\limits_{n = 1}^{\infty} b_{ni} \cos(\delta_{n}) \sin(\lambda_{n} L/2)\right)h t\\ \nonumber &&- \alpha_{i0}(x, t) \beta_{i0}(x, t)^{2}ht, \\ \end{eqnarray}$

(4.26)

$\begin{eqnarray} \alpha_{i0}(x, t) \beta_{i0}(x, t)^{2} & = & \sum\limits_{n = 1}^{\infty} \sum\limits_{m = 1}^{\infty}b_{ni} b_{mi} \cos(\delta_{n}) \sin(\lambda_{n} L/2) \cos(\delta_{m}) \sin(\lambda_{m} L/2) h t \nonumber\\ &&-\sum\limits_{n = 1}^{\infty} \sum\limits_{m = 1}^{\infty} \sum\limits_{r = 1}^{\infty}b_{ni} b_{mi} \cos(\delta_{n}) \sin(\lambda_{n} L/2) b_{mi}\cos(\delta_{m}) \sin(\lambda_{m} L/2) h t \\ \nonumber && \times\mbox{}b_{ri}\cos(\delta_{r}) \sin(\lambda_{r} L/2)h t, \\ \end{eqnarray}$

(2.7)

$\delta_{n} = (0.5 (L-2x)\lambda_{n}), \delta_{m} = (0.5 (L-2x)\lambda_{m}), \delta_{r} = (0.5 (L-2x)\lambda_{r}),$

(4.28)

$\lambda_{m} = \frac{m \pi}{L}, ~\lambda_{r} = \frac{r \pi}{L}.$

(4.29)

And so on, in the same manner the rest of approximations can be obtained using the Mathematica package.

---

4.1. H-curves

To observe the intervals of convergence of the HAM solutions, we plot the $h$-curves of $4$, $5$, $6$ terms of HAM solutions in Figure 1(a)-(d). In Figure 1, we plot $\alpha_{1t}(x, 0)$, $\beta_{1t}(x, 0)$, $\alpha_{2t}(x, 0)$ and $\beta_{2t}(x, 0)$ against $h$ respectively at $k = 0.1, \gamma = 0.2, L = 100, x = 3, a_{n_{1}} = 0.001, a_{n_{2}} = 0.002, b_{n_{1}} = 0.001, b_{n_{2}} = 0.002.$ From these figures, we note that the straight line that parallels the $h$-axis gives the valid region of the convergence ^[14].

---

4.2. Average residual errors

We notice that $h$-curve does not give the best value for the $h$. Therefore, we evaluate its optimal values by the $\min$ of the averaged residual errors ^{[1,3,5,12,15,16,26]}.

$E_{\alpha_{i}}(h) = \frac{1}{N M} \sum\limits_{s = 0}^{N} \sum\limits_{j = 0}^{M} \left[\mathcal{N} \left(\sum\limits_{k = 0}^{m}\alpha_{ik}\left(\frac{100 s}{N}, \frac{30 j}{M}\right)\right)\right]^{2},$

(4.30)

$E_{\beta_{i}}(h) = \frac{1}{N M} \sum\limits_{s = 0}^{N} \sum\limits_{j = 0}^{M} \left[\mathcal{M} \left(\sum\limits_{k = 0}^{m}\beta_{ik}\left(\frac{100 s}{N}, \frac{ 30 j}{M}\right)\right)\right]^{2},$

(4.31)

corresponding to a nonlinear algebraic equations

$\frac{d E_{\alpha_{i}}(h)}{d h} = 0,$

(4.32)

$\frac{d E_{\beta_{i}}(h)}{d h} = 0.$

(4.33)

We represent $E_{\alpha_{i}}(h)$ and $E_{\beta _{i}}(h)$ in Figure 2(a)-(d) and in Tables 1--4. Figure 2 and Tables 1--4 show that the $E_{\alpha_{i}}(h)$ and $E_{\beta _{i}}(h)$ for 2, 3, 4, 5, 6 terms HAM solutions. We set into (4.32)--(4.33) $N = 100$ and $M = 30$ with $k = 0.1, \gamma = 0.2, L = 100, x = 3, a_{n_{1}} = 0.001, a_{n_{2}} = 0.002, b_{n_{1}} = 0.001, b_{n_{2}} = 0.002.$ We use the command FinMinimum of Mathematica to get the optimal values $h$.

---

4.3. Comparison analysis

Now, we compare 5-terms of HAM solutions obtained with a numerical method using the commands with Mathematica 9 for solving the system of partial differential equations numerically. We plot the 5-terms of HAM solutions in Figure 3. Figure 3 shows the comparison of HAM solutions HAM solutions with numerical method for $k = 0.1, \gamma = 0.2, L = 100, x = 3, a_{n_{1}} = 0.001, a_{n_{2}} = 0.002, b_{n_{1}} = 0.001, b_{n_{2}} = 0.002.$ We noted from this figure that the HAM solutions have a good agreement with the numerical method. Figure 4 shows the 3-terms HAM solutions obtained.

---

5. Conclusion

In the present research, the HAM was employed to analytically compute approximate solutions of CIACS. By comparing these approximate solutions with numerical solutions and the averaged residual error were found. We show the convergence region by $h$-curves. The agreement with the numerical solutions are very good. The results show that HAM accurate for solving CIACS. By increasing the number of iterations one can reach any desired accuracy. In this paper, we used Mathematica 9 in all calculations.

---

Acknowledgements

First author thanks Saeid Abbasbandyand, Hossein Jafari and R. V. Gorder for stimulating discussions during the preparation of this article.

---

Conflict of interest

The authors declare that there is no conflict of interests regarding the publication of this paper.

---