How to solve the integral-like energy equation with Sagdeev potential numerically in Python?

Question

I am trying to numerically solve equation (6) of Lakhina 2021 in Python. The equation is

$$\frac{1}{2}\left(\frac{d \phi}{d\xi}\right)^2 + S(\phi, M) = 0\, .$$

The Sagdeev potential expression is given by (7).

What I want is to reproduce the potential profiles in Fig. 3 of Lakhina 2021.

The boundary conditions given in the paper are:

$\phi(0)_{M = 2.55} = 0.023$
$\phi(0)_{M = 2.57} = 0.037$
$\phi(0)_{M = 2.55} = 0.046$

In the code below, I first define a function for the first-order differential equation. Then, set boundary conditions for each mach number, $M$, and finally, I use odeint from the scipy.integrate module in Python to solve the boundary value problem. The plot of the solutions is shown in the last figure.

Here is my attempt, the Python code:

## Importing standard modules
from scipy.integrate import odeint
import numpy as np
import matplotlib.pyplot as plt 


## Reconnection jet plasma parameters
n1 = 0.74 
n2 = 0.26 
sig1 = 0.11 
sig2 = 0.07 
U1 = -1.72
U2 = 1.82 


# Function for Sagdeev potential equation 
def S(phi, M):
    s = (1 - np.exp(phi)) + n1/(6*np.sqrt(3*sig1))*((M - U1 + np.sqrt(3*sig1))**3 -
                                                    ((M - U1 + np.sqrt(3*sig1))**2 - 2*phi)**1.5 -
                                                    (M -  U1 - np.sqrt(3*sig1))**3 + 
                                                    ((M - U1 - np.sqrt(3*sig1))**2 - 2*phi)**1.5) + n2/(6*np.sqrt(3*sig2))*(
                                                    (M - U2 + np.sqrt(3*sig2))**3 -
                                                    ((M - U2 + np.sqrt(3*sig2))**2 - 2*phi)**1.5 -
                                                    (M - U2 - np.sqrt(3*sig2))**3 +
                                                    ((M - U2 - np.sqrt(3*sig2))**2 - 2*phi)**1.5) 
    return s

## Solving the ode

def model(phi, zeta, M):

S = (1 - np.exp(phi)) + n1/(6*np.sqrt(3*sig1))*((M - U1 + np.sqrt(3*sig1))**3 -
                                                ((M - U1 + np.sqrt(3*sig1))**2 - 2*phi)**1.5 -
                                                (M -  U1 - np.sqrt(3*sig1))**3 + 
                                                ((M - U1 - np.sqrt(3*sig1))**2 - 2*phi)**1.5) + n2/(6*np.sqrt(3*sig2))*(
                                                (M - U2 + np.sqrt(3*sig2))**3 -
                                                ((M - U2 + np.sqrt(3*sig2))**2 - 2*phi)**1.5 -
                                                (M - U2 - np.sqrt(3*sig2))**3 +
                                                ((M - U2 - np.sqrt(3*sig2))**2 - 2*phi)**1.5)  
dphi_dzeta = -np.sqrt(-2*S)


return dphi_dzeta


# Boundary conditions
phi0_M255 = 0.023 #For M = 2.55
phi0_M257 = 0.037 #For M = 2.57
phi0_M259 = 0.046 #For M = 2.59



phi_array = np.linspace(-0.01, 0.06, 1000)
zeta_array = np.linspace(-16, 16, 1000)

Phi = odeint(model, phi0, zeta_array, args = (2.57,))

## Plotting

plt.figure(2)
plt.axhline(0, color = 'k', lw = 1)
plt.axvline(0, color = 'k', lw = 1)
plt.plot(zeta_array, Phi, label = "M = 2.55")
plt.xlabel("$\zeta$")
plt.ylabel("S($\phi$, M)")
plt.legend()

Output:

May you please assist? I am really not sure where I am going wrong.

References

1. Lakhina, G. S., Singh, S. V., & Rubia, R. (2021). A mechanism for electrostatic solitary waves observed in the reconnection jet region of the Earth’s magnetotail. Advances in Space Research.

Answer 1

