This system

eqs = {(1/2) Y'[x]^2 == (1 - Log[Y[x]^2]) Y[x]^2, Y[0] == 1}

is known to have a simple solution in terms of Gaussian functions, which can be checked analytically (more precisely, it has two Gaussian solutions, $Y_- = -e^{-x\left(x-\sqrt{2}\right)}$ and $Y_+ = -e^{-x \left(x+\sqrt{2}\right)}$, due to the symmetry $Y \to - Y$ of this ODE).

However, when I try to numerically solve $eqs$ and plot any of its solutions:

nsol = NDSolve[eqs, Y, {x, -5, 5}] // Flatten

Plot[{Y[x] /. nsol[[1]]} // Evaluate, {x, -4, 4}, PlotRange -> {{-3, 3}, All}]

Plot[{Y[x] /. nsol[[2]]} // Evaluate, {x, -4, 4}, PlotRange -> {{-3, 3}, All}]

it produces only half of each of the solutions. Following the warning message of NDSolve, I tried to increase the MaxSteps value but it doesn't work, and it freezes PC for MaxSteps = 10^7 and above (I tried MMA up to 11.2.0).

Any suggestions what options/methods I should use for NDSolve here?

Update 11/27: It has been pointed out by LutzL that the problem might lie in the point where the r.h.s. of my ODE hits a zero value. If it becomes slightly negative (e.g., due to numerical errors/fluctuations) then NDSolve tries to take a square root of a negative value and collapses. If so, how one can prevent this without trying to fit the (Gaussian) solution which is already known to you?

Update 11/28: It was suggested to increase the order of this system. Indeed, a system of 2nd-order ODE

eqs2 = {Y''[x] Y'[x] == -2*Log[Y[x]^2]*Y[x]*Y'[x], Y[0] == 1, Y'[0] == Sqrt[2]}

is equivalent to the original one, but now it doesn't contain powers of derivatives (you can also cancel Y'[x] out, it changes nothing). It also eliminates the issue of a square root of a negative small number mentioned in the previous update and some of Answers. Well, guess what? NDSolve can't correctly integrate $eqs2$ either. Instead of a single Gaussian, it produces some weird quasi-oscillating chain of Gaussians.
P.S. Looking at the existence of so many troubles with numerical solving of such a simple system, I am asking myself a question. How many sophisticated numerical results published in gazillions research articles are actually having hidden errors of a computational origin, let alone the deliberate counterfeiting?... Numerical packages and codes are pretty much black boxes nowadays, therefore I bet that if such errors do happen, 99.99% of reviewers would not be able to spot them, especially for those situations when models are conceptually novel and/or not backed up by analytical calculations or estimates. And we are talking about petabytes of codes and outputs currently used in all branches of science...

  • $\begingroup$ You should probably ask this question on mathematica.stackexchange.com instead. $\endgroup$
    – heropup
    Nov 24 '17 at 7:49
  • $\begingroup$ One issue I can see: due to the square of $Y'(x)^2$ this is not a regular ODE - it is multivalued. If $Y(x)$ is a solution then $-Y(x)$ is also a solution. $\endgroup$
    – Winther
    Nov 24 '17 at 14:49
  • $\begingroup$ Winther, thanks for the idea, but this Z-symmetry is not what causes the problem. This symmetry just leads to the two solutions of the system, being two Gaussians in this case. Mathematica sees both of them (the nsol is a two-item list) but produces/plots only "half" of each of them. $\endgroup$
    – Konstantin
    Nov 27 '17 at 8:22

What you demand of the numerical solver is impossible on several counts

1. Domain of ODE is bounded, solver may violate bound

The immediate problem producing your error picture is that your ODE is not defined outside $|y|\le e^{1/2}$. A numerical solver for $y'=f(x,y)$, especially those with internal steps of adaptive step size, may want to evaluate the ODE function $f$ outside the bounds of the visible solution, that is, with $y$ outside the range of the exact solution and with $x$ outside the given integration interval. There may be options to restrict the domain of $f$, such as enforcing positivity, that may extend to this situation. A stop-gap measure here is to extend the domain of $f$ by constant continuation as in
$$ y'^2=2\max(0,1-\ln(y^2))y^2. $$ The numerical solver will then stay on the constant solution once it is reached, and that is correct as that is an approximation of one of the possible exact solutions to this initial value problems.

1.-2. Full solution set of the implicit ODE

At the non-zero roots of the right side $$0=(1-\ln(y^2))\iff y=\pm e^{1/2},$$ the ODE is not Lipschitz, as the square root function as a vertical tangent. Thus the uniqueness of the Picard-Lindelöf theorem does not apply, non-constant solutions reach these values. One can insert segments of these constant solutions at will in any Gauss bell curve solution to piece together further exact solutions of the implicit ODE.

2. There is no "up'n'down" in order 1 scalar autonomous ODE

In any explicit autonomous ODE $y'=f(y)$ the qualitative behavior is determined by the roots of $f$ and the sign of $f$ in the intervals between the roots. If the sign is positive, then the solution is monotonically increasing, if negative, then monotonically decreasing. If the ODE is differentiable, then the roots of $f$ give the stationary, constant solutions, and no other solution can cross or even reach these. If the derivative is not bounded at some root, then the ODE is stiff at this point, which will give any numerical solver problems with solutions that approach that point. The kind of problems also depends on the solver, fixed-step solvers may step over that point incurring large quantitative and qualitative errors, adaptive-step solvers tend to stall by regulating the step size below the machine epsilon.

3. The sign of the square root is fixed

When preparing the implicit ODE for the numerical solver, the CAS has to make it explicit. Here this involves choosing the sign of the square root. Once that sign is chosen, it stays constant for the full integration process. One would have to introduce some kind of event handling and define suitable events to effect a sign change. How that would work with the given Mathematica syntax I have no idea.

And finally a kind of solution

A) Higher order derivatives can serve as a kind of memory

To select the right sign of the root depends on the behavior of the solution at previous points. If the term under the square root is not zero, then one would chose a positive sign if the solution was previously increasing and a negative sign if it was decreasing. Circumventing all speculation involving Taylor polynomials, let's just take the derivative of the original ODE $$ y'y''=-2yy'+(1-\ln(y^2))2yy'=-2\ln(y^2)yy',\ y(0)=1, \frac12y'(0)^2=1 $$ and excluding constant solutions, even constant segments in a solution, one can divide by $y'$ to find $$ y''=-2\ln(y^2)y,\quad y(0)=1, y'(0)=\pm\sqrt2 $$ Pythons scipy.integrate.odeint for this second order ODE

def odefunc2(y,t): return [ y[1], -2*y[0]*np.log(y[0]**2) ]
sol2 = odeint(odefunc2,[ 1.0, 2**0.5],x)
plt.plot(x,sol2[:,0],'-b',x,sol2[:,1],'.g'); plt.grid(); plt.show()

results in the curve

solution of second order ode

  • $\begingroup$ LutzL, thank you for your answer, it seems robust. However, I'm not sure how your advice to change my original equation to yours $y'^2=2\max(0,1-\ln(y^2))y^2$ can help to restore the missing half of the Gaussian. Your advice replaces the missing half with a constant, which is not exactly the original goal... $\endgroup$
    – Konstantin
    Nov 27 '17 at 15:02
  • $\begingroup$ No, it does not restore that particular solution. Analytically, you can construct a valid solution by breaking apart any Gaussian at the maximum, continuing constant for a while and then going down with the second half of the Gaussian. There is nothing to restore as the solution is valid, the one that you want is just one among an infinity of correct solutions. This is similar and locally equivalent to the standard example $y'=\sqrt{y}$ with $y(0)=0$. Set $y=e^u$ in your problem, then $u'^2=2-4u$, which is exactly that. $\endgroup$ Nov 27 '17 at 15:26
  • $\begingroup$ If it doesn't restore then I don't see much sense in your advice (I still like your observation about the rhs's sign flip). Of course, I can patchwork different solutions along x-axis at my will (matching them somehow to each other), but this is not what I need. I need a numerical algorithm which would fully reproduce the Gaussian solution, even if I do not know that such solution exists... $\endgroup$
    – Konstantin
    Nov 27 '17 at 16:16
  • $\begingroup$ Besides, I don't know why do you call $Y=\sqrt e$ a correct solution here. It does not satisfy a condition $Y(0) =1$ which is a part of $eqs$... $\endgroup$
    – Konstantin
    Nov 27 '17 at 16:33
  • $\begingroup$ You can patch, because you know all the choices and can chose. A numerical algorithm only sees, if at all, the slope go to zero. How should the numerical method chose the correct sign of the square root when that is based on global considerations? $\endgroup$ Nov 27 '17 at 16:39

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.