Homotopy Analysis Method Solution to Time-Fractional Diffusion with a Moving Boundary

Ogugua N. Onyejekwe

Department of Mathematics, Indian River State College, Fort Pierce, FL, USA.

DOI: 10.4236/jamp.2025.136116                

Abstract

It is difficult to obtain exact or analytical solutions to most moving boundary problems. In this paper, we employ the use of Homotopy Analysis Method (HAM) to solve a time-fractional diffusion equation with a moving boundary. HAM is a semi-analytic technique used to solve ordinary, partial, algebraic, delay and fractional differential equations. This method uses the concept of homotopy from topology to generate a convergent series solution for nonlinear systems. The homotopy Maclaurin series is utilized to deal with nonlinearities in the system.

Keywords

Diffusion Equation, Fractional Derivatives, Homotopy Analysis Method, Moving Boundary

1. Introduction

The moving boundary problem is unique because it is a special nonlinear problem whose analytic solution is difficult to obtain. A few numerical and semi analytic methods such as the Homotopy Perturbation Method (HPM), Adomian Decomposition Method (ADM) and Variational Iteration Method (VIM) [1]-[5], have been used to solve some moving boundary problems. In this paper we employ the use of Homotopy Analysis Method (HAM), HAM was first developed by Dr. Shijun Liao in 1992 [6] [7]. He introduced a convergent control parameter $h$ which guarantees convergence for both linear and nonlinear equations. The advantage of using HAM is that we have the freedom to choose the initial guess $u_0\left(x,\tau \right)$ and the value of the convergence control parameter $h$ .

This paper is organized as follows; in Section 2 a short introduction to the mathematical model of drug release and the properties of the Caputo Fractional derivative operator is given. The governing equations to the problem are presented in Section 3. In Section 4, we introduce the HAM solution to the given mathematical model. The comparisons between the approximate and exact solutions are shown in Section 5. Finally, in Section 6, we provide a summary and conclusion to the proposed problem.

2. The Mathematical Model and its Parameters

We consider the model by Liu and Xu [4]. The concentration profile of the drug at time $t$ is shown in Figure 1, where $R$ is the scale of the polymer matrix. Each matrix consists of two regions; the surface zone, $0<\xi <p\left(t\right)$ , in which all solute is dissolved, and the core, $p\left(t\right)<\xi <R$ , where the undissolved solute is contained. The two zones are separated by the diffusion front $x=p\left(t\right)$ , which moves inward as time progresses. $T_0$ and $T_s$ are the initial concentration and the solubility of the drug in the tissue field respectively. In this paper, we will only consider the early stages of loss before the diffusion front moves to $R$ and assume $T_0>T_s$ . A condition of perfect sink is assumed.

Nomenclature

The parameters used in the governing equation are defined as follows:

$p\left(t\right)$ —diffusion front

$T_0$ —initial concentration of drug distributed in matrix

$T_s$ —solubility of drug in the matrix

$T\left(\xi ,t\right)$ —concentration of the drug in the matrix

$\varrho$ —diffusivity of drug in the matrix (assumed to be constant)

${}_0{}^{c}D{}_t^{\alpha }$ —Caputo derivative

$R$ —scale of the polymer matrix

Figure 1. Profile of concentration of drug in the matrix at time t.

3. The Governing Equation

The governing equation and the posed conditions are as follows:

${}_0{}^{c}D{}_t^{\alpha }T\left(\xi ,t\right)=\varrho T_{\xi \xi }\left(\xi ,T\right)$ (1)

with the initial condition

$T\left(\xi ,0\right)=0$ (2)

and the following boundary conditions

$T\left(p\left(1\right),1\right)=k_2,T\left(p\left(t\right),t\right)=T_s,t>0$ (3)

$p\left(1\right)=k_1$ (4)

The values for $p\left(1\right)$ depends on the values of $\alpha$ and $\varrho$ and can be found in the paper written by X.Li and M.Yu [4].

