CODEWITHSUNDEEP
pythonscipymathematical-optimization
heracho
heracho

Reputation: 610

need to improve accuracy in fsolve to find multiples roots

I'm using this code to get the zeros of a nonlinear function. Most certainly, 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 give the 3 expected roots without any problem. But, with c = -1.5 the code miss a root, and with c = -3 it find a non existing root.

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

I appreciate any solution, trick or advice.

Upvotes: 3

Views: 2007

Answers (2)

Dietrich
Dietrich

Reputation: 5531

Between two stationary points (i.e., df/dx=0), you have one or zero roots. In your case it is possible to calculate the two stationary points analytically:

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

So you have three intervals where you need to find a zero. Using Sympy saves you from doing the calculations by hand. Its sy.nsolve() allows to robustly find a zero in an interval:

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

Upvotes: 2

Stelios
Stelios

Reputation: 5521

What you need is an automatic way to obtain good initial estimates of the roots of the function. This is in general a difficult task, however, for univariate, continuous functions, it is rather simple. The idea is to note that (a) this class of functions can be approximated to an arbitrary precision by a polynomial of appropriately large order, and (b) there are efficient algorithms for finding (all) the roots of a polynomial. Fortunately, Numpy provides functions for both performing polynomial approximation and finding polynomial roots.

Let's consider a specific function

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

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

The following code uses polyfit and poly1d to approximate func over the range of interest (-10<x<10) by a polynomial function f_poly of order 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)

As the following plot shows, f_poly is indeed a good approximation of func. Even greater accuracy can be obtained by increasing the order. However, there is no point in pursuing extreme accuracy in the polynomial approximation, since we are looking for approximate estimates of the roots that will be later refined by fsolve

enter image description here

The roots of the polynomial approximation can be simply obtained as 

roots = np.roots(pfit)
roots

array([-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 ])

As expected, Numpy returns 10 complex roots. However, we are only interested for real roots within the interval [-10,10]. These can be extracted as follows:

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

array([-4.9269, 4.7486, 2.9158])

Array x0 can serve as the initialization for fsolve:

fsolve(func, x0)

array([-4.9848, 4.5462, 2.7192])

Remark: The pychebfun package provides a function that directly gives all the roots of a function within an interval. It is also based on the idea of performing polynomial approximation, however, it uses a more sophisticated (yet, more efficient) approach. It automatically chooses the best polynomial order of the approximation (no user input), with the polynomial roots being practically equal to the true ones (no need to refine them via fsolve).

This simple code gives the same roots as those by fsolve

import pychebfun

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

Upvotes: 3

Related Questions