Perturbation analysis generally deals with an unsolvable problem by treating it as a perturbation from a solvable problem. The distinction between regular and singular that in a singular problem there is a qualitative difference in the natures of the solution to the solvable problem and the unsolvable problem.

Asymptotic Orders of Functions:
Big O and Little o

The functions, Big O and Little o, are functions from the set of functions to the set of equivalence classes. Their definitions are as follows:

• Big O: A function f belongs to the same equivalence class as a function g, i.e.,

  f ∈ O(g) for some interval I if there exists a constant k
  such that |f(x)|<k|g(x)| for all x in I.

• Little o: A function f belongs to the same equivalence class as a function g, i.e.,

  f ∈ o(g) if for any k>0 there exists a an interval I_k such that |f(x)|<k|g(x)| for all x in I_k.

The more convenient tests are:

• O(f) = O(g) for some interval I if |f/g| is finite and positive on the interval I
• o(f) = o(g) if lim |f/g| = 0.

Perturbation Theory

Example: Particle Motion
Near a Planet's Surface

Let R be the radius of the planet and y be the distance above its surface. The distance to the center of mass of the planet is thus R+y and hence the force upon a particle of mass m is

-GmM/(R+y)²,

where M is the mass of the planet and G is the gravitational constant.

By Newton's First Law

m(d²/dt²)((R+y)r^) =
-GmM/(R+y)²r^,

where r^ is a unit vector directed radially away from the center of the planet.

The mass of the particle cancels. For convenience let GM/R² equal g and note that R is a constant so the equation for the motion of the particle reduces to:

d²/dt² ((R+y)r^) =
d²/dt²(yr^) =
-[gR²/(R+y)²]r^.

Since all motion takes place in the radial direction the r^ can be dropped from the analysis.

A particle at the surface of the planet with an initial vertical velocity of V and subjected to an acceleration of -g will rise to an maximum altitude of V²/2g. For convenience later let us define the characteristic length h_m as twice this height; i.e.,

h_m = V²/g.

The time taken for the velocity to be reduced from the initial V to zero is

t_m = V/g.

Define new nondimensional variables

x = y/h_m
and s = t/t_m.

The equation of motion is then:

[h_m/t_m²]d²x/ds² = -gR²/[h_mx +R]².

Since [h_m/t_m²] = g the equation of motion can be reduced to:

d²x/ds² = -[(h_m/R)x + 1]^-2.

The quantity h_m/R is small and nondimensional; let us call it ε. Now the equation of motion is:

d²x/ds² = -[1 + εx]^-2.

the initial conditions are now that:

x(0) = 0 and dx(0)/ds = 1.

The function on the righthand side of the equation can be expanded into a power series in εx. Note that

[1+z]^-1 = 1 - z + z² - z³ + ....

If we take the derivative of each side of this equation we get

-[1+z]^-2 = -1 + 2z - 3z² +....

Thus the equation of motion takes the form:

d²x/ds² =
-1 + 2(εx) - 3(εx)² +....

If x(s, ε) has an expansion of the form:

x(s, ε) = Σx_i(s) εⁱ
for i = 0 to ∞.

Then it must be that:

x₀" + x₁"ε + x₂"(ε)² + ^{. . .} = -1 +
2ε[x₀ + x₁ε + x₂(ε)² + ^{. . .}]
- 3ε²[x₀ + x₁ε + x₂(ε)² + ^{. . .}]
+ 4ε³[x₀ + x₁ε + x₂(ε)² + ^{. . .}] + ^{. . .}

The coefficient of each power of ε must be the same on each side of the above equation. Therefore:

ε⁰:   x₀" = -1;
ε¹:   x₁" = 2x₀;
ε²:   x₁" = 2x₁ - 3x₀²;

The initial condition of x(0)=0 implies that

x_i(0) = 0 for all i

and the initial condition of x'(0)=1 implies that

x₀(0) = 1 and
x_i(0) = 0 for all i>0

These initial value problems can be solved successively. The results are:

x₀ = -(1/2)s² + s;
x₁ = - (1/12)s⁴ + (1/3)t³
x₂ = - (11/360)s⁶ + (11/60)s⁵
- (1/4)s⁴

Therefore we have:

x(s, ε) = (-(1/2)s² + s) +
ε(- (1/12)s⁴ + (1/3)s³) +
ε²(- (11/360)s⁶ + (11/60)s⁵
- (1/4)s⁴) + . . .

The analytic solution to:

d²x/ds² = -[1 + εx]^-2
with the initial conditions that
x(0) = 0 and dx(0)/ds = 1

can be achieved conveniently by letting [1 + εx]=z and s=rε so the equation of motion becomes:

z²d²z/dr²=-1

and the initial conditions are now that z(0)=1 and dz/dr(0)=1. Now the equation can be integrated using the technique of integration by parts.

Thus,

z²dz/dr -1 - ∫ 2z(dz/dr)dr = -r

or

z²dz/dr -1 + z² -1 = -r

Regular Perturbation

Consider a Taylor series expansion such as

exp(-εt) -1 = -εt + (-εt)²/2! + . . .

Such a series is said to be uniform or asymptotically regular if

O((εt)²) < O(εt).

But if

O(εt) = O(1)

then the series ordering breaks down and the series is said to be nonuniform and cannot be used for approximations. The region where the series ordering breaks down is called a boundary layer.

Asymptotic Series

Consider the differential equation:

dy/dx + y = 1/x

If one looks for a solution to this equation in the form of series in x that is valid as x→∞ one obtains:

Y(x) = 1/x + 1/x² + 2!/x³ + ... + (n-1)!/xⁿ + . . .

As can be seen from the ratio test this series does not converge.

Perturbation Analysis of a
Damped Linear Oscillator

Consider

m(d²Y/dT²) + b(dY/dT) + kY = I₀δ(T),
Y(0-) = 0, (dY/dT)(0-) = 0
where δ(T) is the Dirac delta function.

An integration over 0- to 0+ produces the result:

