1 Introduction
Diffusion is a physical process that refers to the movement of molecules from a region of high concentration to one of lower concentration. The phenomenon of classical diffusion is well described by the following relationship
(1) 
where is a mean squared displacement of the diffusing molecule in the course of time , is the diffusion coefficient. However, it turns out that in many phenomena occurring in nature, the relation (1) is not sufficient to describe them, which is confirmed by the results of experiments in [32, 27, 12, 15, 16]. A more general process of random motion of molecules is the socalled anomalous diffusion, which is characterized by the nonlinear relationship [19, 20]:
(2) 
where is the generalized diffusion coefficient.
Processes known as moving boundary problems or Stefan like problems (in connection with the early work of Joseph Stefan [28]) involve solution of the diffusion equation in a domain with a free boundary. They are governed by the formula (1). Many monographs have been devoted to the classical Stefan problem, among which we can mention: [8, 11, 25, 10, 21].
The fractional Stefan problem is a generalization of the classical one, in which both the motion of the molecules and the position of the moving boundary are governed by the low described by formula (2). Currently, many analytical solutions (also approximate solutions) of the fractional Stefan problem are known and presented in papers [18, 30, 29, 26, 22, 17, 23, 24, 2, 13] and they constitute the overwhelming majority of all solutions. Numerical solutions and methods dedicated to solutions are not so numerous [9, 31, 6, 7] and still require further research.
The numerical scheme discussed in this paper is an extension of the front fixing method developed in [6]. The most important element of the new method is the much more accurate integration algorithm for the fractional integrodifferential equation. The new approach uses the AlAlaoui operator, which significantly increases the accuracy of numerical integration. The optimal contribution of the fractional trapezoidal rule and the fractional rectangular rule was implemented using an artificial neural network.
2 Preliminaries
This section is devoted to integrals and derivatives of a noninteger order together with some of their properties. We also recall the twoparameter Wright function related to the theory of partial differential equations of fractional order. The definitions in the field of numerical methods will also be presented. At the beginning, let us focus our attention on the definitions of two fractional operators: the leftsided RiemannLiouville integral and the leftsided Caputo derivative.
Definition 1
The leftsided RiemannLiouville integral of order , denoted as , is given by the following formula for :
(3) 
where is the Euler gamma function.
Definition 2
Let . The leftsided Caputo derivative of order is given by the formula:
(4) 
Two of the well known properties of integerorder integral and differential operators are preserved by the following generalizations:
Property 3 (cf. Lemma 2.3 [14])
If , and , then the equation
(5) 
is satisfied at almost every point for where . If , then the above relation holds at any point of .
Property 4 (cf. Lemma 2.22, [14])
Let function . Then, the composition rule for the leftsided RiemannLiouville integral and the leftsided Caputo derivative is given as follows:
(6) 
One of the very commonly used special functions in the theory of partial differential equations of fractional order is twoparameter Wright function:
Definition 5
Let , . The twoparameter Wright function is given as the following series:
(7) 
which is a generalization of the complementary error function.
The numerical scheme proposed in the further part of the paper uses a mesh of nodes defined as follows:
Definition 6
Let be a continuous region of solutions for the partial differential equation and will be the end time. Then the set we call the rectangular regular mesh described by the set of nodes.
The next three definitions will allow us to formulate the artificial neural network in the form of a very simple formula.
Definition 7
Function such that
where
is a vector function, that assigns to each point
a ndimensional vector lying in space.Definition 8
Function such that
where is a vector function, that assigns to each point a ndimensional vector lying in space.
Definition 9
Function such that
where is a vector function, that assigns to each point a n+1dimensional vector lying in space.
3 Mathematical formulation of the problem
The release of the solute is possible by anomalous diffusion through a penetrant fluid in the polymer matrix consists of two regions: first bounded by and where all solute is dissolved and second one bounded by and which contains undissolved solute. The boundary of the domains formed by the diffusion front is described by the unknown function . On the moving boundary local mass conservation law is fulfilled and the concentration is constant and equal to drug solubility . We also assume that: generalized diffusivity of the drug in the matrix is constant; initial concentration of drug is greater than its solubility ; matrix is heterogeneous and nonswellable. The mathematical model described diffusion of the solute through the dissolved drug phase (illustrated by a simple scheme in Figure 1) is formulated by the subdiffusion equation:
(8) 
supplemented with the boundary conditions
(9) 
initial conditions
(10) 
and fractional Stefan condition
(11) 
Model formulated by equations (811) assume that drug release is into a perfect sink, with zero drug concentration i.e. . This condition is a mathematical idealisation that is well approximated if the release medium is exchanged sufficiently rapidly to keep sink conditions, or if the volume of the release medium is so large that drug concentration in the medium is negligible.
Many mathematical models formulated by partial differential equations use dimensionless variables in numerical work. This approach is particularly useful in problems with universal scope applying to different scales [4, 5]. There are some advantages in using dimensionless variables e.g. reduction of parameters in equations has a very positive effect on the time necessary to carry out calculations for large sets of parameters.
Let us note that the above mathematical model depends on four parameters, where constant generalized diffusion coefficient is measured in unit . Thus if we write the group , then we see that is a dimensionless variable. For subdiffusion in a plane sheet we introduce the following dimensionless variables:
which reduces the number of parameters to two. In our case, a whole set of solutions for different values of: generalized diffusion coefficient , initial concentration of drug , drug solubility and matrix thickness can be obtained from a single solution in dimensionless variables by simple scaling. So we can rewrite equations (811) in a dimensionless form:
(12) 
supplemented with the boundary conditions
(13) 
initial conditions
(14) 
and fractional Stefan condition
(15) 
where is fractional Stefan number.
From the point of view of further numerical results presented in this paper, the analytical solution of the onephase fractional Stefan problem obtained by Liu Junyi et al. [13] is very valuable, because it will allow validation of the method proposed in section four. Using the notation introduced in this paper, the closed analytical solution to the problem under consideration is expressed in terms of the Wright’s function:
(16) 
where is the coefficient of the power function describing the position of the diffusion front which can be determined from the transcendental equation
(17) 
4 Solution of the problem
The subject of our considerations is to find a numerical solution of the problem described by the equations (811). The proposed approach takes place in three stages:

