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

Authors: B. A. Jacobs; C. Harley
Journal: Journal of Mathematics
 (2018)

---

## 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 withMathematica’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.


---

## Body

## 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:(1)utα=cduxx+uyy+csu1-uu-awhere (2)utα=∂αu∂tα,uxx=∂2u∂x2,uyy=∂2u∂y2.This model is explored subject to both Dirichlet and Neumann boundary conditions on the bounded domain Ω=[-1,1]×[-1,1] to satisfy the domain required by the Chebyshev polynomials with initial condition u(x,y,0)=f(x,y). 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 α∈(0,1] resolve the diffusion equation to subdiffusion, which preserves symmetry and exhibits some interesting dynamics which are discussed later on in this work. Values of α>1 introduce a transport effect which then breaks symmetry. It is for this reason that we impose the restriction, α∈(0,1].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 Section3 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α>0 of a function u(t) is(3)Jαut=1Γα∫0tt-τα-1uτdτ,x>0.Definition 2.
The fractional derivative ofu(t) according to the Caputo definition with m-1<α≤m,m∈N, is(4)∂αut∂tα=Jm-αDmut=D∗αut≔1Γm-α∫0tt-τm-α-1umτdτ.Ifα∈Z, 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 0<α<1 we have(5)L∂αux,y,t∂tα=sαUx,y-sα-1ux,y,0.This property allows one to treat fractional-order derivatives algebraically.Definition 3.
The Generalized Mittag-Leffler function of the argumentz is(6)Eα,βz=∑k=0∞zkΓkα+β.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 form(7)ut=12πi∫σ-i∞σ+i∞estUsds,σ>σ0,where σ0 is the convergence abscissa. Using the parabolic conformal map(8)s=μiw+12,-∞<w<∞,(7) becomes(9)ut=12πi∫-∞∞eswtUsws′wdwwhich can then be approximated simply by the trapezoidal rule, or any other quadrature technique, as(10)uh;Ntt=h2πi∑k=-NtNteswktUswks′wk,wk=kh.On 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 h 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 as(11)cduxx+cduyy+csa+1u2-csau-csu3-utα=Φuxx,uyy,u,utα=0.The 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 by(12)cduxxn+1+cduyyn+1+X1nun+1+X2nutαn+1=X3n,where n indicates the index of successive approximation. Moreover, u(x,y,t)n is known entirely at the explicit index n. The coefficients are given by(13)X1n=∂Φ∂u=cs2a+1un-csa-3csun2,and(14)X2n=∂Φ∂utα=-1.If the elements in u_=(uxx,uyy,u,utα) are indexed by j, then for the nth iteration in the quasi-linearization we have(15)X3n=∑ju_jn∂Φn∂u_j,=uxx∂Φ∂uxx+uyy∂Φ∂uyy+u∂Φ∂u+utα∂Φ∂utαn,=cduxx+cduyy+2a+1u2-au-3u3-utαn,=cduxx+cduyy+a+1u2-au-u3-utαn+a+1u2-2u3n.Since the first term above satisfies (11), it is replaced with 0 giving(16)X3n=a+1u2-2u3n.Equation (12) can now be transformed by the Laplace transform, a linear operator, to obtain(17)cdUxxn+1+cdUyyn+1+X1nUn+1-sαUn+1=X3ns-sα-1fx,y,where(18)Ux,y,s=Lux,y,t.Equation (17) may be discretized by Chebyshev collocation, described in the following section.


### 3.1. Chebyshev Collocation

