Journal of Mathematics

Volume 2018, Article ID 8924547, 9 pages

https://doi.org/10.1155/2018/8924547

## Application of Nonlinear Time-Fractional Partial Differential Equations to Image Processing via Hybrid Laplace Transform Method

^{1}School of Computer Science and Applied Mathematics, University of the Witwatersrand, Johannesburg, Private Bag 3, Wits 2050, South Africa^{2}DST-NRF Centre of Excellence in Mathematical and Statistical Sciences (CoE-MaSS), South Africa^{3}School of Mathematics and Statistics, UNSW Australia, Sydney, NSW 2052, Australia

Correspondence should be addressed to B. A. Jacobs; az.ca.stiw@sbocaj.noryb

Received 6 April 2018; Accepted 30 August 2018; Published 17 September 2018

Academic Editor: Morteza Khodabin

Copyright © 2018 B. A. Jacobs and C. Harley. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### Abstract

This work considers a hybrid solution method for the time-fractional diffusion model with a cubic nonlinear source term in one and two dimensions. Both Dirichlet and Neumann boundary conditions are considered for each dimensional case. The hybrid method involves a Laplace transformation in the temporal domain which is numerically inverted, and Chebyshev collocation is employed in the spatial domain due to its increased accuracy over a standard finite-difference discretization. Due to the fractional-order derivative we are only able to compare the accuracy of this method with Mathematica’s NDSolve in the case of integer derivatives; however, a detailed discussion of the merits and shortcomings of the proposed hybridization is presented. An application to image processing via a finite-difference discretization is included in order to substantiate the application of this method.

#### 1. Introduction

This work examines the performance of a hybrid Laplace transform–Chebyshev collocation technique applied to the time-fractional diffusion equation in two dimensions with a nonlinear source term:where This model is explored subject to both Dirichlet and Neumann boundary conditions on the bounded domain to satisfy the domain required by the Chebyshev polynomials with initial condition . This method benefits from the analyticity of the Laplace transform and efficient numerical inversion of this transform, an accurate discretization approach through Chebyshev collocation, and a convergent linearization technique, which results in a robust method for solving nonlinear time-fractional partial differential equations on a bounded domain. Following the method detailed by Jacobs and Harley in [1, 2] we also employ a finite-difference discretization in applying this method to images. This is elaborated on further in Section 5.

Recently fractional derivatives and fractional partial differential equations (FPDEs) have received great attention both in analysis and application (see [3–5] and references therein). Despite this large effort, very little attention has been paid to solving FPDEs on a bounded domain through transformation techniques. Agrawal [6] makes use of the Laplace transform and the finite sine transform to obtain an analytic solution to the fractional diffusion-wave equation on a bounded domain. Other techniques such as the variational iteration method, Adomian decomposition method, and differential transform method have all been applied to fractional partial differential equations but with a focus on unbounded domains. More recently Hashemizadeh and Ebrahimzadeh [7] and Zaky et al. [8] present solutions to the linear and nonlinear time-fractional Klein-Gordon equations, respectively. Their solutions were obtained by using an operational matrix approach and Legendre polynomial interpolation.

The linear diffusion model has been applied in many different fields including image processing [9–14]. However, to the best of the authors’ knowledge, the application of a time-fractional partial differential equation has not yet been thoroughly examined within such a context. It is from the standpoint of applying a fractional partial differential equation to an image that we examine the time-fractional diffusion equation in two dimensions. In a recent paper Jacobs and Momoniat [15, 16] show that the diffusion equation with this nonlinear source term is able to binarize a document image with great success. In extending this concept to fractional-order derivatives, we seek to preserve the symmetry of the diffusion equation. Values of resolve the diffusion equation to subdiffusion, which preserves symmetry and exhibits some interesting dynamics which are discussed later on in this work. Values of introduce a transport effect which then breaks symmetry. It is for this reason that we impose the restriction, .

In this work we make use of the Laplace transform which allows us to handle the fractional-order derivative in an algebraic way and incur no error in doing so. Inversion of the Laplace transform is however difficult to obtain analytically. We therefore make extensive use of the numerical inversion procedure described by Weideman and Trefethen in [17] which defines a contour of integration that maps the domain of the Bromwich integral from the entire complex space to the real space, from which we can approximate this integral with a trapezoidal rule.

With a robust method for inverting the Laplace transform we may then hybridize the transform with a discretization technique. The Laplace transform of the temporal variable avoids the need for a time-marching scheme as well as reducing the fractional-order derivative to an algebraic expression. The transformed model is then discretized by the use of Chebyshev collocation due to its superior performance to finite differences as illustrated in [1, 2]. The resulting system is solved and the transform inverted to obtain a semianalytic solution, continuous in time and discrete in space.

In the following section we present some preliminary results which are put to use throughout the paper. In Section 3 the implemented methods are described, including the different cases for boundary conditions. Section 4 presents the solutions obtained by the proposed method as well as an error comparison with Mathematica’s NDSolve for the one- and two-dimensional cases of our model as well as both Dirichlet and Neumann boundary conditions. In Section 5 we provide a real-world application of this method and model in the form of document image binarization, illustrating the ability to obtain a useful result with the present method when coupled with the finite-difference discretization. A discussion of the results and their relationship to work beyond this research is presented in Section 6 and some concluding remarks are made in Section 7.

