Comments on “Homotopy Perturbation Method for Solving Reaction-Diffusion Equations”

Abstract

The paper entitled “Homotopy perturbation method for solving reaction diffussion equation” contains some mistakes and misinterpretations along with a false conclusion. Applying the homotopy perturbation method (HPM) in an incorrect manner, the authors have drawn the false conclusion that this approach is efficient for reaction-diffusion type of equation. We show that HPM in the proposed form is not efficient in most cases, and hence, we will introduce the correct form of HPM.

1. Introduction

In [1] Wang et al. proposed a new approach to apply the HPM for solving the reaction-diffusion problem. We can easily identify in [1] some mistakes, which can mislead the readers, as we show in Section 3 of this paper. The main difficulties in the application of HPM in proposed form are that the corresponding approximate solution is a divergent series on [0,1] and the method of weighted residuals applied in an incorrect form.

We consider a reaction-diffusion process governed by the equation𝑢^′′+𝑢^𝑛=0,0<𝑥<𝐿, (1.1) with boundary conditions𝑢(0)=𝑢(𝐿)=0, (1.2) where 𝑢(𝑥) represents the steady-state temperature for the corresponding reaction-diffusion equation with the reaction term 𝑢^𝑛,𝑛 is the power of the reaction term (heat source), and 𝐿 is the length of the sample (heat conductor). The physical interpretation of (1.1) was given in [2]. For 𝑛=3, the problem (1.1)-(1.2) has a unique positive solution (e.g., as can be seen by phase-plane analysis) which is given in terms of elliptic functions. For an arbitrary real number 𝑛≠3, the general reaction-diffusion process is more complicated. The problem (1.1)-(1.2) was studied by Lesnic using Adomian decomposition method [2] and by Mo [3] using variational method. In this paper, we analyze the HPM used in [1] and demonstrate the correct form of HPM.

2. Homotopy Perturbation Method

In [1] the authors used the HPM in standard form:1−𝑝𝑢^′′−𝑢₀^′′+𝑝𝑢^′′+𝑢^𝑛=0, (2.1) with initial approximation𝑢₀=𝑎𝑥(1−𝑥), (2.2) where 𝑎 is an unknown constant. The HPM supposes that the solution of (1.1)-(1.2) has the form𝑢=𝑢₀+𝑝𝑢₁+𝑝²𝑢₂+⋯. (2.3) And then equating the terms with identical powers of 𝑝, the components 𝑢₁,𝑢₂,… of 𝑢 can be calculated easily. This standard process for 𝑛=2 with 𝐿=1 gives an approximate solution for 𝑢 [1]:𝑢=𝑎𝑥(1−𝑥)+𝑎𝑥²−𝑎²1𝑥30⁶−1𝑥10⁵+1𝑥12⁴−1𝑎−𝑎60²𝑥. (2.4) To find 𝑎 the authors in [1] then used the method of weighted residuals. Substituting (2.4) into (1.1) results in the following residual:𝑅(𝑥,𝑎)=𝑢^′′(𝑥)+𝑢^𝑛(𝑥). (2.5) Then they set 𝑅(1/3,𝑎)=0 and found 𝑎=45.4205 (why 1/3?).

It is not clear why the number 1/3 (but not 1/2) plays some special role for determining 𝑎. Actually, the problem is that the solution corresponding to 1/2 in fact has no relationship with the exact solution; that is error is large enough comparatively with the case 1/3. But it is not difficult to show that the error (in case of 𝑛=2) for 1/3 is also large enough. This solution is not even so good around 1/3. Tables 1 and 2 demonstrate the errors in 𝑢^′′+𝑢^𝑛=0 for some values of 𝑥. If we use one more term (component) in the HPM in [1] applying the method of weighted residuals (for 𝑥=1/3 and 𝑛=2), we obtain 𝑎=23.949 and an approximate solution which in fact has no relationship with normal solution: for example, the “error” for 𝑥=0.9 is 43.780. That is the result received in [1] cannot be improved by taking more terms in approximation.

Table 1 

Errors in approximate solution for Example 3.1.

Table 2 

Errors in approximate solution for Example 3.2.

The failure of the method related also with the convergence rate of the series solution. It is not difficult to show that the coefficients in the highest degree terms of the second derivative are large enough and therefore go to infinity for 𝑥 close to 1.