Chebyshev polynomials form a basis on[-1,1] and hence we dictate the domain of our PDE to be Ω. We note here, however, that any domain in R2 can be trivially deformed to match Ω. We discretize our spatial domain using Chebyshev-Gauss-Lobatto points:(19)xi=cos⁡iπNx,i=0,1,…,Nx,yj=cos⁡jπNy,j=0,1,…,Ny.Given this choice of spatial discretization, we have x0=1, xNx=-1, y0=1, and yNy=-1, 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 Nx=Ny; i.e., we have equal number of collocation points in each spatial direction. We now define a differentiation matrix D(1)=dkl:(20)dkl=ck-1k+lclxk-xl,k≠l,-xk21-x2k,k=l,162Nx2+1,k=l=0,-162Nx2+1,k=l=Nx,whereck=2,k=0,Nx,1,k=1,…,Nx-1.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 D(2)=D(1).D(1), 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 thex direction are(21)D^x1=D1,D^x2=D2,where D(1) is the Chebyshev differentiation matrix of size (Nx+1)×(Ny+1).Because we have assumedNx=Ny, we derive the pleasing property that D^y(1)=D^x(1)T and D^y(2)=D^x(2)T.Writing the discretization of (17) in matrix form yields(22)D^x2Un+1+Un+1D^y2+X1nUn+1-sαUn+1=X3ns-F,where(23)Fij=sα-1fxi,yj.By expanding (22) in summation notation, we have(24)∑k=0Nxdik2Un+1xk,yj+∑k=0Nydkj2Un+1xi,yk+X1i,jnUn+1xi,yj-sαUn+1xi,yj=X3i,jns-sα-1fxi,yj,for i=0,1,…,Nx,j=0,1,…,Ny. By extracting the first and last terms in sums, we obtain(25)di02Un+1x0,yj+diNx2Un+1xNx,yj+d0j2Un+1xi,y0+dNyj2Un+1xi,yNy+∑k=1Nx-1dik2Un+1xk,yj+∑k=1Ny-1dkj2Un+1xi,yk+X1i,jnUn+1xi,yj-sαUn+1xi,yj=X3i,jns-sα-1fxi,yj,for i=1,…,Nx-1,j=1,…,Ny-1. We use the form of (25) to impose the boundary conditions.The solutionU~={U(x1,y1),Ux1,y2,…, U(xNx-1,yNy-1)}, which is the matrix of unknown interior points of U, can be found by solving the system(26)A~-sαIU~+X1n∘U~+U~B~=X3ns-F~,where A~ is the matrix of interior points of D^x(2) and B~ is the matrix of interior points of D^y(2), so that A~ and B~ match the dimensions of U~. We also use ∘ to denote the Hadamard product between two matrices. Also(27)F~ij=sα-1fxi,yj+di02Ux0,yj+diNx2UxNx,yj+d0j2Uxi,y0+dNyj2Uxi,yNy,for i=1,…,Nx-1,j=1,…,Ny-1.


#### 3.1.1. Dirichlet Boundary Conditions

Boundary conditions may be in the form of Dirichlet conditions,(28)u-1,y,t=a,u1,y,t=b,(29)ux,-1,t=c,ux,1,t=d,and hence,(30)U-1,y=La,U1,y=Lb,(31)Ux,-1=Lc,Ux,1=Ld.The parameters a, b, c, and d are potentially functions of the temporal variable and one of the spatial variables; i.e., a=a(y,t). We assume that a, b, c, and d are constant.Dirichlet boundary conditions can be imposed directly by substituting (28) and (29) into (25) and collecting all the known terms in F~.


#### 3.1.2. Neumann Boundary Conditions

Alternatively Neumann boundary conditions give(32)ux-1,y,t=a,ux1,y,t=b,(33)uyx,-1,t=c,uyx,1,t=d,with,(34)Ux-1,y=La,Uy1,y=Lb,(35)Uyx,-1=Lc,Uyx,1=Ld.Neumann boundary conditions given by (30) are discretized as (36)∂U∂xxNx=-1,yj≈∑k=0NxdNxkUnxk,yj=La,(37)∂U∂xx0=1,yj≈∑k=0Nxd0kUnxk,yj=Lb.Similarly for equations (31)(38)∂U∂yxi,yNy=-1≈∑k=0NydkNyUnxi,yk=Lc,(39)∂U∂yxi,y0=1≈∑k=0Nydk0Unxi,yk=Ld.By extracting the first and last terms in the sum, the discretizations can be written as(40)dNx0dNxNxd00d0NxUnx0,yjUnxNx,yj=Lb-∑k=1Nx-1d0kUnxk,yjLa-∑k=1Nx-1dNxkUnxk,yjand(41)d0NydNyNyd00dNy0Unxi,y0Unxi,yNy=Ld-∑k=1Ny-1dk0Unxi,ykLc-∑k=1Ny-1dkNyUnxi,ykand the solutions to these linear systems are then substituted into equation (25).