A simple approach here would be to use the shooting method, by integrating from $\xi$=0 to infinity (some large value) the ODE written as

$ \frac{d}{d\xi} \phi = \sqrt{-2 S(\phi)} $

Due to the symmetry we have the constraint $d_{\xi} \phi$=0 at $\xi$=0, and to find the right value of $\phi$ at $\xi=0$ we'd need to solve $S(\phi)$=0 which is not easy to do. But even if we find the right value $\phi(0)$, there are two solutions here that satisfy the condition $d_{\xi} \phi$=0 at $\xi$=0; one is $\phi$=const and the other one is the correct solution that dies out at $\xi \rightarrow \infty$. In other words, there are two integral curves that touch at $\xi$=0, and integrating the ODE from $\xi=0$ would keep us on the wrong integral curve.

To overcome this difficulty, let's manipulate the equation a little bit. Take the derivative of the equation

$ \frac{d}{d\xi} (\frac{1}{2} (\dot{\phi})^2 + S(\phi) = 0) $

That produces

$ \dot{\phi} \ddot{\phi} + S' \dot{\phi} = 0, $

where the dot stands for $d_{\xi}$ and prime stands for $d_{\phi}$.

Now we can divide by $\dot{\phi}$ to write it as

$ \ddot{\phi} = -S' $

Then introduce variable $\psi=\dot{\phi}$, and we have a system of ODEs amenable to solving by the shooting method

$ \dot{\phi} = \psi \\ \dot{\psi} = -S'(\phi) $

The initial conditions at $\xi$=0 are $\psi$=0 and $\phi$=$\phi_0$, where the parameter $\phi_0$ is found by shooting to satisfy the BC $\phi=0$ at infinity.

Answer 2

Quite hastily, I admit, I wrote the code below and generated the solutions. I did not bother to center them.

One issue you have with your code is that you are trying to integrate

$$\frac{d\phi}{d\xi} \, =\, \sqrt{- \, 2 \, S(\phi, M)}$$

