# The New Approximate Analytic Solution for Oxygen Diffusion Problem with Time-Fractional Derivative

Authors: Vildan Gülkaç
Journal: Mathematical Problems in Engineering
 (2016)

---

## Abstract

Oxygen diffusion into the cells with simultaneous absorption is an important problem and it is of great importance in medical applications. The problem is mathematically formulated in two different stages. At the first stage, the stable case having no oxygen transition in the isolated cell is investigated, whereas at the second stage the moving boundary problem of oxygen absorbed by the tissues in the cell is investigated. In oxygen diffusion problem, a moving boundary is essential feature of the problem. This paper extends a homotopy perturbation method with time-fractional derivatives to obtain solution for oxygen diffusion problem. The method used in dealing with the solution is considered as a power series expansion that rapidly converges to the nonlinear problem. The new approximate analytical process is based on two-iterative levels. The modified method allows approximate solutions in the form of convergent series with simply computable components.


---

## Body

## 1. Introduction

The diffusion of oxygen into absorbing tissue was first studied in [1]. First the oxygen is allowed to diffuse into a medium, some of the oxygen can be absorbed by the medium, and concentration of oxygen at the surface of the medium is maintained constant. This phase of the problem continues until a steady state is reached in which the oxygen does not penetrate any further and is sealed so that no oxygen passes in or out, the medium continues to absorb the available oxygen already in it, and, as a consequence, the boundary in the steady state starts to recede towards the sealed surface.Crank and Gupta [2] also employed uniform space grid moving with the boundary and necessary interpolations are performed with either cube splines or polynomials. Noble [3] suggested repeated spatial subdivision, Reynolds and Dolton [4] also developed the heat balance integral method, and Liapis et al. [5] proposed an orthogonal collocation for solving the partial differential equation of the diffusion of oxygen in absorbing tissue. Gülkaç proposed two numerical methods for solving the oxygen diffusion problem [6]. Mitchell studied the accurate application of the integral method [7]. More references to this problem may be found in [8–17].In recent years, fractional differential equations have drawn much attention. Many important phenomena in physics, engineering, mathematics, finance, transport dynamics, and hydrology are well characterized by differential equations of fractional order. Fractional differential equations play an important role in modelling the so-called anomalous transport phenomena and in the theory of complex systems. These fractional derivatives work more appropriately compared with the standard integer-order models. So, the fractional derivatives are regarded as very dominating and useful tool. For mathematical properties of fractional derivatives and integrals one can consult [18–23].In the present work, we extend a homotopy perturbation method with time-fractional derivatives to obtain solution for oxygen diffusion problem.We give some basic definitions of fractional derivatives as follows.Definition 1.
The Riemann-Liouville fractional integral off
∈
C
α of the order α
≥
0 is defined as (1)
J
t
α
f
t
=
f
t
,
i
f
α
=
0
1
Γ
α
∫
0
t
t
-
τ
α
-
1
f
τ
d
τ
,
i
f
α
>
0
,where Γ denotes gamma function:(2)
Γ
z
=
∫
0
∞
e
-
t
t
z
-
1
d
t
,
z
∈
C
.Definition 2.
The fractional derivatives off
∈
C
α of the order α
≥
0, in Caputo sense, are defined as(3)
D
t
α
f
t
=
J
t
n
-
α
D
t
α
f
t
=
1
Γ
n
-
α
∫
0
t
t
-
τ
n
-
α
-
1
f
n
τ
d
τfor n
-
1
<
α
≤
n, n
∈
N, t
>
0, f
∈
C
α
n, α
≥
-
1.Definition 3.
The Caputo-time-fractional derivative operator of orderα
>
0 is defined as(4)
D
t
α
u
x
,
t
=
J
t
n
-
α
u
x
,
t
=
1
Γ
n
-
α
∫
0
t
t
-
τ
n
-
α
-
1
∂
n
u
∂
τ
n
d
τ
.Lemma 4.
Letn
-
1
<
α
≤
n, n
∈
N, and f
∈
C
α
n, α
≥
-
1; then(5)
D
t
α
J
t
α
f
t
=
f
t
,
J
t
α
D
t
α
=
f
t
-
∑
k
=
0
n
f
k
0
+
t
k
k
!
,
f
o
r
t
>
0
.Lemma 5.
Ifn
-
1
<
α
≤
n, n
∈
N, and k
≥
0 then one has (6)
J
t
α
t
k
α
Γ
k
α
+
1
=
t
k
α
+
α
Γ
k
α
+
α
+
1
.The Mittag-Leffler function plays a very important role in the fractional differential equations, in fact introduced by Mittag-Leffler in 1903 [24]. Mittag-Leffler function E
α
z
=
∑
n
=
0
∞
z
n
/
Γ
α
n
+
1.


