    # Solving coupled Lane-Emden equations by Green's function and decomposition technique

In this paper, the Green's function and decomposition technique is proposed for solving the coupled Lane-Emden equations. This approach depends on constructing Green's function before establishing the recursive scheme for the series solution. Unlike, standard Adomian decomposition method, the present method avoids solving a sequence of transcendental equations for the undetermined coefficients. Convergence and error estimation is provided. Three examples of coupled Lane-Emden equations are considered to demonstrate the accuracy of the current algorithm.

## Authors

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

This paper aims to extend the application of the Adomian decomposition method with Green’s function [1, 2] for solving the following coupled Lane-Emden boundary value problems

 ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩y′′1(x)+α1xy′1(x)=f1(x,y1(x),y2(x)),    x∈(0,1)y′′2(x)+α2xy′2(x)=f2(x,y1(x),y2(x)),y′1(0)=0,  y′2(0)=0,a1y1(1)+b1y′1(1)=c1,  a2y2(1)+b2y′2(1)=c2, (1.1)

where

are real constants. In recent years, singular boundary value problems for ordinary differential equations have been studied extensively

[3, 4, 5, 6, 7, 8, 9, 10, 11, 1, 2, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21] and references therein. However, we find only the following results on coupled Lane-Emden equations.

Recently, in [22, 23, 24] authors studied (1.1) with boundary conditions and that relates the concentration of the carbon substrate and the concentration of oxygen. In [25, 26], authors considered the coupled Lane–Emden equations (1.1) with boundary conditions , and occurs in catalytic diffusion reactions. In [26, 24], the Adomian decomposition method was applied to obtain a convergent analytic approximate solution of (1.1) with . Later, in , the variational iteration method was applied to obtain approximations to solutions of (1.1) for shape factors . In , the Sinc-collocation method was used to obtain the solution of (1.1). In  authors used the reproducing kernel Hilbert space method for solving to obtain the solution of (1.1).

