ODE: Asymptotic inner expansions

I am given the following ODE $$\epsilon \frac{d^2y}{dx^2}+x\frac{dy}{dx}+y=0 \; , \qquad 1<x<2$$

with the BVPs $$y(1)=0 \qquad , \qquad y(2)=1$$

and where $0<\epsilon \ll1$.

I am to use the method of matching inner and outer expansions to solve this.

By letting $y(x,\epsilon) \sim y_0(x)+y_1(x)\epsilon + y_2(x)\epsilon ^2+\cdots$, I have found the leading-order term of the outer expansion to be $2/x$ (and this satisfies only $y(2)=1$).

In computing the inner expansion, I let $x=\epsilon ^a z$ so that we now have $$\epsilon ^{1-2a}\frac{d^2Y}{dz^2}+z\frac{dY}{dz}+Y=0$$

where $Y(z,\epsilon)=y(x,\epsilon)$.

In order to match the orders, we need to choose $a$ such that $1-2a=0$, thus $a=\frac 12$. Plugging this in, we have $$\frac{d^2Y}{dz^2}+z\frac{dY}{dz}+Y=0$$

Letting $Y(z,\epsilon) \sim Y_0(z)+\cdots$, the equation at leading order is simply $$\frac{d^2Y_0}{dz^2}+z\frac{dY_0}{dz}+Y_0=0$$ and solving this gives $$Y_0=e^{-x^2/2}\Bigl(A+B \; \text{erfi} \bigl(\frac{x}{\sqrt 2} \bigr) \Bigr)$$

which doesn't seem to be right. Have I done something wrong? Or, ignoring what I have done, what is the right way to approach this?

• Your change of variable in the inner expansion is wrong, you want x close to 1 not close to 0.
– Ian
Commented Mar 1, 2018 at 20:52
• I dont understand. So what should I have chosen instead? Commented Mar 1, 2018 at 21:04
• $1+\epsilon^a z$.
– Ian
Commented Mar 1, 2018 at 21:04
• Oh, I see. And just to clarify, we want $x$ close to $1$ because the boundary condition is given at $1$? Commented Mar 1, 2018 at 21:06
• Yes, that is correct.
– Ian
Commented Mar 1, 2018 at 21:08

For the fast movements of the solution you in general look for local coordinates $x=x_0+ϵ^aX$ which gives for $Y(X)=y(x_0+ϵ^aX)$ $$ϵ^{1-2a}Y''+(x_0+ϵ^aX)ϵ^{-a}Y'+Y=0\iff ϵ^{1-a}Y''+(x_0+ϵ^aX)Y'+ϵ^{a}Y=0$$ which gives $a=1$ for the width exponent of the inner solution (as the exceptional point $x_0=0$ is not inside the integration interval). As $x_0>0$ the solution to the equation of the leading terms $Y''+x_0Y'=0$ is $Y(X)=Ce^{-x_0X}+D$. As it is not bounded in direction $X\to-\infty$ the boundary layer has to occur at the left interval end $x_0=1$ with first order approximation of the solution as $$y(x)=\frac2x -2e^{-(x-1)/ϵ}$$

$ϵ=0.01$, blue = numerical solution, green = boundary layer approximation

