Perturbation Expansion to the Solution of Differential Equations

Abstract

The main purpose of this chapter is to describe the application of perturbation expansion techniques to the solution of differential equations. Approximate expressions are generated in the form of asymptotic series. These may not and often do not converge but in a truncated form of only two or three terms, provide a useful approximation to the original problem. These analytical techniques provide an alternative to the direct computer solution. Before attempting to solve these problems numerically, one should have an awareness of the perturbation approach.

1. Introduction

The governing equations of physical, biological and economical models often involve features which make it impossible to obtain their exact solution. For instance, problems where we observe “a complicated algebraic equations”, “the occurrence of a complicated integral”, in case of differential equations (DE), “a varying coefficients or nonlinear term” sometimes problems with an awkwardly shaped boundary are tough to solve with the limited methods for finding analytical solutions. The main purpose of this chapter is to describe the application of perturbation expansion techniques to the solution of DE. Approximate expressions are generated in the form of asymptotic series. These may not and often do not converge but in a truncated form of only two or three terms, provide a useful approximation to the original problem. These analytical techniques provide an alternative to the direct computer solution. Before attempting to solve these DE numerically, one should have an awareness of the perturbation approach. An example of this occurs in boundary layer problems where there are regions of rapid change of quantities such as fluid velocity, temperature or concentration. Appropriate scaling of the boundary layer dimension is required before a numerical solution can be generated which will capture the behavior in the rapidly changing region.

When a large or small parameter occurs in a mathematical model of a process there are various methods of constructing perturbation expansions for the solution of the governing equations. Often the terms in the perturbation expansions are governed by simpler equations for which the exact solution techniques are available. Even if exact solutions cannot be obtained, the numerical methods used to solve the perturbation equations approximately are often easier to construct than the numerical approximation for the original governing equation.

First, we consider a model problem for which an exact solution is available against which the perturbation expansion can be compared. A feature of the perturbation expansions is that they often form divergence series. The concept of an asymptotic expansion will be introduced and the value of a truncated divergent series will be demonstrated.

2. Projectile motion

This example studies the effect of small damping on the motion of a particle. Consider a particle of mass $M$ which is projected vertically upward with an initial speed $U_0$. Let $U$ denote the speed at some general time $T$. If air resistance is neglected then the only force acting on the particle is gravity, $-\mathit{Mg}$, where $g$ is the acceleration due to gravity and the minus sign occurs because the upward direction is chosen to be the positive direction. Newton’s second law governs the motion of the projectile, i.e.,

$$M\frac{\mathit{dU}}{\mathit{dT}}=-\mathit{Mg}.$$E1

Integrating (1), we obtain the solution $U=C-\mathit{gT}$. The constant of integration is determined from the initial condition $U\left(0\right)=U_0$, so that

$$U=U_0-\mathit{gT}.$$E2

On defining the non-dimensional velocity $v$, and time $t$, by $v=U/U_0$ and $t=\mathit{gT}/U_0$, the governing equation becomes

$$\frac{\mathit{dv}}{\mathit{dt}}=-1,v\left(0\right)=1,$$E3

with the solution $v\left(t\right)=1-t$.

Taking account of the air resistance, and is included in the Newton’s second law as a force dependent on the velocity in a linear way, we obtain the following linear equation

$$M\frac{\mathit{dU}}{\mathit{dT}}=-\mathit{Mg}-\mathit{KU},$$E4

where the drag constant $K$ is the dimensions of masa/time. In the non-dimensional variables, it becomes

$$\frac{\mathit{dv}}{\mathit{dt}}=-1-\left(\frac{{\mathit{KU}}_0}{\mathit{Mg}}\right)v.$$E5

Let us denote the dimensionless drag constant by $\epsilon$, then the governing equation is

$$\frac{\mathit{dv}}{\mathit{dt}}=-1-\mathit{\epsilon v},v\left(0\right)=1,$$E6

where $\epsilon >0$ is a “small” parameter and the disturbances are very “small”. The damping constant $K$ in (4) is small, since $K$ has the dimensions of mass/time and a small quantity in units of kilograms per second.

2.1 Perturbation expansion

It is possible to solve (6) exactly since it is of variables separable form. Here, we solve by an iterative process, known as perturbation expansion for the solution.

Let $v^{\left(i\right)}$ denotes the $i$th iterate, which is obtained from the equation

$$\frac{{\mathit{dv}}^{\left(i\right)}}{\mathit{dt}}=-1-\epsilon v^{\left(i-1\right)},v^{\left(i\right)}\left(0\right)=1,$$E7

The justification for this iterative scheme is that the term $\mathit{\epsilon v}$ involves the small multiplying coefficient $\epsilon$, and so the term itself may be expected to be small. Thus, the term $\epsilon v^{\left(i\right)}$ which should appear on the RHS of (7) to make it exact, may be replaced by $\epsilon v^{\left(i-1\right)}$ with an error which is expected to be small.

The first iterate is obtained by neglecting the perturbation, thus

$$\frac{{\mathit{dv}}^{\left(0\right)}}{\mathit{dt}}=-1,v^{\left(0\right)}=1.$$

This is known as the unperturbed problem, and direct integration yields

$$v^{\left(0\right)}=1-t.$$

The next iterate $v^{\left(1\right)}$, satisfies

$$\frac{{\mathit{dv}}^{\left(1\right)}}{\mathit{dt}}=-1-\epsilon \left(1-t\right),v^{\left(1\right)}=1.$$

and integration yields

$$v^{\left(1\right)}=1-t\left(1+\epsilon \right)+\frac{1}{2}\epsilon t^2$$

Similarly, $v^{\left(2\right)}$ satisfies

$$\frac{{\mathit{dv}}^{\left(2\right)}}{\mathit{dt}}=-1-\epsilon \left[1-t\left(1+\epsilon \right)+\frac{1}{2}\epsilon t^2\right],v^{\left(2\right)}=1.$$

Direct integration yields the solution

$$v^{\left(2\right)}=1-t\left(1+\epsilon \right)+\epsilon \left(1+\epsilon \right)\frac{t^2}{2}-\frac{1}{6}{\epsilon }^2{t}^3.$$

Rearranging the terms in these iterates in ascending powers of $\epsilon$, we obtain

$$\begin{array}{l}{v}^{\left(0\right)}=1-t,\\ v^{\left(1\right)}=1-t+\epsilon \left(\frac{t^2}{2}-t\right),\\ v^{\left(2\right)}=1-t+\epsilon \left(\frac{t^2}{2}-t\right)+{\epsilon }^2\left(\frac{t^2}{2}-\frac{t^3}{6}\right).\end{array}$$E8

Clearly as the iteration proceeds the expressions are refined by terms which involve increasing powers of $\epsilon$. These terms become progressively smaller since $\epsilon$ is a small parameter. This is an example of a perturbation expansion. It will often be the case that perturbation expansions involve ascending integer powers of the small parameter, i.e., $\left\{{\epsilon }^0,{\epsilon }^1,{\epsilon }^2,\cdots \right\}$. Such a sequence is called an asymptotic sequence. Although this is the most common sequence which we shall meet, it is by no means unique. Examples of other asymptotic sequences are $\left\{{\epsilon }^{1/2},\epsilon ,{\epsilon }^{3/2},{\epsilon }^2,\cdots \right\}$ and $\left\{{\epsilon }^0,{\epsilon }^2,{\epsilon }^4,\cdots \right\}$. In each case the essential feature is that subsequent terms tend to zero faster than previous terms as $\epsilon \to 0$.

An alternative procedure to that of developing the expansion by iteration is to assume the form of the expansion at the outset. Thus, if we assume that the perturbation expansion involves the standard asymptotic sequence $\left\{{\epsilon }^0,{\epsilon }^1,{\epsilon }^2,\cdots \right\}$, then the solution $v$, which depends on the variable $t$, and the parameter $\epsilon$, is expressed in the form

$$v\left(t;\epsilon \right)={\epsilon }^0{v}_0\left(t\right)+{\epsilon }^1{v}_1\left(t\right)+{\epsilon }^2{v}_2\left(t\right)+\cdots .$$E9

The coefficients $v_0\left(t\right),v_1\left(t\right),\cdots$ of powers of $\epsilon$ are functions of $t$ only. Substituting expansion (9) in the governing Eq. (6) yields the following

