Generalized linear multistep formulae
Building up some linear multistep formula of high order for the stiff problem integration we can follow two directions: The low order linear multistep formula proposed in the last chapter can be improved, from the point of view of the error, applying one of the following methods: (i) exponential fitting, (ii) extrapolation.

A second derivative multistep formula uses also the second derivative of the exact solution y"(t)=d(f(t,y(t))/dt. Note g(t):=y"(t)=[(d f/d t)+(d f/d y)f](t,y(t)). The iterative process, which represents the second derivative linear multistep formula, is

S i=0 ka k-i y n-i =h S i=0 kb k-i f n-i +h2S i=0 kg k-i g n-i .

Examples:

A multiderivative linear multistep formula with k steps and m-derivatives is an iterative process of the form

S i=0 mS j=0 k(-1)ia ij hi y n-k+j (i) =0 , n>= k,

where S i=0 ma2 i0 S i=0 ma2 ik # 0. Note lj=max {i| a ij # 0} . We get the standard linear multistep formulae when m=1 and the second derivative formulae when m=2. The onestep methods for k=1 can be rewritten in the form

y n+1 =y n +S i=1 m hi / i! (a i y n+1 (i) +biy n (i) ).

For example, the Taylor 's series correspond to ai=0, i=1, ... ,m, and the formulae

yn=y n-1 +hS k=0 p-1 (-h)k[P (p-k-1) (1)f (k) (yn)-P (p-k-1) (0)f (k) (y n-1 )]

where P(x)= 1 / (2n)! dn / dxn (xn(1-x)n) are A-stable (the Pad\'e 's approximations to the exponential function). The Brown 's formulae generalize the backward differentiation formulae. The general formula is

S i=0 k a i y n+i-k =S j=1 m hj bj fn (i-1) , S i=0 k |a i |>0 , S j=1 m |bj |>0 ,

with b0=-a k , (-1) j+1 bj>= 0, j=0, ... ,m, S i=0 k a i =0, a k # 0, bk # 0. The Brown 's methods have the maximum order p=k+m-1, and the formulae are A0-stable and stiffly-stable for k<= 3 / -2 (m+1).

The block method lies in the iterative process

S j=-m l-1 a ij y kn+j -h S j=-m l-1 b ij f kn+j =0, i=0, ... ,l-1,

where m is the number of starting points, l is the number of points at which we get the solution approximations at one step, and k is the advancing point number with k<= l. We get the standard linear multistep formulae for l=1. The cyclic methods correspond to the case k=l. The set of l methods with m starting values is referred as a l-cyclic m-step method. Examples:

A hybrid method with a single extradivision evaluation has the form

S j=0 k a j y n+j-k =h S j=0 k b j f n+j-k +hb t f n+t ,

where f n+t =f(tn+t h,y n+t ), and y n+t is computed by a predictor formula. For example, the England 's methods are implicit predictor-corrector schemes with a hybrid corrector formula and with stability properties identically to the ones of the Enright 's second derivative formulae. The advantage of using these methods is the replacement of the second derivative evaluation with the approximation of the exact solution at an extradivision point. The general formula is the following:

y n+t = S j=0 k a j y n+j-k +han f(y n ), (P)

y n =y n-1 +h S j=0 k b j f n+j-k +hb f(y n+t ). (C)

The England 's hybrid methods are A-stable for k<= 2, stiffly stable for k<= 7 and of order p=k+2.

The basic idea of the exponential fitting method lies in building up a family of integration formulae with some free parameters, others than the stepsize h. The restrictions imposed by the stability request do not affect the stepsize, only these free parameters. A numerical method is exponentially fitted with the order s at a value hl if it exactly integrates, using a stepsize h, each initial value problem with a solution of the form P(t)e l t , where P is an arbitrary polynomial with the degree deg\, P<= s. The exponential fitting may be accomplished by two techniques:

The basic idea of the extrapolation method is the following. Let a sequence of positive integers n1< n2< n3< ... and the sequence h1>h2>h3> ... related by hi=H/ni. Let a integration step-by-step method of p order. We evaluate, for each i>= 1, the numerical solution of the initial value problem on the interval with the length H, computing ni values at some points of an equidistant division with the stepsize hi. Note T i,0 =y ni . Let the interpolation polynomial

P(h)=e0+ephp+e p+1 h p+1 + ... +e p+k-1 h p+k-1 , P(hi)=T i,0 , i=j-k, ... ,j

We extrapolate the result by evaluating the limit value at ht to 0: P(0)=e0=T j,k . The extrapolation table has the following form:
Tj,0
Tj-1,1
... ... ... Tj,k
Tj-k+1,1
Tj-k,0
The extrapolation method improves the initial approximations with each new column of the above table. The values T j,k are the results of a numerical method of p+k order. Examples: