Adaptive mesh selection asymptotically guarantees a prescribed local error for systems of initial value problems

TL;DR

This paper presents an adaptive mesh selection algorithm for initial value problems that asymptotically guarantees prescribed local errors, improving efficiency over uniform meshes through rigorous analysis and numerical validation.

Contribution

It introduces a theoretically grounded adaptive mesh selection method that guarantees local error bounds and demonstrates its advantages over uniform meshes.

Findings

Adaptive mesh selection asymptotically guarantees prescribed local errors.

The proposed algorithm improves efficiency compared to uniform mesh methods.

Numerical results confirm theoretical advantages of the adaptive approach.

Abstract

We study adaptive mesh selection for the solution of systems of initial value problems. The goal is a rigorous theoretical analysis of potential advantages of adaption. For an optimal method in the sense of the speed of convergence, we propose an algorithm for successive selection of the mesh points. The selection is based on an upper bound on the local error, and it (asymptotically) guarantees the local errors not exceeding a prescribed level. The mesh selection algorithm can be applied to a general class of methods, not only to the chosen one. We rigorously discuss the cost of the proposed algorithm, comparing it to other algorithms equipped with different mesh selection procedures. We specify a quantitative advantage of the adaptive mesh over the uniform one. Adjustment of the mesh points to a local behavior of the solution yields improved efficiency of the algorithm. Some numerical results illustrating theoretical findings are reported.

Equations

$$z^{\prime }(t)=f(t,z(t)),t\in [a,b],z(a)=\eta ,$$

$$\Vert f(t,y)-f(t,\bar{y})\Vert \le L\Vert y-\bar{y}\Vert \text{mbox}fort\in [a,b],y,\bar{y}\in {\mathbf{R}}^d.$$

$$D=\{y\in {\mathbf{R}}^d:\Vert y\Vert \le \underset{t\in [a,b]}{\sup }\Vert z(t)\Vert +1\}.$$

$$\underset{0\le i\le m-1}{\max }(x_{i+1,m}-x_{i,m})\le K\ell (m).$$

$$z_i^{\prime }(t)=f(t,z_i(t)),t\in [x_i,x_{i+1}],z_i(x_i)=y_i.$$

$$\underset{t\in [x_i,x_{i+1}]}{\sup }\Vert z_i(t)-l_i(t)\Vert ,i=0,1,\dots ,m-1.$$

$$q_{i,j}(t)=\sum _{k=0}^{r-1}{g}_{i,j}(t_k)\prod _{p=0,p\ne k}^{r-1}\frac{t-t_p}{t_k-t_p},t\in [x_i,x_{i+1}],$$

$$l_{i,j+1}(t)=y_i+\underset{x_i}{\overset{t}{\int }}{q}_{i,j}(\xi )d\xi ,t\in [x_i,x_{i+1}].$$

$$l_{r+1}(t)=l_{i,r+1}(t),t\in [x_i,x_{i+1}].$$

$$\underset{t\in [a,b]}{\sup }\Vert z(t)-l_{r+1}(t)\Vert \le M\underset{0\le j\le m-1}{\max }{h}_j^r,$$

$$\underset{t\in [a,x_i]}{\sup }\Vert z(t)-l_{r+1}(t)\Vert \le M_i\underset{0\le j\le i-1}{\max }{h}_j^r,$$

$$\Vert z(t)-z_i(t)\Vert \le \exp (L{h}_i)\Vert z(x_i)-y_i\Vert ,t\in [x_i,x_{i+1}].$$

$${\bar{q}}_i(t)=\sum _{k=0}^{r-1}f(t_k,z_i(t_k)))\prod _{p=0,p\ne k}^{r-1}\frac{t-t_p}{t_k-t_p},t\in [x_i,x_{i+1}].$$

$$z_i(t)=y_i+\underset{x_i}{\overset{t}{\int }}f(\xi ,z_i(\xi ))d\xi ,$$

$$\Vert z_i(t)-l_{i,j+1}(t)\Vert \le \underset{x_i}{\overset{t}{\int }}\Vert f(\xi ,z_i(\xi ))-{\bar{q}}_i(\xi )\Vert d\xi +\underset{x_i}{\overset{t}{\int }}\Vert {\bar{q}}_i(\xi )-q_{i,j}(\xi )\Vert d\xi .$$

$$\Vert f(\xi ,z_i(\xi ))-{\bar{q}}_i(\xi )\Vert \le \frac{1}{r!}\underset{\alpha \in [x_i,x_{i+1}]}{\sup }\Vert H_i^{(r)}(\alpha )\Vert \prod _{k=0}^{r-1}|\xi -t_k|,\xi \in [x_i,x_{i+1}].$$

$$\Vert {\bar{q}}_i(\xi )-q_{i,j}(\xi )\Vert \le \sum _{k=0}^{r-1}\Vert f(t_k,z_i(t_k))-f(t_k,l_{i,j}(t_k))\Vert \prod _{p=0,p\ne k}^{r-1}\left|\frac{\xi -t_p}{t_k-t_p}\right|,$$

$$\Vert {\bar{q}}_i(\xi )-q_{i,j}(\xi )\Vert \le L{\hat{C}}_r\underset{t\in [x_i,x_{i+1}]}{\sup }\Vert z_i(t)-l_{i,j}(t)\Vert ,\xi \in [x_i,x_{i+1}],$$

$$e_{i,j+1}\le \frac{1}{r!}\underset{\alpha \in [x_i,x_{i+1}]}{\sup }\Vert H_i^{(r)}(\alpha )\Vert h_i^{r+1}+h_{i}L{\hat{C}}_r{e}_{i,j}.$$

$$e_{i,j}\le \frac{1}{r!}\underset{\alpha \in [x_i,x_{i+1}]}{\sup }\Vert H_i^{(r)}(\alpha )\Vert h_i^{r+1}\frac{1-(h_{i}L{\hat{C}}_r{)}^j}{1-h_{i}L{\hat{C}}_r}+(h_{i}L{\hat{C}}_r{)}^j{e}_{i,0}.$$

$$e_{i,0}\le h_i\underset{\alpha \in [x_i,x_{i+1}]}{\sup }\Vert f(\alpha ,z_i(\alpha ))\Vert ,$$

$$e_{i,j}\le \frac{2}{r!}\underset{\alpha \in [x_i,x_{i+1}]}{\sup }\Vert H_i^{(r)}(\alpha )\Vert h_i^{r+1}+h_i^{j+1}(L{\hat{C}}_r{)}^j\underset{\alpha \in [x_i,x_{i+1}]}{\sup }\Vert H_i(\alpha )\Vert .$$

$$(\alpha ,z_i(\alpha ))\in [a,b]\times D,$$

$$e_{i,r+1}\le \frac{2}{r!}{D}_r{h}_i^{r+1}+h_i^{r+2}(L{\hat{C}}_r{)}^{r+1}{D}_0,$$

$$e_{i,r+1}\le \left(\frac{2}{r!}{D}_r+\frac{1}{2}\right)h_i^{r+1}.$$

