A Collocation Method in Spline Spaces for the Solution of Linear Fractional Dynamical Systems

TL;DR

This paper introduces a collocation method using spline spaces to efficiently solve linear fractional dynamical systems, with proven convergence and demonstrated numerical results.

Contribution

It presents a novel collocation approach in spline spaces specifically designed for fractional dynamical systems, including convergence proof and practical numerical validation.

Findings

Efficient evaluation of collocation matrices using B-spline derivatives

Convergence of the proposed collocation method established

Numerical results demonstrate method effectiveness

Abstract

We used a collocation method in refinable spline space to solve a linear dynamical system having fractional derivative in time. The method takes advantage of an explicit derivation rule for the B-spline basis that allows us to efficiently evaluate the collocation matrices appearing in the method. We proof the convergence of the method. Some numerical results are shown.

Equations

$$\left\{\begin{array}[]{ll}D^{\gamma}_{t}X(t)=A\,X(t)\,,&t>0\,,\quad 0<\gamma<1\,,\\ \\ X(0)=X_{0}\,,\end{array}\right.$$

$$D_{t}^{\gamma}X(t):=\bigl{[}D_{t}^{\gamma}x_{1}(t),D_{t}^{\gamma}x_{2}(t),\ldots,D_{t}^{\gamma}x_{m}(t)\bigr{]}^{T}\,,$$

$$D_{t}^{\gamma}\,x(t):=\bigl{(}{\cal J}^{(k-\gamma)}x^{(k)}\bigr{)}(t)\,,\quad k-1<\gamma<k\,,\quad k\in\mathbb{N}\,,\quad t>0\,,$$

$${\mathcal{J}}^{(\gamma )}x(t):=\frac{1}{\Gamma (\gamma )}{\int }_0^{t}x(\tau )(t-\tau {)}^{\gamma -1}d\tau ,$$

$$X(t)=E_{\gamma ,1}(t^{\gamma },A)X_0,$$

$$E_{\gamma ,\beta }(z,A)=\sum _{k\ge 0}\frac{(zA{)}^k}{\Gamma (\gamma k+\beta )},z\in \mathbb{C},A\in {\mathbb{R}}^{m\times m},$$

$$T_{n}(t):=\bigl{(}\max(0,t)\bigr{)}^{n},\quad n\in\mathbb{N}\cup{0}\,,$$

$${\Delta }^{n}f(t):=\sum _{\ell =0}^n\left(\genfrac{}{}{0ex}{}{n}{\ell }\right)(-1{)}^{l}f(t-\ell ),n\in \mathbb{N}.$$

$$B_n(t):=\frac{1}{n!}{\Delta }^{n+1}{T}_n(t),n\in \mathbb{N}\cup 0.$$

$${\mathcal{B}}_n=\{B_n(t-\ell ),-n\le \ell \},t\in [0,\infty ),$$

$$\text{mbox}\text{calligra}s(t)=\sum _{\ell \ge -n}{c}_k{B}_n(t-\ell ).$$

$$D_t^{\gamma }\text{mbox}\text{calligra}s(t)=\sum _{\ell \ge -n}{c}_k{D}^{\gamma }{B}_n(t-\ell ),k-1<\gamma <k,1\le k\le n-1.$$

$$B_{n,\ell }(t)=B_n(t-\ell ),\ell \ge -n.$$

$$D_t^{\gamma }{B}_n(t)=\frac{1}{\Gamma (n-\gamma +1)}{\Delta }^{n+1}{T}_{n-\gamma }(t),$$

$$T_{\gamma}(t):=\bigl{(}\max(0,t)\bigr{)}^{\gamma},\quad\gamma\geq 0,$$

$$D_t^{\gamma }{B}_{n,\ell }(t)=\frac{1}{\Gamma (n-\gamma +1)}{\Delta }^{n+1}{T}_{n-\gamma }(t-\ell ),\ell \ge 0,$$

$$D_t^{\gamma }{B}_{n,\ell }(t)=\frac{{\Delta }^{n+1}{T}_{n-\gamma }(t-\ell )}{\Gamma (n+1-\gamma )}-\frac{1}{\Gamma (1-\gamma )}{\int }_0^{-\ell }\frac{B_n^{\prime }(\tau )}{(t-\ell -\tau {)}^{\gamma }}d\tau ,-n\le \ell \le -1.$$

$$D_t^{\gamma }{B}_{n,\ell }(t)$$

$$\frac{1}{\Gamma (1-\gamma )}{\int }_{\ell }^0\frac{B_{n,\ell }^{\prime }(\tau )}{(t-\tau {)}^{\gamma }}d\tau =\frac{1}{\Gamma (1-\gamma )}{\int }_{\ell }^0\frac{B_n^{\prime }(\tau -\ell )}{(t-\tau {)}^{\gamma }}d\tau =$$

$$\frac{1}{\Gamma (1-\gamma )}{\int }_0^{-\ell }\frac{B_n^{\prime }(\tau )}{(t-\ell -\tau {)}^{\gamma }}d\tau$$

$$\begin{array}[]{ll}\displaystyle\frac{1}{\Gamma(1-\gamma)}\displaystyle\int_{0}^{-\ell}\frac{B^{\prime}_{n}(\tau)}{(t-\ell-\tau)^{\gamma}}\,d\tau=\\ \\ \displaystyle\rule{28.45274pt}{0.0pt}\frac{1}{\Gamma(n+1-\gamma)}\sum_{r=0}^{-\ell-1}(-1)^{r}\binom{n+1}{r}\bigg{(}(t-\ell-r)^{(n-\gamma)}+\\ \\ \displaystyle\rule{28.45274pt}{0.0pt}t^{1-\gamma}\,\sum_{p=0}^{n-1}\frac{(-1)^{n-p}(-\ell-r)^{n-1-p}(t-\ell-r)^{p}}{(n-1-p)!}\prod_{s=1}^{n-1-p}(\gamma-s)\bigg{)}.\end{array}$$

$$B_n^{\prime }(t)=\frac{1}{n!}{\Delta }^{n+1}{T}_n^{\prime }(t)=\frac{1}{(n-1)!}\sum _{r=0}^{n+1}(-1{)}^r(\genfrac{}{}{0ex}{}{n+1}{r})T_{n-1}(t-r).$$

$$\begin{array}[]{lcl}\displaystyle\frac{1}{\Gamma(1-\gamma)(n-1)!}\sum_{r=0}^{n+1}(-1)^{r}\binom{n+1}{r}\int_{0}^{-\ell}\frac{T_{n-1}(\tau-r)}{(t-\ell-\tau)^{\gamma}}\,d\tau=\\ \\ \displaystyle\frac{1}{\Gamma(1-\gamma)(n-1)!}\sum_{r=0}^{-\ell-1}(-1)^{r}\binom{n+1}{r}\int_{0}^{-\ell-r}\frac{\tau^{n-1}}{(t-\ell-r-\tau)^{\gamma}}d\tau\,.\end{array}$$

$$\begin{array}[]{l}\displaystyle\int_{0}^{-\ell}\frac{\tau^{n-1}}{(t-\ell-\tau)^{\gamma}}\,d\tau=\\ \\ \displaystyle\frac{(n-1)!}{\prod_{s=1}^{n}(\gamma-s)}(t-\ell-\tau)^{(1-\gamma)}\sum_{p=0}^{n-1}\frac{(\ell-t)^{p}\tau^{n-1-p}}{(n-1-p)!}\prod_{s=1}^{n-1-p}(\gamma-s)\bigg{|}_{\tau=0}^{-\ell}=\\ \\ \displaystyle\frac{(n-1)!}{\prod_{s=1}^{n}(\gamma-s)}\left[(t-\ell)^{(n-\gamma)}+t^{1-\gamma}\,\sum_{p=0}^{n-1}\frac{(-1)^{n-p}(-\ell)^{n-1-p}(t-\ell)^{p}}{(n-1-p)!}\prod_{s=1}^{n-1-p}(\gamma-s)\right]\end{array}$$

$$V_j={\mathrm{clos}}_{L_2[0,\infty )}\{{\phi }_{j\ell }(t),\ell \ge -n\},j\in \mathbb{Z},t\ge 0,$$

$${\phi }_{j\ell }(t):=2^{j/2}{B}_n(2^{j}t-\ell ),$$

$$f(t)=\sum _{\ell \ge -n}{c}_{j\ell }{\phi }_{j\ell }(t),$$

$$D_t^{\gamma }f(2^{j}t)=2^{\gamma j}{D}_{2^{j}t}^{\gamma }f(2^{j}t),k-1<\gamma <k,k\in \mathbb{N}.$$

$$\begin{array}[]{lcl}D^{\gamma}_{t}F(t)&=&\displaystyle\frac{1}{\Gamma(k-\gamma)}\int_{0}^{t}\frac{F^{(k)}(\tau)}{(t-\tau)^{\gamma-k+1}}d\tau\\ \\ &=&\displaystyle\frac{1}{\Gamma(k-\gamma)}\,2^{jk}\int_{0}^{t}\frac{f^{(k)}(2^{j}\tau)}{(t-\tau)^{\gamma-k+1}}d\tau\end{array}$$

$$\begin{array}[]{lcl}D^{\gamma}_{t}F(t)&=&\displaystyle\frac{1}{\Gamma(k-\gamma)}2^{j(k-1)}\int_{0}^{2^{j}t}\frac{f^{(k)}(\tau)}{(t-2^{-j}\tau)^{(\gamma-k+1)}}\,d\tau\\ \\ &=&\displaystyle\frac{1}{\Gamma(k-\gamma)}\frac{2^{j(k-1)}}{2^{-j(\gamma-k+1)}}\int_{0}^{2^{j}t}\frac{f^{(k)}(\tau)}{(2^{j}t-\tau)^{(\gamma-k+1)}}d\tau\end{array}$$

$$X_j(t)=\sum _{\ell \ge -n}{C}_{j\ell }{\phi }_{j\ell }(t)\in V_j,C_{j\ell }\in {\mathbb{R}}^m,$$

Full text

Keywords: Fractional differential problem, Projection method, Collocation method, B-Spline

1 Introduction

In recent years fractional differential models were used to describe a great variety of physical phenomena, such as the anomalous diffusion in biological tissues, the viscoelastic properties of smart materials, the growth of population in dynamical systems (see, [8, 20, 22] and references therein). Even if there is a great effort in developing the theory of fractional calculus [5, 13, 19, 20], the analytical solution of fractional differential problems can be obtained in a very few cases. This is why the literature on numerical methods to solve this kind of problems is growing rapidly (see, [2, 11] and the detailed bibliography in the more recent papers [10, 17]).

In [15] two of the authors introduced an efficient collocation method to numerically solve fractional differential problems. In the method the approximating function is assumed to belong to a refinable space and its expression is determined by solving the differential problem in a set of collocation points. Thus, the method is both a projection method and a collocation method and the nonlocal behavior of the fractional derivative can easily be taken into account.

In the present paper we used this method to solve a linear dynamical system having fractional derivative in time. We assume the approximating function belongs to refinable spaces generated by the polynomial splines and we take advantage of an explicit derivation rule for the B-spline basis that allows us to efficiently evaluate the collocation matrix.

The paper is organized as follows. In Section 2 we recall the definition of the Caputo fractional derivative and describe the fractional dynamical system we are interested in. We give also its analytical solution in terms of the matrix Mittag-Leffler function. In Section 3 we describe the B-spline basis we used to construct the approximate solution to the differential system and give the explicit expression of the fractional derivatives of the basis functions. In Section 4 we analyze the refinability properties of the B-spline basis. The collocation method is described in Section 5 where its convergence is also proved. Finally, some numerical tests are provided in Section 6 while some conclusions are drawn in Section 7.

2 Fractional dynamical systems

Let $X(t):\mathbb{R}\to\mathbb{R}^{m}$ be a real-valued vector function, $X_{0}\in\mathbb{R}^{m}$ be a real vector and $A\in\mathbb{R}^{m\times m}$ be a real matrix. We consider the following linear dynamical system:

[TABLE]

having time derivative of fractional, i.e. noninteger, order $\gamma$. In this context, the operator $D_{t}^{\gamma}$ denotes the Caputo fractional derivative with respect to the time $t$. For a sufficiently smooth vector function $X(t)=[x_{1}(t),x_{2}(t),\ldots,x_{m}(t)]^{T}$, the Caputo derivative is defined as

[TABLE]

where

[TABLE]

and

[TABLE]

is the Riemann-Liouville integral operator. Here, $\Gamma(\gamma)$ denotes the Euler’s gamma function. For details on fractional calculus see, for instance, [5, 8, 13, 20].

The existence of a unique solution to (2.1) was proved, for instance, in [5, §7.1]. A detailed analysis of positive linear systems of type (2.1) and of their properties can be found in [7] where the analytical solution in terms of the Mittag-Leffler function is obtained by the Laplace transform. Its explicit expression is

[TABLE]

where

[TABLE]

is the matrix Mittag-Leffler function. We observe that the evaluation of $E_{\gamma,\beta}(z,A)$ is rather cumbersome (cf. [6]). An alternative expression of the analytical solution not involving the matrix Mittag-Leffler function can be found in [5, §7.1].

3 The B-spline basis and its fractional derivatives

In this section we describe the polynomial B-spline basis we will use to approximate the solution to Equation (2.1) and give the analytical expression of its fractional derivative.

The classical cardinal B-splines are piecewise polynomials of integer degree having breakpoints on integer knots (see [3, 21] for details). For our porposes, we define the cardinal B-spline through the truncated power function

[TABLE]

and the backward finite difference operator

[TABLE]

Then, the cardinal B-spline of integer degree $n\geq 0$ is defined as

[TABLE]

The cardinal B-spline $B_{n}$ is a piecewise polynomial of degree $n$ with breakpoints on the integers, compactly supported on $[0,n+1]$ and belonging to $C^{n-1}$.

On the semi-finite interval $[0,\infty)$ the integer translates

[TABLE]

form a function basis for the spline space so that any spline function \calligras   can be represented as

[TABLE]

As a consequence, the fractional derivative of \calligras   can be evaluated as

[TABLE]

Thus, to compute the fractional derivatives of \calligras   we need the fractional derivatives of the functions belonging to the B-spline basis ${\cal B}_{n}$.

Let us denote by $B_{n,\ell}(t)$ the $\ell$-translate of $B_{n}$, i.e.

[TABLE]

First of all, we notice that when $\ell\geq 0$ the functions $B_{n,\ell}$ are interior functions having support $[\ell,\ell+n+1]$ all contained in $[0,\infty)$. Their fractional derivative can be evaluated by the differentiation rule

[TABLE]

where

[TABLE]

is the fractional truncated power function [23]. From (3.7) it follows that the fractional derivative of a polynomial B-spline is a fractional spline, i.e. a piecewice polynomial of noninteger degree. Details on fractional splines can be found in [23].

For $-n\leq\ell\leq-1$, the functions $B_{n,\ell}$ are left edge functions having support $[0,n+\ell+1]$. Their fractional derivative can be explicitly evaluated using definition (2.3) and the differentiation rule (3.7) as the following theorem shows.

Theorem 3.1**.**

For $0<\gamma<1$, the fractional derivative of the B-spline basis functions $B_{n,\ell}$ is given by

[TABLE]

and

[TABLE]

Proof.

The derivation rule (3.9) immediately follows from (3.7). Now, consider the case $-n\leq\ell\leq-1$. From definition (2.3) one has

[TABLE]

The first integral is the fractional derivative of the $\ell$-translate of $B_{n}$ and can be evaluated by the differentiation rule (3.9).

As for the second integral, we get

[TABLE]

so concluding the proof. ∎

In the following theorem we explicitly evaluate the integral appearing in the left hand side of (3.10).

Theorem 3.2**.**

For $-n\leq\ell\leq-1$ the explicitly expression of the integral in (3.10) is

[TABLE]

Proof.

We recall that $B_{n}^{\prime}(t)$ writes:

[TABLE]

Substituting the expression of $B^{\prime}_{n}$ in the integral in the left hand side of (3.10), we obtain

[TABLE]

By a direct computation we get

[TABLE]

and the claim follows. ∎

4 Multiresolution Analysis on $L_{2}[0,\infty)$

The B-spline basis ${\cal B}_{n}$ generates a multiresolution analysis on the semi-infinite interval $[0,+\infty)$ [4]. This means that the sequence of subspaces $\{V_{j}\subset L_{2}[0,+\infty)\}$ defined as

[TABLE]

where

[TABLE]

fulfills the following properties:

[TABLE]

Thus, any function $f\in V_{j}$ can be represented as

[TABLE]

where $\{c_{j\ell}\}\in\ell_{2}(\mathbb{Z})$. Once again, the basis functions $\varphi_{j\ell}$ with $-n\leq\ell\leq-1$ are the left egde functions, while for $\ell\geq 0$ $\varphi_{j\ell}$ are integer translates of $B_{n}(2^{j}\cdot)$. The computation of the fractional derivatives of $f$ requires the evaluation of the fractional derivatives of $\varphi_{j\ell}$. This can be done using Theorem 3.1 and the following lemma.

Lemma 4.1**.**

The Caputo derivative of order $\gamma$ of the $2^{j}$-dilate of a function $f(t)$ is given by:

[TABLE]

Proof.

Let $F(t)=f(2^{j}t)$, then $F^{(m)}(t)=2^{jm}f^{(m)}(2^{j}t)$, $m\in\mathbb{N}$. By definition

[TABLE]

By the change of variables $2^{j}\tau\rightarrow\tau$, we get

[TABLE]

and the claim follows. ∎

In the next section we will describe how to apply the collocation method introduced in [15] to numerically solve the differential problem (2.1) in the refinable spline spaces.

5 The fractional collocation method

We look for an approximating vector function

[TABLE]

that solves the differential problem (2.1) on a set of collocation points. We choose as collocation points the dyadic nodes in which $\varphi_{jk}$ can be efficiently evaluated through well-known recursive algorithms [12].

Let ${\cal I}=[0,T]$ be a finite interval. Without loss of generality we assume $T\in\mathbb{N}$. Since $\varphi_{jk}$ has compact support, for $t\in{\cal I}$ the sum in (5.1) reduces to a finite sum:

[TABLE]

Let us denote by $\{t_{p}=p/2^{s},0\leq p\leq 2^{s}T\}$ the dyadic collocation points in the interval $\cal I$. Substituting (5.2) in (2.1) evaluated on the collocation points gives

[TABLE]

This is a linear algebraic system that can be written in matrix form as follows

[TABLE]

where $I_{m}$ is the identity matrix of order $m$,

[TABLE]

and

[TABLE]

are the collocation matrices of the refinable basis,

[TABLE]

and

[TABLE]

is the unknown vector. We notice that the linear system (5.4) has $m(2^{s}T+1)$ equations and $m(2^{j}T+n)$ unknowns. To guarantee the existence of a unique solution we set $2^{s}T+1\geq 2^{j}T+n$.

Theorem 5.1**.**

For $2^{s}T+1\geq 2^{j}T+n$ the linear system (5.4) has a unique solution.

Proof.

Using definitions (2.2)-(2.4) the differential problem (2.1) can be written as a system of integral equations:

[TABLE]

where $Z(t)=[z_{1}(t)=D^{\gamma}_{t}x_{1}(t),z_{2}(t)=D^{\gamma}_{t}x_{2}(t),\ldots,z_{m}(t)=D^{\gamma}_{t}x_{m}(t)]^{T}$ and $J^{(\gamma)}\,Z=[J^{(\gamma)}\,z_{1}(t),\ldots,J^{(\gamma)}z_{m}(t)]^{T}$. The system above is equivalent to the differential problem (2.1) (cf. [9]) and has a unique solution [24] so that the associated integral operator is invertible. Thus, the linear system (5.4) has a unique solution, too (cf. [1]). ∎

Finally, we proof the convergence of the collocation method (5.3).

Theorem 5.2**.**

The collocation method is convergent, i.e.

[TABLE]

where $\|X(t)\|_{\infty}=\max_{1\leq i\leq m}\left(\max_{t\in[0,T]}|x_{i}(t)|\right)$. Moreover, the approximation order is $\gamma$, i.e.

[TABLE]

where $\kappa$ is a constant independent from $j$.

Proof.

Since the collocation method can be used also to approximate the solution to the system (5.5), the equivalence implies that the approximation error $\|X-X_{j}\|_{\infty}$ is the same as the approximation error $\|Z-Z_{j}\|_{\infty}$ [1, 9]. Now, $Z_{j}$ is a projection operator in the spline space so that the convergence is guarantee with approximation order at least $\gamma$ in the case when $Z$ is sufficiently smooth [3]. ∎

We notice that since the equality $2^{s}T+1=2^{j}T+n$ can be satisfied just in a few special cases, in practice we choose $s$ and $j$ so that $2^{s}T+1>2^{j}T+n$. Thus, the system (5.4) results is an overdetermined linear system that can be solved in the least squares sense.

6 Numerical results

In this section we use the proposed method to solve some test problems. In the tests we used the splines of degree $n=3$ and $n=4$ as approximating functions and set $s=j+1$. The functions of the cubic B-spline basis and their fractional derivatives are shown in Figures 1-2. The ordinary first derivative is also displayed.

To check the accuracy of the approximations obtained by the proposed method, we evaluated the componentwise $L_{\infty}$-norm of the error ${\cal E}_{j}(t)=X(t)-X_{j}(t)$, i.e.

[TABLE]

Moreover, we evaluated the numerical approximation order $\rho_{\gamma,n}(j)$ defined as

[TABLE]

We notice that the matrix Mittag-Leffler function appearing in the analytical solution was evaluated using the procedure proposed in [6].

6.1 Example 1

First of all, we tested the accuracy of the collocation method by solving the following simple fractional differential equation (cf. [5, pg. 137]):

[TABLE]

whose exact solution is $x(t)=t^{2}$. In this case the cubic spline approximation is exact. We numerically solve Equation (6.1) in the interval ${\cal I}=[0,1]$ for $\gamma=0.10$, 0.25, 0.50, 0.75. The table below lists the $L_{\infty}$-norm of the error ${\cal E}_{j}(t)=x(t)-x_{j}(t)$ obtained by the collocation method when $j=7$:

[TABLE]

As expected, the error is in the order of the machine precision.

6.2 Example 2

In the second test we solved the fractional dynamical system

[TABLE]

The exact solution is [5, §7.1]

[TABLE]

where

[TABLE]

is the one-parameter Mittag-Leffler function. We notice that the matrix

[TABLE]

associated with the dynamical system (6.2) has negative eigenvalues, so that the stability of the dynamical system is guaranteed [7].

We solved the differential problem (6.2) by the collocation method described in Section 5 for $\gamma=0.10$, 0.25, 0.50, 0.75, and for different values of $j$. In Figures 3-6 the numerical solution and the approximation error are displayed in the case of the cubic spline approximation and for $j=8$. The numerical solution and the error obtained when solving the classical problem with integer first derivative are displayed in Figure 7. The plots show that the proposed method gives a good accuracy that increases as $\gamma$ increases, i.e. as the smoothness of the analytical solution increases. In Figure 8 the numerical approximation order $\rho_{\gamma,n}(j)$ is displayed as a function of $j$ for different values of $\gamma$ and $n=3$ and $n=4$. The plots show that the numerical approximation order is in accordance with the theoretical one. Moreover, the error is lower for the spline of degree 4. Finally, in Figure 9 the numerical convergence order $\rho_{\gamma,n}(j)$ for $n=3$ and $n=4$ is displayed in the case of ordinary first derivative. We observe that in this case the theoretical convergence order is $n+1$ (cf. [3]).

7 Conclusions

We used a collocation method based on an interpolating projection operator on refinable polynomial spline spaces to approximate the solution of a linear fractional dynamical system. We provide an explicit formula that allows us to evaluate the fractional derivatives of the approximating function in an accurate and easy way. The method can be used to solve several differential problems of fractional order and, in particular, nonlinear problems [17, 18] or boundary value problems [16]. We notice that higher approximation order methods can be obtained by using different types of collocation points as in [1, 9, 14]. This will be the subject of a forthcoming paper.

Acknowledgment

This work was partially supported by: grant of University of Roma “La Sapienza”, Ricerca Scientifica 2017; INdAM-GNCS 2018 Project “Sviluppo di modelli e metodi computazionali per l’elaborazione di segnali e immagini”.