# Analytical Solutions of Boundary Values Problem of 2D and 3D Poisson and Biharmonic Equations by Homotopy Decomposition Method

Authors: Abdon Atangana; Adem Kılıçman
Journal: Abstract and Applied Analysis
 (2013)

---

## Abstract

The homotopy decomposition method, a relatively new analytical method, is used to solve the 2D and 3D Poisson equations and biharmonic equations. The method is chosen because it does not require the linearization or assumptions of weak nonlinearity, the solutions are generated in the form of general solution, and it is more realistic compared to the method of simplifying the physical problems. The method does not require any corrected function or any Lagrange multiplier and it avoids repeated terms in the series solutions  compared to the existing decomposition method including the variational iteration method, the Adomian decomposition method, and Homotopy perturbation method. The approximated solutions obtained converge to the exact solution asN tends to infinity.


---

## Body

## 1. Introduction

The numerical solution of Poisson equations and biharmonic equations is an important problem in numerical analysis. A vast arrangement of investigating effort has been published on the development of numerical solution of Poisson equations and biharmonic equations. The finite difference schemes of second and fourth order for the solution of Poisson’s equation in polar coordinates have been derived by Mittal and Gahlaut [1]. A numerical method to interpolate the source terms of Poisson’s equation by using B-spline approximation has been devised by Perrey-Debain and ter Morsche [2]. Sutmann and Steffen [3] proposed compact approximation schemes for the Laplace operator of fourth and sixth order; the schemes are based on a Padé approximation of the Taylor expansion for the discretized Laplace operator. Ge [4] used fourth-order compact difference discretization scheme with unequal mesh sizes in different coordinate directions to solve a 3D Poisson equation on a cubic domain. Gumerov and Duraiswami [5] developed a complete translation theory for the biharmonic equation in three dimensions. Khattar et al. [6] derived a fourth-order finite difference approximation based on arithmetic average discretization for the solution of three-dimensional nonlinear biharmonic partial differential equations on a 19-point compact stencil using coupled approach. Altas et al. [7] used multigrid and preconditioned Krylov iterative methods to solve three-dimensional nonlinear biharmonic partial differential equations. Jeon [8] derived scalar boundary integral equation formulas for both interior and exterior biharmonic equations with the Dirichlet boundary data. A spectral collocation method for numerically solving two-dimensional biharmonic boundary-value problems has been reported in [9]. An indirect radial-basis-function collocation method for numerically solving biharmonic boundary-value problems has been reported in [10]. A high-order boundary integral equation method for the solution of biharmonic equations has been presented in [11]. A Galerkin boundary node method for solving biharmonic problems was developed in [12]. An integral collocation approach based on Chebyshev polynomials for numerically solving biharmonic equations for the case of irregularly shaped domains has been developed by Mai-Duy et al. [13]. A numerical method, based on neural-network-based functions, for solving partial differential equations has been in [14]. Mai-Duy and Tanner [15] presented a collocation method based on a Cartesian grid and a 1D integrated radial basis function scheme for numerically solving partial differential equations in rectangular domains and Haar wavelet presented in [16]. The aim of this paper is to solve these problems via the homotopy decomposition method.


## 2. Method

In this study we follow the method of [17–20]. In order to illustrate the basic idea of this method we consider a general nonlinear nonhomogeneous partial differential equation with initial conditions of the following form
(1)∂mU(x,t)∂tm=L(U(x,t))+N(U(x,t))+f(x,t),∂mU(x,t)∂tm=L(U(x,t))+N(U2)m=1,2,3,…,
subject to the initial conditions
(2)∂iU(x,0)∂ti=fm(x),∂m-1U(x,0)∂tm-1=0,∂iU(x,0)∂ti=0fm(x),i=0,1,2,…,m-2,
where f is a known function, N is the general nonlinear differential operator, and L represents a linear differential operator. The method’s first step here is to apply the inverse operator ∂m/∂tm of on both sides (1) to obtain
(3)U(x,t)=∑k=0m-1tkk!dku(x,0)dtk+∫0t∫0t1⋯∫0tm-1L(U(x,τ))+N(U(x,τ))+f(x,τ)dτ⋯dt.
The multi-integrals in (3) can be transformed to
(4)∫0t∫0t1⋯∫0tm-1L(U(x,τ))+N(U(x,τ))+f(x,τ)dτ⋯dt1=1(m-1)!∫0t(t-τ)m-1L(U(x,τ))+N(U(x,τ))+f(x,τ)dτ.
So that (3) can be reformulated as
(5)U(x,t)=∑k=0m-1tkk!{dku(x,0)dtk}+1(m-1)!∫0t(t-τ)m-1L(U(x,τ))+N(U(x,τ))+f(x,τ)dτ.
Using the homotopy scheme the solution of the previous integral equation is given in a series form as
(6)U(x,t,p)=∑n=0∞pnUn(x,t),U(x,t)=limp→1U(x,t,p)
and the nonlinear term can be decomposed as
(7)NU(r,t)=∑n=1∞pnℋn(U),
where p∈(0,1] is an embedding parameter. ℋn(U) is He’s polynomials [21] that can be generated by
(8)ℋn(U0,…,Un)=1n!∂n∂pn[N(∑j=0npjUj(x,t))],ℋn(U0,…,Un)Uj(x,t)=1n!∂n∂pn55n=0,1,2….
The homotopy decomposition method is obtained by the graceful coupling of decomposition method with He’s polynomials and is given by
(9)∑n=0∞pnUn(x,t)=T(x,t)+p1(m-1)!∫0t(t-τ)m-1×[f(x,τ)+L(∑n=0∞pnUn(x,τ))+∑n=0∞pnℋn(U)]dτ
with
(10)T(x,t)=∑k=0m-1tkk!{dku(x,t)dtk∣t=0}.
Comparing the terms of the same power of p gives the solutions of various orders. The initial guess of the approximation is T(x,t). Some further related results can be seen in [22–25].Lemma 1 (see [17]).
The complexity of the homotopy decomposition method is of orderO(n).Proof.
The number of computations including product, addition, subtraction, and division are as follows.
In step 2U0
: 0 because it is obtained directly from the initial conditionsU1: 3⋮Un: 3.
Now in step 4 the total number of computations is equal to∑j=0nUj(x,t)=3n=O(n).


## 3. Solutions of the Main Problems