$$\Vert z(t)-l_{r+1}(t)\Vert \le \Vert z(t)-z_i(t)\Vert +\Vert z_i(t)-l_{i,r+1}(t)\Vert \le \exp (L{h}_i)\Vert z(x_i)-y_i\Vert +e_{i,r+1}.$$

$$\Vert z(t)-l_{r+1}(t)\Vert \le \exp (L{h}_i)M_i\underset{0\le j\le i-1}{\max }{h}_j^r+\left(\frac{2}{r!}{D}_r+\frac{1}{2}\right)h_i^{r+1},t\in [x_i,x_{i+1}].$$

$$\underset{t\in [a,x_{i+1}]}{\sup }\Vert z(t)-l_{r+1}(t)\Vert \le M_{i+1}\underset{0\le j\le i}{\max }{h}_j^r,$$

$$e_{i,r+1}\le \frac{2}{r!}\underset{\alpha \in [x_i,x_{i+1}]}{\sup }\Vert H_i^{(r)}(\alpha )\Vert h_i^{r+1}+h_i^{r+2}(L{\hat{C}}_r{)}^{r+1}\underset{\alpha \in [x_i,x_{i+1}]}{\sup }\Vert H_i(\alpha )\Vert .$$

$$e_{i,r+1}\le 2\left(\frac{1}{r!}\underset{\alpha \in [x_i,x_{i+1}]}{\sup }\Vert H_i^{(r)}(\alpha )\Vert +\beta \right)h_i^{r+1}.$$

Full text

**Adaptive mesh selection asymptotically guarantees a prescribed local error for systems of initial value problems 111  This research was partly supported by the Polish NCN grant - decision No. DEC-2017/25/B/ST1/00945 and by the Polish Ministry of Science and Higher Education

**

**Bolesław Kacewicz 222

*AGH University of Science and Technology, Faculty of Applied Mathematics,

Al. Mickiewicza 30, paw. A3/A4, III p., pok. 301,

30-059 Cracow, Poland

E-mail B. Kacewicz:    kacewicz@agh.edu.pl* **

Abstract

We study potential advantages of adaptive mesh point selection for the solution of systems of initial value problems. For an optimal order discretization method, we propose an algorithm for successive selection of the mesh points, which only requires evaluations of the right-hand side function. The selection (asymptotically) guarantees that the maximum local error of the method does not exceed a prescribed level. The usage of the algorithm is not restricted to the chosen method; it can also be applied with any method from a general class. We provide a rigorous analysis of the cost of the proposed algorithm. It is shown that the cost is almost minimal, up to absolute constants, among all mesh selection algorithms. For illustration, we specify the advantage of the adaptive mesh over the uniform one. Efficiency of the adaptive algorithm results from automatic adjustment of the successive mesh points to the local behavior of the solution. Some numerical results illustrating theoretical findings are reported.

Mathematics Subject Classification: 65L05, 65L50, 65L70

1 Introduction

We deal with the solution of systems of initial value problems (IVPs)

[TABLE]

where $a<b$, $\eta\in{\bf R}^{d}$ and $f:[a,b]\times{\bf R}^{d}\to{\bf R}^{d}$ is a $C^{r}$ function. We study how much adaptive mesh points improve efficiency of algorithms for solving (1). For many years, adaption has been a standard tool in numerical packages. Well known examples include the package QUADPACK [8] for numerical integration, or, among other solvers, the DIFSUB procedure by C.W. Gear or the library ODEPACK by A. Hindmarsh for IVPs. Many authors have reported superiority of adaption over nonadaption for (1), based on numerical results for a number of computational examples. One can cite as an example papers such as [3], [6] or [7]. Practical efficiency gives us considerable knowledge about the power of adaption. Conclusions are however not complete; an analysis of theoretical aspects is missing in many cases of step size control strategies. For instance, most often, step size control devices are not supported by cost analysis. The use of variable step size not only improves the efficiency of methods for solving regular problems (1), but it also allows us to manage singularities. Considerable progress has been made in rigorous analysis of adaption in that case, see e.g. [5], [10]. Adaption allows us in many cases to maintain for singular problems the order of convergence known from the regular case.

Advantages of adaption for regular problems are of different type. Integration of scalar $C^{4}$ functions by the Simpson rule have been recently studied in [9]. It is shown in [9] that adaption does not improve the order of convergence, but it can reduce the asymptotic constant of the method. In the similar spirit, adaption has been considered for univariate approximation and minimization in [1]. For scalar autonomous problems (1), adaptive mesh selection has been recently studied in [4]. An adaptive strategy has been proposed and the cost analyzed, based on specific properties of scalar autonomous equations. The particular technique used does not allow us to extend the results from [4] to systems of IVPs.

In the present paper, we consider general systems (1), with a $C^{r}$ right-hand side function $f$. Our contribution can be summarized as follows:

• •

For a maximal order method for solving the system (1), we propose a new mesh selection algorithm that guarantees the local error of the method at a precribed level $\varepsilon$, for sufficiently small $\varepsilon\in(0,1)$. At each time step, we select the step size and compute an approximation to the solution, which requires two runs of the basic method. Information about the function $f$ only consists of function evaluations. We show that the mesh selection algorithm can be applied to a general class of methods for solving (1).

• •

We rigorously analyze the cost of the method equipped with the mesh selection procedure. We show that the cost is minimal among all mesh selection strategies, provided that we accept some absolute constants.

• •

We specify in a quantitative way an advantage of the adaptive mesh over the uniform one.

The adaptive mesh selection for systems (1) does not improve the speed of convergence of algorithms. A potential gain of adaption lies in reducing a coefficient in the error bound.

The paper is organized as follows. We formulate the problem in Section 2. An algorithm $\phi^{*}$ with maximal order is defined in Section 3, and the error analysis is given in Theorem  1. In Section 4, we discuss upper bounds on the local error of $\phi^{*}$, in particular, we show a constructive upper bound in Theorem 2. Section 5 contains the main algorithm ADAPT-MESH which combines $\phi^{*}$ with a new mesh selection algorithm. A generalization of the algorithm ADAPT-MESH is given in Section 6. Section 7 contains a cost analysis of the algorithm ADAPT-MESH compared to other algorithms, see Theorem 3. Possible advantage of the adaptive mesh over the uniform one is discussed in detail. Finally, some results of numerical experiments are reported in Section 8.

2 Problem formulation

We consider problem (1) with a continuous function $f:[a,b]\times{\bf R}^{d}\to{\bf R}^{d}$ such that for some $L\geq 0$

[TABLE]

Here and in what follows $\|\cdot\|$ denotes the maximum norm in ${\bf R}^{d}$. It follows that there exists a unique solution $z$ of (1) defined on $[a,b]$. Let $r\geq 1$. We assume that $f$ is a regular function in a subset of its domain, $f\in C^{r}([a,b]\times D)$, where

[TABLE]

The class of functions $f$ satisfying the above assumptions will be denoted by $F_{r}$.

Let $m\in{\bf N}$. We wish to compute an approximation to the solution $z$ in $[a,b]$. For $m+1$ mesh points $a=x_{0,m}<x_{1,m}<\ldots<x_{m,m}=b$, we do it by computing approximations $l_{i}$ to $z$ in the subintervals $[x_{i,m},x_{i+1,m}]$, $i=0,1,\ldots,m-1$.

