Resolviendo ODE singular no lineal con SciPy odeint / ODEPACK

Quiero resolver la ecuación isotérmica de Lane-Emden [PDF, eq. 15.2.9]

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

con las condiciones iniciales

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

usando SciPyodeint() pero, como se puede ver, la ecuación es singular en el origen. La documentación indica que usa ODEPACK.

Ya conozco la serie de potencia de la solución en un entorno de ( ref ):$\xi = 0$

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

Intenté fijar tcrita np.array([0.0]), pero no funcionó: recibo una advertencia acerca de los valores válidos y luego mi solución es todo NaN. ¿Debo integrar a partir de 0.01 quizás? ¿O hay alguna otra solución?

Muy bien, esta respuesta es un tiro en la oscuridad, pero aquí va.

Primero, transforme el ODE de segundo orden en un sistema de dos ODE. Dejar

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

donde los puntos sobre las funciones corresponden a la diferenciación con respecto a la variable independiente (en este caso, ).$\xi$

Entonces el ODE implícito de segundo orden

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

puede expresarse como el ODE explícito de primer orden

$$\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}$$

Al principio, parecería que no podemos evaluar el lado derecho de este sistema ODE explícito en , como requiere un integrador numérico. Si existe una solución para este sistema, entonces debe ser diferenciable. Suponiendo que existe una solución, tome el límite del lado derecho como .ξ → $\xi = 0$$\xi \rightarrow 0$

Primero, sabemos que

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

porque asumimos que existe una solución, entonces es diferenciable, lo que significa que debe ser continua. El límite de una función continua en un punto es su valor en ese punto, y sabemos el valor de porque es una condición inicial.$\varphi_{2}$$\varphi_2(0)$

También sabemos que

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

por razones similares; hemos asumido que es diferenciable, por lo que es continuo, y porque es una condición inicial.$\varphi_{1}$$\varphi_{1}(0) = 0$

Finalmente,

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

usando la regla de L'Hôpital en la forma indeterminada .$0/0$

Para continuar, tenemos que hacer otra suposición: es continua en . Entonces se sigue que$\dot{\varphi}_{2}$$\xi = 0$

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

Revisando el ODE de primer orden y evaluando el lado derecho en , podemos ver que tenemos:$\xi = 0$

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

de donde se deduce que .$\dot{\varphi}_{2}(0) = 1/3$

Usando este análisis, podría conectar una ifdeclaración que devuelva estos valores de la función del lado derecho en , lo que debería superar la singularidad. Dicho esto, este análisis requiere un par de suposiciones sobre la continuidad que pueden o no mantenerse, así que tome la solución resultante con un grano de sal.$\xi = 0$

Si desea más opciones para su solucionador de ODE, eche un vistazo al paquete Assimulo que implementa enlaces al paquete CVODE (y RADAU y algunos integradores simples para el caso).