To overcome these difficulties we offer the next two versions of the HPM.(1) Since the solution must be symmetric with respect to the line 𝑥=𝐿/2, try to represent the solution as a series of ∑(𝑥−𝐿/2)^𝑘. Since the solution reaches the maximum at 𝐿/2 take 𝑢^′(𝐿/2)=0. We construct a homotopy as in (2.1) but with initial approximation: 𝑢₀=𝑏, (2.6) where 𝑏 is an unknown constant to be further determined. Instead of initial-boundary conditions we consider initial conditions: 𝑢_𝑘𝐿2=0,𝑢^′_𝑘𝐿2=0, (2.7) for each component 𝑢_𝑘,𝑘=1,2,….(2) To get rapidly convergent series solution we offer to take an initial approximation as 𝑢₀=𝑎sin(𝜋𝑥/𝐿) where 𝑎 is an unknown constant to be further determined (by taking 𝑥=𝐿/2 in 𝑢^′′+𝑢^𝑛=0, where 𝑢=𝑢_approx).

3. Applications

Example 3.1. Consider the problem 𝑢^′′+𝑢²𝑢(0)(1)=0,=𝑢=0 (3.1) and apply version 1 of HPM. Substituting 𝑢=𝑢₀+𝑝𝑢₁+𝑝²𝑢₂+𝑝³𝑢₃+⋯ into (2.1) and equating the terms with the identical powers of 𝑝, we obtain 𝑢₀^′′+𝑝𝑢₁^′′+𝑝²𝑢₂^′′+𝑝³𝑢₃^′′+⋯−𝑣₀^′′+𝑝𝑣₀^′′+𝑢₀+𝑝𝑢₁+𝑝²𝑢₂+𝑝³𝑢₃+𝑝⁴𝑢₄+𝑝⁵𝑢₅²𝑢=0,₀𝑢=𝑏,_1′′=−𝑢₀^′′−𝑢²₀=−(𝑏)²,𝑢₁(0.5)=𝑢^′₁()𝑢0.5=0,₁=−0.5𝑏²(𝑥−0.5)²,𝑢₂^′′=−2𝑢₀𝑢₁=−2(𝑏)−0.5𝑏²(𝑥−0.5)²,𝑢₂(0.5)=𝑢^′₂()𝑢0.5=0,₂=1𝑏12³(𝑥−0.5)⁴,𝑢₃^′′=−2𝑢₀𝑢₂−𝑢²₁=−2(𝑏)8.3333×10⁻²𝑏³(𝑥−0.5)⁴−−0.5𝑏²(𝑥−0.5)²²,𝑢₃=−1.3889×10⁻²𝑏⁴(𝑥−0.5)⁶,𝑢₄^′′=−2𝑢₀𝑢₃−2𝑢₁𝑢₂,𝑢₄=0.11111𝑏56⁵(𝑥−0.5)⁸,𝑢₅^′′=−2𝑢₀𝑢₄−2𝑢₁𝑢₃−𝑢²₂,𝑢₅=−2.7558×10⁻⁴𝑏⁶(𝑥−0.5)¹⁰,𝑢₆^′′=−2𝑢₀𝑢₅−2𝑢₁𝑢₄−2𝑢₂𝑢₃,𝑢₆=3.6743×10⁻⁵𝑏⁷(𝑥−0.5)¹²,𝑢₇^′′=−2𝑢₀𝑢₆−2𝑢₁𝑢₅−2𝑢₂𝑢₄−𝑢²₃𝑢₇=−4.7948×10⁻⁶𝑏⁸(𝑥−0.5)¹⁴,𝑢₈^′′=−2𝑢₀𝑢₇−2𝑢₁𝑢₆−2𝑢₂𝑢₅−2𝑢₃𝑢²₄,𝑢₈=6.1407×10⁻⁷𝑏⁹(𝑥−0.5)¹⁶,𝑢₉^′′=−2𝑢₀𝑢₈−2𝑢₁𝑢₇−2𝑢₂𝑢₆−2𝑢₃𝑢₅−𝑢₄²,𝑢₉=−7.7577×10⁻⁸𝑏¹⁰(𝑥−0.5)¹⁸. (3.2) Let us take 𝑢_approx=𝑢₀+𝑢₁+⋯+𝑢₉ and find 𝑏 from 𝑢_approx(0)=0—some numeric method can be used for finding 𝑏 (we used standard LaTex computation): 𝑏=11.79. Thus we get an approximate solution: 𝑢=(11.79)−0.5(11.79)²(𝑥−0.5)²+1()1211.79³(𝑥−0.5)⁴−1.3889×10⁻²(11.79)⁴(𝑥−0.5)⁶+1.9841×10⁻³(11.79)⁵(𝑥−0.5)⁸−2.7558×10⁻⁴(11.79)⁶(𝑥−0.5)¹⁰+3.6743×10⁻⁵(11.79)⁷(𝑥−0.5)¹²−4.7948×10⁻⁶(11.79)⁸(𝑥−0.5)¹⁴+6.1407×10⁻⁷(11.79)⁹(𝑥−0.5)¹⁶−7.7577×10⁻⁸(11.79)¹⁰(𝑥−0.5)¹⁸. (3.3) Now let us apply version 2 of HPM. We take 𝑢₀=𝑎sin𝜋𝑥, and equating the terms with the identical powers of 𝑝, we obtain 𝑢₀(),𝑢=𝑎sin𝜋𝑥_1′′=−𝑢₀^′′−𝑢²₀=𝜋²𝑎sin(𝜋𝑥)−𝑎²sin²(𝜋𝑥),𝑢₁(0)=𝑢^′₁(1)𝑢=0,₁1=−4𝑎²𝑥²+14𝑎²1𝑥−8𝜋²𝑎²cos²1𝜋𝑥+8𝜋²𝑎²sin²1𝜋𝑥+8𝜋²𝑎²𝑢−𝑎sin𝜋𝑥,_2′′=−2𝑢₀𝑢₁((−1=−2𝑎sin𝜋𝑥))4𝑎²𝑥²+14𝑎²𝑎𝑥−²cos²𝜋𝑥8𝜋²+𝑎²sin²𝜋𝑥8𝜋²+𝑎²8𝜋²,𝑢−𝑎sin𝜋𝑥₂1=−2𝜋²𝑎³𝑥²2sin𝜋𝑥−𝜋³𝑎³1𝑥cos𝜋𝑥+2𝜋²𝑎³1𝑥sin𝜋𝑥−24𝜋⁴𝑎³cos²𝑎𝜋𝑥sin𝜋𝑥+³cos𝜋𝑥𝜋³+𝑎³sin³𝜋𝑥72𝜋⁴+278𝜋⁴𝑎³𝑎sin𝜋𝑥+²𝑥²2+𝑎²cos²𝜋𝑥4𝜋²−𝑎²sin²𝜋𝑥4𝜋²−𝑎²𝑥2−𝑎²(𝜋+4𝑎)4𝜋³,𝑢₃^′′=−2𝑢₀𝑢₂−𝑢²₁,𝑢₃=1𝑎160⁴𝑥⁵−1𝑎192⁴𝑥⁴−14𝑎²𝑥²−1𝑎480⁴𝑥⁶−3𝜋³𝑎³cos𝜋𝑥−818𝜋⁴𝑎³2sin𝜋𝑥−𝜋⁵𝑎⁴−3sin𝜋𝑥32𝜋²𝑎⁴𝑥³+364𝜋²𝑎⁴𝑥⁴−435256𝜋⁴𝑎⁴𝑥²−18𝜋²𝑎²cos²𝜋𝑥−451288𝜋⁶𝑎⁴cos²+𝜋𝑥2518432𝜋⁶𝑎⁴cos⁴1𝜋𝑥+8𝜋²𝑎²sin²1𝜋𝑥−24𝜋⁴𝑎³sin³𝜋𝑥+451288𝜋⁶𝑎⁴sin²+𝜋𝑥2518432𝜋⁶𝑎⁴sin⁴1𝜋𝑥+8𝜋⁴𝑎³cos²9𝜋𝑥sin𝜋𝑥+64𝜋⁴𝑎⁴𝑥²cos²9𝜋𝑥−64𝜋⁴𝑎⁴𝑥²sin²+6𝜋𝑥𝜋³𝑎³3𝑥cos𝜋𝑥−2𝜋²𝑎³𝑥sin𝜋𝑥−253072𝜋⁶𝑎⁴cos²𝜋𝑥sin²𝜋𝑥+2532𝜋⁵𝑎⁴−9cos𝜋𝑥sin𝜋𝑥64𝜋⁴𝑎⁴𝑥cos²3𝜋𝑥+2𝜋²𝑎³𝑥²9sin𝜋𝑥+64𝜋⁴𝑎⁴𝑥sin²−𝜋𝑥2516𝜋⁵𝑎⁴1𝑥cos𝜋𝑥sin𝜋𝑥+3840𝜋⁴𝑎²180𝜋²𝑎²+4𝜋⁴𝑎²+960𝜋⁴+6525𝑎²𝑥+96136144𝜋⁶𝑎⁴+3𝜋³𝑎³+18𝜋²𝑎². (3.4) (We used very easy LaTex tools for finding integrals.) Now we calculate the number 𝑎 from the relationship 𝑢^′′+𝑢²=0, (3.5) where 𝑢≈𝑢₀+𝑢₁+𝑢₂+𝑢₃ and 𝑥=0.5. We obtain 𝑎=11.569. Table 1 shows the errors in the solutions and comparison with the errors in [1].

