Solving an ODE beyond existence. What's happening?

Question

As an example for an ODE course I used the ODE $$ y' = \frac{y}{x} + \frac{1}{\cos(\tfrac{y}{x})} $$ to illustrate domains of existence. Standard substitution $z=y/x$ turns the equation to $$ z' = \frac{1}{x\cos(z)} $$ which can be solved by separation of variables to obtain the solution $$ \phi(x) = x\arcsin(\ln(x/x_0) - \sin(y_0/x_0)) $$ with initial value $\phi(x_0)=y_0$ (and the solution ceases to exists at $x = x_0\exp(1-\sin(y_0/x_0))$).

To check that my solution is correct, I compared the solution with the output of some numerical solver (ode45 from MATLAB). Since I was to lazy to figure out the precise upper bound at which the solution ceases to exist, I just plugged in some upper bound and looked at the result. As expected, something weired happened beyond the point of existence. However, when comparing the numerical solution with the computed one (which is correct, by the way), I noticed something strange: Although the solution computed by the exact formula has a problem (when the argument of the $\arcsin$ becomes larger than one and from that point on MATLAB plot the real part and neglects the imaginary part of the solution), the "wrong exact solution" somehow follows that "wrong numerical solution" in a very crude, strange but not deniable way. Here is the plot for $x_0=1$, $y_0=0$:

(the red is the plot of $\phi$ and the blue is the output of ode45).

Why is that? What is going on with an ODE solver beyond existence? Why is it following the real part of the non-existing complex solution?

Answer 1

I'd say the key is when you say "the solution $\phi(x)$", because the solution to ($x_0=1, y_0=0$) $$ \sin(y/x) = \log x $$ is one of $$ y = x\big( \arcsin\log x + 2\pi k \big), \qquad y = x\big( \pi - \arcsin\log x + 2\pi k\big) $$ for an appropriate integer $k$. The ODE itself is well-defined beyond the location of singularity, unlike the solution that starts at $(1,0)$.

First, there is a slight bug in the stepping function — it should be able to correctly detect the singularity, but it doesn't. Instead, it steps just past the singularity.

Second, the singular points of the ODE in the $(x,y)$-plane are defined by $\cos(y/x)=0$, so in this case it is the line $y=\frac\pi2 x$ that matters. When the ODE solver steps "through" a singularity, it ends up at some value that is at $\tilde x$ just after the location of the singularity $x_{\mathrm{sing}}$, and at $\tilde y$ just above (on the other side of) the singular value $\frac\pi2 \tilde x$.

It isn't really correct that it follows "the real part of the complex solution". There is a well-defined real solution that starts at initial conditions $(\tilde x, \tilde y)$, and goes up (as $\arcsin$ above) when $\tilde y<\frac\pi2 \tilde x$, and goes down (as $\pi-\arcsin$ above) when $\tilde y>\frac\pi2\tilde x$. That's why the plot keeps crossing the singular line over and over again.

The way the plot tracks the line $y=\frac\pi2 x$ consists of individual arcs that are ODE solutions themselves. The locations of those arcs are random because they depend on the numerical errors in the ODE solver.

There might also be a separate issue with how it is plotted, because the plot sometimes doesn't cross the red line as expected. It could be that the plotting function isn't sampling the solution often enough, but I can't tell from just the plot. On the other hand it could be that the solution instead crosses the line in some intermediate solver steps only.