First example on boundary layers

Analytic solution compared to boundary layer approximation for
ϵ y'' [x] + (1+ϵ) y' [x] + y [x] = 0,  y [0] = 0, y [1] = 1.

Analytic solution

The matched asymptotic expansion approach does not require equations to be linear. However, since this equation is not only linear, but also has constant coefficients, it can be solved in simple closed form. We store the solution in the variable yexact for display and for later comparison against our asymptotic approximation

ex = 
        DSolve[{ϵ y''[x] + (1 + ϵ) y'[
                    x] +  y[x] == 0, y[0] == 0, y[1] == 1}, y[x], 
          x];
yexact = y[x]  /. ex  [[1]]

Clear[y]

The solutions for ϵ = 10 ^ (-k) ,  {k, 1, 3} look as follows:

\(\(Do\ [\ \ \ pl[k] = 
          Plot[yexact /. ϵ -> 10^\((\(-k\))\), {x, 0, 1}, 
            DisplayFunction -> Identity], 
	\ \ \ \ \ \ \ \ {k, 1, 
          3}];\)\)
Show[pl[1], pl[2], pl[3], 
      DisplayFunction -> $DisplayFunction]

This picture illustrates clearly that a boundary layer is developing at the left boundary. We could also have concluded that directly from the ODE very heuristically as follows:

We first extend the ODE to the time dependent PDE    ∂(y/∂) t =  ϵ (∂^2) (y/∂)(x^2) + (1+ϵ) ∂(y/∂)x+ y   (the equilibrium solution as t->∞ becomes the ODE).  We can now interpret the effect of the different terms of the PDE as follows:

1.    ∂(y/∂) t =  ϵ  (∂^2) (y/∂)(x^2)       This is a well-posed heat equation (thanks to ϵ>0 and our choice of sign when we added  ∂(y/∂) t ),
2.    ∂(y/∂) t =   (1+ϵ) ∂(y/∂)x   This transport information to the left,  and
3.    ∂(y/∂) t =   y.              This causes y to grow away from zero to the extent permitted by the other terms.

The net result of these diverse effects is that the solution is rising between the boundaryes while being pushed to the left boundary. The final time equilibrium (the solution to the ODE) relies on the weak dissipation to stop the solution from getting discontinuous at the left edge.

Outer expansion

We will expand both the inner- and outer solutions in powers of ϵ  up to a certain level, say up to power n = 4;

n = 4;

Since, in the region away from the boundary layer, the limit of ϵ -> 0 appears to agree with the solution for ϵ = 0 (the ϵy''-term is then of significance only in the boundary layer), our first goal is to find a regular outer expansion which satisfies the right boundary condition y [1] = 1. For the ODE

odeouter =  ϵ y''[x] + (1 + ϵ) y'[
            x] + y[x]

