Roots of equations

Solve Linear/Nonlinear Equations

Please also refer to publication.


@_dc.dataclass
class Info:
   """
   err: error (if available)
   iter: number of iterations to reach the root
   conv: whether converged to a root or not
   msg: if convergence is False, a reason is given
   """
   err:float = None
   iter:int = -1
   conv:bool = False
   msg:str =""

bisect

Finds the root of an equation of form f(x)=0.


def bisect(
   f:FunctionType, 
   a:Real, 
   b:Real, 
   tol=1E-5, 
   maxiter=100, 
   method="bf", 
   modified=False)->tuple[float, Info]:
"""
f: A unary function
a, b:	The interval where the root lies in
tol:	tolerance for error 
maxiter: Maximum number of iterations 
method: "bf" for brute-force (halving) 
"rf" for regula falsi (false position) 
modified: True for modified regula falsi method.
"""

Example #1


from scisuit.roots import bisect

f = lambda x: x**2-5
result = bisect(f=f, a = 0, b = 5 )
print(result)

Bisection using brute-force method 
Root=2.23607, Error=9.54e-06 (18 iterations).

Example #2


from scisuit.roots import bisect

result = bisect(lambda x: x**2-5, a = 0, b = 4, method = ("rf", False))
print(result)

Bisection using regula falsi method 
Root=2.23607, Error=2.97e-06 (12 iterations).

itp

Finds the root of an equation of form f(x)=0 using ITP (interpolation, truncation, projection) method.


def itp(
   f:FunctionType, 
   a:Real, 
   b:Real, 
   k1:Real = 0.1,
   k2:Real = 2.5656733089749,
   tol:Real=1E-5, 
   maxiter:int=100)->tuple[float, Info]:
"""
f: A unary function
a, b: The interval where the root lies in
k1, k2: parameters to control interpolation phase
tol: tolerance for error
maxiter: Maximum number of iterations

Reference:
- Oliveira IFD & Takahashi RHC (2020). An Enhancement of the Bisection Method Average 
  Performance Preserving Minmax Optimality. ACM Transactions on Mathematical Software, 47:1
"""

Example


from scisuit.roots import bisect, itp

args = {"f":lambda x: x**2-5, "a":0, "b": 4} 

for func in [bisect, itp]:
    print(func(**args))

Bisection using brute-force method 
Root=2.23607, Error=3.35e-06 (18 iterations).

ITP method
Root=2.23607, Error=3.43e-07 (7 iterations).

brentq

Finds the root of an equation of form f(x)=0. Uses the Brent's method (1973) to find the root of the function using inverse quadratic interpolation.


def brentq(
   f:FunctionType, 
   a:Real, 
   b:Real, 
   tol=1E-5, 
   maxiter=100)->tuple[float, Info]:
"""
f: A unary function whose root is sought after \n
a, b: The interval where root lies in \n
tol: tolerance for error \n
maxiter: Maximum number of iterations during the search for the root
"""

Example


from scisuit.roots import bisect, brentq

args = {"f":lambda x: x**2-5, "a":0, "b": 4} 

for func in [bisect, brentq]:
    print(func(**args))

Bisection using brute-force method 
Root=2.23607, Error=3.35e-06 (18 iterations).

Brent's method (inverse quadratic interpolation)
Root=2.23607, (7 iterations).

fsolve

Solves system of non-linear equations using Newton-Raphson approach.


def fsolve(
   F:list[FunctionType], 
   x0:list[float], 
   tol=1E-5, 
   maxiter=100 )->tuple:
   
"""
F: a list of functions
x0: a list of initial guesses
"""

Example


"""
x^2 + y^2 = 5
x^2 - y^2 = 1

First define the functions, F(x, y) = 0:
"""

def f1(t):
    return t[0]**2 + t[1]**2 - 5

def f2(t):
    return t[0]**2 - t[1]**2 - 1

roots, iter=fsolve( [f1,f2], [1,1] )

print("Root=", roots) 
print("iter:", iter) 

print("f1(roots)=",f1(roots))
print("f2(roots)=", f2(roots))

Root=1.73205     1.41421     
iter:5 

f1(roots)=9.428e-09    
f2(roots)=9.377e-09

muller

Finds the root of an equation of form f(x)=0.


def muller(
   f:FunctionType, 
   x0:Complex, 
   h=None, 
   x1=None, 
   x2=None, 
   tol=1E-5, 
   maxiter=100)->tuple[float, Info]:
"""
f: A unary function
x0, x1, x2:	Initial guesses
h: Step length
tol: tolerance for error
maxiter: Max number of iterations
"""

Notes on usage and implementation:

The function must be comprised of elements all of which can handle complex numbers.
If only x0 is provided then h = 0.5.
If h is provided, x1 and x2 can not be used. Conversely, If x1 and x2 are provided, h can not be used.
If x0 and/or h are provided, then:
x₁ = x₀ + x₀·h
x₂ = x₀ - x₀·h

newton

Finds the root of an equation of form f(x)=0.


def newton(
   f:FunctionType, 
   x0:Real, 
   x1=None, 
   fprime=None, 
   fprime2=None,
   tol=1E-5, 
   maxiter=100)->tuple[float, Info]:
"""
- fprime != None, Newton-Raphson is used,
- fprime2 != None, Halley's method is used,
- fprime == None, x1 must be provided and Secant method is used.

f: A unary function 
fprime: derivative of f 
x0, x1: Initial guesses 
tol: tolerance for error 
maxiter: Max number of iterations
"""

ridder

Finds the root of an equation of form f(x)=0.


def ridder(
   f:types.FunctionType, 
   a:Real, 
   b:Real, 
   tol=1E-5, 
   maxiter=100)->tuple[float, Info]:
"""
f: A unary function 
a, b: The interval where the root lies in 
tol: tolerance for error 
maxiter: Maximum number of iterations
"""

TOMS 748

Finds the root of an equation of form f(x)=0 using TOMS Algorithm 748 by Alefeld, Potra and Shi (1995) (ACMTransactions on Mathematica1 Software, Vol. 21. No. 3. September 1995. Page. 327-344)


def toms748(
   f:FunctionType, 
   a:Real, 
   b:Real, 
   tol=1E-5, 
   maxiter=100)->tuple[float, Info]:
"""
Algorithm TOMS 748: Alefeld, Potra and Shi: Enclosing zeros of continuous functions

f: A unary function whose root is sought after
a, b: The interval where root lies in,
tol: tolerance for error (weakly valid)
maxiter: Maximum number of iterations during the search for the root

## Reference:
https://beta.boost.org/doc/libs/1_82_0/libs/math/doc/html/math_toolkit/roots_noderiv/TOMS748.html
"""


from scisuit.roots import toms748

def func(x):
    return x**2-5

result = toms748(f=func, a=0, b=5, tol=1E-6)

print(result)
"""
(2.236067978, Info(err=1.91e-13, iter=-1, conv=True, msg=''))
"""