Problem 1.
Consider the following equation(11)∂2u∂x2+∂2u∂y2=sin(πx)sin(πy);u(x,y)=0along  the  boundaries,0≤x,y≤1;ux(0,y)=-sin(yπ)2π.
The exact solution of the previous equation is given as
(12)u(x,y)=sin(xπ)sin(πy)-2π2.
In the view of the homotopy decomposition method, (11) can be first transformed to
(13)u(x,y)=u(0,y)-sin(πy)2πx+∫0x(x-τ)[sin(πτ)sin(πy)-uyy(τ,y)],u(x,y,p)=∑n=0∞pnun(x,y).
Following the decomposition techniques, we obtain the following equation
(14)∑n=0∞pnun(x,y)=T(x,y)+p∫0x(x-τ)[[∑n=0∞pnun(x,y)]sin(πτ)sin(πy)∑n=0∞pnun(x,y)-∂2∂y2-∂2∂y2[∑n=0∞pnun(x,y)]].
Comparing the terms of the same power of p leads to
(15)p0:u0(x,y)=-sin(πy)2πx,p1:u1(x,y)=∫0x(x-τ)[sin(πτ)sin(πy)-∂2∂y2[u0]]dτ,p1:u1(x,y)=∫0x(x-τ)u1(x,y)=0alongtheboundaries,p2:u2(x,y)=∫0x(x-τ)[-∂2∂y2[u1]]dτ,p3:u3(x,y)=∫0x(x-τ)[-∂2∂y2[u2]]dτ,pn:un(x,y)=∫0x(x-τ)[-∂2un-1∂y2]dτ,pn:nun(x,y)=0alongtheboundaries.
The following solutions are obtained:
(16)u0(x,y)=-sin(πy)2πx,u1(x,y)=[xπ-πx32×3!]sin(πy)-sin(πτ)sin(πy)π2,u2(x,y)=[-xπ+πx36-π3x5240]sin(πy)+sin(πτ)sin(πy)π2,u3(x,y)=[xπ-πx36+π3x5120-π5x710080]sin(πy)u3(x,y)=-sin(πτ)sin(πy)π2,u4(x,y)=[-xπ+πx36-π3x5120+π5x75040-π7x9725760]sin(πy)+sin(πτ)sin(πy)π2,u5(x,y)=[xπ-πx36+π3x5120-π5x75040+π7x9362880-π9x1179833600]sin(πy)-sin(πτ)sin(πy)π2,u6(x,y)=[-xπ+πx36-π3x5120+π5x75040-π7x9362880+π9x1139916800-π11x1312454041600]sin(πy)+sin(πτ)sin(πy)π2,u7(x,y)=[xπ-πx36+π3x5120-π5x75040+π7x9362880-π9x1139916800=xπ1+π11x136227020800-π13x152615348736000]sin(πy)-sin(πτ)sin(πy)π2.
In the same manner one can obtain the rest of the components. But for eight terms were computed and the asymptotic solution is given by
(17)uN=8(x,y)=[x2π-πx32×3!+π3x52×5!-π5x72×7!+π7x92×9!x2π--π9x112×11!+π11x132×13!-π13x152×15!]sin(πy)-1π2sin(xπ)sin(yπ).
Therefore in general for any N>8 we have
(18)uN=n(x,y)=[12π2∑n=0N(-1)n(xπ)2n+1(2n+1)!]sin(πy)uN=n(x,y)=-1π2sin(xπ)sin(yπ),limN→∞uN(x,y)=12π2sin(πx)sin(πy)limN→∞uN(x,y)-1π2sin(xπ)sin(yπ)limN→∞uN(x,y)=-12π2sin(xπ)sin(yπ).
This is the exact solution of the problem. Figures 1 and 2 show the comparison of the exact solution and the approximated one for N=4. The approximate solution and the exact solution are compared in Figures 1 and 2, respectively.
The numerical errors forN=4 are evaluated in Table 1.Table 1
Evaluation of numerical errors forN=4.
x
Y
u
(
x
,
y
) exact
u
(
x
,
y
)  N=4
Error
0.25
0.25
−0.0253303
−0.0253303
6.27007
·
10
-
11
0.5
−0.0358224
−0.0358224
8.86722
·
10
-
11
0.75
−0.0253303
−0.0253303
6.27007
·
10
-
11
0.95
−0.00560387
−0.00560387
1.38714
·
10
-
11
0.5
0.25
−0.0358224
−0.0358224
1.26904
·
10
-
11
0.5
−0.0506606
−0.0506604
1.79469
·
10
-
7
0.75
−0.0358224
−0.0358223
1.26904
·
10
-
11
0.95
−0.00792506
−0.00792506
2.80752
·
10
-
8
0.75
0.25
−0.0253303
−0.0253195
1.07646
·
10
-
5
0.5
−0.0358224
−0.0358072
1.52235
·
10
-
5
0.75
−0.0253303
−0.0253195
1.07646
·
10
-
5
0.95
−0.00560387
−0.00560148
2.38148
·
10
-
6
0.95
0.25
−0.00560387
−0.00546191
0.000141956
0.5
−0.00792506
−0.00772431
000200756
0.75
−0.00560387
−0.00546191
0.000141956
0.95
−0.00123975
−0.00120835
3.14051
·
10
-
5Figure 1
Exact solution.Figure 2
Approximated solution for the 4 first terms.Problem 2.
Consider 3D Poisson equation:(19)x,,∂2u∂x2+∂2u∂y2+∂2u∂z2=sin(πx)sin(πy)sin(πz),u(x,y,z)=0along  the  boundaries,0≤x,y≤1.
Following the discussion presented earlier we obtain the following set of integral equations:
(20)p0:u0(x,y)=-sin(πy)sin(πz)3πx,p1:u1(x,y)=∫0x(x-τ)[sin(πτ)sin(πy)sin(πz)-∂2u0∂y2]dτ,pn:un(x,y)=∫0x(x-τ)[-∂2un-1∂y2]dτ,un(x,y)=0alongtheboundaries,n≥2.
The following solutions are obtained:
(21)u0(x,y,z)=-sin(πy)sin(πz)3πx,u1(x,y,z)=[1πx-πx39]sin(πy)sin(πz)u1(x,y,z)=-sin(πτ)sin(πy)sin(πz)π2,u2(x,y,z)=[-2πx+πx39-π3x590]sin(πy)sin(πz)u1(x,y,z)=+2sin(πτ)sin(πy)sin(πz)π2,u3(x,y,z)=[4πx-2πx33+π3x530-π5x71890]sin(πy)sin(πz)-4sin(πτ)sin(πy)sin(πz)π2,u4(x,y,z)=[-8πx+4πx33-π3x515+π5x7630]sin(πy)sin(πz)+8sin(πτ)sin(πy)sin(πz)π2,u5(x,y,z)=[16πx-8πx33+2π3x515-π5x7315u5(x,y,z)=+π7x922680-π7x113742200]sin(πy)sin(πz)u5(x,y,z)=-16sin(πτ)sin(πy)sin(πz)π2.
In the same manner one can obtain the rest of the components. But for six terms were computed and the asymptotic solution is given by
(22)u(x,y,z)N=6=[x3π-πx318+π3x5360-π5x715120+π7x91088640-π9x11119750400]sin(πy)sin(πz)-2sin(πτ)sin(πy)sin(πz)π2,u(x,y,z)N=6=13π2[x3π-πx318+π3x5360-π5x715120+π7x91088640-π9x11119750400]sin(πy)sin(πz)-2sin(πτ)sin(πy)sin(πz)3π2,u(x,y,z)N=6=13π2[πx-(πx)33!+(πx)55!-(πx)77!+(πx)99!-(πx)311!]sin(πy)sin(πz)-2sin(πτ)sin(πy)sin(πz)3π2.
Therefore, for any n≥6, the partial sum is given as
(23)uN=n(x,y,z)=13π2[∑k=1N(-1)k(πx)2k+1(2k+1)!]sin(πy)sin(πz)uN=n(x,y,z)=-2sin(πτ)sin(πy)sin(πz)3π2.
Thus
(24)u(x,y,z)=limN→∞uN=n(x,y,z)u(x,y,z)=sin(πτ)sin(πy)sin(πz)3π2u(x,y,z)=-2sin(πτ)sin(πy)sin(πz)3π2u(x,y,z)=-sin(πτ)sin(πy)sin(πz)3π2.
And this is the exact solution to the problem. One can evaluate error committed by choosing the N first terms in the series solutions, in the same manner as in Table 1. The accuracy of the results is estimated by error function
(25)RN(x,y,z)=|uN(x,y,z)-u(x,y,z)|.Problem 3.
Let us consider the following biharmonic equation(26)d4u(x)dx4+4u(x)=0,
for which the exact solution is
(27)u(x)=Exp[1-x]cos[x]cos[1].
The aim of this part is to compare the numerical results obtained via HDM and the method used in [26].
Applying the steps involved in the HDM, we arrive at the following:(28)u0(x)=eSec(1)(1-x+x33),u1(x)=-eSec(1)[x44-x520+x7420],u2(x)=-eSec(1)[-x81120+x910080-x11554400],u3(x)=-eSec(1)[x122217600-x1364864800+x151135134000],u4(x)=-eSec(1)[-x1616144128000+x17274450176000=-eSec(1555)-x1946930980096000],u5(x)=-eSec(1)[x20312873200640000-x216570337213440000=-eSec(1555)+x231662295315000320000].
In the same manner, one can obtain the remaining term by using the following recursive formula:
(29)un+1(x)=-∫0x(x-t)3un(t)dt.
In this paper we consider only the first six terms of the series solution as follows:
(30)uN=6=∑n=05un(x).
To access the accuracy of the method used in paper, we compare in Table 2 the numerical results of the above equation, the solution obtained in [26] with the exact solution.Table 2
Comparison of the HDM and [1] results with the exaction solution for N=6.
x
HDM
Exact
ADM
Err for HDM
Err for ADM
−1.0
7.38906
7.38906
7.38906
6.78
E
-
16
8.88
E
-
16
−0.6
7.56598
7.56598
7.56598
4.76
E
-
12
7.96
E
-
12
−0.2
6.02244
6.02244
6.02244
0.015
E
-
11
1.46
E
-
11
0.2
4.03696
4.03696
4.03696
0.017
E
-
11
1.80
E
-
11
0.6
2.27883
2.27883
2.27883
0.015
E
-
11
1.46
E
-
11
1.0
1.0
1.0
1.0
1.24
E
-
15
2.22
E
-
15Problem 4.
We consider the 2D biharmonic equation(31)∂4u∂x4+2∂4u∂x2∂y2+∂4u∂x4=sin(3πx)sin(3πy),0≤x,y≤1,
subject to the initial conditions:
(32)∂u(x,y)∂x∣(x=0)=sin(3πy)108π3,∂x,xu(0,y)=0,∂x,x,xu(0,y)=-sin(3πy)12π.
In the view of the homotopy decomposition method, the following integral equations are obtained:
(33)p0:u0(x,y)=sin(3πy)108π3x-sin(3πy)12πx3!3,p1:u1(x,y)=∫0x(x-τ)[sin(πτ)sin(πy)-∂2u0∂y2-2∂4u0∂x2∂y2]dτ,pn:un(x,y)=∫0x(x-τ)[-∂2un-1∂y2-2∂4u0∂x2∂y2]dτ,un(x,y)=0alongtheboundaries,n≥2.
It is worth noting that if the zeroth component u0(x,y) is defined, then the remaining components n≥1 can be completely determined such that each term is determined by using the previous terms, and the series solutions are thus entirely determined. Finally, the solution u(x,y) is approximated for n=4:
(34)u(x,y)=sin(3πy)[x108π3-x372π+πx5160-3π3x72240+3π5x917920=sin(3πy555)-27π7x111971200+81π9x13102502400-81π11x157175168000],(35)u(x,y)=sin(3πy)324π4[3πx-(3πx)33!+(3πx)55!-(3πx)77!+(3πx)99!=sin(3πy555)-(3πx)1111!+(3πx)1313!-(3πx)1515!].
Therefore for any N≥4 we have the following:
(36)uN(x,y)=sin(3πy)324π4∑n=0N(3πx)2n+1(2n+1)!.
Thus
(37)limN→∞uN(x,y)=sin(3πy)sin(3πx)324π4.
The exact solution of (31) is given by
(38)sin(3πy)sin(3πx)324π4=u(x,y).
Figures 3 and 4 are the graphical representation of the previous solution. We have plotted the solution for (31) in Figure 3 and showed absolute value of the solution in Figure 4.Figure 3
Analytical solution.Figure 4
Absolute value of the solution.Theorem 2.
Letm be a nonzero natural number and let (x,y)∈[0,1]×[0,1]; then two dimensional biharmonic equation of form
(39)∂4u∂x4+2∂4u∂x2∂y2+∂4u∂y4=sin(mπx)sin(mπy)
with u(x,y)=0 along the boundaries has an exact solution as follows
(40)u(x,y)=sin(mπx)sin(mπy)4m4π4.Proof.
Use the step of the homotopy decomposition method.Problem 5.
We consider the 3D biharmonic equation:(41)∂4u∂x4+2∂4u∂x2∂y2+2∂4u∂x2∂z2+2∂4u∂z2∂y2+∂4u∂y4+∂4u∂z4=sin(πx)sin(πy)sin(πz),0≤x,y,z≤1,u=0,ux,x=uy,y=uz,z=0.
In the view of the homotopy decomposition method, the following integral equations are obtained:
(42)p0:u0(x,y)=sin(3πy)108π3x-sin(3πy)12πx3!3,p1:u1(x,y)=∫0x(x-τ)[∂4u0∂z2∂y2sin(πτ)sin(πy)sin(πz)p1:(x-τ5555555)-2∂4u0∂x2∂y2-2∂4u0∂x2∂z2p1:(x-τ5555555)-2∂4u0∂z2∂y2-∂4u0∂y4-∂4u0∂z4]dτ,pn:un(x,y)=∫0x(x-τ)[-2∂4un-1∂x2∂y2-2∂4un-1∂x2∂z2p1:(x-τ5555555)-2∂4un-1∂z2∂y2-∂4un-1∂y4-∂4un-1∂z4]dτ,un(x,y)=0alongtheboundaries,n≥2.
Solving the previous integral equations, the series solutions for the first N terms are given as
(43)uN(x,y,z)=sin(zπ)sin(πy)9π4∑n=0N(πx)2n+1(2n+1)!.
Therefore taking the limit at N tending to infinity we obtained
(44)u(x,y,z)=limN→∞uN(x,y,z)=sin(xπ)sin(zπ)sin(πy)9π4.


## 4. Conclusion

In this paper the recent homotopy decomposition [18–21] is used to solve the 2D and 3D Poisson equations and biharmonic equations. The method is chosen because it does not require the linearization or assumptions of weak nonlinearity, the solutions are generated in the form of general solution, and it is more realistic compared to the method of simplifying the physical problems. The method does not require any corrected function any Lagrange multiplier and it avoids repeated terms in the series solutions compared to the existing decomposition method including the variational iteration method and the Adomian decomposition method. The approximated solutions obtained converge to the exact solution as N tends to infinity. The numerical values are presented in Table 1 shows that the method is very efficient and accurate.



---

*Source: 380484-2013-09-24.xml*