Resolvendo ODE singular não linear com SciPy odeint / ODEPACK

Quero resolver a equação 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}$$

com as condições iniciais

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

usando SciPyodeint() , mas, como pode ser visto, a equação é singular na origem. A documentação afirma que usa ODEPACK.

Eu já conheço a série de potências da solução em um bairro de ( ref ):$\xi = 0$

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

Tentei configurar tcrita np.array([0.0]), mas não funcionou: eu recebo um aviso sobre valores inválidos e, em seguida, a minha solução é tudo NaN. Devo integrar a partir de 0,01, talvez? Ou existe alguma outra solução?

Tudo bem, essa resposta é um tiro no escuro, mas aqui vai.

Primeiro, transforme o ODE de segunda ordem em um sistema de dois ODEs. Deixei

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

onde os pontos no topo das funções correspondem à diferenciação em relação à variável independente (neste caso, ).$\xi$

Em seguida, a ODE implícita de segunda ordem

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

pode ser expresso como o ODE explícito de primeira ordem

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

A princípio, parece que não podemos avaliar o lado direito desse sistema ODE explícito em , como exige um integrador numérico. Se existe uma solução para esse sistema, ele deve ser diferenciável. Supondo que exista uma solução, considere o limite do lado direito como .ξ → $\xi = 0$$\xi \rightarrow 0$

Primeiro, sabemos que

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

porque assumimos que existe uma solução, então é diferenciável, o que significa que deve ser contínuo. O limite de uma função contínua em um ponto é o seu valor nesse ponto, e sabemos o valor de porque é uma condição inicial. φ 2 ( 0 $\varphi_{2}$$\varphi_2(0)$

Também sabemos que

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

por razões semelhantes; assumimos que é diferenciável, por isso é contínuo e porque é uma condição inicial. φ 1 ( 0 ) = $\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 a regra de L'Hôpital no formulário indeterminado .$0/0$

Para prosseguir, precisamos fazer outra suposição: é contínuo em . Então segue-se queξ=$\dot{\varphi}_{2}$$\xi = 0$

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

Revisitando a ODE de primeira ordem e avaliando o lado direito em , podemos ver que temos:$\xi = 0$

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

a partir do qual se segue que .$\dot{\varphi}_{2}(0) = 1/3$

Usando esta análise, você pode conectar uma ifinstrução que retorne esses valores da função do lado direito em , o que deve levar você além da singularidade. Dito isto, essa análise requer algumas suposições sobre a continuidade que podem ou não se sustentar; portanto, tome a solução resultante com um grão de sal.$\xi = 0$

Se você quiser mais opções para o seu solucionador de ODE, consulte o pacote Assimulo que implementa as ligações ao pacote CVODE (e o RADAU e alguns integradores simples para esse assunto).