# A Recent Development of Computer Methods for Solving Singularly Perturbed Boundary Value Problems

Authors: Manoj Kumar; Parul
Journal: International Journal of Differential Equations
 (2011)

---

## Abstract

This paper contains a surprisingly large amount of material and indeed can serve as an introduction to some of ideas and methods of singular perturbation theory. The work done in this area during the periods 1984–2000 and 2000–2005 has already been surveyed in 2002 and 2007 but our main objective is to produce a collection of important research articles of physical significance. In this paper, the crux of research articles published by numerous researchers during 2006–2010 in referred journals has been presented, and this leads to conclusions and recommendations about what methods to use on singular perturbation problems.


---

## Body

## 1. Introduction

As science and technology develop, many practical problems, such as the mathematical boundary layer theory or approximation of solution of various problems described by differential equations involving large or small parameters, have become increasingly complex and therefore require the use of asymptotic methods. However, the theory of asymptotic analysis for differential operators has mainly been developed for regular perturbations, where the perturbations are subordinate to the unperturbed operator. In some problems the perturbations occur over a very narrow region across which the dependent variable undergoes very rapid changes. These narrow regions are frequently adjacent to the boundaries of the domain of interest because a small parameter multiplies the highest derivative. Consequently, they are usually referred to as boundary layers in fluid mechanics, edge layers in solid mechanics, skin layers in electrical applications, shock layers in fluid, and solid mechanics, transition points in quantum mechanics.The class of boundary value problems which we will consider in this paper has for many years been a source of difficulty and entertainment to applied mathematicians and numerical analysts alike. Mathematically, we may consider a system of ordinary differential equations where one or more of the highest derivatives appearing are multiplied by a small parameterε. If we let ε approach 0, the order of the system reduces, and one cannot in general expect all boundary conditions imposed on the original ordinary differential equation system to be satisfiable. The perturbation is then “singular.” When ε is not 0 but is small, the solution is expected, under certain conditions, to exhibit narrow regions of very fast variation (so-called boundary or interior layers) which connect wider regions where it varies more slowly.Alternatively, we can define a singular perturbation problem as one whose perturbation series either does not take the form of a power series or, if it does, the power series has a vanishing radius of convergence. In singular perturbation theory there is sometimes no solution to the unperturbed problem, that is, the exact solution as a function ofε may cease to exist when ε=0, when a solution to the unperturbed problem does exist, its qualitative features are distinctly different from those of the exact solution for arbitrarily small but nonzero ε. In either case, the exact solution for ε=0 is fundamentally different in character from the “neighboring” solutions obtained in the limit ε→0. If there is no such abrupt change in character, then we would have to classify the problem as a regular perturbation problem.When dealing with a singular perturbation problem, one must take care to distinguish between the zeroth-order solution (the leading term in the perturbation series) and the solution of the unperturbed problem, since the latter may not even exist. There is no difference between these two in a regular perturbation theory, but in a singular perturbation theory the zeroth-order solution may depend onε and may exist only for nonzero ε.Singular perturbations were first described by Prandtl [1] in a seven-page report presented at the Third International Congress of Mathematicians in Heidelberg in 1904. However, the term “Singular Perturbations” was first used by Friedrichs and Wasow [2] in a paper presented at a seminar on nonlinear vibrations at New York University. The solution of singular perturbation problems typically contains layers. Prandtl originally introduced the term “boundary layer,” but this term came into more general use following the work of Wasow [3]. Pearson [4] was among the first to consider the finite-difference method to solve singular perturbation problems. In recent years researchers have used cubic splines and finite-element methods to solve singular perturbation problems, and a large number of papers and books have been published describing the various methods which have been used [5–28].The article [29] features a survey of some recent developments in asymptotic techniques, which are valid for weak nonlinear equations as well as stronger one throughout the whole solution domain. This review article [29] also emphasizes on perturbation method’s limitations that it always assumes small parameter in an equation although most of strong nonlinear problems have no small parameters at all and it is very tricky to choose small parameter which affects the nature of solution. Selected small parameter may be valid only for small values, In [29] author described some new categories of the asymptotic methods, which are variational approaches, parameter-expanding method, parameterized perturbation method, homotopy perturbation method, iteration perturbation method, and Ancient Chinese methods.In 2008, He published another review article [30] describing an elementary introduction to the concept of the recently developed asymptotic methods and new developments. In [30] the author paid particular attention throughout the paper to provide an intuitive grasp for Lagrange multiplier, calculus of variations, optimization, variational iteration method, parameter-expansion method, Exp-function method, homotopy perturbation method and Ancient Chinese mathematics as well. Subsequently, nanomechanics in textile engineering and E-infinity theory in high energy physics, Kleiber’s 3/4 law in biology, possible mechanism in Spider-spinning process, and fractal approach to carbon nanotube are briefly introduced. Bubble-electro spinning for mass production of nanofibers is illustrated.In 2009, Pamuk [31] presents a review of some recent results for the approximate analytical solutions of nonlinear differential equations. To do this, he introduced and compared three methods, namely, Adomian’s decomposition method (ADM), He’s Homotopy perturbation method (HPM), and Variational iteration method (VIM), which are recently studied by researchers in various fields of science and engineering. Finally the author of [31] concluded that numerical results are almost the same; HPM is much easier and more convenient and efficient than ADM and VIM.The area of singular perturbations is a field of increasing interest to applied mathematicians. During the last few years much progress has been made in the theory and in the computer implementation of the numerical treatment of singular perturbation problems. The development of small parameter methods led to the efficient use of boundary layer theory in various fields of applied mathematics, for instance, fluid mechanics, fluid dynamics, elasticity, quantum mechanics, plasticity, chemical-reactor theory, aerodynamics, plasma dynamics, magneto-hydrodynamics, rarefied-gas dynamics, oceanography, meteorology, diffraction theory, reaction-diffusion processes, nonequilibrium and radiating flows and other domains of the great world of fluid motion. In this paper, we give some singular perturbation models which arise in some of the above-mentioned areas. We will omit the details of techniques used to solve these models. Interested readers can see the references cited along with these models, for more details. The present paper deals with the research articles published during 2006–2010 on both linear and nonlinear type of singular perturbation problems.After the introduction and getting an overview of the subject, we illustrate some standard singular perturbation models in Section2 and the results of some numerical techniques used to solve linear singular perturbation boundary value problems are presented in Section 3, Section 4 contains some numerical techniques used to solve nonlinear singular perturbation boundary value problems, and Section 5 reports the conclusion and direction for further developments.


## 2. Some Standard Singular Perturbation Models

In this section, some popular models of singular perturbation problems arising in engineering problems with illustrative figures have been briefly described, and for details we refer to [17, 32].


### 2.1. A Piston Problem

We consider the one-dimensional flow of a gas in a long, open-ended tube. The gas is brought into motion by the action of a piston at one end, which moves forward at a speed which is much less than the speed of sound in the gas as shown in the Figure1. This is usually called the acoustic problem. The gas is modelled by the isentropic law for a perfect gas Pressure∝(Density)γ, 1<γ<2, and is described by the equations(2.1)at+uax+12(γ-1)aux=0,ut+uux+2γ-1aax=0,
where a is the (local) sound speed in the gas and u its speed along the tube. The initial and boundary conditions are u=0 and a=1 at t=0, x>εxp(t), and u=εV(t) on x=εxp(t), t≥0, with xp(t)=∫0tV(t′)dt′ and V(0)=0.Figure 1
Sketch of a piston moving a gas (according tox=εxp(t)) in an open tube.


### 2.2. Child’s Swing

We are all familiar with the child’s swing and the technique for increasing the arc (amplitude) of the swing. The process of swinging the legs (coupled with a small movement of the torso) causes the centre of gravity of the body to be raised and lowered periodically. This can be modeled by treating the swing as a pendulum which changes its lengthl(t) by a small amount. The model equation for this in the absence of damping is(2.2)d2θdt2+2ldldtdθdt+sin⁡θ=0,t≥0
for θ(t), the angle of the swing, given l(t) as in Figure 2. We choose to represent the child’s movement on the swing by l=l0(1+εsinωt), where l0 is a positive constant and ω is a constant frequency to be selected. Thus we approximate (2.2) as
(2.3)d2θdt2+(2εωcos⁡wt)(1+εsin⁡ωt)dθdt+θ=0.Figure 2
Pendulum of variable lengthl(t), swinging through the angle θ(t) (in a vertical plane).


### 2.3. Very Viscous Flow Past a Sphere

Here, we considerRe→0+, since Re=Ud/v, where U is typical speed of the flow, d is typical dimension of the object in the flow and v is the kinematic viscosity. The case of axisymmetric flow produced by a uniform flow at infinity parallel to the chosen axis, past a solid sphere, is illustrated in Figure 3.Figure 3
Coordinates and velocity components for the uniform flow past a sphere.It is convenient to introduce a stream functionΨ usually called a Stokes stream function, eliminate pressure from the Navier-Stokes equation and hence work with the non-dimensional equation:(2.4)1r2sin⁡θ{ψθ∂∂r-ψθ∂∂θ+2(cotθ)ψr-2rψθ}D2ψ=1ReD4ψ,1≤r<∞,0≤θ≤π,
where
(2.5)D2≡∂2∂r2+1r2∂2∂θ2-cotθr2∂∂θwithψ=ψr=0on⁡r=1,ψr2⟶12sin2θasr⟶∞;
this latter condition ensures that there is an axisymmetric flow of speed one at infinity. The velocity components are
(2.6)u=1r2sin⁡θψθ,v=-1rsin⁡θψr.


### 2.4. A Variable-Depth Korteweg-de Vries Equation for Water Waves

In fluid mechanics, we consider the one-dimensional propagation of waves over water (incompressible), which is modelled by an inviscid fluid without surface tension. The water is stationary in the absence of waves, but the local depth varies on the same scale(ε) that is used to measure the weak nonlinearity and dispersive effects in the governing equations. For right-running waves, the appropriate far-field coordinates are ξ=ε-1χ(X)-t and X=εx where χ(X) is to be determined.The nondimensional equations are(2.7)-uξ+ε[u(χ′uξ+εuX)+vuz]=-(χ′pξ+εpX),ε[-vξ+ε[u(χ′vξ+εvX)+vvz]=-pz,χ′uξ+εuX+vz=0
with
(2.8)p=ηandv=-ηξ+u(χ′nξ+εnX)on⁡z=1+ε,v=εuB′(X)on⁡z=B(X),
where (u,v) are the velocity components of the flow, p is the pressure in the fluid relative to the hydrostatic pressure, and z=1+εη is the surface of the water. The bottom is represented by the function z=B(X) (see Figure 4).Figure 4
Wave propagation in stationary water over variable depth.


## 2.1. A Piston Problem

We consider the one-dimensional flow of a gas in a long, open-ended tube. The gas is brought into motion by the action of a piston at one end, which moves forward at a speed which is much less than the speed of sound in the gas as shown in the Figure1. This is usually called the acoustic problem. The gas is modelled by the isentropic law for a perfect gas Pressure∝(Density)γ, 1<γ<2, and is described by the equations(2.1)at+uax+12(γ-1)aux=0,ut+uux+2γ-1aax=0,
where a is the (local) sound speed in the gas and u its speed along the tube. The initial and boundary conditions are u=0 and a=1 at t=0, x>εxp(t), and u=εV(t) on x=εxp(t), t≥0, with xp(t)=∫0tV(t′)dt′ and V(0)=0.Figure 1
Sketch of a piston moving a gas (according tox=εxp(t)) in an open tube.


## 2.2. Child’s Swing

We are all familiar with the child’s swing and the technique for increasing the arc (amplitude) of the swing. The process of swinging the legs (coupled with a small movement of the torso) causes the centre of gravity of the body to be raised and lowered periodically. This can be modeled by treating the swing as a pendulum which changes its lengthl(t) by a small amount. The model equation for this in the absence of damping is(2.2)d2θdt2+2ldldtdθdt+sin⁡θ=0,t≥0
for θ(t), the angle of the swing, given l(t) as in Figure 2. We choose to represent the child’s movement on the swing by l=l0(1+εsinωt), where l0 is a positive constant and ω is a constant frequency to be selected. Thus we approximate (2.2) as
(2.3)d2θdt2+(2εωcos⁡wt)(1+εsin⁡ωt)dθdt+θ=0.Figure 2
Pendulum of variable lengthl(t), swinging through the angle θ(t) (in a vertical plane).


## 2.3. Very Viscous Flow Past a Sphere

Here, we considerRe→0+, since Re=Ud/v, where U is typical speed of the flow, d is typical dimension of the object in the flow and v is the kinematic viscosity. The case of axisymmetric flow produced by a uniform flow at infinity parallel to the chosen axis, past a solid sphere, is illustrated in Figure 3.Figure 3
Coordinates and velocity components for the uniform flow past a sphere.It is convenient to introduce a stream functionΨ usually called a Stokes stream function, eliminate pressure from the Navier-Stokes equation and hence work with the non-dimensional equation:(2.4)1r2sin⁡θ{ψθ∂∂r-ψθ∂∂θ+2(cotθ)ψr-2rψθ}D2ψ=1ReD4ψ,1≤r<∞,0≤θ≤π,
where
(2.5)D2≡∂2∂r2+1r2∂2∂θ2-cotθr2∂∂θwithψ=ψr=0on⁡r=1,ψr2⟶12sin2θasr⟶∞;
this latter condition ensures that there is an axisymmetric flow of speed one at infinity. The velocity components are
(2.6)u=1r2sin⁡θψθ,v=-1rsin⁡θψr.


## 2.4. A Variable-Depth Korteweg-de Vries Equation for Water Waves

In fluid mechanics, we consider the one-dimensional propagation of waves over water (incompressible), which is modelled by an inviscid fluid without surface tension. The water is stationary in the absence of waves, but the local depth varies on the same scale(ε) that is used to measure the weak nonlinearity and dispersive effects in the governing equations. For right-running waves, the appropriate far-field coordinates are ξ=ε-1χ(X)-t and X=εx where χ(X) is to be determined.The nondimensional equations are(2.7)-uξ+ε[u(χ′uξ+εuX)+vuz]=-(χ′pξ+εpX),ε[-vξ+ε[u(χ′vξ+εvX)+vvz]=-pz,χ′uξ+εuX+vz=0
with
(2.8)p=ηandv=-ηξ+u(χ′nξ+εnX)on⁡z=1+ε,v=εuB′(X)on⁡z=B(X),
where (u,v) are the velocity components of the flow, p is the pressure in the fluid relative to the hydrostatic pressure, and z=1+εη is the surface of the water. The bottom is represented by the function z=B(X) (see Figure 4).Figure 4
Wave propagation in stationary water over variable depth.


## 3. Numerical Techniques to Solve Linear Singular Perturbation Boundary Value Problems

In this section, we summarized the research articles published in the last five years of singular linear perturbation problems. Most of the articles collected here involve the equations which arise in physical problems. To maintain a specific chronological order, first we consider self-adjoint differential equations and then differential equation of physical nature.In [33] Kadalbajoo and Arora developed a B-spline collocation method, which leads to a tridiagonal linear system for solving singularly perturbed equations and used the artificial viscosity to capture the exponential features of the exact solution on a uniform mesh. The design of artificial viscosity parameter is confirmed to be a crucial ingredient for simulating the solution of the problem. In this method, they replace the perturbation parameter ε affecting the highest derivative by artificial viscosity. The artificial viscosity is determined using the truncation error of the corresponding scheme which is zero for the boundary layer with constant coefficients. This scheme is shown to be second-order uniformly convergent, and it is also unconditionally stable. B-splines are used because they yield results of higher accuracy as compared with those of polynomial interpolation. In [33] the following class of self-adjoint singularly perturbed two-point boundary value problem is considered:(3.1)L1u≡-ε2u′′(x)+b(x)u(x)=f(x)for0≤x≤1
with boundary conditions(3.2)u(0)=ρ1,u(1)=ρ2,
where ε is a small parameter and f(x),b(x) are smooth functions and satisfy b(x)≥b*>0, for all x∈[0,1] for some constant b*. The boundary value problem (3.1) and (3.2) under these assumptions possesses unique solution u(x). In general, as ε tends to zero, the solution u(x) may exhibit exponential boundary layers at both ends of the interval [0,1]. The computed solution without using viscosity oscillates in the boundary layer regions for smaller ε. To control these disturbances, we can use artificial viscosity, and the results are far better than those without using viscosity.In [34] presented the analysis of an upwind scheme for obtaining the global solution and the normalized flux for a convection diffusion singularly perturbed two-point boundary value problem of the form(3.3)Lu(x)≡-εu′′(x)-a(x)u′(x)+b(x)u(x)=f(x),x∈Ω=(0,1),u(0)=s0,u(1)=s1,
where 0<ε≪1 is a small singular perturbation parameter, the functions a(x), b(x), f(x) are sufficiently smooth, and s0,s1 are given constants. Further, we assume that 2α*≥a(x)>2α>0 and b(x)≥0. Under these assumptions, the above problem (3.3) has a unique solution which exhibits a boundary layer at x=0.The upwind finite difference scheme is applied on a nonuniform mesh formed by equidistributing the arc-length monitor functionM(x)=1+(u′(x))2, which is bounded below by unity. The optimal choice of the monitor function depends on the problem being solved, the numerical discretization being used, and the norm of the error that is to be minimized. The monitor function is often based on a simple function of the derivatives of the unknown solution. It is shown that the discrete solution obtained by the upwind scheme and the global solution obtained via interpolation converges uniformly with respect to the perturbation parameter. Authors proved the uniform first-order convergence of the weighted derivative of the numerical solution on the nonuniform mesh and uniform convergence of the global normalized flux on the whole domain.In [35] authors concluded that the exponential B-spline collocation method on uniform mesh is a satisfactory way for solving the self-adjoint singularly perturbed Dirichlet boundary value problem and method shows second-order uniform convergence. In case of cubic B-spline collocation method as we decrease the value of ε to study the behavior of the solution at the boundary layer, and we need large number of nodal points in that region. For this purpose, Kadalbajoo and Aggarwal [9] considered the fitted mesh approach that consists of more number of nodal points in the boundary layer region. In [35] authors used the uniform mesh that consists of less number of nodal points in the boundary layer region in comparison to the fitted mesh approach described in [9]. They used the exponential B-splines basis that leads to the tridiagonal system which can be solved by the well-known algorithm. This scheme with exponential B-splines on uniform mesh gives more accurate approximations than the scheme using cubic B-splines on fitted mesh as well as the scheme using cubic B-splines on uniform mesh. The essential idea in this method is to use the B-splines basis for the space of exponential spline to approximate the solution of given problems via nodal collocation approach. It is relatively simple to collocate the boundary value problem at the nodal points, to set up the collocation system, and to solve them. This method produces a spline function which is useful to obtain the solution at any point of the interval.Research article [36] deals with tension spline numerical methods developed on a nonuniform mesh for the efficient numerical integration of singularly perturbed two-point singular boundary value problems. This method is valid for problems in rectangular as well as in polar coordinate system, includes the error analysis. Here, they have shown that, by making use of the continuity of the first derivative of the tension spline function, the resulting tension spline difference scheme gives a tridiagonal system which can be solved efficiently by using a tri-diagonal solver. Mohanty and Arora [36] also discussed the application of proposed methods to singular problem and numerical results show that the proposed tension spline methods produce oscillation-free solution for 0<ε≪1 everywhere in the solution region 0<x<1. Test problems have been solved to demonstrate the efficiency of the proposed method when ε is either small or large as compared to the corresponding variable mesh size. The main idea of the proposed nonuniform mesh technique is to consider the geometric location of the grid points, which depends on the choice of mesh ratio parameter σ. This technique may be extended to derive other numerical methods, not necessarily limited to tension spline methods.In [37] a robust computational method is presented for the coupled system of singularly perturbed reaction-diffusion boundary value problems (BVPs). First, they construct a difference scheme using cubic spline on nonuniform mesh for system of reaction-diffusion BVPs. Then, they apply this scheme on a layer resolving piecewise uniform mesh, known as the Shishkin mesh. The Hybrid numerical scheme is applied on a piecewise-uniform Shishkin mesh. In the boundary layer (inner) region the mesh is fine and the cubic spline scheme is stable, whereas in the regular (outer) region the mesh is coarse and the cubic spline scheme is not stable. This newly proposed scheme is uniformly stable throughout the domain and provides second-order uniform convergence results. Here, they deal with the case ε1=ε and ε2=1. Stability analysis, truncation errors, and parameter-uniform error estimates have been derived, which shows that the proposed method is of almost second order up to a logarithmic factor. A numerical experiment is carried out to verify the mathematical results. The maximum pointwise errors and rate of convergence reflect the theoretical results and show the accuracy and efficiency of the method. Nonlinear systems of equations have been handled by this method after linearization.In [38] Rao and Kumar present a higher-order cubic B-spline collocation method for the numerical solution of self-adjoint singularly perturbed boundary value problems that are much easier and more efficient for computing. The new idea in this method is to divide the domain of the differential equation into three nonoverlapping subdomains and solve the regular problems obtained by transforming the differential equation with respective boundary conditions on these sub domains using the present higher-order B-spline collocation method. The boundary conditions at the transition points are obtained using the zeroth-order asymptotic approximation to the solution of the problem. Authors consider the two-step spline approximate method because a special type of tridiagonal system is obtained and method is more efficient than classical finite difference scheme on piecewise uniform Shishkin meshes and reasonably comparable with the quintic spline difference scheme. Quintic splines are considered to obtain fourth order convergence of the scheme. The method [38] gives optimal order convergence with cubic B-splines only.In [39] authors considered a class of singularly perturbed two-point singular boundary value problem of the form (3.4)-εu′′(x)+kxu′(x)+g(x)u(x)=f(x),0≤x≤1,u(0)=γ0,u(1)=γ1,
where 0<ε≤1, f(x), g(x) are bounded continuous functions in (0,1), g(x)>0, and γ0, γ1 are finite constants and developed a numerical technique for a class of singularly perturbed two point singular boundary value problems on an uniform mesh using polynomial cubic spline.To derive the method, the authors used a polynomial functionSΔ(x) of class C2[0,1] in the following form:(3.5)SΔ(x)=ai+bi(x-xi)+ci(x-xi)2+di(x-xi)3,i=0,1,2,…,N-1,Δ:0=x0<x1<x2<⋯<xN-1<xN=1,xi=ih,i=0,1,2,…,N,
where h=1/N and ai,bi,ci, and di are constants.From algebraic manipulation the following expressions are obtained:(3.6)ai=yi,bi=yi+1-yih-h(Mi+1-Mi)6,ci=Mi2,di=Mi+1-Mi6h,i=0,1,2,…,N-1,
where SΔ(xi)=yi,SΔ(xi+1)=yi+1,SΔ′′(xi)=Mi,SΔ′′(xi+1)=Mi+1, and yi is an approximate solution to y(xi), obtained by the splines function SΔ(x) passing through the points (xi,yi) and (xi+1,yi+1).Then, at the grid pointsxi, the proposed equation (3.4) is discretized in the form(3.7)-εui′′+piui′+giui=fiwherepi=kxi,gi=g(xi),fi=f(xi).
Finally, a tridiagonal system is obtained, which gives the approximations u1,u2,…,uN-1 of the solution u(x) at x1,x2,…,xN-1:(3.8)[-ε-14hpi-1+112hpi+1-13hpi+16h2gi-1]ui-1+[2ε+13hpi-1-13hpi+1+23h2gi]ui+[-ε-112hpi-1+14hpi+1+13hpi+16h2gi+1]ui+1=h26(fi-1+4fi+fi+1),withu(0)=γ0,u(1)=γ1,i=2(1)N-1.
As in (3.8) the coefficients of ui, ui-1, and ui+1 become infinite for i=1, so using L’Hospital rule and (3.8) we obtain the following boundary formula: (3.9)[2ε-13hp2+23h2g1]u1+[-ε+13hp1+14hp2+16h2g2]u2=h2(23f1+16f2)-ε6h2(f0-g0γ0k-ε)+(ε+13hp1+112hp2)γ0.
Also, in [39] authors demonstrated this scheme on three numerical examples and presented the root mean square (RMS) errors at the nodal points and compared results with [36, 40]. The scheme derived is of second order and accurate and resulting linear system of equations has been solved by using a tridiagonal solver.In [41] authors design nonstandard finite difference schemes for the following self-adjoint singularly perturbed two-point boundary value problems: (3.10)Ly≡-ε(a(x)y′)′+b(x)y=f(x)on⁡(0,1),y(0)=η0,y(1)=η1,
where η0, η1 are given constants and ε is a small positive parameter. Further, f(x), a(x), and b(x) are sufficiently smooth functions satisfying the conditions a(x)≥a>0 and b(x)≥b>0. The main concern with such problems is the rapid growth or decay of the solution in one or more narrow layer regions. The specific problem under consideration in [41] is called dissipative because the rapidly varying component of the solution decays exponentially (dissipates) away from a localized breakdown or discontinuity point in the layer regions as ε→0. There are essentially two strategies to design schemes which have small truncation errors inside the layer regions. The first approach, which forms the class of fitted mesh methods, consists in choosing a fine mesh in the layer regions. The second approach is in the context of the fitted operator methods in which the mesh remains uniform and the difference schemes reflect the qualitative behavior of the solutions inside the layer regions. The article [41] falls under the second category of strategies and used the Mickens nonstandard finite difference method [42]. Essential physical properties (e.g., dissipativity) of the solutions of such problems are captured in the schemes by an appropriate renormalization of the denominator of the discrete derivative. The proposed method has second-order ε-uniform convergence.In [43] Bawa and Natesan considered the following singularly perturbed self-adjoint boundary value problem:(3.11)Lu(x)≡-εu′′(x)+b(x)u(x)=f(x),x∈D=(0,1)u(0)=A,u(1)=B,
where 0<ε≤1 is a small parameter, b and f are sufficiently smooth functions, such that b(x)≥β>0 on D¯=[0,1]. Under these assumptions, the boundary value problem (3.11) possesses a unique solution (3.12)u(x)∈C2(D)∩C(D̅).
In general, the solution u(x) may exhibit two boundary layers of exponential type at both end points x=0 and x=1. In this article, a parallelizable computational technique for singularly perturbed reaction diffusion problems is analyzed and implemented on parallel computer. In this technique, the domain is decomposed into nonoverlapping subdomains and boundary value problems are posed on each subdomain with suitable boundary conditions. Then, each problem is solved by the adaptive spline-based difference scheme on each subinterval on parallel computer. To analyze and implement this parallel computational technique for singularly perturbed boundary value problems of the form (3.11), firstly, a nonoverlapping domain decomposition of the computational domain is devised, and then independent BVPs on each subdomain are proposed. Suitable boundary conditions have been supplied from the asymptotic approximate solution, which are not difficult to calculate either theoretically or numerically. Stability and truncation error have been discussed, and error estimate has been obtained. It is observed that the implementation of proposed scheme reduces the maximum absolute error and the CPU time to much extent.In [44] authors propose a tailored finite point method for a kind of singular perturbation problems in unbounded domains. First, they introduce two artificial boundaries and use the artificial boundary method [45] to reduce the original problem to a problem on bounded computational domain equivalently or approximately. Then they propose new approach to construct a discrete scheme for the reduced problem, where the finite point method has been tailored to some particular properties or solutions of the problem. The finite point method is a development of finite difference method, in which the meshless technique is emphasized. For one-dimensional singular perturbation problem, the method is very close to the method of “exponential fitting.” The traditional finite element method or finite difference method cannot get the satisfactory numerical results with the same mesh, especially for small ε. Numerical results show that new methods can achieve very high accuracy with very coarse mesh whenever ε is small or large. In contrast, the traditional finite element method or finite difference method cannot get the satisfactory numerical results with the same mesh, especially for small ε.In [46] authors present a numerical method to solve singularly perturbed two point boundary value problems for second-order ordinary differential equations with a discontinuous source term [47] considered on the unit interval Ω=(0,1). A single discontinuity is assumed to occur at a point d∈Ω. An asymptotic initial value method to solve the following class of problems is described:(3.13)Lu(x)≡-εu′′(x)+a(x)u(x)=f(x),x∈Ω-∪Ω+,u(0)=p,u(1)=q,
where ε is a small positive parameter, a(x) are sufficiently smooth functions on Ω̅ and a(x)≥α>0, f(x) is sufficiently smooth function in Ω-∪Ω+∪{0,1}, f and its derivatives have jump discontinuity at d. Because f is discontinuous at d the solution u of (3.13) does not necessarily have a continuous second derivative at the point d. Thus, u need not belong to the class of functions C2(Ω). But the first derivative of the solution exists and is continuous. Under these assumptions, the singularly perturbed problems (3.13) have a solution (3.14)u∈U≡C0(Ω̅)∩C1(Ω)∩C2(Ω-∪Ω+).
Firstly asymptotic expansion approximation of the solution of the boundary value problem is constructed using the WKB perturbation method. Then, some initial value problems and terminal value problems are constructed such that their solutions are the terms of this asymptotic expansion. These initial value problems, and terminal value problems are singularly perturbed problems and therefore fitted mesh method (Shishkin mesh) are used to solve these problems. Since the asymptotic expansion approximation is used and the boundary value problem is converted into initial value problems and terminal value problems, therefore the method is called Asymptotic initial value method (AIVM). In general solving an initial value problem is much easier than solving a boundary value problem.In [48] the Robust eigenvalue clustering in a specified circular region problem that is the robust D-stability problem of linear discrete singular time-delay systems with structured (elemental) parameter uncertainties is investigated. Under the assumptions that the linear nominal discrete singular time-delay system is regular and causal and has all its finite eigenvalues lying inside a specified circular region, a new sufficient condition is proposed to preserve the assumed properties when the structured parameter uncertainties are added into the linear nominal discrete singular time-delay system. When all the finite eigenvalues are just required to locate inside the unit circle of the z-plane, the proposed criterion will become the stability robustness criterion. The presented criterion is mathematically proved to be less conservative than the existing ones reported very recently in the literature. In [48], it is emphasized that the robust regional eigenvalue clustering analysis of singular systems is more complicated than that of standard state space systems, because that the robust regional eigen value-clustering analysis of uncertain singular systems should consider not only regional eigen value-clustering robustness but also system regularity and impulse-free (or causal) behaviours and the research on the regional eigen value clustering robustness analysis for linear discrete singular time-delay systems with parameter uncertainties is considerably rare and almost embryonic.In [49] analytical determination of the fundamental frequencies of vibrating polygonal plates with a concentric circular core is presented to extract the fundamental eigenvalue of the governing biharmonic boundary value problem. The method is then applied to wavy and polygonal plates with clamped and simply supported outer boundary conditions. Periodic plates can be considered with internal core as clamped, simply supported, and free circular cores. Approximate analytical formulations of the fundamental frequency for such plates with core are obtained. Fundamental frequency values of rectangular plates with boundary conditions simply supported-clamped boundary conditions (SC), free-clamped boundary conditions (FC), and simply supported-simply supported boundary conditions (SS) are compared with the existing literature and found in good agreement with results. The comparison is also made with hexagonal plates with simply supported-clamped boundary conditions (SC) and simply supported-simply supported boundary conditions (SS) boundary conditions, and the values are found to be in excellent agreement with the existing literature with the maximum difference being of the order of 2.7 percent. The present analysis can be used to obtain the higher frequencies by simply applying the perturbation method to higher modes and frequencies about the circular state.In [50] an asymptotic solution is presented for the natural frequencies, mode shapes of a cantilevered, coupled beam model. A three-term expansion is obtained for the frequency, which shows good agreement with exact values over a wide range of the stiffness parameter α. The first three mode shapes are also presented for large α and show good agreement for the first and second modes. Agreement for the third mode improves as α gets larger. The modeled structure behaves essentially as a shear beam except in the immediate neighborhood of the ends. Williams in [50] had found that the asymptotic results are a reasonable approximation to the exact results, particularly for large values α>20 of the stiffness parameter. Further, since the asymptotic results for the first mode are indistinguishable from the exact results, they present only the results for the second and third modes. The general character of the motion is essentially that of a shear beam: the mode shape is a sine wave. The effects of the Euler-Bernoulli beam component are to introduce boundary layers near each end of the beam within which the boundary conditions adjust from the exact conditions to those of a shear beam. If they pick x̃=3, x+=-3 as defining the essential width of the boundary layers, then this adjustment width is approximately equal to 3εL.In [51] Steigmann offered some observations on recent efforts to extract models for the stretching and bending of thin plates from three-dimensional finite elasticity. Using an asymptotic argument the author showed that recent work purporting to derive a nonstandard bending theory generates instead a correction to membrane theory of order thickness squared. The problem of deriving two-dimensional models from three-dimensional elasticity to describe the bending and stretching of plates and shells is one of the major open problems of mechanics. In the present work the author adopt the asymptotic method, which affords a systematic analysis of the questions of concern to us here. In particular he showed that recent work [52] based on ideas used in the method of Gamma convergence and purporting to discover a nonstandard bending energy in fact furnishes an order h2 correction to the leading-order membrane energy. Further developments concerned with the structure of a genuine bending energy valid for finite elastic strains are discussed in [53].The article [54] deals with the damage modeling of quasi-brittle interfaces such as the mortar or brick interfaces present in masonry walls. For this purpose, a model is developed, which takes the damage to the mortar joint into account. A quasi-fragile damage interface model is introduced using an asymptotic technique. This model memorizes some of the geometrical and mechanical characteristics of the interface, such as the thickness, elastic coefficients, normal and tangential stress, and the damage variable. Numerical simulations are performed using the Gyptis finite element software: academic cases involving traction and shear loads are presented. Here, authors described a model for studying quasi-brittle interfaces of this kind. The model is based on the study by Gambarotta and Lagomarsino [55] on masonry bricks. An asymptotic analysis was performed to obtain an interface model depending on six parameters, which relates the displacement jump nonlinearly to the stress vector along the interface. The stiffness of this link depends on a damage variable, and the damage is governed by an evolution equation. Two numerical algorithms were studied. The model was implemented in the Gyptis software, and academic tests were performed to test the robustness of the algorithms. In this paper, only academic examples are presented, and the results seem qualitatively in agreement with experimental data [56, 57].In [58] author addressed unsteady aerodynamic forces acting on a wing section oscillating in a steady incompressible (inviscid) uniform flow in presence of a distant flat ground. The problem of an oscillating wing section in (weak) ground effect has four fundamental length scales from which three independent dimensionless ratios can be constructed. One of them is the ratio δ of the transversal displacement of the wing to the wing chord, the other is the ratio ε of the wing (quarter) chord to its average distance from the ground, and the third is the ratio k of the wing (semi) chord to the distance traveled by the wing during one-oscillation period. As δ goes to zero, the wing ceases being a perturbation to the flow, as k goes to zero, the problem reduces to that of a nonoscillating wing section in ground effect, and as ε goes to zero, the problem reduces to that of an oscillating wing section out of ground effect. With the first two serving as small parameters, asymptotic series of the form
(3.15)δf0(k)+δε2f1(k)+δε2g(kε)f2(k)+⋯
have been constructed for the wing lift and pitching moment. In the case of heave oscillations, three-term series for the lift fits nicely the available numerical data for wide range of δ, ε, and k values. The solution for oscillating wing section in weak ground effect was constructed in two steps. In the first one, the problem was linearized with respect to δ, the ratio of the wing transversal displacement to its semichord. This step was essential in separating time-dependent constituents of the aerodynamic loads from time-independent ones, and, ultimately, in reducing the problem to that of an oscillating infinitesimally thin section. In the second step, the integral equation governing the problem of an oscillating infinitesimally thin section in ground effect was solved asymptotically using ε, the ratio of the wing quarter chord to the average distance between the wing trailing edge and the ground, as a small parameter. Formally, the three-term solution is asymptotically accurate to within terms of the order ε4 but without obtaining all terms in the series (or, what is probably simpler, comparing it with numerous numerical simulations) it is impossible to specify its practical applicability limit. Good agreement with numerical simulations we have obtained at the distances of 0.7 and 0.25 chords above the ground seems to support (but certainly not prove) this conjecture. The present theory should be applicable for aircraft wings during ground roll and low flight and even for racing cars high set front wings.In [59] Carpinteri and Paggi discussed the asymptotic analysis of the stress distribution around an elastic wedge-shaped domain which is one of the most fundamental problems in linear elasticity. Here authors talk about the epistemological steps in the relevance of the stress singularities and presented the review of the efforts and research developments towards a full appreciation of the mathematical and engineering relevance of stress-singularities. To keep the scope of the present article within reasonable limits, they have also restricted reader’s attention to the main geometrical configurations and mechanical conditions that can be used to effectively relieve the power of the stress-singularities with respect to the well-known square root singularity typical of linear elastic fracture mechanics. Among them in [59] shown the effect of the notch angle, the influence of the exponent of the stress-strain relationship in power law hardening materials, as well as the effect of the roughness of crack surfaces. Special attention has been devoted to the problem of stress-singularity in multi-material junctions and wedges which is a situation frequently observed in mechanical and composite engineering. Finally, a section on stress-singularities in non-homogeneous materials has been presented, which is a new emerging research field, more and more appealing for the scientific community due to the enormous advances in materials processing achieved during the last few years. Clearly, much work remains to be done in this research area in order to fully elucidate the potentials of the use of functionally graded materials for the removal of stress-singularities in junction problems.In [60] authors present an asymptotic solution for an infinite, adiabatic porous burner considering three different characteristic length scales: the solid-phase diffusion length-scale (lS), where the solid-phase heat conduction, gas-phase convection, and interphase heat transfer dominate the problem, the gas-phase diffusion length-scale (lG), where the gas-phase convection and diffusion dominate the problem, and the reaction length-scale (lR), where reaction and gas phase diffusion dominate the problem. Explicit solutions for the gas and solid temperatures and for the fuel and oxidant consumption were found as function of the problem parameters for thelSand lGcharacteristic length-scales. The description of the problem in the reaction length scale leads to an approximated expression for the flame velocity. The model was used to analyze the influence of the equivalence ratio ϕ, the conductivities ratio Γ, the matrix porosity ε, and the fuel Lewis number (LeF) on the flame structure. The combination among these four parameters defines the heat recirculation induced by the matrix and, consequently, the super adiabatic effect. Since the main focus of this work is the examination of the processes in the outer and the first inner regions, the simplest kinetic mechanism of one global step is adopted to represent the fuel and oxygen consumption. Then, the description of the reaction zone is obtained using the large activation energy asymptotic method. The description of the problem of the order of the gas-phase length scale is obtained using the boundary layer expansion. This work evaluates the influence of the equivalence ratio, the ratio of the solid to the gas thermal conductivities, the porosity of the medium, and the fuel Lewis number on such flames. Due to the simplifications assumed by the model the solution fails for extremely lean mixtures, when the heat transfer parameter N is large, for equivalence ratios near unity, because of the simplified solution for the reaction region, and for very high porosities (ε→1), since in this case the scale separation assumed by the model does not hold.In [61], authors show the efficiency of anisotropic adaptive meshes in the computation of very thin shell problems including boundary and internal layers phenomena. In a first part, they present general results concerning the structure of the singularities appearing inside these layers which depends on the geometrical nature of the shell, then established error estimates which reveal the performance and the efficiency of using anisotropic and adapted meshes inside the layers. Finally, authors used the finite element method for anisotropic elements and adaptive meshes in the computation of very thin shell problems. Classically, when thin shell models are deduced asymptotically from the three dimensional elasticity with the relative thickness ε as small parameter, they obtain either a membrane model or a pure bending model, according to whether the shell is geometrically rigid or not rigid. Using an adaptive anisotropic mesh procedure, they performed computations for different types of shells. The obtained results confirmed the relevance of anisotropic meshes (inside the layers) and also the necessity of an automatic mesh refinement procedure, considering the variety of singularities (i.e., order, propagation, sensitive problem) that may appear when ε↘0. The chosen examples for drastically different geometries exhibit the variety of situations that can be encountered and show the efficiency of anisotropic adaptive meshes for this kind of problems.Remark 3.1.
The study of the articles mentioned in the above section concludes the following points.(i)
Tailored finite point method gives high accurate results on a very coarse mesh in comparison of traditional finite element method or finite difference method.(ii)
B-spline collocation method is second order uniformly convergent and unconditionally stable.(iii)
Exponential B-spline collocation method is also second order uniformly convergent and gives more accurate results comparative to B-spline collocation method.(iv)
The advantage of an upwind scheme [34] for obtaining the global solution and normalized flux for singularly perturbed boundary value problems is that without any prior knowledge of the location of the boundary layer, we can generate an appropriate nonuniform mesh suitable for the layer type problems.(v)
Higher-order cubic B-spline collocation method produces a spline function which is useful to obtain the solution at any point of the interval whereas the finite difference method gives the solution only at the selected nodal points.(vi)
The extension of non-standard finite difference schemes [41] can be applied in different type of singularly perturbed problems, including turning point and nonlinear problems.(vii)
The uniform mesh method based on polynomial cubic spline for solving singularly perturbed boundary value problems gives better result in comparison to non uniform mesh tension spline methods.


## 4. Numerical Techniques to Solve Nonlinear Singular Perturbation Boundary Value Problems

In the present section, we will discuss singular perturbations in nonlinear problems. Due to the large variety of possible nonlinear differential equations arising in science and engineering we are forced to restrict our target to introduce the problems of physical nature. Here, we will discuss briefly the research articles in a sequence involving one-, two- and multidimensional singular non linear perturbation problems.In [62] Mo used a special and simple method; they study a class of singularly perturbed reaction diffusion nonlinear initial boundary value problems with a boundary perturbations. The author consider the following nonlinear singularly perturbed initial boundary value problem:
(4.1)ε∂u∂t-Lu=f(r,ϕ,u,ε),(t,r,ϕ)∈(0,T]×Ωε,(4.2)Bu≡g(u,ε2∂r∂ε,ε)=0,∂Ωε:r=r(ϕ,ε),(4.3)u=h(r,ϕ,ε),t=0,
where ε is a small positive parameter, L is a linear uniformly elliptic operator, Ωε={(r,ϕ)0≤r≤r(ϕ,ε)} denotes a bounded convex domain in R2, ∂Ωε signifies the smooth boundary of Ωε, and ∂/∂r is the outward derivative on ∂Ωε. Without loss of generality, it is useful to assume that the zeroth-order derivative of u is not present in Lu. Problem (4.1)–(4.3) is a nonlinear initial boundary value problem for reaction diffusion equation, and it constructs the asymptotic expansion of the solution and discusses its asymptotic behavior. Under suitable conditions and using the theory of differential inequalities the asymptotic solution of the initial boundary value problems is studied. Author obtained outer solution for nonlinear initial boundary value problem which may not satisfy the initial condition (4.3), so that they construct the initial layer corrective term.In [63] Amiraliyev et al. discussed a finite difference scheme on a special mesh Bakhvalov type for finding the numerical solution of the following form of singular perturbation boundary value problem which depends on a parameter:(4.4)εu′(x)+f(x,u,λ)=0,x∈Ω=(0,l],u(0)=A,u(1)=B,
where ε is a small positive parameter, A and B are any arbitrary constants. In this article the authors assumed that the function f(x,u,λ) is sufficiently differentiable in
(4.5)Ω̅×R2(Ω̅=Ω⋃{x=0}),0<α≤∂f∂u≤a*<∞,0<m1≤|∂f∂λ|≤M1<∞.
The solution of (4.4) indicates a pair {u(x),λ}∈C1(Ω¯)×R for which (4.4) is satisfied and near x=0 the function u(x) has a boundary layer of thickness O(ε) for ε≪1. Here authors analyzed some important properties of the analytical solution and verified these properties for a finite difference scheme developed on a special mesh of the Bakhvalov type and proved the scheme has uniform convergence in discrete maximum norm. This difference scheme has been demonstrated on two test problems. To compute the errors in test problems double mesh principle as the exact solution of test problems is not available. The error estimates are first order uniformly convergent which match the theoretical argument stated by the authors.In [64] a new reliable algorithm based on an adaptation of the standard homotopy perturbation method (HPM) is applied to the Chen system which is a three-dimensional system of ODEs with quadratic nonlinearities. The HPM is treated as an algorithm in a sequence of intervals for finding accurate approximate solutions to the Chen system, this technique as the Multistage Homotopy perturbation method (MHPM). In particular they look at the accuracy of the HPM as the Chen system changes from a nonchaotic system to a chaotic one. Comparisons between the standard HPM solutions and MHPM solution and the fourth-order Runge-Kutta numerical solutions were made. For both nonchaotic and chaotic cases studied and they found that the 15-term MHPM solutions on a larger time step achieved comparable accuracy compared with the fourth-order Runge-Kutta solutions on a much smaller time step. However, note that the orders of magnitude of the errors in the nonchaotic and chaotic cases differ and MHPM solutions were computed via a simple algorithm without any need for perturbation techniques, special transformations, linearization or discretization. Even though the MHPM works well on highly chaotic systems such as the Chen system, still care has to be taken on the choice of time span, time step, and the number of terms used. In [64] authors demonstrate the accuracy of the MHPM against the Maple’s built-in fourth-order Runge-Kutta procedure for the solutions of both nonchaotic and chaotic systems. The MHPM algorithm is coded in the computer algebra package Maple. The Maple environment variable “Digits” controlling the number of significant digits is set to 35 in all the calculations done in [64] and also the values of the parameters are assumed to be a=35 and c=28 with b=12 (nonchaotic) and b=3 (chaotic). The initial conditions are taken as x(0)=-10, y(0)=0 and z(0)=37.In [65] a method for the numerical solution of singularly perturbed two-point boundary value problems is proposed. The given boundary value problem has been converted to an initial value problem for system of two first-order ordinary differential equations, which is more convenient to handle, in the sense that one can use more number of mesh points to integrate the differential equation. In order to get an initial condition for the system, authors used the asymptotic approximate solution. This method is developed for the following singularly perturbed two-point boundary value problem: (4.6)Lu(x)≡εu′′(x)+a(x)u′(x)-b(x)u(x)=f(x),x∈D=(0,1),B0u(0)≡b(0)u(0)-a(0)u′(0)=-f(0),B1u(1)≡u(1)=B,
where ε>0 is a small parameter a,b and f are smooth functions such that a(x)≥a>0, b(x)≥b≥0, x∈D̅=[0,1]. Under these assumptions the boundary value problem (4.6) has a unique solution u(x) exhibiting a less severe boundary layer at x=0[1,10].The “less severe” means the solution u(x) of the boundary value problem (4.6), and its first derivatives are bounded uniformly, for all ε on the interval [0,1]. It may be noted that the reduced problem satisfies the boundary condition at x=0 exactly. Further, in [65] authors also developed two schemes to integrate the singularly perturbed system of initial value problem; the first method is a combination of the classical and exponentially fitted difference scheme, which is of order O(h). The second method is completely of exponential type and is adaptive single-step exponential scheme, which is of order O(h2), and the accuracy of the method mainly depends on the initial approximations and also the efficiency of the method. The accuracy of the second scheme is superior to the first scheme in all the tested cases even in the first iteration of the shooting method. Nonlinear problems can also be solved by the method after linearizing it using Newton quasilinearization method given in [66].In [67] authors studied the maximum number of limit cycles that can bifurcate from the period annulus surrounding the origin of a class of cubic polynomial differential systems using the averaging method. More precisely, they prove that the perturbations of the period annulus of the center located at the origin of the cubic polynomial differential system
(4.7)ẋ=-yf(x,y),ẏ=xf(x,y),
where f(x,y)=0 is a conic such that f(0,0)≠0, by arbitrary cubic polynomial differential systems provide at least six limit cycles bifurcating from the periodic orbits of the period annulus using only the first-order averaging method. There are essentially four usual methods for studying the number of limit cycles which bifurcate from the periodic orbits of a period annulus of a center. The first one is based on the Poincare return map, second is based on the Poincare-Pontrjagin-Melnikov integrals or abelian integrals that are equivalent in the plane. The third is based on the inverse integrating factor. Finally, the last one is based on the averaging theory. The first two methods only give information on the number of periodic orbits of the unperturbed system that turn into limit cycles with the perturbation. The last two methods additionally can give the shape of the bifurcated limit cycle up to any order in the perturbation parameter. Then, the first-order averaging method coincides with the abelian integral method for studying the number of limit cycles which bifurcate from the period annulus of the center located at the origin when ε=0. The averaging method gives a quantitative relation between the solutions of a nonautonomous periodic differential system and the solutions of its averaged differential system, which is an autonomous one. The given theorem-3 in [67] provides a first-order approximation in ε for the periodic solutions of a periodic differential system.The Multidimensional Lindstedt-Poincare (MDLP) method in [68] is extended to the nonlinear vibration analysis of axially moving beams which belong to the gyroscopic system. Galerkin method is used to discretize the governing equations. The forced response of an axially moving beam with internal resonance between the first two transverse modes and fundamental harmonic resonance are studied. The response curves exhibit the same internal resonance characteristics as that of nontransferring thin plates and beams because all these systems possess cubic nonlinearity and similar frequency distribution. Authors show the results of the MDLP method agree reasonably well with those obtained by the Incremental harmonic balance (IHB) method. They considered some typical examples and illustrated those the MDLP method is more straightforward and efficient than other perturbation methods for multiple degree-of-freedom systems such as the method of multiple scales and the Krylov-Bogoliubov-Mitropolsky method. The forced responses of an axially moving beam with the excitation frequency Ω near the first natural frequency ω10 were investigated. When the damping is small, all the response curves exhibit the same internal resonance characteristics as those of non-transferring thin plates and beams because all these systems possess cubic nonlinearity and similar frequency distribution. When the vibration amplitude is small, the predictions of the MDLP method are in good agreement with those of the IHB method. These two methods can complement each other. The method can be considered as a generalization of the Lindstedt-Poincare method to multidimensional systems and will be termed as the Multidimensional Lindstedt-Poincare (MDLP) method. This paper [68] starts with a brief description on the governing equation of the axially moving system followed by an introduction on the essence of the MDLP method. Typical cases of the axially moving beam problem have been investigated. Results obtained are compared with that obtained by the IHB method.In [69] Kadalbajoo and Kumar are devoted to the numerical study of the boundary value problems for nonlinear singularly perturbed differential-difference equations with small delay. Quasilinearization process is used to linearize the nonlinear differential equation. To obtain parameter uniform convergence, a piecewise-uniform mesh is considered, which is dense in the boundary layer region and coarse in the outer region. The method has shown to have almost second-order parameter-uniform convergence. The effect of small shift on the boundary layers has also been discussed. A Differential-difference equation (DDE) is also named as functional differential equation (FDE) which is classified into two categories, retarded and neutral. A delay differential equation is said to be of retarded DDE if the delay argument does not occur in the highest-order derivative term, otherwise it is said to be neutral DDE. In this article, they considered second-order nonlinear retarded DDE and B-spline collocation method with negative shift in the convection term has been carried out. The piecewise-uniform mesh is chosen to grasp the better approximation to the exact solution in the boundary layer region. Instead, if authors use the uniform mesh for solving the singular perturbation problem, not only it requires a larger number of mesh points in the given interval, but also the approximate solution oscillates about the exact solution in the boundary layer region when the value of ε decreases. So, authors considered piecewise-uniform mesh with more mesh points in the boundary layer region. In this article, numerical results are presented in support of the theoretical results and to show the effect of shift on the solution.Authors in [70] explored the idea of deriving a direct higher-order scheme (specially the one which uses variable step sizes because the one with the uniform step size can be derived from this trivially) by considering the following class of general singularly perturbed two-point boundary value problem:(4.8)Ly≡cεy′′(x)+a(x)y′(x)+b(x)y(x)=f(x),x∈[0,1],y(0)=η0,y(1)=η1,
where η0, η1 are given constants and cε=+εor -ε with 0<ε≤1. Further f(x),a(x), and b(x) are sufficiently smooth functions satisfying the conditions a(x)≥α>0, b(x)≥β>0. The indirect methods (those which do not use any acceleration of convergence techniques, e.g., Richardson’s extrapolation or defect correction) for such problems on a mesh of Shishkin type lead the error as ο(n-1Inn), where n denotes the total number of subintervals of [0,1]. In this article, they systematically describe a very simple and direct method which reduces the error to ο(n-2In2n). This method is proved to be ε-uniformly convergent with the above error bounds, on a piecewise uniform mesh of Shishkin type. Using a fitted mesh finite difference method on a mesh of Shishkin type, they described a numerical method for solving general singular perturbation problems and analyzed the method for ε-uniform convergence. Results obtained show that fitted mesh finite difference method is advantageous over the standard finite difference method. This method can also be extended for singularly perturbed non linear problems by applying it to the sequence of quasi-linearized equations. However, in view of the current interest of the researchers, instead of this quasi-linearization approach, they will prefer to have some direct methods for the general nonlinear singularly perturbed problems.A hyperbolic perturbation method is presented in [71] for determining the homoclinic solution of certain strongly quadratic nonlinear autonomous oscillators of the form
(4.9)ẍ+c1x+c2x2=εf(μ,x,ẋ),
where ε is a small positive parameter in which hyperbolic functions can be employed instead of the usual periodic functions in the perturbation procedure. In this article, the homoclinic orbit is not periodic when the time t→+∞ or t→-∞, the phase point on the phase path in the phase portrait approaches the same saddle point. After proving that the infinite periodic solution of the generating equation can be expressed by hyperbolic functions, hyperbolic functions instead of the usual elliptic functions are selected as the basis functions throughout the perturbation procedure. To show the essence and the effectiveness of the proposed method, typical examples of the generalized van der Pol equations are presented.The generalized Vander Pol oscillator is(4.10)f(μ,x,ẋ)=(μ+μ1x-μ2x2)ẋ,
where μ1 and μ2 are constants whilst μ is considered as a control parameter. To illustrate the accuracy of the method, its predictions are compared with those of Runge-Kutta method. All the homoclinic orbits obtained by this method are closed to those obtained by Runge-Kutta method at the critical parameter μ=μc even for moderately large value of ε which shows the hyperbolic perturbation method is reliable. The purpose of this research article is focused on obtaining the homoclinic solution. Similar perturbation procedure and assumption of the elliptic perturbation method can be employed in the present perturbation method. Here, the procedure of using R-K integration method to determine the value of parameter μ of the homoclinic orbit follows that of Merkin and Needham in [72]. The hyperbolic perturbation method presented in this article is an effective method for determining homoclinic solutions of certain strongly nonlinear autonomous oscillators in which hyperbolic functions are the exact homoclinic solution of the generating equation. Based on the functions, the hyperbolic perturbation method can lead to the analytical expression for the homoclinic solutions of the nonlinear autonomous oscillators. This hyperbolic perturbation method can be generalized for determining the heteroclinic solutions of certain strongly nonlinear autonomous oscillators. Comparing with the elliptic Lindstedt-Poincare method, the main distinguishing feature of this method is that the hyperbolic functions are employed instead of the elliptic functions. The differential and integral operations on the hyperbolic functions are easier than those of the elliptic functions.In [73] authors used the theories and methods of mathematical analysis and computer algebra, a reliable algorithm of finite difference method for solving singular perturbation problems was established, and a new MATLAB procedure EHSSP was implemented. Several linear and nonlinear problems are presented to illustrate the implementation of the program. The results show that the procedure EHSSP of patching method had advantages of easy handling, excellent property for operation and powerful competence. In [73] authors considered the nonlinear singular perturbation problem of the form(4.11)εy′′+f(x,y)y′+g(x,y)y=h(x),a≤x≤b,y(a)=α,y(b)=β,
where ε is a small positive parameter (0<ε≪1), a and b are given constants, f(x,y), g(x,y), and h(x) are assumed to be sufficiently continuously differentiable functions, and f(x,y)≥M>0 for a≤x≤b, where M is some positive constant. Now the original problem will be divided into two subproblems, an inner and outer region problem. Then both inner and outer region problem are solved as two-point boundary value problems by employing finite difference method. The proposed method is iterative on the terminal point. Algorithm EHSSP and mechanization of patching method for solving singularly perturbed differential equation using finite difference technique are presented. Authors suggested that construction of the finite difference differentiation matrix (FDDM) can be used to approximate any derivative order with any required accuracy of finite difference with very small execution time in MATLAB environment. This FDDM with the suggested initial values for the nonlinear solver fsolve reduced the execution time in the algorithm EHSSP.In [74] the homotopy perturbation method has been applied to obtain analytical approximate solutions for truly non linear problems that are conservative and periodic. The major conclusion is that this perturbation scheme provides excellent approximations to the solution of these nonlinear systems with high accuracy, and, in particular, the results for the cubic oscillator are valid for the whole solution domain. The analytical representations obtained using the homotopy perturbation technique give excellent approximations to the exact solutions for the whole range of values of oscillation amplitude. The approximate solutions of homotopy perturbation method are better than the approximate solutions obtained using the harmonic balance method for the second-order approximation. The relative error of the analytical approximate frequency obtained using the homotopy perturbation approach for the cubic oscillator is 0.031%, while the relative error is 0.41% when the harmonic balance method’s second-order approximation is considered. An interesting feature considered in [74] is the comparison between the analytical approximate solutions and the Fourier series expansion of the exact solution. This has allowed us to compare the coefficients for the different harmonics. In summary, He’s homotopy perturbation method is very simple in its principle, and it can be used to solve other conservative truly nonlinear oscillators with complex nonlinearities. The homotopy perturbation method can be applied to nonlinear oscillatory problems where a linear term does not exist, the nonlinear terms are not small, and there is no perturbation parameter.In [75] author presents the Extended Lindstedt-Poincare (EL-P) method, which applies multiple time variables to treat nonstationary oscillations arising in dynamical systems with cubic nonlinearities due to the slowly varied excitation parameters. The method is applied extensively in research of non-stationary vibrations of clamped-hinged beams. Recognizing a periodic nature of non-stationary oscillations, the new formulation is described by adding an additional, slow time scale beside time scales of the nonlinear system, which generally correspond to the incommensurate nonlinear frequencies of the response. Using this concept, a generalized approach of the study to the passage through fundamental, superharmonic and subharmonic resonances is introduced. Effects of slowly varying excitation frequency and slowly varying excitation amplitude on the non-stationary oscillations are studied with the computation of deviations from the stationary response. Although the method is formulated for N-dof dynamical systems having weak cubic nonlinearities, it is applied for non-stationary vibrations, where two-mode shape approximation of damped and undamped clamped hinged beam, respectively, is used and the simultaneous appearance of internal resonance is taken into account. Stability analysis of stationary solutions is performed, and comparisons of stationary resonance curves by results obtained with the incremental harmonic balance (IHB) method show good agreement. Branch tracing of the stationary solutions by arc-length continuation is implemented, which leads to the conclusion that the EL-P method is reliable and efficient. Quite so, in the presence of weak nonlinearities, this method has an evident advantage over the IHB method due to its deficiency of computing non-stationary responses, and this method can be easily generalized for treating dynamical systems with other types of nonlinearities, such as quadratic nonlinearities, as well as treating non-stationary combination resonance, which is beyond the scope of the conventional Lindstedt-Poincare method [76, 77].In [78] Öziş and Yıldırım present the iteration perturbation method proposed by He [79, 80] which is used to generate periodic solutions of Vander Pol oscillator with a forcing term, forcing oscillator with quadratic type damping, and Van der Pol oscillator with excitation term for the following second-order differential equation:(4.12)ẍ+ẋ+εf(x,ẋ,t)=0,
where ε need not be small and f(x,ẋ,t) is nonlinear analytic function of the displacement x, the velocity ẋ, and the time t. These nonlinear systems may be solved, with the combination of certain perturbation techniques, by “the method of elimination” similar to the method of elimination for linear differential systems. As it is well known, the objective of the method of elimination for linear differential systems is to eliminate dependent variable in succession until there remains only a single equation containing only one dependent variable. After the solution of remaining equation has been found, the other dependent variables can be found in turn, using the original differential equation or those that have appeared in the elimination process. Authors described how to obtain the amplitude and period for a system described by (4.12). In this article, they applied the iteration perturbation method to variety of forcing Vander Pol’s equations. This method is extremely simple in its principle, easy to apply, gives a very good accuracy even with the first-order approximation and simplest trial functions, and provides a powerful mathematical tool to the determination of limit cycles of more complex nonlinear systems.The parameter expansion method given in [81] is applied to a nonlinear oscillator with discontinuity. In this method only one iteration is required to obtain a highly accurate solution, which is valid for the whole solution domain.Authors considered the following nonlinear oscillator with discontinuity(4.13)u′′+u|u|=0,u(0)=A,u′(0)=0.
There exists no small parameter in the equation. Therefore, the traditional perturbation methods cannot be applied directly. So they used the parameter expansion method and re-wrote (4.13) in the following form and obtained the solution:(4.14)u′′+0.u+1.u|u|=0.
Comparison of the obtained solution with the exact one shows that the method is very effective, convenient, and this method can be applied to many other nonlinear oscillators. Furthermore, the obtained frequency is of high accuracy. The accuracy of frequency can be improved if they continue the solution procedure to a higher order; however, the amplitude obtained by this method is an asymptotic series, not a convergent one.In [82] Lan Xu introduced He’s work on asymptotic techniques and his parameter expanding method including modified Lindstedt-Poincare method and bookkeeping parameter method. They consider the following nonlinear oscillator:(4.15)mu′′+ω02u+εf(u,u′,u′′)=0,u(0)=A,u′(0)=0,
where ε is a small positive parameter and f is a polynomial function of its arguments.Various perturbation methods have been applied frequently to analyze (4.15). The perturbation methods are, however, limited to the case of small ε and mω02>0, that is, the associated linear oscillator must be statically stable so that linear and nonlinear response will be qualitatively similar. Author studied the cases when m=0 or/and ω02=0 and 0<ε<∞ by simple but powerful parameter expanding methods. Some remarkable virtues of the methods are exploited, and their applications for obtaining the higher-order approximate periodic solutions of strongly nonlinear oscillators are illustrated. In [83] authors solved the same equation (4.15) and obtained approximate frequencies by using parameter expansion method. They obtained sufficiently accurate solutions with a first-order approximation that are valid for whole domain. The traditional perturbation methods can give very good approximations when ε has a small value. However, when ε does not have a small value, such methods can yield unacceptably large errors.In [84] author applied He’s parameter-expanding method and obtained the nonlinear frequency-amplitude relationship of nonlinear oscillators. The obtained result is valid, when the amplitude tends to infinite. Traditional perturbation method provides us with a simple approach to the determination of the frequency-amplitude relationship, but the results are valid for weakly nonlinear systems as well as when the amplitude is very small. In 2006, He’s systematically summarized the two methods; the bookkeeping parameter method and He’s modification of Lindstedt-Poincare method, which were proposed by He in 2001 [85] and 2002 [86] respectively, called the technology as the parameter-expanding method [29], which is applied for strongly nonlinear equations. A parameter-expansion technology was embedded successfully with Lindstedt-Poincare method [86]; in case no small parameter or no expanding parameter existed in an equation, the bookkeeping parameter method was introduced [85].In [87] Marinca and Herisanu applied a new perturbation technique coupled with the iteration method. This procedure is obtained by combining the iteration methods into a new iteration procedure such that excellent approximate analytical solutions, valid for small as well as large values of oscillation amplitude, can be determined for nonlinear oscillations with a single degree of freedom. Authors compared the approximate periods obtained by authors’ procedure with the exact known periods. The results show that the approximations are of extreme accuracy. All of the examples presented here illustrate that the validity of this procedure for some nonlinear oscillations does not strongly depend upon the form of the restoring force.He’s Homotopy perturbation method has been adapted in [88], calculating the higher-order approximate periodic solutions for a nonlinear oscillator with discontinuity for which the elastic force term is an antisymmetric and quadratic term. This method works very well for the whole range of initial amplitudes, with excellent agreement of the approximate frequencies, and periodic solutions with the exact ones. Just one iteration leads to high accuracy of the solutions with a maximal relative error for the approximate period of less than 0.73% for all values of oscillation amplitude, while this relative error is as low as 0.040% when the second iteration is considered. The Homotopy perturbation method, a coupling of the traditional perturbation method and Homotopy in topology, deforms continuously a difficult problem to a simple problem which is easily solved. As a conclusion, authors said that the method has great potential and could be applied to other strongly nonlinear oscillators with nonpolynomial terms. Results presented in this article reveal that the method is very effective and convenient for conservative nonlinear oscillators for which the restoring force has a non-polynomial form.Perturbation methods depend on small parameters which are difficult to find for real-life nonlinear problems. He's parameterized perturbation method as discussed in [89], which does not need small parameters, is one of the newest analytical methods to solve nonlinear equations. In this article some nonlinear heat transfer equations are used as examples to illustrate the simple solution procedures of this method. Comparison of the results obtained by the new method with exact solutions reveals that the method is tremendously effective.In [90] Marinca and Herisanu implement a new analytical technique, He’s variational iteration method, valid for weakly and some strongly nonlinear oscillations. A correction functional is constructed by using a general Lagrange multiplier, which can be identified via the variational theory. This procedure is effective and accurate for nonlinear problems with approximations converging rapidly to accurate solutions. In this research article, the variational iteration method is applied to a general nonlinear problem, which can be initially approximated with unknown constants. The variational iteration method (VIM) was proposed by He [91, 92] which is a powerful mathematical tool for various kinds of nonlinear problems. The variational iteration method can also be used for solving three types of nonlinear partial differential equations, such as the coupled Schrodinger-KdV, generalizing KdV and shallow water equations. The amplitude and period of limit cycles for a modified Vander Pol oscillator are calculated by He’s variational method and the Krylov-Bogoliubov-Mitropolsky method.In [93] Öziş and Yıldırım applied He’s modified Lindstedt-Poincare method to determine the periodic solutions of oscillators in a u1/3 force. The elastic restoring forces are non-polynomial functions of the displacement within the nonlinear oscillator equations. These classes of equations represent a new class of nonlinear oscillating systems and were first studied in detail by Mickens [94]. Some new techniques that are extensions of the Lindstedt-Poincare perturbation method to strongly nonlinear systems, so called He’s Modified Lindstedt-Poincare method, were proposed by He [79]. In He’s Modified Lindstedt-Poincare method, a constant, rather than the nonlinear frequency, is expanded in powers of the expanding parameter to avoid the occurrence of secular terms in the perturbation series solution. The results obtained approximate solutions which are uniformly valid throughout solution domain and they are suitable not only for weakly nonlinear systems, but also for strongly nonlinear systems. The comparison of results with analytical solution provides confirmation for the validity of He’s Modified Lindstedt-Poincare method that it is extremely simple and easy to use.In [95] author considered a class of singularly perturbed boundary value problems for semilinear equations of fourth order with two parameters of the form(4.16)εd4ydx4+μp(x)d3ydx3+q(x)d2ydx2=f(x,y),a<x<b,y(r)(a)=Ar,r=0,1,y(r)(b)=Br,r=0,1,
where ε and μ are small positive parameters and Ar and Br(r=1,2) are constants. The problem (4.16) is a singularly perturbed problem with two small parameters.Author needed the following hypotheses:(H1)
μ2/ε→0 as ε→0;(H2)
the functionsp(x), q(x), and f(x,y) possess up to (m+1)th-order partial derivatives with regard to their variables in the corresponding domains, and min(p,q,fy)≥δ>0, where δ is a positive constant;(H3)
ξ=ε1/2, η=μ/ε1/2that is ε=ξ2,μ=ξη.Using hypotheses (H1)–(H3), the solutiony of the singularly fourth-order semilinear perturbed problem (4.16) with two parameters has been obtained and this solution may be expanded into the uniformly valid asymptotic expansion on x∈[a,b].Under suitable conditions, using the method of lower and upper solutions, the existence and the asymptotic behavior of the solution to the boundary value problem are studied and showed that the solution to the original singularly perturbed problem with two parameters has only one boundary layer.In [96] authors considered a new computational method for solving the following third-order singularly perturbed boundary value problem with a boundary layer in the reproducing kernel space: (4.17)εu′′′(x)+a(x)u′′(x)+b(x)u′(x)+c(x)u(x)=f(x),0≤x≤1,u′(0)=q,u(1)=0,u′(1)=0,
where 0<ε≪1, f(x)∈W21[0,1],a(x),b(x),c(x) are sufficiently smooth functions and a(x)≥-α, α>0. In this method, firstly the given problem is transformed into a system of two ODEs subject to suitable initial and boundary conditions. After that a zeroth-order asymptotic expansion is constructed for the solution of the given problem. Then the reduced terminal value problem (TVP) is solved analytically in the reproducing kernel space. Authors claimed that method is effective and easy to implement and results obtained by the method are found to be in good agreement with the exact solution within the boundary layer as well as away from the layer. Problem (4.17) has a boundary layer (at x=0) which is less severe because the boundary conditions are of the Neumann type.In [97] Zhou et al. present a singular perturbation approach to analyze the stability characteristics of the Active Disturbance Rejection Control (ADRC) system for nonlinear time variant plant. The closed-loop system is reformulated to allow the application of singular perturbation method; it enables the decomposition of the original system into a relatively slow and fast subsystem. Based on the decomposed subsystems, the composite Lyapunov function method is used to show that the closed-loop system achieves exponentially stable under certain conditions. Result shows that the system is exponentially stable, upon which a lower bound for the observer bandwidth is established. Control theory is the study of disturbances rejection problem. In regard to external disturbances, it is well known that good disturbance rejection is achieved with a high gain loop together with a high bandwidth, assuming that the plant is linear time invariant (LTI) and accurately described in a mathematical model. Things get a little complicated and interesting when such assumptions do not hold, as in most practical applications, where it is the internal disturbance that is most significant. In many regards, much of the literature on control can be seen as various responses to this dilemma. Two technology upgrades, figuratively speaking, are offered in the modern control framework: adaptive control [98] and robust control theory [99]. The combination of the two offers the third alternative robust adaptive control [100]. Therefore, as improvements of, rather than departure from, the model-based design paradigm, adaptive control and robust control, as solutions to the disturbance rejection problem, did not travel very far from their source, model-based classical control theory. Such approaches are deemed passive as they accept disturbances as they are and merely deal with them as one of the design issues. The main purpose in [97] is to show analytically how to use singular perturbation theory to obtain the sufficient condition of exponential stability of the closed-loop system of Active Disturbance Rejection Control (ADRC) and thus establish a lower bound for the observer bandwidth.In [101] authors introduced an optimal Galerkin method for solving the singularly perturbed high-order elliptic two-point boundary value problem of reaction-diffusion type using Hermite splines with knots adapted to the boundary layer behavior of the solution. Specifically, they present a sufficient condition on the mesh of grid points that ensures the corresponding approximate solution and the optimal order of uniform convergence in the energy norm. They also construct optimal meshes that satisfy the sufficient condition and choose the mesh size for each of the subintervals so that the error of the Hermite spline approximation on each subinterval is bounded in an optimal order uniformly with respect to the perturbation parameter. In this development, authors used estimates of the derivatives of the exact solution to guide the design of the mesh. Here, they describe the Galerkin methods for solving reaction-diffusion-type problems and propose a sufficient condition on the mesh size which ensures the optimal order of uniform convergence of the approximate solution. They construct a specific mesh which satisfies the condition and generates a Hermite spline approximate solution having the optimal order of uniform convergence.In [102] authors deal with the existence of solutions to the following singularly perturbed second-order three-point boundary value problem for nonlinear differential systems:(4.18)ε2x′′=f(t,x,x′),0<t<1,0<ε≪1,
together with the boundary conditions x(0,ε)=0, x(1,ε)=Px(η,ε), where ε is a small positive parameter, x and f are n-dimensional vectors,(4.19)P=diag⁡(p1,…;pn),pi<1(i=1,2,…,n),0<η<1.
They construct an appropriate generalized lower and upper solution pair, concepts defined in [102] and employ the Nagumo conditions and algebraic boundary layer functions to ensure the existence of solutions of the problem. The uniformly valid asymptotic estimate of the solutions is given as well. The differential systems have nonlinear dependence on all order derivatives of the unknown.Remark 4.1.
We have the following observations from the study of the articles mentioned in the above section.(i)
Fitted mesh finite difference method isε-uniform convergent and advantageous over the standard finite difference method.(ii)
The iteration perturbation method [78] provides a powerful mathematical tool to the determination of limit cycles of more complex nonlinear systems.(iii)
He’s variational iteration method is effective and accurate for weakly and strongly nonlinear problems and has some distinct advantages over usual approximation method.(iv)
Algorithm EHSSP and mechanization of patching method [73] for singular perturbed differential equation using finite difference technique had advantages of easy handling, excellent property for operation, and powerful competence.(v)
The multidimensional Lindstedt-Poincare method is reliable analytical method for periodic solutions of multiple degree of freedom systems and gives better result than the incremental harmonic balance method.(vi)
He’s homotopy perturbation method provides conservative and periodic approximate solutions that are better than the harmonic balance method for truly nonlinear problems.(vii)
The parameter-expanding method includes the modified Lindstedt-Poincare method and bookkeeping parameter method and provides sufficiently accurate solution which is valid in the whole solution domain.


## 5. Conclusion and Direction for Further Development

In the present paper computational methods for solving singularly perturbed boundary value problems in differential equations are briefly discussed; it contains and analyzes huge amount of the literature related to problems of differential equations of different orders from various fields. This paper is devoted to crux of various research articles published in refereed journals within last five years to get a better know-how of the state of art covering up the subject and will be an excellent reference for researchers to offer a state of the art of most active recent developments of methods for solving singular perturbation problems with their applications and remaining challenges in singular perturbation theory.From this survey we conclude that most of the authors considered singular perturbation problems in only ordinary differential equations. But it is well known that many phenomena in biology, chemistry, engineering, and physics can be described by boundary value problems associated with various types of partial differential equations or systems. When we associate a mathematical model with a phenomenon, we generally try to capture what is essential, retaining the important quantities and omitting the negligible ones which involve small parameters. The model that would be obtained by maintaining the small parameters is called the perturbed model.We are particularly interested to motivate the researchers for solving nonlinear partial differential equations, which have hardly been examined so far in the literature dedicated to singular perturbations. This paper suggests that the researchers fill this gap, since most of the applications are described by nonlinear models. Their asymptotic analysis is very interesting, but requires special methods and tools.



---

*Source: 404276-2011-11-22.xml*