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

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

• Are you sure you wrote (1) correctly? It is not really an ODE as written. May 17, 2022 at 17:51
• Since you have the exact solution, have you tried expanding it in $\epsilon$?
– lcv
May 17, 2022 at 17:59
• @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. May 17, 2022 at 19:04
• @lcv: Good point. I'll do it, but the road is dreadful. May 17, 2022 at 19:07
• 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. May 17, 2022 at 19:12

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.

• Thanks for your comprehensive answer. Just a final question. In Eqs. (0) and (1), is $\gamma$ the same as your later $g(\epsilon)$? May 19, 2022 at 4:48
• 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.
– lcv
May 19, 2022 at 4:55