ADVANCED STIFF SYSTEMS DETECTION

The paper deals with stiff systems of differential equations. To solve this sort of system numerically is a difficult task. There are many (implicit) methods for solving stiff systems of ordinary differential equations (ODE’s), from the most simple such as implicit Euler method to more sophisticated (implicit Runge-Kutta methods) and finally the general linear methods. The mathematical formulation of the methods often looks clear, however the implicit nature of those methods implies several implementation problems. Usually a quite complicated auxiliary system of equations has to be solved in each step. These facts lead to immense amount of work to be done in each step of the computation. These are the reasons why one has to think twice before using the stiff solver and to decide between the stiff and non-stiff solver. On the other hand a very interesting and promising numerical method of solving systems of ordinary differential equations based on Taylor series has appeared. The potential of the Taylor series has been exposed by many practical experiments and a way of detection and solution of large systems of ordinary differential equations has been found.


INTRODUCTION
This paper is related with computer simulations of continuous systems.The research group HPC ("High performance computing") [8] has been working on extremely exact and fast solutions of homogenous differential equations, nonlinear ordinary and partial differential equations, stiff systems, large systems of algebraic equations, real time simulations and corresponding software and hardware (parallel) implementations since 1980.
The "Modern Taylor Series Method" (MTSM) is used for numerical solution of differential equations.The MTSM is based on a recurrent calculation of the Taylor series terms for each time interval.Thus the complicated calculation of higher order derivatives (much criticised in the literature) need not be performed but rather the value of each Taylor series term is numerically calculated.
An important part of the MTSM is an automatic integration order setting, i.e. using as many Taylor series terms as the defined accuracy requires.Thus it is usual that the computation uses different numbers of Taylor series terms for different steps of constant length.
There are several papers that focus on computer implementations of the Taylor series method in different context "a variable order and variable step" (see, for instance, [1,3]).Another more detailed description of a variable step size version and software implementation of the Taylor series method can be seen in [11].The stability domain for several Taylor methods is presented in [2].Promising Astable combination of implicit Taylor series method with Trapezoidal rule is described in [9,10].

Modern Taylor Series Method
The best-known and most accurate method of calculating a new value of a numerical solution of ordinary differential equation y = f (t, y), y(0) = y 0 is to construct the Taylor series [6].
Methods of different orders can be used in a computation [12].For instance the order method (ORD = n) means that when computing the new value y i+1 uses n Taylor series terms in the form where h is integration step and DY are Taylor series terms.
Similarly we can construct implicit Taylor series method of order n in the form

STIFF SYSTEMS
Generally speaking, a stiff system contains several components, some of them are heavily suppressed while the rest remain almost unchanged.This feature forces the used method to choose an extremely small integration step and the progress of the computation may become very slow.However, we often need to find out the solution in a long range.It is clear that the mentioned facts are troublesome and ways to cope with such problems have to be devised.
One of the most frequently mentioned definition of stiff systems is to use stiffness ratio r [15].Let be a system of k ordinary differential equations.Let J be the Jacobian of the system (5) and λ i the eigenvalues of J .
The eigenvalues λ i are generally time dependent.Let the eigenvalues λ i be arranged in the following way: then the stiffness ratio is in the form The stiffness ratio r is a coefficient that helps to decide whether a problem is stiff or not.A higher r indicates a more stiff system.However, there is no exact value of the stiffness ratio r that would distinguish the non-stiff problems from the stiff-problems.For many problems in common practice the stiffness ratio r is "very high" (say 1 • 10 6 or higher).

Test example 1
Let us examine "school example" system with initial conditions y(0) = 1, z(0) = 1.Note: well known analytic solution of ( 8) is in the form Typically we calculate the Jacobian of the system (8) then we specify the eigenvalues of the system ( 8) We suppose that a > 0.0001 then the stiffness ratio of the system is in the form Many stiff systems solver needs to compute the Jacobian of the ODEs systems to detect the stiffness.Modern Taylor Series Method as implemented in TKSL software needn't compute Jacobian matrix or eigenvalues of the ODEs systems.
Explicit Taylor series solution of ( 8) is in the form similarly Let us analyze Taylor series terms in the first step.The absolute value of explicit Taylor series terms |DZn i | have rapidly decreasing trend.As we can see in Fig. 1 for constant a = 1 (respectively r = 10000) and integration step size h = 1, 15 Taylor series terms are needed to obtain result with absolute error EPS = 10 −10 .When the constant a is increased (the stiffness ratio r is also increased) we need to use more Taylor series terms to keep the stability of numerical method.In Fig. 2 resp.Fig. 3 we can see Taylor series terms for a = 10 (r = 100000) resp.a = 100 (r = 1000000).To obtain local error EPS = 10 −10 for a = 100 we need to use 295 Taylor series terms.
There is a problem when a = 100 (Fig. 3).As we can see in Fig. 3  As we can see in Fig. 3 and Fig. 2 absolute value of Taylor series terms |DY n i | have increasing character.Modern Taylor series method as implemented in TKSL automatically detects (from different and rapidly growing Taylor series terms) the stiffness in system (8) with growing constant a and automatically reduces integration step size h.Tendency of decreasing Taylor series terms after automatic decreasing step size is shown in Fig. 4. Fig. 4 Taylor series terms after automatic step size reduction

