Variable stepsize SDIMSIMs for ordinary differential equations

Abstract

Second derivative general linear methods (SGLMs) have been already implemented in a variable stepsize environment using Nordsieck technique. In this paper, we introduce variable stepsize SGLMs directly on nonuniform grid. By deriving the order conditions of the proposed methods of order p and stage order $q=p$, some explicit examples of these methods up to order four are given. By some numerical experiments, we show the efficiency of the proposed methods in solving nonstiff problems and confirm the theoretical order of convergence.

Introduction

Efficiency of the constructed methods for the numerical solution of initial value problem of ordinary differential equations (ODEs)

$$\{\begin{array}{c}{y}^{\prime }(x)=f(y(x)),x\in [x_0,X],\hfill \\ y(x_0)=y_0,\hfill \end{array}$$

where $f:{\mathbb{R}}^m\to {\mathbb{R}}^m$ is sufficiently smooth function and m is the dimensional of system, does depend on designing of their implementation issues. Indeed, reliable codes are developed based on the methods with desirable features and more sophisticated implementation strategies. Among all the classes of numerical methods, the class of general linear methods (GLMs)—due to the fact of their comprehensive structure—is one of the best. In fact, GLMs are multivalue-multistage methods which include linear multistep and Rungee–Kutta methods as special cases. These methods which were introduced by Butcher [9], [10] take the form

$$\{\begin{array}{c}{Y}_i^{[n+1]}=h\sum _{j=1}^s{a}_{ij}f(Y_j^{[n+1]})+\sum _{j=1}^r{u}_{ij}{y}_j^{[n]},i=1,\dots ,s,\hfill \\ y_i^{[n+1]}=h\sum _{j=1}^s{b}_{ij}f(Y_j^{[n+1]})+\sum _{j=1}^r{v}_{ij}{y}_j^{[n]},i=1,\dots ,r,\hfill \end{array}$$

where $n=0,1,\dots ,N-1$, $Nh=X-x_0$, s is the number of internal stages and r is the number of input and output values. Here, $Y_i^{[n+1]}$, $i=1,2,\dots ,s$, are an approximation of the stage order q to the solution y at the points $x_n+c_{i}h$, with $c={[c_1{c}_2\cdots c_s]}^T$ as the abscissa vector. Moreover, the values $y_i^{[n+1]}$, $i=1,2,\dots ,r$, are approximations of order p to the linear combinations of scaled derivatives of the solution y at the points $x_{n+1}$, i.e.,

$$y_i^{[n+1]}={\alpha }_{i0}y(x_{n+1})+{\alpha }_{i1}h{y}^{\prime }(x_{n+1})+{\alpha }_{i2}{h}^2{y}^{″}(x_{n+1})+\cdots +{\alpha }_{ip}{h}^p{y}^{(p)}(x_{n+1})+O(h^{p+1}),$$

for some real numbers ${\alpha }_{i1},{\alpha }_{i2},\dots ,{\alpha }_{ip}$. These values denote the computed quantities at step $n+1$, which also represent the incoming values for step $n+2$. Diagonally implicit multistage integration methods (DIMSIMs) as a large subclass of GLMs have been studied and discussed by Butcher and Jackiewicz, for instance in [11], [15], [17], [29]. Implementation issues for DIMSIMs have been studied in [8], [14], [16], [18] and efficient Matlab codes dim18.m and dim13s.m based on these methods using the Nordsieck technique have been developed in [12], [28]. Furthermore, another reliable code has been developed based on the Nordsieck GLMs in which the matrix $W:=[{\alpha }_{ij}]$ is equal to the identity matrix, i.e., the input and output vectors $y^{[n-1]}$ and $y^{[n]}$ approximate the Nordsieck vector $z(x,h)$ defined by

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

at the points $x=x_{n-1}$ and $x=x_n$, respectively [7]. Also, in [30], the authors studied a class of variable stepsize DIMSIMs which provides an alternative to the Nordsieck technique of changing the stepsize of integration. This approach is based on the derivation of the methods based directly on nonuniform grid so that the coefficients matrices $A:=[a_{ij}]$, $U:=[u_{ij}]$, $B:=[b_{ij}]$, $V:=[v_{ij}]$, the vector c as well as the matrix W depend, in general, on the ratios of the current stepsize and the past stepsizes.

To increase the chances of finding methods with high order accuracy together with desirable stability properties, higher derivatives, especially second derivative, of the solution can be used. Therefore, several researches have been focused on the construction of the methods incorporating the second derivative of the solution, see, for instance, [20], [21], [22], [26], [27]. Trying to construct second derivative methods in a unified structure and with efficient properties leads to the second derivative general linear methods (SGLMs) where were first introduced by Butcher and Hojjati in [13] and more investigated by Abdi et al. in [1], [2], [3], [4], [5], [6]. Denoting the function $g:=f^{\prime }(\cdot )f(\cdot )$ as the second derivative of the solution y and using of the previous notations, the SGLMs take the form

$$\{\begin{array}{c}{Y}_i^{[n+1]}=h\sum _{j=1}^s{a}_{ij}f(Y_j^{[n+1]})+h^2\sum _{j=1}^s{\bar{a}}_{ij}g(Y_j^{[n+1]})+\sum _{j=1}^r{u}_{ij}{y}_j^{[n]},i=1,\dots ,s,\hfill \\ y_i^{[n+1]}=h\sum _{j=1}^s{b}_{ij}f(Y_j^{[n+1]})+h^2\sum _{j=1}^s{\bar{b}}_{ij}g(Y_j^{[n+1]})+\sum _{j=1}^r{v}_{ij}{y}_j^{[n]},i=1,\dots ,r.\hfill \end{array}$$

The SGLM (1.3) has order p and stage order $q=p$ if and only if [5]

$$\begin{array}{c}\hfill C=ACK+\bar{A}C{K}^2+UW,\\ \hfill WE=BCK+\bar{B}C{K}^2+VW,\end{array}$$

where the matrices C, K, and E are determined by

$$C:=[ec\frac{c^2}{2!}\cdots \frac{c^p}{p!}]\in {\mathbb{R}}^{s\times (p+1)},K:=[0e_1{e}_2\cdots e_p]\in {\mathbb{R}}^{(p+1)\times (p+1)},E:=\exp (K)\in {\mathbb{R}}^{(p+1)\times (p+1)}$$

with $e={[11\cdots 1]}^T\in {\mathbb{R}}^s$ and $e_j$ as the jth vector of the canonical basis in ${\mathbb{R}}^{p+1}$.

Construction and implementation of the Nordsieck SGLMs have been discussed in [1], [6]. Furthermore, in [3], the authors have investigated the implementation of SGLMs in a variable stepsize environment using the Nordsieck technique and the practical Matlab code SGLM4.m based on an L–stable SGLM of order four has been developed.

Second derivative diagonally implicit multistage integration methods (SDIMSIMs) as a subclass of SGLMs have been introduced in [4] in which usually $p=q=r=s$ and the matrix V is a rank-one matrix. The spirit of this paper is to introduce SDIMSIMs directly on nonuniform grid. For such methods, it is not required to update the input vector for the new stepsize; actually, the output vector of the last step can be directly used in the next step as the input vector.

The rest of the paper is organized along the following lines. In Section 2, we formulate the SDIMSIMs (1.3) as variable stepsize methods directly on nonuniform grid. Then, the order conditions of these methods in the case that the stage order is equal to the order of the methods, are derived. Section 3 is devoted to the construction of such methods in a special class up to order four. To show the efficiency of the constructed methods, some numerical results are provided in Section 4.