Recently, many researchers [30, 31, 32, 33, 34, 35, 36] have applied the Adomian decomposition method to deal with many different scientific models. According to the Adomian decomposition method we rewrite (1.1) in a operator form as

 {L1y1(x)=f1(x,y1(x),y2(x)),    x∈(0,1)L2y2(x)=f2(x,y1(x),y2(x)), (2.1)

where and are differential operators and their inverse integral operators are defined as

 ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩L−11[⋅]=x∫0x−α1x∫0xα1[⋅] dxdx,L−12[⋅]=x∫0x−α2x∫0xα2[⋅] dxdx. (2.2)

Operating on (2.1) and using , we get

 {y1(x)=y1(0)+L−11[f1(x,y1(x),y2(x))],y2(x)=y2(0)+L−12[f2(x,y1(x),y2(x))]. (2.3)

According to the ADM, we decompose and as

 ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩y1(x)=∞∑j=0y1j(x),    f1(x,y1(x),y2(x))=∞∑j=0A1j,y2(x)=∞∑j=0y2j(x),    f2(x,y1(x),y2(x))=∞∑j=0A2j, (2.4)

where are Adomian’s polynomials  are given

 Ain=1n!dndλn[fi(x,∞∑j=0y1j λj,  ∞∑j=0y2j λj)]λ=0,   i=1,2. (2.5)

Substituting (2.4) into (2.3), we get

 ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩∞∑j=0y1j(x)=y1(0)+L−11[∞∑j=0A1j],∞∑j=0y2j(x)=y2(0)+L−12[∞∑j=0A2j]. (2.6)

Upon comparing both sides of (2.6), we have

 ⎧⎪ ⎪⎨⎪ ⎪⎩y10(x)=δ1,  y20(x)=δ2,y1j(x,δ1,δ2)=L−11[A1,j−1],y2j(x,δ1,δ2)=L−12[A2,j−1],   j=1,2,3.... (2.7)

where are unknown constants to be determined. The -term series solutions are given as

 ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩ϕ1n(x,δ1,δ2)=n∑j=0y1j(x,δ1,δ2),ϕ2n(x,δ1,δ2)=n∑j=0y2j(x,δ1,δ2). (2.8)

The unknown constants may be obtained by imposing boundary condition at on , that leads to

 {a1ϕ1n(1,δ1,δ2)+b1ϕ′1n(1,δ1,δ2)−c1=0,a2ϕ2n(1,δ1,δ2)+b2ϕ′2n(1,δ1,δ2)−c2=0. (2.9)

Solving above transcendental equations for require additional computational work, and may not be uniquely determined.

To avoid solving the above sequence of difficult transcendental equations, the Adomian decomposition method with Green’s function was introduced in [1, 2]. This technique relies on constructing Green’s function before establishing the recursive scheme for the solution components. Unlike the standard Adomian decomposition method, this avoids solving a sequence of transcendental equations for the undetermined coefficients.

## 3 Green’s function and decomposition technique

In this section, we extend the application of the Adomian decomposition method with Green’s function [2, 12, 37, 38, 39], where we transformed the singular boundary value problem into the integral equation before establishing the recursive scheme for the approximate solution. To apply this technique to coupled Lane-Emden boundary value problems (1.1), we first consider the equivalent integral form of coupled Lane-Emden equation (1.1) as

 ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩y1(x)=c1a1+1∫0 G1(x,s) sα1 f1(s,y1(s),y2(s))ds,y2(x)=c2a2+1∫0 G2(x,s) sα2 f2(s,y1(s),y2(s))ds, (3.1)

where are given by

 Gi(x,s)={lns,x≤s,      αi=1,   i=1,2,lnx,s≤x, (3.2)

and

 Gi(x,s)=⎧⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪⎩s1−αi−11−αi,x≤s,    αi>1,    i=1,2,x1−αi−11−αi,s≤x. (3.3)

Substituting the series (2.4) into (3.1), we obtain

 ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩∞∑j=0y1j(x)=c1a1+1∫0 G1(x,s) sα1 ∞∑j=0A1j ds,∞∑j=0y2j(x)=c2a2+1∫0 G2(x,s) sα2 ∞∑j=0A2j ds. (3.4)

Comparing components from both sides of (3.4) we have the following recursive scheme

 y10(x)=c1a1,   y20(x)=c2a2,y1j(x)=1∫0 G1(x,s) sα1 A1,j−1 ds,y2j(x)=1∫0 G2(x,s) sα2 A2,j−1 ds.⎫⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎬⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎭ (3.5)

Then, we obtain the approximate series solutions as

 ψ1n(x)=n∑j=0y1j(x),ψ2n(x)=n∑j=0y2j(x).⎫⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎬⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎭ (3.6)

Unlike ADM or MADM, the proposed recursive schemes (3.6) do not require any computation of unknown constants.

## 4 Convergence and error analysis

Let be a Banach space with norm

 ∥y∥=max{∥y1∥,∥y2∥},       y∈E, (4.1)

where, and .

From (3.5) and (3.6), we have

 Ψn =n∑j=0yj(x)=ca+n∑j=1[1∫0G(x,s)sαAj−1ds]=ca+1∫0G(x,s)sαn∑j=1Aj−1ds, (4.2)

where

 y =(y1y2),  yj=(y1jy2j),G(x,s)=(G1(x,s)G2(x,s)),  Aj=(A1jA2j) f =(f1f2),  sα=(sα1sα2),  ca=(c1a1c2a2),  Ψn=(ψ1nψ2n).
###### Definition 1.

The function satisfy Lipschitz condition as

 |f(x,y1,y2)−f(x,y∗1,y∗2)|≤2∑j=1lj|yj−y∗j|,  ∀(x,y1,y2), (x,y∗1,y∗2)∈D, (4.3)

where , and , are Lipschitz constants.

###### Theorem 4.1.

Suppose that the nonlinear function satisfy Lipschitz condition (4.3), then the series solution defined by (3.6) is convergent whenever .

###### Proof.

Define

 Ψ0=x0, Ψ1=x0+x1,…,Ψn=n∑k=0yk.

For and using (4.2), we have

 ∥Ψn−Ψm∥ =maxx∈I∣∣∣1∫0G(x,s)s% α(n∑j=1Aj−1−m∑j=1Aj−1)ds∣∣∣.

Using from (), we have

 ∥Ψn−Ψm∥ ≤maxx∈I∣∣∣1∫0G(x,s)sα[f(s,ψ1,n−1,ψ2,n−1)−f(s,ψ1,m−1,ψ2,m−1)]ds∣∣∣.

Applying the Lipschitz condition, we get

 ∥Ψn−Ψm∥ ≤maxx∈I∣∣∣1∫0G(x,s)sαds∣∣∣2∑i=1limaxx∈I|ψi,n−1−ψi,m−1| ≤maxx∈I∣∣∣1∫0G(x,s)sαds∣∣∣2lmax{∥ψ1,n−1−ψ1,m−1∥,∥ψ2,n−1−ϕ2,m−1∥} ≤2ml∥Ψn−1−Ψm−1∥=δ∥Ψn−1−Φm−1∥,

where

 m=max{max∣∣1∫0G1(x,s)sα1ds∣∣,max∣∣1∫0G2(x,s)sα2ds∣∣}, l=max{l1,l2},  γ=2 ml. (4.4)

Thus, we have

 ∥Ψn−Ψm∥≤γ∥Ψn−1−Ψm−1∥. (4.5)

By taking in (4.5), we see that

 ∥Ψm+1−Ψm∥ ≤γ∥Ψm−Ψm−1∥≤γ2∥Ψm−1−Ψm−2∥≤…≤γm∥Ψ1−Ψ0∥.

For all , with , consider

 ∥Ψn−Ψm∥=∥(Ψn−Ψn−1)+(Ψn−1−Ψn−1)+⋯+(Ψm+1−Ψm)∥ ≤∥Ψn−Ψn−1∥+∥Ψn−1−Ψn−2∥+⋯+∥Ψm+1−Ψm∥ ≤[γn−1+γn−2+⋯+γm] ∥Ψ1−Ψ0∥ =γm[1+γ+γ2+⋯+γn−m−1] ∥%y1∥ =γm(1−γn−m1−γ)∥% y1∥.

It follows that as ,

 ∥Ψn−Ψm∥≤ γm1−γ∥y1∥. (4.6)

Letting , we obtain Hence, is a Cauchy sequences in the Banach space . ∎

###### Theorem 4.2.

If the approximate solution converges to , then the maximum absolute truncated error is estimated

 ∥y(x)−Ψm∥≤γm m1−γmaxx∈I∣∣f(x,y10,y20)∣∣. (4.7)
###### Proof.

From (4.6), we have

 ∥Ψn−Ψm∥≤ γm1−γ∥y1∥.

Since as , and the above inequality reduces to

 ∥y(x)−Ψm∥≤γm1−γ∥y1∥. (4.8)

From (3.5), we have , and we find

 ∥y1∥ =maxx∈I∣∣∣1∫0G(x,s)s% αA0ds∣∣∣≤|mmaxx∈I∣∣f% (x,y10,y20)∣∣. (4.9)

Combining (4.8) and (4.9), we obtain error estimate as

 ∥y(t)−Ψm∥≤γm m1−γmaxx∈I∣∣f(x,y10,y20)∣∣. (4.10)

which completes the proof. ∎

## 5 Numerical Results

In this section, we consider three coupled Lane-Emden type boundary value problems to examine the accuracy of the present method. Since the exact solution of the problems is not known, we examine the accuracy and applicability of the present method by the absolute residual error

 (5.1)

where is the absolute residual error and is the present approximate solution. The maximum residual errors are defined as

 ⎧⎪⎨⎪⎩maxr1n:=maxx∈[0,1]r1n(x),maxr2n:=maxx∈[0,1]r2n(x). (5.2)

### 5.1 The catalytic diffusion reactions problem [25, 26]

###### Example 5.1.

Consider the coupled Lane-Emden equation occurs in catalytic diffusion reactions [25, 26] as

 ⎧⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪⎩y′′1(x)+2xy′1(x)=−k1 y21(x)−k2 y1(x)y2(x),     x∈(0,1),y′′2(x)+2xy′2(x)=−k3 y21(x)−k4 y1(x)y2(x),y′1(0)=0,   y′2(0)=0,   y1(1)=1,  y2(1)=2, (5.3)

where the parameters and are the actual chemical reactions. Here , , and .

In Tables 1 and 3, we list the numerical results of the approximate solution and the absolute error obtained by the proposed method of Example 5.1 for () and (), respectively. We also compare the numerical results of the maximum residual error obtained by the present method and the results obtained by the modified ADM  in Table 2. In Table 4, we list the numerical results of the maximum residual error.

#### 5.1.1 When k1=1,k2=2/5,k3=1/2,k4=1

By applying the proposed scheme (3.5) with the initial guesses , we obtain the 5-terms series solutions as

 {ψ15(x)=0.776218+0.199501x2+0.018823x4+0.005706x6−0.0003741x8+0.000125x10,ψ25(x)=1.68423+0.283069x2+0.0258518x4+0.007133x6−0.0004418x8+0.000152x10.

#### 5.1.2 When k1=k2=k3=k4=1/2

On applying the proposed scheme (3.5) with the initial guesses , we obtain the 5-terms series solutions as

 {ψ15(x)=0.80364+0.17704x2+0.017049x4+0.00223x6−0.00001x8+0.0000316x10,ψ25(x)=1.80365+0.17704x2+0.017049x4+0.00223x6−0.00001x8+0.0000316x10.

### 5.2 The concentration of the carbon substrate and the concentration of oxygen problem 

###### Example 5.2.

Consider the coupled Lane-Emden equations, which was used to study the concentration of the carbon substrate and the concentration of oxygen, as

 ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩y′′1(x)+α1xy′1(x)=−b+a y1(x) y2(x)(l1+y1)(m1+y2)+c y1(x) y2(x)(l2+y1(x))(m2+y2(x)),    x∈(0,1),y′′2(x)+α2xy′2(x)=d y1(x) y2(x)(l1+y1(x))(m1+y2(x))+e y1(x) y2(x)(l2+y1(x))(m2+y2(x)),y′1(0)=0,   y1(1)=1,  y′2(0)=0,   y2(1)=1, (5.4)

where the parameters are fixed as given in [27, 24]. Here, .

In Table 5, 6 and 7, we list the numerical results of the approximate solution and the absolute error obtained by the proposed method of Example 5.2 for ), and , respectively. We also compare the numerical results of the maximum residual error obtained by the present method and the results obtained by the modified ADM  in Table 8 for .

#### 5.2.1 When the shape factors α1=α2=1

By applying the proposed scheme (3.5) with the initial guesses , we obtain the 4-terms series solutions as

 {ψ14(x)=2.02484−1.02488x2+0.00006968x4−0.0000308x6+8.57×10−6x8,ψ24(x)=1.0375−0.0374966x2+2.049×10−6x4−9.06×10−7x6+2.52×10−7x8.