## 2. Analysis of Homotopy Perturbation Method with Time-Fractional Derivatives

Let us assume that nonlinear fractional differential equation is as follows:(7)
D
t
α
u
x
,
t
=
A
u
+
f
x
,
t
,
x
,
t
∈
Ω
,with the initial condition u
(
x
,
0
)
=
φ, where A is the operator, f is known functions, and u
x
,
t is sough functions. Assume that operator A can be written as A
u
=
L
(
u
)
+
N
(
u
), where L is the linear operator and N is the nonlinear operator. Hence (7) can be written as follows:(8)
D
t
α
u
x
,
t
=
L
u
+
N
u
+
f
x
,
t
.For solving (7) by homotopy perturbation method, we construct the homotopy (9)
H
V
,
p
=
1
-
p
D
t
α
V
-
u
0
+
p
D
α
V
-
L
V
-
N
V
-
f
x
,
t
=
0or the equivalent one(10)
H
V
,
p
=
D
t
α
V
-
u
0
+
p
u
0
-
L
V
-
N
V
-
f
x
,
t
=
0
,where p
∈
[
0,1
] is an embedding or homotopy parameter, H
x
,
t
;
p
:
Ω
x
[
0,1
]
→
R, and u
0 is the initial approximation for solution equation (8).Clearly, the homotopy equationsH
V
,
0
=
0 and H
V
,
1
=
1 are equivalent to the equations D
t
α
V
-
u
0
=
0 and D
α
V
-
L
V
-
N
V
-
f
x
,
t
=
0, respectively. Thus, a monotonous change of parameter p from 0 to 1 corresponds to a continuous change of the trivial problem D
t
α
V
-
u
0
=
0 to the original problem. Now, we assume that the solution of (8) can be written as a power series in embedding parameter p, as follows:(11)
V
=
V
0
+
p
V
1
,where V
0 and V
1 are functions which should be determined. Now, we can write (11) in the following form:(12)
D
t
α
V
x
,
t
=
u
0
+
p
-
u
0
+
L
V
+
N
V
+
f
x
,
t
.Apply the inverse operator, J
t
α, which is the Riemann-Liouville fractional integral of order α
>
0.On both sides of (12), we have(13)
V
x
,
t
=
V
x
,
0
+
J
t
α
u
0
+
p
J
t
α
-
u
0
+
L
V
+
N
V
+
f
x
,
t
.Suppose that the initial approximation of solutions equation (8) is in the following form:(14)
u
0
=
∑
k
=
0
∞
a
k
x
t
k
α
Γ
k
α
+
1
,where a
k
(
x
) for k
=
1,2
,
… are functions which must be computed. Substituting (11) and (14) into (13) we get (15)
V
0
+
p
V
1
=
V
x
,
0
+
J
t
α
∑
k
=
0
∞
a
k
x
t
k
α
Γ
k
α
+
1
+
p
J
t
α
-
∑
k
=
0
∞
a
k
x
t
k
α
Γ
k
α
+
1
+
L
V
0
+
p
V
1
+
N
V
0
+
p
V
1
+
f
x
,
t
.Synchronizing the coefficients of the same powers leads to(16)
p
0
:
V
0
=
V
x
,
0
+
J
t
α
∑
k
=
0
∞
a
k
x
t
k
α
Γ
k
α
+
1
,
p
1
:
V
1
=
J
t
α
-
∑
k
=
0
∞
a
k
x
t
k
α
Γ
k
α
+
1
+
L
V
0
+
p
V
1
+
N
V
0
+
p
V
1
+
f
x
,
t
.Now, we obtain the coefficients a
k
x, k
=
1,2
,
…. Therefore the exact solution can be obtained as follows:(17)
u
x
,
t
=
V
x
,
t
=
V
x
,
0
+
J
t
α
∑
k
=
0
∞
a
k
x
t
k
α
Γ
k
α
+
1
.Efficiency and reliability of the method are shown.


## 3. Problem Description and Formulation

Crank and Gupta [1] were the first researchers to model oxygen diffusion problem mathematically.The process includes two mathematical levels. At the first level, the stable condition occurs once the oxygen is injected into either the inside or outside of the cell; then the cell surface is isolated.At the second level, tissues start to absorb the injected oxygen. The moving boundary problem is caused by this level. The aim of this process is to find a balance position and to determine the time-dependent moving boundary position. Writing down the time-fractional derivatives of oxygen diffusion problem in [1] is adopted, following [18].


## 4. Solution of Fractional Oxygen Diffusion Problem