\(TraditionalForm
    \`we\ assume\ an\ expansion\ in\ \((yet\ unknown)\) coefficient\ functions
      \ \ y\_k[x]\ \ as\ follows\)

y[x_] := Sum[(Subscript[y, k])[x] ϵ^k,  {k, 0, n}]

With this expansion substituted into the ODE, we get an expansion in ϵ in which every single coefficient function should be equal to zero.

ser = Series[ odeouter,  {ϵ, 0, n}]

We next generate separate equations for the different coefficients

odes = LogicalExpand[ser == 0]

We now need to loop over these ODEs, and solve them in turn with the appropriate boundary conditions.We first move the ODEs over to a usual list format  (else, some things go wrong...)

ode = Table[odes[[i + 1]], {i, 0, n}]

\(TraditionalForm\`Now, \ 
    loop\ over\ the\ different\ ODEs.\ We\ need\ to\ use\ an\ intermediate\ 
      variable\ \ temp\ to\ avoid\ infinite\ recursions\ 
      \((that\ a\ statement	\ y\_k[x_]\  = \ 
          y\_k[x] /. sol[\([1]\)]\ \ otherwise\ would\ have\ produced; \ 
        \(Note\ also\ that\ we\ need\ to\ use\  = \)\ \ , \ 
        not\  := \ \ 
          in\ these\ assignments.
We\ have\ enough\ BC\ for\ each\ of\ the\ 
              ODEs\ \((which\ are\ first\ order\ only)\)\  - \ 
            no\ undetermined\ constants\ will\ be\ \(produced.\)					
        \)\)

\(Do\ [\ 
      sol = DSolve[{ode[\([k + 1]\)]\ , y\_k[1] == \ If[k == 0, 1, 0]}, 
          y\_k[x], x]; 
	temp\  = \ y\_k[x] /. sol[\([1]\)]; \ 
      y\_k[x_] = temp, \ {k, 0, n}]\)

The series expansion for   y [x]  in powers of ϵ  becomes

outerexp = y[x]

In general, we would here have obtained a power series in ϵ with different functions of x as coefficients.In this special case, all coefficients but the first one happened to vanish (as we can see from the analytic solution).

Inner expansion

The equation      odeouter==0    fails to describe the boundary layer. We therefore want next to carry out a stretching  (change of variables)  x = δ X  , y = Y   so that the ODE becomes regular in a neighborhood of   X = 0. Substituting the changes of variable into odeouter amounts only to replacing y'[x] by Y'[X]/δ and y''[x] by Y''[X]/(δ^2) - we do this 'manually' since y[x] is no longer available to us as an unspecified function:

odeinner = Y[X]+(((1 + ϵ) Y'[X])/δ)+((ϵ Y''[X])/((δ^2)))

\(TraditionalForm
    \`Since\ \ δ\  -> \ \(0\ when\ ϵ\  -> \ 0\), \ 
    and\ because\ we\ want\ Y[X]\ and\ its\ derivatives\ to\ remain\ O[1], 
    \ \ the\ two\ dominant\ terms\ which\ need\ to\ balance\ are\ \ 
        \(Y'\)[X]/δ\ \ and\ ϵ/δ\^2\ 
        \(\(Y'\)'\)[x].\ This\ leads\ us\ to\ choose\ δ\  = \ 
      \(ϵ : \)\)

δ = ϵ;

The ODE   odeinner  becomes now  (after multiplying by ϵ to get the coefficient for the highest derivative to be one):

Expand[odeinner  *=  ϵ]

By the same series of steps as for   odeouter, we obtain a regular perturbation expansion for  odeinner

Y[X_] := Sum[(Subscript[Y, k])[X] ϵ^k,  {k, 0, n}];
ser = Series[ odeinner,  {ϵ, 0, n}];
odes = LogicalExpand[ser == 0];
ode = Table[odes[[i]], {i, 1, n + 1}]

For each of these ODEs, we this time get one integration constant, which Mathematica writes as C[1] or C[2] (each equation is second order, but has only one BC specified). This undetermined constant needs to be immediately renamed so as not to interfere with constants of the same name that get generated for the subsequent ODEs. Else all is similar to before

\(Do\ [\ sol = DSolve[{ode[\([k + 1]\)]\ , Y\_k[0] == 0}, Y\_k[X], X]; 

      temp = \ \((Y\_k[X] /. sol[\([1]\)])\) /. C[z_] -> a\_k; \ 
      Y\_k[X_] = temp, \ {k, 0, n}]\)

The inner expansion (in ϵ)  becomes

Y[X]

Since we will 'match' in two ways, we store for now away a second copy of this expansion

Y2[X] = Y[X] /.  (Subscript[a, k])_ -> (Subscript[b, k])

Matching

We have now generated both the outer

y[x]

and inner expansions:

Y[X]

In general, both of these would be power series expansions in ϵ. The outer expansion will normally not contain any free parameters. Our task is to determine the constants in the inner expansion by requiring the two to match in some overlap region. There are two main ways to achieve this is

i.  'van Dykes matching rule' - this is usually the simpler method , but it may fail in tricky cases, or 		
ii.  Introduce a scaling ξ in-between  x and X, and require the two expansions to be identical when expanded in this new variable. This procedure will always work.

We will use both methods below, starting with  'van Dykes matching rule':

van Dyke's matching rule

Change outer expansion to inner (stretched) variable:

y[x]  /.  x -> ϵ  X

Expand for ϵ->0.

ser1 = Series[%, {ϵ, 0, n}]

Recall the series expansion for the inner (stretched) variable:

ser2 = Collect[Expand[Y[X]], ϵ]

Remove exponentials (which will be asymptotically small in the limit of X large)

ser2 = ser2 /. (E^(-X)) -> 0

Reconcile the two expansions

LogicalExpand[ser1 == ser2]

rules = Solve[%]

In this case, it happened (uncharacteristically) that we got perfect matching after just one term - all further terms will vanish. We note that the matching coefficients feature no X-dependence although lots of X-es appeared in the system we solved.  This forms an excellent consistency check on our whole procedure.

We apply now the substitution rules above to get (Subscript[a, k]) to actually assume the values above.

Do[(Subscript[a, k]) = (Subscript[a, k]) /. rules[[1]], {k, 0, n}]

The inner expansion is therefore

Y[X]

The inner expansion is therefore

Matching through an intermediate variable

\(TraditionalForm
    \`By\ choosing\ \ ξ = 
      x/ϵ\^α\ \ with\ 0 < α < 1, \ 
    ve\ can\ create\ an\ intermediate\ variable\ in - 
        between\ x\ and\ X.\ We\ can\ express\ both\ the\ outer\ and\ inner\ 
          expansions\ in\ this\ new\ variable\ by\ substituting\ \ x -> 
      \(\(ϵ\^α\) ξ\ \ and\ \ \ X -> 
        \(ϵ\^\(α - 1\)\) ξ\ \ \(respectively : \)\)\)

OUTER:

Series[y[x], {x, 0, n}] /. 
      x -> (ϵ^α) ξ

INNER:

(Expand[Y2[X]] /. (E^(-X)) -> 0)  /.  
      X -> (ϵ^(α - 1)) ξ

COMPARE:

For every power of ϵ,  
    these two expansions must agree,  
    term - by - 
      term. There are many more terms to match than there are 
        coefficients to choose. The ability to cancel all terms 
        therefore forms an effective consistency check. Here,  
    we note first that (ϵ^0) 
        (i.e. the ϵ - independent coefficient) forces 
        (Subscript[b, 0]) = E. All the powers (ϵ^k)α,  k = 1, 2, 
    ... are then seen to cancel out. There are now no coefficients
       left in the outer expansion,  
    and all remaining coefficients in the inner expansions vanish as
         well if we choose (Subscript[b, k]) = 0 ,   k = 1, 2, 
    ... . We have thus obtained the same matching result as 
      with van Dyke' s rule.

Creation of a uniformly valid expansion

The task is to combine inner and outer expansions into a unifirmly valid expansion

By inspection, the inner expansion is the beginning of the Taylor expansion in ϵ for

(E^(1 - X ϵ)) - (E^(1 - X))

In the outer variable, this becomes

innerexp = %  /.  X -> (ϵ^(-1)) x

The expressions  outerexp and innerexp have in common

match = (E^(1 - x))

The uniform approximation is obtained as

uniform = outerexp + innerexp - match

It is worth noting that, although this expression represents matching to infinite order, and is valid across the full interval [0,1], it still does not agree exaxtly to the analytic solution - the two differ by an exponentially small quantity as ϵ->0. This discrepancy is asymptotically smaller than any power of ϵ, and is undetectable with perturbation analysis.

Comparison between exact and matched asymptotic expansion

We plot the uniform approximations (dashed) against the exact ones (solid) for the same three cases  ϵ = 10 ^ (-k) ,  {k, 1, 3} :

\(\(Do\ [
        pl[k + 3] = 
          Plot[uniform /. ϵ -> 10^\((\(-k\))\), {x, 0, 1}, 
            PlotStyle -> Dashing[{0.05, 0.05}], 
            DisplayFunction -> Identity], 
	\ \ \ \ \ \ \ \ {k, 1, 
          3}];\)\)
Show[Table[pl[k], {k, 1, 6}], 
      DisplayFunction -> $DisplayFunction]

Since the difference is not graphically visible in any of these cases, we next plot the difference in the least accurate case (ϵ = 0.1).

Plot[(yexact - uniform) /. ϵ -> 1/10 , {x, 0, 
        1}]

Our asymptotic result satisfies exactly the left BC, but not the right one.