7
$\begingroup$

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

enter image description here

(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?

$\endgroup$
  • 1
    $\begingroup$ Perhaps I am missing something in your substitution but i get: $$\left(zx\right)'=z+xz'=z+\frac{1}{\cos z}\rightarrow z'=\frac{1}{x\cos z}$$ Dont know if it will help with your question... $\endgroup$ – nluigi Dec 8 '15 at 11:12
  • $\begingroup$ @nluigi Right, thanks! Just a typo as I think the rest of the question is ok (modulo further typos). $\endgroup$ – Dirk Dec 8 '15 at 11:21
6
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ You get an upvote for the explanation of the spurious arcs. However, I don't get the first paragraph. For any initial value there is a unique local solution and only one $k$ is appropriate (for these initial values it is $k=0$). I think the blue plot plots really all sample points of the solution, i.e. the $x$ variable has all sampling points that the stepping function produced. So the last thing I ask myself is, why the real part of the solution lies of the singularity… $\endgroup$ – Dirk Dec 9 '15 at 8:49
  • $\begingroup$ @Dirk $\Re(\arcsin x)=\frac\pi2$ for all $x\geq 1$, so I think it's a property specific to $\arcsin$ and this equation, not something that would hold in general. $\endgroup$ – Kirill Dec 9 '15 at 21:04

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.