## 3.1. Chebyshev Collocation

Chebyshev polynomials form a basis on[-1,1] and hence we dictate the domain of our PDE to be Ω. We note here, however, that any domain in R2 can be trivially deformed to match Ω. We discretize our spatial domain using Chebyshev-Gauss-Lobatto points:(19)xi=cos⁡iπNx,i=0,1,…,Nx,yj=cos⁡jπNy,j=0,1,…,Ny.Given this choice of spatial discretization, we have x0=1, xNx=-1, y0=1, and yNy=-1, 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 Nx=Ny; i.e., we have equal number of collocation points in each spatial direction. We now define a differentiation matrix D(1)=dkl:(20)dkl=ck-1k+lclxk-xl,k≠l,-xk21-x2k,k=l,162Nx2+1,k=l=0,-162Nx2+1,k=l=Nx,whereck=2,k=0,Nx,1,k=1,…,Nx-1.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 D(2)=D(1).D(1), 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 thex direction are(21)D^x1=D1,D^x2=D2,where D(1) is the Chebyshev differentiation matrix of size (Nx+1)×(Ny+1).Because we have assumedNx=Ny, we derive the pleasing property that D^y(1)=D^x(1)T and D^y(2)=D^x(2)T.Writing the discretization of (17) in matrix form yields(22)D^x2Un+1+Un+1D^y2+X1nUn+1-sαUn+1=X3ns-F,where(23)Fij=sα-1fxi,yj.By expanding (22) in summation notation, we have(24)∑k=0Nxdik2Un+1xk,yj+∑k=0Nydkj2Un+1xi,yk+X1i,jnUn+1xi,yj-sαUn+1xi,yj=X3i,jns-sα-1fxi,yj,for i=0,1,…,Nx,j=0,1,…,Ny. By extracting the first and last terms in sums, we obtain(25)di02Un+1x0,yj+diNx2Un+1xNx,yj+d0j2Un+1xi,y0+dNyj2Un+1xi,yNy+∑k=1Nx-1dik2Un+1xk,yj+∑k=1Ny-1dkj2Un+1xi,yk+X1i,jnUn+1xi,yj-sαUn+1xi,yj=X3i,jns-sα-1fxi,yj,for i=1,…,Nx-1,j=1,…,Ny-1. We use the form of (25) to impose the boundary conditions.The solutionU~={U(x1,y1),Ux1,y2,…, U(xNx-1,yNy-1)}, which is the matrix of unknown interior points of U, can be found by solving the system(26)A~-sαIU~+X1n∘U~+U~B~=X3ns-F~,where A~ is the matrix of interior points of D^x(2) and B~ is the matrix of interior points of D^y(2), so that A~ and B~ match the dimensions of U~. We also use ∘ to denote the Hadamard product between two matrices. Also(27)F~ij=sα-1fxi,yj+di02Ux0,yj+diNx2UxNx,yj+d0j2Uxi,y0+dNyj2Uxi,yNy,for i=1,…,Nx-1,j=1,…,Ny-1.


### 3.1.1. Dirichlet Boundary Conditions

Boundary conditions may be in the form of Dirichlet conditions,(28)u-1,y,t=a,u1,y,t=b,(29)ux,-1,t=c,ux,1,t=d,and hence,(30)U-1,y=La,U1,y=Lb,(31)Ux,-1=Lc,Ux,1=Ld.The parameters a, b, c, and d are potentially functions of the temporal variable and one of the spatial variables; i.e., a=a(y,t). We assume that a, b, c, and d are constant.Dirichlet boundary conditions can be imposed directly by substituting (28) and (29) into (25) and collecting all the known terms in F~.


### 3.1.2. Neumann Boundary Conditions

