An Implicit Method for Numerical Solution of Singular and Stiff Initial Value Problems

An implicit method has been presented for solving singular initial value problems. The method is simple and gives more accurate solution than the implicit Eulermethod as well as the second order implicit Runge-Kutta (RK2) (i.e., implicitmidpoint rule)method for some particular singular problems. Diagonally implicit Runge-Kutta (DIRK) method is suitable for solving stiff problems. But, the derivation as well as utilization of this method is laborious. Sometimes it gives almost similar solution to the two-stage third order diagonally implicit Runge-Kutta (DIRK3) method and the five-stage fifth order diagonally implicit Runge-Kutta (DIRK5) method. The advantage of the present method is that it is used with less computational effort.


Introduction
Mathematical models of numerous applications from physics, chemistry, and mechanics take the form of systems of timedependent partial differential equations subject to initial or boundary conditions.For the investigation of stationary solutions, many of these models can be reduced to singular systems of ordinary differential equations, especially when symmetries problem in the geometry and polar, cylindrical, or spherical coordinates can be used.
The leading-edge model describing the avalanches dynamics [1] has the form of a singular initial value problem for a scalar ordinary differential equation.Koch et al. [2,3] applied implicit Euler method (backward) to evaluate the approximate solutions of this type of singular initial value problem and finally used an acceleration technique known as the Iterated Defect Correction (IDeC) to improve the approximations.The second order implicit Runge-Kutta (RK2) method (i.e., implicit midpoint rule) is a higherorder solver than the implicit Euler method for solving singular initial value problems.For some singular initial value problems, the present method gives more accurate solutions than those of implicit Euler method as well as second order implicit Runge-Kutta (RK2) method.
Implicit methods are more suitable than explicit methods for solving stiff problems because of their higher-order accuracy.Stiff differential equations arise in a variety of physical applications, such as network analysis and chemical or nuclear kinetics [4].Stiff differential equations also occur in many kinds of studies, such as biochemistry, biomedical system, weather prediction, mathematical biology, electronics, fluids and heat transfer [5].For example, the stiffness in heat transfer originates physically in one of two ways: sharp changes in the thermal environment or large differences in the rates at which components of the system can transfer heat.Alexander [6] solved stiff problems by using different stage and different order diagonally implicit Runge-Kutta (DIRK) methods.For these methods, the coefficient matrix has lower triangular structure with equal elements on the diagonal, sometimes referred to as singly equal diagonals.Ababneh and Ahmad [5] derived and applied a three-stage third order diagonally implicit Runge-Kutta (AM-DIRK3) method for solving stiff problems.Then, Ababneh et al. [7] evaluated stiff problems by using five-stage fifth order diagonally implicit Runge-Kutta (DIRK5) method.The diagonal elements of these methods (AM-DIRK3 and DIRK5) are equal.In this paper, we use two-stage third order diagonally implicit Runge-Kutta (DIRK3) method whose diagonal elements are not equal.But, the derivation as well as the utilization of the previous diagonally implicit Runge-Kutta (DIRK) methods is laborious.The aim of this paper is to derive a simple implicit method specially for solving singular and stiff initial value problems.
In [9], (2a) and (2b) were solved numerically for unequal interval as ℎ, 3ℎ, 2.3ℎ, 2 2 ⋅ 3ℎ, 2 3 ⋅ 3ℎ, . ... In this regard, step size as well as error gradually increased.To overcome this difficulty, the methods (2a) and (2b) is modified as where  +1 =   + ℎ.This method is useful for solving initial value problems having an initial singular point.Another method [8] can be derived for solving initial value problem (1) having an terminal singular point as Earlier, a two-stage third order diagonally implicit Runge-Kutta (DIRK3) method was used for solving various initial value problems.The formula is given as follow: where  1 =   +(ℎ/3)(  +ℎ/3,  1 ) and The disadvantage of this method is that it is calculated in two stages.On the contrary, (3) is calculated in a single stage.The order of truncation error of the formulae (3), (4), and ( 5) is O(ℎ 4 ).

Examples
In Section 1, it has already been mentioned that the first order singular initial value problem of the form ( 9) is used to represent the leading-edge model describing the avalanche dynamic [1].This model is used in the computation of the run-out length of dry-flowing avalanches.Moreover, the amplitude and phase equations appear in the first order nonlinear equations when Krylov-Bogoliubov-Mitropolski (KBM) method [10] is applied to solve the second order nonlinear differential equation ẍ +  + (, ẋ ) = 0. Recently, ẍ + 2 ẋ / +   = 0,  ̸ = 1 (well known as Lane-Emden equation) is solved by the KBM method.In this case, amplitude and phase equations [11] become first order singular differential equations of the form (9). The method is illustrated by following the first order singular and stiff initial value problems.
Example 1.Consider a first order initial value problem in the form Equation ( 6) can be written as Integrating (7) and applying initial condition, we get the exact solution of (6) which is To compare the present method to other existing methods we have derived the absolute error and plotted it in the figures.First, the absolute error of solution of ( 6) by the present method (i.e., (3)), the implicit Euler method, the second order implicit Runge-Kutta method, and the two-stage third order diagonally implicit Runge-Kutta (DIRK3) method are plotted in Figure 1 with step size ℎ = 0.001 when  = 1/2,  = −1, and  = 2.
Example 2. Consider another form of the first order singular initial value problem which is The exact solution of ( 9) is Journal of Computational Engineering Since the solution of ( 9) is singular at  = 0 for  < 0, so we choose  0 = 1 at  0 = 0.001 as the initial condition.The absolute errors of the solution of ( 9) by the four methods are plotted in Figure 2 with step size ℎ = 0.001 when  = −1/8 and  = 1.
Journal of Computational Engineering

Results and Discussion
A simple method has been presented to solve initial singular value and stiff problems.and 8).But, the DIRK3 and DIRK5 methods take more computational time than the present method.
From Figure 3 it is also observed that the method (i.e., (4)) gives more accurate solutions than those of (3) near the terminal singular point.The solution (3) as well as Euler's solution and DIRK3 solution is useless when singularity arises at terminal point (i.e., (4)) (see Figure 4).Even if these solutions are used near the singular point, the errors become significant.In this regard, solution (3) has been modified to (4).
Seeing all the above figures, we conclude that the present method (i.e., solution (3) or solution (4)) is more suitable than the implicit Euler and the second order implicit Runge-Kutta methods.However, the present method gives almost similar results to those obtained by DIRK3 and DIRK5 methods.The disadvantage of DIRK3 and DIRK5 methods is that they are laborious.

Figure 3 :Figure 4 :
Figure 3: The absolute error of the four methods with step size ℎ = 0.001.

Figure 5 :Figure 6 :
Figure 5: The absolute error of the three methods with step size ℎ = 0.01.

Figure 7 :Figure 8 :
Figure 7: The absolute error of the three methods with ℎ = 0.002.