$\left(T_0-T_s\right){}_0{}^{c}D{}_t^{\alpha }p\left(t\right)=\varrho T_{\xi }\left(\xi ,t\right),\xi =p\left(t\right),t>0$ (5)

${}_0{}^{c}D{}_t^{\alpha }$ is defined as the Caputo derivative

${}_0{}^{c}D{}_t^{\alpha }f\left(t\right)={\displaystyle {\int }_0^t\frac{{\left(t-\tau \right)}^{n-\alpha -1}}{\text{Γ}\left(n-\alpha \right)}{f}^{\left(n\right)}\left(\tau \right)\text{d}\tau },\alpha >0$ (6)

for $n-1<\alpha <n,n\in N$ and $\text{Γ}(.)$ represents the Gamma function.

4. The HAM Solution to the Equations

The reduced dimensionless variables are defined as

$x=\frac{\xi }{R},\tau ={\left(\frac{\varrho }{R^2}\right)}^{\frac{1}{\alpha }}t,u=\frac{T}{T_s},P\left(\tau \right)=\frac{p\left(t\right)}{R}$ (7)

The governing Equation (1) subjected to conditions (2)-(5) can be reduced to the dimensionless forms

${}_0{}^{c}D{}_t^{\alpha }u\left(x,\tau \right)=\frac{{\partial }^{2}u\left(x,\tau \right)}{\partial x^2},0<x<P\left(\tau \right)$ (8)

$u\left(0,\tau \right)=0$ (9)

$u\left(P\left(\tau \right),\tau \right)=1,\tau >0$ (10)

$\frac{\partial u\left(x,\tau \right)}{\partial x}=\eta {}_0{}^{c}D{}_t^{\alpha }P\left(\tau \right),x=P\left(\tau \right),\tau >0$ (11)

where $\eta =\frac{T_0-T_s}{T_s}$ .

To solve Equation (11) by homotopy analysis method (HAM), the initial guess for $u\left(x,\tau \right)$ is chosen as

$u_0\left(x,\tau \right)={\left(a_0\right)}^{-1}x{\tau }^{\gamma }$ (12)

where $a_0={\left(\frac{\Gamma \left(1-\frac{\alpha }{2}\right)}{\eta \Gamma \left(1+\frac{\alpha }{2}\right)}\right)}^{\frac{1}{2}},\gamma =-\frac{\alpha }{2}$ .

The auxiliary linear operator is

$ℒ\left[\varphi \left(x,\tau ;q\right)\right]={\varphi }_{xx}\left(x,\tau \right)$ (13)

with the property

$ℒ\left[c\right]=0$ (14)

where $c$ is the integral constant, $\varphi \left(x,\tau ;q\right)$ is an unknown function.

The nonlinear operator is given as

$\text{Ν}\left[\varphi \left(x,\tau ;q\right)\right]={\varphi }_{xx}\left(x,\tau ;q\right)-\frac{{\partial }^{\alpha }\varphi \left(x,\tau ;q\right)}{\partial {\tau }^{\alpha }}$ (15)

By means of HAM, defined by S. Liao [8] [9], we construct a zeroth-order deformation

$\left(1-q\right)ℒ\left[\varphi \left(x,\tau ;q\right)-u_0\left(x,\tau \right)\right]=qh\text{Ν}\left[\varphi \left(x,\tau ;q\right)\right]$ (16)

where $q\in \left[0,1\right]$ is the embedding parameter, $h\ne 0$ , is the convergence-control parameter and $u_0\left(x,\tau ;q\right)$ is the initial/best guess $u_0\left(x,\tau \right)$ .

Expanding $\varphi \left(x,\tau ;q\right)$ in Taylor series with respect to $q$ we obtain

$\varphi \left(x,\tau ;q\right)=u_0\left(x,\tau \right)+{\displaystyle \sum }_{m=1}^{\infty }{u}_m\left(x,\tau \right)q^m$ (17)

Clearly, we see that when $q=0$ and $q=1$ respectively, Equation (17) becomes

