Qualitative behavior of stiff ODEs through a stochastic approach

ABSTRACT


Introduction
Differential equations are used to model real-life systems by conserving their physical structures. There are different types of differential equations which have been named by according to their characteristics. Stiff differential equations are one of those. While developing a model of a system, it is necessary to consider suddenly occurred reactions with small time steps without neglecting that the system continue to behave over the whole-time interval. Stiff equations represent unstable behaviors for very small values. In other words, a model contains a point which decays or grows very rapidly than others. Despite natural restrictions of physical systems represented by stiff Ordinary Differential Equations (ODEs), they are commonly used in modelling various problems, through chemical reactions, while creating electrical circuits or studying in control theory etc. Not only modelling a stiff behavior but also solving the model accurately play a key role for capturing real-life behavior. Stiffness was firstly named by Curtiss and Hirschfelder [1] in 1952. Although this explanation leads to be realized that almost all real-life problems include stiff property, the first efficient algorithm for solving the model equations was suggested relatively late, in 1976 by Shampine and Gear [2]. Finding exact solution for stiff problems is generally limited to simple cases and conventional numerical methods have to be reconstructed with small time steps for these types of problems. However, the increased number of steps might possibly cause an accumulation of error. This fact gives rise to a necessity of alternative approaches for stiff equations. In the last few decades, various implicit and explicit methods related to stiffness have been developed. The explicit methods find a solution by using the current time information to produce later time information. However, implicit ones use the current and later time information at the same time. While analyzing stiff behavior, it should be taken into consideration how much small changes in the current time information affects the later time. Explicit methods generally do not work efficiently for catching the changing behavior in small step sizes or if they do, it converges very slowly than expected [3]. If the initial conditions cause a divergence in the solution, an explicit method requires impractically small step sizes to control the convergence. Although the implicit ones need more computation and requires sensitive implementations, they are properly applied to many stiff problems.
Even though an application area of numerical methods has a broad range, they are occasionally suffering from their restrictions. They may be seen to be efficient for the aim of the solving the problems iteratively, but these methods cannot be a first choice considering their computational time, computational error or effort spent for construction of a stiff structure. At this point, new approaches such that simulation techniques emerge by paying attention to these corresponding issues [4][5]. The Monte Carlo Method (MCM) is one of the basic simulation techniques [6][7][8]. It has been generally defined as a random sampling method for solving any model. Since this method uses basically random variables to represent the behavior of physical processes, it is classified as a stochastic approach. The MCMs can be applied to a wide range of problems in three different ways; sampling, estimation and optimization [9][10]. This classification depends on aim and a way of building algorithm. If a researcher wants to use simulation to mimic the nature of the system by creating objects or unreal systems, sampling methods are more useful than the rest. Therefore, random sampling and estimation techniques are used in this study to observe the behavior of the stiff differential equations.

Implementation of the method
The main intention of this study is to capture the exact behavior of stiff differential equations by using simulation techniques. To achieve this, differential equations are described by using integrals since the Monte Carlo integration is based on random sampling [11][12][13]. This randomness comes from uniformly distributed pseudorandom numbers selected by a sample space. The method is named by rejection sampling which is used for generating random variables with density function . The main advantage of using rejection sampling is that sampling can be used even if the density function cannot be integrated analytically. Let us then consider any first order differential equation in an implicit form: where function represent an arbitrary function with variables. After modifying the equation in this form, the algorithm needs a reference number which is chosen to do comparison in the related steps of the algorithm. The first reference number is generated by using initial conditions 0 and 0 and this number should be revised for each iteration. The step size is determined by dividing uniformly the interval to points. Let us call this reference number as Classification Number (CN) defined as follows where = 0,1, … , .
Next step, determination of upper and lower bounds for generating random numbers is expected to lead to more accurate estimation. The estimation can be made under the consideration of the physical realities of the problem. These bounds are determined by initial conditions. After determining an upper and a lower bound, random numbers can be created according to these bounds for making a comparison with the CN. To create random numbers, rand function of MATLAB can be used. This function generates different pseudorandom numbers between 0 and 1. They are known as pseudorandom since even if they act as a random number they are generated according to some artificial algorithm by the function. These random numbers between 0 and 1 are extended to the interval which determined by upper and lower bounds.
Then the comparison starts with created N positive random variables and N negative random variables by using rand function with respect to the CN of the algorithm. The way of implementation is given in the following pseudocode. Qualitative behavior of stiff ODEs through a stochastic approach 183 installed on a computer which has the properties of 2.3 GHz intel core i5 and 16 GB ram.

