Perturbation Methods
Perturbation Methods
Perturbation Methods
GM01
Dr. Helen J. Wilson
Autumn Term 2008
2
1 Introduction
3
1.3 A real research example
This comes from my own research2 . I will not present the equations or the working here:
but the problem in question is the stability of a polymer extrusion flow. The parameter
varied is wavelength: and for both very long waves (wavenumber k ≪ 1) and very short
waves (k −1 ≪ 1) the system is much simplified. The long-wave case, in particular, gives
very good insight into the physics of the problem.
If we look at the plot of growth rate of the instability against wavenumber (inverse
wavelength):
we can see good agreement between the perturbation method solutions (the dotted lines)
and the numerical calculations (solid curve): this kind of agreement gives confidence in
the numerics in the middle region where perturbation methods can’t help.
1.4 References
The principal reference for this course is the book by Hinch. In particular, many of the
examples and exercises are taken from it. The others are also very good – your choice is
really a question of style preference. Also see http://www.ucl.ac.uk/~ucahhwi/GM01/.
4
2 Regular perturbation expansions
x2 + εx − 1 = 0. (1)
If ε is small, ε ≪ 1, we can then expand for small ε (using the Taylor series for the square
root3 ) to have
1 1 1 4
⎧
⎨ 1 − ε + ε2 −
⎪ ε + ···
x= 2 8 128
⎩ −1 − 1 ε − 1 ε2 + 1 ε4 + · · ·
⎪
2 8 128
x = −1 + εx1 + ε2 x2 + ε3 x3 + · · ·
3
Actually, in this case the expansion is valid for |ε| < 2; but we will usually require ε ≪ 1
4
This is a massive assumption. If it is true, then a regular perturbation expansion will do the job: if
not, we need a singular expansion. We’ll see them soon enough.
5
If we substitute this into equation (1) we get:
We equate powers of ε:
ε0 : 1−1 =0
ε1 : −2x1 − 1 = 0 , x1 = − 12
ε2 : −2x2 + x21 + x1 = 0 , x2 = − 18
ε3 : −2x3 + 2x1 x2 + x2 = 0 , x3 = 0
Note that the equation at ε0 was automatically satisfied. This will always happen if we
have solved the ε = 0 equation correctly.
The series we have found is
1 1
x = −1 − ε − ε2 + O(ε4 ).
2 8
This matches with the exact solution. Remember the notation O(ε4 ): this means that
the missing terms from the equation tend to zero at least as fast as ε4 , as ε → 0.
df (x)
+ f (x) − εf 2 (x) = 0, f (0) = 2. (2)
dx
We can’t solve this directly. However, let’s look at ε = 0:
df (x)
+ f (x) = 0, f (0) = 2,
dx
which has solution
f (x) = 2e−x .
So in the same way as for the algebraic equation, let us try a regular perturbation
expansion in ε:
f = 2e−x + εf1 (x) + ε2 f2 (x) + ε3 f3 (x) + · · ·
where in order to satisfy the initial condition f (0) = 2, we will have f1 (0) = f2 (0) =
f3 (0) = · · · = 0. Substituting into (2) gives
6
and we can collect powers of ε:
ε0 : −2e−x + 2e−x =0
ε1 : f1 (x) + f1 (x) − 4e−2x
′
=0
ε2 : f2′ (x) + f2 (x) − 4e−x f1 (x) =0
ε3 : f3′ (x) + f3 (x) − f12 (x) − 4e−x f2 (x) =0
has solution
f1 (x) = −4e−2x + c1 e−x
and the boundary condition f1 (0) = 0 gives c1 = 4:
f2′ (x) + f2 (x) = 4e−x f1 (x) =⇒ f2′ (x) + f2 (x) = 16e−x (e−x − e−2x )
with solution
f2 (x) = 8(−2e−2x + e−3x ) + c2 e−x
and the boundary condition f2 (0) = 0 gives c2 = 8:
which becomes
f3′ (x) + f3 (x) = [4(e−x − e−2x )]2 + 4e−x [8(e−x − 2e−2x + e−3x )]
= 48(e−2x − 2e−3x + e−4x )
7
This is an example of a case where carrying out a perturbation expansion can give us an
insight into the full solution. Notice that, for the terms we have calculated,
fn (x) = 2n+1 e−x (1 − e−x )n ,
suggesting a guessed full solution
∞ ∞
( ( 2e−x
f (x) = εn 2n+1 e−x (1 − e−x )n = 2e−x [2ε(1 − e−x )]n = .
n=0 n=0
1 − 2ε(1 − e−x )
We can check this solution, and it is indeed the correct solution to the ODE of equa-
tion (2).
Exercise 1: Find the first three terms of an expansion for each root of the following
equation:
x3 − (2 − ε)x2 − x + 2 + ε = 0.
8
3 Rescaling
Now let us look at the true solution to see what’s gone wrong.
√
−1 ± 1 + 4ε
x=
2ε
As ε → 0, the leading-order terms of the two roots are
1
x = 1 + O(ε); and − + O(1).
ε
The first of these is amenable to the simplistic approach; we haven’t seen the second root
because it → ∞ as ε → 0.
9
and collecting powers of ε gives:
ε−1 : x2−1 + x−1 = 0 ; x−1 = 0 , −1
0
ε : 2x−1 x0 + x0 − 1 = 0 ; x0 = 1 , −1
1 2
ε : x0 + 2x−1 x1 + x1 = 0 ; x1 = −1 , 1
Note that we can now get the expansions for both of the roots using the same method.
εx2 + x − 1 = 0.
We put in the new form, and then look at the different possible values of δ.
LHS = εδ 2 X 2 + δX − 1 = 0
δ≪1 LHS = small + small − 1 = 0 ℵ
δ=1 LHS = small + X − 1 = 0 regular root
1 LHS
1≪δ≪ = small + X − small = 0 ℵ
ε δ
1 LHS
δ= = X2 + X − small = 0 singular root
ε δ
1 LHS
δ≫ = X2 + small − small = 0 ℵ
ε εδ 2
Note that X = 0 cannot be a solution because we require X = ord(1). The scalings that
work (in this case δ = 1 and δ = ε−1 ) are called distinguished scalings.
ε3 x3 + x2 + 2x + ε = 0.
and find the first two nonzero terms in the expansion of each root.
εx3 + x2 + (2 − ε)x + 1 = 0.
and find the first two terms in the expansion of each root. [Hint: you may find an exact
root – then there will be no more terms in the expansion.]
10
3.3 Non-integral powers
x = 1 + εx1 + ε2 x2 + · · ·
Substituting in gives
then we get
1 + 2ε1/2 x1/2 + ε(x21/2 + 2x1 ) + ···
− ε + ···
− 2 − 2ε1/2 x1/2 − 2εx1 + ···
+ 1 = 0
2x1/2 − 2x1/2 = 0
giving two solutions x1/2 = ±1. Both of these are valid and will lead to valid expansions
if we continue.
We could have predicted that there would be trouble when we found the double root:
near a quadratic zero of a function, a change of order ε1/2 in x is needed to change the
11
function value by ε:
in which xi is strictly order 1 as ε → 0 (i.e. tends neither to zero nor infinity) and the
series of functions δi (ε) has δ0 (ε) ≫ δ1 (ε) ≫ δ2 (ε) · · · for ε ≪ 1.
Then at each order we look for a distinguished scaling. Let us work through an example:
√ ) π* 1 1
2 sin x + − 1 − x + x2 = − ε.
4 2 6
In this case there is a solution near x = 0, which we will investigate.
First let us sort out the trigonometric term, expanding it as a Taylor series about x = 0:
√ ) π* √ + )π * ) π *,
2 sin x + = 2 sin x cos + cos x sin =
-4 .4 4
√ 1 1 x2 x3 x4 x5
2 √ sin x + √ cos x = sin x + cos x = 1 + x − − + + + ···
2 2 2! 3! 4! 5!
x3 x4 x5 1
− + + + O(x6 ) = − ε.
6 24 120 6
x4 x5
x3 − − + O(x6 ) = ε.
4 20
12
We pose a series
x = x0 δ0 (ε) + x1 δ1 (ε) + · · ·
and substitute it. The leading term on the left hand side is x30 δ03 , and on the right hand
side is ε. So we set δ0 = ε1/3 and x0 = 1.
Now we have
x = ε1/3 + x1 δ1 (ε) + · · ·
which we substitute into the governing equation. Remembering that δ1 ≪ ε1/3 and
keeping terms up to order ε2/3 δ1 and ε4/3 (neglecting only terms which are guaranteed to
be smaller than one of these), we have
ε4/3
3x1 ε2/3 δ1 − = 0.
4
To make this work, we need δ1 = ε2/3 and then x1 = 1/12.
The first two terms of the solution are:
1 2/3
x = ε1/3 + ε + ···
12
εx4 − x2 − x + 2 = 0.
3.5 Logarithms
HERE BE DRAGON S
There is a worse case than fractional power of ε: logarithms. These are beyond the scope
of this course. If you ever have a problem in which a quantity like ln (1/ε) appears to be
important, either give up or go to a textbook.
13
14
4 Scalings with differential equations
f = A0 exp [−x/ε].
This gives us the clue that what we should have done was change to a stretched variable
z = x/ε.
Let us ignore the full solution and simply make that substitution in our governing equa-
tion. Note that df /dx = df /dz dz/dx = ε−1 df /dz.
df df
εε−1 +f =0 + f = 0.
dz dz
Now the two terms balance: that is, they are the same order in ε. Clearly the solution
to this equation is now A0 exp [−z] and we have found the result.
This is a general principle. For a polynomial, we look for a distinguished scaling of the
quantity we are trying to find. For a differential equation, we look for a stretched version
of the independent variable.
The process is very similar to that for a polynomial. We use a trial scaling δ and set
x = a + δ(ε)X.
Then we vary δ, looking for values at which the two largest terms in the scaled equation
balance.
Let’s work through the process for the following equation:
d2 f df
ε 2
+ − f = 0.
dx dx
Again, we note that if x = a + δX then d/dx = d/dX dX/dx = δ −1 d/dX. We substitute
in these scalings, and then look at the different possible values of δ.
LHS = εδ −2 d2 f /dX 2 + δ −1 df /dX − f = 0
−1 2
δ≪ε ε δ LHS = f ′′ + small − small = 0 ℵ
δ=ε εLHS = f ′′ + f′ − small = 0
ε≪δ≪1 δLHS = small + f′ − small = 0 ℵ
δ=1 LHS = small + f′ − f = 0
δ≫1 LHS = small + small − f = 0 ℵ
15
Here the two distinguished stretches are δ = ε and δ = 1.
For δ = 1 we can treat this as a regular perturbation expansion:
For δ = ε we use our new variable X = ε−1 (x − a) and work with the new governing
equation:
d2 f df
2
+ − εf = 0
dX dX
Now we try a regular perturbation expansion:
f = f0 + εf1 + ε2 f2 + · · ·
and so on. Of course, without boundary conditions to apply, this process spawns large
numbers of unknown constants. Rescaling to our original variable completes the process:
x−a /
-" # .
−(x−a)/ε −(x−a)/ε
0 −(x−a)/ε
f (x) ∼ A0 + B0 e +ε A0 − B0 e + A1 + B1 e + ···
ε
Note that this expansion is only valid where X = (x − a)/ε is order 1: that is, for x close
to the (unknown) value a.
Exercise 8: Find the distinguished stretches for the following differential equation:
d3 f d2 f df
ε3 3
+ ε 2
+ +f =0
dx dx dx
Find the leading-order term of each solution.
16
4.2 Nonlinear differential equations: scaling and stretching
Recall that for a linear differential equation, if f is a solution then so is Cf for any
constant C. So if f (x; ε) is a solution as an asymptotic expansion, then Cf is a valid
asymptotic solution even if C is a function of ε.
The same is not true of nonlinear differential equations. Suppose we are looking at the
equation:
d2 f df
2
+ εf (x) + f 2 (x) = 0
dx dx
There are two different types of scaling we can apply: we can scale f , or we can stretch
x. To get all valid scalings we need to do both of these at once.
Let us take f = εα F where F is strictly ord(1), and x = a + εβ z with z also strictly
ord(1). Then a derivative scales like d/dx ∼ ε−β d/dz and we can look at the scalings of
all our terms:
d2 f df
2
+ εf (x) + f 2 (x) = 0
dx dx
εα ε−2β εε2α ε−β ε2α
As always with three terms in the equation, there are three possible balances.
• For terms I and II to balance, we need α−2β = 2α+1−β. This gives α+β+1 = 0, so
that terms I and II scale as ε2+3α , and term III scales as ε2α . We need the balancing
terms to dominate, so we also need 2α > 2 + 3α which gives α < −2.
• For terms I and III to balance, we need α − 2β = 2α. This gives α = −2β, so
that terms I and III scale as ε2α and term II scales as ε1+5α/2 . Again, we need the
non-balancing term to be smaller than the others, so we need 1 + 5α/2 > 2α, i.e.
α > −2.
• Finally, to balance terms II and III, we need 2α − β + 1 = 2α which gives β = 1.
Then terms II and III scale as ε2α and term I scales as εα−2 , so to make term I
smaller than the others we need α − 2 > 2α, giving α < −2.
We can plot the lines in the α–β plane where these balances occur, and in the regions
between, which term (I, II or III) dominates:
We can see that there is a distinguished scaling α = −2, β = 1 where all three terms
balance. If we apply this scaling to have z = (x − a)/ε and F = ε2 f then the governing
17
ODE for F (z) (after multiplication of the whole equation by ε4 ) becomes
d2 F dF
+ F + F 2 = 0.
dz 2 dz
This is very nice: but it may not always be appropriate: the boundary conditions may
fix the size of either f or x, in which case the best you can do may be one of the simple
balance points (i.e. a point (α, β) lying on one of the lines in the diagram).
d2 f df
ε + f −f =0
dx2 dx
Find the scalings f = εα F and stretches x = a + εβ z at which two dominant terms
balance, and sketch these balance scalings in the α–β plane. Hence determine the critical
values of α and β for which all three terms balance. Give also the possible values of β if
we are constrained by the boundary conditions to have α = 0.
18
5 Matching: Boundary Layers
Consider the following equation (rather similar to exercise 7):
d2 f df
ε 2
+ +f =0
dx dx
There are two solutions. One is regular:
f0′ + f0 = 0 =⇒ f0 = a0 e−x .
At order ε we have
The second solution is singular, and the distinguished scaling (to balance the first two
terms) is δ = ε. We introduce a new variable z = (x − a)/ε to have
d2 f df
2
+ + εf = 0
dz dz
with solution
f = F0 (z) + εF1 (z) + · · ·
At order 1 we have
F0′′ + F0′ = 0 =⇒ F0′ = −B0 e−z
F0 (z) = A0 + B0 e−z .
At order ε we have
Question: Will we ever need to use both of these in the same problem?
Answer: The short answer is yes. This is a second-order differential equation, so we are
entitled to demand that the solution satisfies two boundary conditions.
Suppose, with the differential equation above, the boundary conditions are
df
f = e−1 at x = 1 and = 0 at x = 0.
dx
19
We will start by assuming that the unstretched form will do, and apply the boundary
condition at x = 1 to it:
so at x = 0,
f ′ (0) = −1 − 2ε + · · ·
This is where we have to use the other solution. If we fix a = 0 in the scaling for z, then
the strained region is near x = 0. We can re-express the boundary condition in terms of
z:
df
= 0 at z = 0.
dz
Now applying this boundary condition to our strained expansion gives:
We now have two perturbation expansions, one valid at x = 1 and for most of our region,
the other valid near x = 0. We have not determined all our parameters. How will we do
this? The answer is matching.
• One scales normally and satisfies a boundary condition somewhere away from the
tricky region: we will call this the outer solution.
• The other is expressed in terms of a scaled variable, and is valid in a narrow region,
(probably) near the other boundary. We will call this the inner solution.
20
In order to make sure that these two expressions both belong to the same real (physical)
solution to the problem, we need to match them.
In the case where the outer solution is
x = εα ξ =⇒ z = εα−1 ξ,
Here the exponential terms become very small indeed so we neglect them and have
F (z) = A0 − A0 εα ξ + εA1 + · · ·
1 = A0
5
However, occasionally you may find it quicker to pick a value of α = 1/2, say. Be warned: sometimes
there is only a specific range of α which works.
21
and at order εα ,
−ξ = −A0 ξ
which is automatically satisfied if A0 = 1. If we fix α > 1/2 then the next term is order
ε, giving
1 = A1 .
The next term in the outer expansion is order ε2α , but to match that we would have to
go to order ε2 in the inner expansion.
We have now determined all the constants to this order: so in the outer we have
Note: The beauty of the intermediate variable method for matching is that it has so
much structure. If you have made any mistakes in solving either inner or outer equation,
or if (by chance) you have put the inner region next to the wrong boundary, the structure
of the two solutions won’t match and you will know something is wrong!
22
and an inner solution
with z = (x − a)/ε.
This time we will try to fit the outer solution to the boundary condition at x = 0. We
have
df
∼ −a0 e−x + ε[a0 x − a0 − a1 ]e−x + · · ·
dx
so the condition is
df
=0 at x = 0.
dx
0 = −a0 + ε[−a1 − a0 ] + · · ·
which gives a0 = 0, a1 = 0 and so on. It is clear that we’re not going to get a solution
this way!
However, there is another problem, which appears when we try to fit the inner solution
at the other boundary. We are setting a = 1 and trying to fit F (z) = e−1 at z = 0. This
gives:
e−1 = A0 + B0 + ε[A1 + B0 + B1 ] + · · ·
so A0 = e−1 − B0 and A1 = −B0 − B1 . This seems fine, but look at the solution we get:
Remember that, now the boundary layer is at the top, the outer limit of the inner solution
will be for large negative z: in other words, all of these exponentials will be growing!
This can never match onto a well-behaved outer solution.
Key fact: The boundary layer is always positioned so that any exponentials in the inner
solution decay as you move towards the outer.
23
5.3 A worse example
This example comes from Cole’s book (available in the department library). I would not
expect you to solve equations of this difficulty without hints. The governing equation is
d2 f df
ε 2
+f −f =0
dx dx
with boundary conditions f (0) = −1, f (1) = 1.
Note that you have seen this equation before in exercise 9; you should have found the
distinguished scaling was f = ord(ε1/2 ). Here we can’t use that scaling as the boundary
conditions fix f to be ord(1).
Outer
Let us look first at the outer solution. We pose a series
Note that for both of these, d2 f0 /dx2 = 0 and so f0 is an exact solution of the equation,
and f1 = f2 = · · · = 0.
Clearly the branch f0 = 0 can’t match either of the boundary conditions, so we know our
outer solution must be
f (x) = x + C.
We have not yet found where the boundary layer will be; since the outer is so simple, we
might as well work out the constant for both possibilities now.
If the outer meets x = 1 then we have C = 0:
fouter,1 (x) = x.
fouter,0 (x) = x − 1.
Inner
What stretch do we expect for the inner? Note that the boundary conditions mean we
can’t scale f , we can only stretch x. In your solutions to exercise 9, this is equivalent
to fixing α = 0. There were two balancing scalings that crossed the axis α = 0: β = 0
(which is the outer) and β = 1, which will be our boundary layer scaling.
24
[If you didn’t do that example, simply find the scaling for x = a + εb z which balances the
first two terms of the ODE.]
We introduce z = (x − a)/ε and rewrite our differential equation:
d2 f df
2
+f − εf = 0
dz dz
Now we pose an inner expansion:
d2 F0 dF0
2
+ F0 = 0.
dz dz
We can integrate this directly once:
dF0 1 2
+ F0 = C.
dz 2
Now remember that for a boundary layer solution, we are going to need solutions which
decay exponentially to some fixed value out of the layer. This means that as z → ±∞
(but not necessarily both), we need dF0 /dz → 0 and so C ≥ 0. Let us set C = 2k 2 for
convenience.
This ODE for F0 has three different possible forms of solution. If k = 0 the solution is
2
F0 = ,
z+C
which is no good as it only decays algebraically. For k > 0 there are three solutions, two
of which work.
First we look at the possibility that |F0 | = 2k. In that case
dF0
=0 F0 = ±2k.
dz
This doesn’t have the exponential decay we need either. So we move on to the two other
cases: |F0 | < 2k and |F0 | > 2k.
25
In the two cases |F0 | < 2k and |F0 | > 2k we can solve the ODE by partial fractions.
If |F0 | < 2k then we have:
dF0
2 = 4k 2 − F02 .
dz
1 " #
2 1 1 1
1 1
dz = dF0 = + dF0
4k 2 − F02 2k 2k − F0 2k + F0
2kz + 2B = − ln (2k − F0 ) + ln (2k + F0 )
2k + F0 4k
exp [2(kz + B)] = = −1 +
2k − F0 2k − F0
4k
= 2k − F0
exp [2(kz + B)] + 1
4k 4k exp [−(kz + B)]
F0 = 2k − = 2k −
exp [2(kz + B)] + 1 exp [kz + B] + exp [−(kz + B)]
exp [kz + B] − exp [−(kz + B)] sinh [kz + B]
F0 = 2k = 2k = 2k tanh [kz + B].
exp [kz + B] + exp [−(kz + B)] cosh [kz + B]
26
Now we have two solutions which decay exponentially to some limit as z → ∞:
We can see that the tanh solution moves smoothly from one value to another over the
width of the boundary layer, whereas the coth profile cannot be given a value z = 0.
This means that the coth profile can only be used if the boundary layer is at one end or
other of the region, whereas the tanh profile can be used anywhere.
Matching
Let us try first to put the boundary layer near x = 0. The outer solution must match
the boundary condition at x = 1 so
fouter = x.
fouter = x − 1.
27
Finally, let us try having the “boundary layer” in the middle, at some general position a
between 0 and 1. This time we have two different branches of the outer solution:
F (z → −∞) = a − 1 F (z → ∞) = a.
The only profile we are allowed is the tanh profile, which goes from −2k to 2k over the
width of the layer. This fixes
28
Exercise 12: Consider the advection-diffusion equation for f (weak diffusion):
∇ · [f V ] − ε∇2 f = 0.
Vx = κx/y Vy = −κ.
(a) Substitute the velocity field into the governing equation and boundary conditions.
(b) Setting ε = 0, find a solution f0 which matches the upper boundary condition at
y = 2. [Hint: try f0 (x, y) = g(y).]
(c) Substitute your solution back into the full governing equation. What can you say
about the corrections to f0 for ε ̸= 0?
(d) The inner boundary condition at y = 1 is not satisfied. Assume that there is a
boundary layer close to y = 1. How does the size of this layer scale with ε? Assume
that derivatives with respect to x remain order 1.
(e) Introduce a scaled variable z to replace y near y = 1. Replace y in the governing
equation and give your new PDE to two orders of magnitude. Do the same for the
inner boundary condition at y = 1.
(f) Using your new PDE and boundary condition alone, calculate the first term of a
perturbation expansion for F (x, z) = f (x, y), valid within the boundary layer near
y = 1. You will not be able to determine all the constants (or even all the functions
of x) at this stage.
[Hint: if you have a PDE in which all the derivatives are with respect to z, you can
solve it like an ODE in z but all the “constants of integration” must be functions of
x.]
(g) Can your solution be made to match onto the outer solution as z → ∞ and y → 1?
(h) Continue to the next order correction in your inner expansion. What PDE must be
satisfied by the next term? What boundary conditions will be applied to it?
(i) Solve your PDE for the second term in the expansion of the inner solution. Apply
the boundary condition to get a relation between the different unknowns.
(j) Carry out matching between your outer and inner solutions. Thus find an ODE in
x for the unknown function from the calculation of (f).
(k) Solve the ODE to complete the calculation of the leading-order term in the inner
expansion. You may assume that the function is well-behaved at x = 0.
29