Peculiar second order autonomous DE

Question

I discovered this equation that is baffling me how to solve. The physical nature of it is important, so i will discuss. Suppose you have a tank with a pipe coming out of it. There is an item stuck in this pipe. To remove this item, you begin pressurizing the tank with air until the item breaks free. Once it breaks free, you stop increasing the pressure in the tank and you wish to know the equations of motion of the item as it goes down the pipe. If the pipe is long enough, the item will eventually get stuck again due to friction, and that's the end of it.

For sheer simplicity, take all constants to be one. Then we eventually end up with the equation (using newton force law $f=ma$ and ideal gas law $p_1v_1=p_2v_2$)

$x"+x'=\frac{1}{x}$

Where x is one plus the distance the item has travelled down the pipe.

The interesting part is if you do the standard $y=x'$ and $x"=y\frac{dy}{dx}$ then you get the equation

$y\frac{dy}{dx}+y=\frac{1}{x}$

With initial condition $y=0$ when $x=1$. This equation is elegant because solving it allows one to find the velocity $y$ as a function of distance down the pipe ($x$ shifted by 1). But it can't be solved numerically using first order method! The very first iteration has no solution, try it! However, the original second order equation can be solved numerically, but you have the awkward time independent variable to interpert! This phenomenon holds even we don't assume all constants are one!

Strange, isn't it? Not really sure what to make of this.... Can anyone tell me what's wrong here?

Answer 1

$$\frac{d^2x}{dt^2}+\frac{dx}{dt}=\frac{1}{x} \tag 1$$ I would like to add a few considerations to the pertinent answer from Claude Leibovici.

The second equation : $$y\,\frac{dy}{dx}+y=\frac{1}{x}\quad \text{with}\quad y(1)=0 \tag 2$$ poses a particular problem because $\quad 0\,\frac{dy}{dx}+0=\frac{1}{1}\quad$ implies $\quad\frac{dy}{dx}=\infty\quad$ at $(x=1\:,\:y=0)$.

The difficulty for numerical solving comes from the infinite slope at the starting point.

All this is due to the transformation of the first equation to the second with the change $\quad \frac{dx}{dt}=y(x) \quad\to\quad \frac{d^2x}{dt^2}=y\frac{dy}{dx}\quad$ which introduces the product $0$ by $\infty$ at the starting point.

On the other hand, there is no difficulty with Eq.$(1)$ where $\frac{dx}{dt}=0$ and $\frac{d^2x}{dt^2}=1$ at the starting point.

If we really want to use Eq.$(2)$ for numerical calculus, we must chose another starting point than $(x=1\:,\:y=0)$ , not far, but not exactly at this point.

In order to find the initial condition at a point $(x=x_0\:,\:y=y_0)$ we have to solve Eq.$(2)$ with $y$ small compared to $\frac{1}{x}\simeq 1$. In first approximation $\quad y\,\frac{dy}{dx}\simeq\frac{1}{x}\quad\to\quad y\simeq \sqrt{2\ln(x)} $

So, for numerical calculus, let $\quad\begin{cases}x_0=1+\epsilon \\y_0=\sqrt{2\ln(x_0)}\simeq \sqrt{2\epsilon} \end{cases}$

$$y\,\frac{dy}{dx}+y=\frac{1}{x}\quad \text{with}\quad y(x_0)=y_0 $$

I cannot say what will be the deviation due to the use of a first approximate only for the starting point. This could be investigated by comparison to the numerical calculus from Eq.$(1)$.

This is in the same vein than Ian's comment.

Answer 2

This is probably not an answer.

If we consider the second equation $$y\,\frac{dy}{dx}+y=\frac{1}{x}\qquad \text{with}\qquad y(1)=0$$ you end needing to solve for $y$ the following equation $$x=\frac{2\, e^{\frac{1}{2} (x+y)^2}}{\sqrt{2 \pi } \left(\text{erfi}\left(\frac{x+y}{\sqrt{2}}\right)-\text{erfi}\left(\frac{1}{\sqrt{2}}\right)\right)+2 \sqrt{e}}\tag 1$$ which is not the most pleasant one.

You can visualize $y(x)$ using a contour plot.

Being unable to produce correct plots (I am almost blind), I give you below a few values $$\left( \begin{array}{cc} x & y \\ 1.0 & 0.000000 \\ 1.1 & 0.371958 \\ 1.2 & 0.475887 \\ 1.3 & 0.533912 \\ 1.4 & 0.568128 \\ 1.5 & 0.587369 \\ 1.6 & 0.596296 \\ 1.7 & 0.597751 \\ 1.8 & 0.593637 \\ 1.9 & 0.585306 \\ 2.0 & 0.573773 \\ 2.1 & 0.559823 \\ 2.2 & 0.544083 \\ 2.3 & 0.527065 \\ 2.4 & 0.509195 \\ 2.5 & 0.490831 \\ 3.0 & 0.400678 \\ 3.5 & 0.326873 \\ 4.0 & 0.273999 \\ 4.5 & 0.237019 \\ 5.0 & 0.209913 \\ 6.0 & 0.171945 \\ 7.0 & 0.146052 \\ 8.0 & 0.127090 \\ 9.0 & 0.112557 \\ 10.0 & 0.101043 \end{array} \right)$$

Edit

If we put constant coefficients and consider $$ay\,\frac{dy}{dx}+by=\frac{1}{x}\qquad \text{with}\qquad y(c)=0$$ the equation will be $$x=\frac{2 \sqrt{a} c\, e^{\frac{(a y+b x)^2}{2 a}}}{2 \sqrt{a} e^{\frac{b^2 c^2}{2 a}}+\sqrt{2 \pi } b c \left(\text{erfi}\left(\frac{a y+b x}{\sqrt{2a} }\right)-\text{erfi}\left(\frac{b c}{\sqrt{2a} }\right)\right)}\tag 2$$

Update after JJacquelin's answer

Using the implicit function theorem for the computation of $\frac{dy}{dx}$ from equation $(2)$ and using Taylor expansion around $x=c$, we get around this point $$\frac{dy}{dx}=\frac{1}{b c (x-c)}-\frac{b}{2 a}+O\left(x-c\right)$$ which could give some further ideas for the selection of the initial step size.

In any manner, my personal choice would be to stay with equation $(2)$ to be solved using Newton.