Conclusion:
The TKSL automatically detects stiff system (8) using Taylor series terms and automatically reduces integration step size until the strategy in Fig. 4 is obtained.After detection of stiffness (using explicit MTSM), implicit Taylor series method must be used as presented in the Test example 2 as follows.

Test example 2
Stiff systems in some literature [7] are defined as systems of ODEs where explicit numerical methods don't work and implicit numerical methods must be used.
Well known analytic solution of ( 15) is in the form The system (15) becomes stiff for b 0 and stiffness ratio is r = b.
Let's try to find the solution of ( 15) with explicit Taylor series method.Absolute error of numerical solution which is defined as difference between numerical z i and analytical z(t i ) solution where   TKSL automatically detects stiffness in system (15) (when b is growing) from Taylor series terms and TKSL uses automatically smaller step size.Implicit Taylor series method (as implemented in iTKSL software) has prosperous properties to solve stiff systems especially implicit Taylor series method has bigger absolute stability domain than those of explicit Taylor series method.
Let's solve the system (15) with implicit Taylor series method.Implicit Taylor series is in the form where higher derivations are in the form After substitution higher derivations into implicit Taylor series form (18) we obtain numerical solution in the form 1+h 2 b+hb+h , for implicit Taylor series order 1 (ORD = 1) that is implicit Euler method.
Similarly for ORD = 2 we obtain formula in the form Absolute error of numerical solution using implicit Taylor series method of different order is shown in Tab. 4. Note that constant b has no influence on error of computation.There is a problem with general formulation of y i+1 , z i+1 from implicit Taylor series formula (18) -other implicit numerical methods have the same problem.We have to use some iteration method to compute y i+1 , z i+1 in implicit form.Implicit Taylor series method with recurrent calculation of Taylor series terms and Newton iteration method (ITMRN) was implemented.Absolute errors of ITMRN of ORD = 1, 2, 3, 4 are the same as explicit calculations presented in Tab. 4 and only two Newton iterations are needed.Absolute errors of ITMRN ORD = 5, 6, 7 and number of Newton iterations j which are needed to obtain numerical results with tolerance T OL = 10 −10 are shown in Tab. 5. Problem with increasing the number of Newton iterations j with increasing order of implicit Taylor series method is shown in Tab. 5. We should use multiple arithmetic with growing constant a and order of ITMRN or we should reduce integration step size to obtain better stability of ITMRN (Tab.6).Absolute errors of ITMRN ORD = 8 and number of Newton iterations using in each step are shown in Tab. 6.There are two integration step size used: h 1 = 0.1 with absolute error |Error 1 (y)| and number of Newton iterations j 1 resp.h 2 = 0.05 with absolute error |Error 2 (y)| and number of Newton iterations j 2 .The same arithmetic (double precision) is used in both cases.

CONCLUSION
A very interesting and promising numerical method of solving systems of ordinary differential equations based on Taylor series has appeared.The question was how to harness the said "Modern Taylor Series Method" for solving of stiff systems.The potential of the Taylor series has been exposed by many practical experiments and a way of detection and solution of large systems of ordinary differential equations has been found.
The paper was presented during the INFORMATICS 2011 conference [21].

Fig. 2 10 nFig. 3
Fig.2Taylor series terms, a = 10 Multiple arithmetic: with growing constant a multiple arithmetic is needed.Only 9 Taylor series terms are used for integration step h = 0.1 and local error EPS = 10 −20 .Corresponding word length of Taylor series terms for large constant b are shown in Tab. 3.

Table 3
Multiple arithmetic