we introduce to subdiffusion equation (12) new spatial variable, which immobilizes dissolution boundary for all . Next we convert the fractional partial differential equation to the equivalent fractional integrodifferential equation. We are using two methods of numerical integration: the fractional rectangle rule and the fractional trapezoidal rule to discretization of the integrodifferential equation. The share of the trapezoidal rule in the integration scheme is determined by the parameter , whose optimal value at this stage is unknown.

in the second stage, we determine the optimal value of the parameter. This task is implemented as follows: for fixed values of the model and grid parameters (the value of we get from the equation (17)) we minimize some functional expressing the sum of the absolute errors generated by the numerical scheme received in the previous point. The values obtained in this way and the associated values of the model and mesh parameters create training set for the artificial neural network. Next, we create the artificial neural network to predict the values of the parameter.

At this stage, the value of the parameter is already known and can be used in the scheme from the first point. To make a numerical scheme independent of the analytical solution, we formulate an iterative algorithm that allows us to determine the value of the parameter regardless of the equation (17).
4.1 Front fixing method
The basic idea of the front fixing method is a suitable choice of a new space variable for the governing equation. The transformation:
(18) 
fixes dissolution boundary at the point for all . By using the following relationships
(19) 
(20) 
(21) 
the subdiffusion equation (12) can be written
(22) 
Now, we introduce the auxiliary function
(23) 
So, finally, we get the integrodifferential equation
(24) 
supplemented by the boundary conditions
(25) 
and initial condition
(26) 
It should be noted that due to further numerical considerations, we assume that is a very small positive real number.
Here we should explain our new approach. We can see that the integrodifferential equation discussed above, and in particular the last integral term on the right hand side of this equation can be discretized in different ways. In [6], the fractional trapezoidal rule was used to approximate the integral. However, it turns out that the use of the AlAlaoui operator ([1],[2],[3]) significantly improves the integration accuracy, which will be shown in examples later in the paper.
The second component on the right hand side of equation (24) can be approximated using the central differential quotient to discretize the derivative of function with respect to variable and the rectangle rule for calculating the integral:
(27) 
where
(28) 
Now let us apply the integration concept proposed by AlAlaoui. First, we are using the fractional trapezoidal rule for the leftsided RiemannLiouville integral and the differential quotient to discretize the second derivative
(29) 
where
(30) 
Then we calculate the same integral, but this time using the fractional rectangle rule
(31) 
where
(32) 
Finally, the approximation of the leftsided RiemannLiouville integral can be expressed in the following form
(33) 
Using formulas (2733) and (24), we get a numerical scheme that can be written in the matrix form
(34) 
where A is an matrix of coefficients of system of linear equations

