Find a solution to , using x0 as starting point or interval for x.
Multidimensional overdetermined systems are supported. You can specify them using a function or a list of functions.
If the found root does not satisfy , an exception is raised (this can be disabled with verify=False).
Arguments
solver has to be callable with (f, x0, **kwargs) and return an generator yielding pairs of approximative solution and estimated error (which is expected to be positive). You can use the following string aliases: ‘secant’, ‘mnewton’, ‘halley’, ‘muller’, ‘illinois’, ‘pegasus’, ‘anderson’, ‘ridder’, ‘anewton’, ‘bisect’
See mpmath.calculus.optimization for their documentation.
Examples
The function findroot() locates a root of a given function using the secant method by default. A simple example use of the secant method is to compute as the root of closest to :
>>> from mpmath import *
>>> mp.dps = 30; mp.pretty = True
>>> findroot(sin, 3)
3.14159265358979323846264338328
The secant method can be used to find complex roots of analytic functions, although it must in that case generally be given a nonreal starting value (or else it will never leave the real line):
>>> mp.dps = 15
>>> findroot(lambda x: x**3 + 2*x + 1, j)
(0.226698825758202 + 1.46771150871022j)
A nice application is to compute nontrivial roots of the Riemann zeta function with many digits (good initial values are needed for convergence):
>>> mp.dps = 30
>>> findroot(zeta, 0.5+14j)
(0.5 + 14.1347251417346937904572519836j)
The secant method can also be used as an optimization algorithm, by passing it a derivative of a function. The following example locates the positive minimum of the gamma function:
>>> mp.dps = 20
>>> findroot(lambda x: diff(gamma, x), 1)
1.4616321449683623413
Finally, a useful application is to compute inverse functions, such as the Lambert W function which is the inverse of , given the first term of the solution’s asymptotic expansion as the initial value. In basic cases, this gives identical results to mpmath’s built-in lambertw function:
>>> def lambert(x):
... return findroot(lambda w: w*exp(w) - x, log(1+x))
...
>>> mp.dps = 15
>>> lambert(1); lambertw(1)
0.567143290409784
0.567143290409784
>>> lambert(1000); lambert(1000)
5.2496028524016
5.2496028524016
Multidimensional functions are also supported:
>>> f = [lambda x1, x2: x1**2 + x2,
... lambda x1, x2: 5*x1**2 - 3*x1 + 2*x2 - 3]
>>> findroot(f, (0, 0))
[-0.618033988749895]
[-0.381966011250105]
>>> findroot(f, (10, 10))
[ 1.61803398874989]
[-2.61803398874989]
You can verify this by solving the system manually.
Please note that the following (more general) syntax also works:
>>> def f(x1, x2):
... return x1**2 + x2, 5*x1**2 - 3*x1 + 2*x2 - 3
...
>>> findroot(f, (0, 0))
[-0.618033988749895]
[-0.381966011250105]
Multiple roots
For multiple roots all methods of the Newtonian family (including secant) converge slowly. Consider this example:
>>> f = lambda x: (x - 1)**99
>>> findroot(f, 0.9, verify=False)
0.918073542444929
Even for a very close starting point the secant method converges very slowly. Use verbose=True to illustrate this.
It is possible to modify Newton’s method to make it converge regardless of the root’s multiplicity:
>>> findroot(f, -10, solver='mnewton')
1.0
This variant uses the first and second derivative of the function, which is not very efficient.
Alternatively you can use an experimental Newtonian solver that keeps track of the speed of convergence and accelerates it using Steffensen’s method if necessary:
>>> findroot(f, -10, solver='anewton', verbose=True)
x: -9.88888888888888888889
error: 0.111111111111111111111
converging slowly
x: -9.77890011223344556678
error: 0.10998877665544332211
converging slowly
x: -9.67002233332199662166
error: 0.108877778911448945119
converging slowly
accelerating convergence
x: -9.5622443299551077669
error: 0.107778003366888854764
converging slowly
x: 0.99999999999999999214
error: 10.562244329955107759
x: 1.0
error: 7.8598304758094664213e-18
ZeroDivisionError: canceled with x = 1.0
1.0
Complex roots
For complex roots it’s recommended to use Muller’s method as it converges even for real starting points very fast:
>>> findroot(lambda x: x**4 + x + 1, (0, 1, 2), solver='muller')
(0.727136084491197 + 0.934099289460529j)
Intersection methods
When you need to find a root in a known interval, it’s highly recommended to use an intersection-based solver like 'anderson' or 'ridder'. Usually they converge faster and more reliable. They have however problems with multiple roots and usually need a sign change to find a root:
>>> findroot(lambda x: x**3, (-1, 1), solver='anderson')
0.0
Be careful with symmetric functions:
>>> findroot(lambda x: x**2, (-1, 1), solver='anderson')
Traceback (most recent call last):
...
ZeroDivisionError
It fails even for better starting points, because there is no sign change:
>>> findroot(lambda x: x**2, (-1, .5), solver='anderson')
Traceback (most recent call last):
...
ValueError: Could not find root within given tolerance. (1 > 2.1684e-19)
Try another starting point or tweak arguments.
1d-solver generating pairs of approximative root and error.
Needs starting points x0 and x1 close to the root. x1 defaults to x0 + 0.25.
Pro:
Contra:
1d-solver generating pairs of approximative root and error.
Needs starting points x0 close to the root.
Pro:
Contra:
1d-solver generating pairs of approximative root and error.
Needs starting point x0 close to the root. Uses modified Newton’s method that converges fast regardless of the multiplicity of the root.
Pro:
Contra:
1d-solver generating pairs of approximative root and error.
Needs a starting point x0 close to the root. Uses Halley’s method with cubic convergence rate.
Pro:
Contra:
1d-solver generating pairs of approximative root and error.
Needs starting points x0, x1 and x2 close to the root. x1 defaults to x0 + 0.25; x2 to x1 + 0.25. Uses Muller’s method that converges towards complex roots.
Pro:
Contra:
1d-solver generating pairs of approximative root and error.
Uses bisection method to find a root of f in [a, b]. Might fail for multiple roots (needs sign change).
Pro:
Contra:
1d-solver generating pairs of approximative root and error.
Uses Illinois method or similar to find a root of f in [a, b]. Might fail for multiple roots (needs sign change). Combines bisect with secant (improved regula falsi).
The only difference between the methods is the scaling factor m, which is used to ensure convergence (you can choose one using the ‘method’ keyword):
Pro:
Contra:
1d-solver generating pairs of approximative root and error.
Uses Pegasus method to find a root of f in [a, b]. Wrapper for illinois to use method=’pegasus’.
1d-solver generating pairs of approximative root and error.
Uses Anderson-Bjoerk method to find a root of f in [a, b]. Wrapper for illinois to use method=’pegasus’.
1d-solver generating pairs of approximative root and error.
Ridders’ method to find a root of f in [a, b]. Is told to perform as well as Brent’s method while being simpler.
Pro:
Contra:
EXPERIMENTAL 1d-solver generating pairs of approximative root and error.
Uses Newton’s method modified to use Steffensens method when convergence is slow. (I.e. for multiple roots.)
Find the root of a vector function numerically using Newton’s method.
f is a vector function representing a nonlinear equation system.
x0 is the starting point close to the root.
J is a function returning the Jacobian matrix for a point.
Supports overdetermined systems.
Use the ‘norm’ keyword to specify which norm to use. Defaults to max-norm. The function to calculate the Jacobian matrix can be given using the keyword ‘J’. Otherwise it will be calculated numerically.
Please note that this method converges only locally. Especially for high- dimensional systems it is not trivial to find a good starting point being close enough to the root.
It is recommended to use a faster, low-precision solver from SciPy [1] or OpenOpt [2] to get an initial guess. Afterwards you can use this method for root-polishing to any precision.
[1] http://scipy.org