# scipy.optimize.root not converging and RuntimeWarning

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!

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).

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.$$

• I completely forgot to reassign my boundary conditions in the function, now my output makes a lot more sense, thanks! I still get the runtime warning though, could that be causing sol.success to remain as false?
– HWIK
Commented Oct 15, 2022 at 17:49
• In my version, out of 4614 calls to the function I get in call 507 suddenly the value 2122500.386 in place 501, falling until in call 530 unproblematic values 695.49411 and below are reached. This is up to the internals of the root finder, I have no idea how this happens. Commented Oct 15, 2022 at 18:25
• I tried defining a Jacobian to see if it would improve at all, but I'm getting odd results, would you mind giving me your thoughts on it
– HWIK
Commented Oct 16, 2022 at 12:23
• You would have to define AJ as global if you want to change this out-of-scope variable inside the function. You also need to treat again the (first and) last line as special case. Commented Oct 16, 2022 at 12:35
• Taking zero for the initial guess, and the corrected Jacobian, I get to "The iteration is not making good progress, as measured by the improvement from the last five Jacobian evaluations." after 141 function calls, but without any overflow problems. The problem appears to be that the solution, exact as numerical, makes a jump of more than $70\%$ towards zero between y[0] and y[1], this does not lead to numerical stability. This gets a little better with N=700, but also markedly slower. Commented Oct 16, 2022 at 16:01