On a multiwavelet spectral element method for integral equation of a generalized Cauchy problem

In this paper we deal with construction and analysis of a multiwavelet spectral element scheme for a generalized Cauchy type problem with Caputo fractional derivative. Numerical schemes for this type of problems, often suffer from the draw-back of spurious oscillations. A common remedy is to render the problem to an equivalent integral equation. For the generalized Cauchy type problem, a corresponding integral equation is of nonlinear Volterra type. In this paper we investigate wellposedness and convergence of a stabilizing multiwavelet scheme for a, one-dimensional case (in [a, b] or [0, 1]), of this problem. Based on multiwavelets, we construct an approximation procedure for the fractional integral operator that yields a linear system of equations with sparse coefficient matrix. In this setting, choosing an appropriate threshold, the number of non-zero coefficients in the system is substantially reduced. A severe obstacle in the convergence analysis is the lack of continuous derivatives in the vicinity of the inflow/ starting boundary point. We overcome this issue through separating a J (mesh)-dependent, small, neighborhood of a (or origin) from the interval, where we only take L2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L_2$$\end{document}-norm. The estimate in this part relies on Chebyshev polynomials, viz. As reported by Richardson( Chebyshev interpolation for functions with endpoint singularities via exponential and double-exponential transforms, Oxford University, UK, 2012) and decreases, almost, exponentially by raising J. At the remaining part of the domain the solution is sufficiently regular to derive the desired optimal error bound. We construct such a modified scheme and analyze its wellposedness, efficiency and accuracy. The robustness of the proposed scheme is confirmed implementing numerical examples.

