On the development of a type-insensitive code

Advances in Engineering Software 1992, 14, 119-127

On the development of a type-insensitive code K.C. Wade*, M.G. Everett & C.W. Richards~ Centre for Numerical Modelling and Process Analysis, Thames Polytechnic, London, SEI8, UK

In this paper the need for a numerical algorithm that can automatically select, not only the step size and order, but also the most suitable numerical Ordinary Differential Equation (ODE) integrator from a prescribed set, is demonstrated. Such an integrator based upon explicit and implicit Runge-Kutta methods is described along with all the required switching strategies. Finally, the performance of the new algorithm, SARK (Switching Algorithm for Runge-Kutta methods) is assessed over a number of test cases.

stability function, f(qi = 2ih), i = l(l)m is generated. This stability function can then be used to test for the property of A-stability or A(~)-stability which a numerical method must possess for it to be able to solve stiff problems. This places a severe restriction on the numerical method as it must be implicit. 2 As most stiff systems have a mixture of fast and slow transients, the initial steps must be taken with a small step size to ensure that the fast decaying transients are followed faithfully. Thus an implicit integration method which utilizes a larger number of function evaluations on every step, will be uneconomical over the initial range. As the integration continues, however, and the fast transients decay, a larger step size can be employed by the implicit methods whereas the explicit method must continue, for stability restrictions, to use the smaller step size. For oscillatory problems, it has been widely upheld that a precisely A-stable method should be employed, as these types of problems are usually characterized by imaginary eigenvalues of the Jacobain matrix Of/~y.3.4,7.10 However, Wade ~ has shown that explicit methods can be more economical and as accurate as implicit methods for this type of problem.

INTRODUCTION This paper is concerned with the numerical integration of a system of m ODEs, y' =

f(x,y)

y(0) =

x >~ 0

(1)

a

at a series of discrete points x~, i = I(1)N. It is often convenient to classify the above system as either a: (classically) stiff, b: oscillatory or, c: non-stiff, depending upon its characteristics. 5 Problems classified as being non-stiff are generally numerically the simplest to integrate. Stiff problems generally occur when the system of equations possess a mixture of fast and slow decaying transient components and the step size of the numerical method is restricted due to the fast transients. The category of oscillatory cases includes such problems as simple harmonic motion, damped harmonic motion, etc. It is generally recognized that to solve stiff problems numerically a method must provide the correct behaviour when applied to the standard test case, y' = y(0) =

Ay

x ~> 0

NUMERICAL INTEGRATORS

(2)

In solving a particular problem efficiently, care must be taken to select the correct integration method. Unfortunately, it is not always possible to know in advance the nature of the solution to a particular problem and as mentioned in the previous section it may be beneficial to use different methods for different regions of the integration range. The aim of this work is to develop an algorithm capable of detecting the nature of the solution (stiff, highly oscillatory or non-stiff) and automatically

1

After applying the numerical method to the system (2), where A is an m*m diagonal matrix of eigenvalues, 2~, a Now at *Baker Jardine and Associates, 19 Heathmans Road, London SW6 4TJ, UK ~:Fluid Power Centre, University of Bath, Bath, Avon BA2, UK. Advances in Engineering Software 0965-9978/92/$05.00 ~ 1992 Elsevier Science Publishers Ltd. 119

120

K.C. Wade, M.G. Everett, C.W. Richards

select the most appropriate integration method. Such methods are known as type-insensitive codes and have obvious attractions as general ODE solvers. The fact that highly oscillatory problems can be solved adequately using explicit methods leads to a useful simplification of the task in hand. The methods chosen for the non-stiff/oscillatory problems were the highly popular one-step explicit Runge-Kutta methods (ERK), of the general form, Y,+l

(3a)

= y, + h ~ ciki

exact)] ~ The mechanism chosen, therefore, is the well proven and reliable, although computational expensive, Richardson extrapolation, i,e. taking three steps, two of step size h and one repeated step of size 2h. Although this method is expensive in terms of function evaluation per step, the use of it in this case does have the additional advantage that it allows the switching strategies described later to work efficiently.

OVERALL SWITCHING A L G O R I T H M

i-I

k~ = f

x. + b~h,y. + h ~ a u

i = l(1)s