Example 1
Let us consider a first order stiff ODE with an initial condition, The exact solution of the differential equation is: The Monte Carlo based algorithm is applied to Equation (3) by dividing the time axis uniformly. Qualitative results including the solution produced by the proposed algorithm, the exact solution and the absolute errors have been illustrated in Figure 1. The corresponding numerical results can be seen in Table 1. It can be easily observed in Figure 2 that the stiffness occurs between the points 0 and 0.006, near to initial value. The initial deviation dampens fast due to the large value of coefficients. Despite the fact that each trial uses different set of random numbers to predict the results, each trial indicates the common feature at this stiff point. So the quantitative results are also close to each other. Even though quantitative results have slightly little deflections, qualitative results can be seen in good agreement with the exact results.  (3) Computational time of the proposed algorithm is 0.5548 s for this set of trial. Moreover, the accuracy is expected to be improved by reasonably decreasing the step size of the interval. However, the small step size leads a large number of comparisons in the algorithm, so that increasing the computational cost. Even if there are higher computational costs for some complex problems, it is seen that the proposed algorithm is the accurate solver as one of the simulation techniques.
The proposed algorithm is applied to Equation (5) by dividing the time axis uniformly. The comparison of predicted results with the analytical solution of Equation (5) is given Figure 3 and the corresponding absolute errors are shown in Figure 4. Quantitative results of Equation (5) are exhibited in Table 2.
Unlike Example 1, from Figure 3 it can be seen that the solutions already computed are much closer to the exact solution for all points in the range. It can be seen from the absolute errors that the results are very accurately predicted by the current algorithm. Computational cost of the present algorithm is 0.2913 s for this set of trial. Even if the computational time of the algorithm seems to be higher than its rivals depending on random numbers, its accuracy level is relatively in good agreement with available analytical solution.

Example 3
Now take a first order stiff differential equation system with initial conditions   (5) The Monte Carlo based algorithm is applied to Equation (7) by dividing the time axis uniformly. The comparison of the results of the current algorithm with the ode23s results and a closer view can be seen in Figures 5 and 6, respectively. The corresponding differences are illustrated in Figure 7. Quantitative results of Equation (7) are exhibited in Table 3.
As seen in the corresponding figures, the deviations originating near to the initial conditions have arisen rapidly. Two different equations behave separately; one is increasing while the other one is decreasing. Though the solution remains close to the referenced solution curves in a large scale of vertical axis, the deviations may occur. However, the qualitative and quantitative results can be seen in good agreement with the ode23s results. Even though ode23s is commonly accepted as one of the most suitable methods for properly capturing stiff behavior, the current method is seen to be as suitable as the ode23s. To support this, another suitable method, RK4 can be applied to this example. The difference between results of the simulation technique and the RK4 results are seen to be relatively small. Therefore, it has been claimed that the approach has ability to capture the stiff behavior. The predicted results have reasonable agreement with the results of ode23s and RK4 according to the Figure 8 and Table 3 and the differences between the predicted and ode23s solutions.
The computational costs of the current algorithm, ode23s and RK4 are 0.4214, 0.0180 and 0.1835 s, respectively. Despite the relatively higher computational cost of the proposed algorithm, the rest of its advantages is taken us to see attractiveness of the approach. In this respect, the computational cost can be sacrificed in simulation techniques in case of especially discrete and continuous methods have serious lack of accuracy or not existing solution for intricate problems.

Conclusions and recommendations
In this study, a Monte Carlo based stochastic algorithm has been developed to discover the behavior of realworld processes governed by stiff differential equations. All qualitative and quantitative results produced by the present algorithm have been seen to be in good agreement with the real environment. Despite the effect of randomness to error, the current procedure has been seen to produce highly acceptable results. Even if reconstructing the conventional methods with small time steps for stiff problems is affordable, simulation techniques can be better choices for challenging problems. In real-life problems, when there are sudden deviations in the consequences of random movements, it is necessary to consider the current stochastic approach that can handle the rapidly changing behavior.