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

Authors: M. Kamrul Hasan; M. Suzan Ahamed; M. S. Alam; M. Bellal Hossain
Journal: Journal of Computational Engineering
 (2013)

---

## Abstract

An implicit method has been presented for solving singular initial value problems. The method is simple and gives more accurate solution than the implicit Euler method as well as the second order implicit Runge-Kutta (RK2) (i.e., implicit midpoint 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.


---

## Body

## 1. Introduction

Mathematical models of numerous applications from physics, chemistry, and mechanics take the form of systems of time-dependent 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 higher-order 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.


## 2. Derivation of the Formula

Recently, an implicit formula [8, 9] has been derived for solving singular initial value problems
(1)y′(x)=f(x,y),y(x0)=y0
having an initial singular point in the form as(2a)y1=y0+3h4[3f(x0+h,y0)+f(x0+3h,y1)],(2b)yi=yi-2+3h4[3f(x3i-3,yi-1)+f(x3i,yi)]fori≥2,where, x3i=x3i-3+2(i-1)(3h),i≥2.In [9], (2a) and (2b) were solved numerically for unequal interval as h,3h,2.3h,22·3h,23·3h,…. In this regard, step size as well as error gradually increased. To overcome this difficulty, the methods (2a) and (2b) is modified as
(3)yi+1=yi+h4[3f(xi+h3,(yi+(yi+1-yi)3))hhhhhhhh+f(xi+1,yi+1)(xi+h3,(yi+(yi+1)3))],i=0,1,2,…,
where xi+1=xi+h.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
(4)yi+1=yi+h4[f(xi,yi)+3f(xi+2h3,(yi+2(yi+1-yi)3))],i=0,1,2,….
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:
(5)yi+1=yi+h4[3f(xi+h3,Y1)+f(xi+1,Y2)],hhhhhhhhhhhhhhhhhhhhhhi=0,1,2,…,
where Y1=yi+(h/3)f(xi+h/3,Y1) and Y2=yi+(h/2)[f(xi+h/3,Y1)+f(xi+1,Y2)].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(h4).


## 3. Examples

In Section1, 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  x¨+x+εf(x,x˙)=0. Recently,  x¨+2x˙/t+xp=0, p≠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(6)y′(x)=qyrxp,x>0,y(0)=1,0<p<1,-1≤q<0.
Equation (6) can be written as
(7)dyyr=qdxxp.
Integrating (7) and applying initial condition, we get the exact solution of (6) which is
(8)y=(x1-pq(-1+r)+(-1+p)(-1+p))1/(1-r),r≠1=e(qx1-p/(1-p)),r=1.
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 h=0.001 when p=1/2, q=-1, and r=2.Figure 1
The absolute error of the four methods with step sizeh=0.001.Example 2.
Consider another form of the first order singular initial value problem which is(9)y′(x)=qyrx,y(x0)=y0,q<0.
The exact solution of (9) is
(10)y=(y01-r+(q-qr)log(x)+q(-1+r)log(x0)y01-r)1/(1-r),r≠1=y0(xx0)q,r=1.
Since the solution of (9) is singular at x=0 for q<0, so we choose y0=1 at x0=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 h=0.001 when q=-1/8 and r=1.Figure 2
The absolute error of the four methods with step sizeh=0.001.Example 3.
Consider another first order initial value problem in the form as(11)y′(x)=qyr(1-x),0≤x<1,y(0)=13,0<p<1,q>0.
Clearly, (11) is singular at x=1. The exact solution of (11) with  p=1/2, q=1, and r=2 is 1/(1+21-x). The absolute error of the solution of (11) by the present methods (i.e., (3) and (4)), the second order implicit Runge-Kutta method, and the two-stage third order diagonally implicit Runge-Kutta (DIRK3) method are plotted in Figure 3 in the range 0.9≤x<1, and also the present method (i.e., (4)) and the second order implicit Runge-Kutta method are shown in Figure 4 in the range 0.998≤x≤1 with step size h=0.001 when  p=1/2, q=1, and r=2.Figure 3
The absolute error of the four methods with step sizeh=0.001.Figure 4
The absolute error of the two methods with step sizeh=0.001.Example 4.
Consider the stiff initial value problem [7]
(12)y′(x)=-200(y(x)-cos(x)),y(0)=0
with the exact solution y=(1/40001)(200(-200e-200x+200cos(x)+sin(x))), which is considered over the range 2≤x≤5, using step size h=0.01. The absolute error of the solution of (12) by the implicit Euler method, the second order implicit Runge-Kutta method, and the present method (i.e., (3)) are shown in Figure 5, and also the present (i.e., (3)), the DIRK3, and the DIRK5 methods are shown in Figure 6.Figure 5
The absolute error of the three methods with step sizeh=0.01.Figure 6
The absolute error of the three methods with step sizeh=0.01.Example 5.
Consider the highly stiff initial value problem [7]
(13)y′(x)=-1000y(x)+e-2x,y(0)=0
with the exact solution y=(1/998)(e-1000x(-1+e998x)), which is considered over the rang 0≤x≤1, using step size h=0.002. The absolute error of the solution of (13) by the implicit Euler method, the second order implicit Runge-Kutta method, and the present method (i.e., (3)) are shown in Figure 7, and also the present (i.e., (3)), the DIRK3, and the DIRK5 methods are shown in Figure 8.Figure 7
The absolute error of the three methods withh=0.002.Figure 8
The absolute error of the three methods withh=0.002.


## 4. Results and Discussion

A simple method has been presented to solve initial singular value and stiff problems. From the Figures1, 2, 3, 5, and 7 it is observed that the present method (i.e., (3)) gives better results than those of implicit Euler and the RK2 methods and almost similar results to DIRK3 and DIRK5 methods (Figures 1, 2, 3, 6, and 8). But, the DIRK3 and DIRK5 methods take more computational time than the present method.From Figure3 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.



---

*Source: 720812-2013-09-26.xml*