It is necessary to improve accuracy in fsolve to find multiple roots

Question

I use this code to get the zeros of a nonlinear function. Most likely, the function should have 1 or 3 zeros

import numpy as np
import matplotlib.pylab as plt
from scipy.optimize import fsolve

[a, b, c] = [5, 10, 0]

def func(x):
    return -(x+a) + b / (1 + np.exp(-(x + c)))

x = np.linspace(-10, 10, 1000)

print(fsolve(func, [-10, 0,  10]))
plt.plot(x, func(x))
plt.show()

In this case, the code gives 3 expected roots without any problems. But with c = -1.5, the code skips the root, and with c = -3 it finds a nonexistent root.

I want to calculate the roots for many combinations of parameters, so manually changing the seeds is not a practical solution.

I appreciate any solution, trick or tip.

Answer 1

Between two fixed points (i.e. df/dx=0) you have one or zero root. In your case, analytically you can calculate two fixed points:

[-c + log(1/(b - sqrt(b*(b - 4)) - 2)) + log(2),
 -c + log(1/(b + sqrt(b*(b - 4)) - 2)) + log(2)]

, , . Sympy . sy.nsolve() :

import sympy as sy

a, b, c, x = sy.symbols("a, b, c, x", real=True)

# The function:
f = -(x+a) + b / (1 + sy.exp(-(x + c)))
df = f.diff(x)  # calculate f' = df/dx
xxs = sy.solve(df, x)  # Solving for f' = 0 gives two solutions

# numerical values:
pp = {a: 5, b: 10, c: .5}  # values for a, b, c
fpp = f.subs(pp)
xxs_pp = [xpr.subs(pp).evalf() for xpr in xxs]  # numerical stationary points
xxs_pp.sort()  # in ascending order

# resulting intervals:
xx_low = [-1e9,      xxs_pp[0], xxs_pp[1]]
xx_hig = [xxs_pp[0], xxs_pp[1],       1e9]

# calculate roots for each interval:
xx0 = []
for xl_, xh_ in zip(xx_low, xx_hig):
    try:
        x0 = sy.nsolve(fpp, (xl_, xh_), solver="bisect")  # calculate zero
    except ValueError:  # no solution found
        continue
    xx0.append(x0)

print("The zeros are:")
print(xx0)

sy.plot(fpp)  # plot function

Answer 2

. , , . , , () , () () . , Numpy , .

[a, b, c] = [5, 10, -1.5]

def func(x):
    return -(x+a) + b / (1 + np.exp(-(x + c)))

polyfit poly1d func (-10<x<10) f_poly 10.

x_range = np.linspace(-10,10,100)
y_range = func(x_range)

pfit = np.polyfit(x_range,y_range,10)

f_poly = np.poly1d(pfit)

, f_poly func. , . , , fsolve

roots = np.roots(pfit)
roots

([- 10.4551 + 1.4893j, -10.4551-1.4893j, 11.0027 + 0.j,          8.6679 + 2.482j, 8.6679-2.482j, -5.7568 + 3.2928j,         -5.7568-3.2928j, -4.9269 + 0.j, 4.7486 + 0.j, 2.9158 + 0.j])

, Numpy 10 . [-10,10]. :

x0 = roots[np.where(np.logical_and(np.logical_and(roots.imag==0, roots.real>-10), roots.real<10))].real
x0

([- 4.9269, 4.7486, 2.9158])

x0 fsolve:

fsolve(func, x0)

([- 4.9848, 4.5462, 2.7192])

. pychebfun , . , ( ) . ( ), ( fsolve).

, fsolve

import pychebfun

f_cheb = pychebfun.Chebfun.from_function(func, domain = (-10,10))
f_cheb.roots()