# First-term approximation for singular perturbation of ODE (with two turning points)

I'm reading "Introduction to Perturbation Methods" by Mark Holmes, and I came across an exercise that I don't know how to approach. As I decided to independently read this book, I have no friends/classmates/teachers with whom I can discuss this subject.

The exercise reads as follows (ex. 2.47):

Consider the problem $$\epsilon y'' -(x-a)(x-b)y' - x(y-1) = 0, \text{ for } 0<x<1$$ with the boundary conditions $y(0)=-2$ and $y(1)=2$. The numerical solution in the case $a=1/4$, $b=3/4$ and $\epsilon = 10^{-4}$ is shown below. Based on this information, derive a first-term approximation of the solution for arbitrary $0<a<b<1$.

I'll be more than happy if you just tell me/direct me on how to solve the problem for the particular case $a=1/4$, $b=3/4$.

-EDIT- Looking at the graph and after skimming through a lot of books, I think there's an interior boundary layer at $x=1/4$ and another boundary layer at $x=1$.

I found the outer solution, valid in $0 \leq x<1/4$, to be: $$\tag{1} y_{\text{out}}=1 - 3 \left( \frac{3}{3-4x} \right)^{\frac{3}{2}} \sqrt{1-4x}$$

I still have to find/guess the boundary layer thickness at $x=1/4$ and compute an inner solution. Then, it somehow seems that this boundary layer will "connect" the outer solution (1) to the solution $y(x)=1$, which will then be connected to the boundary layer at $x=1$. I'd appreciate any help/insight on how to proceed.

You don't say where you went wrong, so I here is an outline of what I think the exercise requires you to do. This is exactly as it is described in the chapter on matched asymptotic expansions.

You already have the outer solution $1-A|\frac14-x|^{1/2}|\frac34-x|^{-3/2}$. Stylistic note: it is easier to write down the general solution in the form $$-3\left|\frac{a-x}{a}\right|^{a/\alpha} \left|\frac{b-x}{b}\right|^{b/\alpha},$$ where $\alpha=b-a$.

If the first boundary layer, is at $x=a$, then you introduce a new variable $\bar x=(x-\frac14)/\epsilon^\alpha$, and substitute $$y = 1 + Y(\bar x),$$ together with the change of independent variable to $\bar x$. (Adding $1$ makes the equation linear.) Using $\frac{d}{dx}=\epsilon^{-\alpha}\frac{d}{d\bar x}$ and $x=a + \bar x\epsilon^\alpha$, the resulting equation (in $\bar x$) is $$\epsilon^{1-2\alpha}Y'' -\bar x \epsilon^{\alpha}(a+\bar x\epsilon^{\alpha}-b)\epsilon^{-\alpha}Y'-(a+\bar x\epsilon^\alpha)Y = 0.$$ Requiring $\epsilon^{1-2\alpha}Y''$ to be the same order in $\epsilon$ as the other highest order terms, which are $\epsilon^0$, gives the condition $\alpha=\frac12$, and the equation to leading order: $$Y'' + \frac12\bar x Y'-\frac14 Y = 0,$$ with $\alpha=\frac12$, and you can solve this using Kummer functions.

The book also helpfully gives asymptotics for Kummer functions. The solution is of the form $$Y = A_0 M(-1/4, 1/2, -\eta) + B_0 \bar x M(1, 3/2, -\eta),$$ where $\eta = -\frac14\bar x^2$. Then as $\bar x\to\pm\infty$, the asymptotic form of the inner solution is $$Y \sim \sqrt\pi\left(\frac{A_0}{\Gamma(\frac34)}\pm \frac{B_0}{\Gamma(\frac54)}\right) \eta^{1/4}.$$ This asymptotic needs to be matched to $Y\sim y$ ($\bar x\to-\infty$), and to $Y\sim 0$ ($\bar x\to+\infty$). When matching the solution to $y$, you need to introduce an intermediate variable $x_\rho = (x-\frac14)/\epsilon^\rho$, with $0<\rho<\frac12$, so that $x=\frac14+x_\rho\epsilon^\rho$, and $\bar x = x_\rho\epsilon^{\rho-1/2}$. These two matches give you $A_0$ and $B_0$ in terms of the constant in the outer solution $A$ and also $\epsilon$.

It is similar for the boundary layer at $x=1$. Expanding in $\bar x=(x-1)/\epsilon$, and using $y=1+Y(\bar x)$ gives $$Y''-\frac3{16}Y' = 0,$$ which has solution $Y=A_0\bar x+B_0e^{3\bar x/16}$, which requires $A_0=0$ and $B_0=1$.

Here is what this looks like (red is the numerical solution, orange, green and blue are analytical):

This is the code used to generate the plot. Mathematica's NDSolve won't solve the boundary value problem itself, because it uses a shooting method, and the problem posed that way is incredibly numerically ill-conditioned (e.g., solutions with $y(1)=2$ and with $y(2)=2$ are almost the same), so this uses a simple relaxation scheme. It's not very efficient.

