where f _n-i =f(t _n-i ,y _n-i ) and M>= n>= k.

If (a₀²+b₀²)(a²_k+b²_k)# 0, then k is the formula stepnumber.

The formula is an explicit one if b_k=0, else it is an implicit one. If a_k=0, b_k#0, then the formula is an extended one. At the step n, knowing the (k-1)^th previous approximate values, we get the k^th value by writing the formula in the form a_ky_n-hb_kf_n=constant . In the explicit case, when b_k=0, solving this equation means only a simple division of a constant expression with a_k. In the implicit case and a N-dimensional system function f we must, generally, solve a nonlinear system of N equations, so that it is necessary to use an iterative procedure for this special problem.

The classical procedure for solving the implicit equation or the system of implicit equations is the predictor-corrector scheme. The predictor linear multistep formula, usually an explicit one, give a starting value for the iterative process for solving the implicit corrector formula. The approximation delivered by the predictor formula is improved applying the corrector formula. The standard iterative process for solving the implicit corrector equation is the method of simple iterations, in the nonstiff case, and the Newton's method, in the stiff case.

The linear multistep formula can be applied only if we know the values y₀ , ... ,y _{k - 1} , referred to as the starting values. These are delivered by a procedure which is independent from the basic formula, a starting procedure. Naturally, the stepnumber of such a procedure is smallest that of the multistep formula. In almost all present numerical codes, the starting procedure includes some onestep methods.

Let the linear operator associated with a linear multistep formula L=S _j=0 ^ka_jH^j-hS _j=0 ^kb_jH^j d/dt, where H is the displacement operator defined as Hy(t)=y(t+h). L has the precision degree p, if L~0 for any t^j, j=0,...,p, and p is the maximum integer that satisfies this property. A given linear multistep formula has p order if the associated linear operator has a precision degree p. The linear multistep formula is consistent if its order satisfies p>= 1.

The algebraic system , which is equivalent with the condition of p-order, is the following:

Moreover,

For a linear multistep with k steps, the inequality p<= 2k is referred to as the order barrier.

The methods with k steps are generated by polynomials pairs (r,s) with real coefficients, where

The consistency condition can be written by r(1)=a_k=0, r'(1)=s(1)=b_k=a_{k - 1} /2# 0. The local error of a p order formula can be represented in the form y(t_n)-y_n=c^*h ^p+1 y_n ^(p+1) + O (h ^p+2 ). The value c^*=-c _p+1 /s(1) is referred to as the error constant of the linear multistep formula (r,s) of p>1 order.

The stability study of the linear multistep methods is based on the concept of root condition. A polynomial satisfies the root condition if all its roots are inside the complex closed unitary disk and those from the boundary are not roots of the derivative polynomial. If a linear multistep method (r,s) is stable, then r(x) satisfies the root condition. The reverse implication is true when the linear multistep formula is also consistent. If r satisfies the root condition, we say that the formula is zero-stable.

In the construction of efficient numerical codes, the Dahlquist 's barriers are very important. These are the theorems which impose some restrictions on the stable formula class. The first Dahlquist 's barrier is the following theorem: the order of a consistent and zero-stable linear multistep formula can not exceed the following values

The absolute stable methods are particular cases of zero-stable methods. The idea is to request the decrease of the secondary numerical solution with the increase of the stepnumber. The study is made on the scalar test equation y'=l y, Re l<0. The absolute stability is the minimum property that must be imposed on any numerical integration method.

The behaviour of the method on the test equation is a prevision model of the behaviour in the case of some nonlinear systems. Let the equation produced applying the method with constant stepsize h to the test equation:

The characteristic equation can be get from this formula replacing the values y _n-k+j with x^j: S _j=0 ^ka_jx^j-qS _j=0 ^kb_jx^j=0, q=l h. The polynomial r-qs is referred to as the characteristic polynomial of the linear multistep formula.

A method is absolute stable for a value q=hl if, for the value q, the characteristic polynomial satisfies the root condition. A consistent linear multistep formula is absolute stable at q=hl in C iff all the solutions y_n of the above equation, produced applying the formula to the test equation y'=l y with constant stepsize h, are bounded when n goes to infinity.