and the square root may get in the way, because you have to choose whether it is the one above or this one, especially when the derivative becomes zero (and that's why you are not seeing the maximum when you integrate the problem):

$$\frac{d\phi}{d\xi} \, =\, - \,\sqrt{-\,2 \, S(\phi, M)}$$

Because I do not want to deal with this issue, I integrated the Hamiltonian equations instead

\begin{align} &\frac{d\phi}{d\xi} = p\\ &\frac{dp}{d\xi} = - \,\frac{\partial S}{\partial \phi}(\phi, M)\\ \end{align}

In addition to that, I used a python solver for initial conditions (I think that was another issue in your code), so I had to integrate forward in time $[0, \, 16]$ and backward in time $[-16, \, 0]$ and plot both solutions to get the solution from $[-16, \, 16]$

Finally, since the system is time-invariant, if you have a solution $\phi(\xi)$ then $\phi(\xi-\xi_0)$ is also a solution. So I assume that in the paper, the solutions were shifted to have their peaks at $\xi = 0$.

import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
 
def _grad_U(phi, A, a):
    return - 3 * np.sqrt( A**2 - 2*phi ).dot(a) + np.exp(phi)

def f(t, y, A, a):
    # a = [a1, -a1, a2, -a2]
    return np.array([ y[1], _grad_U(y[0], A, a) ])
   
def initialize(phi, A, a):
    S = (A**3 - np.sqrt( A**2 - 2*phi )**3).dot(a) + 1 - np.exp(phi)
    return np.sqrt(-2*S)

def set_parameters(M, U1, U2, sigma1, sigma2, n1, n2):
    sqrt_3sigma1 = np.sqrt(3*sigma1)
    sqrt_3sigma2 = np.sqrt(3*sigma2)
    A1  = M - U1 + sqrt_3sigma1
    A1_ = M - U1 - sqrt_3sigma1
    A2  = M - U2 + sqrt_3sigma2
    A2_ = M - U2 - sqrt_3sigma2
    a1  = n1 / (6*sqrt_3sigma1)
    a2  = n2 / (6*sqrt_3sigma2)
    return np.array([A1, A1_, A2, A2_]), np.array([a1, -a1, a2, -a2])    

def solve_equation(phi_in, M, U1, U2, sigma1, sigma2, n1, n2, zeta_left, zeta_right, n_zetas):
    A, a = set_parameters(M, U1, U2, sigma1, sigma2, n1, n2)
    v_in = initialize(phi_in, A, a)
    zeta_array = np.linspace(0, zeta_right, n_zetas)
    y_in = np.array([phi_in, v_in])
    sol  = solve_ivp(fun = lambda t, y :  f(t, y, A, a), t_span=[0, zeta_right], y0=y_in, t_eval=zeta_array, method='Radau')
    zeta_array = np.linspace(0, zeta_left, n_zetas)
    sol_ = solve_ivp(fun = lambda t, y : f(t, y, A, a), t_span=[0, zeta_left], y0=y_in, t_eval=zeta_array, method='Radau')
    return sol, sol_


## Reconnection jet plasma parameters
n1 = 0.74 
n2 = 0.26 
sig1 = 0.11 
sig2 = 0.07 
U1 = -1.72
U2 = 1.82 

# Boundary conditions
phi0_M255 = 0.023 #For M = 2.55
phi0_M257 = 0.037 #For M = 2.57
phi0_M259 = 0.046 #For M = 2.59


zeta_left =  -16
zeta_right = 16
n_zetas = 5000

sol_255, sol_255_ = solve_equation(phi0_M255, 2.55, U1, U2, sig1, sig2, n1, n2, 
                     zeta_left, zeta_right, n_zetas)

sol_257, sol_257_ = solve_equation(phi0_M257, 2.57, U1, U2, sig1, sig2, n1, n2, 
                     zeta_left, zeta_right, n_zetas)

sol_259, sol_259_ = solve_equation(phi0_M259, 2.59, U1, U2, sig1, sig2, n1, n2, 
                     zeta_left, zeta_right, n_zetas)

#phi_array = np.linspace(-0.01, 0.06, 1000)
#zeta_array = np.linspace(-16, 16, 1000)


## Plotting

plt.figure(2)
plt.axhline(0, color = 'k', lw = 1)
plt.axvline(0, color = 'k', lw = 1)
plt.plot(sol_255.t, sol_255.y[0,:], color = 'b', label = "M = 2.55")
plt.plot(sol_255_.t, sol_255_.y[0,:], color = 'b')
plt.plot(sol_257.t, sol_257.y[0,:], color = 'g', label = "M = 2.57")
plt.plot(sol_257_.t, sol_257_.y[0,:], color = 'g')
plt.plot(sol_259.t, sol_259.y[0,:], color = 'r', label = "M = 2.59")
plt.plot(sol_259_.t, sol_259_.y[0,:], color = 'r')
plt.xlabel("$\zeta$")
plt.ylabel("$\phi$")
plt.legend()

Edit. The first function, _grad_U() calculates $ - \, \frac{\partial S}{\partial \phi}(\phi, M)$, taking as input $\phi$, as well as the parameters featured in the formula of $S$. I wrote it in a vectorized form, that is why it looks shorter (it is easier to keep track of things and make corrections if necessary, plus numpy is faster with vectorized quantities) The second one, f(), takes in $y = (\phi, \, p)$ and outputs $\Big(p, \, - \, \frac{\partial S}{\partial \phi}\partial(\phi, M)\Big )$. Additional parameters are the time t, which is there only because solve_ivp() needs it, and the rest are the parameters for the calculation of the potential. Basically, f() is the right hand side of the differential equations I wrote above. The third function, initialize() takes in the initial state $\phi(0)$ and calculates the corresponding initial momentum $p(0)$ so that they satisfy the total energy equation: $$\frac{1}{2}p(0)^2 \, + \, S(\phi(0), M) \, =\, 0 $$ That is why, the latter function calculates $$p(0) = \sqrt{2 \, S(\phi(0), M)}$$
The forth function, set_parameters() takes in the original parameters that are included in the calculation of the potential $S$ and reorganzies them into the new sets of more convenient parameters A and a, that are fed to the two functions before it.