2010 Mathematics Subject Classification : $65\mathrm{L}04,65\mathrm{L}05,65\mathrm{H}05$$65\mathrm{L}04,65\mathrm{L}05,65\mathrm{H}05$65L04,65L05,65H0565 \mathrm{~L} 04,65 \mathrm{~L} 05,65 \mathrm{H} 05$65\mathrm{L}04,65\mathrm{L}05,65\mathrm{H}05$
Keywords: initial value problems, stiff equations, Steffensen method, Newton method, convergence order

The mathematical modeling of many problems in physics, engineering, chemistry, biology, etc. gives rise to ordinary differential equations or systems of ordinary differential equations.

It is known that a high-order initial value problem (IVP) for differential equations or systems of equations can be rewritten as a first-order IVP system (see e.g. [4, [5]) so that the standard IVP can be written in the form:

where: $I\subseteq \mathbb{R},y_0\in {\mathbb{R}}^m,f:I\times {\mathbb{R}}^m\to {\mathbb{R}}^m$$I\subseteq \mathbb{R},y_0\in {\mathbb{R}}^m,f:I\times {\mathbb{R}}^m\to {\mathbb{R}}^m$I subeR,y_(0)inR^(m),f:I xxR^(m)rarrR^(m)I \subseteq \mathbb{R}, y_{0} \in \mathbb{R}^{m}, f: I \times \mathbb{R}^{m} \rightarrow \mathbb{R}^{m}$I\subseteq \mathbb{R},y_0\in {\mathbb{R}}^m,f:I\times {\mathbb{R}}^m\to {\mathbb{R}}^m$ and $a\in I$$a\in I$a in Ia \in I$a\in I$.
A solution is sought on the interval $[a,b]\subset I$$[a,b]\subset I$[a,b]sub I[a, b] \subset I$[a,b]\subset I$, where $a,b$$a,b$a,ba, b$a,b$ are finite. In this paper we consider only the scalar case, i.e., $m=1$$m=1$m=1m=1$m=1$.

In practice, the number of cases where an exact solution can be found by analytical means is very limited, so that one uses numerical methods for the approximation of the solution.

Integrating (1), for $m=1$$m=1$m=1m=1$m=1$, using an implicit linear multistep method with step-size $h$$h$hh$h$, leads to the solving at each step of an equation of the form:

Here $A$$A$AA$A$ is a constant determined by the numerical method and $\psi$$\psi$psi\psi$\psi$ is a known value.
This equation can be solved by the fixed point iteration

which converges to the unique solution of (2) provided that:

where $L$$L$LL$L$ is the Lipschitz constant of $\varphi$$\varphi$phi\phi$\varphi$ with respect to the second variable.
Condition (4) becomes too restrictive for stiff problems. Thus, if we use an explicit method to solve a stiff equation, we have to use an excessively small step-size to avoid instability; if we use an implicit method with an absolute stability region large enough to impose no stability restriction, we can choose a step-size as large as we want, but we will not be able to solve the implicit equation (2) by the iteration (33) unless the step-size is excessively small.

In order to overcome this difficulty one uses the Newton iteration instead of the fixed point iteration. Newton iteration applied to the equation:

where $F:[c,d]\to \mathbb{R},c,d\in \mathbb{R},c<d$$F:[c,d]\to \mathbb{R},c,d\in \mathbb{R},c<d$F:[c,d]rarrR,c,d inR,c < dF:[c, d] \rightarrow \mathbb{R}, c, d \in \mathbb{R}, c<d$F:[c,d]\to \mathbb{R},c,d\in \mathbb{R},c<d$, has the form:

When applied to the equation (2), where $F(y)=y-hA\varphi (x,y)-\psi$$F(y)=y-hA\varphi (x,y)-\psi$F(y)=y-hA phi(x,y)-psiF(y)=y-h A \phi(x, y)-\psi$F(y)=y-hA\varphi (x,y)-\psi$, we get:

i.e.,

One step of Newton iteration requests considerably more computing time than one step of fixed point iteration. Each step of the latter costs one function evaluation, whereas each step of the former calls for the updating of the derivative.

In this paper we approximate the solution of equation (5) using the Steffensen type method:

where $g:[c,d]\to [c,d]$$g:[c,d]\to [c,d]$g:[c,d]rarr[c,d]g:[c, d] \rightarrow[c, d]$g:[c,d]\to [c,d]$ is an auxiliary function such that the equation:

is equivalent to (5), and $[u,v;F]$$[u,v;F]$[u,v;F][u, v ; F]$[u,v;F]$ represents the first order divided difference of $F$$F$FF$F$ at the points $u,v\in [c,d]$$u,v\in [c,d]$u,v in[c,d]u, v \in[c, d]$u,v\in [c,d]$. This method does not require the calculation of the derivative of the function $F$$F$FF$F$.

Let $y^{\ast }\in (c,d)$$y^{\ast }\in (c,d)$y^(**)in(c,d)y^{*} \in(c, d)$y^{\ast }\in (c,d)$ be the root of equation (5). If the elements of the sequence ${\left(y^{(\nu )}\right)}_{\nu \ge 0}$${\left(y^{(\nu )}\right)}_{\nu \ge 0}$(y^((nu)))_(nu >= 0)\left(y^{(\nu)}\right)_{\nu \geq 0}${\left(y^{(\nu )}\right)}_{\nu \ge 0}$ belong to the interval $[c,d]$$[c,d]$[c,d][c, d]$[c,d]$ then from Newton identity and (9) we obtain:

where $[u,v,w;F]$$[u,v,w;F]$[u,v,w;F][u, v, w ; F]$[u,v,w;F]$ represents the second order divided difference of $F$$F$FF$F$ at the points $u,v,w\in [c,d]$$u,v,w\in [c,d]$u,v,w in[c,d]u, v, w \in[c, d]$u,v,w\in [c,d]$. If $g$$g$gg$g$ is Lipschitz on $[c,d]$$[c,d]$[c,d][c, d]$[c,d]$ with constant $L$$L$LL$L$ and if we assume that there exist the real numbers $M,m>0$$M,m>0$M,m > 0M, m>0$M,m>0$ such that:

for all $u,v,w\in [c,d]$$u,v,w\in [c,d]$u,v,w in[c,d]u, v, w \in[c, d]$u,v,w\in [c,d]$, then:

which shows that the $q$$q$qq$q$-convergence order for the method (9) is 2 , i.e., the same as for the Newton method.
In 77 are given conditions for the convergence of the sequences generated by relation (9), and the function $g$$g$gg$g$ is defined such that the sequences ${\left(y^{(\nu )}\right)}_{\nu \ge 0}$${\left(y^{(\nu )}\right)}_{\nu \ge 0}$(y^((nu)))_(nu >= 0)\left(y^{(\nu)}\right)_{\nu \geq 0}${\left(y^{(\nu )}\right)}_{\nu \ge 0}$ and ${\left(g\left(y^{(\nu )}\right)\right)}_{\nu \ge 0}$${\left(g\left(y^{(\nu )}\right)\right)}_{\nu \ge 0}$(g(y^((nu))))_(nu >= 0)\left(g\left(y^{(\nu)}\right)\right)_{\nu \geq 0}${\left(g\left(y^{(\nu )}\right)\right)}_{\nu \ge 0}$ approximate bilaterally the exact solution $y^{\ast }$$y^{\ast }$y^(**)y^{*}$y^{\ast }$. Thus, we have an a posteriori error control.

For the functions $F$$F$FF$F$ and $g$$g$gg$g$ we suppose the following hypothesis:
( $\alpha$$\alpha$alpha\alpha$\alpha$ ) the equations (5) and (10) are equivalent;
$(\beta )$$(\beta )$(beta)(\beta)$(\beta )$ the function $g$$g$gg$g$ is continuous and decreasing on $[c,d]$$[c,d]$[c,d][c, d]$[c,d]$;
( $\gamma$$\gamma$gamma\gamma$\gamma$ ) the equation (5) has a unique solution $y^{\ast }\in (c,d)$$y^{\ast }\in (c,d)$y^(**)in(c,d)y^{*} \in(c, d)$y^{\ast }\in (c,d)$.
The following theorem holds (see [7]):
Theorem 1 If the functions $F$$F$FF$F$ and $g$$g$gg$g$ satisfy the conditions $(\alpha )-(\gamma )$$(\alpha )-(\gamma )$(alpha)-(gamma)(\alpha)-(\gamma)$(\alpha )-(\gamma )$ and moreover the following conditions hold:
(i) $F$$F$FF$F$ is increasing and convex on $[c,d]$$[c,d]$[c,d][c, d]$[c,d]$;
(ii) $F\left(y_0\right)<0$$F\left(y_0\right)<0$F(y_(0)) < 0F\left(y_{0}\right)<0$F\left(y_0\right)<0$;
(iii) $g\left(y_0\right)\le d$$g\left(y_0\right)\le d$g(y_(0)) <= dg\left(y_{0}\right) \leq d$g\left(y_0\right)\le d$,
then the elements of the sequences ${\left(y^{(\nu )}\right)}_{\nu \ge 0}$${\left(y^{(\nu )}\right)}_{\nu \ge 0}$(y^((nu)))_(nu >= 0)\left(y^{(\nu)}\right)_{\nu \geq 0}${\left(y^{(\nu )}\right)}_{\nu \ge 0}$ and ${\left(g\left(y^{(\nu )}\right)\right)}_{\nu \ge 0}$${\left(g\left(y^{(\nu )}\right)\right)}_{\nu \ge 0}$(g(y^((nu))))_(nu >= 0)\left(g\left(y^{(\nu)}\right)\right)_{\nu \geq 0}${\left(g\left(y^{(\nu )}\right)\right)}_{\nu \ge 0}$ belong to the interval $[c,d]$$[c,d]$[c,d][c, d]$[c,d]$ and the following properties hold:
( $j$$j$jj$j$ ) the sequence ${\left(y^{(\nu )}\right)}_{\nu \ge 0}$${\left(y^{(\nu )}\right)}_{\nu \ge 0}$(y^((nu)))_(nu >= 0)\left(y^{(\nu)}\right)_{\nu \geq 0}${\left(y^{(\nu )}\right)}_{\nu \ge 0}$ is increasing and convergent;
(jj) the sequence ${\left(g\left(y^{(\nu )}\right)\right)}_{\nu \ge 0}$${\left(g\left(y^{(\nu )}\right)\right)}_{\nu \ge 0}$(g(y^((nu))))_(nu >= 0)\left(g\left(y^{(\nu)}\right)\right)_{\nu \geq 0}${\left(g\left(y^{(\nu )}\right)\right)}_{\nu \ge 0}$ is decreasing and convergent;
(jjj) $y^{(\nu )}\le y^{\ast }\le g\left(y^{(\nu )}\right),\nu =0,1,\dots$$y^{(\nu )}\le y^{\ast }\le g\left(y^{(\nu )}\right),\nu =0,1,\dots$y^((nu)) <= y^(**) <= g(y^((nu))),quad nu=0,1,dotsy^{(\nu)} \leq y^{*} \leq g\left(y^{(\nu)}\right), \quad \nu=0,1, \ldots$y^{(\nu )}\le y^{\ast }\le g\left(y^{(\nu )}\right),\nu =0,1,\dots$.
(jv) $\underset{\nu \to \mathrm{\infty }}{lim}{y}^{(\nu )}=\underset{\nu \to \mathrm{\infty }}{lim}g\left(y^{(\nu )}\right)=y^{\ast }$$\underset{\nu \to \mathrm{\infty }}{lim} y^{(\nu )}=\underset{\nu \to \mathrm{\infty }}{lim} g\left(y^{(\nu )}\right)=y^{\ast }$lim_(nu rarr oo)y^((nu))=lim_(nu rarr oo)g(y^((nu)))=y^(**)\lim _{\nu \rightarrow \infty} y^{(\nu)}=\lim _{\nu \rightarrow \infty} g\left(y^{(\nu)}\right)=y^{*}$\underset{\nu \to \mathrm{\infty }}{lim}{y}^{(\nu )}=\underset{\nu \to \mathrm{\infty }}{lim}g\left(y^{(\nu )}\right)=y^{\ast }$;
(vj) $|y^{\ast }-y^{(\nu )}|\le |g\left(y^{(\nu )}\right)-y^{(\nu )}|,\nu =0,1,\dots$$\left|y^{\ast }-y^{(\nu )}\right|\le \left|g\left(y^{(\nu )}\right)-y^{(\nu )}\right|,\nu =0,1,\dots$|y^(**)-y^((nu))| <= |g(y^((nu)))-y^((nu))|,quad nu=0,1,dotsquad\left|y^{*}-y^{(\nu)}\right| \leq\left|g\left(y^{(\nu)}\right)-y^{(\nu)}\right|, \quad \nu=0,1, \ldots \quad$|y^{\ast }-y^{(\nu )}|\le |g\left(y^{(\nu )}\right)-y^{(\nu )}|,\nu =0,1,\dots$.
In the above theorem the auxiliary function $g$$g$gg$g$ can be taken as: $g(y)=y-\frac{F(y)}{F^{\prime }(c)}$$g(y)=y-\frac{F(y)}{F^{\prime }(c)}$g(y)=y-(F(y))/(F^(')(c))g(y)=y-\frac{F(y)}{F^{\prime}(c)}$g(y)=y-\frac{F(y)}{F^{\prime }(c)}$. Similar results have been obtained in [7] if $F$$F$FF$F$ verifies the properties:
-F is increasing and concave; $g$$g$gg$g$ can be taken as $g(y)=y-\frac{F(y)}{F^{\prime }(d)}$$g(y)=y-\frac{F(y)}{F^{\prime }(d)}$g(y)=y-(F(y))/(F^(')(d))g(y)=y-\frac{F(y)}{F^{\prime}(d)}$g(y)=y-\frac{F(y)}{F^{\prime }(d)}$;
-F is decreasing and concave; $g$$g$gg$g$ can be taken as $g(y)=y-\frac{F(y)}{F^{\prime }(c)}$$g(y)=y-\frac{F(y)}{F^{\prime }(c)}$g(y)=y-(F(y))/(F^(')(c))g(y)=y-\frac{F(y)}{F^{\prime}(c)}$g(y)=y-\frac{F(y)}{F^{\prime }(c)}$;
-F is decreasing and convex; $g$$g$gg$g$ can be taken as $g(y)=y-\frac{F(y)}{F^{\prime }(d)}$$g(y)=y-\frac{F(y)}{F^{\prime }(d)}$g(y)=y-(F(y))/(F^(')(d))g(y)=y-\frac{F(y)}{F^{\prime}(d)}$g(y)=y-\frac{F(y)}{F^{\prime }(d)}$.
The interval $[a,b]$$[a,b]$[a,b][a, b]$[a,b]$ is partitioned by the point set $\left\{x_n\right\}$$\left\{x_n\right\}${x_(n)}\left\{x_{n}\right\}$\left\{x_n\right\}$ defined by $x_n=a+nh,n=0,1,\dots ,N,h=(b-a)/N,$$x_n=a+nh,n=0,1,\dots ,N,h=(b-a)/N,$x_(n)=a+nh,n=0,1,dots,N,h=(b-a)//N,quadx_{n}=a+n h, n=0,1, \ldots, N, h=(b-a) / N, \quad$x_n=a+nh,n=0,1,\dots ,N,h=(b-a)/N,$ and $y_n$$y_n$quady_(n)quad\quad y_{n} \quad$y_n$ denotes an approximation to the exact solution $y$$y$yy$y$ of (1) at $x_n$$x_n$x_(n)x_{n}$x_n$.

If we use an implicit linear multistep method then $y_n,n=1,\dots ,N$$y_n,n=1,\dots ,N$y_(n),n=1,dots,Ny_{n}, n=1, \ldots, N$y_n,n=1,\dots ,N$, are the solutions of the equation:

where ${\psi }_n={\psi }_n(a,h,y_{n-1},y_{n-2},\dots ,y_0)$${\psi }_n={\psi }_n\left(a,h,y_{n-1},y_{n-2},\dots ,y_0\right)$psi_(n)=psi_(n)(a,h,y_(n-1),y_(n-2),dots,y_(0))\psi_{n}=\psi_{n}\left(a, h, y_{n-1}, y_{n-2}, \ldots, y_{0}\right)${\psi }_n={\psi }_n(a,h,y_{n-1},y_{n-2},\dots ,y_0)$. We call this equation as approximant equation and we denote by $y_n^{\ast }\in (c,d),n=1,\dots ,N$$y_n^{\ast }\in (c,d),n=1,\dots ,N$y_(n)^(**)in(c,d),n=1,dots,Ny_{n}^{*} \in(c, d), n=1, \ldots, N$y_n^{\ast }\in (c,d),n=1,\dots ,N$, the exact solution.

For each $n=1,\dots ,N$$n=1,\dots ,N$n=1,dots,Nn=1, \ldots, N$n=1,\dots ,N$ let $F_n:[c,d]\to \mathbb{R}$$F_n:[c,d]\to \mathbb{R}$F_(n):[c,d]rarrRF_{n}:[c, d] \rightarrow \mathbb{R}$F_n:[c,d]\to \mathbb{R}$ be defined by

Then equation (11) can be rewritten in the form $F_n(y)=0$$F_n(y)=0$F_(n)(y)=0F_{n}(y)=0$F_n(y)=0$.
To approximate bilaterally the solution $y_n^{\ast },n=1,\dots ,N$$y_n^{\ast },n=1,\dots ,N$y_(n)^(**),n=1,dots,Ny_{n}^{*}, n=1, \ldots, N$y_n^{\ast },n=1,\dots ,N$, we generate the sequence ${\left(y_n^{(\nu )}\right)}_{\nu \ge 0}$${\left(y_n^{(\nu )}\right)}_{\nu \ge 0}$(y_(n)^((nu)))_(nu >= 0)\left(y_{n}^{(\nu)}\right)_{\nu \geq 0}${\left(y_n^{(\nu )}\right)}_{\nu \ge 0}$, by:

or, using (12),

From Theorem 1.1, if $F_n$$F_n$F_(n)F_{n}$F_n$ is increasing and convex, and the initial guess $y_n^{(0)}$$y_n^{(0)}$y_(n)^((0))y_{n}^{(0)}$y_n^{(0)}$ satisfy $F_n\left(y_n^{(0)}\right)<0$$F_n\left(y_n^{(0)}\right)<0$F_(n)(y_(n)^((0))) < 0F_{n}\left(y_{n}^{(0)}\right)<0$F_n\left(y_n^{(0)}\right)<0$, then the sequences ${\left(y_n^{(\nu )}\right)}_{\nu \ge 0},{\left(g\left(y_n^{(\nu )}\right)\right)}_{\nu \ge 0}$${\left(y_n^{(\nu )}\right)}_{\nu \ge 0},{\left(g\left(y_n^{(\nu )}\right)\right)}_{\nu \ge 0}$(y_(n)^((nu)))_(nu >= 0),(g(y_(n)^((nu))))_(nu >= 0)\left(y_{n}^{(\nu)}\right)_{\nu \geq 0},\left(g\left(y_{n}^{(\nu)}\right)\right)_{\nu \geq 0}${\left(y_n^{(\nu )}\right)}_{\nu \ge 0},{\left(g\left(y_n^{(\nu )}\right)\right)}_{\nu \ge 0}$ converge to $y_n^{\ast }$$y_n^{\ast }$y_(n)^(**)y_{n}^{*}$y_n^{\ast }$ and we have the inequalities:

The rest of the cases can be treated in a similar fashion.

We consider the trapezoidal rule to integrate the initial value problem (1), for $m=1$$m=1$m=1m=1$m=1$, and the Steffensen method described above to solve the approximant equation (11).

The trapezoidal rule is a 1 -step Adams-Moulton method (an implicit method), and for (1) is defined by:

It is known that the trapezoidal rule is an $A$$A$AA$A$-stable method and has order 2 (see [5]).
For any point $x_n,n=1,\dots ,N$$x_n,n=1,\dots ,N$x_(n),n=1,dots,Nx_{n}, n=1, \ldots, N$x_n,n=1,\dots ,N$, we have:

and in this case $F_n(y)=y-\frac{h}{2}f(x_n,y)-\frac{h}{2}f(x_{n-1},y_{n-1})-y_{n-1}$$F_n(y)=y-\frac{h}{2}f\left(x_n,y\right)-\frac{h}{2}f\left(x_{n-1},y_{n-1}\right)-y_{n-1}$F_(n)(y)=y-(h)/(2)f(x_(n),y)-(h)/(2)f(x_(n-1),y_(n-1))-y_(n-1)F_{n}(y)=y-\frac{h}{2} f\left(x_{n}, y\right)-\frac{h}{2} f\left(x_{n-1}, y_{n-1}\right)-y_{n-1}$F_n(y)=y-\frac{h}{2}f(x_n,y)-\frac{h}{2}f(x_{n-1},y_{n-1})-y_{n-1}$. Thus, in (11) we have $A=\frac{1}{2}$$A=\frac{1}{2}$A=(1)/(2)A=\frac{1}{2}$A=\frac{1}{2}$, $\varphi (x_n,y)=f(x_n,y)$$\varphi \left(x_n,y\right)=f\left(x_n,y\right)$phi(x_(n),y)=f(x_(n),y)\phi\left(x_{n}, y\right)=f\left(x_{n}, y\right)$\varphi (x_n,y)=f(x_n,y)$ and ${\psi }_n=\frac{h}{2}f(x_{n-1},y_{n-1})+y_{n-1},n=1,\dots ,N$${\psi }_n=\frac{h}{2}f\left(x_{n-1},y_{n-1}\right)+y_{n-1},n=1,\dots ,N$psi_(n)=(h)/(2)f(x_(n-1),y_(n-1))+y_(n-1),n=1,dots,N\psi_{n}=\frac{h}{2} f\left(x_{n-1}, y_{n-1}\right)+y_{n-1}, n=1, \ldots, N${\psi }_n=\frac{h}{2}f(x_{n-1},y_{n-1})+y_{n-1},n=1,\dots ,N$.

For simplicity we consider only the autonomous case, i.e. $f=f(y)$$f=f(y)$f=f(y)f=f(y)$f=f(y)$, and in this case equation (15) becomes:

and $F_n(y)=y-\frac{h}{2}f(y)-{\psi }_n,{\psi }_n=\frac{h}{2}f\left(y_{n-1}\right)+y_{n-1},n=1,\dots ,N$$F_n(y)=y-\frac{h}{2}f(y)-{\psi }_n,{\psi }_n=\frac{h}{2}f\left(y_{n-1}\right)+y_{n-1},n=1,\dots ,N$F_(n)(y)=y-(h)/(2)f(y)-psi_(n),psi_(n)=(h)/(2)f(y_(n-1))+y_(n-1),n=1,dots,NF_{n}(y)=y-\frac{h}{2} f(y)-\psi_{n}, \psi_{n}=\frac{h}{2} f\left(y_{n-1}\right)+y_{n-1}, n=1, \ldots, N$F_n(y)=y-\frac{h}{2}f(y)-{\psi }_n,{\psi }_n=\frac{h}{2}f\left(y_{n-1}\right)+y_{n-1},n=1,\dots ,N$.
Using the fact that

and