The absolute stability domain of a linear multistep formula is the set

The stability domains of the classical methods, like the Adams-Bashforth, Adams-Moulton, Nystrom or Milne-Simpson 's formulae, are almost all bounded. The stability requirement leads to the condition that the values |hl_i|, where l_i are the eigenvalues of the Jacobian 's matrix at the moment t_n, must be inside the stability domain. A large value |l_i| implies a small value h. If the stability domain is unbounded, for example includes the complex left half-plane, then the condition hl_i in A does not restricted the stepsize.

A formula is A-stable if the absolute stability domain includes the complex left half-plane, i.e. those values z for which Re (z)< 0.

For an onestep method, the A-stability condition is derived from a property of the stability function R(q)=y _n+1 /y_n. An onestep method is A-stable if the stability function satisfies |R(z)|<= 1, \forall Re(z)<0..

A linear multistep formula is stable at infinity if lim _{hl to -infinity} |y _n+1 /y_n|=0, where y_n is the numerical solution of the scalar test equation. The formula is L-stable or stiffly A-stable if it is A-stable and it is stable at infinity. A linear multistep formula is strongly A-stable if it is A-stable and there is a value w>0 such that sup _hl<-w |R(hl)|<1.

The class of A-stable linear multistep formulae is not very large, since such a formula must be implicit and with an order of accuracy lowest than three. A linear multistep formula, with a stability domain which includes the negative real half-axis, must be an implicit one. The following result, due to Dahlquist , gives a full characterization of the restrictions impose by the A-stability condition. (it is referred to as the second Dahlquist 's barrier): the order of an A-stable linear multistep formula can not exceed two. The principal ways to break the Dahlquist 's second barrier are:

• weakling the stability condition: stiff stability, A(a) stability and so on;
• using the implicit Runge-Kutta 's process;
• using the high derivatives of the solution:
• multiderivative multistep formulae;
• multiderivative Runge-Kutta 's process;
• using some nonlinear interpolators:
• nonlinear multistep methods;
• rational Runge-Kutta 's process;
• adding some new stages, extradivision points or future points:
• cyclic linear multistep methods and composed methods;
• hybrid methods (extradivision points) and multistep Runge-Kutta methods;
• extended methods (future points);
• reducing the approximation error:
• exponential fitting;
• extrapolation methods;
• exploiting the special structure of the system: stiff separability.

A linear multistep formula is A(a)-stable , with 0<a<p/2, if all the solutions of the difference equation, produced applying the formula to the scalar test equation y'=l y, converge to zero when n\to \infty, for any fixed stepsize h>0, and for any l# 0 so that hl in C _a = {q||arg(-q)|<a,q#0 } holds. The algebraic condition is that, for any q in C _a , all the roots of the characteristic equation satisfy |x_i|<1, i=1, ... ,k.

The formula is A⁰-stable if, applying it to the test equation with l a negative real number, we get some approximations y_n with the property y_n goes to 0 when n goes to infinity, for any constant stepsize h. The formula is A₀-stable if the characteristic polynomial satisfies the root condition for any negative real value q.

Let a linear multistep formula and x_i, the roots of the characteristic polynomial. Then

is referred to as the relative stability domain, and also as the star region.

A linear multistep formula is stiffly-stable with respect to the D, a, t parameters, if R₁ union with R₂ is subset of A, R₃ is subset of S where

Any numerical method must satisfy the condition of convergence to demonstrate its efficiency in the integration of an initial value problem, since the convergence supposes a control of the error produced by the numerical process. The global error (accumulated error) of a numerical method is e_n=y(t_n)-v_n, where y(t_n) is the exact solution, and v_n is the approximate solution at t_n, produced by the numerical code.

A method is convergent if lim _{h goes to 0} max_n||e_n||=0 holds for all the initial value problems with the standard conditions for the solution unicity.

Note y_h(t) the interpolation function of the approximate values, produced applying a linear multistep formula at an integration interval division with constant stepsize h. A linear multistep formula is convergent if ||y(t)-y _h (t)|| goes to 0, when h goes to 0, t in [t₀,t₀+T], holds for any initial value problem with the standard conditions for the unicity of the exact solution, independent from the starting values, which only satisfy ||y(t₀+ih)-y_h(t₀+ih)|| goes to 0, h goes to 0, i= 0,...,k-1 . The algebraic conditions for convergence are the zero-stability and the consistency.

There are three ways to build-up a linear multistep formula:

• by numerical integration; for example, we can mention the Adams-Bashforth (for a_ka _k-1 # 0, a _k-2 =...=a₀=b_k=0), Adams-Moulton (for a_ka _k-1 b_k# 0, a _k-2 =...=a₀=0), Nystrom (for a_ka _k-2 #0, a _k-1 =a _k-3 =...=a₀=b_k=0), Milne-Simpson (for a_ka _k-2 b_k#0, a _k-1 =a _k-3 =...=a₀=0) or Newton-Cotes (for a_ka₀#0, a _k-1 =...=a₁=0) formulae;
• by numerical diferentiation; for example, we can mention the backward differentiation formulae (for a_kb_k# 0, b _k-1 =...=b₀=0);
• by combination of the coefficients of some known formulae.

A bounded stability domain is a characteristical property for the methods of the above mentioned first way. Consequently, these formulae are not suitable for the integration of stiff problems. Most adequate are the following methods:

• the backward differentiation formulae (short noted BDF). For k<= 6, the methods with k steps and maximum order k are known as Gear 's formulae. The Gear 's formulae are A-stable for k<= 2 and stiffly stable for k<= 6. The maximum order for a stiffly stable formula is 6.
• the Varah 's formulae reduce with an unity the method order, in the favour of a greatest stability angle than that of Gear 's formulae. The Varah 's formulae are stiffly stable for k<= 12. The maximum order for a stiffly stable formula is 8.
• the Cash 's extended BDF suppose an implementation in a predictor-corrector scheme: the predictor formula is one of the Gear 's methods, and the corrector formula is an extended one:
  (P) S _j=0 ^k a'_jy _n+j-k =hb'_kf_n, (C) S _j=0 ^ka_jy _n+j-k =hb_kf_n+hb_k+1 f_n+1 .

  The scheme has k+1 order of accuracy. The purpose of using this scheme is to improve the maximum order of the stiffly stable formulae, and to maximize the stability angle. The Cash's extended BDF schemes are A-stable for k=1,2 and stiffly stable for k<= 8. The maximum order for a stiffly stable formula is 9.

• the Cryer 's methods are defined by the polynomial s(e)=(e+d)^k, d in (-1,1) and the condition of maximum order of accuracy. The Cryer 's method with k steps is A(0)-stable and stiffly stable if -1 < d <= -1+2/(1+2 ^k+1 ), where k=1,...,7;
• the Grigorieff-Schroll 's methods are computed using the polynomial s(e)=(e+d) ^k-1 (e+c) and the condition of maximum order. The numerical tests have proved that, for some values d and c, the methods are stiffly stable for k<= 12 steps.
• The split linear multistep methods have been proposed by Cash. A split formula with k steps has the form
  S _j=0 ^k a _jy _n+j-k -hS _j=0 ^k-1 b _jf _n+j-k =ht f(t, z_n)+h(b_k-t)f(t_n,y_n)

  where z_n is a value generated by a predictor formula of the form

  (P) S _j=0 ^k a'_jy _n+j-k =hS _j=0 ^kb'_jf _n+j-k

  Cash has studied the effect on the stability of splitting the BDFs. The Cash 's split BDF schemes are A-stable for k=3 and some t values, and stiffly stable for k<= 6 and some t values. The maximum order for a stiffly stable formula is 6.

• The Dill 's formulae lie in the iterative process
  S _j=0 ^ka_jy _n+j-k =h(b_kf_n+b _n-u f _n-u ), b_k=1, u in {1, ... ,k}

  The purpose of building these formula is to increase the order of stiffly stable formulae. The Dill 's class of formulae of p=k-1 order includes some methods of 7 and 8 order which are stiffly stable. The maximum order for a stiffly stable formula is 8.