B is a column vector with entires
and is a column vector with entires. We define the elements of matrix A and B as follows:
The formulas below allow us to recover the original spatial variable and function .
(35) 
(36) 
Particular attention should be given to the fact that abovederived numerical scheme depends on two unknown parameters: and . In the next two subsections we will demonstrate how to determine their values.
4.2 Prediction of the parameter value using an artificial neural network
Simultaneous determination of the optimum value of the parameter and is very problematic. Therefore, this task will be divided into two stages. First, we determine the optimal value. We will use an artificial neural network for this purpose, because this mathematical construct can generate predictions for very complex problems. Let us assume that the parameter depends on three variables: , and (second grid parameter depends on the spatial step as follows in definition 6 adopted
). The proposed network is a fourlayer, unidirectional and a full neural network. The first network layer (input layer) consists of three neurons loading variables:
and bias neuron with a fixed value of 1. The first hidden layer of the network contains five neurons applyingas an activation function and bias neuron with a fixed value of 1. The next hidden layer also contains five neurons using a sigmoid function as an activation function and bias neuron with a fixed value of 1. The output layer contains one neuron applying sigmoid function as an activation function. Two criteria were used to select the number and type of neurons: the accuracy of the results generated by the network and it’s as simple structure as possible. As we know, a large number of neurons in the hidden layer increases the computing power and flexibility of the neural network while learning complicated patterns, but also contributes to good matching of the network to the training set and loss of generalizing ability on the test set. The learning process (adjustment of the weights to every connection between neurons) of the neural network was carried out using a backpropagation algorithm.
The neural network learning process requires a training set, which is a group of sample inputs we can feed into the neural network in order to train the model. The training set used in this paper was created as follows. We first determine the possible values of loading variables: :
Next, we solve a socalled inverse problem for set . In an inverse problem, some information is not known, in our case we do not know parameter . By minimizing the following functional
(37) 
we find the approximate value of the parameter By and we denoted numerical solution obtained in mesh node by the scheme proposed in subsection 4.1, where the value of is known obtained from formula (17) and closed analytical solution obtained from formulas (16), (17) respectively.
As a result of the neural network learning process, a set of weights was obtained that we wrote in the matrix form:
The matrix contains the weight values connecting the network input layer with the first hidden layer. The row number of matrix means the number of the neuron in the input layer, while the number of the matrix column means the number of the neuron in the first hidden layer. Matrixes and have an analogous interpretation as in the case of the matrix, but refer to a subsequent network layers. Using the above matrix notation, we can write our neural network in the form of the following function of three variables:
(38) 
whose values are illustrated in the Figure 2.
An important step in creating artificial neural networks is the process of validation, in which we use another group of sample inputs which were not included in training and are different from the samples in the training set. The 5th section of this paper is devoted to the validation process of the proposed artificial neural network.
4.3 Determination of the parameter
A boundary surface, on which dissolution occurs, moves across the slab separating a region of a dissolved drug from an undissolved core as in Figure 1. The position of the dissolution front is described by the function (dimensionless form) whose two values are known: at the beginning of the dissolution process and when dissolution is complete . Based on this simple observation, we formulate the convergence criterion
(39) 
where
(40) 
The proposed algorithm is a simple and more elegant version of the algorithm from [6], that can be written on the following block diagram
The block diagram presented in Figure 3 requires some additional detailed explanations, special attention should be given to the superscript of the coefficient . The difference in notation has its interpretation. The values: , , contribute to the formation of a new spatial variable (18) used in numerical scheme (34)they allow us to determine the end times , while the superscript means that the values: , , have been calculated using the formula (40). The iterative algorithm ends when and are almost identical, when condition (39) is fulfilled.
5 Numerical examples
The method proposed in the paper is a modification of the numerical method described in [6], therefore we adopted the same values of the model parameters for validation of the new version of the method. Two variants of the meshes were used in the calculations: grid with a division into and nodes. The smaller mesh was used to generate the results i Tables 2 and 3, while the mesh with higher nodes density allowed to generate the data illustrated in Figures 418. Let us pay attention to an important fact, in the validation process of the neural network
, we assumed sample imputes different from the samples in the training set.
In Table 1 we collected the values of the coefficient obtained from transcendental equation (17) using the Newton algorithm.
1/3  0.546438  0.598238  0.669592  0.77614 

