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

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.
• Welcome to SciComp.SE. If you add a description of your attempt it might be easier to help you. Aug 24 at 11:17
• Is my code not clear? I used Python to solve the differential equation. Aug 24 at 13:27
• I didn't check the code since I don't know what you are trying to do. Also, my guess is that you want to solve a boundary value problem but you are not listing the boundary conditions. If it is a boundary value problem, I suppose that you would like to use a shooting method since odeint is targeted to initial value problems. Aug 24 at 16:06
• @Linda, I have edited your question a bit. I suggest that you add the details mentioned above using MathJax. Aug 24 at 17:04
• Interesting problem, $\phi$=const is also a solution matching the boundary condition $\phi '$=0 at $\xi$=0; but it would not satisfy the BC at infinity; I've seen something like that. You cannot solve it by shooting the way the equation is written; but I would rewrite the equation as a 2nd order ODE and then solve it by shooting. Aug 25 at 5:29

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.

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

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.

• Hi @Futurologist, thank you for the answer. What do the first three functions in your Python code represent? I am a bit confused by those functions. Aug 25 at 18:04
• @Linda I added an edit to my question, explaining the role of the functions. You can replace them with a different more direct implementation. I personally try to have my code modulated and organized in smaller functions that can be corrected more easily. Aug 25 at 20:34