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


  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$.

  • $\begingroup$ Are you sure you wrote (1) correctly? It is not really an ODE as written. $\endgroup$
    – copper.hat
    May 17, 2022 at 17:51
  • $\begingroup$ Since you have the exact solution, have you tried expanding it in $\epsilon$? $\endgroup$
    – lcv
    May 17, 2022 at 17:59
  • $\begingroup$ @copper.hat: It's a non-linear first-order ODE written implicit, $f(y',x)=0$. We can solve for $y'$ and obtain the explicit representation, $y' = g(x)$. But, since it's a third-grade polynomial equation in $y'$, the right-hand side of that rearranging deems to be ugly. $\endgroup$ May 17, 2022 at 19:04
  • $\begingroup$ @lcv: Good point. I'll do it, but the road is dreadful. $\endgroup$ May 17, 2022 at 19:07
  • 1
    $\begingroup$ The coefficients of the perturbation series grow very fast. $ϵ$ must at least be within the radius of convergence. The exact expression has a singularity at $-5/1458$, which is a bound on the radius of convergence. $\endgroup$ May 17, 2022 at 19:12

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.

  • $\begingroup$ Thanks for your comprehensive answer. Just a final question. In Eqs. (0) and (1), is $\gamma$ the same as your later $g(\epsilon)$? $\endgroup$ May 19, 2022 at 4:48
  • $\begingroup$ Yes it's (essentially) the same. $g$ is one solution of the polynomial whose variable I denoted with $\gamma$. $\gamma_1$ would probably have been a better notation but I was lazy about typing. $\endgroup$
    – lcv
    May 19, 2022 at 4:55

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .