Jacques Gaudin - 2 years ago 162
Python Question

# Scipy root-finding method

I am writing a class that has an mathematical function as an attribute, say f.

f is:

• Defined on a real segment [-w;+w]

• Positive and bounded above by a real H

• even (for all x in [-w;+w], f(x)=f(-x)) and f(w)=f(-w)=0

• Differentiable over [-w;+w] and its derivative is positive and continuous over [-w;0]

My class looks like :

from scipy.misc import derivative
from math import cosh, sqrt

class Function(object):

w = 1.
PRECISION = 1e-6

def f(self, x):
'''This is an example but f could be
any math function matching requirments above.
'''
return 0.5+1.07432*(1-cosh(x/1.07432))

def deriv_f(self, x):
return derivative(self.f, x, self.PRECISION)

def x_to_arc_length(self, x):
def func(x):
return sqrt(1+self.deriv_f(x)**2)

def arc_length_to_x(self, L):
bound = [-self.w, self.w]
while bound[1]-bound[0] > self.PRECISION:
mid= sum(bound)/2
bound[(self.x_to_arc_length(mid)-L > 0)] = mid
return sum(bound)/2

I use bisection to inverse the arc length method, but I was looking at changing this for one of the
scipy.optimize
root-finding method to gain speed.
I am new to scipy and must admit that my math are a bit rusted...
Scipy gives me the choice between
brentq
,
brenh
,
ridder
,
bisect
and
newton
.

Could anyone point me to the best-suited method for a case like this ? Or maybe there is a better library for this ?