Our objective is to study the spectral element method for the nonlinear differential equation where c D α a + denotes the Caputo fractional differential operator defined via Riemann-Liouville fractional derivatives and the integral representation: x a u (n) (t)dt (x − t) α−n+1 =: I n−α a + D n (u)(x), D := d dx . (1.2) In this paper we consider the study of the one-dimensional case where, for simplicity, we shall let [a, b] := [0, 1] = Ω, and denote f on the right hand side of (1.1) by For u ∈ AC n (Ω) ( space of absolutely continuous complex-valued functions with continuous derivatives up to order n − 1 in Ω) the Caputo fractional derivative exists almost everywhere in Ω, see, e.g. [22]. Below are some related, previous, studies of the problem (1.1). Kilbas et al. [21] derived existence of unique L 1 solution for Cauchy type problem with the Riemann-Liouville fractional derivative: (1.4) under the assumption that f (x, y) ∈ L 1 [a, b] satisfies a Lipschitz condition in y.
Their investigation is based on rendering (1.4) into the equivalent Volterra integral equation: (1.5) and Banach fixed point theorem. This problem was first considered by Picher and Swell [32], for α ∈ (0, 1) and a bounded, Lipschitz, function f (x, y). System of such problems is studied by Bonilla et al. [6]. Existence of a unique solution for (1.4), with Caputo fractional derivative, is given in [11,13,14]. In [14] a "shadowing-like" approach reduces, both in linear and nonlinear cases, the non-rational fractional orders to rational ones and then follow Grönwall type argument. Kilbas and Marzan [21] studied (1.4), in a similar integral equation as (1.5), with Caputo fractional derivative: (1.6) Some other approaches are Laplace transform [17], Adomian decomposition method [9]. System of Caputo fractional differential equations: investigated in [8] and [14], contain full wellposedness proofs. This approach was generalized to nonlinear Cauchy type problem with the Riemann-Liouville fractional derivative of order α ∈ C ( (α) > 0) in [22]: where 0 < (α 1 ) < · · · < (α σ ) < (α), σ ≥ 2, giving conditions for a unique L 1 (a, b) solution of (1.6) with α ∈ C. An explicit approach, via multivariate Mittag-Leffler functions, for the linear version is given in [26]. In [25] an operational method is introduced, where using B-spline functions the multi-order fractional differential equation: is solved combining the operational matrix of the Caputo fractional derivatives with the collocation method, and without using the equivalent Volterra integral equation form.
Other related studies are considering initial boundary value problems of subdiffusion type and combine the spectral Galerkin in time with finite element [18] or finite difference schemes [42] in space. They construct efficient schemes tackling singularities that substantially raise the convergence rate.
A main strength of wavelets lies on the property of removing/reducing the singularities, see, e.g. [27]. Wavelet bases have been widely used for alternative representations of integral and differential operators [3,28], leading to more thorough study of the underlying equations. As for the numerical studies, wavelet bases have significant advantage, e.g., in the study of adaptivity of discontinuous Galerkin schemes [10,19,[38][39][40][41], collocation method [29], etc. In this setting, Rehman et al. [34], approach (1.4) with the Caputo fractional derivative, using Legendre wavelet method and convert the problem to a system of algebraic equations, hence reducing it to finding unknown coefficients of the system. In an other study, Patera [31], introduced the spectral method into the finite element schemes that uses higher degree piecewise polynomial bases. This procedure may lead to higher order of accuracy, where convergence is achieved either by rasing the spectral order or making mesh refinements (as hp approach). Compared to the standard finite elements, using hp-mesh, faster convergence can be achieved with fewer degrees of freedom. The disadvantage of hp approach is that it is not easily applicable in complex geometries.
Fractional calculus and fractional differential equations have applications in numerous fields of science and engineering, such as material science in mechanical engineering, anomalous diffusion, control and robotics, signal processing and system identifications, friction modeling, wave propagation, turbulence, seepage in fractal media, etc. [4,7,20,23,30,45]. There are several equivalent definitions given for the fractional derivatives, e.g., Grüwald-Letnikov, Riemann-Liouville, and the Caputo fractional derivative [33].
In this paper we consider a multiwavelet approach rather than scalar wavelets. This would allow higher vanishing moments without extending the supports of the involved functions. Thus a smooth function gets negligible projections on most of the bases and, hence, can be locally approximated by lower-order polynomials [1]. The advantageous interpolating property of the scaling functions makes multiwavelet more suitable in solving differential equations: a property that helps to find the coefficients of expansions of a solution. Our objective is to solve the generalized Cauchy type problem (1.1) with Caputo fractional derivative. We represent the fractional integral operator in multiwavelet bases that leads to a sparse representation of the solution operator on a finite interval. Doing so we reduce the problem to the equivalent Volterra integral equation and then apply the multiwavelet spectral element method to reach faster discretization scheme. We investigate existence, uniqueness and derive convergence for the proposed scheme.
An outline of the remaining part of the paper is as follows. In Sect. 2 we prove the existence of a unique solution to our model problem. Section 3 is devoted to the properties of multiwavelets and related projections. In Sect. 4 we introduce multiscale transformation and give a representation of fractional integrals in multiwavelet bases in a sparse matrix form. In Sect. 5 we define the multiwavelets spectral element method and prove error estimates, via a splitting technique, where we employ Chebyshev polynomials, cf Richardson [35], to tackle the lack of continuous derivatives in the vicinity of the origin. Finally, in our concluding Sect. 6, we present numerical results justifying the robustness of our constructed scheme.

Wellposedness
We show the wellposedness for the solution of the differential Eq. (1.1) by verifying the existence of a unique solution for the, equivalent, nonlinear Volterra integral equation: For ordinary differential equations, equivalence of Cauchy type problem, of Riemann-Liouville fractional order, and its corresponding nonlinear Volterra integral equation is proved in Kilbas et al. [22], where also the existence of a unique solution for this problem is proved. In this study, however, we consider the Caputo fractional order rather than the Riemann-Liouville. Since, then verifying the equivalence between (1.1) and (2.1), though similar to the proof of Theorem 3.24 in [22], is much easier. Hence, we establish existence of a unique solution for the Cauchy problem (1.1) in To proceed, for a locally integrable function f ∈ C α n−α (Ω), we define the fractional integrals of order α ∈ R: We shall use the following result by [22]. , where d(Ω) := |Ω| is the size of Ω.
We rewrite the Eq. (2.1), introducing Volterra integral operator K , in the form where To apply Banach fixed point theorem (see Theorem 1.9 [22]), we need to show that for any u 1 , u 2 ∈ C n−1 [a, x 1 ], (satisfying same initial data), Here, the distance in the complex metric space C n−1 (Ω) is given by .
Differentiating (K u)(x), κ times, κ = 1, · · · , n − 1, and using (D κ I α a + )(u)(x) = (I α−κ a + )(u)(x), which holds for α > κ and any sufficiently regular u(x) ∈ C(Ω), we have, using (2.7), The first term on the right hand side: To show continuity for the second term, we use Lemma 2.1 for γ = 0, and a new α, by relabeling: α → α − κ − 1. Then we can prove that the second term is also continuous on [a, x 1 ] and for any κ = 0, 1, · · · , n − 1, and To prove (2.8), by Lipschitz condition (2.5), inequality (2.9), and successive estimates, Hence, the assumption ( By (2.6), 0 < η < 1, hence using Banach fixed point theorem, there exists a unique solutionû ∈ C n−1 [a, Following theorem 1.19 [22], this solution is the limit of a convergence sequence K l (û 0 ) (l → ∞) whereû 0 is any function in C n−1 (Ω 1 ). We can takeû 0 = u 0 , if at least one b κ = 0, put u l = K lû 0 and use (2.7), to get the recursion Note that here the method of successive approximation is used to obtain a unique solution to the nonlinear Volterra integral Eq. (2.1): first for x ∈ (a, x 1 ) and then, in order to establish uniqueness on Ω, we choose at the next step and set Due to the previous step, the first integral at the right hand side of (2.12) is a known function. Hence, by the same argument as above, there exists a unique solutionû( Using (2.11), we get , and the proof is complete.
and B j := {0, 1, · · · , 2 j − 1} for j ∈ Z + ∪ {0}. Further, define the dilation and translation operators: As a property of the Multi-Resolution Analysis (M R A), we can introduce the nested subspaces: where ⊕ denotes orthogonal sum, and a family of bases functions ψ k j,b that generates W r j : The functions ψ k = ψ k 0,0 are called multiwavelets. Alpert's multiwavelets are constructed using dual basis ψk˜j ,b with the bi-orthogonality condition: A corresponding duality for the scaling functions is defined replacing ψ by ϕ in (3.3).
Due to the fact that V r j ⊂ V r j+1 and W r j ⊂ V r j+1 , the vector functions Φ r ,0 matrices with elements obtained from the inner products using the orthonormality property of wavelets and scaling functions. To identify the closed form for multiwavelets, from (3.5), there are 2r 2 unknown coefficients to be found. This is achieved using (i) The orthonormality that yields 2r of the unknown coefficients of (3.5) using: (ii) The number of, N k ψ = k + r − 1, vanishing moments defined by gives the remaining 2r (r − 1) coefficients.

Projection onto V r J
For a fixed integer J ≥ 0, we define an orthonormal projection operator P r J ( f ) :

Multiscale transformation
Using (3.2), we write the multiscale decomposition as V r . Thus, one can approximate any function f ∈ L 2 (Ω) by ISFs of the space V r 0 and multiwavelets of the spaces W r j , j = 0, 1, · · · , J − 1. To proceed, we introduce multiscale operator M r J : i.e., by a multiscale transformation, with .
The single-scale coefficients f k 0,0 , k = 0, 1, · · · , r − 1 are computed using (3.8) and the interpolation property of interpolating scaling functions. But, to evaluate the multiwavelets coefficientsf k j,b , we do not have this property available. Therefore, these integrals are computed numerically, by introducing the N ×N wavelet transform matrix T J , obtained using matrix refinement Eqs. (3.4) and (3.5). Then, the multiwavelets are obtained as the result of T J acting on the multiscaling functions: Below is a brief way to construct T J . In general, for ISFs, the refinement equation between neighboring scales is given by and I 2 j is the identity matrix of order 2 j . Let now the vector function Then, we readily identify the wavelet transform matrix T J as The elements of the matrices G 0 , G 1 , H 0 and H 1 , are, implicitely, determined using (3.4) and (3.5). Further, one can find the reconstruction formula (see also [3]): Relations (3.4), (3.5) and (4.7) yield algorithms for transition between different scales.

Thresholding
By (3.6), each multiwavelet provides vanishing moments of order N k ψ = k + r − 1, k = 0, 1, · · · , r − 1. The Alpert's multiwavelets are uniformly bounded both in L ∞ and L 1 , up to a constant: (4.8) The vanishing moments and normalizations (4.8), imply that the detail coefficients f k J ,b become smaller when the underlying function is locally smooth, and we have, cf [19] , (4.9) where N k ψ is the space of polynomial functions of degree less than N k ψ and W 1,N k ψ (Ω) is the Sobolev space of real valued functions. Thus the detail coefficients decay at the rate 2 −J N k ψ , N k ψ = k + r − 1. Using higher vanishing moments, and increasing the refinement level J , more detail coefficients may be discarded in smooth regions. This yields thresholding with operator T D ε :  The thresholding operator T D ε , acts on the detail coefficients leaving the coarse scale coefficients unaffected. The approximation error due to the thresholding can be estimated similar to that of the classical wavelets, e.g., for the approximation operator

Representation of fractional integral in multiwavelet bases
Recall the fractional integrals of order α ∈ R defined in (2.3): Although the Alpert's multiwavelet bases are discontinuous, by construction, they are locally integrable. Hence, the fractional integral operator I α acting on the vector function Φ r J can be expressed by the projection operator P r J as, To find the entries of the matrix I α φ , we use the following auxiliary result: Let now i = br + k + 1, j = b r + k + 1 for k, k = 0, · · · , r − 1, and b, b ∈ B J , then the coefficients [I α φ ] i, j are given by To compute the integral in (4.16), correlating the integration interval and the support of φ r J ,b , we get the following three cases: Here, by the change of variable x = 2 J t − b and with λ := τ k + b − b, (4.16) can be rewritten as Using definition of the multi-scaling function φ k in (3.1) and Lemma 4.1, we have By a further change of variable: x = λt, we can rewrite (4.19) as where B is the β function and 2 F 1 is the hypergeometric function defined by the power series representation below (see [5]), Here (·) m is the Pochhammer symbol. Since l + 1 − r is a non-positive integer, this series terminates at a finite sum and the function 2 F 1 reduces to the polynomial

and the integral in (4.16) can be rewritten as
(4.22) Using the same argument as in Case (2), we represent (4.22) in form of (4.18). Now, once again, using the definition of scaling function (3.1), Lemma 4.1 and the change of variable x = λy we end up with (4.23) Finally, applying the binomial expansion of (2λy − 1) r −1−l : and further simplification, we get Here the last integral is the incomplete β function which is defined as Hence, letting σ := r − l − m, we can write In this way, the elements of the block upper triangular matrix I α φ are specified in all possible cases.
where D ψ := T J D φ T T J (herein, D φ is the matrix representation of differential operator for scaling functions Φ r J ) is the matrix representing the differential operator D for Alpert's multiwavelets and I α ψ := T J I α φ T T J (see [3]). Approximating (2.1), using (5.2)-(5.3), (4.14) and the wavelet transform matrix T J , the residual in approximating (2.1) is written as Minimizing (5.4) yields the unknown coefficients. The spectral method is implemented either by the Galerkin-or the collocation-method. In the wavelet Galerkin method we use the orthogonality relation r J , Ψ r J = 0 to find the Fourier coefficients for r J associated with Ψ r J . This yields the nonlinear system where the orthonormality property of the Alpert's multiwavelets has been used. Then, the Newton's method is used to solve this nonlinear system.

Error estimate
There are some concerns about the approximation Lemma 3.1. The common wisdom knows that one cannot expect the exact solution to have continuous derivatives in the vicinity of the origin. A remedy cf [35] reads: let f (x) satisfy the condition and set L = ln 2 J . Then one can show that there exists A 1 > 0, such that for a sufficiently small x L ∈ (0, 1), Suppose that the Chebyshev interpolant C n (y) provides an approximation for the function 1], and x L depends on and decreases, exponentially, to zero, as L → ∞. Now we consider C n L , interpolating f (x) on [0, 1], as Under these assumptions, the approximation error is given by Then, similar to the Theorem 3.2 in [35], we can derive the following estimates: and Proof Using (5.9), (5.7), and proof of lemma 3.1 (see [1]), we have Then (5.10) follows with the constant C = max 2M r To estimate P r J f − f , we use the error estimate for the polynomial based on Chebyshev nodes for interpolation, and (5.10), to get (5.11), viz.

Theorem 5.2 Under the assumptions of Theorem 2.1, and for u(x) and u J (x) the exact and multiwavelets Galerkin solutions of Eq. (2.1), respectively, we have that: If
then, for any u(x) ∈ C(Ω), the following estimate holds true , , for someC > 0, and There are two possible cases here: 1. If r < n, then we have and it follows from (2.4) of Lemma 2.1, that 2. If r ≥ n. Combining Lemma 2.21 of [22] and Theorem 3.14 in [15], it is easy to write (5.14) Thus using Lemma 5.1, we have where m = [r − α] + 1.
In general, these two cases can be summarized in one, as where Using (1.3), the hypothesis and proof of the Theorem 2.1 (Lipschitz condition 2.5), and (2.9), we also have the following estimate, Now, we can write Taking C(Ω)-norm of (5.17) and using (5.13)-(5.16) and triangle inequality, we get and hence, sinceη = Λ σ σ j=0 where .

Numerical examples
Here we present some canonical examples to illustrate the validity of the proposed scheme through numerical implementations for the following two configurations: Three linear Cauchy type initial value problems Three nonlinear Cauchy type initial value problems Here if α < 1, then we suppress the data u (0) = 0, and recall the exact solution given in [25]: .
is the Bagley-Torvik equation [13] with exact solution x 7/3 and The equation E 1 (c) demonstrate ability of the proposed method in solving problems whose solutions do not have continuous derivatives near the origin. The exact solution of this problem is given by u(x) = 1 − e π x erfc( √ π x).
The exact solutions for the configuration E 2 (E 2 (a), E 2 (b), and E 2 (c)) all are of the form u(x) = x β , and the source functions f 1 (x) and f 2 (x) are obtained via The coefficient a 1 (x) is a given smooth function.
As for the equation E 2 (c), the exact solution is u(x) = x 3 , for which one can observe that the corresponding function F α 1 ,··· ,α σ [x, u] in (1.3) has not a continuous derivative at zero. In this problem a 2 (x) := x 2 and Using the multiwavelet spectral element method of Sect. 5, the equivalent Volterra integral equation for the configuration E 1 is reduced to a system of linear algebraic equations in matrix form as AU = U 0 . According to (4.9), the entries of the coefficient matrix A decay at the rate of 2 −J N ψ . Thus we expect that by increasing the refinement level J and the vanishing moment of multiwavelets (equivalently r ), the coefficients tend to zero for the sufficiently smooth underlying functions. Using a proper threshold, one may set most coefficients to zero and thus obtain a sparse coefficient matrix A.
Obviously this process increases the computational speed but it may also increase the error when choosing a larger threshold. The error is controlled by a given tolerance ξ : To preserve the desired accuracy, the optimal threshold value (6.1) must be used. Let N ε be the number of nonzero elements of the coefficient matrix after thresholding. Then, the compression percentage (percentage of matrix sparsity) P ε is defined by To check efficiency and accuracy of the proposed method for the configuration E 1 , we illustrate the effects of the refinement level J and the multiplicity parameter r in Figs. 1, 5 and 6. We observe that by increasing these parameters the L 2 errors decrease and thus resulting in higher CPU times. Since our goal is to decrease the computational cost for linear Cauchy type problem by thresholding, we plot the effect of varying threshold values on the coefficient matrix A and the CPU time (see Figs. 2 and 3). In Fig. 3, the percentage of compression is also presented. One can see that with increasing threshold value, the number of elements of the coefficients matrix decreases which results in lower CPU time and a higher percentage of compression. Figure 4 shows the effect of the thresholding on L 2 error with different values of r , J and ε. Note that with increasing threshold value the error increases. Table 1 shows the effect of thresholding on L 2 -error , for the configuration E 1 (a) and with α = 0.85, r = 5 and J = 4. Table 2 shows the effect of thresholding, for the refinement parameter J (= 2, 3) and the multiplicity r (= 3, 2), on L 2 -error for the configuration E 1 (b) and with different ε-values. Table 5 gives the L 2 error for the configuration E 2 (a), for different α, β, and for J = 2, 3, and r = 2, 3, 4, 5, 6. We observe that, by increasing the refinement level J and the multiplicity parameter r , the L 2 error decreases. As for the configuration E 2 (b): choosing α = 2.5, α 1 = 1.5 α 2 = 0.5, different values for β, and various functions for a 1 (x), we can see the effect of parameters J and r in Fig. 7.
We return to the critical case of non-differentiability in the vicinity of the origin, i.e. the problems whose corresponding "F"-functions do not have a continuous derivative near the origin.
In the linear case, Tables 3-4 demonstrate the effect of thresholding on the L 2 error and the percentage of matrix sparsity for E 1 (c). As for the nonlinear Example E 2 (c): Fig. 8 is a further justification for robustness of the error analysis, where it can be seen that, by increasing J , the error decreases exponentially, and the harm near x = 0 is indeed removable. To verify the accuracy of our approach we compare Implementing E 1 (a), with the existing results in the literature, e.g. in [12,24,36]: In Table 6, the absolute error for the proposed method at x = 1 is compared with those of the methods presented in [12,24]. As seen, we get a substantially better convergence with, relatively, larger mesh size h (lines 5 and 6 of the table compared with the lines 1 and 2 in [12], and lines 3 and 4 in [24]). In Table 7 a comparison is made with the results of the method in [36], for different values of both α and x. Once again, the results confirm that our method with lower spectral order (degree of approximating polynomilal) gives better accuracy than in [36]. For α = 1 and α = 2, the exact solutions for E 1 (a) is reported in [36] as u(x) = e −x and u(x) = cos(x), respectively. Our approximate solutions for different values of α, and with r = 3 and J = 3, plotted in Fig. 9, shows that the numerical solutions converge to the analytical ones as α approaches an integer order. This comfirms that the solution of the fractional differential equation approaches to that of the integer-order differential equation. This is yet a further test to justisy the merit of our approach.
All examples are carried out with the combined use of Maple and Matlab softwares. Altogether confirms the advantageous effects of thresholding discussed in the paper.     . 8 Effects of the refinement level J (left) and the multiplicity parameter r (right) on L 2 error for configuration E 2 (c) Table 6 Comparison of absolute value error reported in [12,24] and proposed method at x = 1 for E 1 (a) Method h α = 0.5 α = 1.5 Predictor-Corrector approach [12]

Conclusion
We propose a reliable and efficient scheme based on multiwavelet spectral element method for numerical solution of a generalized Cauchy type problem with Caputo fractional derivative. Under certain assumptions, including Lipschitz continuity of the right hand side, we prove existence of a unique solution for the problem. To this approach, we start transforming the equation into a Volterra integral equation and then reduce it to an algebraic system. Our proposed scheme is based on representing the fractional integral operator in the multiwavelet bases as a sparse coefficient matrix. For the resulting equation, in nonlinear form, the numerical computations are more challenging. We propose an adequate numerical scheme to deal with this problem. In the linear case, selecting an appropriate threshold, the number of non-zero coefficients in the system is substantially reduced, thus obtaining a fast convergence with lower cost. More specifically, selecting the high values of the refinement level J and multiplicity parameter r , the error will decrease correspondingly. Furthermore, due to this observation, it is obvious that, by increasing the threshold the CPU time decreases and the error increases. A main obstacle in error estimates is the lack of continuous derivatives of the solution near the origin. This, however, is tackled by a splitting technique, separating a local neighborhood of the origin, where L 2 -estimate based on Chebyshev polynomials is employed. This approach can be extended to two-dimensional domains, with no reentrants, using 2D integral equation representation, and is the subject of a forthcoming study (see also [10]). This, however, requires excessive CPU time to implement. Implementing, goal oriented, numerical examples show the robustness of this approach. The scheme yields a desired converges for an appropriate choice of threshold and is cost-effective where we may, substantially, reduce the number of non-zero coefficients of the system. Here we mention that the constructed wavelet transform matrix T j has an epical role in decreasing the number of nonzero coefficients, hence, reducing the computational load. More specifically, T j converts the coefficients in the expasions based on scaling functions to those in multiwavelets. Here, the coefficients of the expansion with respcet to the scaling functions are detemined using interpolation, and no integration is performed. The matrix I α φ , which defines the fractional integral operator, and as a spars matrix, reduces the computational load, is not considered elsewhere.