Clear[s1];
s1[\[Epsilon]_?NumericQ, a_: 1/4, b_: 3/4, n_: 1000] :=

Module[{eq, y, x, k, i, h, eqs, arr, sol, yd},
h = 1/n;
eq = \[Epsilon] y''[x] - (x - a) (x - b) y'[x] - x (y[x] - 1) == 0;
eq = eq /. {
y''[x] -> yd[2, k]
, y'[x] -> yd[1, k]
, y[x] -> Subscript[y, k], x -> Subscript[x, k]};
(*yd[1,k_Integer]/;0<=k-2&&k+2<=n:=1/h{1/12,-2/3,0,2/3,-1/12}.Table[
Subscript[y, i],{i,k-2,k+2}];*)

yd[1, k_Integer] := (Subscript[y, k + 1] - Subscript[y, k - 1])/(2 h);
(*yd[2,k_Integer]/;0<=k-2&&k+2<=n:=1/h^2{-1/12,4/3,-5/2,4/3,-1/
12}.Table[Subscript[y, i],{i,k-2,k+2}];*)

yd[2, k_Integer] := (Subscript[y, k + 1] - 2 Subscript[y, k] + Subscript[y, k - 1])/h^2;
eqs = Table[eq, {k, n - 1}];
eqs = eqs /. {Subscript[y, k_ /; k <= 0] -> -2, Subscript[y, k_ /; k >= n] -> 2, Subscript[x, k_] :> k/N[n]};
arr = CoefficientArrays[eqs, Table[Subscript[y, k], {k, n - 1}]];
sol = Join[{-2}, LinearSolve[arr[[2]], -arr[[1]]], {2}];
Print@
Show[ListPlot[sol, DataRange -> {0, 1}, PlotRange -> {-2, 2},
PlotMarkers -> None, Joined -> True, PlotStyle -> Red],
Plot[1 -
3 Abs[(a - x)/a]^(a/(b - a)) Abs[(b - x)/b]^(-b/(b - a)), {x,
0, (1 + 1.*^-2) a}, PerformanceGoal -> "Quality",
Exclusions -> None, PlotStyle -> Green],
Plot[1 -
9/2 Sqrt[3/\[Pi]] \[Epsilon]^(1/4)
Gamma[3/4] Hypergeometric1F1[-1/4,
1/2, -1/4 ((x - 1/4)/Sqrt[\[Epsilon]])^2] +
9/2 Sqrt[3/\[Pi]] \[Epsilon]^(1/4)
Gamma[5/4] (x - 1/4)/Sqrt[\[Epsilon]] Hypergeometric1F1[1/4,
3/2, -1/4 ((x - 1/4)/Sqrt[\[Epsilon]])^2]
, {x, 0, b}, Exclusions -> None, PlotStyle -> Orange]
,
Plot[1 + E^((1 - a) (1 - b) (x - 1)/\[Epsilon]), {x, b, 1.1},
PlotStyle -> Blue, Exclusions -> None]
]
]
s1[0.1]
s1[0.01]
s1[0.001]
s1[0.0001]

• thanks a lot for the outline. I'm curious how to justify the (in)existence of the boundary layers. I'd expect a boundary layer at $x=3/4$ since $p(3/4)=0$, where $p(x) = -(x-1/4)(x-3/4)$. also, how did you determine the thickness of the boundary layers? the number of possible combinations to consider for dominant balance is huge! – Maya Jun 20 '14 at 23:05
• @MayaLda. I didn't try to find where the layers are; I just took them to be at $x=a$ and $x=1$ from the numerical solution (he says that in the exercise as well). Near $x=3/4$, $y$ is almost constant $1$, so $y'=0$ and a vanishing coefficient doesn't matter; $x=\frac34$ isn't special. What do you mean by thickness? You introduce $\bar x=(x-a)/\epsilon^\alpha$, substitute using $\frac{d}{dx}=\epsilon^{-\alpha}\frac{d}{d\bar x}$, and require all the leading terms to be the same order in $\epsilon$, giving an equation in $\alpha$. – Kirill Jun 20 '14 at 23:27
• nice looking plot! the thickness of a boundary layer is the value of $\alpha$. then, how do you determine which terms should be leading? supposedly you can do it using dominant balance, considering all the combinations between all the terms in the resulting ODE. do you have another technique, or is it just a good guess? PS: I wasn't able to plot the numerical solution using NDSolve in Mathematica, but you were! how did you do it? can you please post the code? thanks in advance. – Maya Jun 20 '14 at 23:31
• @MayaLda. I reworded the computation of $\alpha$. There are terms of order $O(\epsilon^0)$, $O(\epsilon^{1-2\alpha})$ and $O(\epsilon^\alpha)=o(\epsilon^0)$. You set the first two types to be the same order ($\alpha=\frac12$) and ignore the third type. I wasn't able to get NDSolve to work either. – Kirill Jun 21 '14 at 0:45
• amazing work! thanks for your detailed answer and patience. – Maya Jun 21 '14 at 0:58