Implementation of general linear methods for Volterra integral equations

Abstract

General linear methods (GLMs) are a large family of methods which have been already introduced for the numerical solution of Volterra integral equations of the second kind. In this paper, various issues relevant to the implementation of these methods in a variable stepsize environment are discussed, and the corresponding code is developed.

Introduction

In this paper we describe a code for the numerical solution of Volterra integral equations (VIEs) of the second kind

$$y\left(t\right)=g\left(t\right)+{\int }_{t_0}^{t}K(t,s,y\left(s\right))ds,t\in I≔[t_0,\mathrm{T}],$$

where the forcing function $g:I\to {\mathbb{R}}^m$ and the kernel $K:\Delta \times {\mathbb{R}}^m\to {\mathbb{R}}^m$ with $\Delta ≔\{(t,s):t_0\le s\le t\le \mathrm{T}\}$ are assumed to be sufficiently smooth. Various versions of numerical methods have been investigated to approximate the solution of VIEs (1.1); see, for instance, [1], [2], [3], [4], [5], [6], [7], [8], [9]. Moreover, some efficient codes for VIEs based on various class of numerical methods have been already introduced: variable step size predictor–corrector schemes for nonlinear second kind VIEs [10]; a variable stepsize one-step method of collocation type for general VIEs of the second kind [11]; a variable stepsize algorithm for VIEs of convolution type based on an embedded pair of Runge–Kutta methods of order $p=5$ and $p=4$ [12]; variable stepsize two-step continuous methods for VIEs [13]; a variable stepsize code for VIEs of the second kind [14] based on the natural Volterra Runge–Kutta method of order $p=4$ derived in [15].

In this paper we develop a variable stepsize code based on general linear methods (GLMs) of order four for VIEs (1.1). These methods were introduced by Izzo et al. [9] and investigated more by Abdi et al. in [16] and [17]. For a given uniform grid

$$t_0<t_1<\cdots <t_N,$$

on the interval $I$ with the fixed stepsize $h≔t_{j+1}-t_j$, $j=0,1,\dots ,N-1$, in a GLM with the abscissa vector $c$ and the coefficients matrices

$$A=\left[a_{ij}\right]\in {\mathbb{R}}^{s\times s},U=\left[u_{ij}\right]\in {\mathbb{R}}^{s\times r},$$$$B=\left[b_{ij}\right]\in {\mathbb{R}}^{r\times s},V=\left[v_{ij}\right]\in {\mathbb{R}}^{r\times r},$$

for the numerical solution of VIEs (1.1), the quantities imported into and evaluated in step number $n$ are related by

$$\left\{\begin{array}{c}{Y}_i^{\left[n\right]}=\sum _{j=1}^s{a}_{ij}\left(F_h^{\left[n\right]}\left(t_{n-1,j}\right)+{\Phi }_h^{\left[n\right]}\left(t_{n-1,j}\right)\right)+\sum _{j=1}^r{u}_{ij}{y}_j^{[n-1]},i=1,2,\dots ,s,\hfill \\ y_i^{\left[n\right]}=\sum _{j=1}^s{b}_{ij}\left(F_h^{\left[n\right]}\left(t_{n-1,j}\right)+{\Phi }_h^{\left[n\right]}\left(t_{n-1,j}\right)\right)+\sum _{j=1}^r{v}_{ij}{y}_j^{[n-1]},i=1,2,\dots ,r,\hfill \end{array}\right.$$

for $n=1,2,\dots ,N$, in which $s$ is the number of internal stages, $r$ is the number of external stages, $t_{n-1,j}=t_{n-1}+c_{j}h$, $j=1,2,\dots ,s$, and $F_h^{\left[n\right]}\left(t_{n-1,j}\right)$ and ${\Phi }_h^{\left[n\right]}\left(t_{n-1,j}\right)$ are respectively approximations of sufficiently high order to the lag term $F^{\left[n\right]}\left(t_{n-1,j}\right)$ defined by

$$F^{\left[n\right]}\left(t_{n-1,j}\right)=g\left(t_{n-1,j}\right)+{\int }_{t_0}^{t_{n-1}}K(t_{n-1,j},s,y\left(s\right))ṣ,$$

and the increment term ${\Phi }^{\left[n\right]}\left(t_{n-1,j}\right)$ defined by

$${\Phi }^{\left[n\right]}\left(t_{n-1,j}\right)={\int }_{t_{n-1}}^{t_{n-1,j}}K(t_{n-1,j},s,y\left(s\right))ṣ.$$

Here the stage values $Y_i^{\left[n\right]}$, $i=1,2,\dots ,s$, are approximations of stage order $q$ to the solution $y$ of (1.1) at the stage points $t_{n-1}+c_{i}h$, $i=1,2,\dots ,s$, i.e., 

$$Y_i^{\left[n\right]}=y\left(t_{n-1,i}\right)+\mathcal{O}\left(h^{q+1}\right),$$

Also, the vector of external values $y^{[n-1]}≔{\left[y_i^{[n-1]}\right]}_{i=1}^r$ denotes the vector of approximations imported into step $n$ and $y^{\left[n\right]}≔{\left[y_i^{\left[n\right]}\right]}_{i=1}^r$ is the vector of the quantities computed in this step and exported for use in the subsequent step. Assuming the solution $y$ has derivatives up to order $p$, the input and output vectors $y^{[n-1]}$ and $y^{\left[n\right]}$ are respectively approximations of order $p$ to the linear combinations of the scaled derivatives of the solution $y$ at the points $t_{n-1}$ and $t_n$, i.e., 

$$y^{[n-1]}≔(W\otimes I_m)z(t_{n-1},h)+\mathcal{O}\left(h^{p+1}\right),$$

and

$$y^{\left[n\right]}≔(W\otimes I_m)z(t_n,h)+\mathcal{O}\left(h^{p+1}\right),$$

where $I_m$ is the identity matrix of dimension $m$, the vector $z(t,h)$ is the Nordsieck vector defined by

$$z(t,h)=\left[\begin{array}{c}\hfill y\left(t\right)\hfill \\ \hfill h{y}^{\prime }\left(t\right)\hfill \\ \hfill ⋮\hfill \\ \hfill h^p{y}^{\left(p\right)}\left(t\right)\hfill \end{array}\right],$$

and the matrix $W$ is a $r\times (p+1)$ matrix defined by

$$W=\left[\begin{array}{cccc}\hfill {\alpha }_{10}\hfill & \hfill {\alpha }_{11}\hfill & \hfill \cdots \hfill & \hfill {\alpha }_{1p}\hfill \\ \hfill {\alpha }_{20}\hfill & \hfill {\alpha }_{21}\hfill & \hfill \cdots \hfill & \hfill {\alpha }_{2p}\hfill \\ \hfill ⋮\hfill & \hfill ⋮\hfill & \hfill \ddots \hfill & \hfill ⋮\hfill \\ \hfill {\alpha }_{r0}\hfill & \hfill {\alpha }_{r1}\hfill & \hfill \cdots \hfill & \hfill {\alpha }_{rp}\hfill \end{array}\right],$$

for some real parameters ${\alpha }_{ik}$. The integer $p$ is referred to as the order of the method (1.2). We assume that the approximations $F_h^{\left[n\right]}\left(t_{n-1,j}\right)$ and ${\Phi }_h^{\left[n\right]}\left(t_{n-1,j}\right)$ take the following forms

$$F_h^{\left[n\right]}\left(t_{n-1,j}\right)=g\left(t_{n-1,j}\right)+h\sum _{\nu =1}^{n-1}\sum _{\ell =1}^s{b}_{\ell }K(t_{n-1,j},t_{\nu -1,\ell },Y_{\ell }^{\left[\nu \right]}),$$

and

$${\Phi }_h^{\left[n\right]}\left(t_{n-1,j}\right)=h\sum _{\ell =1}^s{w}_{j,\ell }K(t_{n-1,j},t_{n-1,\ell },Y_{\ell }^{\left[n\right]}),$$

with the quadrature weights $b_{\ell }$ and $w_{j,\ell }$ depending on the abscissae $c_i$, $i=1,2,\dots ,s$, and precomputed in advance.

The GLM (1.2) is usually represented by its coefficients matrices in a partitioned $(s+r)\times (s+r)$ matrix of the form (cf. [9], [16], [17])

together with the quadrature weights appearing in (1.6), (1.7) in the vector $b={\left[b_1{b}_2\cdots b_s\right]}^T\in {\mathbb{R}}^s$ and the matrix

$$\overline{W}=\left[\begin{array}{cccc}\hfill w_{11}\hfill & \hfill w_{12}\hfill & \hfill \cdots \hfill & \hfill w_{1s}\hfill \\ \hfill w_{21}\hfill & \hfill w_{22}\hfill & \hfill \cdots \hfill & \hfill w_{2s}\hfill \\ \hfill ⋮\hfill & \hfill ⋮\hfill & \hfill \ddots \hfill & \hfill ⋮\hfill \\ \hfill w_{s1}\hfill & \hfill w_{s2}\hfill & \hfill \cdots \hfill & \hfill w_{ss}\hfill \end{array}\right].$$

Pre-consistency, consistency, and zero-stability of the GLMs (1.2) have been discussed in [9], [17]. Here, we recall that zero-stability of these methods is governed by the matrix $V$; in fact, the GLM (1.2) is zero-stable if the matrix $V$ is power bounded, i.e., $sup\{\Vert V^n\Vert :n=1,2,\dots \}<\infty$. Now, we state the fundamental theorem on GLMs.

Theorem 1 [9]

A zero-stable and consistent GLM (1.2) together with the quadrature rules (1.6) and (1.7), satisfying

$$V^{\bar{k}}B\mathrm{e}=0,\text{for some}\bar{k}\in \mathbb{N},$$

is convergent.

Introducing the matrices $D$, and $C\in {\mathbb{R}}^{s\times (p+1)}$ as

$$C≔{\left[\mathrm{e}c\cdots c^p\right]}^T,D=diag\{0!,1!,2!,\dots ,p!\},$$

with $\mathrm{e}$ as $s$-dimensional all-ones vector, and $c^j$ as the component-by-component power of the vector $c$, and also the matrix $P=\left[p_{ij}\right]\in {\mathbb{R}}^{(p+1)\times (p+1)}$ as an upper triangular matrix with

$$p_{ij}=\frac{(j-1)!}{(j-i)!},j\ge i,$$

necessary and sufficient conditions for GLMs to be of order $p$ and stage order $q=p$ are given in the following theorem whose proof was presented in [17].

Theorem 2 [17]

Assume that $y^{[n-1]}$ satisfies (1.4). Then, the GLM (1.2) is of order $p$ and stage order $q=p$, i.e. it satisfies (1.3) and (1.5), if and only if

$$AC+UWD=C,$$$$BC+VWD=WP.$$

Examples of GLMs of order $p$ and stage order $q=p$ up to order four with desirable stability properties have been given in [9], [16], [17]. The constructed methods in [9] are in the Nordsieck form, i.e., the matrix $W$ is the identity matrix while in [17], by introducing the order conditions of the methods in general form, rather than the Nordsieck form, methods of order $p$ and stage order $q=p$ up to order four have been constructed. The methods of orders three and four in [17] have been constructed by setting arbitrary values for free parameters of the methods. These methods have bounded stability regions. In [16], the methods of orders three and four have been constructed by minimizing the objective function for the negative area of the intersection of the stability region of the method with the negative half plane. The stability region of the obtained methods in this way are much larger than those of the methods in [17]. In this paper, we show that the proposed methods in [16], [17] are not useable in a variable stepsize environment in the sense of zero-stability property. Then, we construct a Nordsieck GLM of order and stage order four which is zero-stable in a variable stepsize environment for any stepsize pattern and develop a variable stepsize code for VIEs of the second kind based on that.

The rest of the paper is organized along the following lines. Transformation of the GLMs (1.2) corresponding to a nonsingular matrix $T=\left[t_{ij}\right]\in {\mathbb{R}}^{r\times r}$ is discussed in Section 2. In Section 3, implementation of Nordsieck GLMs with $r=p+1$ in which the matrix $W$ is equal to the identity matrix $I_r$, in a variable stepsize environment by using Nordsieck technique is investigated. In Section 4, it is shown that the constructed GLMs in [16], [17] are not useable in a variable stepsize environment in the sense of zero-stability property, and then an example of Nordsieck GLM of order and stage order four is given which is unconditionally zero-stable for any step size pattern. This method is utilized in our developed Matlab code which is described in Section 5. Some numerical experiments confirming the robustness and efficiency of the designed code are presented in Section 6. Finally, some concluding remarks are given in Section 7 where some plans for future work are also presented.