$$\left\{\begin{array}{l}\frac{{\mathit{dv}}_0}{\mathit{dt}}+\epsilon \frac{{\mathit{dv}}_1}{\mathit{dt}}+{\epsilon }^2\frac{{\mathit{dv}}_2}{\mathit{dt}}+\cdots =-1-\epsilon v_0-{\epsilon }^2{v}_1-\cdots \\ v_0\left(0\right)+\epsilon v_1\left(0\right)+{\epsilon }^2{v}_2\left(0\right)+\cdots =1.\end{array}\right.$$E10

Thus, the coefficients of powers of $\epsilon$ can be equated on the left– and right–hand sides of (10):

$$\left\{\begin{array}{ll}{\epsilon }^0:& \frac{{\mathit{dv}}_0}{\mathit{dt}}=-1,v_0\left(0\right)=1,\\ {\epsilon }^1:& \frac{{\mathit{dv}}_1}{\mathit{dt}}=-v_0,v_1\left(0\right)=0\\ {\epsilon }^2:& \frac{{\mathit{dv}}_2}{\mathit{dt}}=-v_1,v_2\left(0\right)=0,\mathrm{etc}.\end{array}\right.$$E11

The proof of validity of this fundamental procedure can be developed by first setting $\epsilon =0$ in (10) which yields the first equation of (11). This result allows the first member of the left– and right–hand side of Eq. (10) to be removed. Then, after dividing the remaining terms by $\epsilon$ we obtain the equation

$$\frac{{\mathit{dv}}_1}{\mathit{dt}}+\epsilon \frac{{\mathit{dv}}_2}{\mathit{dt}}+\cdots =-v_0-\epsilon v_1-\cdots$$

This is valid for all nonzero values of $\epsilon$ so that on taking the limit as $\epsilon \to 0$ we obtain the second equation of (11). Repeating the procedure, we obtain the other equations.

Integrating the equations in (11), we obtain

$$v_0=1-t,v_1=t^2/2-t,v_2=t^2/2/-t^3/6.$$

Using these values in (9), we obtain that

$$v\left(t;\epsilon \right)=1-t+\epsilon \left(t^2-t\right)+{\epsilon }^2\left(t^2/2-t^3/6\right)+\cdots$$E12

This is the same as the expansion (8) which is generated by iteration.

The IVP (6) can be solved exactly as

$$v\left(t\right)=\left[\left(1+\epsilon \right)e^{-\mathit{\epsilon t}}-1\right]{\epsilon }^{-1}.$$

The perturbation expansion can be obtained from (12) by replacing the exponential function by its Maclaurin expansion, i.e.,

$$v\left(t\right)=\frac{1}{\epsilon }\left[1-\mathit{\epsilon t}+\frac{{\epsilon }^2{t}^2}{2}-\frac{{\epsilon }^3{t}^3}{6}+\cdots +\epsilon -{\epsilon }^{2}t+\frac{{\epsilon }^3{t}^2}{2}+\cdots -1\right]$$E13

$$=\left(1-t\right)+\epsilon \left(\frac{t^2}{2}-t\right)+{\epsilon }^2\left(\frac{t^2}{2}-\frac{t^3}{6}\right)+\cdots$$E14

This is the same as the expansion (12). Thus, the perturbation expansion approach is justified in this case. One can refer the books [1, 2].

3. Asymptotics

The letters $O$ and $o$ are order symbols. They are used to describe the rate at which functions approach limit values. We will consider the types of limit values, namely zero, a finite number but nonzero and infinite.

If a function $f\left(x\right)$ approaches a limiting value at the same rate of another function $g\left(x\right)$ as $x\to x_0$, then we write

$$f\left(x\right)=O\left(g\left(x\right)\right),\mathrm{as}x\to x_0$$E15

The functions are said to be of the same order as $x\to x_0$. The test for this is the limit of the ratio. Thus, if $\underset{x\to x_0}{lim}\frac{f\left(x\right)}{g\left(x\right)}=C,$ where $C$ is finite, then we say (15) holds.

For example, we have the following functions:

$$\begin{array}{ll}{x}^2=O\left(x\right),& \mid x\mid <2,\\ sin\left(x\right)=O\left(\sqrt{x}\right),& x\to 0,\\ sin\left(x\right)=O\left(x\right),& -\infty <x<\infty .\end{array}$$

The expression

$$f\left(x\right)=o\left(g\left(x\right)\right),\mathrm{as}x\to x_0$$E16

means that $\underset{x\to x_0}{lim}\frac{f\left(x\right)}{g\left(x\right)}=0.$ This is a stronger assertion that the corresponding $O$–formula. The relation (16) implies the relation (15), as convergence implies boundedness from a certain point onwards.

We have the following functions satisfy the $o$–relation:

$$\begin{array}{ll}cos\left(x\right)=1+o\left(x\right),& \mid x\mid <2,\\ e^x=1+o\left(x\right),& x\to 0\\ n!=e^{-n}\cdot n^n\sqrt{2\mathit{\pi n}}\left(1+o\left(1\right)\right),& n\to \infty .\end{array}$$

3.1 Asymptotic expansions

Consider the expansion

$$f\left(x\right)=a_0+\frac{a_1}{x}+\frac{a_2}{x^2}+\cdots +\frac{a_N}{x^N}+R_N,$$E17

is an asymptotic expansion as $x\to \infty$, if, for any $N$,

$$R_N=O\left(\frac{1}{x^{N+1}}\right),\mathrm{as}x\to \infty$$E18

The following expansion is used when (17) and (18) hold,

$$f\left(x\right)\sim \sum _{n=0}^{\infty }\frac{a_n}{x^n},\mathrm{as}x\to \infty$$E19

Here, $\underset{n\to \infty }{lim}{R}_N=0$, for any value of N.

The sequence $\left\{1,1/x,1/x^2,\cdots \right\}$ is an asymptotic sequence as $x\to \infty$. The characteristic feature of such sequences is that each member is dominated by the previous member. In constructing examples it is easier to deal with the limit zero than any other. Thus, for the case $x\to \infty$, we let $\epsilon =1/x$, which for $x\to x_0$, we let $\epsilon =x-x_0$ so that without loss of generality we may confirm our attention to the limit $\epsilon \to 0$. The standard asymptotic sequence is $\left\{1,\epsilon ,{\epsilon }^2,\cdots \right\}$ as $\epsilon \to 0$. If we let ${\delta }_n\left(\epsilon \right)$ represent members of an asymptotic sequence $\left\{{\delta }_0\left(\epsilon \right),{\delta }_1\left(\epsilon \right),\cdots \right\}$ as $\epsilon \to 0$, then the following condition must hold

$${\delta }_{n+1}\left(\epsilon \right)=o\left({\delta }_n\left(\epsilon \right)\right),\mathrm{as}\epsilon \to 0.$$

Some examples of asymptotic sequences are

1. $\left\{1,sin\left(\epsilon \right),{\left(sin\left(\epsilon \right)\right)}^2,{\left(sin\left(\epsilon \right)\right)}^3,\cdots \right\}$, here we have

   $$\underset{\epsilon \to 0}{lim}\frac{{\delta }_{n+1}}{{\delta }_n}=\underset{\epsilon \to 0}{lim}sin\left(\epsilon \right)=0.$$

2. $\left\{1,ln\left(1+\epsilon \right),ln\left(1+{\epsilon }^2\right),ln\left(1+{\epsilon }^3\right),\cdots \right\}$, with ${\delta }_0=1,{\delta }_n=ln\left(1+{\epsilon }^n\right)n\ge 1$, we have

$$\begin{array}{l}\underset{\epsilon \to 0}{lim}\frac{{\delta }_1}{{\delta }_0}=\underset{\epsilon \to 0}{lim}ln\left(1+\epsilon \right)=0,\\ \underset{\epsilon \to 0}{lim}\frac{{\delta }_{n+1}}{{\delta }_n}=\underset{\epsilon \to 0}{lim}\frac{ln\left(1+{\epsilon }^{n+1}\right)}{ln\left(1+{\epsilon }^n\right)}=\underset{\epsilon \to 0}{lim}\frac{{\epsilon }^{n+1}+O\left({\epsilon }^{2n+2}\right)}{{\epsilon }^n+O\left({\epsilon }^{2n}\right)}=0.\end{array}$$

The general expression for an asymptotic expansion of a function $f\left(\epsilon \right)$, in terms of an asymptotic sequence ${\delta }_n\left(\epsilon \right)$ is

$$f\left(x\right)\sim \sum _{n=0}^{\infty }{a}_n{\delta }_n\left(\epsilon \right),\mathrm{as}\epsilon \to 0,$$E20

where the coefficients $a_n$ are independent of $\epsilon$. The expression (20) involving the symbol $\sim$, means that for all $N$,

$$f\left(x\right)=\sum _{n=0}^N{a}_n{\delta }_n\left(\epsilon \right)+R_N,$$E21

where

$$R_N=O\left[{\delta }_{N+1}\left(\epsilon \right)\right],\mathrm{as}\epsilon \to 0,$$E22

$$a_n=\underset{\epsilon \to 0}{lim}\left(\frac{f\left(\epsilon \right)-{\sum }_{n=0}^{N-1}{a}_n{\delta }_n\left(\epsilon \right)}{{\delta }_N\left(\epsilon \right)}\right).$$E23

If a function possesses an asymptotic expansion involving the sequence $\left\{{\delta }_0\left(\epsilon \right),{\delta }_1\left(\epsilon \right),\cdots \right\}$ then the coefficients $a_n$ of the expansion (21) given by the expression (24) are unique. However, another function may share the same set of coefficients. Thus, while functions have unique expansions, an expansion does not correspond to a unique function.

Consider a function $f\left(x;\epsilon \right)$, which depends on both an independent variable $x$, and a small parameter $\epsilon$. Suppose that $f\left(x;\epsilon \right)$ is expanded using an asymptotic sequence $\left\{{\delta }_n\left(\epsilon \right)\right\}$,

$$f\left(x;\epsilon \right)=\sum _{n=0}^N{a}_n\left(x\right){\delta }_n\left(\epsilon \right)+R_N\left(x;\epsilon \right).$$E24

The coefficients of the gauge functions ${\delta }_n\left(\epsilon \right)$ are functions of $x$, and the remainder after $N$ terms is a function of both $x$ and $\epsilon$. For this to be an asymptotic expansion, we require

$$R_N\left(x;\epsilon \right)=O\left[{\delta }_{N+1}\left(\epsilon \right)\right],\mathrm{as}\epsilon \to 0.$$E25

Refer [3, 4] for more details. For (24) to be a uniform asymptotic expansion the ultimate proportionality between $R_N$ and ${\delta }_{N+1}$ must be bounded by a number independent of $x$, i.e.,

$$\mid R_N\left(x;\epsilon \right)\mid \le K\mid {\delta }_{N+1}\left(\epsilon \right)\mid ,$$E26

for $\epsilon$ in the neighborhood near zero, where $K$ is a fixed constant.

An example of a uniform asymptotic expansion is $f\left(x;\epsilon \right)=\frac{1}{1-\epsilon sin\left(x\right)}.$

An example of a nonuniform expansion is

$$f\left(x;\epsilon \right)\sim \sum _{n=0}^N{x}^n{\epsilon }^n+R_N\left(x;\epsilon \right),\mathrm{as}\epsilon \to 0.$$E27

Here, one cannot find a fixed $K$ which satisfy $\mid R_N\mid \le K\mid {\epsilon }^{N+1}\mid ,$ because for any choice of $K$, $x$ can be chosen so that $x^{N+1}$ exceeds this value.

3.2 Nonuniformity

The expansion (27) becomes nonuniform when subsequent terms are no longer small corrections to previous terms. This occurs when subsequent terms are of the same order or of dominant order than previous terms. Subsequent terms dominate previous terms for larger $x$, for example, when $x=O\left(1/{\epsilon }^2\right)$. The expansion is valid for $x=O\left(1\right)$ since then subsequent terms decrease by a factor of $\epsilon$. The expansion remains valid for large $x$, provided $x$ is not as large as $1/\epsilon$. For instance, the expansion is valid for $x=O\left(1/\sqrt{\epsilon }\right)$, as $\epsilon \to 0$.

The critical case is such that subsequent terms are of the same order. This determines the region of nonuniformity. In (27), the region of nonuniformity occurs when $\mathit{\epsilon x}=O\left(1\right)$, i.e., $x=O\left({\epsilon }^{-1}\right)$, as $\epsilon \to 0$.

3.2.1 Sources of nonuniformity

There are two common reasons for nonuniformities in asymptotic expansions, they are

1. Infinite domains which allow long-term effects of small perturbations to accumulate.

2. Singularities in governing equations which lead to localized regions of rapid change.

Consider the nonlinear Duffing equation

$$\left\{\begin{array}{l}\frac{d^{2}u}{{\mathit{dt}}^2}+u+\epsilon u^3=0,t\in \left[0,\infty \right)\\ u\left(0\right)=a,\frac{\mathit{du}}{\mathit{dt}}\left(0\right)=0.\end{array}\right.$$E28

Suppose the solution may be expanded using the standard asymptotic sequence

$$u\left(t;\epsilon \right)\sim u_0\left(t\right)+\epsilon u_1\left(t\right)+{\epsilon }^2{u}_2\left(t\right)+\cdots .$$E29

On substituting this in (28) and in the initial conditions, we get

$$\left\{\begin{array}{l}\frac{d^2{u}_0}{{\mathit{dt}}^2}+\epsilon \frac{d^2{u}_1}{{\mathit{dt}}^2}+\cdots +u_0+\epsilon u_1+\cdots +\epsilon u_0^3+\cdots \sim 0,\\ u_0\left(0\right)+\epsilon u_1\left(0\right)+\cdots =a+0\cdot \epsilon +\cdots ,\\ \frac{{\mathit{du}}_0}{\mathit{dt}}\left(0\right)+\epsilon \frac{{\mathit{du}}_1}{\mathit{dt}}\left(0\right)+\cdots =0+0\cdot \epsilon +\cdots .\end{array}\right.$$

Equating like of powers of $\epsilon$ on both sides, we get

$$\left.\begin{array}{ll}O\left(1\right):& \frac{d^2{u}_0}{{\mathit{dt}}^2}+u_0=0,\\ & u_0\left(0\right)=a,\frac{{\mathit{du}}_0}{\mathit{dt}}\left(0\right)=0,\end{array}\right\}$$E30

and

$$\left.\begin{array}{ll}O\left(\epsilon \right):& \frac{d^2{u}_1}{{\mathit{dt}}^2}+u_1=-u_0^3,\\ & u_1\left(0\right)=0,\frac{{\mathit{du}}_1}{\mathit{dt}}\left(0\right)=0.\end{array}\right\}$$E31

Solving Eqs. (30) and (31), we obtain

$$u\sim acos\left(t\right)+\epsilon \left[\frac{a^3}{32}\left(cos\left(3t\right)-cos\left(t\right)\right)-\frac{3a^3}{8}tsin\left(t\right)\right]+\cdots .$$E32

The term $tsin\left(t\right)$ in the expansion (32) is called a secular term. It is an oscillating term of growing amplitude. All other terms are oscillating of fixed amplitude. The secular term leads to a nonuniformity for large $t$. The region of nonuniformity is obtained by equating the order of the first and second terms,

$$cos\left(t\right)=O\left(\mathit{\epsilon t}sin\left(t\right)\right),\mathrm{as}\epsilon \to 0.$$

The trigonometric functions are treated as $O\left(1\right)$ terms. Thus, the region of nonuniformity is $t=O\left(1/\epsilon \right)$, as $\epsilon \to 0$.

The second common source of nonuniformities is associated with the presence of singularities. Consider, the following initial-value problem:

$$\left\{\begin{array}{l}\epsilon \frac{\mathit{dy}}{\mathit{dx}}+y=e^{-x},x>0\\ y\left(0\right)=2,\end{array}\right.$$E33

where $\epsilon >0$ is a small parameter. Suppose $y$ has the expansion

$$y\sim y_0\left(x\right)+\epsilon y_1\left(x\right)+{\epsilon }^2{y}_2\left(x\right)+\cdots .$$E34

Substituting (34) in (33), we have

$$\left\{\begin{array}{l}\epsilon \left(\frac{{\mathit{dy}}_0}{\mathit{dx}}+\epsilon \frac{{\mathit{dy}}_1}{\mathit{dx}}+\cdots \right)+\left(y_0+\epsilon y_1+\cdots \right)=e^{-x},\\ y_0\left(0\right)+\epsilon y_1\left(0\right)+\cdots =2.\end{array}\right.$$E35

Equating coefficients of like powers of $\epsilon$ on both sides, we get

$$\begin{array}{ll}O\left(1\right):& y_0=e^{-x},y_0\left(0\right)=2,\\ O\left(\epsilon \right):& y_1=-\frac{{\mathit{dy}}_0}{\mathit{dx}}=e^{-x},y_1\left(0\right)=0,\\ O\left({\epsilon }^2\right):& y_2=-\frac{{\mathit{dy}}_1}{\mathit{dx}}=e^{-x},y_2\left(0\right)=0.\end{array}$$

Clearly, $y_0$ cannot satisfy the boundary condition $y_0\left(0\right)=2$ as no constant of integration is available because the equation determining $y_0$ is an algebraic equation not a differential equation, and no additional conditions are required. Thus, we have obtain the expression

$$y\sim e^{-x}+\epsilon e^{-x}+{\epsilon }^2{e}^{-x}+\cdots ,$$E36

but the initial condition $y\left(0\right)=2$ has not been satisfied.

The unperturbed problem, obtained by setting $\epsilon =0$ is not a DE, but an algebraic equation $y=e^{-x}$. This cannot satisfy an arbitrarily imposed condition at $x=0$. For any nonzero value of $\epsilon$, (33) becomes a first-order DE which can satisfy an initial condition. This is an example of a singular perturbation problem (SPP), where the behavior of the perturbed problem is very different from that of the unperturbed problem.

Thus, the perturbation expansion (36) is a good approximation of the exact solution away from the region $x=0$. To see this, let us compare (36) with the following exact solution:

$$y_{\mathit{ex}}=\frac{1-2\epsilon }{1-\epsilon }{e}^{-x/\epsilon }+\frac{e^{-x}}{1-\epsilon }=\left[\left(1-\epsilon -{\epsilon }^2-\cdots \right)e^{-x/\epsilon }\right]+\left[\left(1+\epsilon +{\epsilon }^2+\cdots \right)e^{-x}\right].$$E37

The perturbation expansion (36) generates the second member of (37), but not the first member. The coefficient $e^{-x/\epsilon })$ is a rapidly varying function which takes the value of unity at $x=0$, and rapidly decays to zero for $x>0$. Clearly, $y_0$ provides a good approximation away from the region $x=0$. The region near $x=0$ is called the boundary layer. These regions usually occur when the highest order derivative of a DE is multiplied by a small parameter. The unperturbed problem, obtained by setting $\epsilon =0$ is of lower order and consequently cannot satisfy all the boundary conditions. This leads to boundary layer regions where the solution varies rapidly in order to satisfy the boundary condition.

Boundary layers are regions of nonuniformity in perturbation expansions of the form (36).

4. Boundary layer

Boundary layers are regions in which a rapid change occurs in the value of a variable. Some physical examples include “the fluid velocity near a solid wall”, “the velocity at the edge of a jet of fluid”, “the temperature of a fluid near a solid wall.” Ludwig Prandtl pioneered the subject of boundary layer theory in his explanation of how a quantity as small as the viscosity of common fluids such as water and air could nevertheless play a crucial role in determining their flow. The viscosity of many fluids is very small and yet taking account of this small quantity is vital. The essential point is that the viscous term involves higher order derivatives so that its omission necessitates the loss of a boundary condition. The ideal flow solution allow slip to occur between a solid and fluid. In reality the tangential velocity of a fluid relative to a solid is zero. The fluid is brought to rest by the action of a tangential stress resulting from the viscous force.

Mathematically the occurrence of boundary layers is associated with the presence of a small parameter multiplying the highest derivative in the governing equation of a process. A straightforward perturbation expansion using an asymptotic sequence in the small parameter leads to differential equations of lower order than the original governing equation. In consequence not all of the boundary and initial conditions can be satisfied by the perturbation expansion. This is an example of what is commonly referred to as a singular perturbation problem. The technique for overcoming the difficulty is to combine the straightforward expansion, which is valid away from the layer adjacent to the boundary. The straightforward expansion is referred to as the outer expansion. The inner expansion associated with the boundary layer region is expressed in terms of a stretched variable, rather than the original independent variable, which takes due account of the scale of certain derivative terms. The inner and outer expansions are matched over a region located at the edge of the boundary layer. The technique is called the method of matched asymptotic expansions.

Consider the following two-point boundary value problem:

$$\left\{\begin{array}{l}\epsilon \frac{d^{2}u}{{\mathit{dx}}^2}+\frac{\mathit{du}}{\mathit{dx}}=2x+1,x\in \left(0,1\right)\\ u\left(0\right)=1,u\left(1\right)=4,\end{array}\right.$$E38

where $\epsilon >0$ is a small parameter. If we assume that $u$ possesses a straightforward expansion in powers of $\epsilon$,

$$u\left(x;\epsilon \right)\sim u_0\left(x\right)+\epsilon u_1\left(x\right)+{\epsilon }^2{u}_2\left(x\right)+\cdots ,$$E39

then the equations associated with powers of $\epsilon$ leads to

$$O\left(1\right):\frac{{\mathit{du}}_0}{\mathit{dx}}=2x+1,$$E40

$$O\left({\epsilon }^n\right):\frac{{\mathit{du}}_n}{\mathit{dx}}=-\frac{d^2{u}_{n-1}}{{\mathit{dx}}^2},\text{for}n=\mathrm{1,2,3},\cdots .$$E41

and the boundary conditions require

$$\begin{array}{l}{u}_0\left(0\right)+\epsilon u_1\left(0\right)+\cdots \sim 1+\epsilon \cdot 0+\cdots ,\\ u_0\left(1\right)+\epsilon u_1\left(1\right)+\cdots \sim 4+\epsilon \cdot 0+\cdots ,\end{array}$$

which leads to

$$\begin{array}{l}{u}_0\left(0\right)=1,u_0\left(1\right)=4,\\ u_n\left(0\right)=0,u_n\left(1\right)=0,\text{for}n=1,2,\cdots .\end{array}$$E42

Equation (42) require that each $u_n\left(x\right)$ satisfy two boundary conditions. This is in general impossible since Eqs. (41) and (42) governing each $u_n$ are of first-order. Now the question is which one of the boundary condition has to be taken into account. We will find out that the boundary condition at $x=0$ must be abandoned and consequently the expansion (39) is invalid near $x=0$.

The general solution of (42) is $u_0\left(x\right)=x^2+x+C$, using the boundary condition $u_0\left(1\right)=4$, we obtain

$$u_0\left(x\right)=x^2+x+2.$$

From (42), we obtain the equations

$$\begin{array}{ll}\frac{{\mathit{du}}_1}{\mathit{dx}}=-2,& u_1\left(1\right)=0,\\ \frac{{\mathit{du}}_2}{\mathit{dx}}=0,& u_2\left(1\right)=0,\end{array}$$

and its solutions are

$$u_1\left(x\right)=-2\left(x-1\right),u_n\left(x\right)=0,n\ge 2.$$

Therefore, the outer expansion is

$$u^{\text{out}}\left(x;\epsilon \right)=\left(x^2+x+2\right)+\epsilon 2\left(1-x\right),$$E43

where ‘out’ label is used to indicate that the solution is valid away from the region near $x=0$. Clearly $u^{\text{out}}$ fails to satisfy the boundary condition at $x=0$. The reason why the outer solution is of use is that it closely follows the exact solution of the problem except in a narrow region near $x=0$, where the exact solution changes rapidly in order to satisfy the boundary condition.

The exact solution of the BVP (38) can be obtained as

$$u\left(x\right)=A+{\mathit{Be}}^{-x/\epsilon }+x^2+x\left(1-2\epsilon \right).$$E44

The constants $A$ and $B$ are determined from the boundary conditions:

$$\left\{\begin{array}{l}A+B=1,\\ A+{\mathit{Be}}^{-1/\epsilon }+2-2\epsilon =4\end{array}\right.$$E45

We know that $e^{-1/\epsilon }=o\left({\epsilon }^N\right)$, as $\epsilon \to 0$, for all $N$. This means that the exponential term tends to zero faster than any power of $\epsilon$, as $\epsilon \to 0$. It is called a transcendentally small term (T.S.T.) and can always be neglected since its contribution is asymptotically always less than any power of $\epsilon$. Thus, (45) gives

$$A=2\left(1+\epsilon \right),B=-\left(1+2\epsilon \right),$$

and the exact solution is

$$u^{\mathrm{ex}}\left(x\right)=2\left(1+\epsilon \right)-\left(1+2\epsilon \right)e^{-x/\epsilon }+x^2+x\left(1-2\epsilon \right),$$E46

after rearranging the terms in asymptotic order, we obtain

$$u^{\mathrm{ex}}\left(x\right)=\left(x^2+x+2\right)-e^{-x/\epsilon }+\epsilon \left[2\left(1-x\right)-2e^{-x/\epsilon }\right].$$E47

Comparing the exact solution with the outer expansion shows that the terms involving $e^{-x/\epsilon }$ are absent. The effect of these terms is negligible when $x=O\left(1\right)$. But, when $x=O\left(\epsilon \right)$, then $e^{-\mathit{x\epsilon }}=O\left(1\right)$. It is clear that as $\epsilon \to 0$ the region in which the outer solution departs from the exact solution becomes arbitrarily close to $x=0$ with a thickness $O\left(\epsilon \right)$. This region is called the boundary layer.

The behavior of the exact solution and the zeroth-order term of the outer expansion are plotted in Figure 1 for various values of $\epsilon$.

Figure 1.

Exact solution of (38) for various values of $\epsilon$.

By differentiating the leading order term $u_0^{\mathrm{ex}}$, of the exact solution, we have

$$\begin{array}{l}{u}_0^{\mathrm{ex}}=x^2+x+2-e^{-x/\epsilon }\\ \frac{{\mathit{du}}_0^{\mathrm{ex}}}{\mathit{dx}}=2x+1+\frac{1}{\epsilon }{e}^{-x/\epsilon }\\ \frac{d^2{u}_0^{\mathrm{ex}}}{{\mathit{dx}}^2}=2-\frac{1}{{\epsilon }^2}{e}^{-x/\epsilon }\end{array}$$

Outside the boundary layer, i.e., for $x=O\left(1\right)$, we have $e^{-x/\epsilon }=o\left({\epsilon }^N\right),\forall N$, so ${\epsilon }^{-1}{e}^{-x/\epsilon }$ and ${\epsilon }^{-2}{e}^{-x/\epsilon }$ are also transcendentally small. Within the boundary layer when $x=O\left(\epsilon \right)$, we have $e^{-x/\epsilon }=O\left(1\right)$. The order of $u_0^{\mathrm{ex}}$ and its derivatives are given below:

$$\begin{array}{ccc}& \text{Outside}\mathrm{BL}& \text{Inside}\mathrm{BL}\\ u_0^{\mathrm{ex}}& O\left(1\right)& O\left(1\right)\\ \frac{{\mathit{du}}_0^{\mathrm{ex}}}{\mathit{dx}}& O\left(1\right)& O\left(\frac{1}{\epsilon }\right)\\ \frac{d^2{u}_0^{\mathrm{ex}}}{{\mathit{dx}}^2}& O\left(1\right)& O\left(\frac{1}{{\epsilon }^2}\right)\end{array}$$

This indicates that $x$ is the appropriate independent variable outside the boundary layer where $u_0^{\mathrm{ex}}$ and its derivatives are of $O\left(1\right)$ quantities. However, within the boundary layer the appropriately scaled independent variable is $s=x/\epsilon$, then

$$\frac{\mathit{du}}{\mathit{dx}}={\epsilon }^{-1}\frac{\mathit{dv}}{\mathit{ds}},\frac{d^{2}u}{{\mathit{dx}}^2}={\epsilon }^{-2}\frac{d^{2}v}{{\mathit{ds}}^2},$$

so that within the boundary layer

$$\frac{\mathit{du}}{\mathit{dx}}=O\left(1\right),\text{and}\frac{d^{2}u}{{\mathit{dx}}^2}=O\left(1\right).$$

The variable $s=x/\epsilon$ is called a stretched variable. The differential equations becomes

$$\frac{d^{2}v}{{\mathit{ds}}^2}+\frac{\mathit{dv}}{\mathit{ds}}=\epsilon +2{\epsilon }^{2}s.$$E48

We assume a boundary layer expansion, called the inner expansion of the form

$$v\left(s;\epsilon \right)\sim v_0\left(s\right)+\epsilon v_1\left(s\right)+\cdots .$$E49

The inner expansion will satisfy the boundary condition at $x=s=0$ namely $v_0\left(s=0\right)=1$ giving $v_0\left(0\right)=1$, and $v_n\left(0\right)=0,n=1,2,\cdots$. Substituting (49) into the DE (48), we obtain the following set of equations:

$$\left\{\begin{array}{ll}O\left(1\right):& \frac{d^2{v}_0}{{\mathit{ds}}^2}+\frac{{\mathit{dv}}_0}{\mathit{ds}}=0,v_0\left(0\right)=1\\ O\left(\epsilon \right):& \frac{d^2{v}_1}{{\mathit{ds}}^2}+\frac{{\mathit{dv}}_1}{\mathit{ds}}=1,v_1\left(0\right)=0\\ O\left({\epsilon }^2\right):& \frac{d^2{v}_2}{{\mathit{ds}}^2}+\frac{{\mathit{dv}}_2}{\mathit{ds}}=2s,v_1\left(0\right)=0\\ O\left({\epsilon }^n\right):& \frac{d^2{v}_n}{{\mathit{ds}}^2}+\frac{{\mathit{dv}}_n}{\mathit{ds}}=0,v_n\left(0\right)=0,n=3,4,\cdots \end{array}\right.$$E50

with solutions

$$\left\{\begin{array}{l}{v}_0=A+\left(1-A\right)e^{-s}\\ v_1=B-{\mathit{Be}}^{-s}+s\\ v_2=C-{\mathit{Ce}}^{-s}+s^2-2s\\ v_n=D_n-D_n{e}^{-s},n=3,4,\cdots \end{array}\right.$$E51

The boundary condition at $x=1$ cannot be used to determine the constants appearing in these solutions because the DEs (50) are only valid in the boundary layer. The constants in (51) are determined by matching the inner and outer expansions. We shall first restrict our attention to matching the leading order expansions $u_0$ and $v_0$. The method which we shall apply is Prandtl’s matching condition.

The leading order terms in the ‘inner’ and ‘outer’ expansions are to be matched at the ‘edge of the boundary layer’. Of course there is no precise edge of the boundary layer, we simply know that it has thickness of order $O\left(\epsilon \right)$. A plausible matching procedure would be to equate $u_0$ and $v_0$ at a value of $x$ such that the region of rapid change has passed. We might choose to equate the terms at the point $x=5\epsilon$. The leading order expansions are

$$u_0=x^2+x+2v_0=A+\left(1-A\right)e^{-s}.$$

Equating at $x=5\epsilon$ gives the following:

$$A=\frac{2+5\epsilon +25{\epsilon }^2-e^{-5}}{1-e^{-5}}.$$

If, instead we choose to match at $x=6\epsilon$, then we obtain

$$A=\frac{2+6\epsilon +36{\epsilon }^2-e^{-6}}{1-e^{-6}}.$$

These two expressions differ in the argument of the exponential and differ algebraically with $5\epsilon$ replaced by $6\epsilon$. The exponential functions are approaching transcendentally small values so that their contribution can be neglected. The algebraic difference is of $O\left(\epsilon \right)$. Thus, the arbitrariness in the decision of the point at which we choose to equate the expansions leads to a difference of $O\left(\epsilon \right)$. But we are only dealing with leading order expansions anyway. The difference between the exact solution and the leading order expansions will of $O\left(\epsilon \right)$ so that an arbitrariness in $v_0$ and $u_0$ of $O\left(\epsilon \right)$ is immaterial. Rather than choose between, for example, $5\epsilon$ and $6\epsilon$ as the value of $x$ to evaluate $u_0$ we may take the value at $x=0$, since

$$u_0\left[x=O\left(\epsilon \right)\right]=u_0\left(0\right)+O\left(\epsilon \right),$$

where the remainder is uniformly $O\left(\epsilon \right)$ since the gradient of $u_0$ is $O\left(1\right)$. For the inner expansion we are to ensure that the rapidly varying function has achieved its asymptotic value at the edge of the boundary layer. This means that the term $e^{-x/\epsilon }$ should be replaced by zero. This can be achieved by taking the limit $s\to \infty$. Thus, rather than choosing a specific point to equate the inner and outer terms er are led to the following Prandtl’s matching condition:

$$\underset{x\to 0}{lim}{u}_0\left(x\right)=\underset{s\to \infty }{lim}{v}_0\left(s\right).$$E52

The limit $s\to \infty$ may appear rather dangerous since although it certainly removes the exponential term it could lead to an algebraically unbounded term. For example, if $v_0=\mathit{As}+\left(1-A\right)e^{-s}$, then the first member would be unbounded as $s\to \infty$. This possibility can be eliminated since the inner expansion must be of a form which varies rapidly for $x=O\left(\epsilon \right)$ but not for $x=O\left(1\right)$, i.e., not for $s\to \infty$. In practice, if the boundary layer has been properly located and the correct inner variable is used then Prandtl’s matching condition is valid and elegantly avoids the need to choose an arbitrary ‘edge’ of the boundary layer.

Applying these conditions to the current example leads to

$$\underset{x\to 0}{lim}\left(x^2+x+2\right)=\underset{s\to \infty }{lim}\left[A+\left(1-A\right)e^{-s}\right],$$

which yields $A=2$. Thus the leading order terms in the expansion solutions are

$$\begin{array}{ll}\text{Outer}\text{region}:& u_0=x^2+x+2,\text{for}x=O\left(1\right)\\ \text{Inner}\text{region}:& v_0=2-e^{-x/\epsilon },\text{for}x=O\left(\epsilon \right)\end{array}$$

To prove that these are valid leading terms we consider $u^{\mathrm{ex}}$:

$$\begin{array}{ll}\text{If}x=O\left(1\right),& \text{then}{u}_0^{\mathrm{ex}}=x^2+x+2+\mathrm{T}.\mathrm{S}.\mathrm{T}.\\ \text{If}x=O\left(\epsilon \right),& \text{then}{u}_0^{\mathrm{ex}}=2-e^{-x/\epsilon }+O\left(\epsilon \right)\end{array}$$

We conclude that the matching condition has correctly predicted the leading order terms.

4.1 Composite expansion

As single composite expression for these leading order terms can be constructed using the combination

$$u_0^{\text{comp}}=u_0+v_0-u_0^{\text{match}},$$E53

where $u_0^{\text{match}}$ is given by (52). Then,

$$\begin{array}{ll}\text{for}x=O\left(1\right),v_0=u_0^{\text{match}}+\mathrm{T}.\mathrm{S}.\mathrm{T}.,& \mathrm{so}\text{that}{u}_0^{\text{comp}}=u_0^{\text{match}}+\mathrm{T}.\mathrm{S}.\mathrm{T}.\\ \text{for}x=O\left(\epsilon \right),u_0=u_0^{\text{match}}+O\left(\epsilon \right),& \mathrm{so}\text{that}{u}_0^{\text{comp}}=v_0+O\left(\epsilon \right)\end{array}$$

For the current example, $u_0^{\text{match}}=2$, so the composite expansion is

$$u_0^{\text{match}}=x^2+x+2-e^{-x/\epsilon }.$$E54

Prandtl’s matching condition can only be used for the leading order terms in the asymptotic expansions.

The outer, inner and composite expansions of the BVP (38) are presented in Figures 2 and 3 for different values of $\epsilon$. From these figures, one can easily identify the need and efficiency of the composite expansion.

Figure 2.

Outer, inner and composite expansions. (a) For ε = 0.2; (b) For ε = 0.1.

Figure 3.

Outer, inner and composite expansions. (a) For ε = 0.05; (b) For ε = 0.025.

4.2 Boundary layer location

Consider the following linear DE

$$\epsilon \frac{d^{2}u}{{\mathit{dx}}^2}+a\left(x\right)\frac{\mathit{du}}{\mathit{dx}}+b\left(x\right)u=c\left(x\right),x\in \left(x_1,x_2\right).$$E55

The following general statements can be made about the boundary layer location and the nature of the inner expansion.

Case I. If $a\left(x\right)>0$ throughout $\left(x_1,x_2\right)$, then the boundary layer will occur at $x=x_1$. The stretching transformation will be $s=\left(x-x_1\right)/\epsilon$, and the one-term inner expansion will satisfy

$$\frac{d^2{v}_0}{{\mathit{ds}}^2}+a\left(x_1\right)\frac{{\mathit{dv}}_0}{\mathit{ds}}=0.$$

The solution of this equation is

$$v_0=A+{\mathit{Be}}^{a\left(x_1\right)\left(x-x_1\right)/\epsilon },$$

where $A+B=u\left(x=x_1\right)$. The other condition to determine the constants $A$ and $B$ is obtained by matching with the value of the outer expansion at $x=x_1$.

Case II. If $a\left(x\right)<0$ throughout $\left(x_1,x_2\right)$, then the boundary layer will occur at $x=x_2$. The stretching transformation will be $s=\left(x_2-x\right)/\epsilon$, and the one-term inner expansion will involve the rapidly decaying function $e^{a\left(x_2\right)\left(x_2-x\right)/\epsilon }$.

Case III. If $a\left(x\right)$ changes sign in the interval $x_1<x<x_2$, then a boundary layer occurs at an interior point $x_0$, where $a\left(x_0\right)=0$ and boundary layers may also occur at both ends $x_1$ and $x_2$.

4.3 Boundary layer thickness and the principle of least degeneracy

The boundary layers which we have met so far have all had thickness $O\left(\epsilon \right)$. By this we mean that a variation of $O\left(\epsilon \right)$ in the independent variable will encompass the region of rapid change in the dependent variable. The associated stretched independent variable $s$, appropriate for the boundary layer is related to $x$ by a linear transformation involving division by $\epsilon$.

There are practical situations where the boundary layer thickness will be of $O\left({\epsilon }^p\right)$. This means that if the boundary layer is located at $x=x_0$, then the appropriate stretching transformation is $s=\left(x-x_0\right)/{\epsilon }^p$. More generally, the choice of the function $\delta \left(\epsilon \right)$ to use in the stretching transformation $s=\left(x-x_0\right)/\delta \left(\epsilon \right)$ is determined by the need to represent the region of rapid change correctly. We must ensure that the boundary layer solution contains rapidly varying functions. The form of the governing equation in the boundary layer region must have sufficient structure to allow such solutions.

Consider the example

$$\left\{\begin{array}{l}\epsilon \frac{d^{2}u}{{\mathit{dx}}^2}+\frac{\mathit{du}}{\mathit{dx}}+u=x,\left(0,1\right)\\ u\left(0\right)=1,u\left(1\right)=2.\end{array}\right.$$E56

Since the signs of the first and second derivatives are the same, and the boundary layer will occur at $x=0$. We are not going to assume at the outset that the boundary layer thickness is $O\left(\epsilon \right)$. Our intension is to deduce that the appropriate stretching variable is $s=x/\epsilon$.

The one-term outer expansion $u_0$ satisfies $\frac{{\mathit{du}}_0}{\mathit{dx}}+u_0=x,u_0\left(1\right)=2$ The solution is

$$u_0\left(x\right)=2e^{1-x}+x-1.$$E57

To determine the inner expansion we first wrongly assume that the boundary layer thickness is $O\left({\epsilon }^{1/2}\right)$. The stretching transformation $s=x/{\epsilon }^{1/2}$ changes the original DE (56) into the following one:

$$\frac{d^{2}v}{{\mathit{ds}}^2}+\frac{1}{{\epsilon }^{1/2}}\frac{\mathit{dv}}{\mathit{ds}}+v={\epsilon }^{1/2}s$$E58

If the appropriate stretching transformation has been used for the boundary layer then $\mathit{dv}/\mathit{ds}$ and $d^{2}v/{\mathit{ds}}^2$ will be of $O\left(1\right)$ within it. The leading order expansion $v_0$ will satisfy the dominant part of (58), i.e., the component of $O\left({\epsilon }^{-1/2}\right)$

$$\frac{{\mathit{dv}}_0}{\mathit{ds}}=0,v_0\left(0\right)=1.$$E59

The solution is $v_0\left(s\right)=1$. This of course does not have the rapidly varying behavior which we anticipate in the boundary layer. Prandtl’s matching condition cannot be satisfied since

$$\underset{x\to 0}{lim}\left(2e^{1-x}+x-1\right)=2e-1\ne \underset{s\to \infty }{lim}{v}_0\left(s\right)=1.$$

Thus, we reject the assumption of a boundary layer of thickness $O\left({\epsilon }^{1/2}\right)$.

Next, suppose that the boundary layer thickness is $O\left({\epsilon }^2\right)$ and again we will discover that this is incorrect because the corresponding inner expansion cannot be matched to the outer expansion. Proceeding with the analysis we introduce the stretching transformation $s=x/{\epsilon }^2$ which leads to the equation

$$\frac{1}{{\epsilon }^3}\frac{d^{2}v}{{\mathit{ds}}^2}+\frac{1}{{\epsilon }^2}\frac{\mathit{dv}}{\mathit{ds}}+v={\epsilon }^{2}s.$$

Again we argue that if the appropriate stretching has been used then all derivatives are of $O\left(1\right)$ so that the governing equation for the leading term in $O\left({\epsilon }^{-3}\right)$, namely

$$\frac{d^2{v}_0}{{\mathit{ds}}^2}=0,v_0\left(0\right)=1.$$E60

The solution is $v_0\left(s\right)=1+\mathit{As}$, where the constant $A$ is to be determined from matching. This solution is rapidly varying but the rapidity does not decay at the edge of the boundary layer (i.e., as $s\to \infty$). Indeed, we cannot match $v_0$ to the outer expansion because the term $\mathit{As}$ becomes arbitrarily large as $s\to \infty$.

The correct choice of stretching transformation is $s=x/\epsilon$ showing that the boundary layer thickness is $O\left(\epsilon \right)$. The boundary layer equation becomes

$$\frac{1}{\epsilon }\frac{d^{2}v}{{\mathit{ds}}^2}+\frac{1}{\epsilon }\frac{\mathit{dv}}{\mathit{ds}}+v=\mathit{\epsilon s}.$$

The dominant equation satisfied by $v_0$ is $O\left(1/\epsilon \right)$, namely

$$\frac{d^2{v}_0}{{\mathit{ds}}^2}+\frac{{\mathit{dv}}_0}{\mathit{ds}}=0,v_0\left(0\right)=1.$$E61

The solution is $v_0\left(s\right)=1-A+{\mathit{Ae}}^{-s}$. The last member provides the necessary rapid decay away from the point $x=s=0$. Prandtl’s matching condition requirest

$$\underset{x\to 0}{lim}\left(2e^{1-x}+x-1\right)=\underset{s\to \infty }{lim}\left(1-A+{\mathit{Ae}}^{-s}\right),$$

which leads to $A=2-2e$, and

$$v_0\left(x\right)=2e-1+2\left(1-e\right)e^{-x/\epsilon }.$$

The one-term composite expansion is

$$u^{\text{comp}}=\left(2e^{1-x}+x-1\right)+\left(2e-1\right)+2\left(1-e\right)e^{-x/\epsilon }-\left(2e-1\right).$$E62

The leading order boundary layer equation associated with the stretching transformation $s=x/\epsilon$, (61) involves more terms than (59), associated with $s=x/{\epsilon }^{1/2}$, and (60) associated with $s=x/{\epsilon }^2$. The extra term in (61) allows sufficient structure in the solution to produce the required boundary layer behavior. An aid for choosing the boundary layer thickness is to seek a stretching transformation which retains the largest number of terms in the dominant equation governing $v_0$. This referred to as the principle of least degeneracy by Van Dyke.

The composite expansion (62) can be verified by comparing with the exact solution of (56). The general solution of (56) is

$$u^{\mathrm{ex}}=C_1{e}^{m_{1}x}+C_2{e}^{m_{2}x}+\left(x-1\right),$$

where

$$m_1=\frac{-1+\sqrt{1-4\epsilon }}{2\epsilon },m_2=\frac{-1-\sqrt{1-4\epsilon }}{2\epsilon }.$$

We expand $\sqrt{1-4\epsilon }$ using the binomial series, $\sqrt{1-4\epsilon }=1-2\epsilon +O\left({\epsilon }^2\right)$, then

$$m_1=-1+O\left(\epsilon \right),\text{and}{m}_2=-\frac{1}{\epsilon }+1+O\left(\epsilon \right),$$

so that

$$u^{\mathrm{ex}}=C_1{e}^{-x}+C_2{e}^{-x/\epsilon }\cdot e^x+\left(x-1\right)+O\left(\epsilon \right).$$E63

Using the boundary conditions and by neglecting the transcendentally small term $e^{-1/\epsilon }$, we have $C_1=2e,C_2=2\left(1-e\right)$. Then, (63) becomes

$$u^{\mathrm{ex}}=2e^{1-x}+2\left(1-e\right)e^{-x/\epsilon }\cdot e^x+\left(x-1\right)+O\left(\epsilon \right).$$E64

There is an apparent discrepancy between (64) and the composite expansion (62) in the coefficient of the $e^{-x/\epsilon }$ term. There is an extra term only contributes in the boundary layer where $x=O\left(\epsilon \right)$ so that the coefficient $e^x$ may to leading order, be replaced by unity. Thus, the leading order composite expansion and the leading order term in the exact solution are in complete agreement.

4.4 Boundary layer of thickness of $O\left(\sqrt{\epsilon }\right)$

Consider the following two-point BVP:

$$\left\{\begin{array}{l}\epsilon \frac{d^{2}u}{{\mathit{dx}}^2}+x^2\frac{\mathit{du}}{\mathit{dx}}-u=0,\left(0,1\right)\\ u\left(0\right)=1,u\left(1\right)=2.\end{array}\right.$$E65

We seek a one-term composite expansion for the above BVP. We will tentatively assume that a boundary layer occurs at $x=0$ although the vanishing of the coefficient of the first derivative suggests the possibility of nonstandard behavior.

The one term outer expansion satisfies

$$x^2\frac{{\mathit{du}}_0}{\mathit{dx}}-u_0=0,u_0\left(1\right)=2.$$

Its exact solution is $u_0\left(x\right)=2e^{\left(1-1/x\right)}$.

Let us assume that the boundary layer thickness is of $O\left({\epsilon }^p\right)$, where $p$ is to be determined from the principle of least degeneracy. The stretched variable is $s=x/{\epsilon }^p$, and (65) becomes

$${\epsilon }^{1-2p}\frac{d^{2}v}{{\mathit{ds}}^2}+{\epsilon }^p{s}^2\frac{\mathit{dv}}{\mathit{ds}}-v=0.$$

The second-term is always dominated by the third, so the principle of degeneracy requires the first term to be of the same order as the third term (i.e., $O\left(1\right)$). Thus, $p=1/2$, and the one-term inner expansion satisfies

$$\frac{d^2{v}_0}{{\mathit{ds}}^2}+\frac{{\mathit{dv}}_0}{\mathit{ds}}=0,v_0\left(0\right)=1.$$

The solution of the above problem is $v_0\left(s\right)={\mathit{Ae}}^s+\left(1-A\right)e^{-s}$. Prandtl’s matching condition requires

$$\underset{x\to 0}{lim}2e^{\left(1-1/x\right)}=\underset{s\to \infty }{lim}\left[{\mathit{Ae}}^s+\left(1-A\right)e^{-s}\right]$$

which yields $A=0$. This example is rather special in that $A$ will be zero for all boundary conditions.

The on-term composite expansion is

$$u_0^{\text{comp}}=2e^{\left(1-x\right)}+e^{-x/\sqrt{\epsilon }}.$$

We conclude this example with the observation that a choice for the value of the index $p$ other than $p=1/2$ leads to boundary layer equations with insufficient structure to generate the required rapidly decaying behavior.

Thus, if $p>1/2$, the dominant equation becomes

$$\frac{d^2{v}_0}{{\mathit{ds}}^2}=0,v_0\left(0\right)=1,$$

which gives $v_0\left(s\right)=1+\mathit{As}$. It is obvious that Prandtl’s matching condition cannot be used to determine $A$. Whereas, if $p<1/2$ the dominant equation degenerates to $v_0\left(s\right)=0$ which does not satisfy the boundary condition at $s=0$.

4.5 Interior layer

Consider the BVP:

$$\left\{\begin{array}{l}\epsilon \frac{d^{2}u}{{\mathit{dx}}^2}+x\frac{\mathit{du}}{\mathit{dx}}+\mathit{xu}=0,\left(-1,1\right)\\ u\left(-1\right)=e,u\left(1\right)=2e^{-1}.\end{array}\right.$$E66

The coefficient of the first derivative (convective term) is positive in $\left(0,1\right)$ which indicates the occurrence of a boundary layer at the left hand limit $x=0$. While the corresponding coefficient is negative in the range $-1<x<0$ indicates a boundary layer located at the right-hand limit which is again is $x=0$. Thus, we are led to expr = ect two outer expansions for positive and negative $x$ respectively and an inner expansion in the interior layer located at $x=0$. We denote the leading term in the outer expansion for positive $x$ by $u_0^{+}$, it satisfies

$$\frac{{\mathit{du}}_0^{+}}{\mathit{dx}}+u_0^{+}=0,u_0^{+}\left(1\right)=2e^{-1}$$E67

with the solution $u_0^{+}\left(x\right)=2e^{-x}$.

The outer expansion for negative $x$, $u_0^{-}$ satisfies

$$\frac{{\mathit{du}}_0^{-}}{\mathit{dx}}+u_0^{-}=0,u_0^{+}\left(-1\right)=e$$E68

with the solution $u_0^{-}\left(x\right)=e^{-x}$.

We suppose the boundary layer at $x=0$ has thickness $O\left({\epsilon }^p\right)$ and determine the index $p$ using the principle of least degeneracy. Let $s=x/{\epsilon }^p$ so that the DE becomes

$${\epsilon }^{1-2p}\frac{d^{2}v}{{\mathit{ds}}^2}+s\frac{\mathit{dv}}{\mathit{ds}}+{\epsilon }^p\mathit{sv}=0.$$

The third term is dominated by the second term. The first term has the same order as the second term if $p=1/2$. For this choice of $p$ the leading term of the inner expansion $v_0$ satisfies

$$\frac{d^2{v}_0}{{\mathit{ds}}^2}+s\frac{{\mathit{dv}}_0}{\mathit{ds}}=0.$$

Its solution can be given by

$$v_0\left(s\right)=Berf\left(s/\sqrt{2}\right)+v_0\left(0\right),$$

Prandtl;s matching condition applied to the region $x>0$ is

$$\underset{s\to +\infty }{lim}{v}_0\left(s\right)=\underset{x\to 0^{+}}{lim}{u}_0^{+}\left(x\right)$$

and corresponding for $x<0$, we have

$$\underset{s\to -\infty }{lim}{v}_0\left(s\right)=\underset{x\to 0^{-}}{lim}{u}_0^{-}\left(x\right)$$

Using the limiting values $\mathit{erf}\left(\pm \infty \right)=\pm 1$ yields $v_0\left(0\right)=1.5$ and $B=0.5$. The leading order terms over the whole region are

$$\begin{array}{ll}{u}_0^{+}\left(x\right)=2e^{-x},& x>O\left(\sqrt{\epsilon }\right)\\ v_0=0.5erf\left(x/\sqrt{2\epsilon }\right)+1.5,& x=O\left(\sqrt{\epsilon }\right)\\ u_0^{-}\left(x\right)=e^{-x},& x<-O\left(\sqrt{\epsilon }\right)\end{array}$$

A composite expansion cannot be formed in the standard way when there is more than one outer solution. However, the behavior of $v_0$ for $\mid x\mid >O\left(\sqrt{\epsilon }\right)$ is as follows:

$$\begin{array}{l}{v}_0\left[x>O\left(\sqrt{\epsilon }\right)\right]=0.5+1.5+\mathrm{T}.\mathrm{S}.\mathrm{T}\\ v_0\left[x<-O\left(\sqrt{\epsilon }\right)\right]=-0.5+1.5+\mathrm{T}.\mathrm{S}.\mathrm{T}\end{array}$$

Utilizing this enables a uniformly valid one-term composite expansion to be constructed which yields the correct coefficient of $e^{-x}$ outside the boundary layer and the correct leading order behavior within the boundary layer. It is

$$u_0^{\text{comp}}=\left[0.5erf\left(x/\sqrt{2\epsilon }\right)+1.5\right]e^{-x}.$$

4.6 Nonlinear differential equation

Consider the following semilinear

$$\left\{\begin{array}{l}\epsilon \frac{d^{2}u}{{\mathit{dx}}^2}+\frac{\mathit{du}}{\mathit{dx}}+u^2=0,\left(0,1\right)\\ u\left(0\right)=2,u\left(1\right)=1/2.\end{array}\right.$$E69

The coefficient of the first and second order derivatives have the same sign, so the boundary layer will occur at the left boundary $x=0$. The one-term outer expansion satisfies

$$\frac{{\mathit{du}}_0}{\mathit{dx}}+u_0^2=0,u_0\left(1\right)=1/2,$$

and the solution is $u_0\left(x\right)=1/\left(1+x\right)$. The stretching transformation for the inner region will be $s=x/\epsilon$ and therefore, the inner expansion satisfies

$$\frac{d^{2}v}{{\mathit{ds}}^2}+\frac{\mathit{dv}}{\mathit{ds}}+\epsilon v^2=0,v\left(0\right)=2.$$

The one-term inner expansion $v_0$ satisfies the dominant part of this equation, i.e.,

$$\frac{d^2{v}_0}{{\mathit{ds}}^2}+\frac{{\mathit{dv}}_0}{\mathit{ds}}=0,v_0\left(0\right)=2,$$

which gives $v_0\left(s\right)=A+\left(2-A\right)e^{-s}$. Prandtl’s matching condition yields $A=1$, and the composite one-term uniformly valid expansion is

$$u_0^{\text{comp}}=\frac{1}{1+x}+e^{-x/\epsilon }.$$

Next, consider the quasilinear problem

$$\left\{\begin{array}{l}\epsilon \frac{d^{2}u}{{\mathit{dx}}^2}+2u\frac{\mathit{du}}{\mathit{dx}}-4u=0,\left(0,1\right)\\ u\left(0\right)=0,u\left(1\right)=4.\end{array}\right.$$E70

The nonlinearity is associated with the first derivative term. The location of the boundary layer depends on the relative sign of the first and second derivative coefficients. If we assume that the dependent variable is nonnegative throughout the interval $0<x<1$, then the boundary layer will occur at $x=0$. The one-term outer expansion satisfies

$$2u_0\frac{{\mathit{du}}_0}{\mathit{dx}}-4u_0=0,u_0\left(1\right)=4,$$

with the solution $u_0\left(x\right)=2x+2$.

Assuming that the boundary layer thickness is $O\left(\epsilon \right)$, therefore, the dominant-order equation for the one-term inner expansion becomes

$$\frac{d^2{v}_0}{{\mathit{ds}}^2}+2v_0\frac{{\mathit{dv}}_0}{\mathit{ds}}=0,v_0\left(0\right)=0.$$

Its solution is $v_0\left(s\right)=atanh\left(\mathit{as}\right)$. Prandtl’s matching condition yields $a=2$. Thus, $v_0\left(s\right)=2tanh\left(2s\right)$, and the uniformly valid one-term composite expansion is

$$u_0^{\text{comp}}=2x+2+2tanh\left(2s\right)-2.$$

Application of perturbation techniques to partial differential equations, and other types of problems can be seen in the books [5, 6].

Conflict of interest

The authors declare no conflict of interest.

Notes/thanks/other declarations

I owe a great debt to my mentor Prof. S Natesan, Department of Mathematics, IIT Guwahati, who introduced me to this topic. The chapter was discussed during my stay at IIT Guwahati.