$\varphi \left(x,\tau ;0\right)=u_0\left(x,\tau \right),\varphi \left(x,\tau ;1\right)=u\left(x,\tau \right)$ (18)

If the auxiliary linear operator $ℒ$ , the initial guess $u_0\left(x,\tau \right)$ and the convergence-control parameter $h$ are properly chosen so that the series described in equation (17) converges at $q=1$ , then $u\left(x,\tau \right)$ will be one of the solutions of the problem we have considered.

Using the following property

$u_m\left(x,\tau \right)={\frac{1}{m!}\frac{{\partial }^m}{\partial x^m}\varphi \left(x,\tau ;q\right)|}_{q=0}$ (18)

On Equation (16), we obtain an mth-order deformation equation

$ℒ\left[u_m\left(x,\tau \right)-{\chi }_m{u}_{m-1}\left(x,\tau \right)\right]=h{V}_m\left[u_{m-1}\left(x,\tau \right)\right]$ (19)

where

${\chi }_m=\{\begin{array}{l}0,m\le 1\\ 1,m>1\end{array}$

and

$V_m\left[u_{m-1}\left(x,\tau \right)\right]=\frac{{\partial }^2{u}_{m-1}\left(x,\tau \right)}{\partial x^2}-\frac{{\partial }^{\alpha }{u}_{m-1}\left(x,\tau \right)}{\partial {\tau }^{\alpha }}$ (20)

We have

$u_m\left(x,\tau \right)={\chi }_m{u}_{m-1}\left(x,\tau \right)+h{ℒ}^{-1}{V}_m\left[u_{m-1}\left(x,\tau \right)\right]+c_1$ (21)

and the integration constant $c_1$ is determined by the boundary condition given in Equation (10).

Looking at Equation (21), the values for $u_m\left(x,\tau \right)$ for $m=1,2,3,\cdots$ can be obtained and the series solution obtained.

The approximate analytic solution is gained by truncating the following series

$u\left(x,\tau \right)={\displaystyle \sum }_{i=0}^m{u}_i\left(x,\tau \right)$ (22)

Equation (22) contains the convergence-control parameter $h$ , which determines the convergence region and convergence rate of the homotopy-series solution.

In this paper, the convergence-control parameter $h$ is obtained by setting $u{\left(P\left(1\right),1\right)}_{exact}=u{\left(P\left(1\right),1\right)}_{HAM}$ unlike other papers [8] [9] where $h$ is determined by the minimum of residual square of the original governing equation. The approximation for the diffusion front $u{\left(P\left(\tau \right),\tau \right)}_{HAM}=1$ .

5. Comparison and Analysis of Results

The exact solutions for $u\left(x,\tau \right)$ and $P\left(\tau \right)$ [2] [3] respectively are given as follows

$u\left(x,\tau \right)=H{\displaystyle \sum }_{n=0}^{\infty }\frac{{\left(\frac{x}{{\tau }^{\alpha /2}}\right)}^{2n+1}}{\left(2n+1\right)!\text{Γ}\left(1-\frac{2n+1}{2}\alpha \right)}$ (23)

where the values for $H$ are calculated from the following equation

$H={\left(\frac{\text{Γ}\left(1-\frac{\alpha }{2}\right)}{\eta \text{Γ}\left(1+\frac{\alpha }{2}\right)}\right)}^{-\frac{1}{2}}\left(\text{Γ}\left(1-\frac{\alpha }{2}\right)\right)$ (24)

$P\left(\tau \right)=\overline{p}{\tau }^{\alpha /2}$ (25)

In which $\overline{p}$ and $H$ satisfy the following

$H{\displaystyle \sum }_{n=0}^{\infty }\frac{{\left(\overline{p}\right)}^{2n+1}}{\left(2n+1\right)!\text{Γ}\left(1-\frac{2n+1}{2}\alpha \right)}=1$ (26)

and

$H{\displaystyle \sum }_{n=0}^{\infty }\frac{{\left(\overline{p}\right)}^{2n}}{\left(2n\right)!\Gamma \left(1-\frac{2n+1}{2}\alpha \right)}=\overline{p}\eta \left(\frac{\text{Γ}\left(1+\frac{\alpha }{2}\right)}{\text{Γ}\left(1-\frac{\alpha }{2}\right)}\right)$ (27)

To evaluate the accuracy of HAM, the Error Analysis is shown for both $u\left(x,\tau \right)$ and $P\left(\tau \right)$ respectively. It can be observed that $u\left(x,\tau \right)$ and $P\left(\tau \right)$ are most accurate when $\alpha =0.75$ . The values of both Table 1, and Table 2 are not affected by the values of $x$ and $\tau$ . When looking at Table 2, we see that for any value of $\alpha$ the Relative Error (RE) decreases as $\eta$ increases.

Table 1. Error analysis for $u\left(x,\tau \right)$ for $x\ge 0,\tau \ge 0$ and $\eta >0$ where $h=5\times {10}^{-5}$ .

$\alpha$	Relative Error (RE)
0.5	0.0514
0.75	0.0295
1	0.0839

Table 2. Error analysis for $P\left(\tau \right)$ for $\tau \ge 0$ where $h=5\times {10}^{-5}$ .

$\eta$	$\alpha =0.5$ Relative Error (RE)	$\alpha =0.75$ Relative Error (RE)	$\alpha =1.0$ Relative Error (RE)
0.5	0.1000	0.1080	0.2490
1	0.0640	0.0505	0.1403
3	0.0243	0.0156	0.0519
5	0.0150	0.0089	0.0319
7	0.0110	0.0064	0.0230
9	0.0085	0.0049	0.0179

5. Conclusion

We have shown that when the initial value for $u_0\left(x,\tau \right)$ is well chosen, Homotopy Analysis Method (HAM) can be used to accurately predict drug distribution $u\left(x,\tau \right)$ for different values of $\alpha$ and $\eta$ . Computational efficiency and a strong rate of convergence can be observed.

Conflicts of Interest

The author declares no conflicts of interest regarding the publication of this paper.

References

[1]	Zheng, M., Liu, F., Liu, Q., Burrage, K. and Simpson, M.J. (2017) Numerical Solution of the Time Fractional Reaction-Diffusion Equation with a Moving Boundary. Journal of Computational Physics, 338, 493-510.[CrossRef]
[2]	Das, S. and Rajeev (2010) Solution of Fractional Diffusion Equation with a Moving Boundary Condition by Variational Iteration Method and Adomian Decomposition Method. Zeitschrift für Naturforschung, 65, 793-799.
[3]	Gao, X.L., Jiang, X.Y. and Chen, S.Z. (2015) The Numerical Method for the Moving Boundary Problem with Space-Fractional Derivative in Drug Release Devises. Applied Mathematical Modelling, 39, 2385-2391.
[4]	Li, X.C., Xu, M.Y. and Jiang, X.Y. (2009) Homotopy Perturbation Method of Time-Fractional Diffusion Equation with a Moving Boundary Condition. Applied Mathematics and Computation, 208, 434-439.
[5]	Li, X., Wang, S. and Zhao, M. (2013) Two Methods to Solve a Fractional Single Phase Moving Boundary Problem. Central European Journal of Physics, 11, 1387-1391.
[6]	Liao, S. (2012) Homotopy Analysis Method in Nonlinear Equations. Springer.
[7]	Liao, S. (2014) Advances in the Homotopy Analysis Method. World Scientific.
[8]	Kumar, A. and Rajeev (2020) A Moving Boundary Problem with Space-Fractional Diffusion Logistic Population Model and Density-Dependent Dispersal Rate. Applied Mathematical Modelling, 88, 951-965.[CrossRef]
[9]	Chmielowska, A. and Słota, D. (2022) Fractional Stefan Problem Solving by the Alternating Phase Truncation Method. Symmetry, 14, Article No. 2287.[CrossRef]