Lösen nichtlinearer singulärer ODE mit SciPy-Odeint / ODEPACK

Ich möchte die Lane-Emden-Isothermengleichung lösen [PDF, Gl. 15.2.9]

$$\frac{d^2 \!\psi}{d \xi^2} + \frac{2}{\xi} \frac{d \psi}{d \xi} = e^{-\psi}$$

mit den Anfangsbedingungen

$$\psi(\xi = 0) = 0 \quad \left.\frac{d\psi}{d \xi}\right|_{\xi = 0} = 0$$

mit SciPy,odeint() aber wie zu sehen ist, ist die Gleichung am Ursprung singulär. In der Dokumentation wird angegeben, dass ODEPACK verwendet wird.

Ich kenne bereits die Potenzreihen der Lösung in einer Nachbarschaft von $\xi = 0$ ( ref ):

$$\psi(\xi) \simeq \frac{\xi^2}{6} - \frac{\xi^4}{120} + \frac{\xi^6}{1890}$$

Ich habe versucht Einstellung tcritzu np.array([0.0]), aber nicht Arbeit: Ich habe eine Warnung über ungültige Werte erhalten und dann meine Lösung ist alles NaN. Sollte ich vielleicht ab 0,01 integrieren? Oder gibt es eine andere Lösung?

Gut, diese Antwort ist ein Schuss in die Dunkelheit, aber hier geht es weiter.

Transformieren Sie zunächst die ODE zweiter Ordnung in ein System aus zwei ODEs. Lassen

$$\begin{align} \varphi_{1} &= \psi, \\ \varphi_{2} &= \dot{\psi}, \end{align}$$

$\xi$

Dann die implizite ODE zweiter Ordnung

$$\begin{align} \ddot{\psi}(\xi) + 2\xi^{-1}\dot{\psi}(\xi) &= e^{-\psi(\xi)} \\ \psi(0) &= 0 \\ \dot{\psi}(0) &= 0 \end{align}$$

kann als explizite ODE erster Ordnung ausgedrückt werden

$$\begin{align} \dot{\varphi}_{1}(\xi) &= \varphi_{2}(\xi) \\ \dot{\varphi}_{2}(\xi) &= -2\xi^{-1}\varphi_{2}(\xi) + e^{-\varphi_{1}(\xi)} \\ \varphi_{1}(0) &= 0 \\ \varphi_{2}(0) &= 0. \end{align}$$

$\xi = 0$$\xi \rightarrow 0$

Erstens wissen wir das

$$\begin{align} \lim_{\xi \rightarrow 0} \varphi_{2}(\xi) = 0, \end{align}$$

Da wir davon ausgegangen sind, dass eine Lösung existiert, ist differenzierbar, was bedeutet, dass sie kontinuierlich sein muss. Die Grenze einer stetigen Funktion an einem Punkt ist ihr Wert an diesem Punkt, und wir kennen den Wert von da es sich um eine Anfangsbedingung handelt.$\varphi_{2}$$\varphi_2(0)$

Das wissen wir auch

$$\begin{align} \lim_{\xi \rightarrow 0} e^{-\varphi_{1}(\xi)} = 1 \end{align}$$

aus ähnlichen Gründen; Wir haben angenommen, dass differenzierbar ist, also stetig, und da es sich um eine Anfangsbedingung handelt.$\varphi_{1}$$\varphi_{1}(0) = 0$

Schließlich,

$$\begin{align} \lim_{\xi \rightarrow 0} \frac{-2\varphi_{2}(\xi)}{\xi} = \lim_{\xi \rightarrow 0} -2\dot{\varphi}_{2}(\xi), \end{align}$$

unter Verwendung der L'Hôpital-Regel auf der unbestimmten Form .$0/0$

Um weiter fortzufahren, müssen wir eine andere Annahme treffen: ist stetig bei . Dann folgt daraus$\dot{\varphi}_{2}$$\xi = 0$

$$\begin{align} \lim_{\xi \rightarrow 0} -2\dot{\varphi}_{2}(\xi) = -2\dot{\varphi_{2}}(0). \end{align}$$

Wenn wir die ODE erster Ordnung erneut betrachten und die rechte Seite bei auswerten , können wir sehen, dass wir:$\xi = 0$

$$\begin{align} \dot{\varphi}_{1}(0) = 0 \\ \dot{\varphi}_{2}(0) = -2\dot{\varphi}_{2}(0) + 1, \end{align}$$

woraus folgt, dass .$\dot{\varphi}_{2}(0) = 1/3$

Mit dieser Analyse können Sie eine ifAnweisung einfügen, die diese Werte der Funktion auf der rechten Seite bei zurückgibt , wodurch Sie die Singularität überwinden sollten. Diese Analyse erfordert jedoch einige Annahmen über die Kontinuität, die möglicherweise zutreffen oder nicht. Nehmen Sie die resultierende Lösung mit einem Salzkorn.$\xi = 0$

Wenn Sie mehr Optionen für Ihren ODE-Solver wünschen, schauen Sie sich das Assimulo- Paket an, das Bindungen an das CVODE-Paket (und RADAU und einige einfache Integratoren) implementiert.