m(d²Y/dT²) + b(dY/dT) + kY = 0
Y(0+ = 0, (dY/dT)(0+)= I₀/m

This equation can be put into nondimensional form by making the substitutions

t = [k/m]^1/2T,
ε = b/(2[km]^1/2
y = Y[km]^1/2/I₀

The nondimensional form of the equation is then:

(d²y/dt²) + 2ε(dy/dt) + y = 0
y(0+ = 0, (dy/dt)(0+)= 1

The exact solution of this system is known and it is:

y(t) =
exp(-εt)sin([1-εt]^1/2t)/ [1-εt]^1/2

This exact solution is useful for comparison with the approximation derived using perturbation theory.

Regular perturbation theory makes the assumption that the solution can be expression in a series of the form:

y(t,ε) = f₀(t) + εf₁(t) + ε² f₂(t) + . . .

When this series form is substituted into the differential equation the result is:

f"₀(t) + εf"₁(t) + ε² f"₂(t) + . . .
     2εf'₀(t) + 2ε²f'₁(t) +
                                             2ε³ f'₂(t) + . . .
f₀(t) + εf₁(t) + ε² f₂(t) + . . .
= 0

Equating the coefficients of the poweres of ε to zero gives:

ε⁰: f"₀(t) + f₀(t)= 0;
ε¹: f"₁(t) + 2f'₀(t) + f₁(t) = 0

the initial conditions are that:

f₀(0) = 0, f'₀(0) = 1
f₁(0) = 0, f'₁(0) = 0

The solutions are:

f₀(t) = sin(t), f₁(0) = -t·sin(t),

Thus the approximate solution is:

y(t,ε) (1 - εt)sin(t)

This approximation is good only for εt being small. For sufficiently large t the approximate solution differs very drastically from the exact solution.

Singular Perturbation Analysis

The previous example illustrates regular perturbation ananlysis where the damping coefficient is a vanishingly small parameter. But the damped linear oscillator also provides a case of singular perturbation analysis when the mass is a vanishingly small parameter.

The difference between the two cases is that mass is the coefficient of the second derivative in the equation so its vanishing reduces the order of the differential equation from two to one. This reduction in order makes it impossible to simultaneously satisfy both initial conditions. This problem leads to boundary layer phenomena.

First consider the limiting case of m=0. The original system is then

b(dY/dT) + kY = I₀δ, Y(0-)=0.

Integration across 0- to 0+ gives

bY(0+) = I₀

Now the system for t>0 is

b(dY/dT) + kY = 0 and Y(0+) = I₀/b.

The solution to this initial value problem is:

Y(T) = (I₀/b)[1 - exp(-bt/m)].

The original problem can be transformed into a nondimensional form by letting

y =(b/I₀)Y,
t =(k/b)T,
and ε = mk/b².

The result of these transformations is:

ε(d²y/dt²) +(dy/dt) + y = 0
y(0) = 0, y'(0) = 1/ε.

The initial condition of y'(0) = 1/ε precludes the series used in regular perturbation analysis because such an initial condition cannot be satisfied by that series. An alternatet series is:

y(t, ε) =
ν₀(ε)h₀(t) + ν₁(ε)h₁(t) + ν₂(ε)h₂(t) +
where ν_n+1 is o(ν_n)
so ν_n+1/ν_n → 0 as ε →0.

Substitution of this form into the nondimensional differential equation results in

          εν₀(ε)h"₀(t) + εν₁(ε)h"₁(t) + εν₂(ε)h"₂(t) +
ν₀(ε)h'₀(t) + ν₁(ε)h'₁(t) +
ν₂(ε)h'₂(t) +
ν₀(ε)h₀(t) + ν₁(ε)h₁(t) +
ν₂(ε)h₂(t) = 0

Since the ν_n(ε) are of different orders of magnitude we get the conditions

h'₀(t) + h₀(t) = 0
and
h'₁(t) + h₁(t) + h"₀(t) = 0
if h₀(t) = εν₀(ε) = ν₁(ε)
or
h'₁(t) + h₁(t) = 0
if εν ₀(ε)/ ν₁(ε) →0 as ε →0.

and likewise for n>1.

The initial conditions for the h_n(t) functions are uncertain.

Leaving the initial conditions unspecified, the solutions for h₀(t) and h₁(t) are:

h₀(t) = A₀exp(-t)
and
h₁(t) = A₁exp(-t) - A₀t exp(-t)
if εν₀ /ν₁=1
or
h₁(t) = A₁exp(-t)
if εν₀ /ν₁=0

This is the outer solution.

Inner Solution

An inner solution is a solution which is valid in the vicinity of the boundary or intitial condition. In this an inner solution is one that is valid for small values of time.

Consider a general transformation of the time variable; i.e.,

τ = t/ φ (ε)
where φ(ε) →0 as ε→0

In terms of τ the differential equation is:

(ε/φ²(ε)(d²y/dτ²) + (1/φ(ε))(dy/dτ) + y = 0.

Now consider the various possibilities for φ (ε).

Case I: ord(φ) < ord(ε)

In this case (ε/φ) → 0 as ε → 0.

Since

(d²y/dτ²) + (φ(ε)/ε)(dy/dτ) + (φ²/ε)y = 0
as ε → 0
the equation goes to
(d²y/dτ²) = 0.

Since this is a second order equation the initial conditions of y(0)= 0 and y'(0) = 1/ε can be satisfied.

Case I: φ(ε) = ε

In this case the differential equation, in terms of τ, becomes:

(d²y/dτ²) + (dy/dτ) + εy = 0
and hence in the limit
as ε → 0
the equation goes to
(d²y/dτ²) + (dy/dτ) = 0
which has the solution
y(τ) = 1 - exp(-τ)
= 1 - exp(-t/ε)

(to be continued)