Low convergence of a perturbation solution to a non-linear ODE

Question

I have the following first-order non-linear ODE:

$$ \frac{1}{12}y' + \frac{1}{40} \epsilon \left( y' \right)^3 = -1 \tag{1}\label{1} $$

where $y=y(x), \enspace x\in \left[0,1\right], \enspace \epsilon \ll 1$ is the perturbation parameter, and Eq. \eqref{1} is subject to,

$$ y(1) = 0. \tag{2}\label{2} $$

I sought a perturbation solution in powers of $\epsilon$ expanding $y$ as follows,

$$ y = y_0 + y_1 \epsilon + y_2 \epsilon^2 + \cdots. \tag{3}\label{3} $$

Substituting \eqref{3} in \eqref{1} and \eqref{2} produces the following problems:

$ O(1): $ $$ \frac{1}{12}y_0' = -1 \\ y_0(1) = 0 $$

$ O(\epsilon): $ $$ \frac{1}{12}y_1' + \frac{1}{40} \left( y_0' \right)^3 = 0 \\ y_1(1) = 0 $$

$ O(\epsilon^2): $ $$ \frac{1}{12}y_2' + \frac{3}{40} \left( y_0' \right)^2 y_1' = 0 \\ y_2(1) = 0 $$

The corresponding solutions are:

$$ y_0 = 12(1-x), \tag{4}\label{4} $$

$$ y_1 = -\frac{2592}{5}(1-x), \tag{5}\label{5} $$

$$ y_2 = \frac{1679616}{25}(1-x). \tag{6}\label{6} $$

Then, inserting Eqs. \eqref{4}-\eqref{6} in \eqref{3} yields the perturbation solution up to $ O(\epsilon^2) $. Finally, the exact solution of \eqref{1} is presented to measure the accuracy of \eqref{3}, namely

$$ y = \frac{ \sqrt[3]{100} \epsilon - \sqrt[3]{10 \alpha^2} }{ 3 \epsilon \sqrt[3]{\alpha} } (1-x), \tag{7}\label{7} $$

where $ \alpha = \sqrt{ 2 \epsilon^3 \left( 1458 \epsilon + 5 \right) } - 54 \epsilon^2 $.

Now, the exact solution \eqref{7} is plotted with its corresponding numerical solution (via MATLAB ode15i function). Figure 1 shows a perfect agreement between the exact and the numerical solution.

Figure 1. Exact and numerical solutions

Here comes my problem (after all this wordiness). Figure 2 shows that the perturbation solution \eqref{3} only converges to the exact solution \eqref{7} for really small values of $\epsilon$, namely, $\epsilon \sim O(10^{-3})$. The title in each subplot shows the $\epsilon$ value and the relative error of the perturbation solution to the exact solution. This means that the contributions of $O(\epsilon)$ and higher-order terms are not important, and all the physics of the problem is driven by the leading-order solution $y_0$ \eqref{4}.

Figure 2. Perturbation and exact solutions

Questions:

1. Did I achieve a correct perturbation solution?
2. Is this non-linear equation not suitable for a regular expansion approach?
3. Why does the perturbation parameter need to be so small?

A final comment. I proposed $ y = y_0 + y_1 \epsilon^\beta + \cdots,$ and I got $ \beta = 1$.

Answer 1

Writing the solution as

$$ y(x) = \gamma(\epsilon) (1-x) \ \ \ \quad \mathrm(0) $$

we see that $\gamma$ must satisfy

$$ \gamma^3\frac{\epsilon}{40}+\frac{\gamma}{12} -1=0. \ \ \ \quad \mathrm(1) $$

This is a cubic equation so it always has one real solution which is

$$ g(\epsilon) = \frac{1}{3} \left(\frac{10^{2/3}}{\sqrt[3]{\sqrt{2} \sqrt{1458 \epsilon ^4+5 \epsilon ^3}-54 \epsilon ^2}}-\frac{\sqrt[3]{10} \sqrt[3]{\sqrt{2} \sqrt{1458 \epsilon ^4+5 \epsilon ^3}-54 \epsilon ^2}}{\epsilon }\right). $$

Now, the polynomial in (1) changes order when $\epsilon=0$ so you don't even have guarantee that the solution is continuous in $\epsilon$ at zero. Indeed $g(\epsilon)$ is complex for $\epsilon<0$. Computing the derivatives as $\epsilon \to 0^+$ one obtains

$$ g(\epsilon) = 12-\frac{2592 \epsilon }{5}+\frac{1679616 \epsilon ^2}{25}+O\left(\epsilon ^{3}\right), \ \ \ \quad \mathrm(3) $$

so your calculation is correct. However the coefficients $a_n$ of this expansion defined as

$$ g(\epsilon) = \sum_{n=0} a_n \epsilon^n \ \ \ \quad \mathrm(4) $$

grow exponentially $a_n \sim b a^n$ (with $a\approx 186.9$, $b\approx3.14$) which means that the series (4) diverges for $\epsilon a >1$. I think this explains what you are observing. You could make it even more evident by plotting the exact form of $g(\epsilon)$ together with one of its Taylor approximants.