Let ${\ell}(m)$ be any nonincreasing sequence convergent to [math] as $m\to\infty$. We consider for any $f$ a class of partitions of $[a,b]$ defined as follows. We assume that there exist $K=K(f,a,b,\eta)$ and $k_{0}=k_{0}(f,a,b,\eta)$ such that for all $m\geq k_{0}$ and any partition it holds

[TABLE]

Thus, the partitions under consideration are uniformly normal. Note that we always have $\max\limits_{0\leq i\leq m-1}(x_{i+1,m}-x_{i,m})\geq(b-a)/m$ for any $m\geq 1$. Thus, the condition (4) implies that ${\ell}(m)$ cannot go to zero faster than $1/m$. Note that the convergence of ${\ell}(m)$ to zero can be arbitrarily slow, and the constant $K$ can be arbitrarily large. We shall omit in the sequel the second subscript $m$, remembering that the choice of points $x_{i}$ can be different for varying $m$.

For a given $y_{i}\in{\bf R}^{d}$, we denote by $z_{i}$ the solution of a local problem

[TABLE]

If $l_{i}$ is an approximation to $z_{i}$ given by a certain method, then local errors of the method are given by

[TABLE]

Our aim is to select possibly small $m$ and mesh points $\{x_{i}\}_{i=0}^{m}$ such that the local errors remain at a precribed level $\varepsilon>0$.

3 The basic method $\phi^{*}$ and its error

The basic method makes use of the approximate Picard iteration, an idea that turned out useful in several contexts, see e.g. [2], [5]. Let $m\in{\bf N}$ and $x_{0}=a<x_{1}<\ldots<x_{m}=b$ be mesh points satisfying (4). Let $y_{0}=\eta$. For a given $y_{i}$, we define approximations $l_{i,j}$ in $[x_{i},x_{i+1}]$ as follows.

We set $l_{i,0}(t)\equiv y_{i}$. Let $l_{i,j}$ be given. Denote by $t_{0},t_{1},\ldots,t_{r-1}$ the equidistant nodes in $[x_{i},x_{i+1}]$, with $t_{0}=x_{i}$ for $r=1$ and $t_{0}=x_{i}$, $t_{r-1}=x_{i+1}$ for $r\geq 2$. (The points $t_{k}$ depend on $i$; we shall omit this index in the notation.) We define $q_{i,j}$ to be the Lagrange interpolation polynomial of degree at most $r-1$ for the function $g_{i,j}(t)=f(t,l_{i,j}(t))$ based on the nodes $t_{0},t_{1},\ldots,t_{r-1}$,

[TABLE]

where $\prod\limits_{p=0,p\neq k}^{0}=1$. An approximation $l_{i,j+1}$ is given by

[TABLE]

The final approximation in $[x_{i},x_{i+1}]$ is given by $l_{i,r+1}$. To complete the definition, we set $y_{i+1}=l_{i,r+1}(x_{i+1})$.

For $t\in[a,b]$ we define a continuous approximation to $z$ by

[TABLE]

The transformation that assignes to $f$ the approximation $l_{r+1}$ will be denoted by $\phi^{*}$, $\phi^{*}(f)(t)=l_{r+1}(t)$, $t\in[a,b].$

The following theorem provides error analysis of the method $\phi^{*}$. The proof follows usual lines of the analysis of approximate Picard iteration, it is however focused on our specific requirements. We shall need in the next sections the error bound (9) for a non-uniform mesh, as well as specific local error bounds derived in the body of the proof. Let $h_{i}=x_{i+1}-x_{i}$.

**Theorem 1  ** Let $f\in F_{r}$. There exists $m_{0}$ such that for all $m\geq m_{0}$ and any $\{x_{i}\}_{i=0}^{m}$ satisfying (4), the global error of $\phi^{*}$ at $f$ satisfies

[TABLE]

*where $M=\exp(L(b-a))(b-a)\left(2D_{r}/r!+1/2\right)$, and the number $D_{r}$, only dependent on $f$, $r$, $a$, $b$, is defined below before the inequality (17). *

**Proof  ** Let $M_{0}=0$ and $M_{i+1}=\exp(Lh_{i})M_{i}+\left(2D_{r}/r!+1/2\right)h_{i}$, $i=0,1,\ldots,m-1$. One can check that $M_{i}\leq M$, $i=0,1,\ldots,m$. We shall show by induction on $i$ that

[TABLE]

where $\max\limits_{0\leq j\leq-1}=1$.

For $i=0$, (10) holds true. Let (10) hold for some $i$. Consider the interval $[x_{i},x_{i+1}]$. We have that

[TABLE]

We shall now study the local error in $[x_{i},x_{i+1}]$ given by $e_{i,j}=\sup\limits_{t\in[x_{i},x_{i+1}]}\|z_{i}(t)-l_{i,j}(t)\|$. Denoting $H_{i}(t)=f(t,z_{i}(t))\;(=z_{i}^{\prime}(t))$, we let $\bar{q}_{i}$ be the Lagrange interpolation polynomial for $H_{i}$,

[TABLE]

Since

[TABLE]

we have that

[TABLE]

By the Lagrange interpolation error formula applied component by component, we have that

[TABLE]

Furthermore,

[TABLE]

which yields that

[TABLE]

where $\hat{C}_{r}$ only depends on $r$. From (13) we get for $j=0,1,\ldots$

[TABLE]

By solving (14) we get for $j=0,1,\ldots$

[TABLE]

Since

[TABLE]

for $m$ sufficiently large (such that $h_{i}L\hat{C}_{r}\leq 1/2$) we have that

[TABLE]

Note that $H_{i}^{(r)}(\alpha)$ and $H_{i}(\alpha)$ can be expressed in terms of partial derivatives of the function $f$ of order $0,1,\ldots,r$, evaluated at $(\alpha,z_{i}(\alpha))$. Due to (11) and the inductive assumption, for sufficiently large $m$ we have that $\|z(\alpha)-z_{i}(\alpha)\|\leq 1$ for $\alpha\in[x_{i},x_{i+1}]$. Hence,

[TABLE]

where the set $D$ is given in (3). This yields that $\sup\limits_{\alpha\in[x_{i},x_{i+1}]}\|H_{i}^{(r)}(\alpha)\|$ and $\sup\limits_{\alpha\in[x_{i},x_{i+1}]}\|H_{i}(\alpha)\|$ are bounded from above independently of $i$ and $m$ by some numbers $D_{r}$ and $D_{0}$, respectively. For the final approximation $l_{i,r+1}$ we have from (15) that

[TABLE]

which yields for sufficiently large $m$ that

[TABLE]

For the final global approximation $l_{r+1}$, we have for $t\in[x_{i},x_{i+1}]$

[TABLE]

By the inductive assumption, we get that

[TABLE]

Hence,

[TABLE]

where $M_{i+1}=\exp(Lh_{i})M_{i}+\left(2D_{r}/r!+1/2\right)h_{i}$. The induction is finished.