Alternatively Neumann boundary conditions give(32)ux-1,y,t=a,ux1,y,t=b,(33)uyx,-1,t=c,uyx,1,t=d,with,(34)Ux-1,y=La,Uy1,y=Lb,(35)Uyx,-1=Lc,Uyx,1=Ld.Neumann boundary conditions given by (30) are discretized as (36)∂U∂xxNx=-1,yj≈∑k=0NxdNxkUnxk,yj=La,(37)∂U∂xx0=1,yj≈∑k=0Nxd0kUnxk,yj=Lb.Similarly for equations (31)(38)∂U∂yxi,yNy=-1≈∑k=0NydkNyUnxi,yk=Lc,(39)∂U∂yxi,y0=1≈∑k=0Nydk0Unxi,yk=Ld.By extracting the first and last terms in the sum, the discretizations can be written as(40)dNx0dNxNxd00d0NxUnx0,yjUnxNx,yj=Lb-∑k=1Nx-1d0kUnxk,yjLa-∑k=1Nx-1dNxkUnxk,yjand(41)d0NydNyNyd00dNy0Unxi,y0Unxi,yNy=Ld-∑k=1Ny-1dk0Unxi,ykLc-∑k=1Ny-1dkNyUnxi,ykand the solutions to these linear systems are then substituted into equation (25).


## 3.1.1. Dirichlet Boundary Conditions

Boundary conditions may be in the form of Dirichlet conditions,(28)u-1,y,t=a,u1,y,t=b,(29)ux,-1,t=c,ux,1,t=d,and hence,(30)U-1,y=La,U1,y=Lb,(31)Ux,-1=Lc,Ux,1=Ld.The parameters a, b, c, and d are potentially functions of the temporal variable and one of the spatial variables; i.e., a=a(y,t). We assume that a, b, c, and d are constant.Dirichlet boundary conditions can be imposed directly by substituting (28) and (29) into (25) and collecting all the known terms in F~.


## 3.1.2. Neumann Boundary Conditions

Alternatively Neumann boundary conditions give(32)ux-1,y,t=a,ux1,y,t=b,(33)uyx,-1,t=c,uyx,1,t=d,with,(34)Ux-1,y=La,Uy1,y=Lb,(35)Uyx,-1=Lc,Uyx,1=Ld.Neumann boundary conditions given by (30) are discretized as (36)∂U∂xxNx=-1,yj≈∑k=0NxdNxkUnxk,yj=La,(37)∂U∂xx0=1,yj≈∑k=0Nxd0kUnxk,yj=Lb.Similarly for equations (31)(38)∂U∂yxi,yNy=-1≈∑k=0NydkNyUnxi,yk=Lc,(39)∂U∂yxi,y0=1≈∑k=0Nydk0Unxi,yk=Ld.By extracting the first and last terms in the sum, the discretizations can be written as(40)dNx0dNxNxd00d0NxUnx0,yjUnxNx,yj=Lb-∑k=1Nx-1d0kUnxk,yjLa-∑k=1Nx-1dNxkUnxk,yjand(41)d0NydNyNyd00dNy0Unxi,y0Unxi,yNy=Ld-∑k=1Ny-1dk0Unxi,ykLc-∑k=1Ny-1dkNyUnxi,ykand 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 is(42)ux,0=e-x2,with Dirichlet boundary conditions consistent with this initial condition,(43)u-1,t=u1,t=e-1.Table 1illustrates the maximum absolute error between the solution obtained by the present method and the solution obtained by NDSolve in Mathematica 9 for α=1. To the best of the authors’ knowledge, no exact solution exists for the fractional case and, hence, no comparison can be made.Table 1
Maximum absolute error in the presented method’s solution of the problem described by Example 1 at timet=0.5.
Nx Error vsNDSolve 5 0.029145584018683945 10 0.00015607552473428932 15 0.00014992336248420557 20 0.00014949081096204964 25 0.0001494889958276735 30 0.00014948859870700382However, we present Figure1 which illustrates the behaviour of the mid-point of the solution as time evolves for various values of α.Figure 1
Plot ofu(x=0,t) for various values of α.


### 4.2. Example 2

We now consider the one-dimensional case with initial condition(44)ux,0=e-x+1and Neumann Boundary conditions(45)ux-1,t=-1(46)ux1,t=-1e2.Once again our solution converges with N moving from 5 to 10 but does not continue to improve. Table 2 presents the errors in our method when compared to NDSolve for α=1. We present the behaviour of the mid-point u(x=0,t) as time evolves for different values of α in Figure 2.Table 2
Maximum absolute error in the presented method’s solution of the problem described by Example 2 at timet=0.5.
Nx Error vsNDSolve 5 0.18888552094993993 10 0.06402298769357995 15 0.06716933046894086 20 0.06886100702660541 25 0.0685232103205754 30 0.06782287252683938Figure 2
Plot ofu(x=0,t) for various values of α.


