A new stepsize change technique for Adams methods

Consider the numerical solution of the initial value problem

y′(t)=f(t,y(t)),t∈[0,T]y(0)=y0$$\left\{\begin{array} {}y'(t) = f(t,y(t)), & \quad t \in [0,T] \\ y(0) = y_0 & \end{array}\right. $$(1)

where T > 0, and y₀ are given real numbers and f:[0,T]×R→R$\begin{array}{} \displaystyle f:[0,T]\times \mathbf{R} \to \mathbf{R} \end{array}$ is a continuously differentiable function in a neighborhood {(t, y) | t ∊ [0, T ], |y − y(t)| ≤ a} of the solution (t, y(t)), t ∊ [0, T] of (1). Although the following results can be also applied to systems of ODE’s, for simplicity we restrict our considerations to a single equation.

Assume that approximations yj,y′j,j=0,…,n−1,n≥k$\begin{array}{} \displaystyle {y_j},{y'_j},j = 0, \ldots ,n - 1,n \ge k \end{array}$ to the true solution and its derivative are known at the grid points t_j, j = 0,...,n−1 and we want to advance the numerical solution computing new approximations at t_n = t_n−1 +h_n with a k-step Adams-Moulton method (AMk). Then, denoting by P_n(t) the polynomial of degree ≤ k + 1 such that

Pn(tn)=yn,Pn(tn−1)=yn−1,P′n(tn−j)=f(tn−j,Pn−1(t¯n−j)),j=1,…,k,$$\begin{array}{} \displaystyle {P_n}({t_n}) = {y_n},\quad {P_n}({t_{n - 1}}) = {y_{n - 1}},\quad {P'_n}({t_{n - j}}) = f({t_{n - j}},{P_{n - 1}}({\bar t_{n - j}})),\quad j = 1, \ldots ,k, \end{array}$$(2)

where t¯n−j$\begin{array}{} \displaystyle \overline{t}_{n-j} \end{array}$ are back points that define the class of AMk method and the new approximation y_n is determined so that P_n(t) satisfies the differential equation at t_n, i.e.,

Pn'(tn)=f(tn,yn).$$\begin{array}{} \displaystyle P_n^{^\prime }({t_n}) = f({t_n},{y_n}). \end{array}$$(3)

By substituting the Lagrange form of the interpolatory polynomial (2) into (3), one gets the difference equations of the AMk method. However, for our purposes it turns out more convenient to deal with the polynomial form of the AMk method. To define the method in this way, note that polynomials P_n and P_n−1 associated to t_n and t_n−1 respectively satisfy

Pn(t)=Pn−1(t)+γqn(t),qn(t)=∫tn−1t(s−t¯n−1)…(s−t¯n−k)ds$$\begin{array}{} \displaystyle P_n(t) = P_{n-1}(t) + \gamma q_n(t), \quad q_n(t) = \int_{t_{n-1}}^t (s-\overline{t}_{n-1}) \ldots (s-\overline{t}_{n-k})\,\text{d}s \end{array}$$(4)

and γ is a real constant which may be determined from the two conditions

Pn(tn)=yn,P′n(tn)=f(tn,yn).$$\begin{array}{} \displaystyle P_n(t_n) = y_n, \quad\quad P'_n(t_n) = f(t_n , y_n). \end{array}$$(5)

The equations (4) imply that P_n and P_n−1 are such that their derivatives agree at the k back points t¯n−1,…,t¯n−k$\begin{array}{} \displaystyle {\bar t_{n - 1}}, \ldots ,{\bar t_{n - k}} \end{array}$ whereas y_n is chosen so that P_n(t) satisfies the ODE at the advanced point t_n.

On a uniform grid we take as back points t¯n−j=tn−jh$\begin{array}{} \displaystyle \overline{t}_{n-j}=t_n-jh \end{array}$ where h is the fixed stepsize of the grid. For non uniform grids there are many possibilities about the choice of back points. Firstly, taking as t¯n−j,j=1,…,k$\begin{array}{} \displaystyle {\bar t_{n - j}},j = 1, \ldots ,k \end{array}$ the past grid points we get the AMk method with the so called variable coefficient technique for stepsize changing. On the other hand, taking the k equally spaced back points t¯n−j=tn−jhn,j=1,…,k$\begin{array}{} \displaystyle {\bar t_{n - j}} = {t_n} - j{h_n},j = 1, \ldots ,k \end{array}$, we have the AMk method with the interpolation technique for stepsize changing. Other techniques have been proposed by Skeel [14], Gupta and Wallace [8] and Jackson and Sack-Davis [10].

The above techniques are the most usual in practical codes for non-stiff IVP’s. Thus, STEP written by Shampine and Gordon [12], DVDQ by Krogh [11] and EPISODE by Byrne and Hindmarsh [1] employ AM formulas with the variable coefficient technique while DIFSUB by Gear [5] and LSODE by Hindmarsh [9] use the interpolation technique.

The stability properties of these techniques have been studied by several authors (e.g. Gear and Tu [6], Crouzeix and Lisbona [4], Grigorieff [7], Skeel and Jackson [13] and Calvo, Lisbona and Montijano [2, 3]). From these papers it follows that the variable coefficient technique is more stable than the interpolation one. However, since the computational cost of the first is higher than the second one, for those problems with stability requirements not too strong, the interpolation technique can be advantageously used.

In this paper, a new technique for stepsize change in Adams methods is proposed. Its computational cost is similar to the interpolation technique but it has better stability properties and smaller leading error term. The paper is organized in the following manner: In section 2 the new technique is introduced in the frame of the polynomial formulation of AM methods. In section 3 the leading term of the local error is obtained. Finally, in section 4 the stability (0-stability) of this technique is studied and some comparisons with the interpolation technique are presented.

As it was remarked earlier, when AMk methods are written in polynomial form, the difference between interpolation and variable coefficient techniques is the choice of the back points t_n−j on which agree the derivatives Pn′=Pn−1′$\begin{array}{} \displaystyle P_n^{^\prime } = P_{n - 1}^{^\prime } \end{array}$ of two consecutive polynomials. Moreover, when the back points are equally spaced, the calculations required to advance one step are considerably simplified.

In view of these facts we will take the back points t¯n−j=tn−jh¯,j=1,…,k$\begin{array}{} \displaystyle \overline t_{n-j}=t_n - j \overline{h} , \quad j=1,\ldots, k \end{array}$ equally spaced with a stepsize h¯$\begin{array}{} \displaystyle \overline{h} \end{array}$ adjusted to improve the stability and leading error term over the corresponding ones of the interpolation technique.

Denoting by r_n = h_n/h_n−1 the stepsize ratio, we have considered h¯$\begin{array}{} \displaystyle \overline{h} \end{array}$ given by

h¯=hnϕ(rn)$$\begin{array}{} \displaystyle \overline{h} = h_n \phi(r_n) \end{array}$$

where ϕ(r) is a bounded function for 0 ≤ r_* < r < r^* with ϕ(1) = 1. In particular, the following functions have been studied:

ϕ(r)=αk+(1−αk)r−1,$$\begin{array}{} \displaystyle \phi(r) = \alpha_k + (1-\alpha_k) r^{-1}, \end{array}$$(T1)

ϕ(r)={αk+(1−αk)r−1,if r>11,if r≤1,$$\begin{array}{} \phi(r) = \begin{cases} \alpha_k + (1-\alpha_k) r^{-1}, & \text{if } r>1 \\ 1 , & \text{if } r \le 1, \end{cases} \end{array}$$(T2)

ϕ(r)=αk,if r>11,if r≤1.$$\begin{array}{} \phi(r) = \begin{cases} \alpha_k , & \text{if } r>1 \\ 1, & \text{if } r \le 1. \end{cases} \end{array}$$(T3)

where α_k is, for the AMk method, a constant that will be determined by imposing that the stability of the technique and the leading error term behave in an optimal way in a sense that will be made precise later.

Note that if α_k < 1, then ϕ(r) < 1 for all r for the T3 technique, and the same holds for the T1 and T2 techniques if α_k ∊ (0, 1), for r > 1, i.e., when the current stepsize is decreased, and therefore h¯<hn$\begin{array}{} \displaystyle \overline{h} <\ h_n \end{array}$. This means that our interpolation points t¯n−j$\begin{array}{} \displaystyle \bar t_{n-j} \end{array}$ are closer to t_n than those used in the interpolation technique. The choice ϕ(r) = 1 for r ≤ 1 in T2 and T3 is mainly due to the fact the stability in the interpolation technique is satisfactory for r ≤ 1.

On the other hand, in the techniques T1 and T2, the interpolation stepsize h¯$\begin{array}{} \displaystyle \overline{h} \end{array}$ is taken as a convex combination of the two last stepsizes h_n−1 and h_n, i.e.

h¯=αkhn+(1−αk)hn−1.$$\begin{array}{} \displaystyle \overline{h}= \alpha_k h_n + (1-\alpha_k) h_{n-1}. \end{array}$$

Next, to obtain the propagation matrices we write the AMK methods in matrix form by means of Nordsieck vectors (see [15]). Let

qn(t)=∫t¯n−1t(s−t¯n−1)…(s−t¯n−k)ds.$$\begin{array}{} \displaystyle q_n(t) = \int_{\overline t_{n-1}}^t (s-\overline t_{n-1}) \ldots(s-\overline t_{n-k})\,\text{d}s. \end{array}$$(6)

the modifier polynomial of the new technique and N_n(t, h) and C_n(t, h) the (k + 2)-Nordsieck vectors at the point t, with stepsize h of P_n(t) and q_n(t) respectively, given by

Nn(t,h)=(Pn(t),hP′n(t),…,hk+1(k+1)!Pn(k+1)(t))T,Cn(t,h)=(qn(t),hq′n(t),…,hk+1(k+1)!qn(k+1)(t))T.$$\begin{array}{} \displaystyle \begin{array}{*{20}{c}} {{N_n}(t,h){\rm{ }} = {{({P_n}(t),h{{P'}_n}(t), \ldots ,\frac{{{h^{k + 1}}}}{{(k + 1)!}}P_n^{(k + 1)}(t))}^T},} \\ {{C_n}(t,h){\rm{ }} = {{({q_n}(t),h{{q'}_n}(t), \ldots ,\frac{{{h^{k + 1}}}}{{(k + 1)!}}q_n^{(k + 1)}(t))}^T}.} \\ \end{array} \end{array}$$(7)

With these notations, the first equation of (4) can be written equivalently in the form

Nn(tn,hn)=Nn−1(tn,hn)+γCn(tn,hn).$$\begin{array}{} \displaystyle {N_{n}(t_{n},h_{n}) = N_{n-1}(t_{n},h_{n}) + \gamma\, C_n(t_{n},h_{n})}. \end{array}$$(8)

Now, using the Taylor formula to expand each component of N_n−1(t_n,h_n) at the point t_n−1 and substituting h_n by r_nh_n−1 we have

Nn−1(tn,hn)=PD(rn)Nn−1(tn−1,hn−1),$$\begin{array}{} \displaystyle {N_{n - 1}}\left( {{t_n},{h_n}} \right) = PD\left( {{r_n}} \right){N_{n - 1}}\left( {{t_{n - 1}},{h_{n - 1}}} \right), \end{array}$$(9)

where D(r) = diag(1,...,r^k+1) and P=(pij)i,j=0k+1∈R(k+2)×(k+2)$\begin{array}{} \displaystyle P=\left(p_{ij}\right)_{i,j=0}^{k+1}\in \mathbf{R}^{(k+2)\times(k+2)} \end{array}$ is the Pascal matrix whose entries are

pij=ji,if k+1≥j≥i≥00,otherwise.$$\begin{array}{} p_{ij}= \begin{cases} {j\choose i}, & \text{if } k+1\ge j \ge i \ge0 \\ \strut \ 0, & \text{otherwise.} \end{cases} \end{array}$$

On the other hand, putting in (6)s=tn+xh¯$\begin{array}{} \displaystyle s=t_n + x \overline{h} \end{array}$, the Nordsieck vector of q_n(t) at t_n becomes

Cn(tn,hn)=h¯k+1D(r¯n)c,$$\begin{array}{} \displaystyle C_n(t_n, h_n) = \overline{h}^{k+1} D(\bar r_n)\, c, \end{array}$$(10)

where c = (c₀, c₁,...,c_k+1)^T ∊ R^k+2 is a constant vector whose components are the coefficients of the polynomial

Λ(x)=∫−1x(s+1)…(s+k)ds=∑j=0k+1cjxj,$$\begin{array}{} \displaystyle \Lambda(x) = \int_{-1}^x (s+1) \ldots (s+k)\,\text{d}s = \sum_{j=0}^{k+1}c_j x^j, \end{array}$$(11)

and r¯n=hn/h¯=1/ϕ(r)$\begin{array}{} \displaystyle \bar r_n = h_n/\overline{h}=1/\phi(r) \end{array}$.

Substituting (9) and (10) into (8) and putting N_j = N_j(t_j, h_j) we have

Nn=PD(rn)Nn−1+γh¯k+1D(r¯n)c.$$\begin{array}{} \displaystyle N_n = P D(r_n) N_{n-1} + \gamma \overline{h}^{k+1} D(\bar r_n) c. \end{array}$$(12)

Finally, the value of constant γ is determined from the second component of (12) and the equation (5). Thus, we arrive at the matrix equation of the AMk method with the new technique of stepsize change

Nn=Ω(rn)Nn−1+hnc1r¯nf(tn,yn)D(r¯n)c$$\begin{array}{} \displaystyle N_n = \Omega(r_n) N_{n-1} + {h_n\over c_1 \bar r_n}f(t_n, y_n)D(\bar r_n) c \end{array}$$(13)

where the propagation matrix Ω(r) is given by

Ω(rn)=(I−ω(r¯n)−1D(r¯n)ce2T)PD(rn)ω(r¯n)=e2TD(r¯n)c=k!r¯n,e2=(0,1,0,…)T∈Rk+2,r¯n=1/ϕ(rn).$$\begin{array}{} \displaystyle \begin{array}{*{20}{c}} {\Omega ({r_n}) = (I - \omega {{({{\bar r}_n})}^{ - 1}}D({{\bar r}_n})ce_2^T)PD\left( {{r_n}} \right)} \hfill \\ {\omega ({{\bar r}_n}) = e_2^TD({{\bar r}_n})c = k!{\kern 1pt} {{\bar r}_n},\quad {e_2} = {{(0,1,0, \ldots )}^T} \in {{\bf{R}}^{k + 2}},} \hfill \\ {{{\bar r}_n} = 1/\phi ({r_n}).} \hfill \\ \end{array} \end{array}$$(14)

Although the matrix form (13)-(14) of the AMk method will be useful to study the stability of the methods, for practical purposes it is more convenient to write the equations in the equivalent form

Nn,0=PD(rn)Nn−1=(yn,0,hnyn,0′,…,hnk+1(k+1)!yn,0(k+1))Tyn=yn,0+1ℓ1r¯n(hf(tn,yn)−hyn,0′)Nn=Nn,0+(yn−yn,0)D(r¯n)ℓ$$\begin{array}{} \displaystyle \begin{array}{*{20}{c}} {{N_{n,0}} = PD({r_n}){N_{n - 1}} = {{({y_{n,0}},{h_n}y_{n,0}^\prime , \ldots ,\frac{{h_n^{k + 1}}}{{(k + 1)!}}y_{n,0}^{(k + 1)})}^T}} \hfill \\ {{y_n} = {y_{n,0}} + \frac{1}{{{\ell _1}{{\overline r }_n}}}(hf({t_n},{y_n}) - hy_{n,0}^\prime )} \hfill \\ {{N_n} = {N_{n,0}} + ({y_n} - {y_{n,0}})D({{\overline r }_n})\ell } \hfill \\ \end{array} \end{array}$$(15)

where ℓ is the (k + 2)-vector defined by ℓ=1c0c=(1,ℓ1,…,ℓk+1)T$\begin{array}{} \displaystyle \ell={\displaystyle 1\over c_0} c = (1, \ell_1, \ldots, \ell_{k+1})^T \end{array}$. The first equation of (15) is the prediction stage in which the Nordsieck vector at t_n is predicted from the corresponding one at t_n−1. The second equation of (15) is a nonlinear algebraic equation which must solved by some iterative method (e.g. fixed point) to get y_n. In the last equation (15) the predicted N_n,0 is corrected to the final value N_n.

To define the truncation error EL(t_n) at the grid point t_n, we assume that the true solution y(t) is sufficiently differentiable and then

EL(tn)=y(tn)−yn*$$\begin{array}{} \displaystyle EL(t_n) = y(t_n) - y_n^* \end{array}$$

where yn*$\begin{array}{} \displaystyle y_n^* \end{array}$ is the solution computed by using (15) with an exact Nordsieck vector at t_n−1, i.e.

Nn−1*=(y(tn−1),hn−1y′(tn−1),…hn−1k+1(k+1)!y(k+1)(tn−1))T.$$\begin{array}{} \displaystyle N_{n-1}^* = \big( y(t_{n-1}), h_{n-1} y'(t_{n-1}), \ldots {h_{n-1}^{k+1}\over (k+1)!} y^{(k+1)}(t_{n-1}) \big)^T. \end{array}$$

Next, let us calculate the leading term of EL(t_n) when it is expanded in powers of h_n. We write it in the form

EL(tn)=(y(tn)−yn,0)+(yn,0−yn*).$$\begin{array}{} \displaystyle EL(t_n) = (y(t_n) - y_{n,0}) + ( y_{n,0} - y_n^*). \end{array}$$(16)

Taking into account the first equation of (15), the first term in the right hand side of (16) is

y(tn)−yn,0=hnk+2(k+2)!y(k+2)(tn−1)+O(hnk+3).$$\begin{array}{} \displaystyle y(t_n)-y_{n,0} = {h_n^{k+2}\over (k+2)!} \, y^{(k+2)}(t_{n-1}) + {\mathscr O}(h_n^{k+3}). \end{array}$$

Similarly, from the second equation of (15) we get

yn*−yn,0=1ℓ1r¯nhnk+2(k+1)!y(k+2)(tn−1)+O(hnk+3),$$\begin{array}{} \displaystyle y_n^*- y_{n,0} ={1\over \ell_1 \bar r_n} {h_n^{k+2}\over (k+1)!} \, y^{(k+2)}(t_{n-1}) + {\mathscr O}(h_n^{k+3}), \end{array}$$

and therefore,

EL(tn)=Ck+2hnk+2y(k+2)(tn−1)+O(hnk+3),$$\begin{array}{} \displaystyle EL(t_n) = C_{k+2} h_n^{k+2}y^{(k+2)}(t_{n-1}) + {\mathscr O}(h_n^{k+3}), \end{array}$$

with

Ck+2=1(k+2)!(1−k+2r¯nℓ1).$$\begin{array}{} \displaystyle {C_{k + 2}} = \frac{1}{{(k + 2)!}}\left( {1 - \frac{{k + 2}}{{{{\bar r}_n}{\ell _1}}}} \right). \end{array}$$

The values of (k + 2)/ℓ₁, for 2 ≤ k ≤ 8 are given in Table 1.

Table 1

Since all these values are greater than unity, it follows that

|Ck+2(1)|>|Ck+2(r¯)|$$\begin{array}{} \displaystyle |C_{k+2}(1)| > |C_{k+2}(\bar r)| \end{array}$$

whenever 1<r¯<(k+2)/ℓ12−(k+2)/ℓ1$\begin{array}{} \displaystyle \displaystyle 1<\bar r<{(k+2)/\ell_1 \over 2- (k+2)/\ell_1} \end{array}$, i.e., when 1>ϕ(r)>2ℓ1−k−2k+2$\begin{array}{} \displaystyle \displaystyle 1>\phi(r)>{2\ell_1 - k-2 \over k+2} \end{array}$. This means that for both T2 and T3 techniques with α_k ∊ (0, 1) the coefficient of EL(t_n) is smaller than the corresponding to the IT. On the other hand, for the T1 technique, this coefficient is smaller only if r_n > 1.

In this section we study the 0-stability of AMk methods with the techniques T1 and T2 for the stepsize changing introduced in section 2. Taking into account that f (t, y) satisfies a uniform Lipschitz condition in a neighborhood of the true solution and assuming that the stepsize ratios are bounded, it may be proved (e.g. [4,7]) that the stability of (13) is independent of the non-linear term and therefore the method (13) is stable if and only if there exist numbers K, r*, h* > 0 such that for any grid {t_n} ⊂ [0, T ] with h_n = t_n − t_n−1 ≤ h, and r_n ≤ r*, the inequality

‖Ωk(rn)Ωk(rn−1)…Ωk(rm)‖<K<∞$$\begin{array}{} \displaystyle { \Vert \Omega_k(r_n)\Omega_k(r_{n-1})\ldots \Omega_k(r_m)\Vert < K} <\ \infty \end{array}$$(17)

holds for all m ≤ n

Since the propagation matrix Ω_k(r) depends only on the stepsize ratio and on function ϕ(r), following the ideas of [2] we introduce

Next, we turn our attention to the choice of the α-parameter in the new techniques. We consider first the two steps AM method. In this case, after some elementary calculations it is found that Ω¯2(r)$\begin{array}{} \displaystyle \overline{\Omega}_2(r) \end{array}$ is a 2 × 2 matrix given by

Ω¯2(r)=r2(1−32r¯(3−94r¯)r−13r¯2(1−12r¯2)r)$$\begin{array}{} \displaystyle {\bar \Omega _2}(r) = {r^2}\left( {\begin{array}{*{20}{c}} {1 - \frac{3}{2}\bar r{\rm{ }}} & {(3 - \frac{9}{4}\bar r)r} \\ { - \frac{1}{3}{{\bar r}^2}} & {(1 - \frac{1}{2}{{\bar r}^2})r} \\ \end{array}} \right) \end{array}$$

where r¯=1/ϕ(r)$\begin{array}{} \displaystyle \bar r=1/\phi(r) \end{array}$. By using Schur’s criteria, it follows that the spectral radius ρ(Ω¯2(r))<1$\begin{array}{} \displaystyle \rho(\overline{\Omega}_2(r))< 1 \end{array}$ if and only if the following inequalities hold

r5|(r¯−1)(r¯−2)|<2r2|(2−3r¯)+r(2−r¯2)|<|2+r5(r¯−1)(r¯−2)|.$$\begin{array}{} \displaystyle \begin{array}{*{20}{c}} {{r^5}|(\bar r - 1)(\bar r - 2)| < 2} \hfill \\ {{r^2}|(2 - 3\bar r) + r(2 - {{\bar r}^2})| <\ |2 + {r^5}(\bar r - 1)(\bar r - 2)|.} \hfill \\ \end{array} \end{array}$$(18)

Thus, starting with the technique T1, for each value α = α₂ we may found the set of all r ≥ 0 such that ρ(Ω¯2(r))<1$\begin{array}{} \displaystyle \rho(\overline{\Omega}_2(r))< 1 \end{array}$ which is given by (18). This set has the form [0, φ₁(α)). Figure 1 shows the function φ₁(α), and a remarkable fact is that it attains a maximum for α=α2*=0.7677$\begin{array}{} \displaystyle \alpha=\alpha_2^* = 0.7677 \end{array}$. In view of this, a reasonable choice for the AM2 method with the T1 technique would be α₂ = 0.7677. Note that in the interpolation technique α₂ = 1 and then the spectral radius ρ(Ω¯2(r))<1$\begin{array}{} \displaystyle \rho(\overline{\Omega}_2(r))<1 \end{array}$ only for r ∊ [0, 1.695) while for the new technique, the same condition is satisfied for r ∊ [0, 1.803). A similar study can be carried out for the remaining techniques as well as for k > 3. In Table 2 the optimal values of the α-parameters for the T1 and T3 techniques and k ≤ 7 are given. Moreover, Table 2 includes the greatest value of r_k such that ρ(Ω¯k(r))<1$\begin{array}{} \displaystyle \rho(\overline{\Omega}_k(r))<1 \end{array}$ for r ∊ [0, r_k) for the T1, T3 and IT techniques respectively. In all cases, it is clear that the stability of the new techniques is better than that of the IT.

Fig. 1

Table 2

Although the sets Ik={r|ρ(Ω¯k(r))<1}$\begin{array}{} \displaystyle I_k= \{ r\; \vert\; \rho(\overline{\Omega}_k(r))<1 \} \end{array}$, see Figures 2-7, are good guides to select the optimal values of the parameters, it is clear that I_k is not in general a stability set (see e.g. Th. 4.1 [3]). One approach to get stability sets consist in finding a suitable matrix H such that the set M_k defined by

Fig. 2

Fig. 3

Fig. 4

Fig. 5

Fig. 6

Fig. 7

Mk={r∈R+;||H−1Ω¯k(r)H||1<1}⊂Ik$$\begin{array}{} \displaystyle M_k = \{ r\in R^+ \ ; \ \vert \vert H^{-1} \overline{\Omega}_k(r) H \vert \vert_1 <\ 1 \} \subset I_k \end{array}$$(19)

be as large as possible. Moreover, it is convenient for practical applications for M_k to have the form [0, d_k) with d_k > 1. It has been proved in [3], that all closed set of M_k is a stability set for the AMk method.

Finally, we briefly outline the choice of H for the new techniques and the corresponding stability sets. The approach is purely numerical and it follows essentially the same ideas used in [3]. For a value r ∊ I_k such that Ω¯k(r)$\begin{array}{} \displaystyle \overline{\Omega}_k(r) \end{array}$ possesses a complete set of eigenvectors, we take these vectors as columns of H(r). Then, we get the corresponding set M_k(H(r)). By a numerical search we obtain an optimal r_k ∊ I_k. In Table 3 the optimal values of r_k for k ≤ 6 are given. It must be pointed that for the T1 technique M_k = I_k, k = 2,...,5 and for the T3 technique M_k = I_k only for k = 2 and k = 4. Therefore, optimal norms have obtained in these cases.

Table 3