Solution of a simple ODE by a singular perturbation method

This notebook has been written in Mathematica by 

Mark J. McCready
Professor of Chemical Engineering
University of Notre Dame
Notre Dame IN 46556
USA

mjm@nd.edu

http://www.nd.edu/~mjm/

It is copyrighted to the extent allowed by whatever laws pertain to the World Wide Web and the Internet.

I would hope that as a professional courtesy that this notice remain visible to other users. 
There is no charge for copying and dissemination   

Version:  4/20/16

Reference:  E. J. Hinch, Pertubation Methods, Cambridge, 1991.

Example Problem from Chapter 5.

As a very simple model of the Boundary-Layer equation flow past a flat plate (the plate is at x =0) we could write.

Here is the exact answer

DSolve[{ ϵ y''[x] +    y'[x]==-Exp[-x], y[0]==0,y[1] ==1},y[x],x]

why=y[x]/.%[[1]]

Plot[Re[why/.{ϵ->.01}],{x,0,1},PlotRange->All]

Plot[why/.ϵ->.1,{x,0,1},PlotRange->All]

Plot[{why/.ϵ->.02,why/.ϵ->.01,why/.ϵ->.05,why/.ϵ->.1},{x,0,1},PlotRange->All]

So there is a boundary-layer at x -->1.  It gets thinner as ϵ increases

Plot[{why/.ϵ->.02,why/.ϵ->.01,why/.ϵ->.05,why/.ϵ->.1},{x,0,.1},PlotRange->{0,.5}]

Numerical solutions.

We start off well

AbsoluteTiming[ans =NDSolve[{ .01 y''[x] +  y'[x]  ==- Exp[-x], y[0]==0,y[1] ==1},y[x],x]]

Plot[y[x]/.ans,{x,0,1}]

Plot[y[x]/.ans,{x,.9,1}]

AbsoluteTiming[ans =NDSolve[{ .001 y''[x] +  y'[x]  ==- Exp[-x], y[0]==0,y[1] ==1},y[x],x]]

Plot[y[x]/.ans,{x,.99,1},PlotRange->All]

But it gets harder “wander” begins to creep in...

AbsoluteTiming[ans =NDSolve[{ .0001 y''[x] +  y'[x]  ==-Exp[-x], y[0]==0,y[1] ==1},y[x],x]]

Plot[y[x]/.ans,{x,0,.001}]

Plot[y[x]/.ans,{x,0,.1}]

Even if it looks OK on the whole range

Plot[y[x]/.ans,{x,0,1},PlotRange->All]

For higher Re = 1/ϵ

AbsoluteTiming[ans =NDSolve[{ .000001 y''[x] +  y'[x]  ==-Exp[-x], y[0]==0,y[1] ==1},y[x],x]]

Plot[y[x]/.ans,{x,0,1}]

Plot[y[x]/.ans,{x,0,.001}]

Plot[y[x]/.ans,{x,0,1},PlotRange->All]

I think that Mathematica is cheating!

AbsoluteTiming[ans =NDSolve[{ .00000001 y''[x] +  y'[x]  ==-Exp[-x], y[0]==0,y[1] ==1},y[x],x]]

Plot[y[x]/.ans,{x,0,1}]

Let’s see what we can do with an expansion

Outer solution

We expect a boundary layer at x =0.  Hence an expansion over the entire domain will likely not work.  We will have an “outer” expansion in unscaled variable, x, and then do a second “inner” expansion in a rescaled variable, ξ.

Here is the equation.  I will use a “general” function, h[x] and express it as a derivative.  This actually makes the analysis go a bit easier.  In the end, h[x] = Exp[-x]

eq =  ϵ y''[x] +  y'[x]  -D[ h[x],x]

The “outer” expansion is

perturbexpansion = Y0[x] + ϵ Y1[x] + ϵ^2 Y2[x]+ ϵ^3 Y3[x]

result1 =eq/.{y[x] ->perturbexpansion,Derivative[a1_][y][x]:>  (∂_({x,a1}))perturbexpansion}

result2 =Expand[result1]

result3 = Collect[result2 ,ϵ]

So at lowest order we would have

eqϵm1 = Coefficient[result3,ϵ,0]

DSolve[{eqϵm1 ==0,Y0[1] ==1},Y0[x],x]

why0 = Y0[x]/.%[[1]]

The next order...

eqϵ0 = Coefficient[result3,ϵ,1]

DSolve[{(eqϵ0 ==0)/.Y0''[x]->D[1-h[1]+h[x],{x,2}],Y1[1] ==0},Y1[x],x]

why1 =Y1[x]/.%[[1]]

Another term is:

eqϵ1 = Coefficient[result3,ϵ,2]

DSolve[{(eqϵ1 ==0)/.Y1''[x]->D[h'[1]-h'[x],{x,2}],Y2[1] ==0},Y2[x],x]

why2 = Y2[x]/.%[[1]]

We can see the pattern but one more...

eqϵ2 = Coefficient[result3,ϵ,3]

DSolve[{(eqϵ2 ==0)/.Y2''[x]->D[-h''[1]+h''[x],{x,2}],Y3[1] ==0},Y3[x],x]

why3 = Y3[x]/.%[[1]]

So we see how this goes and could get as many terms as we want.  But this does not help us close to x=0

Inner solution

We need to rescale the “x” close to x =0

We try a powerlaw scaling in ϵ with an as yet undertimed value

yee = (x)/ϵ^α

Solve[yee ==ξ,x]

rescaledeq = eq/.{x-> (ϵ^α) ξ , y[x]->ψ[ξ],y'[x]-> 1/(ϵ^α) ψ'[ξ], y''[x]->1/(ϵ^(2 α))ψ''[ξ]}

Since we don’t know the function h (and in this case we want it to be an exponential which could change value very fast, we need to expand this around ξ = 0

hpexpansion = Normal[Series[h'[(ϵ^α) ξ],{ξ,0,3}]]

rescaledeq1 = Expand[rescaledeq/.h'[(ϵ^α) ξ]->hpexpansion]

The inner solution will be:

innerexpansion = ψ0[ξ] + (ϵ^α) ψ1[ξ] + (ϵ^(2 α)) ψ2[ξ]+ (ϵ^(3 α)) ψ3[ξ]

inner1 =rescaledeq1/.{ψ[ξ] ->innerexpansion,Derivative[a1_][ψ][ξ]:>  (∂_({ξ,a1}))innerexpansion}

inner2 = Expand[inner1]

inner3 = Collect[inner2,ϵ]

We see that we need ϵ ((ϵ^(-2 α)) (ψ0)''[ξ]) and (ϵ^(-α)) (ψ0)'[ξ] to be of the same order...

So α is chosen to make this happen

Solve[-α == 1 - 2 α,α]

We continue with α = 1

inner4 = inner3/.α->1

inner5 =Collect[inner4,ϵ]

The inner ode at 0 order is:

ieqϵm1 = Coefficient[inner5,ϵ,-1]

We solve this to match ψ = 0 at ξ =0 (i.e., x->0)

innersolution1 = DSolve[{ieqϵm1==0,ψ0[0] ==0},ψ0[ξ],ξ]

psi0 = Expand[ψ0[ξ]/.innersolution1[[1]]]

Let’s get another term of the inner solution.

ieqϵ0 = Coefficient[inner5,ϵ,0]

ieqϵ01 = ieqϵ0

innersolution2 = DSolve[{ieqϵ01==0,ψ1[0] ==0},ψ1[ξ],ξ]

ieqϵ01 = ieqϵ0/.h'[ϵ ξ]-> D[Exp[-ϵ ξ],ξ]

innersolution2 = DSolve[{ieqϵ01==0,ψ1[0] ==0},ψ1[ξ],ξ]

psi1 =( ψ1[ξ]/.innersolution2[[1]])/.C[1]->C[2]

psi11 = PowerExpand[psi1,Assumptions->{ξ>0}]

Let’s get another term of the inner solution.

ieqϵ1 = Coefficient[inner5,ϵ,1]

innersolution3 = DSolve[{ieqϵ1==0,ψ2[0] ==0},ψ2[ξ],ξ]

psi2 =( ψ2[ξ]/.innersolution3[[1]])/.C[1]->C[3]

psi21 = PowerExpand[psi2,Assumptions->{ξ>0}]

So three terms of inner solution are:

inner3 = psi0+ ϵ psi11  + ϵ^2 psi21

Matching

We have an inner and an outer solution.  You can see that the outer solution has no unknown constants, but the inner solution does.  We find these by matching the inner and outer solutions at each order.  That is, each term of the inner solition should match the outer solution as x->1 (ξ->Infinity) as well as this can be made to occur.

One way to do this is to re-scale the other solution to be in an intermediate variable and to do the same for inner variable with the inverse scaling.  I pick “1/2” as this intermediate scaling, ϵ^1/2

innerscaled = inner3/.ξ->η/ϵ^(1/2)

innerscaled1 =ExpandAll[innerscaled]

As it turns out, these will not match exactly.  The “exponentially small” terms, (e^(-(η/(Sqrt[ϵ])))) C[1] can be neglected.    So we set these to 0.

innerscaled2 =innerscaled1/.Exp[-η/Sqrt[ϵ]]->0

Four terms of the outer solution are:

outer3 =  perturbexpansion/.{Y0[x]->1-h[1]+h[x], Y1[x]->h'[1]-h'[x],Y2[x]->-h''[1]+h''[x],Y3[x]-> h'''[1]-h'''[x]}

outerscaled =outer3/.{x-> η ϵ^(1/2)}

As this is an expansion in small epsilon, we can expand h around ϵ=0,

outerscaled1 = Normal[Series[outerscaled,{ϵ,0,2}]]

The first term is easy to compare.  Just subtract the two expansions.  We get the coefficient for ψ0[ξ]->(e^(-ξ)) (-1+(e^ξ)) C[1]

Limit[innerscaled2-outerscaled1,ϵ->0,Assumptions->{η>0}]

Solve[%==0,C[1]]

Let’s get the next one...

term2 = ExpandAll[(innerscaled2-outerscaled1)/.C[1]->1+h[0]-h[1]]

Coefficient[term2,ϵ,1]

Solve[%==0,C[2]]

term3 = ExpandAll[((innerscaled2-outerscaled1)/.C[1]->1+h[0]-h[1])/.C[2]->-h'[0]+h'[1]]

Coefficient[%,ϵ,2]

I think it is safe to say that

Solve[C[3]-h''[0]+h''[1]==0,C[3]]

This is for:  (1/2) (e^(-ξ)) (-2 C[1]+2 (e^ξ) C[1]-2 (e^ξ) ξ h''[0]+(e^ξ) (ξ^2) h''[0])

voila the answer:

solution1 = ((innerexpansion/.{α->1})/.ψ0[ξ]->(E^(-ξ)) (-1+(E^ξ)) C[1])/.{C[1]->1+h[0]-h[1]}

solution2 = ((solution1/.ψ1[ξ]->(E^(-ξ)) (-C[2]+(E^ξ) C[2]+(E^ξ) ξ h'[0]))/.C[2]->-h'[0]+h'[1])

solution3 = ((solution2/.ψ2[ξ]->(1/2) (E^(-ξ)) (-2 C[3]+2 (E^ξ) C[3]-2 (E^ξ) ξ h''[0]+(E^ξ) (ξ^2) h''[0]))/.C[3]->h''[0]-h''[1])

solution4=solution3/.ψ3[ξ]->0

We change back to “x” as the indepedent variable.

solution5 = solution4/.ξ->x/ϵ

We can now replace h[x] and its derivatives.

plott1 = solution5/.{h[0] -> Exp[-0], h[1] -> Exp[-1],h'[0]->-1, h'[1]-> - Exp[-1], h''[0] ->1, h''[1]-> Exp[-1]}

innerplot = Plot[plott1/.ϵ->.1,{x,0,1},PlotRange->All, PlotLegends->{"ϵ=.1"}]

The outer solution is

1-h[1]+h[x];

outerplot =Plot[1 - Exp[-1] + Exp[-x],{x,0,1}, PlotStyle->Red,PlotLegends->{"outer solution"}]

Compare the full inner expansion and the other solution

Show[innerplot,outerplot]

Here is the analytical solution

exactplot = Plot[Re[why/.{ϵ->.1}],{x,0,1},PlotRange->All,PlotStyle->Dashed, PlotLegends->{"analytical solution"}]

Show[innerplot,outerplot,exactplot]

innerplot = Plot[plott1/.ϵ->.001,{x,0,1},PlotRange->All, PlotLegends->{"ϵ=.001"}]

The outer solution is

1-h[1]+h[x];

outerplot =Plot[1 - Exp[-1] + Exp[-x],{x,0,1}, PlotStyle->Red,PlotLegends->{"outer solution"}]

Compare the full inner expansion and the other solution

Show[innerplot,outerplot]

Here is the analytical solution

exactplot = Plot[Re[why/.{ϵ->.001}],{x,0,1},PlotRange->All,PlotStyle->Dashed, PlotLegends->{"analytical solution"}]

Show[innerplot,outerplot,exactplot]

innerplot = Plot[plott1/.ϵ->.00001,{x,0,.0001},PlotRange->All, PlotLegends->{"ϵ=.0001"}]

The outer solution is

1-h[1]+h[x];

outerplot =Plot[1 - Exp[-1] + Exp[-x],{x,0,.0001}, PlotStyle->Red,PlotLegends->{"outer solution"}]

Compare the full inner expansion and the other solution

Show[innerplot,outerplot]

Here is the analytical solution

exactplot = Plot[Re[why/.{ϵ->.00001}],{x,0,.0001},PlotRange->All,PlotStyle->Dashed, PlotLegends->{"analytical solution"}]

Show[innerplot,outerplot,exactplot]

Let’s see what the numerical solution does at this epsilon

AbsoluteTiming[ans =NDSolve[{ .00001 y''[x] +  y'[x]  ==- Exp[-x], y[0]==0,y[1] ==1},y[x],x]]

numerical = Plot[y[x]/.ans,{x,0,.0001},PlotRange->All,PlotStyle->Green,PlotLegends->{"numerical"}]

Show[innerplot,outerplot,exactplot,numerical]