i=[

(3b) for an s-stage pth order method (p ~< s). To complement the explicit Runge-Kutta methods, suitable, implicit methods capable of integrating stiff problems are required, Cash ~ devised a class of RungeKutta method known as mono-implicit Runge-Kutta methods (MIRK), and it is a class of these that we shall focus upon as they possess the correct properties that we require. We shall refer to these methods as Backward Runge-Kutta methods (BRK). These BRK methods have the general form, Y,+I

= y, + h ~ cik i

ki = f

i=l

(4)

x,+, - bih, y,+. - h Y, a u

i = l(l)s

Often in engineering the use of numerical integrators is restricted to that of a 'black box' whereby the ODE to be solved is 'plugged-in' and the answer generated. This means that the engineer does not want to go through the laborious exercise of calculating Jacobain matrices, etc. For this reason the algorithm described works on the numerical evaluation of the Jacobain matrix by finite differences. Also, all the switching strategies are performed numerically rather than, say, by calculating eigenvalues during the integration range. This may appear to be detrimental to the overall scheme, but it does allow it to be used in a more flexible manner and, as will be shown below, the overheads incurred are very small. To ensure that a wide variety of problems at different accuracy levels can be solved the algorithm is also capable of switching order, a complete schema is presented in Fig. 1, showing all the possible switching opportunities.

j=l

for an s-stage pth order method. The main advantage of using BRK and E R K methods is the close coupling between the two methods; we note that the coefficients, (a and b) and weights e r are the same for an s-stage ERK method and an s-stage BRK method, e.g. Euler (ERK) and Backward Euler (BRK). The methods stability functions are also complementary. ~ If the stability function of the explicit method is Fe(2) then the stability function of the corresponding backward method is I/F~(- 2).

EBRK4

V

4,-[

ERK 1/2 ~.

-4"

~BRK3 ~-;tep Fail ! ~Fail l; ~ -! ~ EBRK3- ~t__ . ! ~ IERK2/3

~SStep

ERROR CONTROL The most commonly employed error control mechanism used on ERK methods is the embedding method and this is used here. Thus approximations are taken with an order p method and an order p + i method to estimate the error produced in the lower order method, which is, providing that the step is successful, carried forward. Most practitioners, however, carry forward the approximation produced by the higher order method and the same practice is adopted here. Unfortunately, this type of error estimate can not be used for the BRK methods as the error that is predicted is at Yn (assuming that y,+~ is

RK1/2

i

i /

Fig. I. Schema of the overall algorithm.

On the development of a type-insensitive code

121

and for the (6, 1) method

SWITCHING FROM EXPLICIT METHOD TO IMPLICIT METHOD

6

)',,+1 = Y, + h ~ e~ki The algorithm starts the integration with an explicit method using a very small step, as this is computational the cheapest and the state of the equations are initially undetermined. The explicit method will be efficient providing that the system remains non-stiff or oscillatory. If the equations are, or become, stiff then the numerical integrator must be changed to a stiff solver, hence some means of detecting stiffness is needed. The approach taken here is to use a stiffness detection scheme based on an idea by Shampine & Hiebert. 9 Assuming that the integration is being performed by an s-stage explicit method of order p, (s, p), and the error controlled by an embedded (s, p + 1) method, then after a step has been taken the error control mechanism will select the next step size based upon the two estimates generated and the requested accuracy level. This new step will be the largest possible in order to meet the requested accuracy requirements and may correspond to a q~ value being on the boundary of the methods stability region, in which case, the step size is restricted for stability reasons. Alternatively, all qg values may be well within the stability region and the step size is restricted by accuracy. It is on this basis that the stiffness detection scheme works. The main requirement, apart from reliability of the stiffness detection scheme, is economy. As the scheme is to be employed after every successful step the cost of it must be minimal. Consider taking a step with the (6, 4) and (6, 5) Fehlberg methods (eqn (3) with the coefficients and weights given by the Butcher matrix, Table 13. At each step six function evaluations at discrete points between xo and x,+~ are available. From these a (6, 1) and a (6, 2) method can be constructed, i.e. 1st and 2nd order methods with 6 stages. This is clearly inefficient in a normal explicit scheme (as s-order s-stage methods can be found for explicit method with s <~ 5), but in this case the function evaluations and k values are already available from the (6, 5) method and hence no extra work is involved. Thus, a different linear combination of them can be taken to form the lower order methods. For the (6, 2) method this will be,

where the weights are calculated below and taken with the coefficients from Table 1 to form the required order method. The two equations to be satisfied for a 6-stage 2nd order method are

6

y.+,

=

(5)

y. + h Z diki

(6)

i=l

t=]

6

Z di i

=

1

(7)

1

and 6

Z dibi

=

1/2

i=1

By incorporating the coefficients from the main (6, 5) method, the six free parameters can be reduced to four with

d2 and 5

a,=l-Z4 i=2

In a similar manner one of the free parameters can be removed for the 1st order method. The only requirement of this lower order pair is that their stability regions are uniformly larger than that of the main integration method. A computer search was made for a (6, 13 and a (6, 2) pair which has a stability region uniformly larger than the (6, 5) method. The weights produced are given in Table 2. The boundary of the absolute stability regions of the 5th order related methods are shown in Fig. 2. Clearly, the stability region (the enclosed finite regions) of the embedded explicit methods, ERK 2 and ERK 1, used for the switching of numerical methods, are uniformly larger than the main explicit methods. Similarly for a (3, 2) method with a (3, 3) method for error control, embedded (3, 1) and (3, 2) methods for stiffness detection can be found. This idea of embedding a 2nd and 1st order pair within the higher order methods for the detection of stiffness,

Table 1. Coefficients for (6, 5) and (6, 4) Fehiberg methods

0 1/4 3/8 12/13 1

1/2

0 1/4 3/32 I 932/2 197 439/216 - 8/27 16/135 25/216

9/32 - 7 200/2 197 - 8 2 0 0

7 296/2 197 3 680/513 - 3 544/2 565 6 656/12 825 1 408/2 565

- 845/4 104 1 859/4 104 28 561/56 430 2 197/4 104

- 11/40

- 9/50 --1/15

2/55 0

(6, 5) (6, 4)

122

K.C. Wade, M.G. Everett, C.W. Richards

5th

order

related

stability

PSPLOT Version 3.58b

regions

PIPESIM Plotting Utilitz -

.

o~l I"" "

4~ .%.,%, | If$#"

k

i

o.° oO

ERK 1 .,"

ERR 5

-"

....... ......

• •

. .. ,.:':

,o oO ,o oO

tm

oO°

o,l,,ot

" "'"ERK 2 ,'.-

.*e°°''°

•

o

•

.."."" ERK 4

oO • oOJo. .s

/J oO

2: "s"

g ,%

oo

4:

li

• • o" ;

O[ . . . . 6

-5

4

-3

t

Real

J

-2

-1

q~

Fig. 2. Stability regions of 5th order related methods. can be extended to higher orders. It cannot, however, be used in conjunction with any E R K pair lower than 3rd order. For this scheme to work a 2nd and 1st order embedded pair are required with stability regions uniformly larger than the main method, and if the main method is 2nd order, in 2 stages, then the embedded 2nd order method must have exactly the same stability region as the main method. Thus, a switch from E R K 2 to BRK 2 cannot be provided in this manner. The main integrator pair will select a step size that just meets the accuracy requirements. I f the step size is restricted on the grounds of accuracy rather than stability, then the lower order embedded pair are unlikely to satisfy the accuracy requirements. We know that the lower order pair are failing for accuracy reasons, and not for stability reasons, as its stability region is larger than that of the main method. Thus, the system is not stiff at the current integration point. If the lower order embedded method is repeatedly able to meet the same accuracy requirements, then the step size of the higher order method must be restricted for stability reasons. This implies that the problem is becoming stiff. If the embedded method is successful for 50% of the time over 50 consecutive steps then the

problem is deemed stiff and a switch to the implicit method advocated.

SWITCHING FROM IMPLICIT METHOD TO EXPLICIT METHOD

In order to ensure that the code is competitive, the implicit to explicit switching must also be computationally inexpensive. One approach, 6'8 is to use the approximation to the Jacobian matrix of the problem, or the exact Jacobian matrix if this is supplied to the code, to form some estimate of the Lipschitz constant. However, in the case of BRK methods this matrix is not computed. (An approximation to the iteration matrix is computed, but this is different from the jacobian). The approach taken here is to generate a non-stiff detection scheme based on the function evaluations available from the implicit method. These will be the k values after the implicit equations have been solved to an acceptable tolerance by the quasi-Newton process. Thus, the detection scheme will be a numerical method in its own right. The only requirements of this new method is that it

Table 2. Weights for (6, 2) and (6, 1) embedded pair

-0-139682 0084 227

--0-198633 - 0" 163 140

0.724462 0.761 013

0.428953 0.405 846

--0"141 85 - 0"319 70

0-047041 0.044 024

(6. 2) (6. 1)

On the development of a type-insensitive code Table 3. Weights for embedded 3rd order EBRK

uses the k values from the implicit method and must have characteristics similar to an explicit method. In particular it must possess a finite stability region. Taking a linear combination of the k values produced by an implicit BRK method, clearly cannot result in a method with a finite stability region. If this were so then the corresponding E R K method would have an infinite stability region! After a successful step with any Runge-Kutta method the k values are discarded and a new collection generated on the next step. In fact, the only value carried forward is the solution y,+~, which becomes y, on the next step. By using previous solution values, ( y values), the resulting method becomes akin to linear multistep methods and will inherit some of their deficiencies. If, however, the k values computed on a step are stored and used on the next step then a new method, an Extended BRK method (EBRK), is generated, which can possess a finite stability region. Thus,

y,+, = y , + h £ k i + h ~ k ; i=l

1 2 3

(9)

is used to advance the step from y, to y,+~ instead of (3). In (9) the k~ are selected from the k i, i = l(1)s values from the previous step. Thus, the k values in generating a solution at y° are stored and used in the generation of the solution at y,+~. Clearly, this technique must be employed on the second h step of the Richardson error control process. The solution produced by the EBRK can only be

Order

ci

di

1-01000 -2.0320 - 0.676 7

4.2167 - 1.5180

acceptable providing that qi, i = I(1)N is within its stability region for all i such that Re(qi) < 0. The solution of the EBRK method is thus compared with that obtained by the BRK over the second h step and is regarded as successful if it satisfies the same error tolerance as is employed for the BRK method. By being successful it implies that an ERK method, which has a similar stability region, would be capable of continuing with the integration of this non-stiff problem. If the solutions agree for five consecutive steps then q~ is deemed to be 'small' and a switch to the cheaper explicit method performed. The weights for the EBRK method were found numerically by a computer search. Clearly, the weights are constrained by the order equations that must be satisfied for a particular order, i.e. for a 3rd order method there are four equations to satisfy. The successful weights being those that would produce stability regions that were finite and comparable with the stability regions of the corresponding E R K method. With t = 2 a 3rd order EBRK method with stability region shown in Fig. 3 is

i=l

3rd

123

related

stability

2.5 PSPLOTVersion 3.58b

regions

PIPESIM Plotting Utility •

•.,

•."

..!::~. t.

2.0 • ••• • •

t~

. • •

o

~

•

-"

:

t

ERK 2

•.

BI~

° ~

•

,,

1.5

• •°l°s

|

3

"•-t ..

~

~

1.0

0.5 • EBRK

3

:

*

O.

-:3.0

-2.5

-2.0

.5

.0

Real

.5

ql

Fig. 3. Stability regions of 3rd order related methods.

0.5

K.C. Wade, M.G. Everett, C.W. Richards

124

Table 4. Weights for embedded 4th order EBRK

l 2 3 4 5 6

0.000 00 -0.91475 3.37465 0.000 00 - 0.299 77 - 1.793 42

0.453 71 0.07037 0.167 89 - 0.058 68

produced, which is comparable with the stability region of the 3rd order E R K method. The weights for this are given in Table 3. For the 5th order method a 4th order E B R K was initially tested to determine if the same idea would work for higher orders. Suitable weights were found, using t = 4, Table 4. Initial testing was performed with this 4th order EBRK. The method proved very reliable and was always able to detect when the problem became non-stiff. Ideally a 5th order E B R K method should be used. However, no suitable method has yet been found and since the 4th order method proved so effective it was retained. Numerical testing with a 2nd order E B R K and the 3rd order E B R K for use with the 3rd order B R K demonstrated that the order of the E B R K need not match exactly the order of the B R K method. However, as a 3rd order E B R K method could be generated it ws retained for use with the 3rd order B R K method. Similarly an E B R K for the 2nd order BRK method is required. This can be generated by the use of one past value, t = 1. One suitable choice of weights is given in Table 5, with the resulting method being 2nd order.

SWITCHING ORDER

To meet the requirements of the user and to make the code more flexible it should ideally be able to select the most appropriate order. The orders employed for this code are 2, 3 and 5 for the explicit methods and orders 1, 2, 3 and 5 for the implicit ones. The main orders in the code are 3 and 5 with the others playing a supporting role in case of serious failures. The 2nd order E R K method being employed mainly so that the algorithm can recover from the use of lower order B R K methods, e.g. being forced to low order if the iteration matrix is regularly computed as singular. This is a problem that occasionally appears when solving some linear problems. Table 5. Weights for embedded 2nd order EBRK

1

2

Ci

~Ji

0.5 -0'5

1-0

Table 6. Weights for embedded 2nd and 3rd order methods

I 2 3 4 5 6

3rd order

2nd order

0-083 1292 - 0.002 969 8 0.618 7300 0.1847600 0.080 000 0 0-036 363 6

0.003 1290 0.007 032 0 0.128 5200 0.1234860 - 0-065 661 0 0.803 4940

ORDER REDUCTION

As the decision to switch order must be inexpensive, the k values produced by the main integrator must be used. If the integrator has order p then clearly, a method of order less than p can be embedded within this method such that the same function evaluations are used. When either the E R K or B R K 5th order method is being used a decision must be made as to whether to switch down to the 3rd order method. By embedding a 3rd order method within the 6-stages of the 5th order the comparison can be made. For the E R K method an embedded 3, 2 pair will be required for the error estimation while for the B R K method the same 3rd order method used in backward mode can be employed. The weights were again determined numerically, the final chosen values given in Table 6. The embedding error test, for the E R K method, is performed by generating a 3rd order and a 2nd order approximation at y,+~ after the 5th order E R K method has taken a successful step. A Fehlberg type test is then carried out to determine if the lower order method can match the required accuracy. For the BRK method the soluton obtained from the lower order method, over the second h step, is compared directly with the B R K solution. The same embedding process can be used for a 2nd order, plus 1st for the E R K method, within the 3rd order methods.

INCREASING ORDER

All of the switching strategies developed so far rely on information readily available at the end of a successful step. The only extra expense involved is a few multiplications and additions. Various ways of increasing order were considered that involved computing extra function evaluations. Only two possibilities came to light that would not require any additional work. The first was to use another extended method, i.e. using previous k values. However, it proved difficult to generate the E B R K methods with the required stability and accuracy properties. When employing the 3rd order method the proposed change is to 5th order, hence the extended method will require at least 6-stages. Even by taking all

On the development of a type-insensitive code

125

T a b l e 7. R e s u l t s f o r n o n - s t i f f p r o b l e m p l

Method

Order

Log TOL

FE

JE

Steps

Sig Fig. Accuracy

CPU Time

ERK

3

ERK

5

BRK

3

BRK

5

3 - 6 - 3 - 6 - 3 - 6 -3 -- 6 - 3 - 6

108 849 126 354 468 1374 1566 170 4 126 354

14 10 14 13 0 0

36 283 21 59 68 282 106 90 21 59

4.16 7.08 4.85 8.31 4.37 6.61 3.71 7"09 4"85 8'31

0.21 1-07 0.17 0.40 0.63 1.98 1.73 1.74 0'20 0.54

SARK

-

the k values from the previous step a 5th order method with the correct stability region is not guaranteed. The E B R K 4th order method required 4 previous values, 10-stages in total to produce the required stability region. The much simpler option is to continue with the current order method until the step fails for accuracy reasons. At that point an order increase is advocated. This method was employed for both the E R K and B R K methods. If the E R K method is being used, then a step size failure must be due to accuracy if the stiffness detection scheme was not triggered. Thus, the stiffness prediction scheme must be given priority over order change. This very simple order increase scheme is justified purely on the grounds that numerical tests indicate that it works reasonably well.

NUMERICAL RESULTS The numerical results are split into two sections to test the correct behaviour of all the switching strategies described earlier. The headings for each table of results depict: Order of the numerical method; Logz0 of the requested tolerance; - N u m b e r of function evaluations; - N u m b e r of Jacobain evaluations; -

- Number of steps required; Significant figures accuracy of the components Y3, Y3, Y4 and y~ for problems pl, p2, p3 and p4 respectively at the end of the integration range; - Total C P U time in seconds on a Prime 550 minicomputer running under PRIMOS. -

INTEGRATING

NON-STIFF

PROBLEMS

In this section the ability of the code to select methods for non-stiff problems is tested. Two problems are considered, problem pl and problem p2. Clearly in these problems stiffness should not be detected and thus the explicit method should be used throughout the integration range. The results of integrating these two problems with ERK3, ERK5, (explicit methods of order 3 and 5), BRK3, BRK5 and S A R K are shown in Tables 7 and 8. The most efficient method for p 1, Table 7, is clearly the high order explicit method with the implicit methods being inefficient. Clearly, S A R K chooses the correct method of integration by staying with the high order explicit method. The C P U times produced by S A R K are higher than the E R K 5 because of the slight overheads incurred in the switching strategies. For problem p2 a similar pattern emerges with the explicit methods being far superior to the implicit ones. Again, S A R K can correctly deduce that the problem can

T a b l e 8. R e s u l t s f o r n o n - s t i f f p r o b l e m p2

Method

Order

Log TOL

FE

JE

Steps

Sig Fig. Accuracy

CPU Time

ERK

3

ERK

5

BRK

3

BRK

5

- 6 -9 - 6 - 9 - 6 - 9 - 6 -9 -9

429 4 101 126 366 843 2 937 I 308 1 980 366

13 13 19 15 0

143 1 367 21 61 90 414 68 21 61

2.16 5.58 1.79 4.51 1-21 4.52 3.16 5.33 4.51

0-70 6.77 0.20 0.53 1.23 4.70 1-80 2.66 0.55

126

K.C. Wade, M.G. Everett, C.W. Richards Table 9. Comparison of BRK 5 and SARK (starting on BRK method) on problem 1

Method

Order

BRK SARK*

Log TOL

FE

JE

Steps

Sig Fig. Accuracy

CPU Time

% ERK

- 3 -6 - 3 - 6

1566 1 704 702 1290

14 13 14 13

106 90 42 71

3.71 7.09 7'60 9.59

1.73 1.74 108 169

90.4 561

Table 10. Comparison of BRK 5 and SARK (starting on BRK method) on problem p2

Method

Order

Log TOL -

BRK SARK*

3 6 3 6

FE

JE

Steps

Sig Fig. Accuracy

CPU Time

% ERK

1260 1308 576 744

25 19 13 10

44 68 31 53

0-72 3.16 1-14 3.22

1.72 1.80 0.99 1.26

90.5 81.8

be solved more efficiently by an explicit method and thus no switching occurs. (Table 8). The other occasion when a non-stiff integrator is required is when a stiff problem becomes non-stiff. This can be simulated by commencing the integrating of a non-stiff problem with a stiff integrator. The detection scheme should detect that the problem is non-stiff and switch accordingly. Tables 9 and 10 display results that compare BRK5 with SARK, run in this manner, on problems pl and p2 respectively. By commencing the integration of SARK with the 5th order implicit method it clearly switches over to the 'cheaper' explicit method when the problem is deemed to be non-stiff. For pl the ERK methods are used for 90% of the integration range at a tolerance of 1.e-3 and for 56% at a tolerance of 1.4-6. On problem p2 these percentages are even higher at 91% (l.e-3) and 82% (1.e-6), clearly showing that lack of stiffness is detected and more importantly acted upon.

ODEs. The switching of order is also tested. Clearly no explicit method is able to integrate a stiff problem efficiently and hence S A R K must switch over to an implicit method. The results of running the stiff problems p3 and p4 are shown in Tables 11 and 12. For problem p3 the low order implicit method is the most efficient, hence SARK should switch over to this method, which it does, producing comparable accuracy in a faster time. Problem p4 on the other hand is integrated quicker by the high order method at a tight tolerance and by the lower order method at a high tolerance. S A R K manages to switch to the appropriate method and produce superior results to either of the implicit methods at both tolerances. For both problems the BRK methods were used for 99.9% of the integration. These results show that, not only is the switching of numerical methods handled correctly, but also the most appropriate order is selected.

INTEGRATING STIFF PROBLEMS

CONCLUSION

This section examines the ability of the code to select the correct method when confronted with a stiff system of

These four problems demonstrate that S A R K is capable of correctly selecting the most suitable method available

Table 1 I. Results for stiff problems p3

Method

Order

Log TOL

FE

JE

Steps

Sig Fig. Accuracy

CPU Time

% BRK

BRK

3

BRK

5

-3 - 6 -3 -6 - 3 - 6

516 1950 2 136 2958 585 1 521

12 27 36 36 12 12

58 346 132 222 113 301

3.18 6.45 3.77 4.57 3-04 6.43

0"75 2.96 229 338 150 2"66

99.9 99.9

SARK

127

On the development o f a type-insensitive code Table 12. Results for stiff problems p4 Method

Order

Log TOL

FE

JE

Steps

Sig Fig. Accuracy

CPU Time

% BRK

BRK

3

BRK

5

-4 - 7 -4 - 7 -4 -7

6 153 31329 50496 67158 1 545 11 316

33 39 87 63 14 45

1288 6850 2 690 2 828 245 1691

1.31 3.25 1.68 4.26 1.47 3.90

7.68 69.19 43.60 54.93 1.82 11.35

99.9 99.9

SARK

to it and that it is often superior to specialized integrators. Clearly, when the characteristics o f the problem are such that it can be easily categorized as stiff or non-stiff, then a specialized integrator m a y be superior. However, it is often difficult to categorize the problem and furthermore m a n y users o f codes are not interested in attempting to do so. In such cases S A R K is to be recommended.

il. Wade, K.C. SARK - - A type-insensitive Runge-Kutta Code, PhD Thesis, CNAA, Thames Polytechnic, 1988.

PROBLEMS

CONSIDERED

Problem: pl

(I + x):

ACKNOWLEDGEMENT y~

The authors would like to thank D r David Sayers o f N A G Ltd, Oxford for his assistance in making available relevant subroutines from the N A G library required for the preparation o f this paper.

=

y2+y~(l

y' 3 ~ x

y~

=

y(O)

+x) 2

[o, lo]

=

=

[1, l , l ] T I

REFERENCES

Problem: p2 1. Cash, J.R. A class of implicit Runge-Kutta method for the numerical integration of stiff ODE, J. ACM, 1975, 22, 504-511. 2. Dahlquist, G. A special stability problem for linear multistep methods, BIT, 1963, 3, 27-43. 3. Gear, C.W. Numerical solution of ODEs: Is there anything left to do?, SIAM Rev., 1981, 23, 19-24. 4. Jeltsh, A.A. Stability on the imaginary axis and A-stability of linear multistep methods, BIT, 1978, 18, 170-174. 5. Lambert, J.D. Computational Methods in ODE, Wiley, London, 1973. 6. Norsett, S.P., & Thomsen, P.G. Switching between modified and fix-point iteration for implicit ODE solvers BIT, 1986, 26, 339-345. 7. Prothero, A., & Robinson, A. On the stability and accuracy of one-step methods for solving stiff systems of ODE, Maths. Comput., 1974, 25, 145-162. 8. Shampine, L.F. Stiffness and non-stiffdifferential equation solvers, Numerisehe Behlanglung v a o n Differetialgleiehungen, ed. L. Collatz, Inst. series Numer. Maths, 27, 287-301, Birkhausser Basel, Switz., 1975. 9. Shampine, L.F., & Hiebert, K.L. Detecting stiffness with the Fehlberg (4, 5) formula, Comp. & Maths. with Appls., 1977, 3, 41-46. 10. Singhal, A. Implicit Runge-Kutta formulae for the integration of ODE, PhD Thesis, University of London, 1980.

Y~ =

YrY3e' Y2 I+X

)"2 -)'~ =

y2(1 + x)e '

x

[0, lO]

=

y(O) =

[l.e-2, - 1 ,

-1] 1

Problem: p3 Problem: pl with/3

=

Problem: p4 Yl =

3'2

y~ =

10001yl - 2), 2

3'~ =

y~(I + x)e ~

x

[0, 10]

=

y(0)

=

[l, l] T

l.e6