Example of linear differential equation with boundary layer

We will investigate the behavior of solutions to the following boundary value problem:

(ϵy)''+2y'+2y=0
y(0)=0
y(1)=1

Numerical Experiment

We begin by using Mathematica’s built in ODE solver to gain some intuition about the from the problem.

ϵ=.1;
sol1=NDSolve[{ϵ*y''[x]+2*y'[x]+2*y[x]==0,y[0]==0,y[1]==1},y[x],x];
plotN1=Plot[Evaluate[y[x]/.sol1],{x,0,1},AxesLabel->{"x","y"},TicksStyle->Large,PlotStyle->Thick,PlotRange->{0,3}];
ϵ=.05;
sol2=NDSolve[{ϵ*y''[x]+2*y'[x]+2*y[x]==0,y[0]==0,y[1]==1},y[x],x];
plotN2=Plot[Evaluate[y[x]/.sol2],{x,0,1},AxesLabel->{"x","y"},TicksStyle->Large,PlotStyle->{Thick,Red},PlotRange->{0,3}];
ϵ=.01;
sol3=NDSolve[{ϵ*y''[x]+2*y'[x]+2*y[x]==0,y[0]==0,y[1]==1},y[x],x];
plotN3=Plot[Evaluate[y[x]/.sol3],{x,0,1},AxesLabel->{"x","y"},TicksStyle->Large,PlotStyle->{Thick,Green}];
Show[plotN1,plotN2,plotN3]

From the numerics we see that the near x=0 there is a boundary layer where the second derivative is very large. The rest of the domain consists of an outer layer where the function behaves like an exponential function.

O(1) uniform asymptotic approximation.

We now construct a O(1) unfirom asymptotic approximation by pieciing together outer and inner approximations.

O(1) Outer Approximation

We assume an asymptotic expansion of the following form for the outer layer approximation: (Subscript[y, out])=(Subscript[y, 0])+(Subscript[ϵy, 1])+(ϵ^2)(Subscript[y, 2])+… . The equation for (Subscript[y, 0])is given by:

(Subscript[y, 0])'+(Subscript[y, 0])=0
y(1)=0

DSolve[{y0'[x]+y0[x]==0,y0[1]==1},y0[x],x]

y0[x_]:=Exp[1-x];
y0[x]

The lowest order outer solution is highlighted above in light green. Lets compare how well the outer expansion matches with the numerical solution.

plotOut1=Plot[y0[x],{x,0,1},AxesLabel->{"x","y"},TicksStyle->Large,PlotStyle->Thick,PlotRange->{0,3}];
Show[plotN1,plotN2,plotN3,plotOut1]

As expected the outer approximation does an excellent job throughout most of the domain except near x=0.

O(1) Inner Approximation

To determine the precise width of the boundary layer we rescale the independant varaible by X=(ϵ^(-α))x. This gives us the following differential equation

(ϵ^(1-2α))(((d^2)y)/(dX^2))+2(ϵ^(-α))(dy/dX)+2y=0

To balance terms we must have that the first two terms balance. This gives us that α=1. We now expand as (Subscript[y, in])=(Subscript[Y, 0])+(Subscript[ϵY, 1])+(ϵ^2)(Subscript[Y, 2])+…. At the lowest order this gives us the following equation.

(((d^2)(Subscript[Y, 0]))/(dX^2))+2((Subscript[dY, 0])/dX)=0
Y(0)=0.

DSolve[{Y0''[X]+2*Y0'[X]==0,Y0[0]==0},Y0[X],X]

Y0[X_]:=(1/2)*A*(1-Exp[-2*X]);

Note that we still have undetermined contstant of integration, namely A. This constant must be determined from matching the inner and outer layers.

Matching

To match the inner and outer boundary layers we select an intermediate scale: (Subscript[x, η])=(ϵ^(-β))x where 0<β<1. We now expand both the inner and outer solutions in terms of (Subscript[x, η]) and select A so that they match to lowest order.

Clear[ϵ]
y0[(ϵ^β)*xη]

Series[Exp[1-y],{y,0,3}]/.y->(ϵ^β)*xη

y0∼e-e xη (ϵ^β)+(1/2) e (((xη (ϵ^β)))^2)-(1/6) e (((xη (ϵ^β)))^3)+((O[xη (ϵ^β)])^4)

The above expression highlighted in light pink is the asymptotic expansion for the outer solution y0 as ϵ->0 in terms of the rescaled coordinate (Subscript[x, η]). Now lets look at the inner approximation.

Y0[x/ϵ]/.x->xη*(ϵ^β)

Since β-1<0 the term (e^(-2 xη (ϵ^(-1+β)))) is transcendentally small. Therefore, in this scaled region.

Y0∼(1/2)A

Consequently to ensure that the inner and outer approxmations match to lowest order we must have that A=2e. This gives us the following

A=2*E;
Y0[x/ϵ]

The lowest order inner solution is highlighted above in green. Lets see how the inner approximation matches with the exact solution for ϵ=.05. Do not forget the Y0 is defined in terms of the rescaled coordinate X=x/ϵ.

ϵ=.05;
PlotInner1=Plot[Y0[x/ϵ],{x,0,1},AxesLabel->{"x","y"},TicksStyle->Large,PlotStyle->Thick,PlotRange->{0,3}];
Show[plotN2,PlotInner1]

The inner layer does an excellent job matching the boundary layer but of course does recover the outer solution. Now, lets see how to combine them into a composite expansion.

Composite Expansion

The idea of a composite expansion is that we should be able to simply add the two expansion together and subtract the part where they overlap. The overlap comes from exactly the terms that we matched in above. Recall that the term that matched in the inner and outer layer was simply e.

y0Com[x_]:=y0[x]+Y0[x/ϵ]-E;

Now lets compare with the numerical solution.

ϵ=.05;
PlotCom1=Plot[y0Com[x],{x,0,1},AxesLabel->{"x","y"},TicksStyle->Large,PlotStyle->{Thick,Dashed},PlotRange->{0,3}];
Show[plotN2,PlotCom1]

This is so amazing that this works so well!! Lets look how does for more modest values of ϵ.

ϵ=.1;
PlotCom1=Plot[y0Com[x],{x,0,1},AxesLabel->{"x","y"},TicksStyle->Large,PlotStyle->{Thick,Dashed},PlotRange->{0,3}];
Show[plotN1,PlotCom1]

Lets take stock on what we have done. We have broken up a differential equation into two parts that could be solved relatively painlessly and then pieced them together to give an excellent approximation to the solution to differential equation. If we can extend this to nonlinear problems we will have a very powerful tool! Lets look at solving the equation to higher orders. Below is highlighted in green the composite expansion.

Clear[ϵ]
Expand[y0Com[x]]

y~(e^(1-x))-(e^(1-((2 x)/ϵ)))

O(ϵ) uniform asymptotic approximation.

We now construct a O(ϵ) unfirom asymptotic approximation by pieciing together outer and inner approximations.

O(1) Outer Approximation

Recall that we are assuming yin~ (Subscript[y, 0])+(Subscript[ϵy, 1])+…. Lets first determine the equation for (Subscript[y, 1]). We will simply substitute the full expansion into the differential equations and look at the coefficient of the (ϵ^1)term.

yin[x_]:=y0[x]+ϵ*y1[x];
CoefficientList[ϵ*yin''[x]+2*yin'[x]+2*yin[x],ϵ]

DSolve[{(E^(1-x))+2 y1[x]+2 y1'[x]==0,y1[1]==0},y1[x],x]

Clear[ϵ];
y1[x_]:=(1/2)Exp[1-x](1-x);
y0[x]+ϵ*y1[x]

Highlighted in green above is the first two terms in the outer layer expansion. Lets compare how well the outer expansion matches with the numerical solution

ϵ=.1;
plotOut2=Plot[{y0[x]+ϵ*y1[x],y0[x]},{x,0,1},AxesLabel->{"x","y"},TicksStyle->Large,PlotStyle->{{Thick,Red},{Thick,Green}},PlotRange->{0,3}];
Show[plotN1,plotOut2]

The red curve is our improved outer approximation while the green curve is the original approximation.

O(1) Inner Approximation

Lets now do a similar calculation for the next term in the inner solution.

Clear[Y1]
Clear[ϵ];
yin[X_]:=Y0[X]+ϵ*Y1[X];
CoefficientList[yin''[X]+2*yin'[X]+2*ϵ*yin[X],ϵ]

DSolve[{2 E (1-(E^(-2 X)))+2 Y1'[X]+Y1''[X]==0,Y1[0]==0},Y1[X],X]

The above solution is a bit of a mess that we should be able to write in a more simple form.

Expand[-(1/2) (E^(-2 X)) (E-(E^(1+2 X))+2 E X+2 (E^(1+2 X)) X+C[1]-(E^(2 X)) C[1])]

Simplify[(E/2)-(1/2) (E^(1-2 X))]
Simplify[-E X-(E^(1-2 X)) X]
Simplify[((C[1])/2)-(1/2) (E^(-2 X)) C[1]]

Notice we have the repeated term (1-(e^(-2 X))) multiplied times a constant. Therefore, we can express the solution in the following form:

Y1[X_]:=B*(1-(E^(-2 X)))-E (1+(E^(-2 X))) X;

Again, we have an undetermined coefficient B which must be determined from matching the inner and outer layers.

Matching

To match the inner and outer boundary layers we select an intermediate scale: (Subscript[x, η])=(ϵ^(-β))x where 0<β<1. We now expand both the inner and outer solutions in terms of (Subscript[x, η]) and select A so that they match to lowest order.

Clear[ϵ]
y0[(ϵ^β)*xη]+ϵ*y1[(ϵ^β)*xη]

Series[y0[y]+ϵ*y1[y],{y,0,3}]/.y->(ϵ^β)*xη

y0+ϵy1 ~(e+((e ϵ)/2))+(-e-e ϵ) xη (ϵ^β)+((e/2)+((3 e ϵ)/4)) (((xη (ϵ^β)))^2)+(-(e/6)-((e ϵ)/3)) (((xη (ϵ^β)))^3)+((O[xη (ϵ^β)])^4)

The above expression highlighted in light pink is the asymptotic expansion for the outer solution y0 as ϵ->0 in terms of the rescaled coordinate (Subscript[x, η]). Now lets look at the inner approximation.

Expand[Y0[x/ϵ]+ϵ*Y1[x/ϵ]/.x->xη*(ϵ^β)]

Since β-1<0 terms (e^(-2 xη (ϵ^(-1+β)))) are transcendentally small. Therefore, we have that

Y0[x/ϵ]+ϵ*Y1[x/ϵ] ~ e+B ϵ-e xη (ϵ^β)

Consequently to ensure that the inner and outer approxmations match to lowest order we must have that B=(e/2). This gives us the following

B=E/2;
Y0[x/ϵ]+ϵ*Y1[x/ϵ]

e (1-(e^(-((2 x)/ϵ))))-e (1+(e^(-((2 x)/ϵ)))) x+ (ϵ/2) e (1-(e^(-((2 x)/ϵ))))

Above is the two terms in the inner approximation. Lets see how the inner approximation matches with the exact solution for ϵ=.05. Do not forget the Y0 and Y1 are defined in terms of the rescaled coordinate X=x/ϵ.

ϵ=.05;
PlotInner2=Plot[{Y0[x/ϵ],Y0[x/ϵ]+ϵ*Y1[x/ϵ]},{x,0,1},AxesLabel->{"x","y"},TicksStyle->Large,PlotStyle->{{Thick,Blue},{Thick,Green}},PlotRange->{0,3}];
Show[plotN2,PlotInner2]

Above the blue curve is the O(1) inner approximation while the green is the O(ϵ) inner approximation.

Composite Expansion

To creat the composite expansion we again should be able to simply add the two expansion together and subtract the part where they overlap. The overlap comes from the terms that matched in the two expansions above.

y1Com[x_]:=y0[x]+ϵ*y1[x]+Y0[x/ϵ]+ϵ*Y1[x/ϵ]-E-((E ϵ)/2)+E*x;

Now lets compare with the numerical solution.

ϵ=.05;
PlotCom2=Plot[{y0Com[x],y1Com[x]},{x,0,1},AxesLabel->{"x","y"},TicksStyle->Large,PlotStyle->{{Thick,Dashed},{Thick,Dashed}},PlotRange->{0,3}];
Show[plotN2,PlotCom2]

Notice that the second dashed curve essentially matches the solution perfectly! This is really amazing! Again lets check what happens when ϵ=.1.

ϵ=.1;
PlotCom2=Plot[{y0Com[x],y1Com[x]},{x,0,1},AxesLabel->{"x","y"},TicksStyle->Large,PlotStyle->{{Thick,Dashed},{Thick,Dashed}},PlotRange->{0,3}];
Show[plotN1,PlotCom2]

Lets now write out the full two term composite expansion.

Clear[ϵ]
CoefficientList[Expand[y1Com[x]],ϵ]

Simplify[(E^(1-x))-(E^(1-((2 x)/ϵ)))-(E^(1-((2 x)/ϵ))) x]

y~(e^(1-x))-(e^(1-((2 x)/ϵ)))(1+x)+(ϵ/2)((e^(1-x))(1-x)-(e^(1-((2 x)/ϵ))))