In this section, the Trotter OS scheme [53] is applied to Eq. (1.1) and converts it to some simpler PDEs. Then, the time dimension is discretized by a backward Euler approach. Finally, the obtained equations are solved by the CS–RBF method.
Temporal discretization
In this part, we apply a divide-and-conquer approach which is known as the OS method, to divide Eq. (1.1) into two subproblems in separate directions. Then, the Caputo time derivative is discretized in each subproblem.
Trotter splitting scheme
First, in order to reduce the complexity of problem (1.1), a first-order OS method called Trotter splitting approach is employed on Eq. (1.1) to solve FPMM along separate directions. Using this scheme reduces the size of obtained matrices and thus substantially reduces the computational complexity. Consider the time interval I which is divided into (M+1) equally spaced steps (t_{0},t_{2},dots ,t_{M}), (t_{k}=kDelta t), (k=0,1,dots ,M), and (Delta t=frac{T}{M}). Equation (1.1) can be written as follows:
$$ _{0}^{C}mathcal{D}_{t}^{alpha}u= mathcal{L}_{x}(u)+mathcal{L}_{y}(u), $$
(2.1)
where (mathcal{L}_{x}) and (mathcal{L}_{y}) are suboperators on function u as follows:
$$begin{aligned}& mathcal{L}_{x}(u)=frac{partial}{partial x}g bigl( Vert nabla G_{ sigma}*u Vert bigr)frac{partial}{partial x}u+g bigl( Vert nabla G_{sigma}*u Vert bigr)frac{partial ^{2} }{partial x^{2}}u, end{aligned}$$
(2.2)
$$begin{aligned}& mathcal{L}_{y}(u)=frac{partial}{partial y}g bigl( Vert nabla G_{ sigma}*u Vert bigr)frac{partial}{partial y}u+g bigl( Vert nabla G_{sigma}*u Vert bigr)frac{partial ^{2}}{partial y^{2}}u. end{aligned}$$
(2.3)
By employing the Trotter splitting method, the approximate solution of (u(x,y,t_{k+1})), which can simply be written as (u^{k+1}), is obtained as (u^{k+1}=[mathcal{S}_{x}^{Delta t} mathcal{S}_{y}^{Delta t}]u^{k}), in which (mathcal{S}_{x}^{Delta t}) and (mathcal{S}_{y}^{Delta t}) are ({}_{0}^{C}mathcal{D}_{t}^{alpha}u^{*}=mathcal{L}_{x}(u^{*})) and ({}_{0}^{C}mathcal{D}_{t}^{alpha}u^{**}=mathcal{L}_{y}(u^{**})), respectively. Finally, we have (u^{k+1}=u^{**^{k+1}}). The interested readers can see [43, 54–56] for more details about the OS methods.
Now, a proper scheme should be proposed to discretize the Caputo time derivative operator ({}_{0}^{C}mathcal{D}_{t}^{alpha}) in the subproblems.
Discretization of Caputo derivative
At each time step (t_{k}), (k=0,2,dots ,M), using the backward Euler approach, the time fractional derivative is approximated as follows:
$$ {}_{0}^{C}mathcal{D}_{t}^{alpha}u(x,y,t_{k+1})= frac{1}{Gamma (1-alpha )}sum_{j=0}^{k} int _{t_{j}}^{t_{j+1}} frac{1}{(t_{k+1}-tau )^{alpha}} frac{partial}{partial tau}u(x,y, tau ),dtau . $$
(2.4)
So, we have
$$ begin{aligned} &{}_{0}^{C} mathcal{D}_{t}^{alpha}u(x,y,t_{k+1}) \ &quad =frac{1}{Gamma (1-alpha )}sum_{j=0}^{k} frac{u(x,y,t_{j+1})-u(x,y,t_{j})}{Delta t} int _{t_{j}}^{t_{j+1}} frac{dtau}{(t_{k+1}-tau )^{alpha}}+R_{k+1}, end{aligned} $$
(2.5)
in which (R_{k+1}) is the truncation error and we have
$$ R_{k+1}leq mathcal{C} Biggl[frac{1}{Gamma (1-alpha )}sum _{j=0}^{k} int _{t_{j}}^{t_{j+1}} frac{t_{j+1}+t_{j}-2tau}{(t_{k+1}-tau )^{alpha}},dtau + mathcal{O} bigl(Delta t^{2} bigr) Biggr], $$
(2.6)
where (mathcal{C}) is a constant that depends only on u. Based on Lemma 3.1 in [57] and the fact that (Gamma (2-alpha )geq 2) for all (alpha in [0,1]), it can be concluded that
$$ Bigglvert frac{1}{Gamma (1-alpha )}sum_{j=0}^{k} int _{t_{j}}^{t_{j+1}} frac{t_{j+1}+t_{j}-2tau}{(t_{k+1}-tau )^{alpha}},dtau Biggrvert leq 2Delta t^{2-alpha}. $$
So
$$ vert R_{k+1} vert leq mathcal{C}Delta t^{2-alpha}. $$
(2.7)
Moreover, from Eq. (2.4), the following equations can be computed:
$$ begin{aligned} &frac{1}{Gamma (1-alpha )}sum _{j=0}^{k} int _{t_{j}}^{t_{j+1}} frac{1}{(t_{k+1}-tau )^{alpha}} frac{partial}{partial tau}u(x,y, tau ),dtau \ &quad =frac{1}{Gamma (1-alpha )}sum_{j=0}^{k} frac{u(x,y,t_{j+1})-u(x,y,t_{j})}{Delta t} int _{t_{k-j}}^{t_{k+1-j}} frac{dt}{t^{alpha}} \ &quad =frac{1}{Gamma (1-alpha )}sum_{j=0}^{k} frac{u(x,y,t_{k+1-j})-u(x,y,t_{k-j})}{Delta t} int _{t_{j}}^{t_{j+1}} frac{dt}{t^{alpha}} \ &quad =frac{1}{Gamma (2-alpha )}sum_{j=0}^{k} frac{u(x,y,t_{k+1-j})-u(x,y,t_{k-j})}{Delta t^{alpha}} bigl((j+1)^{1- alpha}-j^{1-alpha} bigr). end{aligned} $$
(2.8)
The discrete fractional differential operator (D^{alpha}_{t}) is as follows:
$$ D^{alpha}_{t}=frac{1}{Gamma (2-alpha )}sum _{j=0}^{k} frac{u(x,y,t_{k+1-j})-u(x,y,t_{k-j})}{Delta t^{alpha}} bigl((j+1)^{1- alpha}-j^{1-alpha} bigr). $$
(2.9)
Thus, we can write
$$ _{0}^{C}mathcal{D}_{t}^{alpha}u(x,y,t_{k+1})=D^{alpha}_{t}+R_{k+1}. $$
(2.10)
So (D^{alpha}_{t}) is an approximation to ({}_{0}^{C}mathcal{D}_{t}) and we have
$$ begin{aligned} &{}_{0}^{C} mathcal{D}_{t}^{alpha}u(x,y,t_{k+1}) \ &quad approx frac{1}{Gamma (2-alpha )(Delta t)^{alpha}}sum_{j=0}^{k} bigl((j+1)^{1-alpha}-j^{1-alpha} bigr) bigl(u(x,y,t_{k-j+1})-u(x,y,t_{k-j}) bigr). end{aligned} $$
(2.11)
For the sake of simplicity, we define (mathcal{C}_{Delta t}^{alpha}= frac{1}{Gamma (2-alpha )(Delta t)^{alpha}}), (mathcal{A}_{j}^{alpha}=(j+1)^{1-alpha}-j^{1-alpha}), and (u(x,y,t_{k})=u^{k}).
Remark 2.1
From Theorems 3.1 and 3.2 in [57], it can be deduced that the proposed temporal discretization scheme is unconditionally stable for (Delta t>0) and has globally ((2-alpha ))-order accuracy for (0<alpha <1).
Now, by using Eq. (2.11), the subproblems are discretized along the temporal dimension as follows:
$$begin{aligned} & begin{aligned} mathcal{C}_{Delta t}^{alpha} Biggl[u^{*^{k+1}}-u^{k}+sum_{j=1}^{k} bigl(u^{{k+1-j}}-u^{{k-j}} bigr)mathcal{A}_{j}^{alpha} Biggr]&= frac{partial}{partial x}g bigl( biglVert nabla G_{sigma}*u^{k} bigrVert bigr) frac{partial}{partial x}u^{*^{k+1}} \ &quad {} +g bigl( biglVert nabla G_{sigma}*u^{k} bigrVert bigr) frac{partial ^{2}}{partial x^{2}}u^{*{k+1}}, end{aligned} end{aligned}$$
(2.12)
$$begin{aligned} &begin{aligned} mathcal{C}_{Delta t}^{alpha} Biggl[u^{**^{k+1}}-u^{*^{k}}+sum_{j=1}^{k} bigl(u^{*^{k+1-j}}-u^{*^{k-j}} bigr)mathcal{A}_{j}^{alpha} Biggr]&= frac{partial}{partial y}g bigl( biglVert nabla G_{sigma}*u^{*^{k}} bigrVert bigr) frac{partial}{partial y}u^{**^{k+1}} \ &quad {} +g bigl( biglVert nabla G_{sigma}*u^{*^{k}} bigrVert bigr) frac{partial ^{2}}{partial y^{2}}u^{**{k+1}}. end{aligned} end{aligned}$$
(2.13)
Spatial discretization
Before going to the space discretization of problems (2.12) and (2.13), we focus on the convolution term (G_{sigma}*u^{k} (G_{sigma}*u^{*^{k}})). Since the kernel (G_{sigma}) is smooth, the term (G_{sigma}*u^{k} (G_{sigma}*u^{*^{k}})) can be replaced by a heat equation with initial condition (u^{k}(u^{*^{k}})). These heat equations are solved by the OS method in conjunction with the Crank–Nicolson and CS–RBF approaches. The solutions of heat equations replaced in Eqs. (2.12) and (2.13) are represented by (u^{c}) and (u^{*^{c}}), respectively.
By substituting of the solution of the heat equations, Eqs. (2.12) and (2.13) can be considered in the following forms:
$$begin{aligned} &begin{aligned} &frac{partial}{partial x}g bigl( biglVert nabla u^{c} bigrVert bigr) frac{partial}{partial x}u^{*^{k+1}}+g bigl( biglVert nabla u^{c} bigrVert bigr) frac{partial ^{2}}{partial x^{2}}u^{*{k+1}}- mathcal{C}_{Delta t}^{ alpha}u^{*^{k+1}} \ &quad =mathcal{C}_{Delta t}^{alpha}sum_{j=1}^{k} bigl(u^{{k+1-j}}-u^{{k-j}} bigr)mathcal{A}_{j}^{alpha}- mathcal{C}_{Delta t}^{alpha}u^{k}, end{aligned} end{aligned}$$
(2.14)
$$begin{aligned} &begin{aligned} &frac{partial}{partial y}g bigl( biglVert nabla u^{*^{c}} bigrVert bigr) frac{partial}{partial y}u^{**^{k+1}}+g bigl( biglVert nabla u^{*^{c}} bigrVert bigr)frac{partial ^{2}}{partial y^{2}}u^{**{k+1}}- mathcal{C}_{ Delta t}^{alpha}u^{**^{k+1}} \ &quad =mathcal{C}_{Delta t}^{alpha}sum_{j=1}^{k} bigl(u^{*^{k+1-j}}-u^{*^{k-j}} bigr)mathcal{A}_{j}^{alpha}- mathcal{C}_{Delta t}^{alpha}u^{*^{k}}. end{aligned} end{aligned}$$
(2.15)
By obtaining (u^{c}) and (u^{*^{c}}), functions (g(Vert nabla u^{c}Vert )) and (g(Vert nabla u^{*^{c}}Vert )) are computed as follows:
$$ g bigl( biglVert nabla u^{c} bigrVert bigr)= frac{1}{1+frac{(frac{partial u^{c}}{partial x})^{2}+(frac{partial u^{c}}{partial y})^{2}}{K^{2}}}, qquad g bigl( biglVert nabla u^{*^{c}} bigrVert bigr)= frac{1}{1+frac{(frac{partial u^{*^{c}}}{partial x})^{2}+(frac{partial u^{*^{c}}}{partial y})^{2}}{K^{2}}}.$$
Moreover, (frac{partial g(Vert nabla u^{c}Vert )}{partial x}) and (frac{partial g(Vert nabla u^{*^{c}}Vert )}{partial y}) can be obtained by the chain rule.
Now, CS–RBF method is applied on Eqs. (2.14) and (2.15) which has two main benefits. First, as there is no need for any mesh generation for constructing the shape functions, the CS–RBFs scheme is considered truly meshless, and the coefficient matrix of this approach is well-conditioned and can be solved by common methods easily. Moreover, the local RBFs like Wendland CS–RBFs do not have a shape parameter and are more stable than global ones.
At the first step of the space discretization, the nodal points over the domain Ω should be selected at each time step (t_{k}). We consider pixels of the images, i.e., (Xtimes Y) where (X=lbrace x_{1},dots ,x_{N_{x}}rbrace ) and (Y=lbrace y_{1},dots ,y_{N_{y}}rbrace ), as the mesh nodes on the boundary and inside of the domain. In order to approximate (u^{k}) for all (kgeq 0), (u^{0}) is computed by the interpolation, then (u^{k}) for (k>0) is obtained by solving Eq. (2.12) using the CS–RBFs scheme. In addition, CS–RBFs are symmetric with respect to their center points by the following definition:
Definition 2.2
([45])
A function (phi :mathbb{R}^{s}rightarrow mathbb{R}) is called radial provided there exists a univariate function (varphi :[0,infty )rightarrow mathbb{R}) such that
$$ phi (mathbf{x})=varphi (r),quad text{where } r= Vert mathbf{x} Vert $$
(2.16)
and (Vert cdot Vert ) is some norm on (mathbb{R}^{s}), usually the Euclidean norm.
Among different choices of existing radial basis functions, in this paper, Wendland’s compactly supported radial basis functions (WCS–RBFs) with (C^{2}) smoothness are selected because they do not include any shape parameters that make these functions easier to use. These functions are defined as follows:
$$ phi _{i}(x)=(1-delta r_{i})^{4}_{+}(1+4 delta r_{i}), $$
(2.17)
in which (r_{i}=Vert x-x_{i}Vert ) is the distance from node (x_{i}) to x, and δ is the size of support for the radial function (phi _{i}(x)). Moreover, function ((1-delta r_{i})^{4}_{+}) is defined such that for (0leq delta r_{i} <1) it is equal to ((1-delta r_{i})^{4}), otherwise is zero.
To approximate functions (u^{k}), (u^{*^{k+1}}), and (u^{**^{k+1}}) in Ω over center points ((x_{i},y_{j})), (i=1,2,dots ,N_{x}), (j=1,2,dots ,N_{y}), the CS–RBF interpolation approximations (hat{u}^{k}), (hat{u}^{*^{k+1}}), and (hat{u}^{**^{k+1}}) are defined as
$$begin{aligned}& hat{u}^{k}(x,y_{j}) =sum _{i=0}^{N_{x}}theta _{ij}^{k} phi _{i}(x),quad j= 1, 2, dots ,N_{y}, end{aligned}$$
(2.18)
$$begin{aligned}& hat{u}^{*^{k+1}}(x,y_{j}) =sum _{i=1}^{N_{x}}lambda _{ij}^{k+1} phi _{i}(x),quad j= 1, 2, dots ,N_{y}, end{aligned}$$
(2.19)
$$begin{aligned}& hat{u}^{**^{k+1}}(x_{i},y) =sum _{j=1}^{N_{y}}gamma _{ij}^{k+1} phi _{j}(y), quad i= 1, 2, dots ,N_{x}, end{aligned}$$
(2.20)
in which (phi _{i}) and (phi _{j}) are RBFs based on centers (lbrace x_{i}rbrace _{i=0}^{N_{x}}) and (lbrace y_{j}rbrace _{j=0}^{N_{y}}). Additionally, (lambda _{ij}) and (gamma _{ij}) are unknown coefficients which should be specified.
Finally, substituting the approximation function of Eq. (2.19) into Eq. (2.14) associated to the x-direction, for all nodes in Ω, the matrix forms of their discrete equations are obtained as follows:
$$ mathbf{A}_{j}Lambda ^{k+1}_{j}= mathbf{B}_{j},quad j=1, 2, dots , N_{y}, $$
(2.21)
where (Lambda ^{k+1}_{j}= [ lambda _{1j}^{k+1},lambda _{2j}^{k+1}, dots ,lambda _{N_{x}j}^{k+1} ]^{T}), (mathbf{A}_{j}) and (mathbf{B}_{j}) are (N_{x} times N_{x}) matrix and (N_{x} times 1) vector, respectively, defined as
$$begin{aligned}& mathbf{A}_{j} =mathbf{D} bigl(mathbf{G}_{x}^{j} Upsilon _{x}+ mathbf{G}^{j}Psi _{x} – mathcal{C}_{Delta t}^{alpha}Phi _{x} bigr)+ mathbf{H}_{x}, end{aligned}$$
(2.22)
$$begin{aligned}& mathbf{B}_{j} =mathbf{D} Biggl(mathcal{C}_{Delta t}^{alpha} sum_{p=1}^{k} mathcal{A}_{k}^{alpha} bigl(Phi _{x}Theta ^{p+1-k}_{j}-Phi _{x} Theta ^{p-k}_{j} bigr)- mathcal{C}_{Delta t}^{alpha}Phi _{x} Theta ^{k}_{j} Biggr), end{aligned}$$
(2.23)
where
$$begin{aligned}& Phi _{x} = begin{bmatrix} phi _{1}(x_{1}) & phi _{2}(x_{1}) & ldots & phi _{N_{x}-1}(x_{1})& phi _{N_{x}}(x_{1}) \ phi _{1}(x_{2})& phi _{2}(x_{2}) & ldots & phi _{N_{x}-1}(x_{2})& phi _{N_{x}}(x_{2}) \ vdots & vdots &ddots &vdots & vdots \ phi _{1}(x_{N_{x}-1})& phi _{2}(x_{N_{x}-1}) & ldots & phi _{N_{x}-1}(x_{N_{x}-1})& phi _{N_{x}}(x_{N_{x}-1}) \ phi _{1}(x_{N_{x}})& phi _{2}(x_{N_{x}}) & ldots & phi _{n}(x_{N_{x}})& phi _{N_{x}}(x_{N_{x}}) end{bmatrix}, \& Upsilon _{x} = begin{bmatrix} phi ‘_{1}(x_{1}) & phi ‘_{2}(x_{1}) & ldots & phi ‘_{N_{x}-1}(x_{1})& phi ‘_{N_{x}}(x_{1}) \ phi ‘_{1}(x_{2})& phi ‘_{2}(x_{2}) & ldots & phi ‘_{N_{x}-1}(x_{2})& phi ‘_{N_{x}}(x_{2}) \ vdots & vdots &ddots &vdots & vdots \ phi ‘_{1}(x_{N_{x}-1})& phi ‘_{2}(x_{N_{x}-1}) & ldots & phi ‘_{N_{x}-1}(x_{N_{x}-1})& phi ‘_{N_{x}}(x_{N_{x}-1}) \ phi ‘_{1}(x_{N_{x}})& phi ‘_{2}(x_{N_{x}}) & ldots & phi ‘_{n}(x_{N_{x}})& phi ‘_{N_{x}}(x_{N_{x}}) end{bmatrix}, \& Psi _{x} = begin{bmatrix} phi ”_{1}(x_{1}) & phi ”_{2}(x_{1}) & ldots & phi ”_{N_{x}-1}(x_{1})& phi ”_{N_{x}}(x_{1}) \ phi ”_{1}(x_{2})& phi ”_{2}(x_{2}) & ldots & phi ”_{N_{x}-1}(x_{2})& phi ”_{N_{x}}(x_{2}) \ vdots & vdots &ddots &vdots & vdots \ phi ”_{1}(x_{N_{x}-1})& phi ”_{2}(x_{N_{x}-1}) & ldots & phi ”_{N_{x}-1}(x_{N_{x}-1})& phi ”_{N_{x}}(x_{N_{x}-1}) \ phi ”_{1}(x_{N_{x}})& phi ”_{2}(x_{N_{x}}) & ldots & phi ”_{n}(x_{N_{x}})& phi ”_{N_{x}}(x_{N_{x}}) end{bmatrix}, \& mathbf{D}=operatorname{diag}(0,1,1,1,dots ,1,0), \& begin{aligned} mathbf{G}^{j}={}&operatorname{diag} bigl(g bigl( biglVert nabla u^{c}(x_{1},y_{j}) bigrVert bigr),g bigl( biglVert nabla u^{c}(x_{2},y_{j}) bigrVert bigr), \ & dots ,g bigl( biglVert nabla u^{c}(x_{N_{x}-1},y_{j}) bigrVert bigr),g bigl( biglVert nabla u^{c}(x_{N_{x}},y_{j}) bigrVert bigr) bigr), end{aligned} \& begin{aligned} mathbf{G}_{x}^{j}={}&operatorname{diag} biggl( frac{partial}{partial x}g bigl( biglVert nabla u^{c}(x_{1},y_{j}) bigrVert bigr),frac{partial}{partial x}g bigl( biglVert nabla u^{c}(x_{2},y_{j}) bigrVert bigr), \ & dots ,frac{partial}{partial x}g bigl( biglVert nabla u^{c}(x_{N_{x}-1},y_{j}) bigrVert bigr),frac{partial}{partial x}g bigl( biglVert nabla u^{c}(x_{N_{x}},y_{j}) bigrVert bigr) biggr), end{aligned} end{aligned}$$
(mathbf{H}_{x}) is an (N_{x}times N_{x}) matrix whose first and last rows are first and last rows of (Upsilon _{x}) and the other rows contain only zero elements. Moreover,
$$ Theta ^{m}_{j}= bigl(theta _{1j}^{m}, theta _{2j}^{m},dots ,theta _{N_{x}j}^{m} bigr)^{T}, $$
that are known coefficients computed in previous iterations. Indeed, there are (N_{y}) linear systems of (N_{x}) algebraic equations which are well-posed and can be solved by a standard approach such as LU factorization. By solving these systems, the coefficient vectors (Lambda _{j}^{k+1}) are obtained and the approximation of (hat{u}^{*^{k+1}}) is computed as follows:
$$ mathbf{U}^{*^{k+1}}=Phi bigl[Lambda _{1}^{k+1}, Lambda _{2}^{k+1}, dots ,Lambda _{N_{y}}^{k+1} bigr], $$
(2.24)
in which
$$ mathbf{U}^{*^{k+1}}= begin{bmatrix} u^{*^{k+1}}(x_{1},y_{1}) & u^{*^{k+1}}(x_{1},y_{2}) & ldots & u^{*^{k+1}}(x_{1},y_{N_{y}-1})& u^{*^{k+1}}(x_{1},y_{N_{y}}) \ u^{*^{k+1}}(x_{2},y_{1}) & u^{*^{k+1}}(x_{2},y_{2}) & ldots & u^{*^{k+1}}(x_{2},y_{N_{y}-1})& u^{*^{k+1}}(x_{2},y_{N_{y}}) \ vdots & vdots &ddots &vdots & vdots \ u^{*^{k+1}}(x_{N_{x}-1},y_{1}) & u^{*^{k+1}}(x_{N_{x}-1},y_{2}) & ldots & u^{*^{k+1}}(x_{N_{x}-1},y_{N_{y}-1})& u^{*^{k+1}}(x_{N_{x}-1},y_{N_{y}}) \ u^{*^{k+1}}(x_{N_{x}},y_{1}) & u^{*^{k+1}}(x_{N_{x}},y_{2}) & ldots & u^{*^{k+1}}(x_{N_{x}},y_{N_{y}-1})& u^{*^{k+1}}(x_{N_{x}},y_{N_{y}}) end{bmatrix} , $$
The same procedure is applied on Eq. (2.15); as a result, the following system is obtained:
$$ mathbf{A}_{i}Gamma ^{k+1}_{i}= mathbf{B}_{i}, quad i=1, 2, dots , N_{x}, $$
(2.25)
where (Gamma ^{k+1}_{i}= [ gamma _{i1}^{k+1},gamma _{i2}^{k+1},dots , gamma _{iN_{y}}^{k+1} ]^{T}), (mathbf{A}_{i}) and (mathbf{B}_{i}) are (N_{y} times N_{y}) matrix and (N_{y} times 1) vector, respectively, defined as follows:
$$begin{aligned}& mathbf{A}_{i} =mathbf{D} bigl(tilde{mathbf{G}}_{y}^{i} Upsilon _{y}+ tilde{mathbf{G}}^{i}Psi _{y} -mathcal{C}_{Delta t}^{alpha}Phi _{y} bigr)+mathbf{H}_{y}, end{aligned}$$
(2.26)
$$begin{aligned}& mathbf{B}_{i} =mathbf{D} Biggl(mathcal{C}_{Delta t}^{alpha} sum_{p=1}^{k} mathcal{A}_{k}^{alpha} bigl(Phi _{y}overline{Lambda}^{p+1-k}_{i}- Phi _{y}overline{Lambda}^{p-k}_{i} bigr)-mathcal{C}_{Delta t}^{ alpha}Phi _{y} overline{Lambda}^{k}_{i} Biggr). end{aligned}$$
(2.27)
Vector (overline{Lambda}^{m}_{i}) is calculated by the multiplying (Phi ^{-1}) in the mth column of ((mathbf{U}^{*^{k}})^{T}); (Phi _{y}), (Upsilon _{y}), and (Psi _{y}) are (N_{y}times N_{y}) matrices which are constructed same as (Phi _{x}), (Upsilon _{x}), and (Psi _{x}), respectively, with (lbrace y_{1},y_{2},dots ,y_{N_{y}}rbrace ). Moreover, (mathbf{H}_{y}) is an (N_{y} times N_{y}) matrix like (mathbf{H}_{x}) whose first and last rows are equal to the first and the last rows of matrix (Upsilon _{y}); (tilde{mathbf{G}}^{i}) and (tilde{mathbf{G}}_{y}^{i}) have the following structures:
$$begin{aligned}& begin{aligned} tilde{mathbf{G}}^{i}={}& operatorname{diag} bigl(g bigl( biglVert nabla u^{*^{c}}(x_{i},y_{1}) bigrVert bigr),g bigl( biglVert nabla u^{*^{c}}(x_{i},y_{2}) bigrVert bigr), \ & dots ,g bigl( biglVert nabla u^{*^{c}}(x_{i},y_{N_{y}-1}) bigrVert bigr),g bigl( biglVert nabla u^{*^{c}}(x_{i},y_{N_{y}}) bigrVert bigr) bigr), end{aligned} \& begin{aligned} tilde{mathbf{G}}_{y}^{i}={}& operatorname{diag} biggl( frac{partial}{partial y}g bigl( biglVert nabla u^{*^{c}}(x_{i},y_{1}) bigrVert bigr), frac{partial}{partial y}g bigl( biglVert nabla u^{*^{c}}(x_{i},y_{2}) bigrVert bigr), \ & dots ,frac{partial}{partial y}g bigl( biglVert nabla u^{*^{c}}(x_{i},y_{N_{y}-1}) bigrVert bigr)g(x_{i},y_{N_{y}-1}), frac{partial}{partial y}g bigl( biglVert nabla u^{*^{c}}(x_{i},y_{N_{y}}) bigrVert bigr) biggr), end{aligned} end{aligned}$$
By solving (N_{x}) systems (2.25), (hat{u}^{**^{k+1}}) is calculated as
$$ mathbf{U}^{**^{k+1}}=Phi bigl[Gamma _{1}^{k+1}, Gamma _{2}^{k+1}, dots ,Gamma _{N_{x}}^{k+1} bigr], $$
(2.28)
in which
$$ mathbf{U}^{**^{k+1}}= begin{bmatrix} u^{**^{k+1}}(x_{1},y_{1}) & u^{**^{k+1}}(x_{1},y_{2}) & ldots & u^{**^{k+1}}(x_{1},y_{N_{y}-1})& u^{**^{k+1}}(x_{1},y_{N_{y}}) \ u^{**^{k+1}}(x_{2},y_{1}) & u^{**^{k+1}}(x_{2},y_{2}) & ldots & u^{**^{k+1}}(x_{2},y_{N_{y}-1})& u^{**^{k+1}}(x_{2},y_{N_{y}}) \ vdots & vdots &ddots &vdots & vdots \ u^{**^{k+1}}(x_{N_{x}-1},y_{1}) & u^{**^{k+1}}(x_{N_{x}-1},y_{2}) & ldots & u^{**^{k+1}}(x_{N_{x}-1},y_{N_{y}-1})& u^{**^{k+1}}(x_{N_{x}-1},y_{N_{y}}) \ u^{**^{k+1}}(x_{N_{x}},y_{1}) & u^{**^{k+1}}(x_{N_{x}},y_{2}) & ldots & u^{**^{k+1}}(x_{N_{x}},y_{N_{y}-1})& u^{**^{k+1}}(x_{N_{x}},y_{N_{y}}) end{bmatrix} . $$
Finally, (u^{k+1}(x_{i},y_{j})), (i=1,2,dots ,N_{x}), (j=1,2,ldots ,N_{y}) are obtained from ((mathbf{U}^{**^{k+1}} )^{T}).
We summarize the proposed approach in Algorithm 1.
The developed numerical algorithm to solve FPMM
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
Disclaimer:
This article is autogenerated using RSS feeds and has not been created or edited by OA JF.
Click here for Source link (https://www.springeropen.com/)
