2nd order differential equation coupled to integro-differential equation in python

Question

I'm trying to solve the following equations numerically in python

$$\begin{align} 12\pi\int_0^\infty drf(r)\phi(r)r^4&=E\\ f(r)-\frac{1}{2\mu}\bigg(\frac{d^2\phi(r)}{dr^2}+\frac{2}{r}\frac{d\phi(r)}{dr} \bigg)+m_\pi\phi(r)&=E\phi(r) \end{align} $$

where $f(r)=S\exp(-r^2/b^2)$, while $\mu$ and $m_π$ are constants.

I have tried rewriting the 2nd order differential equation into two first order equations and solving the system of equations with root finding.

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
from scipy.integrate import trapz
from scipy.optimize import root
b=1    
S=20    
m=135   
mn = 939.5  
mu = m*mn/(mn+m)

g = (2*mu)

def f(r): 
    return S*np.exp(-r**2/b**2)

def diff(phi,r,E):
    return (phi[1],g*(-E+m)*phi[0]-2/r*phi[1]+g*f(r))

phi0 = [b/m,b/m] #Initial

def phi_fun(E):
    rs = np.linspace(1e-5,50,1000)
    ys = odeint(lambda phi,r: diff(phi,r,E), phi0, rs)
    integral = 12*np.pi*trapz(ys[:,0]*f(rs)*rs**4,rs)
    return integral - E

E_true = root(phi_fun, -2).x

rs = np.linspace(1e-5,50,1000)
ys = odeint(lambda phi,r: diff(phi,r,E_true), phi0, rs)

phi_true = ys[:,0]
plt.plot(rs, phi_true,linewidth=2,label=r'$\phi(r)$')
print("Minimum found at E =",E_true)

This does not produce a result I would expect and maybe I have done something wrong. I was thinking about using IDEsolver but I'm not sure how to implement this. Any suggestions? Thanks in advance!

Answer 1

The next transformation step would see the system reformulated as a boundary value problem with the additional equations and conditions replacing the integral $$ I'(r)=12πf(r)ϕ(r)r^4,~~I(0)=0,~~ I(\infty)=E $$ scipy provides integrate.solve_bvp, which has a mechanism to handle free parameters like $E$.

Now the equation proceeds from a singular point to an infinite point. At both ends approximations are required if the solution is to stay finite.

• Close to zero the equation is approximately an Euler-Cauchy equation $$r^2ϕ''(r)+2rϕ'(r)=0$$ which has basis solutions $1$ and $r^{-1}$. The second one is to be avoided for finite solutions, so that $ϕ'(a)=0$ is the condition to demand for $a\approx0$, so use that as the left boundary condition with $a\gtrapprox0$.

• For large $r$, other terms in the equation become negligible, the equation is approximately $$ -ϕ''(r)+2μ(m_π-E)ϕ(r)=0 $$ You seem to be expecting a negative value for $E$. In that case, the basis solutions are $e^{\pm\sqrt{2μ(m_π+|E|)}r}$. To avoid a diverging behavior also on this end, the positive sign has to be avoided. The basis solution with negative exponents satisfies $$ ϕ'(r)+\sqrt{2μ(m_π+|E|)}ϕ(r)=0. $$ This is then also a suitable boundary condition on the right boundary.

def sys(r,u,E):
    y,v,I = u
    dy = v
    dv = g*(-E+m)*y-2/r*v+g*f(r)
    dI = f(r)*r**4*y
    return dy,dv,dI

def bc(ua, ub,E):
    ya,va,Ia = ua
    yb,vb,Ib = ub
    return va, vb+(g*(m+abs(E)))**0.5*yb, Ia, Ib-E

r = np.logspace(-5,0,20)*5
u = [ 0*r,0*r,E*r/r[-1]]
E = -2

res = solve_bvp(sys,bc,r,u,p=[E],tol=1e-4)
print(res.message,", E: ",res.p[0])

plt.plot(res.x,res.y.T,'-+'); 
plt.grid(); plt.legend(r"$\phi$ $\phi'$ $I$".split()); 
plt.show()

This works well, reports

The algorithm converged to the desired accuracy. , E:  -0.3472387364434603

and gives a plot