We consider the following oxygen diffusion problem:(18)
D
t
α
c
x
,
t
=
c
x
x
-
1
,
x
,
t
∈
Ωwith the following initial and boundary conditions: (19)
c
x
,
0
=
0.5
1
-
x
2
,
0
≤
x
≤
1
,
t
=
0
,
(20)
∂
c
∂
x
=
0
,
x
=
0
,
t
≥
0
,
(21)
c
=
∂
c
∂
x
=
0
,
x
=
s
t
,
t
≥
0
w
i
t
h
s
0
=
1
,where 0
<
α
≤
1. To solve (18)–(21) by present method, we construct the following homotopy:(22)
H
V
,
p
=
1
-
p
D
t
α
V
-
c
0
+
p
D
α
V
-
L
V
-
1
=
0or(23)
H
V
,
p
=
D
t
α
V
-
c
0
+
p
c
0
-
L
V
-
1
=
0
,where p
∈
[
0,1
] and 0
<
α
≤
1. Consider(24)
D
t
α
c
x
,
t
=
c
0
x
,
t
-
p
c
0
x
,
t
-
c
x
x
x
,
t
-
1
.Assume that the initial approximation of solutions equation (18) is in the following form:(25)
c
0
x
,
t
=
∑
n
=
0
∞
a
n
x
t
n
α
Γ
n
α
+
1
,where a
n
(
x
) for n
=
1,2
,
… are functions which must be computed.Applying the inverse operatorJ
t
α of D
t
α on both sides of (24) we obtain(26)
c
x
,
t
=
c
x
,
0
+
J
t
α
c
0
x
,
t
-
p
J
t
α
c
0
x
,
t
-
c
0
x
x
x
,
t
-
1
.Suppose the solution of (26) has the following form: (27)
c
x
,
t
=
c
0
x
,
t
+
p
c
1
x
,
t
,where c
x
,
t is functions which should be determined. Substituting (27) into (26), collecting the same powers of p, and equating each coefficient of p to zero yield(28)
c
0
x
,
t
=
0.5
1
-
x
2
+
J
t
α
c
0
x
,
t
,
(29)
p
0
:
c
0
x
,
t
=
0.5
1
-
x
2
+
∑
n
=
0
∞
a
n
x
t
n
α
+
α
Γ
n
α
+
α
+
1
,
(30)
p
1
:
c
1
x
,
t
=
-
J
t
α
c
0
x
,
t
-
c
0
x
x
x
,
t
-
1
,where(31)
c
0
x
x
x
,
t
=
t
α
Γ
α
+
1
+
t
2
α
Γ
2
α
+
1
+
⋯
+
t
n
α
Γ
n
α
+
1
+
⋯
,
c
0
x
x
x
,
t
=
E
α
t
=
∑
n
=
0
∞
t
n
α
Γ
n
α
+
1
;then (32)
c
1
x
,
t
=
-
∑
n
=
0
∞
a
n
x
t
n
α
+
α
Γ
n
α
+
α
+
1
+
∑
n
=
0
∞
t
n
α
+
α
Γ
n
α
+
α
+
1
,
c
1
x
,
t
=
∑
n
=
0
∞
1
-
a
n
x
t
n
α
+
α
Γ
n
α
+
α
+
1and from (29) and (30) we obtain(33)
a
0
x
=
0.5
1
-
x
2
,
a
0
x
=
a
1
x
=
a
2
x
=
⋯
=
a
n
-
1
x
=
a
n
x
=
⋯
;therefore, we obtain solution of (26) for p
=
1:(34)
c
x
,
t
=
0.5
1
-
x
2
+
∑
n
=
0
∞
a
n
x
t
n
α
+
α
Γ
n
α
+
α
+
1
+
∑
n
=
0
∞
1
-
a
n
x
t
n
α
+
α
Γ
n
α
+
α
+
1or(35)
c
x
,
t
=
0.5
1
-
x
2
+
∑
n
=
0
∞
t
n
+
1
Γ
n
+
1
i
f
n
+
1
=
m
,
n
,
m
∈
N
,
0
<
α
≤
1
,
f
o
r
α
=
1and we obtain the following solution:(36)
c
x
,
t
=
0.5
1
-
x
2
+
e
t
.If α
=
1
/
2 then(37)
c
x
,
t
=
0.5
1
-
x
2
+
e
t
2
erf
⁡
-
t
.If α
=
2 then(38)
c
x
,
t
=
0.5
1
-
x
2
+
s
i
n
t
.We can now obtain an expression for the location of the moving boundary s
(
t
). We can write(39)
d
s
d
t
=
-
∂
c
∂
x
x
=
s
t
a
t
x
=
s
,
s
0
=
1following [18]:(40)
D
t
α
s
=
-
∂
c
∂
x
x
=
s
tand following initial and boundary conditions (41)
s
0
=
1
,
∂
c
∂
x
=
0
,
x
=
0
,
c
=
∂
c
∂
x
=
0
,
x
=
s
t
,
t
≥
0
.We construct the following homotopy for moving boundary as(42)
H
W
,
p
=
1
-
p
D
t
α
W
-
c
0
+
p
D
α
W
-
L
W
=
0or(43)
H
W
,
p
=
D
t
α
W
-
c
0
+
p
c
0
-
L
W
=
0
,where p
∈
[
0,1
] and 0
<
α
≤
1. Consider(44)
D
t
α
s
x
,
t
=
c
0
x
,
t
-
p
c
0
x
,
t
-
c
x
x
,
t
;applying the inverse operator J
t
α of D
t
α on both sides of (44), we obtain(45)
s
x
,
t
=
c
x
,
0
+
J
t
α
c
0
x
,
t
-
p
J
t
α
c
0
x
,
t
-
c
0
x
x
,
t
;suppose the solution of (45) has the following form:(46)
s
x
,
t
=
c
0
x
,
t
+
p
c
1
x
,
t
.Substituting (46) into (45), collecting the same powers of p, and equating each coefficient of p to zero yield(47)
p
0
:
c
0
x
,
t
=
1
+
∑
n
=
0
∞
a
n
x
t
n
α
+
α
Γ
n
α
+
α
+
1
,
p
1
:
c
1
x
,
t
=
-
J
t
α
c
0
x
,
t
-
c
0
x
x
x
,
t
-
1
,where (48)
c
0
x
x
,
t
=
a
0
x
+
∑
n
=
1
∞
a
n
x
t
n
α
+
α
Γ
n
α
+
α
+
1
,
J
t
α
c
0
x
x
,
t
=
J
t
α
∑
n
=
1
∞
t
n
α
Γ
n
α
+
1
=
t
n
α
+
α
Γ
n
α
+
α
+
1
,so(49)
c
1
x
,
t
=
-
a
0
x
t
α
Γ
α
+
1
.Then we have from (46) if we let p
=
1
(50)
c
1
x
,
t
=
a
1
x
t
α
Γ
α
+
1
+
a
2
x
t
2
α
Γ
2
α
+
1
+
⋯
,
a
1
x
=
-
a
0
x
,
a
2
x
=
a
3
x
=
⋯
=
a
n
x
=
⋯
=
0
.Therefore, we obtain the solutions of moving boundary condition as(51)
s
x
,
t
=
∑
n
=
0
1
a
n
x
t
n
α
Γ
n
α
+
1
+
a
1
x
t
α
Γ
α
+
1or (52)
s
x
,
t
=
0.5
1
-
x
2
1
-
t
α
Γ
α
+
1
.