Example 3.2. Consider the problem 𝑢^′′+𝑢³𝑢(0)(1)=0,=𝑢=0. (3.6) Substituting 𝑢=𝑢₀+𝑝𝑢₁+𝑝²𝑢₂+𝑝³𝑢₃+⋯ into (2.1), where 𝑢₀=𝑎sin(𝜋𝑥), and equating the terms with the identical powers of 𝑝, we obtain from 𝑢₀^′′+𝑝𝑢₁^′′+𝑝²𝑢₂^′′+𝑝³𝑢₃^′′+⋯−𝑣₀^′′+𝑝𝑣₀^′′+𝑢₀+𝑝𝑢₁+𝑝²𝑢₂+𝑝³𝑢₃+𝑝⁴𝑢₄+𝑝⁵𝑢₅³=0 (3.7) that 𝑢₀(),𝑢=𝑎sin𝜋𝑥_1′′=−𝑢₀^′′−𝑢³₀=𝜋²𝑎sin(𝜋𝑥)−𝑎³sin³(𝜋𝑥),𝑢₁(0)=𝑢^′₁(1)𝑢=0,₁1=−112𝜋𝑎3𝜋³9sin3𝜋𝑥−𝜋𝑎³,𝑢sin𝜋𝑥+12𝜋𝑎sin𝜋𝑥_2′′=−3𝑢²₀𝑢₁=−3(𝑎sin(𝜋𝑥))²−𝑎³cos²𝜋𝑥sin𝜋𝑥12𝜋²+𝑎³sin³𝜋𝑥36𝜋²+3𝑎³sin𝜋𝑥4𝜋²,𝑢−𝑎sin𝜋𝑥₂=11440𝜋³2460𝜋𝑎⁵sin𝜋𝑥−290𝑎3𝜋⁵sin3𝜋𝑥+6𝑎⁵sin5𝜋𝑥5𝜋−3240𝜋𝑎³sin𝜋𝑥+120𝜋𝑎³,𝑢sin3𝜋𝑥_3′′=−3𝑢²₀𝑢₂−3𝑢₀𝑢²₁,𝑢₃1=−756000𝜋⁵63000𝜋³𝑎³sin3𝜋𝑥−3925250𝜋𝑎⁷sin𝜋𝑥−1701000𝜋³𝑎³+sin𝜋𝑥470470𝑎3𝜋⁷sin3𝜋𝑥−12614𝑎5𝜋⁷sin5𝜋𝑥+130𝑎7𝜋⁷sin7𝜋𝑥+6457500𝜋𝑎⁵sin𝜋𝑥−253750𝜋𝑎⁵sin3𝜋𝑥+3150𝜋𝑎⁵sin5𝜋𝑥. (3.8) To find 𝑎 we take 𝑢=𝑢₀+𝑢₁+𝑢₂+𝑢₃ in 𝑢^′′+𝑢³=0 and replace 𝑥=0.5. Thus we obtain 𝑎=3.6033. Table 2 demonstrates the errors in the solution and the comparison with the errors in [1].

4. Conclusions

We proved that Wang et al. who intended to develop in [1] a new approach to apply the HPM for solving reaction-diffusion equations, included in their paper some mistakes and misinterpretations. They used the method of weighted residuals in an incorrect form, and initial approximation is not so good for the convergence of the series solution. We introduced two versions of the HPM, which are very efficient for initial-boundary type of problems.

The main strength of HPM is in fact its fast convergence, but the form of initial approximation and initial-boundary conditions is very important for the efficiency of the HPM.