#### 2. Preliminaries

In this work we employ Caputo’s definition of a fractional derivative over the Riemann-Liouville derivative due to the fact that the Caputo derivative makes use of the physical boundary conditions, whereas the Riemann-Liouville derivative requires fractional-order boundary conditions.

*Definition 1. *The Riemann-Liouville integral of order of a function is

*Definition 2. *The fractional derivative of according to the Caputo definition with , is

If , the Caputo fractional derivative reduces to the ordinary derivative or integral. Podlubny [4] illustrates the pleasing property of the Laplace transform of a Caputo derivative, as can be seen in (5). In our case where we haveThis property allows one to treat fractional-order derivatives algebraically.

*Definition 3. *The Generalized Mittag-Leffler function of the argument is

The use of Laplace transform allows one to circumvent the problems that arise in the time-domain discretization. However, using the Laplace transform for the fractional-order derivative presents the problem of inverting the transform to find a solution. Analytic inversion of the transform is infeasible and hence the numerical scheme for evaluating the Bromwich integral presented by Weideman and Trefethen in [17] is put to extensive use. The Bromwich integral is of the formwhere is the convergence abscissa. Using the parabolic conformal map(7) becomeswhich can then be approximated simply by the trapezoidal rule, or any other quadrature technique, asOn the parabolic contour the exponential factor in (7) forces rapid decay in the integrand making it amenable to quadrature. All that is left is to choose the parameters and optimally which is done by asymptotically balancing the truncation error and discretization errors in each of the half-planes. The methodology described by Weideman and Trefethen achieves near optimal results (see [17] and figures therein) and the interested reader is directed there for a thorough description of the implementation of the method. In this work we make use of the parabolic contour due to the ease of use and the hyperbolic contour only exhibits a slight improvement in performance over the parabolic contour.

#### 3. Methods

This section introduces the methodologies used for the two-dimensional model previously presented. We may write our model asThe quasi-linearization technique can be viewed as a generalized Newton-Raphson method in functional space. An iterative scheme is constructed creating a sequence of linear equations that approximate the nonlinear equation (11) and boundary conditions. Furthermore, this sequence of solutions converges quadratically and monotonically [18–20]. The quasi-linear form is given bywhere indicates the index of successive approximation. Moreover, is known entirely at the explicit index . The coefficients are given byandIf the elements in are indexed by , then for the th iteration in the quasi-linearization we haveSince the first term above satisfies (11), it is replaced with 0 givingEquation (12) can now be transformed by the Laplace transform, a linear operator, to obtainwhereEquation (17) may be discretized by Chebyshev collocation, described in the following section.

##### 3.1. Chebyshev Collocation

Chebyshev polynomials form a basis on and hence we dictate the domain of our PDE to be . We note here, however, that any domain in can be trivially deformed to match . We discretize our spatial domain using Chebyshev-Gauss-Lobatto points:Given this choice of spatial discretization, we have , , , and , indicating that the domain is in essence reversed and one must exercise caution when imposing the boundary conditions.

In mapping our domain to we may assume that ; i.e., we have equal number of collocation points in each spatial direction. We now define a differentiation matrix :Bayliss et al. [21] describe a method for minimizing the round-off errors incurred in the calculations of higher order differentiation matrices. Since we write , we implement the method, described in [21], in order to minimize propagation of round-off errors for the second derivative in space.

The derivative matrices in the direction arewhere is the Chebyshev differentiation matrix of size .

Because we have assumed , we derive the pleasing property that and .

Writing the discretization of (17) in matrix form yieldswhereBy expanding (22) in summation notation, we havefor . By extracting the first and last terms in sums, we obtainfor . We use the form of (25) to impose the boundary conditions.

The solution , , which is the matrix of unknown interior points of , can be found by solving the systemwhere is the matrix of interior points of and is the matrix of interior points of , so that and match the dimensions of . We also use to denote the Hadamard product between two matrices. Alsofor .

###### 3.1.1. Dirichlet Boundary Conditions

Boundary conditions may be in the form of Dirichlet conditions,and hence,The parameters , , , and are potentially functions of the temporal variable and one of the spatial variables; i.e., . We assume that , , , and are constant.

Dirichlet boundary conditions can be imposed directly by substituting (28) and (29) into (25) and collecting all the known terms in .

###### 3.1.2. Neumann Boundary Conditions

Alternatively Neumann boundary conditions givewith,

Neumann boundary conditions given by (30) are discretized as Similarly for equations (31)By extracting the first and last terms in the sum, the discretizations can be written asandand the solutions to these linear systems are then substituted into equation (25).

#### 4. Results

In this section we consider only the results obtained by Chebyshev collocation due to the enormous increase in accuracy obtained over the finite-difference method that was presented by the authors in [1, 2].

##### 4.1. Example 1

The first example we consider iswith Dirichlet boundary conditions consistent with this initial condition,Table 1illustrates the maximum absolute error between the solution obtained by the present method and the solution obtained by NDSolve in Mathematica 9 for . To the best of the authors’ knowledge, no exact solution exists for the fractional case and, hence, no comparison can be made.