## 5. Numerical Simulations

In this section numerical results for the solution of the oxygen diffusion problem using the constructed homotopy perturbation method with the time-fractional derivative are presented. These proposed homotopy perturbation methods are applied and figures present solutions are presented using different values for the derivative orderα. Figures 1(a) and 1(b) show the surface concentration c
(
x
,
t
) for α
=
1
/
2, Figures 2(a) and 2(b) show the surface concentration c
(
x
,
t
) for α
=
1, and Figures 3(a) and 3(b) show the surface concentration c
(
x
,
t
) for α
=
2. Figures 4(a) and 4(b) show position of moving boundary s
(
x
,
t
) for α
=
1, Figures 5(a) and 5(b) show position of moving boundary s
(
x
,
t
) for α
=
2, and finally Figures 6(a) and 6(b) show position of moving boundary s
(
x
,
t
) for α
=
3.Figure 1
Surface concentrationc
(
x
,
t
) for α
=
1
/
2.
(a)
(b)Figure 2
Surface concentrationc
(
x
,
t
) for α
=
1.
(a)
(b)Figure 3
Surface concentrationc
(
x
,
t
) for α
=
2.
(a)
(b)Figure 4
Position of moving boundarys
(
x
,
t
) for α
=
1.
(a)
(b)Figure 5
Position of moving boundarys
(
x
,
t
) for α
=
2.
(a)
(b)Figure 6
Position of moving boundarys
(
x
,
t
) for α
=
3.
(a)
(b)


## 6. Conclusion

In this study, we extended homotopy perturbation method with time-fractional derivative to find the exact solution of oxygen diffusion problem with moving boundary. It is effortless and also easy to apply and we can say that the present method is an effective method and has appropriate technique to find the exact solution to many complex problems.



---

*Source: 8409839-2016-06-15.xml*