To complete the proof we recall that $M_{i}\leq M$ for $i=0,1,\ldots,m,$ where $M$ is given in the statement of the theorem.

Given $m$ and a mesh $\{x_{i}\}_{i=0}^{m}$, the method $\phi^{*}$ describes the construction of the approximations $l_{i,r+1}$ in $[x_{i},x_{i+1}]$ for $i=0,1,\ldots,m-1$. Theorem 1 provides the bound on the global error of $\phi^{*}$. A selection of the mesh $\{x_{i}\}_{i=0}^{m}$ still remains an open question; we will study this issue in the next sections.

4 **Local error bounds for $\phi^{*}$ **

We now extract from the proof of Theorem 1 bounds on the local error $e_{i,r+1}$ of the method $\phi^{*}$. By (15)

[TABLE]

Let $\beta>0$. From (22), for sufficiently large $m$ we have

[TABLE]

The function $H_{i}(\alpha)=f(\alpha,z_{i}(\alpha))=z_{i}^{\prime}(\alpha)$ is not known. We now show how the term ’$\sup$’ above can be (asymptotically) replaced by a known quantity.

We take $\bar{x}_{i+1}$, $x_{i}<\bar{x}_{i+1}\leq b$. For $\bar{h}_{i}=\bar{x}_{i+1}-x_{i}$ we assume that

[TABLE]

where $\gamma\geq 1$ is a given number which may depend on $f$, but is independent of $i$ and $m$.

Let $\bar{t}_{k}$, $k=0,1,\ldots,r$ be equidistant points from $[x_{i},\bar{x}_{i+1}]$, $\bar{t}_{0}=x_{i}$, $\bar{t}_{r}=\bar{x}_{i+1}$ (we omit the index $i$). We construct an auxillary approximation $\bar{l}_{i,r+1}$ in the interval $[x_{i},\bar{x}_{i+1}]$ in the same way as we did in (6)–(8) in the case of the approximation $l_{i,r+1}$ in the interval $[x_{i},x_{i+1}]$, using now as interpolation nodes the points $\bar{t}_{0},\bar{t}_{1},\ldots,\bar{t}_{r-1}\in[x_{i},\bar{x}_{i+1}]$.

Let $\tilde{H}_{i}(\alpha)=f(\alpha,\bar{l}_{i,r+1}(\alpha))$ and $\tilde{H}_{i}[\bar{t}_{0},\bar{t}_{1},\ldots,\bar{t}_{r}]$ be the divided difference, computed component by component, for $\tilde{H}_{i}$. We shall need the bounds stated in the following two lemmas. Recall that $H_{i}(t)=f(t,z_{i}(t))$.

**Lemma 1  ** Let $f\in F_{r}$, $\beta>0$ and $\varphi\in(0,1)$. There exists $m_{0}$ such that for any $m\geq m_{0}$, for any $\{x_{i}\}_{i=0}^{m}$ satisfying (4) and $i=0,1,\ldots,m-1$ we have

[TABLE]

*where $S_{i}=\left\|H_{i}[\bar{t}_{0},\bar{t}_{1},\ldots,\bar{t}_{r}]\right\|$. *

**Proof  ** Let $\bar{\alpha}$ be a point from $[x_{i},x_{i+1}]$ for which

[TABLE]

For the $l$th component $H_{i}^{l}$ of the function $H_{i}$, we express the divided difference as

[TABLE]

where $\hat{\alpha}^{l}$ is some point from $[x_{i},\bar{x}_{i+1}]$. For any $l=1,2,\ldots,d$ it holds

[TABLE]

where

[TABLE]

Similarly to what we have already noticed, the quantity $\left(H_{i}^{l}\right)^{(r)}(t),\;t\in[x_{i},\max\{x_{i+1},\bar{x}_{i+1}\}],$ can be expressed by values of a continuous function defined by partial derivatives of $f$, evaluated at $(t,z_{i}(t))$, where the argument $(t,z_{i}(t))$ belongs to the compact set $[a,b]\times D$. By the uniform continuity of this function, we have that

[TABLE]

Hence $|\kappa_{i}^{l}|\leq\varphi$ for $m\geq\bar{m}_{0}$, which leads to (25).

**Lemma 2  ** Let $f\in F_{r}$. There exist $\bar{C}$, $\tilde{C}$, $m_{0}$ such that for any $m\geq m_{0}$, for any $\{x_{i}\}_{i=0}^{m}$ satisfying (4) it holds

[TABLE]

*for $i=0,1,\ldots,m-1$.

(Here $\bar{t}_{k}$ are, as above, the equidistant nodes in $[x_{i},\bar{x}_{i+1}]$; the index $i$ is omitted.) *

**Proof  ** The proof follows from the fact that

[TABLE]

and from (18) applied to the approximation $\bar{l}_{i,r+1}$ in the interval $[x_{i},\bar{x}_{i+1}]$ with $h_{i}$ replaced by $\bar{h}_{i}$.

From (23), Lemmas 1 and 2 we get the following computable (asymptotic) upper bound on the local error of the method $\phi^{*}$.

**Theorem 2  ** Let $f\in F_{r}$, $\beta>0$ and $\varphi\in(0,1)$. There exists $m_{0}$ such that for all $m\geq m_{0}$, for any $\{x_{i}\}_{i=0}^{m}$ satisfying (4), and any $\bar{x}_{i+1}$ satisfying (24) it holds

[TABLE]

where

[TABLE]

**Proof  ** We successively use (23), (25) and (27). We first get

[TABLE]

and next

[TABLE]

For sufficiently large $m$ such that $\tilde{C}\bar{h}_{i}\leq\beta/3$ we get (28).

For given $x_{i}<b$ and $y_{i}$, after selecting a point $\bar{x}_{i+1}$, we are able to construct an auxillary approximation $\bar{l}_{i,r+1}$ and compute $\tilde{H}_{i}[\bar{t}_{0},\bar{t}_{1},\ldots,\bar{t}_{r}]$. Hence, for given $\beta$ and $\varphi$, the coefficient $G_{i}$ can be effectively computed. For further purposes, note that

[TABLE]

On the other hand, due to (27), we have in terms of $H_{i}$ that

[TABLE]

for sufficiently large $m$. Due to (26) and the observation made in the proof of Lemma 1 regarding the derivatives of $H_{i}$, we have the bound

[TABLE]

where $N(f)$ is independent of $i$ and $m$.

5 Algorithm with guaranteed local error

Let $\varepsilon\in(0,1)$. Our aim is to select mesh points $\{x_{i}^{*}\}$ in the method $\phi^{*}$ in order to guarantee that in all steps the local error is at most $\varepsilon$,

[TABLE]

Given $x_{i}^{*}$ and $y_{i}^{*}$, we set

[TABLE]

