Second example on boundary layers

Numerical solution compared to boundary layer approximation for
ϵ y'' [x] + (1+x) y' [x] + y [x] = 0,  y [0] = 1, y [1] = 1.  (cf. Example 1 in B&O page 431)

Analytic solution

Remove["Global`*"]

Mathematica manages to find a solution in closed form, but it is not simple; with a small change in the equation, it would not be likely to find one:

ode = ϵ y''[x] + (1 + x) y'[x] +  y[x] == 
        0

DSolve[{ode, y[0] == 1, y[1] == 1}, y[x], x]

Clear[y]

Numerical solution

nsol = 
      NDSolve[{ode  /.  ϵ -> 0.1, y[0] == 1, y[1] == 1}, 
        y[x], {x, 0, 1}]

The output of NDSolve takes the form of an interpolating function - essentially a table of intermediate values, from which the function values at arbitrary points can be obtained by standard interpolation  (for this, Mathematica internally uses divided differences to construct Lagrange- or Hermite interpolating polynomials) .

\(\(\(
\)\(Plot[Evaluate[y[x]\  /. \ nsol], {x, 0, 1}]\)\)\)

The numerical algorithm that Mathematica uses for boundary value problems is a version of 'shooting' (technically known as the Gelfand-Lokutsiyevskii chasing method) - this is ill suited for singular perturbation problems. In this example, it will fail for ex. if ϵ is 0.05 (or less). Still, even with a rather large value of ϵ, we are getting an idea of what the solution looks like.

Clear[y]

Heuristic argument of what we would expect the solution to look like

Just like in the previous example, we can conclude that ODE has a boundary layer at left:

We first extend the ODE to the time dependent PDE    ∂(y/∂) t =  ϵ (∂^2) (y/∂)(x^2) + (1+x) ∂(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+x) ∂(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 boundaries 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;

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

\(TraditionalForm
    \`Assume\ an\ expansion\ and\ find\ successive\ perturbation\ equations\ 
      as\ before\)

y[x_] := Sum[(Subscript[y, k])[x] ϵ^k,  {k, 0, n}];
ser = Series[ odeouter,  {ϵ, 0, n}];
odes = LogicalExpand[ser == 0];
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]

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 + δ X) 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] == \ If[k == 0, 1, 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]

Matching  (using intermediate variable)

We have now generated both the outer

y[x]

and inner expansions:

Y[X]

Matching through an intermediate variable

\(TraditionalForm
    \`By\ choosing\ \ ξ = 
      x/ϵ\^α\ \ with\ 0 < α < 1, \ 
    we\ 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:

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

INNER:  (Matching will involve increasingly large X, so remove (e^(-X))-terms)

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

DIFFERENCE:

diff = outα - innα

diff /. {ϵ -> 0, ξ -> 0}

(Subscript[a, 0]) = (Subscript[a, 0]) /. Solve[% == 0][[1]]

diff

\(TraditionalForm
    \`Notice\ that\ when\ \ we\ eliminated\ the\ constant\ term, \ 
    \(we\ also\ got\ rid\ of\ all\ terms\ of\ form\ \ \ \ 
        \(ϵ\^α\) ξ\  - \), 
    \(\(ϵ\^\(2  α\)\) ξ\^2 - \), \ 
    \(ϵ\^\(3  α\)\) ξ\^3 - \ etc.\ Since\ \ α > 0, 
    \ next\ term\ along\ in\ size\ is\ therefore\ ϵ\^1.\ So\ we\ \ 
      divide\ by\ ϵ\ and\ isolate\ constant\ factor\ \(again : \)\)

Expand[diff/ϵ]  /.  {ϵ -> 0, ξ -> 0}

(Subscript[a, 1]) = (Subscript[a, 1]) /. Solve[% == 0][[1]]

diff

Not only did we get rid of our  ϵ-terms, we this time also lost terms    
      (ϵ^(1 + α)) ξ  - , 
      (ϵ^(1 + 2  α)) (ξ^2) - ,  
      (ϵ^(1 + 3  α)) (ξ^3) -  etc.  The term (59/2) (ϵ^(1 + 4 α)) is left only because we ran out of terms of that kind in one of the expansions. So it should be gone too - easily verified by increasing n. So

Expand[
        diff/(ϵ^2)]  /.  {ϵ -> 0, ξ -> 0}

(Subscript[a, 2]) = (Subscript[a, 2]) /. Solve[% == 0][[1]]

diff

\(TraditionalForm
    \`etc.\ \ This\ matching\ can\ be\ continued\ to\ give\)

\(\(\(
\)\(a\_3 = 447/16\ ; \ a\_4 = 12855/64\ ;\)\)\)

After that, we would have needed more terms initially. Our inner expansion is now

Y[X]

and thus, both inner and outer expansions are now determined to O[ϵ]^5.

Determining a uniformly valid approximation

The expansion for Y[X] contains both exponentially decaying terms and polynomial terms. The former describe the structure of the boundary layer whereas the latter are beginnings of the Taylor expansions of the successive terms of the outer solution. To get the uniform expansion, we first remove the matching part

YXremovematch =  
      Collect[Simplify[
          Y[X] - (Expand[ Y[X]]  /.  (E^(-X)) -> 
                  0)], ϵ]

and we then add this to the outer expansion

uniform = Collect[YXremovematch + y[x], ϵ]

\(TraditionalForm\`\(B &\) 
      O\ gives\ this\ same\ expansion\ at\ the\ bottom\ of\ page\ 433\ up\ to\ order\ ϵ\^3. \ In\ that\ term, \ 
    they\ have\ a\ coefficient\ \  - \ 
      1949\/72\ where\ we\ have\  - \(1341\/48\)\ \  . \ \ Strangely, \ 
    Mathematica\ \ seems\ to\ come\ up\ with\ a\ different\ answer\ for\ this\ coefficient\ \((but\ not\ for\ any\ others)\)\ dependent\ of\ which\ \ version\ of\ Mathematica\ I\ use . \ 
        So\ this\ smells\ like\ a\ Mathematica\ \(\(bug\)\(.\)\(\ \ \)\)\)

Comparison to numerical solution

The uniform expansion above has both x and X as independent variables. In order to plot it - and to compare with the numerical solution to the problem - we convert first X to x:

uniformx  =  uniform  /.  X -> x/ϵ ;

Plot[uniformx  /.  ϵ -> 0.1, {x, 0, 1}]

This looks almost identically the same as the numerical solution, so we plot the difference

Plot[(Evaluate[
            y[x]  /.  nsol]) - (uniformx  /.  ϵ -> 
              0.1), {x, 0, 1}]