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.