$\newcommand{\bbx}[1]{\,\bbox[15px,border:1px groove navy]{\displaystyle{#1}}\,} \newcommand{\braces}[1]{\left\lbrace\,{#1}\,\right\rbrace} \newcommand{\bracks}[1]{\left\lbrack\,{#1}\,\right\rbrack} \newcommand{\dd}{\mathrm{d}} \newcommand{\ds}[1]{\displaystyle{#1}} \newcommand{\expo}[1]{\,\mathrm{e}^{#1}\,} \newcommand{\ic}{\mathrm{i}} \newcommand{\mc}[1]{\mathcal{#1}} \newcommand{\mrm}[1]{\mathrm{#1}} \newcommand{\pars}[1]{\left(\,{#1}\,\right)} \newcommand{\partiald}[3][]{\frac{\partial^{#1} #2}{\partial #3^{#1}}} \newcommand{\root}[2][]{\,\sqrt[#1]{\,{#2}\,}\,} \newcommand{\totald}[3][]{\frac{\mathrm{d}^{#1} #2}{\mathrm{d} #3^{#1}}} \newcommand{\verts}[1]{\left\vert\,{#1}\,\right\vert}$

Lets $\ds{\,\mrm{y}\pars{x} \equiv \exp\pars{\mrm{S}\pars{x,\delta} \over \delta}\,,\quad \delta \equiv \root{\epsilon}.\quad}$

$$\mbox{Then,}\quad\left\{\begin{array}{rcl} \ds{\mrm{y}'\pars{x}} & \ds{=} & \ds{\mrm{y}\pars{x}{\mrm{S}'\pars{x,\delta} \over \delta}} \\ \ds{\mrm{y}''\pars{x}} & \ds{=} & \ds{\mrm{y}\pars{x}{\bracks{\mrm{S}'\pars{x,\delta}}^{\,2} \over \delta^{2}} + \mrm{y}\pars{x}{\mrm{S}''\pars{x,\delta} \over \delta}} \end{array}\right.$$ and

\begin{align} &\delta^{2}\braces{\mrm{y}\pars{x}{\bracks{\mrm{S}'\pars{x,\delta}}^{\,2} \over \delta^{2}} + \mrm{y}\pars{x}{\mrm{S}''\pars{x,\delta} \over \delta}} + x\bracks{\mrm{y}\pars{x}{\mrm{S}'\pars{x,\delta} \over \delta}} + \mrm{y}\pars{x} = 0 \\[5mm] & \bbx{\bracks{\mrm{S}'\pars{x,\delta}}^{\,2}\,\delta + \mrm{S}''\pars{x,\delta}\delta^{2} + x\,\mrm{S}'\pars{x,\delta} + \delta = 0} \end{align}

Now, you can expand $\ds{\,\mrm{S}\pars{x,\delta}}$ in powers of $\ds{\delta}$.

Namely, $$\mrm{S}\pars{x,\delta} = \mrm{S}_{0}\pars{x} + \mrm{S}_{1}\pars{x}\delta + \mrm{S}_{2}\pars{x}\delta^{2} + \cdots$$ $$\left\{\begin{array}{rcl} \ds{\mrm{S}_{0}'\pars{x}} & \ds{=} & \ds{0} \\[1mm] \ds{x\,\mrm{S}_{1}'\pars{x} + 1} & \ds{=} & \ds{0} \\[1mm] \ds{\mrm{S}_{2}'\pars{x}} & \ds{=} & \ds{0} \\[1mm] \ds{\mrm{S}_{1}'^{2}\pars{x} + \mrm{S}_{1}''\pars{x} + x\,\mrm{S}_{3}'\pars{x}} & \ds{=} & \ds{0} \\[1mm] \ds{\vdots} & \ds{\vdots} & \ds{\vdots} \end{array}\right.$$

This approach is the WKB Approximation.

• With this approach the $1/δ$ term reduces to a constant that is a constant factor that gets renormalized away. No boundary layer happens. I think you should have $δ=ϵ$. Commented Mar 2, 2018 at 6:19
• @LutzL $\color{red}{No}$. See the discussion in the above link. I didn't elaborate further because it's too long. So, the best option was to show a useful link. Commented Mar 2, 2018 at 19:48
• Moved previous comment to additional answer for more space and graphics. Commented Mar 3, 2018 at 8:24

Doing the WKB approximation right

As Felix Marin cited, one starts by representing the solution as $y(x)=\exp\left[\frac1δS(x;δ)\right]$ where later on $S(x;δ)$ is developed as power series in $δ$, $S(x;δ)=S_0(x)+S_1(x)δ+S_2(x)δ^2+\dots$. Inserting the derivatives $y'=\frac1δyS'$ and $y''=\frac1δyS''+\frac1{δ^2}yS'^2$ into the differential equation gives $$0=ϵ\left(\frac1δyS''+\frac1{δ^2}yS'^2\right)+\frac1δxyS'+y =\frac1δy\left(\fracϵδS'^2+ϵS''+xS' +δ\right).$$ In the last form one sees that the dominant terms in the last factor are $$\fracϵδS'^{\,2}+xS'.$$ To get meaningful results from the approximation method these terms have to be balanced in magnitude. This happens for $δ=ϵ$ with the resulting equations for the perturbation series \begin{align} &&(S'^2+xS')+δ(S''+1)&=0\\[1em]\hline &δ^0:& S_0' (S_0' +x)&=0 \\ &δ^1:& 2S_0' S_1' +xS_1' +S_0''+1&=0 \\ &δ^2:& S_1'^2+2S_0'S_2'+xS_2'+S_1'' &=0 \\ &&&\vdots \end{align} The first equation has two different solutions for $S_0'$ that give rise to two basis solutions for $y$. With $S_0'=-x$ the following equations result in $S_1' =S_2' =0,...$. The second solution has $S_0' =0$, thus no fast term, and leads to $S_1' =-\frac1x$, $S_2' =-\frac2{x^3},...$.

Integration of the terms leads to the assembled general solution approximation $$y(x)=A\exp\left(-\frac{x^2-1}{2ϵ}+O(ϵ^2)\right)+B\frac1x\exp\left(\frac{ϵ}{x^2}-ϵ+O(ϵ^2)\right),$$ where integration constants were chosen to give $0=y(1)=A+B$. The second boundary condition then gives $$1=y(2)=A\exp\left(-\frac{3}{2ϵ}\right)+\frac{B}2\exp\left(-\frac34ϵ\right) \implies -A=B=2\exp\left(\frac34ϵ\right)$$ ignoring the very small term $\exp(-\frac{3}{2ϵ})$.

Comparing against the numerical solution confirms a good fit even for $ϵ$ as large as $0.1$, while moving away visibly for larger $ϵ$. $ϵ=0.1$, blue = numerical solution, red = $O(ϵ^2)$ WKB approximation

$ϵ=0.2$, blue = numerical solution, red = $O(ϵ^2)$ WKB approximation