NDSolve löst diese gewöhnliche Differentialgleichung nur „halbwegs“

Was Sie vom numerischen Löser verlangen, ist in mehrfacher Hinsicht unmöglich

1. Der Bereich von ODE ist begrenzt, der Solver kann die Grenze verletzen

Das unmittelbare Problem, das Ihr Fehlerbild erzeugt, besteht darin, dass Ihre ODE nicht außerhalb definiert ist$|y|\le e^{1/2}$. Ein numerischer Löser für$y'=f(x,y)$, insbesondere diejenigen mit internen Schritten mit adaptiver Schrittgröße, möchten möglicherweise die ODE-Funktion auswerten$f$außerhalb der Grenzen der sichtbaren Lösung, das heißt mit$y$außerhalb des Bereichs der exakten Lösung und mit$x$außerhalb des angegebenen Integrationsintervalls. Möglicherweise gibt es Optionen zum Einschränken der Domäne von$f$, wie z. B. das Erzwingen von Positivität, die sich auf diese Situation erstrecken können. Eine Notlösung ist hier die Erweiterung des Geltungsbereichs von$f$durch ständige Fortsetzung wie in

$$y'^2=2\max(0,1-\ln(y^2))y^2.$$ Der numerische Löser bleibt dann bei der konstanten Lösung, sobald sie erreicht ist, und das ist richtig, da dies eine Annäherung an eine der möglichen exakten Lösungen für dieses Anfangswertproblem ist.

1.-2. Vollständiger Lösungssatz der impliziten ODE

An den Nullwurzeln der rechten Seite

$$0=(1-\ln(y^2))\iff y=\pm e^{1/2},$$ die ODE ist nicht Lipschitz, da die Quadratwurzel als vertikale Tangente funktioniert. Damit gilt die Eindeutigkeit des Satzes von Picard-Lindelöf nicht, nicht konstante Lösungen erreichen diese Werte. Man kann beliebig Segmente dieser konstanten Lösungen in jede Gauß-Glockenkurven-Lösung einfügen, um weitere exakte Lösungen der impliziten ODE zusammenzusetzen.

2. Es gibt kein "up'n'down" in Ordnung 1 skalare autonome ODE

In jeder expliziten autonomen ODE$y'=f(y)$das qualitative Verhalten wird durch die Wurzeln bestimmt$f$und das Zeichen von$f$in den Zwischenräumen zwischen den Wurzeln. Bei positivem Vorzeichen ist die Lösung monoton steigend, bei negativem Vorzeichen monoton fallend. Wenn die ODE differenzierbar ist, dann sind die Wurzeln von$f$geben die stationären, konstanten Lösungen, und keine andere Lösung kann diese kreuzen oder auch nur erreichen. Wenn die Ableitung nicht an einer Wurzel begrenzt ist, ist die ODE an diesem Punkt steif, was jedem numerischen Löser Probleme mit Lösungen gibt, die sich diesem Punkt nähern. Die Art der Probleme hängt auch vom Löser ab, Löser mit fester Schrittweite können diesen Punkt überschreiten und große quantitative und qualitative Fehler verursachen, Löser mit adaptiver Schrittweite neigen zum Stillstand, indem sie die Schrittgröße unter das Maschinen-Epsilon regulieren.

3. Das Vorzeichen der Quadratwurzel ist fest

Beim Vorbereiten der impliziten ODE für den numerischen Löser muss der CAS sie explizit machen. Hier geht es darum, das Vorzeichen der Quadratwurzel zu wählen. Sobald dieses Vorzeichen gewählt ist, bleibt es für den vollständigen Integrationsprozess konstant. Man müsste eine Art Ereignisbehandlung einführen und geeignete Ereignisse definieren, um einen Vorzeichenwechsel zu bewirken. Wie das mit der gegebenen Mathematica-Syntax funktionieren würde, habe ich keine Ahnung.

Und schließlich eine Art Lösung

A) Ableitungen höherer Ordnung können als eine Art Gedächtnis dienen

Die Wahl des richtigen Vorzeichens der Wurzel hängt vom Verhalten der Lösung an den vorherigen Punkten ab. Wenn der Term unter der Quadratwurzel nicht Null ist, würde man ein positives Vorzeichen wählen, wenn die Lösung zuvor steigend war, und ein negatives Vorzeichen, wenn sie fallend war. Um alle Spekulationen mit Taylor-Polynomen zu umgehen, nehmen wir einfach die Ableitung der ursprünglichen ODE

$$y'y''=-2yy'+(1-\ln(y^2))2yy'=-2\ln(y^2)yy',\ y(0)=1, \frac12y'(0)^2=1$$ und unter Ausschluss konstanter Lösungen, sogar konstanter Segmente in einer Lösung, kann man durch dividieren $y'$finden $$y''=-2\ln(y^2)y,\quad y(0)=1, y'(0)=\pm\sqrt2$$ Pythons scipy.integrate.odeintfür diese ODE zweiter Ordnung

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

Ergebnisse in der Kurve