where $h:(0,1)\to{\bf R}_{+}$ is any function such that $h(\varepsilon)=O(\varepsilon^{1/(r+1)})$, as $\varepsilon\to 0^{+}$, with a constant in the ${}^{\prime}O^{\prime}$-notation possibly dependent on $f$ (but not on $i$ or the number of subintervals). We discuss the form of $h(\varepsilon)$ below. We then compute the auxillary approximation $\bar{l}_{i,r+1}$ defined in Section 4 before Lemma 1, and the coefficient $G_{i}$ from (29), with $[x_{i},\bar{x}_{i+1}]$ replaced by $[x_{i}^{*},\bar{x}_{i+1}]$. The mesh point $x_{i+1}^{*}\leq b$ is now selected such that

[TABLE]

that is, we put

[TABLE]

With this mesh point $x_{i+1}^{*}$, the approximation $l_{i,r+1}$ defined in the algorithm $\phi^{*}$ satisfies, due to Theorem 2, the condition (32).

**Remark 1  ** We comment on the choice of $h(\varepsilon)$. We can choose $\bar{x}_{i+1}$ ’large’ by taking $h(\varepsilon)=\varepsilon^{1/(r+1)}$, so that

[TABLE]

The condition (24) holds for $\bar{x}_{i+1}$ and $x_{i+1}^{*}$, since $x_{i+1}^{*}-x_{i}^{*}=\Theta(\bar{x}_{i+1}-x_{i}^{*})$, which follows from (30) and (31). If we take $h(\varepsilon)=O(\varepsilon)$, then the condition (24) holds for $\bar{x}_{i+1}$ and $x_{i+1}^{*}$, since $\bar{x}_{i+1}\leq x_{i+1}^{*}$ for sufficiently small $\varepsilon$. In this case $\bar{x}_{i+1}$ can be arbitrarily close to $x_{i}^{*}$. Our results hold for such $h(\varepsilon)$. It seems that the faster $h(\varepsilon)$ goes to zero with $\varepsilon\to 0$, the earlier the asymptotics shows up.

We have arrived at the following algorithm.

Algorithm ADAPT-MESH

**1.  ** Choose $\varepsilon\in(0,1)$, $\beta>0$ and $\varphi\in(0,1)$. Set $x_{0}^{*}=a$ and $y_{0}^{*}=\eta$.

**2.  ** Given $x_{i}^{*}$ and $y_{i}^{*}$, compute $\bar{x}_{i+1}$ from (33). Compute the equidistant points $\bar{t}_{0}=x_{i}^{*},\bar{t}_{1},\ldots,\bar{t}_{r}=\bar{x}_{i+1}$ from $[x_{i}^{*},\bar{x}_{i+1}]$.

**3.  ** Compute the approximation $\bar{l}_{i,r+1}$ in $[x_{i}^{*},\bar{x}_{i+1}]$ from (6) and (7), based on $\bar{t}_{0},\bar{t}_{1},\ldots,\bar{t}_{r-1}$.

**4.  ** Compute (component by component) the divided difference $\tilde{H}_{i}[\bar{t}_{0},\bar{t}_{1},\ldots,\bar{t}_{r}]$, where $\tilde{H}_{i}(t)=f(t,\bar{l}_{i,r+1}(t))$.

**5.  ** Compute $G_{i}$ from (29).

**6.  ** Compute $x_{i+1}^{*}$ from (35), and the equidistant points from $[x_{i}^{*},x_{i+1}^{*}]$ with $t_{0}=x_{i}^{*}$ for $r=1$ and $t_{0}=x_{i}^{*},t_{1},\ldots,t_{r-1}=x_{i+1}^{*}$ for $r\geq 2$.

**7.  ** Compute the approximation $l_{i,r+1}$ in $[x_{i}^{*},x_{i+1}^{*}]$ from (6) and (7), based on $t_{0},t_{1},\ldots,t_{r-1}$.

**8.  ** Set $y_{i+1}^{*}=l_{i,r+1}(x_{i+1}^{*})$. If $x_{i+1}^{*}<b$, then go to 2 with $i:=i+1$.

STOP

Steps 2–5 can be viewed as a ’prediction’ stage, while the final approximation is computed in steps 6 and 7.

We denote by $m^{*}=m^{*}(\varepsilon)$ the number of subintervals defined by the algorithm ADAPT-MESH. The local error of the approximation $l_{i,r+1}$ in the interval $[x_{i}^{*},x_{i+1}^{*}]$ is guaranteed to be at most $\varepsilon$,

[TABLE]

for sufficiently small $\varepsilon$, see (32). The cost of each step of the method $\phi^{*}$, when applied on a given mesh, is $\kappa^{*}(r)=r^{2}+\Theta(r)$ evaluations of $f$ which are needed to produce $l_{i,r+1}$. The algorithm ADAPT-MESH is additionally equipped with the mesh selection procedure which makes the cost twice as large. The cost of ADAPT-MESH equals $2r^{2}+\Theta(r)$ function evaluations per step, which is roughly $2m^{*}\kappa^{*}(r)$ in total. In Section 7 we compare it with the cost of other methods and mesh selection procedures.

By the definition of $\{x_{i}^{*}\}$, we have that $x_{i+1}^{*}-x_{i}^{*}=\left(\varepsilon/G_{i}\right)^{1/(r+1)}$, $i=0,1,\ldots,m^{*}-2$, and $x_{m^{*}-1}^{*}+\left(\varepsilon/G_{m^{*}-1}\right)^{1/(r+1)}\geq b$. This yields that $m^{*}$ is the minimal number $m\in{\bf N}$ such that

[TABLE]

where

[TABLE]

Taking into account the bounds on $G_{i}$, we have that

[TABLE]

We see that the mesh selection procedure in the algorithm ADAPT-MESH does not reduce the speed of growth of the cost as $\varepsilon\to 0$, with respect to the equidistant mesh. As in the latter case (see Section 7), we have that

[TABLE]

A potential gain of adaption lies in reducing the coefficient.

Note that the condition (4) is satisfied for $\{x_{i}^{*}\}$. Indeed, we have using (39) that

[TABLE]

for sufficiently small $\varepsilon$. Hence, (4) holds with any

[TABLE]

6 **Mesh selection for a general class of methods **

The mesh selection procedure described above can be applied to a class of methods $\phi$ for solving (1), not only for $\phi^{*}$. We assume that for any discretization $\{x_{i}\}_{i=0}^{m}$ satisfying (4), a method $\phi$ successively computes in each interval $[x_{i},x_{i+1}]$ an approximation $l_{i}$ to $z_{i}$, with $l_{i}(x_{i})=y_{i}$ and $y_{i+1}=l_{i}(x_{i+1})$, starting from $x_{0}=a$, $y_{0}=\eta$. Global approximation $l$ computed by $\phi$ in $[a,b]$ is composed of the approximations $l_{i}$ in $[x_{i},x_{i+1}]$, $\phi(f)(t)=l(t)=l_{i}(t)$, $t\in[x_{i},x_{i+1}]$, $i=0,1,\ldots,m-1$. We assume that the computation of $l_{i}$ requires a certain number of evaluations of some functionals on $f$ (information about $f$). The total number of evaluations in a single interval $[x_{i},x_{i+1}]$ is $\kappa_{\phi}(r)$, where $\kappa_{\phi}(r)$ is independent of $i$ and $m$. For instance, for the method $\phi^{*}$ the functionals are defined by evaluations of $f$, and $\kappa_{\phi^{*}}(r)=2r^{2}+\Theta(r)$. We assume about $\phi$ that:

**A.  ** There are $\bar{\beta},\beta>0$ such that for any $f\in F_{r}$ there is $m_{0}$ such that for all $m\geq m_{0}$, for any $\{x_{i}\}_{i=0}^{m}$ satisfying (4)

[TABLE]

Assumption A has been verified for $\phi=\phi^{*}$ in (23). Of course, the method $\phi^{*}$ is not the only example of $\phi$. It can also be defined in many different ways, e.g., by Taylor’s approximation.

**Remark 2  ** It is easy to see that Theorem 1 (with slightly different constant $M$), Lemma 1, Lemma 2 and Theorem 2 hold for $\phi$, with $\beta$ given in assumption A, with $l_{i,r+1}$ replaced by $l_{i}$ and $\bar{l}_{i,r+1}$ replaced by $\bar{l}_{i}$. The coefficient $G_{i}$ is now given by

[TABLE]

We have that

[TABLE]

where $N(f)$ depends on $\bar{\beta}$, and it is independent of $i$ and $m$.

We now discuss yet another local error bound for $\phi$. We show that the local solution $z_{i}$ in the bound (42) can be replaced, at cost of changing a constant, by the global solution $z$.

**Lemma 3  ** For any $f\in F_{r}$, $\beta>0$ and $\varphi\in(0,1)$ there is $m_{0}$ such that for all $m\geq m_{0}$, for any $\{x_{i}\}_{i=0}^{m}$ satisfying (4) it holds

[TABLE]

$i=0,1,\ldots,m-1$*, for some $\kappa_{i}$, where $|\kappa_{i}|\leq\varphi$. *

**Proof  ** Let the two sup above be achieved in points $t_{1},t_{2}\in[x_{i},x_{i+1}]$, respectively. We have

[TABLE]

where $\kappa_{i}=(1/r!)\left(\|z_{i}^{(r+1)}(t_{1})\|-\|z^{(r+1)}(t_{2})\|\right)/\left((1/r!)\|z^{(r+1)}(t_{2})\|+\beta\right)$. Hence,

[TABLE]

The last term is bounded by $\sup\limits_{|t_{1}-t_{2}|\leq h_{i}}\|z^{(r+1)}(t_{1})-z^{(r+1)}(t_{2})\|$, so that, due to the uniform continuity of $z^{(r+1)}$, it tends to [math] (uniformly with respect to $i$) as $m\to\infty$. Note that $z_{i}^{(r+1)}(t_{1})$ and $z^{(r+1)}(t_{1})$ can be expressed by partial derivatives of $f$ of order $0,1,\ldots,r$ evaluated in $(t_{1},z_{i}(t_{1}))$ and $(t_{1},z(t_{1}))$, respectively. Due to Theorem 1 for $\phi$, we have

[TABLE]

Hence, $z(t_{1}),z_{i}(t_{1})\in D$, see (3), for sufficiently large $m$. This and the uniform continuity of the partial derivatives of $f$ in $[a,b]\times D$ yield that the first term also tends to [math] as $m\to\infty$, uniformly with respect to $t_{1}$ and $i$. Hence, $|\kappa_{i}|$ tends to [math] as $m\to\infty$, uniformly with respect to $t_{1}$, $t_{2}$ and $i$. This proves the lemma.

We now list upper bounds on the local error of $\phi$ that appeared so far. The basic one is given in assumption A

[TABLE]

where

[TABLE]

The bound (45) involves the local solution $z_{i}$; it usually appears in the error analysis of a method $\phi$. The second one follows from Lemma 3 and has the form (45) with $c_{i}$ replaced by

[TABLE]

We observe that the bounds $\bar{c}_{i}\,h_{i}^{r+1}$, $i=0,1,\ldots,m-1$, only depend on the local behavior of $f$ and on the mesh $\{x_{i}\}$. They hold for any method $\phi$ satisfying A, and are useful for theoretical reasons. Note that the function $p(x_{i},x_{i+1})=\bar{c}_{i}\,h_{i}^{r+1}$ is an increasing function with respect to $x_{i+1}$ (for fixed $x_{i}$), and a decreasing function with respect to $x_{i}$ (for fixed $x_{i+1}$).

The third bound is constructive and it will be used in the algorithm ADAPT-MESH-GEN below. It is given by (45) with $c_{i}$ replaced by $G_{i}$ from (43). Note that the coefficients (47) and (43) are not overestimated compared to (46); they are equivalent up to a coefficient only dependent on $\varphi$, for sufficiently large $m$. It follows from Lemmas 1, 2 and 3 that for any $m\geq m_{0}$, any $\{x_{i}\}_{i=0}^{m}$ satisfying (4) and $i=0,1,\ldots,m-1$ it holds

[TABLE]

Hence, for a given $\varphi$, the three error bounds

[TABLE]

are equivalent up to absolute constants (for a fixed $\varphi$). In particular, they all reflect a local behavior of $f$.

As in the case of the method $\phi^{*}$, we are free to choose the mesh points for $\phi$. The following algorithm, very much similar to ADAPT-MESH, describes the mesh selection for $\phi$ that allows us to keep the local error at level $\varepsilon$.

Algorithm ADAPT-MESH-GEN

**1.  ** Choose $\varepsilon\in(0,1)$, and $\varphi\in(0,1)$. Set $x_{0}=a$ and $y_{0}=\eta$.

**2.  ** Given $x_{i}$ and $y_{i}$, compute

[TABLE]

**3.  ** Compute an approximation $\bar{l}_{i}$ in $[x_{i},\bar{x}_{i+1}]$ using $\phi$.

**4.  ** For $\tilde{H}_{i}(t)=f(t,\bar{l}_{i}(t))$, compute (component by component) the divided difference $\tilde{H}_{i}[\bar{t}_{0},\bar{t}_{1},\ldots,\bar{t}_{r}]$, where $\bar{t}_{0}=x_{i},\bar{t}_{1},\ldots,\bar{t}_{r}=\bar{x}_{i+1}$ are the equidistant points from $[x_{i},\bar{x}_{i+1}]$.

**5.  ** Compute $G_{i}=(4/3)\bar{\beta}\left(\|\tilde{H}_{i}[\bar{t}_{0},\bar{t}_{1},\ldots,\bar{t}_{r}]\|+\beta\right)(1+\varphi)$.

**6.  ** Compute

[TABLE]

**7.  ** Compute an approximation $l_{i}$ in $[x_{i},x_{i+1}]$ using $\phi$.

**8.  ** Set $y_{i+1}=l_{i}(x_{i+1})$. If $x_{i+1}<b$, then go to 2 with $i:=i+1$.

STOP

Since $c_{i}\,h_{i}^{r+1}\leq G_{i}\,h_{i}^{r+1}\leq\varepsilon$, we see that the local error of the approximation $l_{i}$ computed in the step 7 is guaranteed to be at most $\varepsilon$,

[TABLE]