### 4.3. Example 3

This example considers the two-dimensional case with initial condition(47)ux,y,0=x-1x+1y-1y+1and consistent Dirichlet boundary conditions described as(48)u-1,y,t=u1,y,t=ux,-1,t=ux,1,t=0.Table 3 presents the errors in the present method when compared to NDSolve. Figure 3 describes the evolution of the mid-point in two dimensions with time for various α values.Table 3
Maximum absolute error in the presented method’s solution of the problem described by Example 3 at timet=0.5.
Nx=Ny Error vsNDSolve 5 0.000050581035499991776 10 0.00006478718535577835 15 0.00006310249073780343 20 0.00006478739642413641 25 0.0000642415234639226 30 0.00006478739639844785Figure 3
Plot ofu(x=0,y=0,t) for various values of α.


### 4.4. Example 4

Finally we consider the case of a two-dimensional initial condition with Neumann conditions,(49)ux,y,0=cos⁡xcos⁡yand(50)ux-1,y,t=sin⁡1cos⁡y(51)ux1,y,t=-sin⁡1cos⁡y(52)uyx,-1,t=cos⁡xsin⁡1(53)uyx,1,t=-cos⁡xsin⁡1.The errors for this example when compared to NDSolve are presented in Table 4. Figure 4 describes the evolution of the mid-point in two dimensions with time for various α values.Table 4
Maximum absolute error in the presented method’s solution of the problem described by Example 4 at timet=1.
Nx=Ny Error vsNDSolve 5 0.00875011619006949 10 0.002380333748528196 15 0.002658963151695559 20 0.0027658847054590208 25 0.002819390822498269 30 0.0028497574705675377Figure 4
Plot ofu(x=0,y=0,t) for various values of α.


## 4.1. Example 1

The first example we consider is(42)ux,0=e-x2,with Dirichlet boundary conditions consistent with this initial condition,(43)u-1,t=u1,t=e-1.Table 1illustrates the maximum absolute error between the solution obtained by the present method and the solution obtained by NDSolve in Mathematica 9 for α=1. To the best of the authors’ knowledge, no exact solution exists for the fractional case and, hence, no comparison can be made.Table 1
Maximum absolute error in the presented method’s solution of the problem described by Example 1 at timet=0.5.
Nx Error vsNDSolve 5 0.029145584018683945 10 0.00015607552473428932 15 0.00014992336248420557 20 0.00014949081096204964 25 0.0001494889958276735 30 0.00014948859870700382However, we present Figure1 which illustrates the behaviour of the mid-point of the solution as time evolves for various values of α.Figure 1
Plot ofu(x=0,t) for various values of α.


## 4.2. Example 2

We now consider the one-dimensional case with initial condition(44)ux,0=e-x+1and Neumann Boundary conditions(45)ux-1,t=-1(46)ux1,t=-1e2.Once again our solution converges with N moving from 5 to 10 but does not continue to improve. Table 2 presents the errors in our method when compared to NDSolve for α=1. We present the behaviour of the mid-point u(x=0,t) as time evolves for different values of α in Figure 2.Table 2
Maximum absolute error in the presented method’s solution of the problem described by Example 2 at timet=0.5.
Nx Error vsNDSolve 5 0.18888552094993993 10 0.06402298769357995 15 0.06716933046894086 20 0.06886100702660541 25 0.0685232103205754 30 0.06782287252683938Figure 2
Plot ofu(x=0,t) for various values of α.


## 4.3. Example 3

This example considers the two-dimensional case with initial condition(47)ux,y,0=x-1x+1y-1y+1and consistent Dirichlet boundary conditions described as(48)u-1,y,t=u1,y,t=ux,-1,t=ux,1,t=0.Table 3 presents the errors in the present method when compared to NDSolve. Figure 3 describes the evolution of the mid-point in two dimensions with time for various α values.Table 3
Maximum absolute error in the presented method’s solution of the problem described by Example 3 at timet=0.5.
Nx=Ny Error vsNDSolve 5 0.000050581035499991776 10 0.00006478718535577835 15 0.00006310249073780343 20 0.00006478739642413641 25 0.0000642415234639226 30 0.00006478739639844785Figure 3
Plot ofu(x=0,y=0,t) for various values of α.


