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.
However, we present Figure 1 which illustrates the behaviour of the mid-point of the solution as time evolves for various values of .
4.2. Example 2
We now consider the one-dimensional case with initial conditionand Neumann Boundary conditionsOnce again our solution converges with moving from 5 to 10 but does not continue to improve. Table 2 presents the errors in our method when compared to NDSolve for . We present the behaviour of the mid-point as time evolves for different values of in Figure 2.
4.3. Example 3
This example considers the two-dimensional case with initial conditionand consistent Dirichlet boundary conditions described asTable 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.
4.4. Example 4
Finally we consider the case of a two-dimensional initial condition with Neumann conditions,andThe 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.
5. Image Processing Application
Despite the impressive accuracy obtained by Chebyshev collocation for small values of , the method suffers from severe round-off error for values of 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.
(a)
(b)
(c)
(d)
(e)
(a)
(b)
(c)
(d)
If we take, for example, then the stability criteria for a standard time-marching scheme would be . As decreases from 1, becomes so small that the number of iterations required to attain a final time of 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. Figure 7 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 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 results in an optimally processed image, supporting the application of fractional partial differential equations in image processing.
(a)
(b)
(c)
(d)
(e)
(f)
6. Discussion
The results above indicate a strong convergence to a solution with increasing . The similar results obtained for 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 . Despite this, the present method attains results similar to that of an industry standard, NDSolve, which is encouraging.
Figures 1, 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 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 increasing [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 for 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 [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, , 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 our solution obtains a result similar to that of NDSolve in Mathematica 9 and extends trivially to fractional-order temporal derivative of order . 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.
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Disclosure
Opinions expressed and conclusions arrived at are those of the authors and are not necessarily to be attributed to the CoE-MaSS. The work contained in this paper also forms part of B. A. Jacobs’ Doctoral thesis.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
B. A. Jacobs acknowledges the support of the NRF under Thuthuka grant 94005. The DST-NRF Centre of Excellence in Mathematical and Statistical Sciences (CoE-MaSS) is acknowledged for funding.