for sufficiently small $\varepsilon$. The cost of each step of ADAPT-MESH-GEN, measured by the number of evaluations of functionals on $f$ needed to produce $l_{i}$, is doubled with respect to the cost of $\phi$ applied on a given mesh.

7 Adaptive mesh – the cost analysis

Consider an arbitrary method $\phi$ satisfying A, based on a mesh $x_{0}=a<x_{1}<\ldots<x_{m}=b$ for which (4) holds. We measure the cost of $\phi$, ${\rm cost}(\phi,m)$, by the total number of evaluations of functionals on $f$ needed for computing $l$ in all subintervals $[x_{i},x_{i+1}]$ , i.e.,

[TABLE]

Our goal is to keep the local error at a prescribed level $\varepsilon$, see (51). Assuming that the goal is achieved by some $\phi$ with some mesh $x_{0},x_{1},\ldots,x_{m}$, we wish to compare ${\rm cost}(\phi,m)$ with the cost of the algorithm ADAPT-MESH.

To compare the costs of algorithms, we use the local error bound $c_{i}\,h_{i}^{r+1}$ from assumption A. We wish to assure that

[TABLE]

which implies (51). We define the reference quantity $\hat{m}(\varepsilon)$ as follows. Let

[TABLE]

Then we define

[TABLE]

Thus, $\hat{m}(\varepsilon)$ is the minimal number of subintervals $m$ for which there exists a mesh with $m+1$ points such that $c_{i}\,h_{i}^{r+1}\leq\varepsilon$, $i=0,1,\ldots,m-1$.

Define similarly to (54) and (55) the following, technically useful, quantities

[TABLE]

and

[TABLE]

Since $\bar{c}_{i}\,h_{i}^{r+1}$ is an increasing function of $x_{i+1}$ (for fixed $x_{i}$) and a decreasing function of $x_{i}$ (for fixed $x_{i+1}$), for sufficiently small $\varepsilon>0$ the quantity $\bar{m}(\varepsilon)$ can be computed as follows. We start with $\hat{x}_{0}=a$, and for a given $\hat{x}_{i}$, we compute $\hat{x}_{i+1}$ as the unique solution of $\bar{c}_{i}\,(x_{i+1}-\hat{x}_{i})^{r+1}=\varepsilon$. Then $\bar{m}(\varepsilon)$ is the minimal $i$ such that $\hat{x}_{i}\geq b$.

Note further that for any $m\geq m_{0}$ and $\{x_{i}\}_{i=0}^{m}$ satisfying (4), it follows from (48) that

[TABLE]

which yields for $m\geq m_{0}$ that

[TABLE]

Hence, for any $\varepsilon\in(0,\varepsilon_{0}]$ (and fixed $\phi$, $\bar{\beta}$)

[TABLE]

We now compare the cost of ADAPT-MESH with the cost of other algorithms $\phi$ equipped with any mesh selection procedure. The number of subintervals computed by ADAPT-MESH is $m^{*}(\varepsilon)$ and the cost of producing $l_{i,r+1}$ in each subinterval is $\kappa^{*}(r)$. Since the cost is doubled due to the mesh selection, it holds

[TABLE]

The quantities $\hat{k}(m)$ and $\hat{m}(\varepsilon)$ depend on $\phi$ (and $\bar{\beta}$). In the following result, we shall use for clarity the notation $\hat{k}_{\phi}(m)$ and $\hat{m}_{\phi}(\varepsilon)$. We compare $m^{*}(\varepsilon)$ from ADAPT-MESH (where $\bar{\beta}=2$) with the minimal number of intervals for any other method $\phi$ with $\bar{\beta}=2$, equipped with any mesh selection strategy. We have

**Theorem 3  ** Let $f\in F_{r}$, $\varphi\in(0,1)$ and $\phi$ be any method satisfying A with $\bar{\beta}=2$. Then there exists $\varepsilon_{0}\in(0,1)$ such that for any $0<\varepsilon\leq\varepsilon_{0}$ it holds

[TABLE]

*with $p_{1}=(1+\varphi)/(1-\varphi)$ and $p=(1-\varphi)^{2}/(2(1+\varphi)^{2})$.

Hence,*

[TABLE]

**Proof  ** The algorithm ADAPT-MESH defines $m^{*}+1$ points $x_{i}^{*}$ such that for $h_{i}=x_{i+1}^{*}-x_{i}^{*}$ and $G_{i}$ for $\phi^{*}$ we have

[TABLE]

We show the lower bound in (61). For $\{x_{i}^{*}\}_{i=0}^{m^{*}}$ and $\varepsilon$ sufficiently small, we have that $c_{i}\,h_{i}^{r+1}\leq G_{i}\,h_{i}^{r+1}\leq\varepsilon$ for $i=0,1,\ldots,m^{*}-1$. This yields that $\hat{k}_{\phi^{*}}(m^{*})\leq\varepsilon$, which implies that $\hat{m}_{\phi^{*}}(\varepsilon)\leq m^{*}(\varepsilon)$. The further lower bound follows from (59).

We now show the upper bound. By (48), for any $m\geq m_{0}$, any $\phi$, any $\{x_{i}\}_{i=0}^{m}$ satisfying (4) and any $i$, we have for $h_{i}=x_{i+1}-x_{i}$ (and $G_{i}$ for $\phi$) that

[TABLE]

For $\phi=\phi^{*}$, $m=m^{*}$ and the mesh $\{x_{i}^{*}\}$ given by ADAPT-MESH, due to (63), we have for $h_{i}=x_{i+1}^{*}-x_{i}^{*}$ and $\varepsilon$ sufficiently small that

[TABLE]

that is,

[TABLE]

The observation made after (57) yields that the points $\hat{x}_{i}$ computed for the accuracy $(1-\varphi)/(2(1+\varphi))\varepsilon$ satisfy in the light of (64) the inequalities $\hat{x}_{i}\leq x_{i}^{*}$ for $i=0,1,\ldots,m^{*}-1$. This implies that

[TABLE]

Finally, (59) and (65) give us the desired inequality (61)

[TABLE]

The inequalities for the cost follow immediately.

Theorem 3 says that the cost of ADAPT-MESH can only exceed by the constant $2\kappa^{*}(r)/\kappa_{\phi}(r)$ the cost of any algorithm $\phi$ with $\bar{\beta}=2$, with any mesh selection strategy such that the local error is at most $(1-\varphi)^{2}/(2(1+\varphi)^{2})\,\varepsilon$. This accuracy is more demanding than $\varepsilon$. For instance, if we take $\varphi=1/2$, then $\varepsilon$ in the accuracy demand is replaced by $\varepsilon/18$.

Observe that adaptive mesh cannot reduce the speed of growth of the cost as $\varepsilon\to 0$. It follows from (61) and (40) that for the best choice of points we have, similarly as for the equidistant mesh, that

[TABLE]

The asymptotics is thus the same. Possible advantage of the adaptive mesh selection is hidden in the size of the quantity $\hat{m}(\varepsilon)$, see the definition (55). To illustrate this, we now discuss possible advantage of the algorithm ADAPT-MESH with respect to another algorithm $\phi$ with $\bar{\beta}=2$, based on the uniform mesh. Note that for any $\phi$, any $m$ and the uniform mesh we have

[TABLE]

where $\bar{N}(f)=\bar{\beta}\left((1/r!)\sup\limits_{t\in[a,b]}\|z^{(r+1)}(t)\|+\beta\right)(1+\varphi)$. Hence,

[TABLE]

It follows from (48) that we also have

[TABLE]

and

[TABLE]

where constants in the $\Theta$-notation only depend on $\varphi$.

We want to compare $m^{*}(\varepsilon)$ (which decides about the cost of ADAPT-MESH) with $\hat{m}^{{\rm equid}}(\varepsilon)$ (which decides about the cost of $\phi$ with the equidistant mesh). We have by (65) the following sequence of inequalities

[TABLE]

[TABLE]

where again constants in the $\Theta$-notation only depend on $\varphi$.

The second inequality in (68) allows us to understand when the adaption pays off. Just below the definition (57), we gave a comment on how to compute $\bar{m}(\varepsilon)$. The comment yields that $\bar{m}(\varepsilon)$ is the minimal number $m\in{\bf N}$ such that

[TABLE]

where

[TABLE]

To see this, consult similar reasoning leading to (37) and (38). In the case of $\bar{m}^{{\rm equid}}\left(\varepsilon\right)$ an analogous condition to (69) is given in the first equality in (67):

[TABLE]

Comparing (69) and (71), we see that the second term in (68) is much less than the third term if $\bar{S}(m)^{-1}$ is much smaller than $\bar{N}(f)^{1/(r+1)}$. Since

[TABLE]

and

[TABLE]

we can identify cases when the gain of adaption is significant. Adaption pays off for functions for which the size of $\|z^{(r+1)}(t)\|$ changes significantly in parts of the interval $[a,b]$. Of course, the second inequality in (68) can also turn into equality. For such functions $f$ there is no gain of adaption. Translating the above discussion to similar properties of the cost of the algorithm is straightforward.

8 **Numerical example **

We illustrate the performance of the mesh selection mechanism in ADAPT-MESH by an example (other test examples in ${\rm C}^{++}$ are in progress), see [11]. We consider a scalar test problem from [4] with a parameter $\delta>0$

[TABLE]

with the global solution given by $z(t)=\left(\frac{15}{8}t+\delta^{5/2}\right)^{2/5}+1$. The solution with the initial condition $z(x)=y\;\;$ ($x\geq 0$, $y>1$) is given by

[TABLE]

The right-hand side $f(t,y)=\frac{3}{4}(y-1)^{-3/2}$ is a $C^{\infty}$ function for $y>1$.

The problem (72) is a typical test problem whose computational difficulty can be controlled by $\delta$; it grows as $\delta$ tends to zero. We use the algorithm ADAPT-MESH with $r=1$, which corresponds to the Euler method equipped with the mesh selection algorithm, and with $r=2$. For the global solution $z$, we have that $|z^{(r+1)}(t)|$ for $t$ close to [math] behaves like $1/\delta^{4}$ for $r=1$ and $1/\delta^{6.5}$ for $r=2$. For $t$ away from [math], $|z^{(r+1)}(t)|$ is essentially a constant. That is, for small $\delta$ we should observe a significant advantage of adaptive mesh points over the equidistant points.

The computer precision is $macheps=10^{-15}$. Obviously, since $macheps$ is fixed and computing time is limited, we cannot verify the asymptotic behavior of the algorithm as $\varepsilon\to 0$; we are only able to see results for some number of values of $\varepsilon$.

Let us now briefly discuss a practical choice of $h(\varepsilon)$ in step 2 of the algorithm. In fixed precision computation, the crucial point is accuracy of computing the divided difference in (27) of Lemma 2. Due to round off errors in computing both $f(\bar{t}_{k},\bar{l}_{i,r+1}(\bar{t}_{k}))$ and the divided difference, the bound (27) changes to $\tilde{C}(\bar{h}_{i}+macheps/\bar{h}_{i}^{r})$, for some $\tilde{C}$ dependent on $f$. The minimum of the function of $\bar{h}_{i}$ is achieved for $\bar{h}_{i}=(r\cdot macheps)^{1/(r+1)}$. Thus, in step 2, neglecting the coefficient dependent on $r$, we fix $h(\varepsilon)$ independently of $\varepsilon$ to be $h(\varepsilon)=10^{-15/(r+1)}$.

In step 5, we set

[TABLE]

(Note that for $r=1$ we have $l_{i,r+1}=\bar{l}_{i,r+1}$.) The following table shows results computed by ADAPT-MESH for a number of values of $\delta$ and $\varepsilon$. We denote

[TABLE]

where $z_{i}(x_{i}^{*})=y_{i}^{*}$, and $m^{*}$ is the number of intervals computed in ADAPT-MESH. The value EQUIDIST is the maximal local error of the respective method for $r=1$ or $r=2$ applied on the equidistant mesh $x_{i}=i/m^{*}$, $i=0,1,\ldots,m^{*}$, with the same number of subintervals equal to $m^{*}$. In the successive columns we show the values $m^{*}$ , MAXERR/$\varepsilon$ and EQUIDIST/$\varepsilon$.

[TABLE]

According to the theory, for sufficiently small $\varepsilon$ the values in the 4th column (for $r=1$) and 7th column (for $r=2$) should be at most 1. This is the case for $\delta=10^{-1},10^{-2}$ for $r=1,2$, and $\delta=10^{-3}$ for $r=1$. For small values of $\delta$, the round off errors do not allow us to observe the asymptotic behavior of the algorithm, since the value of $\varepsilon$ is too large. Comparison of columns 4 and 5 for $r=1$ and 7 and 8 for $r=2$ shows the gain of the adaptive mesh selection algorithm applied in ADAPT-MESH over the equidistant points. In the test we have computed results for the equidistant mesh with the same number of points. We may wish to compare the behavior of adaption with nonadaption using the same number of evaluations of $f$. For $r=1$, the adaptive method uses 2 function evaluations, while the nonadaptive one only one value. Hence, in this case the value in the 5th column should be divided by 4. For $r=2$, the respective numbers are 10 and 4 evaluations, that is the result in the 8th column should be divided roughly by 16. This does not change the picture – in both cases, for small $\delta$ the tests show a very significant advantage of the adaption over nonadaption.

We shortly comment on comparison between the algorithm defined in [4] for scalar autonomous problems and the current algorithm designed for systems of IVPs, for the test problem (72). As it can be expected, the algorithm from [4] allows us to better treat small values of $\delta$. This follows from the fact that, roughly speaking, the step size control in [4] was based on two-sided estimates of local errors. Specific properties of scalar autonomous problems were used in [4]; they cannot be extended to systems of initial value problems. In order to handle systems of IVPs, the present algorithm uses upper local error bounds, see (28) and (29).

9 Conclusions

We have proposed a mesh selection algorithm for systems of IVPs that (asymptotically) guarantees a given level of the local error. The algorithm only requires evaluations of the right-hand side $f$. Rigorous analysis of the cost has been given, including comparison with the best choice of the mesh points, as well as with the uniform mesh.