2/3  0.736836  0.808016  0.90623  1.05134 
1  0.871649  0.956298  1.07232  1.24014 
Table 2 shows the values of the coefficient obtained by the old version of the method. A comparison of the values in Tables 2 and 3 with Table 1 clearly indicates that the scheme supported by the neural network gives much more accurate results.
1/3  0.584814  0.618749  0.674536  0.774755 

2/3  0.78916  0.822729  0.895849  1.050512 
1  0.928442  0.957251  1.041601  1.240087 
1/3  0.543311  0.597387  0.660986  0.774755 

2/3  0.730993  0.803931  0.892553  1.050512 
1  0.862646  0.948828  1.041601  1.240087 
In the further part of the paper, we will show detailed results obtained for and , because according to Figure 2 and data collected in the tables above, the new numerical scheme has a large advantage over the old one, especially for small values of and .
In Figure 4, 5 and 6, we presented the drug concentration profiles obtained by the old version of the front fixing method. In all cases, the graphs were plotted for the time when the solid core of the drug was completely dissolved.
Similar results were presented in graph 7, 8 and 9, but were received by the method supported by an artificial neural network.
A boundary surface on which dissolution occurs, moves according to the graphs shown in Figure 10, 11 and 12. An analytical solution of the function has been marked with a black line. The green and red lines represent the numerical solution obtained with the old version of the front fixing method and the new version respectively. Graph analysis in all three cases clearly indicates the advantage of the method supported by an artificial neural network.
Figure 13, 14 and 15 shows the distribution of absolute errors in the domain generated by the old version of the front fixing method.
In the case of graphs in Figure 16, 17 and 18 we notice that the new method generates much smaller absolute errors. The analysis of Figures 1318 also leads to an important observation. In all considered cases, we note the highest absolute error values for small values of variable .
The essence of the front fixing method is the use of a new spatial variable (18) which has singularity at and is not defined at this point. Using a numerical approach consisting in taking as a starting point a very small positive number we generate some error at the beginning of calculations. This initial error is transferred to subsequent time layers of the mesh, indirectly by the preceding time layers, and directly through the first time layer in proportion to a certain weight. It seems that the numerical method supported by an artificial neural network allows us to limit the adverse effect of the results from the initial time layers.
6 Conclusions
In this paper we proposed a new version of the front fixing method supported by an artificial neural network. The new scheme allows us to obtain accurate results especially in cases where the previous version ot the method generated relatively large errors, i.e. for
. It should be noted that the neural network presented in the paper is only a simple example of the application of artificial intelligence in the optimization process of numerical integration. The use of a more complex network could hypothetically improve results, but at the expense of simplicity and ease of use. The presented approach can also be extended to the twophase fractional LaméClapeyronStefan problem.
References
 [1] (1993) Novel digital integrator and differentiator. Electronics Letters 29 (4), pp. 376–378. Cited by: §4.1.
 [2] (2006) Alalaoui operator and the approximation for discretization of analog systems. Facta Universitatis Ser.: Elec. Energ. 19 (1), pp. 143–146. Cited by: §1, §4.1.
 [3] (2008) Alalaoui operator and the new transformation polynomials for discretization of analogue systems. Electrical Engining 90 (6), pp. 455–467. Cited by: §4.1.
 [4] (1996) Scaling, selfsimilarity, and intermediate asymptotics. Cambridge University Press, Cambridge. Cited by: §3.
 [5] (2003) Scaling. Cambridge University Press, Cambridge. Cited by: §3.
 [6] (2015) Numerical solution of the one phase 1d fractional Stefan problem using the front fixing method. Mathematical Methods in the Applied Sciences 38 (15), pp. 3214–3228. Cited by: §1, §1, §4.1, §4.3, §5.
 [7] (2014) Numerical scheme for onephase 1d fractional stefan problem using the similarity variable technique. Journal of Applied Mathematics and Computational Mechanics 13 (1), pp. 13–21. Cited by: §1.
 [8] (1984) Free and moving boundary problems. Clarendon Press, Oxford. Cited by: §1.
 [9] (2015) The numerical method for the moving boundary problem with spacefractional derivative in drug release devices. Applied Mathematical Modelling 39, pp. 2385–2391. Cited by: §1.
 [10] (2003) The classical stefan problem. basic concepts, modeling and analysis. Elsevier, Amsterdam. Cited by: §1.
 [11] (1987) Onedimensional stefan problems: an introduction. Longman Scientific and Technical, New York. Cited by: §1.
 [12] (2010) Environmental context explains Lévy and Brownian movement patterns of marine predators. Nature 465, pp. 1066–1069. Cited by: §1.
 [13] (2009) Some exact solutions to stefan problems with fractional differential equations. Journal of Mathematical Analysis and Applications 351, pp. 536–542. Cited by: §1, §3.
 [14] (2006) Theory and applications of fractional differential equations. Elsevier, Amsterdam. Cited by: Property 3, Property 4.
 [15] (2005) How to measure subdiffusion parameters. Physical Review Letters 94, pp. 170602. Cited by: §1.
 [16] (2005) Measuring subdiffusion parameters. Physical Review E 71, pp. 041105. Cited by: §1.
 [17] (2008) Scaleinvariant solutions to partial differential equations of fractional order with a moving boundary condition. Journal of Physics A: Mathematical and Theoretical 41, pp. 155202. Cited by: §1.
 [18] (2004) An exact solution to the moving boundary problem with fractional anomalous diffusion in drug release devices. Zeitschrift fur Angewandte Mathematik und Mechanik 84, pp. 22–28. Cited by: §1.
 [19] (2000) The random walk:s guide to anomalous diffusion: a fractional dynamics approach. Physics Reports 339, pp. 1–77. Cited by: §1.
 [20] (2004) The restaurant at the end of the random walk: recent developments in the description of anomalous transport by fractional dynamics. Journal of Physics A: Mathematical and General 37, pp. 161–208. Cited by: §1.
 [21] (1993) Heat conduction. 2 edition, Wiley. Cited by: §1.
 [22] (2013) Homotopy perturbation method for a limit case stefan problem governed by fractional diffusion equation. Applied Mathematical Modelling 37, pp. 3589–3599. Cited by: §1.
 [23] (2013) Two equivalent stefan?s problems for the time fractional diffusion equation. Fractional Calculus and Applied Analysis 16 (4), pp. 802–815. Cited by: §1.
 [24] (2015) Hopf lemma for the fractional diffusion operator and its application to a fractional freeboundary problem. Journal of Mathematical Analysis and Applications 434 (1), pp. 125–135. Cited by: §1.
 [25] (1971) The stefan problem. American Mathematical Society. Cited by: §1.
 [26] (2011) Homotopy perturbation method to spacetime fractional solidification in a finite slab. Applied Mathematical Modelling 35, pp. 1937–1945. Cited by: §1.
 [27] (1993) Observations of anomalous diffusion and Lévy flights in a 2dimensional rotating flow. Physical Review Letters 71, pp. 3975–3979. Cited by: §1.
 [28] (1891) Uber die theorie der eisbildung, insbesondere uber die eisbildung im polarmeere. Annalen der Physik und Chemie 278 (2), pp. 269–286. Cited by: §1.
 [29] (2013) Two exact solutions of a stefan problem with varying diffusivity. International Journal of Heat and Mass Transfer 58, pp. 80–85. Cited by: §1.
 [30] (2010) An exact solution of a limit case stefan problem governed by a fractional diffusion equation. International Journal of Heat and Mass Transfer 53, pp. 5622–5625. Cited by: §1.
 [31] (2014) Fractional stefan problems. International Journal of Heat and Mass Transfer 74, pp. 269–277. Cited by: §1.
 [32] (1996) Anomalous diffusion in asymmetric random walks with a quasigeostrophic flow example. Physica D: Nonlinear Phenomena 97, pp. 291–310. Cited by: §1.
Comments
There are no comments yet.