An Algorithm to Optimize the Calculation of the Fourth Order Runge-Kutta Method Applied to the Numerical Integration of Kinetics Coupled Differential Equations ()
1. Introduction
Adirovitch proposed a description of the luminescence in crystal phosphor through a kinetic process involving electrons trapped in a crystal defect, their subsequent thermally or optically induced promotion to a conduction band, and their recombination in another defect [1]. The use of this model allows the inclusion of physically meaningful quantities [2,3]. To date, the model has been applied and further extended to analyse, for example, to the optical absorption decay of Fe impurity oxidation [4], the thermal decay of the optical absorption of irradiation-induced color centers in spodumene [5,6], the luminescence of irradiated spodumene [7], the EPR of atomic hydrogen in glass and beryl [8-10], defects in diamond films [11], charge transport in thin-layer photodielectric systems [12], thermally stimulated processes in dielectrics [13], stimulated luminescence in insulators and semiconductors [14], soliton in biological systems [15], photo-induced-degradation of amorphous hydrogen silicon (a-Si:H) [16] and the degradation of a polymer [17].
In Adirovitch’s model, the concentration of trapped charges (
) and free charges in the conduction band (
) are related by the rate equations:
(1a)
(1b)
where
, 
and 
are positive adjustable parameters associated with the release of charges from traps, retrapping and the charge and electron-hole recombinetion process, respectively. In general the energy dissipation of the recombination process occurs through a radiation process. Thus, the function 
presented in Equation (1) is proportional to the luminescence in crystal phosphor and 
is the concentration of traps [1].
Previous research findings have shown that the solutions for the non-linear stability analysis of equations 1a and 1b are stable, following a hyperbolic path [18]. Mizukami et al. [19] observed that these solutions can be represented by the sum of two terms. The first term leads to slow first order decay and the other to a fast decay. Based on these observations they have proposed an approximate method to reach stable solutions. To obtain the exact solutions, iterative methods can be applied. A family of explicit and implicit iterative methods that has been frequently used are the Runge-Kutta methods [20-22].
In this paper we initially report numerical analyses of the fourth order Runge-Kutta method as applied to the solution of Adirovitch model Equations (1a) and (1b). Even though the RK method is stable, we identified a disconcerting property that emerges from the stiffness of the method when solving these equations: the numerical results are reliable and efficient only for a limited range of parameter values. To solve this problem, this paper proposes an algorithm to determine the best values for parameters in the Runge-Kutta method in order to guarantee the reliability of the results while using the smallest numbers of steps to reach a solution. This algorithm was implemented and has been used in different situations with different parameters. Thus far, to the best of our knowledge, the proposed algorithm has not produced any incorrect solutions.
The paper is organized as follows: in the first section we present the numerical analysis of the Adirovitch model; next, we present the problem of calculating the smallest number of steps to find the solution of this model using a Runge-Kutta Method; and finally we present an algorithm, together with its implementation and tests, to optimize the calculation of the smallest noncritical number of steps to reach the solution.
2. Numerical Analysis
To numerically solve the Equations (1a) and (1b) for the Adirovitch model presented in the introduction of this work it is necessary to find the correct values for each adjustable parameter
, 
, and
. Powerful genetic [23] and statistical [24] methods were developed recently and can be applied in fitting these parameters. Also, other classical methods for fitting can be adopted [25,26]. Nevertheless, fit methods may use a high number of iterations, and the number of steps used for the integration method may increase the computational cost. They may also require large adjustable parameter values. In this context, the solution of Equations (1) using the RungeKutta method, according to the parameters given, may require a large number of steps. Furthermore, the stiffness of the method can create difficulties in the numerical integration of Equations (1) resulting in abnormal and numerically unstable solutions. In this section we will examine the chaotic behaviours of these solutions. In the next section, we propose an algorithm to overcome them.
Equations (1) are made more tractable by introducing the transformations:
(2)
where 
is the total concentration of trapped charges at the initial instant
, and
. Since the parameters
, 
and 
enter linearly in the right hand side of Equations (1), mathematically only two of them are really independent. Therefore, it is natural to make an additional scaling. Then, using 
, we deal with only two parameters. Thus, Equations (1) are written as:
(3a)
(3b)
Numerical methods for the solution of Equations (3) can be divided into two categories [27]. One is the initial value problem, in which the values of 
are known at an initial time
, with
, requiring the evaluation of
. The other is the boundary condition for two points, in which boundary conditions are given for 
and
, typically.
The numerical solution of Equations (3) belonging to the first category are mostly grouped as Runge-Kutta [21, 22,27-30], Bulirsch-Stoer [31], predictor-corrector [32] solutions among others [33,34].
A widely used form of the Runge-Kutta method is of the fourth order. Using a vector notation typical in differential equations, 
and 
, the advancing formula 
from time 
to 
is given by:
(4)
where 
is the time increment for N steps in which the interval is divided,
The precision of 
is proportional to
. For orders greater than four, it is necessary to evaluate M to M + 2 more functions, increasing the computational cost. The equilibrium between computational efficiency and cost has been found for the fourth order Runge-Kutta method [13,23].
As the error is proportional to
, the Runge-Kutta method has a strong dependence on the number of steps. It predicts that the solution is invariant for large numbers of steps and that errors accumulate with any decrease in the number of steps.
Because of the dynamic nature of the Equations (3), it is possible that the accumulation of error makes the solutions erratic. Therefore, we investigated the types of errors that can happen depending on the number of steps and the values of the parameters of the Equations (3).
We have found at least four types of erratic solutions for Equations (3), which are shown in Figure 1. Here we assume that all charge traps are filled (
) and have no free charges at first (
). The time interval used for all calculations in the present work is
. Figure 1(a) shows an error in which the solutions behave as expected for the exact solutions but whose values are wrong. Figure 1(b) shows deformation errors in which the solution has anomalous depressions or elevations. Figure 1(c) shows an error in which solutions are duplicated, suggesting chaotic behavior. Finally, Figure 1(d) shows erratic solutions leading to divergence. Erratic solutions such as divergence, doubling and deformations are easily identified. Nevertheless, wrong solutions are more difficult to identify.
The existence of incorrect solutions can be evaluated by analysing the correction for the n-th step given by
. Figure 2 shows 
as a function of 
in the interval
, with
, initial parameters
, 
and (a)
, (b)
, (c)
, and (d)
. Solid lines linking the values obtained at each i-th step makes the visualization of the properties of 
easy.
In Figure 2(a) (
) it is possible to see that the modulus of differences decreases at each step and grows slowly with the increase of
. This behaviour changes starting at
, where oscillations can be seen. In Figure 2(b) (
) a stranger fin-like formation appears at approximately
, at the same time that the differences become erratic. In Figure 2(c) (
) the fin-like formation appears around 
and doubling appears in
. In the interval between the beginning of the fin-like formation and the beginning of the doubling, the solutions are incorrect. In Figure 2(d) (
), we observed three critical points, one at the fin-like formation (
), the second at the doubling (
), and the other at a point where the solution diverges completely (
). The differences show the appearance of incorrect fin-like, doubling, and divergence type solutions, but they do not allow for prediction
Figure 1. Solution y1 of Equations (3) using the Runge-Kutta method.
of the result of the sum of these differences. In fact we verified through an exact numerical solution of Equations (3) that wrong solutions for N = 170, 200, 230 and 270 begin at approximately
, respectively. Therefore, 
is not capable of revealing cumulative errors and a new algorithm must be considered in order to do so.
3. An Algorithm to Optimize the Calculation of the Number of Steps in the Runge-Kutta Method
The analysis given in Section 2 reveals that the Equations (3) behave chaotically depending on the number of steps. Previous research findings have tried to provide mathematical explanations for why such phenomena occur in nonlinear systems [35]. Our work does not intend to build upon these mathematical explanations. Instead, this work intends to create an algorithm that calculates the smallest number of steps necessary to solve the Equations (3), thereby preventing users/programs from encountering these chaotic behaviours. Along these lines, Hagebeuk and Kivits have presented an algorithm in the form of expansion of a very small parameter 
to overcome the problem of stiffness [35]. More recently, a new improved class of RK methods of the fourth order were applied with success in several sample problems [36-38]. Nevertheless, the problem of finding the smallest number of steps in which the solutions of RK methods are stable, numerically meaningful and efficiently calculated still remains a subject of great interest [39-42].
3.1. How to Determine the Smallest Non-Critical Number of Steps?
The Runge-Kutta method numerically solves Equations (3) if the solution in a given interval becomes independent of the number of steps (i.e. the number of steps is large enough). However, as discussed in Section 2, we find that 
can detect fin-like, doubling and divergence (Figure 2), but cannot show incorrect solutions. This observation is disturbing because it reminds us of the need to perform additional tests to ensure that the number of steps is large enough to give reliable results for the numerical integration. Therefore, we propose an algorithm that executes this task automatically as the solution to find the smallest number of steps to finding a correct solution to Equations (3).
In Section 3, we observed that incorrect solutions may occur from the beginning of the calculations. For example, considering the same interval of integration, for
, 
, which gives wrong results, 
in the first step is 0.954288... and for 
, 
, which gives correct results, it is 0.952203... after 10 steps (the interval is the same in both calculations). This result shows that there is a difference from the beginning of the calculation. On the other hand for N = 2000, p1 = 1, p2 = p3 = 65, the result of the calculation in the first step is 0.995014870675078 and for N = 20,000 after 10 steps the same level of numerical precision is achieved.
By examining these results, we conclude that the result of the application of the RK method using the interval of the first step and dividing it into 10 steps may reveal whether the number of steps will lead to incorrect results. So we will build the algorithm considering only the interval of the first step. We calculate the result of the application of the RK method with a single step in this interval and subtract the result of the calculation obtained by dividing this same interval into M steps. We write this algorithm as:
(5)
Let us consider the difference 
where 
and 
are calculated with one and 10 steps in the interval
, respectively. Then, since we are dealing with an expansion of the fourth order, we expect that to produce a correct solution this difference must be smaller than
. In case of incorrect solutions, the difference does not maintain proportional to 
and will increase.
Figure 3 shows 
as a function of the parameter
, with 
and the number steps (a) N = 270, (b)
, (c)
, and (d)
. In Figure 3(a) (
) the arrow shows a sharp increase in
, at approximately
, which coincides with the region where the the results are wrong. In Figure 3(b) (
) the arrow shows an abrupt increase in
, at about
, which coincides with the fin-like region. In Figure 3(c) (
), the first arrow at approximately
, shows the beginning of a plateau-like region that extends to
, where the fin-like region indicated by a second arrow begins. In Figure 3(d) (
), the first arrow shows the beginning of a plateau-like region extending to 
at about
, where the fin-like region indicated by a second arrow begins and the third arrow indicates divergence. The plateau-like region is not shown in the plot of 
(Figure 2).
Figure 4 shows the values of 
and 
as a function of the parameter
, with p2 = 65 and
. The features of 
and 
are very similar, until the second arrow; above the second arrow there are differences in the structure of the results. Therefore, the best description of the error in the application of the Runge-Kutta method in the solution of Equations (3) may include both 
and 
features.
To include the behavior of 
and
, as suggested by the analysis of Figure 4, we chose the last term of equation 1b, because it contains the variables y1 and y2. Therefore we write the new algorithm 
as:
(6)
Figure 5 shows 
as a function of
, with 
and (a)
, (b)
, (c)
, and (d)
. Figure 5(a) (
) shows agreement with the exact solution until 
, where a sharp increase begins. Figure 5(b)
(
) shows the exact solution region until
. Figure 5(c) (
) shows the plateau-like region between two arrows. Figure 5(d) (
) shows the plateau-like, fin-like and divergence regions successively.
The results of 
and 
are similar for small values of the number of steps (
and
). However when the number of steps increases closer to the number ideal of steps, we observed that 
reveals more details. For this reason, we adopted the algorithm 
to determine the smallest number of steps with which it is possible to reliably solve Equations (3) using the RK method.
Table 1 shows the number of steps and the values of 
at the end of the time interval determined using the algorithm
. For values of parameters 
, we see that the solution for
, is bifurcated for values equal to or greater than 0.004, but is convergent below 0.0039. We note that in the interval of values of delta between 0.0039 and E-06 the values of steps are equal and the values of 
match up to the third decimal place when compared to the more accurate results. The accuracy in the fifth decimal place is obtained in the interval of E-08 to E-10. For delta equal to E-08, the number of steps is 245 and for E-10 it is 340. For values of parameters
, we see that the solution is bifurcated for 
values equal to or greater than 0.02. Between 0.01 and 0.019, the solution is wrong. Between 0.004 and 0.001 accuracy is limited to the first decimal place, between 0.0001 and 0.00001 it extends to the third decimal place and for values less than 0.000001 to the fifth decimal place. For delta equal to E-06, the number of steps is 161 and 677 for E-10. In any case, we find that, even for very small values of delta, the magnitude of the number of steps is in an acceptable range of values.