## 4.4. Example 4

Finally we consider the case of a two-dimensional initial condition with Neumann conditions,(49)ux,y,0=cos⁡xcos⁡yand(50)ux-1,y,t=sin⁡1cos⁡y(51)ux1,y,t=-sin⁡1cos⁡y(52)uyx,-1,t=cos⁡xsin⁡1(53)uyx,1,t=-cos⁡xsin⁡1.The errors for this example when compared to NDSolve are presented in Table 4. Figure 4 describes the evolution of the mid-point in two dimensions with time for various α values.Table 4
Maximum absolute error in the presented method’s solution of the problem described by Example 4 at timet=1.
Nx=Ny Error vsNDSolve 5 0.00875011619006949 10 0.002380333748528196 15 0.002658963151695559 20 0.0027658847054590208 25 0.002819390822498269 30 0.0028497574705675377Figure 4
Plot ofu(x=0,y=0,t) for various values of α.


## 5. Image Processing Application

Despite the impressive accuracy obtained by Chebyshev collocation for small values ofNx, the method suffers from severe round-off error for values of Nx>100 due to finite-precision arithmetic [22, 23]. In applying this scheme to an image, the dimensions of the input image dictate the resolution of the scheme, and in most cases the input image dimensions will exceed 100. This forces us to employ a finite-difference discretization scheme instead of the Chebyshev scheme as described by Jacobs and Harley in [2]. Jacobs and Momoniat [15] present the results of applying the integer order model to an image for the purpose of document image binarization. We present here some examples of the results obtained using the hybrid-transform method. The power of the hybrid-transform method lies in the ability to immediately determine the solution at any point in time as opposed to needing to iterate to a given point in time. Figure 5 illustrates the solutions obtained, given an input image, at varying points in time. Figure 6 shows the effects of the fractional term on the resulting image, wherein the steady state has been altered as well as the transient phase being executed much more quickly. This mirrors the effects illustrated in the examples above, Figures 1, 2, 3, and 4, where the steady state of the equation is altered with changing values of α as well as the equation reaching a steady state more rapidly. This suggests the order of the derivative may be used as an extra dimension of control in image processing, where a smaller α value actually brightens the background of the image; it also reduces diffusivity of the model by altering the steady state away from a uniformly diffused image.Figure 5
Figures depicting results obtained at increasing time using the current method andα=1. (a) The original image, (b) binarized original image at t≈0, (c) processed binary image at t=0.003, (d) processed binary image at t=0.008, (e) overdiffused image at t=0.01.
(a) (b) (c) (d) (e)Figure 6
Figures depicting results obtained at increasing time using the current method andα=0.2. (a) The original image, (b) processed image at t=0.01, (c) processed binary image at t=0.001, (d) processed binary image at t=0.01.
(a) (b) (c) (d)If we take, for example,Nx=Ny=100 then the stability criteria for a standard time-marching scheme would be △t=1/1002/α. As α decreases from 1, △t becomes so small that the number of iterations required to attain a final time of t=0.003 is unfeasibly large. This further justifies the use of a hybrid-transform method in the application of image processing.


### 5.1. Application Example

We present here a document image binarization application, where the input image is corrupted by noise as well as ink bleed-through from the reverse side of the page. The solution image we seek has black pixels which represent text and white pixels everywhere else. Figure7 illustrates the input image and the resulting binary images for different values of α as well as the ground truth image against which we may quantitatively measure our result. The L1 norm is used as an error measure and is presented at various points in time in Figure 8. Due to the changing α values we must normalize time, to accurately compare the performance of the present method over various α values. Due to the accelerated convergence to the steady state with varying α, as illustrated in the examples in Section 4, we note that α=0.5 results in an optimally processed image, supporting the application of fractional partial differential equations in image processing.Figure 7
Figures illustrating processed binary images obtained for various values ofα. (a) The original image, (b) α=1, (c) α=0.75, (d) α=0.5, (e) α=0.25, (f) ground truth image.
(a) (b) (c) (d) (e) (f)Figure 8
PlotL1 errors against number of iterations for various values of α.


