scipy.optimize.root not converging and RuntimeWarning

Question

I am trying to solve the following problem: $$ \frac{d^2y}{dx^2}=\sinh(y) $$ Where the boundary conditions are: $y(0)=-1$, and $ \frac{dy(x\rightarrow \infty)}{dx}=0 $. Through central difference approximation, this can be simplified to: $$ y(x_{i+1})-2\cdot y(x_i)+y(x_{i-1})=h^2\cdot \sinh(y(x_i)) $$ Where $h$ is the step in x, which approaches 0.

I am trying to solve it in scipy.optimize.root (not my choice, have been told to do so...). And the code I have come up with is:

import numpy as np
import matplotlib.pyplot as plt
import scipy
N=500
x = np.linspace(0,5,N+1)
h = x[1]-x[0]
kb = 1.38e-23 # m^2 kg s^-2 K^-1
T = 298 # K
ec = 1.6e-19 # Coulombs
p0 = -1/(kb*T/ec)

# Prealocate matrices
Ah1 = np.zeros([N+2,N+2]) # expand Ah to include "ficticious solution at N+2 (to account for differential BC)
f1 = np.zeros([N+2,1]) # same as in Ah
# initial boundary conditions:
Ah1[0,0] = 1 
f1[0] = p0
# final BC: central approximation for phi'(x)=(phi(x_i+1)-phi(x_i-1))/2h
# applied to our BC: phi(x_N+1)-phi(x_N-1)=0, assuming h-->0, phi(x_N+1)-phi(x_N-1)=phi(x_N+2)-phi(x_N)=0
# f[N+2] stays as 0 (due to RHS of equation showed in last line)
Ah1[-1,-3:] = [1, 0, -1]
# Finite approximation: y(x_i+1)-2*y(x_i)+y(x_i-1)=h**2*sinh(y)
for i in range(1,N+1): # 1 to N+1 as Ah solution has been expanded 
    Ah1[i,i-1] = 1
    Ah1[i,i] = -2
    Ah1[i,i+1] = 1
def idk(y):
    f1 = h**2*np.sinh(y)
    return np.reshape(np.abs(np.dot(Ah1,y)-f1),(f1.size,))

f0 = np.ones((f1.size,1))*p0
f0[-1] = 0
sol = scipy.optimize.root(idk,f0)
phi_f1 = sol.x*(kb*T/ec)
phi_f1 = phi_f1[:-1] # subtract last ("ficticious") term
print("--Ah1 matrix, dimensions:", np.shape(Ah1))
print(Ah1)
plt.figure(dpi=100)
plt.plot(x,phi_f1)
plt.xlabel(r' $ \frac{x}{\lambda} $ ', fontsize=13)
plt.ylabel(r'$\phi$ (V)', fontsize=13)
plt.title(r'Non-linear Poisson-Boltzman $(\phi_0=-1 V)$')
plt.grid()

But I get the following error and it doesn't converge: RuntimeWarning: overflow encountered in sinh Any idea why that could be? Thanks!

EDIT: Jacobian added:

AJ = Ah1[:]
def jac(y):
    for i in range(1,N+1):
        AJ[i,i] = -2-h**2*np.cosh(y[i])
    return AJ

New call for root optimizer root(idk,f0,jac=jac).

Answer 1

The only stationary point is the saddle point at $(y,y')=(0,0)$. Linearization around that point gives that approximately at infinity $y''=y$. The stable solution satisfies $y'=-y$, which would also give a more realistic right boundary condition.

In this case one can integrate once to $$ y'^2=2\cosh(y)+C. $$ For a solution converging to the saddle point one has $C=-2$, so $$ y'^2=4\sinh^2(y/2). $$ The stable solution follows again the negative feedback path $$ y'=-2\sinh(y/2), $$ this would give an improved boundary condition.

---

Most problematic is

def idk(y):
    f1 = h**2*np.sinh(y)
    return np.reshape(np.abs(np.dot(Ah1,y)-f1),(f1.size,))

You do not set the specific right sides for the boundary conditions. You would have to insert in the middle

f1[0] = p0
f1[-1] = +4*h*np.sinh(y[-2]/2) # or +2*h*y[-2],    or 0 if you insist

Note that by construction of Ah, the last equation is y[-3]-y[-1]=..., that is, y(b-h)-y(b+h)=...

---

The error message means that the root-finding procedure has diverged to large values. It takes not much to cause overflow in the exponential.

---

The equation with negative feedback can be solved, set $u=e^{y/2}$, then $$ 2u'=uy'=-(u^2-1)\implies \frac{u-1}{u+1}=-Ce^{-x},~~ u=\frac{1-Ce^{-x}}{1+Ce^{-x}} \\~\\ y(x)=2\Bigl(\ln(1-Ce^{-x})-\ln(1+Ce^{-x})\Bigr),~~ C=\frac{1-e^{y_0/2}}{1+e^{y_0/2}}>0. $$