## 5.1. Application Example

We present here a document image binarization application, where the input image is corrupted by noise as well as ink bleed-through from the reverse side of the page. The solution image we seek has black pixels which represent text and white pixels everywhere else. Figure7 illustrates the input image and the resulting binary images for different values of α as well as the ground truth image against which we may quantitatively measure our result. The L1 norm is used as an error measure and is presented at various points in time in Figure 8. Due to the changing α values we must normalize time, to accurately compare the performance of the present method over various α values. Due to the accelerated convergence to the steady state with varying α, as illustrated in the examples in Section 4, we note that α=0.5 results in an optimally processed image, supporting the application of fractional partial differential equations in image processing.Figure 7
Figures illustrating processed binary images obtained for various values ofα. (a) The original image, (b) α=1, (c) α=0.75, (d) α=0.5, (e) α=0.25, (f) ground truth image.
(a) (b) (c) (d) (e) (f)Figure 8
PlotL1 errors against number of iterations for various values of α.


## 6. Discussion

The results above indicate a strong convergence to a solution with increasingN. The similar results obtained for N>5 are attributed to the comparison being drawn between two numerical methods which may be different from the true exact solution. Results obtained in [1] indicate that, when the present method is compared with an exact solution, the accuracy increases dramatically with increasing N. Despite this, the present method attains results similar to that of an industry standard, NDSolve, which is encouraging.Figures1, 2, 3, and 4 indicate the dynamic behaviour of the subdiffusive process. Interestingly the diffusive process becomes increasingly more aggressive as α decreases from 1 to 0; however, the final “steady state” achieved by these processes is typically less severe than the standard diffusion model. In image processing terms we could achieve a semidiffused result extremely quickly by choosing α to be small, rather than using α=1 and diffusing over a longer period.In both the one- and two-dimensional cases, the Neumann boundary conditioned problems result in an accuracy of two orders of magnitude worse than that of the Dirichlet cases. In the linear case, where an exact solution exists, the accuracy of the present method increases dramatically with increasingN [2]; however, due to the comparison being drawn with another numerical method, NDSolve, the free ends affect the comparative solution, rather than reinforcing a condition on the boundary as in the Dirichlet case.


## 7. Conclusions

In addition to obtaining a high accuracy, the present method is robust for0<α≤1 and is in fact exact in the transformation from a time-fractional partial differential equation to partial differential equation or differential equation depending on the dimensionality. The errors incurred are then discretization error and numerical inversion of the Bromwich integral which is O(10.2-N) [17].The application of this method to images introduces a computational hurdle since the resolution of the image dictates the spatial discretization resolution. After the transformation one must solve a linear system that is of the same dimension as the input image. This is exacerbated by the use of Chebyshev collocation, since the derivative matrices are full and not tridiagonal or banded as they are when using a standard finite-difference scheme, where the Thomas algorithm can be put to use as in [18–20]. In such a case, a numerical method should be implemented to improve the computational time required for large input data. The Chebyshev collocation method is also known to have large round-off errors for large problems, Nx>100, due to finite-precision [22, 23]. Alternatively, for large systems a finite-difference discretization may be used to speed up computational time at the cost of accuracy in the solution. We have shown, for example, the images obtained by employing the hybrid-transform method to illustrate that the method does produce the expected results in a real-world application.We have presented a method that is robust for a time-fractional diffusion equation with a nonlinear source term for one- and two-dimensional cases with both Dirichlet and Neumann boundary conditions. We have shown through numerical experiments that in the case ofα=1 our solution obtains a result similar to that of NDSolve in Mathematica 9 and extends trivially to fractional-order temporal derivative of order α∈(0,1]. The application to image processing substantiates this method for real-world problems, highlighting the effectiveness of the fractional ordered derivative as a further dimension of control which has a dramatic effect both on the transition time of the model and on the final state attained. Within the realm of image processing this opens new avenues of research allowing more control over physically derived processes to construct methods that are physically substantial.To the best of the authors’ knowledge this hybridization of the Laplace transform, Chebyshev collocation, and quasi-linearization scheme has not yet been applied to FPDEs in one or two dimensions. This collaboration yields a method that is accurate and robust for time-fractional derivatives.



---

*Source: 8924547-2018-09-17.xml*