# Root finding strategies in Numerical Analysis

Charbel Tannous

► To cite this version:

Charbel Tannous. Root finding strategies in Numerical Analysis. Master. Analyse Numérique, UBO Brest, France. 2017, pp.22. hal-04691429

HAL Id: hal-04691429

https://hal.science/hal-04691429v1

Submitted on 8 Sep 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Root finding strategies in Numerical Analysis

C. Tannous

Université de Brest, Lab-STICC, CNRS-UMR 6285, F-29200 Brest, FRANCE

(Dated: September 8, 2024)

Algorithms for finding roots of non-linear functions in 1D and higher dimensions are described and analyzed from the accuracy, speed and application points of view.

PACS numbers: 84.30.Vn, 87.55.de, 02.60.Ed, *43.28.Js, 43.25.+y, 07.07.Tw

Keywords: Electronic Filters, optimization techniques, Interpolation methods, numerical methods, nonlinear systems, Control devices

## Contents

I. Introduction 1

II. General 1D methods 2

A. Bisection method 2

B. Regula-falsi method 3

C. Secant method 3

D. Newton-Raphson method 4

1. Variants of Newton-Raphson 6

E. Fixed point iteration 7

1. Oscillation and Chaos 8

III. Polynomial roots 8

A. Root counting 8

1. Sturm chain 8

2. Routh array 10

B. Real roots 12

C. Complex roots 13

IV. Extremum Problems in Higher dimensions 14

A. Powell Method and minimization 14

B. Broyden method and fitting 15

V. Applications 16

A. Evaluation of mathematical constants 16

B. Arithmetic operations and evaluation of functions 17

C. Signal reconstruction from its zeroes 18

D. Filter design from zeroes and poles 19

E. Stability of control and feedback systems 20

References 21

## I. INTRODUCTION

Finding roots of non-linear functions or systems of non-linear functions is a very important problem in Science, Engineering, Optimization, Social Sciences...

1. Crossing point of two functions in ballistics, missile guiding, economics... can be simply described mathematically by the equation satisfied by $x^*$ such that:

$$f_1(x^*) = f_2(x^*) \rightarrow f(x^*) = f_1(x^*) - f_2(x^*) \equiv 0 \quad (1)$$
2. Level crossing point $x^*$ of a function in hydraulics, stock market share value... can be simply viewed as a solution to:

$$f_1(x^*) = A \rightarrow f(x^*) = f_1(x^*) - A \equiv 0 \qquad (2)$$

3. Extremum problems in Optimization, maximization, minimization... can be described by:

$$f'_1(x^*) = 0 \rightarrow f(x^*) = f'_1(x^*) \equiv 0 \qquad (3)$$

4. Eigenvalue problems in $n \times n$ matrix $A$ diagonalization ... can be described by:

$$\det(A - \lambda \mathbf{1}) = 0 \rightarrow f(x^*) = \det(A - x^* \mathbf{1}) \equiv 0 \qquad (4)$$

5. Signal construction from its zeroes [1, 2]

6. Filter design from its zeroes and poles

7. Inverse problems [3]: Given $F(x) = a$, find $x = F^{-1}(a)$ by solving $f(x) = F(x) - a = 0$

8. Fixed-point problems: $x_{n+1} = F(x_n)$ then solve $f(x) = F(x) - x = 0$

## II. GENERAL 1D METHODS

<figure><img src="image_2.png" /><figcaption>Fig.1: (Color on-line) Graph of a typical function $f(x)$ displaying a single zero $x^*$ over a bracketing interval $[a, b]$ with $f(a) \cdot f(b) \le 0$. In the bisection method, the root is approximated by $(a + b)/2$ at mid-interval whereas in the Regula-Falsi method, $f(x)$ is replaced by the chord (in blue) intersecting the $x$ axis at the approximate root $x_c$.</figcaption></figure>

### A. Bisection method

Given a function $f(x)$ (cf. Fig. 1) displaying a single zero in a bracketing interval $[a, b]$ with condition $f(a) \cdot f(b) \le 0$, the root is estimated at mid-interval $(a + b)/2$. After first root is approximated, bisection is performed repeatedly to narrow the bracketing interval and estimate root at mid interval following the algorithm:

1. Given $f(x)$ continuous over $[a, b]$ with $f(a) \cdot f(b) \le 0$
2. Bisect the interval to estimate root and redefine the root bracketing interval

3. At step $n + 1$, root $x_{n+1}$ is always bracketed by interval $[a_n, b_n]$ and estimated to be $x_{n+1} = (a_n + b_n)/2$ with $f(a_n) \cdot f(b_n) \le 0$

4. Estimated error: $\epsilon_n = \max(b_n - a_n)$

Thus we infer that $\epsilon_{n+1} = \epsilon_n/2$ whose solution is $\epsilon_n = \epsilon_0/2^n$. This means that starting with $\epsilon_0$, then 20 iterations are needed to reduce the error to $\epsilon_{20} \approx 10^{-6}\epsilon_0$.

Example functions given below:

1. $f(x) = -0.9x^2 + 1.7x + 2.5$ in interval $[2.8, 3.0]$

2. $f(x) = \frac{1-0.61x}{x}$ in interval $[1.5, 2.0]$

3. $f(x) = x^2 |\sin x| - 4.1$ in interval $[0, 4]$

show that after three iterations, in the case of $f(x) = -0.9x^2 + 1.7x + 2.5$, the root is $x^* = 2.8601$, whereas for $f(x) = \frac{1-0.61x}{x}$, the root is $x^* = 1.6393$. In the case of $f(x) = x^2 |\sin x| - 4.1$, $x^* = 3.485775$ with respective errors given by $\epsilon_3 = \epsilon_0/8$ where $\epsilon_0$ is the initial error given by $|x_0 - x^*|$.

## B. Regula-falsi method

This method which is a variant of bisection, consists of replacing the mid-interval $[x_1, x_2]$ point $x_m$ by the intersection $x_m^*$ of the line (called chord) joining points $(x_1, f(x_1))$ and $(x_2, f(x_2))$ and the x axis, thus the conditions:

1. Given $f(x)$ continuous over $[a, b]$ with $f(a) \cdot f(b) \le 0$

2. Potential root should always be bracketed

3. Estimated error: $\epsilon_{n+1} \sim \epsilon_n$

In Fig. 1, it is shown how $f(x)$ is replaced by its chord to get an approximate root as $x_c$.

Applying this approach to $f(x) = -0.9x^2 + 1.7x + 2.5$ no convergence is obtained, whereas for $f(x) = \frac{1-0.61x}{x}$, the root is $x^* = 1.6393$ as in bisection. No convergence is obtained for $f(x) = x^2 |\sin x| - 4.1$ as well.

## C. Secant method

The secant method is a root-finding algorithm that uses a succession of approximate roots obtained with chord intersections (as in the regula-falsi) but with a different strategy. The conditions to apply secant method are:

1. Given $f(x)$ continuous over $[a, b]$ with $f(a) \cdot f(b) \le 0$

2. $f'(x^*) \ne 0$ and $f''(x)$ non-singular

3. $x^*$ bracketed only initially

4. Start point $x_0$ should be close to $x^*$

In order to find a zero of a function $f$, the secant method is defined by a recursion relation over an interval $[a_n, b_n]$, exploiting chord intersection with x axis as $x_{n+1}$ given by:

$$x_{n+1} = a_n - \frac{f(a_n)(b_n - a_n)}{f(b_n) - f(a_n)} \quad (5)$$

A major paradigm shift is then used which is to relax the bracketing of the root during the iterations. Thus the root is bracketed only at start, to initiate the iteration mechanism and no longer during following iterations.
<figure><img src="image_3.png" /><figcaption>**Fig.2:** (Color on-line) Secant method illustrated by a recursion $x_1, x_2 \to x_3$ where $x_1, x_2$ are bounds not necessarily bracketing the real zero but used to derive from the chord replacing the function $f(x)$, a point $x_3$, the chord intersection point with x axis. By inspection, we infer $x_3 = x_1 - \frac{f(x_1)}{f(x_2)-f(x_1)}(x_2 - x_1)$.</figcaption></figure>

This implies that given an interval limited by two points $x_1, x_2$ the intersection of the chord replacing the function is $x_3$ such that

$$x_3 = x_1 - \frac{f(x_1)}{f(x_2) - f(x_1)}(x_2 - x_1) \quad (6)$$

Considering this relationship as generalizable to the $(n+1)$-st iterate depending on $x_n, x_{n-1}$, we get the correspondence $x_1 \to x_{n-1}, x_2 \to x_n, x_3 \to x_{n+1}$. This can be written as:

$$x_{n+1} = x_{n-1} - \frac{f(x_{n-1})}{f(x_n) - f(x_{n-1})}(x_n - x_{n-1}) \quad (7)$$

Thus $x_{n+1} = F(x_n, x_{n-1})$ allows to Taylor expand $F$ about $x_n \sim x^*, x_{n-1} \sim x^*$ such that:

$$x_{n+1} = F(x^*, x^*) + \epsilon_n \left(\frac{\partial F}{\partial x_n}\right)_{(x^*, x^*)} + \epsilon_{n-1} \left(\frac{\partial F}{\partial x_{n-1}}\right)_{(x^*, x^*)} \\ + \frac{1}{2} \epsilon_n^2 \left(\frac{\partial^2 F}{\partial x_n^2}\right)_{(x^*, x^*)} + \frac{1}{2} \epsilon_{n-1}^2 \left(\frac{\partial^2 F}{\partial x_{n-1}^2}\right)_{(x^*, x^*)} + \epsilon_n \epsilon_{n-1} \left(\frac{\partial^2 F}{\partial x_n \partial x_{n-1}}\right)_{(x^*, x^*)} \dots \quad (8)$$

Writing $x_{n+1} = x^* + \epsilon_{n+1}$ and evaluating derivatives, we get:

$$\epsilon_{n+1} = \epsilon_n \epsilon_{n-1} \left( \frac{\partial^2 F}{\partial x_n \partial x_{n-1}} \right)_{(x^*, x^*)} \quad (9)$$

The recursion relation $\epsilon_{n+1} \sim \epsilon_n \epsilon_{n-1}$ is solved by taking the logarithm $u_n = \ln \epsilon_n$ that yields $u_{n+1} \sim u_n + u_{n-1}$, the Fibonacci recurrence. The latter is solved by taking $u_n \sim r^n$ to get the characteristic equation: $r^2 = r + 1$ whose solution is the Golden mean $\phi = (1 + \sqrt{5})/2$.

Thus the error evolves as $\epsilon_{n+1} \sim \epsilon_n^\phi$ which is a superlinear convergence since $\phi \approx 1.6$

## D. Newton-Raphson method

Newton-Raphson method is based on the function $f(x)$ and its derivative $f'(x)$ since the algorithm is given by:
<figure><img src="image_4.png" /><figcaption>Fig.3: (Color on-line) Newton-Raphson algorithm is to find the intersection $x_{n+1}$ with the x axis of the line (in blue) tangent to $f(x)$ at $x_n$.</figcaption></figure>

1. Start from an interval $[a, b]$ bracketing the root $x^*$

2. Fig. 3 shows that the line tangent to $f(x)$ at $x_n$ intersects the x axis such that $x_{n+1} = x_n - f(x_n)/f'(x_n)$

Let $f(x) = x^3 - x - 3$ the goal is to find $x^*$ such that $f(x^*) = 0, x^* \in \mathbb{R}$.

Newton-Raphson method converges starting with $x_0 = 1$ and moving toward $x^* = 1.67169988165734$ in double precision after 6 iterations. It diverges when started at $x_0 = 0$.

<figure><img src="image_5.png" /><figcaption>Fig.4: Graph of $f(x) = x^3 - x - 3$ displaying a single zero.</figcaption></figure>

In order to estimate convergence speed, let $g(x) = x - f(x)/f'(x)$, then Newton-Raphson iteration is $x_{n+1} = g(x_n)$ with a fixed point $x^* = g(x^*)$.

Writing $x^* + \epsilon_{n+1} = g(x^* + \epsilon_n)$ and expanding we get:

$$ \epsilon_{n+1} = \epsilon_n \left( \frac{dg}{dx} \right)_{x^*} + \frac{1}{2} \epsilon_n^2 \left( \frac{dg}{dx} \right)_{x^*} \dots \qquad (10) $$

Derivative evaluation yields: $\left(\frac{dg}{dx}\right) = \frac{f(x)f''(x)}{[f'(x)]^2}$ resulting in $\epsilon_{n+1} = \frac{1}{2}\epsilon_n^2 \left(\frac{dg}{dx}\right)_{x^*}$ or $\epsilon_{n+1} \sim \epsilon_n^2$ i.e. superlinear quadratic convergence. In order to realize convergence speed, if one starts with $\epsilon_0 = 0.1$ after one iteration $\epsilon_1 = 10^{-2}$, then $\epsilon_2 = 10^{-4}$ and after three iterations $\epsilon_3 = 10^{-8}$.
# 1. Variants of Newton-Raphson

1. We consider three variants of Newton-Raphson method: with division, without division and chord iteration.

Let $a$ be a real positive number. Consider the following algorithm:

Starting at $x_0$, iterate $x_{n+1} = x_n + \frac{1}{2}(a - x_n^2)$, for $n = 0, 1, 2, ...$

Results of Newton-Raphson method:

With division: $x_{n+1} = x_n + \frac{1}{2}(a - x_n^2)$ converges to $\sqrt{a}$ after 3 iterations:

<table><thead><tr><th>Iterations</th><th>x<sup>*</sup> value</th></tr></thead><tbody><tr><td>0</td><td>1.700000000000000</td></tr><tr><td>1</td><td>1.73235294117647</td></tr><tr><td>2</td><td>1.73205083391591</td></tr><tr><td>3</td><td>1.73205080756888</td></tr></tbody></table>

Better method without division: $x_{n+1} = \frac{x_n}{2}(3 - ax_n^2)$ (it converges to $1/\sqrt{a}$ after 4 iterations and requires a multiplication by $a$ at the end of iterations.)

<table><thead><tr><th>Iterations</th><th>x<sup>*</sup> value</th></tr></thead><tbody><tr><td>0</td><td>0.588235294117647</td></tr><tr><td>1</td><td>1.73112151434968</td></tr><tr><td>2</td><td>1.73205005981531</td></tr><tr><td>3</td><td>1.73205080756839</td></tr><tr><td>4</td><td>1.73205080756888</td></tr></tbody></table>

Validation results:

1. Test with division: $(x^*)^2 - a \sim 10^{-16}$

2. Test without division: $(x^*)^2 - a \sim 10^{-16}$

2. Double iteration:

Another variety of Newton-Raphson method is to consider two mathematically equivalent expressions and study their convergence properties.

An example is given by a real number $\lambda > 0$ and requiring to find $\bar{x}$ such that $\bar{x} = \lambda e^{\bar{x}}$.

Fig. 5 shows that when $\lambda < 1/e$, there are two positive solutions $\bar{x}_1$ et $\bar{x}_2$, whereas for $\lambda > 1/e$, there is no solution and when $\lambda = 1/e$, there is one solution.

Thus we perform two types of iterations:

$$ (a) \text{ Type 1: } x_{n+1} = \lambda e^{x_n}, \quad n = 0, 1, 2, 3, \dots $$

$$ (b) \text{ Type 2: } x_{n+1} = \ln x_n - \ln \lambda, \quad n = 0, 1, 2, 3, \dots $$

<figure><img src="image_6.png" /><figcaption>Fig.5: Graphs showing that for $\lambda = 1/e$ there are two solutions whereas for $\lambda < 1/e$, only one solution is possible.</figcaption></figure>
As displayed in Fig. 5, type 1 method converges toward $\bar{x}_1$ when $x_0$ is close to $\bar{x}_1$. Type 2 method converges toward $\bar{x}_2$ when $x_0$ is close to $\bar{x}_2$.

### 3. Chord Newton-Raphson iteration:

Replacing the iteration $x_{n+1} = x_n - f(x_n)/f'(x_n)$ by $x_{n+1} = x_n - f(x_n)/f'(x_0)$ makes Newton-Raphson a lot faster since the derivative $f'(x_0)$ is calculated only once. Example below evaluates $\sqrt{2} = 1.414213562373$ in double precision with chord and standard Newton-Raphson with a start point $x_0 = 1$:

Chord Newton-Raphson results:

<table><thead><tr><td>Iterations</td><td>x* value</td><td>Error</td></tr></thead><tbody><tr><td>1</td><td>1.50000000000000</td><td>0.08578643762690</td></tr><tr><td>2</td><td>1.37500000000000</td><td>0.03921356237309</td></tr><tr><td>3</td><td>1.42968750000000</td><td>0.01547393762690</td></tr><tr><td>4</td><td>1.4076843261718</td><td>0.00652923620122</td></tr><tr><td>5</td><td>1.4168967450968</td><td>0.00268318272379</td></tr></tbody></table>

Standard Newton-Raphson results:

<table><thead><tr><td>Iterations</td><td>x* value</td><td>Error</td></tr></thead><tbody><tr><td>1</td><td>1.50000000000000</td><td>0.0857864376269049</td></tr><tr><td>2</td><td>1.41666666666666</td><td>0.0024531042935716</td></tr><tr><td>3</td><td>1.41421568627450</td><td>0.0000021239014148</td></tr><tr><td>4</td><td>1.41421356237468</td><td>1.5948210634053E-12</td></tr><tr><td>5</td><td>1.41421356237309</td><td>9.6710833785707E-17</td></tr></tbody></table>

Convergence of chord Newton-Raphson is slower than standard Newton-Raphson but execution time is much smaller since the derivative is not calculated at every step.

In order to estimate convergence speed of chord Newton-Raphson, Writing $x_{n+1} = x^* + \epsilon_{n+1}$, $x_n = x^* + \epsilon_n$, substituting in $x_{n+1} = x_n - f(x_n)/f'(x_0)$ and expanding we get:

$$ \epsilon_{n+1} = \epsilon_n \left[ 1 - \frac{1}{f'(x_0)} \left( \frac{df}{dx} \right)_{x^*} \right] - \frac{1}{2} \frac{\epsilon_n^2}{f'(x_0)} \left( \frac{d^2 f}{dx^2} \right)_{x^*} \dots \quad (11) $$

This means that in general $\epsilon_{n+1} \sim \epsilon_n$ i.e. linear convergence unless one selects the start point $x_0$ such that $f'(x_0) = \left(\frac{df}{dx}\right)_{x^*}$ yielding quadratic convergence $\epsilon_{n+1} \sim \epsilon_n^2$ like standard Newton-Raphson.

However, when comparing with secant method where $\epsilon_{n+1} \sim \epsilon_n^\phi$ and no derivative evaluation, a fair comparison should be such that one Newton-Raphson step is equivalent to two secant steps.

This implies that Newton-Raphson $\epsilon_{n+1}$ is equivalent to secant $\epsilon_{n+2} = \epsilon_{n+1}^\phi = [\epsilon_n^\phi]^\phi = \epsilon_n^\phi = \epsilon_n^{\phi+1} \approx \epsilon_n^{2.6}$. This means, the secant method achieves super-quadratic convergence when true function evaluations are accounted for.

## E. Fixed point iteration

Consider root finding of $f_1(x^*) = 0$. Reformulating the problem as $f_1(x^*)+x^* = x^*$, then writing $f(x) = f_1(x)+x$, yields: $f(x^*) = x^*$, a fixed point problem in the sense it can be solved iteratively as: $x_{n+1} = f(x_n)$ such that $\lim_{n\to\infty} x_n = x^*$.

Fixed point iteration has the following conditions tied to its performance:

1. Given $f(x)$ continuous with $|1 + f'(x)| < 1$ to warrant a unique root

2. $f'(x)$ smooth and non-singular

3. $x^*$ bracketed only initially

For illustration, let us analyze of the following functions, along with respective stability of their fixed points:

1. $f(x) = 4x - x^2$

2. $f(x) = \sqrt{x}$
3. $f(x) = \sin^{-1} x$

Starting from the iteration $x_{n+1} = f(x_n)$ and expanding about the fixed point $x^* = f(x^*)$, we get:

$$x^* + \epsilon_{n+1} = f(x^* + \epsilon_n) \approx f(x^*) + \epsilon_n \left(\frac{df}{dx}\right)_{x^*} \qquad (12)$$

Thus the error $\epsilon_{n+1} \sim \epsilon_n$ induces linear convergence. We infer that the fixed point is attractive (stable) when $|a| \le 1$ and repulsive (unstable) when $|a| > 1$. Thus, fixed point $x^* = 3$ of $f(x) = 4x - x^2$ is unstable since $a = -2$, whereas $x^* = 1$ of $f(x) = \sqrt{x}$ is stable since $a = \frac{1}{2}$. Same applies to fixed point $x^* = 0$ of $f(x) = \sin^{-1} x$ which is stable since $a = 1$.

## 1. Oscillation and Chaos

An important drawback of fixed point iteration, is a possibility of getting into an oscillatory or chaotic regime that diverts away from the root being sought. A straightforward criterion to tell whether a recursion function $f(x)$ induces chaos or not is given by the Schwarzian derivative [4] sign:

$$Sf(x) = \frac{f'''(x)}{f'(x)} - \frac{3}{2} \left( \frac{f''(x)}{f'(x)} \right)^2, \quad f'(x) = \frac{df}{dx}, f''(x) = \frac{d^2 f}{dx^2}, f'''(x) = \frac{d^3 f}{dx^3} \qquad (13)$$

The condition is $Sf(x) < 0 \quad \forall x \neq x^*$.

For instance, in the logistic recursion $x_{n+1} = ax_n(1-x_n)$ case defined over the interval $[0, 1]$, $Sf(x)$ negative over the interval $[0, 1]$.

The fixed point $x^*$ obtained in the limit $n \to \infty$, is $x^* = ax^*(1-x^*)$ i.e. $x^* = 1 - 1/a$, assuming $x^* \neq 0$.

$x_{n+1} = ax_n(1-x_n)$ recursion possesses subtle behavior when $a$ is changed over a small interval $[3,4]$ as displayed in Fig.6, thus it is important from the numerical point of view to avoid being trapped in an oscillatory or chaotic behavior and miss the root being sought.

# III. POLYNOMIAL ROOTS

Polynomials deserve a separate treatment since they require a number of particular operations, such as deflation that consists of dividing a polynomial $p(x)$, to reduce its order, by $(x-r)$ where $r$ is a root such that $q(x) = p(x)/(x-r)$. Moreover there are rules for root counting as described below.

## A. Root counting

### 1. Sturm chain

Consider a polynomial of degree $n$ with real coefficients:

$$p(x) = a_n x^n + a_{n-1} x^{n-1} + \cdots + a_1 x + a_0, \quad a_n \neq 0 \qquad (14)$$

The fundamental theorem of algebra states that $p(x)$ has $n$ real or complex roots, counting multiplicities. If the coefficients $a_0, a_1, \cdots, a_n$ are all real, then the complex roots occur in conjugate pairs, that is, in the form $c \pm di$, where $c, d \in \mathbb{R}$ and $i$ the unit complex number.

Using Descartes rules of sign [1], it is possible to count the number of real positive zeros of $p(x)$ a polynomial with real coefficients. More specifically, let $v$ be the number of sign change of $a_n, a_{n-1}, \cdots, a_0$ coefficients (ignoring zero coefficients). Let $n_p$ be the number of real positive zeros. Then we have two rules:

(i) $n_p \le v$,

(ii) $v - n_p$ is an even integer.

A negative zero of $p(x)$, if it exists, is a positive zero of $p(-x)$. The number of real negative zeros of $p(x)$ is related to the number of sign changes in $p(-x)$ coefficients.
<figure><img src="image_7.png" /><figcaption>Fig.6: (Color on-line) Logistic recursion $x_{n+1} = ax_n(1-x_n)$ displaying fixed point, oscillatory and chaotic behavior when $a= 2.9$, 3.1 and 3.92 respectively. Diagram $[x_{n+1}, x_n]$ is shown at left whereas sequence $x_n$ versus $n$ is shown at right. All iterations start from 0.1.</figcaption></figure>

Example:

Consider the polynomial $p(x) = x^4 + 2x^2 - x - 1$. Then $v = 1$, so $n_p$ is either 0 or 1 by rule (i). In addition, by rule (ii) $v - n_p$ must be even. Hence $n_p = 1$.

When we consider $p(-x) = x^4 + 2x^2 + x - 1$, odd coefficients have a sign change, so $p(-x)$ has one positive zero. In other words, $p(x)$ has one negative zero.

To summarize, simply by looking at $p(x)$ coefficients, we conclude that $p(x)$ has one positive real root, one negative real root, and two complex roots as a conjugate pair.

Descartes rule of sign still leaves an uncertainty as to the exact number of real zeros. A rigorous test was given by Sturm, who showed how to count the real roots within an arbitrary range.

Let $f(x)$ be a real polynomial. Denote it by $f_0(x)$ and its derivative $f'(x)$ by $f_1(x)$. Perform polynomial division successively according to Euclid algorithm [5] such as: $f_0(x) = q_1(x)f_1(x) + r(x)$ where $q_1(x)$ is the quotient of $f_0(x)/f_1(x)$ and $r(x)$ is the remainder of the division. Then replace $r(x)$ by $-f_2(x)$ to induce a recurrence mechanism. Continue this procedure until a constant is obtained.

Thus the recursion algorithm:

$$ 
\begin{aligned}
f_0(x) &= q_1(x)f_1(x) - f_2(x), \\
f_1(x) &= q_2(x)f_2(x) - f_3(x), \\
\vdots \\
f_{k-2}(x) &= q_{k-1}(x)f_{k-1}(x) - f_k(x), \\
f_{k-1}(x) &= q_k(x)f_k(x),
\end{aligned} 
\quad (15) $$

where $deg(f_i(x)) < deg(f_{i-1}(x))$, for $1 \le i \le k$.
Note that the divisor $f_k$ that yields zero remainder is a GCD (greatest common divisor) of $f(x)$ and $f'(x)$. The sequence $f_0, f_1, \dots, f_k$ is called a Sturm chain for the polynomial $f(x)$.

**Sturm Theorem:** The number of distinct real zeros of a polynomial $f(x)$ with real coefficients in $(a, b)$ is equal to the excess of the number of sign change in the sequence $f_0(a), \dots, f_{k-1}(a), f_k(a)$ over the number of sign change in the sequence $f_0(b), \dots, f_{k-1}(b), f_k(b)$.

As an example, we illustrate Sturm theorem by isolating the real roots of:

$$x^3 - 3x^2 - 4x + 13 = 0 \qquad (16)$$

Using polynomial division with MAXIMA symbolic software, the Sturm chain is given by:

$$ \begin{aligned} f_0(x) &= x^3 - 3x^2 - 4x + 13, \\ f_1(x) &= 3x^2 - 6x - 4, \\ f_2(x) &= 14x/3 - 35/3, \\ f_3(x) &= 1/4, \end{aligned} \qquad (17) $$

Taking the entire real interval $x =] -\infty, \infty[$, the following table displays the procedure:

<table><caption>**Table 1:** Sturm method to count zeros in different intervals over the real axis. Note that sign change is along a horizontal direction whereas root counting is related to sign change alteration along a vertical direction. In effect, going vertically from -3 to -2, sign change goes from 3 to 2 thus there is a single root. In contrast, going vertically from 2 to 3, sign change goes from 2 to 0 thus there are two roots.</caption><thead><tr><th>x</th><th>f<sub>0</sub></th><th>f<sub>1</sub></th><th>f<sub>2</sub></th><th>f<sub>3</sub></th><th>Sign change</th></tr></thead><tbody><tr><td>-∞</td><td>-</td><td>+</td><td>-</td><td>+</td><td>3</td></tr><tr><td>-3</td><td>-</td><td>+</td><td>-</td><td>+</td><td>3</td></tr><tr><td>-2</td><td>+</td><td>+</td><td>-</td><td>+</td><td>2</td></tr><tr><td>0</td><td>+</td><td>-</td><td>-</td><td>+</td><td>2</td></tr><tr><td>2</td><td>+</td><td>-</td><td>-</td><td>+</td><td>2</td></tr><tr><td>3</td><td>+</td><td>+</td><td>+</td><td>+</td><td>0</td></tr><tr><td>+∞</td><td>+</td><td>+</td><td>+</td><td>+</td><td>0</td></tr></tbody></table>

<figure><img src="image_8.png" /><figcaption>**Fig.7:** (Color on-line) (Left) Polynomial $x^3 - 3x^2 - 4x + 13$ with one zero in (-3, -2) interval, and two in (2,3). (Right) Zoomed graph of the (2,3) interval showing clearly the roots.</figcaption></figure>

From Table 1, one negative root is in interval (-3, -2) and two positive roots are in (2,3) as confirmed by Fig. 7.

## 2. Routh array

Another method for counting and determination of roots is provided by Routh array algorithm [6]. This method is used in Control engineering and particularly in determination of system stability [6].
Consider the Linear Time Invariant (LTI) system whose transfer function denominator polynomial leads to the characteristic equation

$$s^n + a_1 s^{n-1} + a_2 s^{n-2} + \dots + a_n = 0 \quad (18)$$

To determine the Routh array, we first arrange the coefficients of the characteristic polynomial in two rows, containing respectively even and odd-numbered coefficients:

$$\begin{array}{l} s^n: \quad 1 \quad a_2 \quad a_4 \quad \dots \\ s^{n-1}: \quad a_1 \quad a_3 \quad a_5 \quad \dots \end{array} \quad (19)$$

We then add subsequent rows to complete the Routh array [6] having $n/2 + 1$ columns and $n + 1$ rows:

$$\begin{array}{llll} \text{Row} & n & s^n: & 1 \quad a_2 \quad a_4 \quad \dots \\ \text{Row} & n-1 & s^{n-1}: & a_1 \quad a_3 \quad a_5 \quad \dots \\ \text{Row} & n-2 & s^{n-2}: & b_1 \quad b_2 \quad b_3 \quad \dots \\ \text{Row} & n-3 & s^{n-3}: & c_1 \quad c_2 \quad c_3 \quad \dots \\ \vdots & & & \\ \text{Row} & 2 & s^2: & * \quad * \quad * \quad \dots \\ \text{Row} & 1 & s: & * \quad * \quad * \quad \dots \\ \text{Row} & 0 & s^0: & * \quad * \quad * \quad \dots \end{array} \quad (20)$$

The $b_i, c_i$ elements are calculated from the coefficients of the two previous rows as [6]:

$$b_1 = \frac{a_1 a_2 - a_3}{a_1}, \quad b_2 = \frac{a_1 a_4 - a_5}{a_1}, \quad b_3 = \frac{a_1 a_6 - a_7}{a_1}, \quad c_1 = \frac{b_1 a_3 - a_1 b_2}{b_1}, \quad c_2 = \frac{b_1 a_5 - a_1 b_3}{b_1}, \quad c_3 = \frac{b_1 a_7 - a_1 b_4}{b_1} \quad (21)$$

The rest of the elements are generated similarly by cross multiplying the coefficients of the two previous rows and proceeding downward until all subsequent coefficients are zero. Moreover missing coefficients in every row are replaced by zeros. The column beginning with $a_1$ is considered to be the first column of the array. The Routh array is obtained as a triangular matrix with non-zero elements followed by zeroes to complete the rows.

It can be shown that a necessary and sufficient condition for roots to have positive real parts i.e. located in the right hand plane [6] is by counting the number of sign changes in the first column of the array,

Example with ($n = 3$):

$$s^3 + 14s^2 + 41s - 56 = 0 \quad (22)$$

Routh table:

$$\begin{array}{llll} \text{Row} & 3 & s^3: & 1 \quad 41 \\ \text{Row} & 2 & s^2: & 14 \quad -56 \\ \text{Row} & 1 & s^1: & 45 \quad 0 \\ \text{Row} & 0 & s^0: & -56 \quad 0 \end{array} \quad (23)$$

One sign change in first column, thus we have a single zero in the right hand plane as verified by the factorization:

$$s^3 + 14s^2 + 41s - 56 = (s - 1)(s + 7)(s + 8) \quad (24)$$

Example with ($n = 4$):

$$s^4 - 10s^3 + 35s^2 - 50s + 24 = 0 \quad (25)$$

Routh table:

$$\begin{array}{llll} \text{Row} & 4 & s^4: & 1 \quad 35 \quad 24 \\ \text{Row} & 3 & s^3: & -10 \quad -50 \quad 0 \\ \text{Row} & 2 & s^2: & 30 \quad 24 \quad 0 \\ \text{Row} & 1 & s^1: & -42 \quad 0 \quad 0 \\ \text{Row} & 0 & s^0: & 24 \quad 0 \quad 0 \end{array} \quad (26)$$

There are four sign changes in first column, thus we have four zeroes in the right hand plane as verified by polynomial factorization:

$$s^4 - 10s^3 + 35s^2 - 50s + 24 = (s - 1)(s - 2)(s - 3)(s - 4) \quad (27)$$
## B. Real roots

Vibrations of mechanical systems provide situations leading to determination of oscillation modes linked to eigenvalue problems with real roots. As an example, consider vibrational planar modes of a triangular spring-mass system representing a molecule.

Planar mass displacements [7, 8] shown in Fig. 8 and given by $\mathbf{u}_i = (x_i, y_i)$, $i = 1, 2, 3$ yield the potential energy:

$$V = \frac{1}{2} \kappa \left[ \frac{1}{2}(x_2 - x_1) + \frac{\sqrt{3}}{2}(y_2 - y_1) \right]^2 + \frac{1}{2} \kappa \left[ -\frac{1}{2}(x_1 - x_3) + \frac{\sqrt{3}}{2}(y_1 - y_3) \right]^2 + \frac{1}{2} \kappa \left[ \frac{1}{2}(x_2 - x_3) \right]^2 \quad (28)$$

$\kappa$ is a mechanical constant characterizing the spring linking neighboring masses.

This mechanical system possesses six degrees of freedom with six modes of vibration.

<figure><img src="image_9.png" /><figcaption>**Fig.8:** (Color on-line) Spring-mass model of a triangular molecule [7, 8]. All atoms have mass $m$ and neighbors are linked pairwise with springs characterized by $\kappa$ along $\mathbf{r}_{23}$, $\mathbf{r}_{31}$ and $\mathbf{r}_{21}$. Small excursions of masses are along: $\mathbf{u}_i = (x_i, y_i)$, $i = 1, 2, 3$.</figcaption></figure>

The kinetic energy is given by:

$$T = \frac{1}{2} m [\dot{x}_1^2 + \dot{y}_1^2 + \dot{x}_2^2 + \dot{y}_2^2 + \dot{x}_3^2 + \dot{y}_3^2] \quad (29)$$

Euler-Lagrange equations [7, 8] of motion are given by:

$$\frac{d}{dt} \left( \frac{\partial T}{\partial \dot{q}_i} \right) + \frac{\partial V}{\partial q_i} = 0, \quad q_i = x_i, i = 1, 2, 3, \quad q_i = y_i, i = 4, 5, 6 \quad (30)$$

Substituting $q_i = a_i \exp(j\omega t)$ with $\omega = \sqrt{\frac{\lambda\kappa}{m}}$, Euler-Lagrange equations 30 yield the secular equation (obtained with MAXIMA):

$$D(\lambda) = \begin{vmatrix} \frac{1}{2} - \lambda & 0 & -\frac{1}{4} & -\frac{1}{4}\sqrt{3} & -\frac{1}{4} & \frac{1}{4}\sqrt{3} \\ 0 & \frac{3}{2} - \lambda & -\frac{1}{4}\sqrt{3} & -\frac{3}{4} & \frac{1}{4}\sqrt{3} & -\frac{3}{4} \\ -\frac{1}{4} & -\frac{1}{4}\sqrt{3} & \frac{5}{4} - \lambda & \frac{1}{4}\sqrt{3} & -1 & 0 \\ -\frac{1}{4}\sqrt{3} & -\frac{3}{4} & \frac{1}{4}\sqrt{3} & \frac{3}{4} - \lambda & 0 & 0 \\ -\frac{1}{4} & \frac{1}{4}\sqrt{3} & -1 & 0 & \frac{3}{4} - \lambda & -\frac{1}{4}\sqrt{3} \\ \frac{1}{4}\sqrt{3} & -\frac{3}{4} & 0 & 0 & -\frac{1}{4}\sqrt{3} & \frac{3}{4} - \lambda \end{vmatrix} = 0 \quad (31)$$
Expanding the determinant with MAXIMA, we get: $D(\lambda) = \lambda^6 - 6\lambda^5 + 45\lambda^4/4 - 27\lambda^3/4$. $D(\lambda) = 0$ roots are vibrational mode eigenvalues $\lambda = 0, \lambda = 3/2$ and $\lambda = 3$.

Confirmation of eigenvalues and their multiplicities is given by the plot of $D(\lambda)$ versus $\lambda$ displayed in fig. 9. The roots are: triple $\lambda = 0$, double $\lambda = 3/2$ and simple $\lambda = 3$ implying that $D(\lambda)$ can be factored as $D(\lambda) = \lambda^3(\lambda - 3/2)^2(\lambda - 3)$ as also checked with MAXIMA.

<figure><img src="image_10.png" /><figcaption>**Fig.9:** (Color on-line) Inspection of the $D(\lambda)$ curve locally around crossings with $\lambda$ axis provides roots location and multiplicities as: $\lambda = 0$, triple (cubic aspect), $\lambda = 3/2$, double (parabolic) and $\lambda = 3$, simple (linear)</figcaption></figure>

An alternative to the graphic solution is Routh criterion that determines the number of roots of a polynomial and their values.

Let us apply this method to $D(\lambda) = \lambda^6 - 6\lambda^5 + 45\lambda^4/4 - 27\lambda^3/4 = 0$ that is transformed after multiplying by 4 into: $4s^6 - 24s^5 + 45s^4 - 27s^3$

The Routh array is given by:

$$ 
\begin{array}{l c c c c c}
\text{Row 6 } s^6: & 4 & 45 & 0 & 0 & 0 \\
\text{Row 5 } s^5: & -24 & -27 & 0 & 0 & 0 \\
\text{Row 4 } s^4: & 40.5 & 0 & 0 & 0 & 0 \\
\text{Row 3 } s^3: & -27 & 0 & 0 & 0 & 0 \\
\text{Row 2 } s^2: & -81 & 0 & 0 & 0 & 0 \\
\text{Row 1 } s^1: & -162 & 0 & 0 & 0 & 0 \\
\text{Row 0 } s^0: & 0 & 0 & 0 & 0 & 0
\end{array}
\qquad (32) $$

There are three sign changes in column 1 (Row 0 does not count in this case), implying we have three zeroes in the right hand plane agreeing with polynomial factorization $s^3(s - 3/2)^2(s - 3)$ that yields triple $s=0$, double $s = 3/2$ and simple $s = 3$.

### C. Complex roots

Moving on to the case of complex roots, consider the following $n = 4$ polynomial with real coefficients:

$$ s^4 + 2s^3 + 3s^2 + 4s + 5 = 0 \qquad (33) $$

The Routh array is given by:
$$
\begin{array}{l}
\text{Row 4 } s^4: 1 \ 3 \ 5 \\
\text{Row 3 } s^3: 2 \ 4 \ 0 \\
\text{Row 2 } s^2: 1 \ 5 \ 0 \\
\text{Row 1 } s^1: -6 \ 0 \ 0 \\
\text{Row 0 } s^0: 5 \ 0 \ 0
\end{array}
\qquad (34)
$$

There are two sign changes in column 1, implying presence of two zeroes in the right hand plane.

According to the fundamental theorem of algebra stating that if all coefficients of a polynomial are real, then the complex roots occur in conjugate pairs, meaning that the polynomial can be factored as $(z-z_1)(z-z_2)(z-z_3)(z-z_4)$ with $z_1, z_2$ in the right hand plane and $z_3, z_4$ in the left hand plane with: $z_4 = z_1^*$ and $z_3 = z_2^*$.

The determination of $z_1, z_2$ can be done numerically in double precision with Jenkins-Traub algorithm [9] (JTA). JTA consists of three stages with a first, called zero-shift, operating to accentuate smaller roots. A second stage implements a fixed-shift operation which isolates one zero followed by a deflation operation in order to reduce the polynomial order. A third stage is a variable-shift operation which converges very quickly toward a root. The third stage may also be viewed as a Newton-Raphson iteration performed on a sequence of rational functions. Another interpretation of the third stage is its equivalence to a generalized Rayleigh quotient [5] iteration applied to the companion matrix of the polynomial.

In the ordinary version of matrix inverse iteration [5], a constant value $\alpha$ is used at each step to approximate some eigenvalue $\lambda$. In Rayleigh quotient iteration method, the current iterate $\mathbf{x}_n$ is used to improve the value of $\alpha$ at each step by setting $\alpha = \mathbf{x}_n^T A \mathbf{x}_n / \mathbf{x}_n^T \mathbf{x}_n$.

Function $R(\mathbf{x}) = \mathbf{x}^T A \mathbf{x} / \mathbf{x}^T \mathbf{x}$ is called the Rayleigh quotient. It can be shown generally that if $\mathbf{x}$ is a good approximation to an eigenvector, then $R(\mathbf{x})$ is a good approximation of the associated eigenvalue. The range of all Rayleigh quotients is called the numerical range of $A$ that contains the spectrum of $A$ (i.e., all eigenvalues of $A$).

Additionally, Jenkins and Traub [9] have shown their algorithm to be globally convergent and, furthermore, that its convergence toward a root is faster than second-order.

Thus, the complex roots in double precision are obtained as:

$$
\begin{align*}
z_1 &= 0.287815479557648 + 1.416093080171908 i, & z_4 &= z_1^* \\
z_2 &= -1.2878154795576482 + 0.8578967583284902 i, & z_3 &= z_2^*
\end{align*}
$$

IV. EXTREMUM PROBLEMS IN HIGHER DIMENSIONS

Finding zeroes $\mathbf{x}^*$ of the partial derivatives $f_i' = (\frac{\partial f}{\partial x_i})$ of a multidimensional function $f(\mathbf{x})$ amounts to search for its extremum. The latter is maximum or minimum depending on whether the Hessian matrix $H(f)$ (or local convexity) whose components are defined by the second derivatives: $H_{ij} = (\frac{\partial^2 f}{\partial x_i \partial x_j})$ is positive or negative at $\mathbf{x}^*$.

A. Powell Method and minimization

Powell method is a procedure [10, 11] for carrying out minimization in N-dimensional space with no derivatives. It performs repeated application of a 1D optimization procedure. Thus Powell method is considered as a direction-set method for minimization [10] which means that if we start at a point P in N-dimensional space, and proceed from there along a vector direction n, then any function of N variables $f(P)$ can be minimized along the line n by the 1D optimization method. Powell method belongs to a type of procedures that optimize $f(\mathbf{x})$ through iterative updating of a candidate vector $\mathbf{x}^*$, where updates are given by $\mathbf{x}^* \leftarrow \mathbf{x}^* + \lambda^* e_j$, with $\lambda^*$ defined as a 1D optimum of $f(\mathbf{x}^* + \lambda^* e_j)$. Vectors $\{e_j\}$ make a basis for $\mathcal{R}^d$, and cycles are performed through basis $j = 1, \dots, d$ as $\mathbf{x}^*$ is being updated.

A major drawback occurs when {$e_j$} basis is not well-aligned with function $f$ contours. The algorithm converges very slowly making erratic walks on $f$ surface. The optimal basis for a given starting point $\mathbf{x}^*$ is obtained from the eigenvectors [10, 11] of $H(f)_x^*$, where $H(F)_u$ is the Hessian of $f$ evaluated at $u$ (i.e., the $d \times d$ matrix with value $(\frac{\partial^2 f}{\partial x_i \partial x_j})_u$). The set of directions is optimal when $f$ is quadratic, the exact minimum being reached in at most $d$ iterations [10, 11].
Since Powell method [10, 11] does not use first nor second derivatives of $f$, the optimal basis cannot be used. Powell method uses instead an approximation to the optimal basis updated after every full cycle of the algorithm. Let $x_0^i$ be an initial guess. The $i^{th}$ cycle of Powell method is as follows:

For $j = 1, \dots, d$:

$$ \lambda^* = \max_{\lambda} f(x_{j-1}^i + \lambda e_j). $$

$$ \text{If } j \le d, \quad x_{j+1}^i = x_{j-1}^i + \lambda^* e_j. $$

$$ \text{If } j < d, \text{ set } e_j = e_{j+1}. $$

$$ \text{If } j = d, \text{ set } x_0^{i+1} = x_d^i, e_d = x_d^i - x_0^i. $$

The main step is that basis vectors $e_j$ are re-written with the increment in $x$ for every cycle through a basis. This vector points in a direction where the function is known to increase substantially, so it should be mostly aligned with, rather than perpendicular to the contours of $f$ [10, 11]. Note that after every set of $d$ iterations of the above procedure, $\{e_j\}$ basis will have been completely re-written (this occurs after $d^2$ 1D optimizations are achieved). Moreover, it is important to orthonormalize $\{e_d, \dots, e_1\}$ before the next set of $d$ iterations begins. Moreover, orthonormalizing in reverse order, induces the next set to search along the direction found in the final cycle of the preceding one.

Using Powell algorithm [10], we solve two 3D examples [12]:

1. $f(x_1, x_2, x_3) = 3 - 4x_1 + 2x_1^2 - 4x_2 + 2x_1x_2 + 2x_2^2 + 4x_3 - 2x_1x_3 - 2x_2x_3 + 2x_3^2$ with a global minimum: $x^* = (0.500000, 0.500000, -0.500000), f(x^*) = 0$.

The steps to minimize $f(x_1, x_2, x_3)$ are:

Start point: 1.5,1.5,2.5

Iterations: 250 (the maximum allowed by the algorithm [10])

Minimum found at: $x^* = (0.484187, 0.496851, -0.529027), f(x^*) = 0.001204$

In order to have faster convergence, the start point is changed, following a hint of negating $x_3$ from the above results:

Start point: 1, 1, -5

Iterations: 8

Minimum found at: $x^* = (0.500000, 0.500000, -0.500000), f(x^*) = 0.000000$

2. Hartmann function: $H_{3,4}(x) = - \sum_{i=1}^{4} \alpha_i \exp \left( - \sum_{j=1}^{3} A_{ij}(x_j - P_{ij})^2 \right)$

$$ \text{with } \alpha = [1, 1, 2, 3, 3, 2]^T, A = \begin{bmatrix} 3.0 & 10 & 30 \\ 0.1 & 10. & 35 \\ 3.0 & 10.0 & 30.0 \\ 0.1 & 10 & 35 \end{bmatrix}, P = 10^{-4} \begin{bmatrix} 6890 & 1170 & 2673 \\ 4699 & 4387 & 7470 \\ 1091 & 8732 & 5547 \\ 381 & 5743 & 8828 \end{bmatrix} $$

The global minimum is: $x^* = (0.114614, 0.555649, 0.852547), H_{3,4}(x^*) = -3.86278$.

The results for $H_{3,4}(x)$ are:

Start point: 1.5,1.5,2.5

Iterations: 4

Minimum found at: $x^* = (0.114690, 0.555646, 0.852544), H_{3,4}(x^*) = -3.862782$

## B. Broyden method and fitting

Broyden method is similar to Newton-Raphson, in the sense it updates approximations to both solution and Jacobian over iteration runs. However, it is considered also as a multidimensional version of the secant method.

Broyden update for solving a nonlinear equation is given by:

$$ F(x) = 0, $$

where $F: \mathbb{R}^n \rightarrow \mathbb{R}^n$ is differentiable. We denote its Jacobian as $J(x) \in \mathbb{R}^{n \times n}$. We need to iteratively update the approximate Jacobian $\{B_k\}$ in order to approach the true Jacobian matrix, while updating parameters $\{x_k\}$ as:

$$ x_{k+1} = x_k - B_k^{-1} F(x_k). $$
The alternative matrix $B_k$ should satisfy the secant condition, that is

$$B_k u_k = y_k, \text{ with } u_k = x_{k+1} - x_k, y_{k+1} = F(x_{k+1}) - F(x_k). \quad (35)$$

However, we could not obtain a definite solution with only eq. (35) because it is an under-determined linear system with an infinite number of solutions, since there are less constraints than intrinsic parameters.

The most common procedure is to select a Jacobian approximation $B_{k+1}$ as close as possible to $B_k$ in the absolute norm sense, i.e. $||B_{k+1} - B_k||$ should be as small as possible. Broyden update [10, 13] defined as follows

$$B_{k+1} = B_k + \frac{(y_k - B_k u_k) u_k^\top}{u_k^\top u_k}, \quad (36)$$

is obtained by solving the following problem:

$$\min_B \|B - B_k\|_F \quad \Rightarrow \quad B u_k = y_k.$$

Updating the inverse of matrix $G_k := B_k^{-1}$ is as follows [10, 13]:

$$G_{k+1} = G_k - \frac{(G_k y_k - u_k) u_k^\top G_k}{u_k^\top G_k y_k}.$$

Broyden updating operations are done in $O(n^2)$ time [13] (rather than $O(n^3)$). For large problems, the complexity of matrix $B_k$ gets more involved as one proceeds, imposing memory limit to keep only a few recent updates to the initial (often diagonal) Jacobian approximation, might resolve this issue.

It is interesting to realize that Broyden method converges in a superlinear fashion [10]. In addition, Broyden method works in over-determined (more constraints than parameters) and under-determined problems [1, 13] (less constraints than parameters) finding always an solution $x_{opt}$ that is close to the ideal one $x_{ideal}$ in the absolute norm sense, i.e. $||x_{opt} - x_{ideal}||$ is minimized.

In order to test the robustness of Broyden method, we consider a synthetic model, where an analytic function with a number of parameters is chosen and perturbed by an additive White noise of variable amplitude.

The challenge for Broyden method, is determination of the right set of parameters in presence of White noise.

Let us consider a model having 3 parameters: $a, b, c$ such that: $f(x) = \frac{a}{[1+b(1-x/c)^2]}$ and a noise source given by an amplitude $wU(0,1)$ with $U(0,1)$ a uniform random number generator over the interval $[0,1]$. This means every function value $f(x_i)$ is transformed into $f(x_i) + wU(0,1)$.

Since we have a model with three parameters, the solution is $y(x)$ is subject to three relations to be satisfied by the algorithm in order to find a fitting solution $y(x)$ as close as possible to $f(x)$ without over or under determination.

These relations are:

1. First condition: Least sum of moduli between noisy and model points: $f_1 = \sum_i |y(x_i) - f(x_i)|$

2. Second condition: Fitting to a value at $x = x_1$: $f_2 = |y(x = x_1) - f(x = x_1)|$

3. Third condition: Fitting to a value at $x = x_2$: $f_3 = |y(x = x_2) - f(x = x_2)|$

Model parameters are chosen as: $a = 4, b = 5, c = 6$.

In the noiseless case, Broyden algorithm delivered $a = 2.8479817, b = 1.9756421, c = 5.837635$ with a start point (1,1,1). Comparison between model and fit found by Broyden is shown in Fig.10.

In the noisy case with $w = 0.1$, Broyden algorithm delivered $a = 2.9574127, b = 1.9596522, c = 5.9582214$ with a start point (3,2,6). Comparison between noisy model and fit found by Broyden is shown in Fig.11.

In the noisy case with $w = 0.5$, Broyden algorithm delivered $a = 3.8551593, b = 2.6939712, c = 6.057486$ with a start point (3,2,1). Comparison between noisy model and fit found by Broyden is shown in Fig.12.

# V. APPLICATIONS

## A. Evaluation of mathematical constants

### Evaluation of $\pi$

Consider the root equation $\cos x = -1$ whose solution is $x = \pi$. Thus any equation whose solution is $\pi$ can be used in order to determine $\pi$ by root finding methods. For instance $\sin x = 0$ or $\tan x = 0$...
<figure><img src="image_11.png" /><figcaption>Fig.10: Comparison between model $f(x) = \frac{a}{[1+b(1-x/c)^2]}$ and fit in absence of noise with a start point (1,1,1).</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Fig.11: Comparison between noisy model and fit in presence of noise with an amplitude $w = 0.1$</figcaption></figure>

• Evaluation of Euler e

Consider the root equation $\ln x = 1$ whose solution is $x = e$. Thus any equation whose solution is $e$ can be used in order to determine $e$.

• Evaluation of $\ln 2$

Consider the root equation $e^x = 2$ whose solution is $x = \ln 2$. Thus any similar equation whose solution is $\ln 2$ can be used to determine it.

• Evaluation of the Golden mean $\phi$

Consider the root equation $x^2 = x + 1$ whose solutions are $x = \phi$ and $x = -1/\phi$. Thus any equation whose solution is $\phi$ can be used in order to determine it. For instance, $\phi$ can be expressed by an infinitely iterated form [14]:

$$ \phi = \sqrt{1 + \sqrt{1 + \sqrt{1 + \sqrt{1 + \dots}}}} \qquad (37) $$

which can be rewritten as a fixed point iteration: $x_{n+1} = \sqrt{1 + x_n}$ and solved by: $\lim_{n \to \infty} x_n = x^* \text{ resulting in } (x^*)^2 = x^* + 1$.

## B. Arithmetic operations and evaluation of functions

• Division operation
<figure><img src="image_13.png" /><figcaption>Fig.12: Comparison between noisy model and fit in presence of noise with an amplitude w = 0.5</figcaption></figure>

Given a real constant A consider the root equation $f(x) = 1/x-A$ whose solution is $x^* = 1/A$ then $f'(x) = -1/x^2$ and Newton-Raphson iteration is $x_{n+1} = x_n(2 - Ax_n)$ containing only multiply and subtract operations.

* **Square root operation**

Given a real constant A consider the root equation $f(x) = x^2 - A$ whose solution is $x^* = \sqrt{A}$ then $f'(x) = 2x$ and Newton-Raphson iteration is $x_{n+1} = \frac{1}{2}(x_n + A/x_n)$ containing add, multiply and divide operations.

* **Inverse cosine function**

Given a real constant A consider the root equation $f(x) = \cos x - A$ whose solution is $x^* = \cos^{-1} A$ then $f'(x) = -\sin x$ and Newton-Raphson iteration is $x_{n+1} = x_n + (\cos x_n - A/\sin x_n)$ containing no inverse trigonometric operations.

* **Inverse tangent function**

Given a real constant A consider the root equation $f(x) = \tan x - A$ whose solution is $x^* = \tan^{-1} A$ then $f'(x) = 1 + \tan^2 x$ and Newton-Raphson iteration is $x_{n+1} = x_n + (\tan x_n - A)/(1 + \tan^2 x_n)$ that can be rewritten as: $x_{n+1} = x_n + A \cos^2 x_n - \sin(x_n) \cos(x_n)$ containing no inverse trigonometric operations.

C. Signal reconstruction from its zeroes

In power systems, acoustics, mechanical vibrations, signal processing... zero detection of a signal can be performed by electronic circuits or special sensors from its sign change [15, 16]. Afterwards, interpolation methods are employed to construct the signal from its zeroes according to:

* Lagrange method requiring signal values only:

Given a set of $n + 1$ points $a \le x_0 < x_1 < ... < x_{n-1} < x_n \le b$, Lagrange interpolation polynomial is:

$$p_L(x) = \sum_{k=0}^{n} f(x_k)\ell_k(x), \quad \ell_k(x) = \prod_{j=0, j \neq k}^{n} \left( \frac{x - x_j}{x_k - x_j} \right). \qquad (38)$$

* Hermite method requiring signal values and its derivatives

Hermite method requires a function $f(x)$ and its derivatives $f_i^{(l)}$ at the given points $x_i, i = 0, ..m$ with $x_0 < x_1.. < x_m$ and $k = 0...n_i - 1$ such that the interpolating polynomial writes: $p_H(x) = \sum_{i=0}^{m} \sum_{k=0}^{n_i-1} f_i^{(k)} L_{ik}(x)$ where $L_{ik}(x)$ are generalized Lagrange [5] polynomials and defined from the auxiliary polynomials:

$$\ell_{ik}(x) = \frac{(x - x_i)^k}{k!} \prod_{j=0, j \neq i}^{m} \frac{(x - x_j)}{x_i - x_j}, \quad 0 \le i \le m, \quad 0 \le k \le n_i \qquad (39)$$

to get recursively, starting from $L_{i,n_i-1}(x) = \ell_{i,n_i-1}, i = 0...m$, and: $L_{ik}(x) = \ell_{ik}(x) - \sum_{\nu=k+1}^{n_i-1} \ell_{ik}^{(\nu)}(x) L_{i\nu}(x)$, $k = n_i - 2, n_i - 3, 0$ obeying the constraints [5]:

$$p_H(x_i) = f_i, \quad p_H^{(1)}(x_i) = f_i^{(1)}, \quad p_H^{(m)}(x_i) = f_i^{(m)}, \quad i = 0, ..k \qquad (40)$$
where the upper index (i), i = 0...k indicates derivation order.

* Trigonometric interpolation (Fourier series) requiring signal values only

* Spline method requiring signal values only

* B-spline method requiring signal values and other properties

## D. Filter design from zeroes and poles

In signal processing, telecommunication... an input signal (message) $x(t)$ at time $t$ goes through a channel (or medium, apparatus...) and emerges from it as an output $y(t)$ (cf. Fig. 13).

<figure><img src="image_14.png" /><figcaption>Fig.13: Impulse response of a signal $x(t)$ (input voltage) and output signal $y(t)$ with $h(t)$ the impulse response relating them.</figcaption></figure>

In a linear channel, input and output [2] are related by:

$$y(t) = \int_{-\infty}^{\infty} x(t')h(t, t')dt' \qquad (41)$$

Adopting time-invariance of the system through relation $h(t, t') = h(t - t')$, the operator $O[\cdot]$ representing the system effect on its input is transformed into a convolution operation (*) [2] of $x(t)$ and $h(t)$ functions as:

$$y(t) = \int_{-\infty}^{\infty} x(\tau)h(t - \tau)d\tau = x(t) * h(t), \quad \text{or} \quad y(t) = \int_{-\infty}^{\infty} h(\tau)x(t - \tau)]d\tau = h(t) * x(t) \qquad (42)$$

In other words, a message $x(t)$ transmitted through a Linear Time Invariant (LTI) medium is convoluted with the channel impulse response $h(t)$ resulting in an output message $y(t)$.

More generally, this process can be described by a time domain Ordinary Differential Equation (ODE) relating input $x(t)$ to output $y(t)$ according to Ref. [2, 17]:

$$\sum_{k=0}^{n} a_k \frac{d^k}{dt^k} y(t) = \sum_{k=0}^{m} b_k \frac{d^k}{dt^k} x(t) \qquad (43)$$

Taking Laplace [2] transform of eq. 43, it can be represented as an algebraic equation in the s domain:

$$Y(s) \left[ \sum_{k=0}^{n} a_k s^k \right] = X(s) \left[ \sum_{k=0}^{m} b_k s^k \right] \qquad (44)$$

and its transfer function $H(s)$ the Laplace transform of impulse response $h(t)$ is rational

$$H(s) = \frac{Y(s)}{X(s)} = \frac{\sum_{k=0}^{m} b_k s^k}{\sum_{k=0}^{n} a_k s^k} = \frac{b_m \prod_{k=0}^{m} (s - s_{z_k})}{a_n \prod_{k=0}^{n} (s - s_{p_k})} \qquad (45)$$

where $s_{z_k}$, ($k = 1, 2, \dots, m$) are the roots of the numerator and $s_{p_k}$, ($k = 1, 2, \dots, n$) are the roots of the denominator.

Some example filters are given below:

* **One-Pole Filter:** Choosing a pole $s_p = -\omega_c$, the cutoff angular frequency, we write:

$$H(s) = \frac{\omega_c}{s + \omega_c}, \quad H(j\omega) = \frac{\omega_c}{j\omega + \omega_c} = \frac{1}{1 + j\omega/\omega_c}$$

$$|H(j\omega)| = \frac{1}{\sqrt{1 + (\omega/\omega_c)^2}}$$
## Two-Pole Filter

Assuming a repeated pole $s_{p_1} = s_{p_2} = -\omega_c$, we get:

$$H(s) = \frac{\omega_c^2}{(s + \omega_c)^2}, \quad H(j\omega) = \frac{\omega_c^2}{(j\omega + \omega_c)^2} = \frac{1}{(1 + j\omega/\omega_c)^2} \\ |H(j\omega)| = \frac{1}{1 + (\omega/\omega_c)^2}$$

Alternatively, if we choose a pair of complex conjugate poles $s_{p_1,p_2} = \omega_c e^{\pm j3\pi/4}$, then we have

$$H(s) = \frac{\omega_c^2}{s^2 - 2\omega_c \cos(3\pi/4)s + \omega_c^2}, \quad H(j\omega) = \frac{1}{1 - (\omega/\omega_c)^2 + \sqrt{2}(j\omega/\omega_c)} \\ |H(j\omega)| = \frac{1}{\sqrt{1 + (\omega/\omega_c)^4}}$$

## n-order Butterworth Filter

The magnitude of the frequency response function of an n-th order Butterworth filter takes the following form:

$$|H(j\omega)|^2 = \frac{1}{1 + (\omega/\omega_c)^{2n}} = \frac{1}{1 + (\omega^2/\omega_c^2)^n}$$

The corresponding impulse response function $h(t) = \mathcal{F}^{-1}[H(j\omega)]$ where $\mathcal{F}^{-1}$ is the inverse Fourier transform. In order for $h(t)$ to be real, $H(j\omega)$ must satisfy

$$H^*(j\omega) = H(-j\omega)$$

Thus

$$|H(j\omega)|^2 = H(j\omega)H^*(j\omega) = H(j\omega)H(-j\omega) = \frac{1}{1 + (\omega^2/\omega_c^2)^n}$$

Replacing $j\omega$ by $s$; i.e., $\omega^2 = -s^2$, the transfer function $H(s)$ of the filter is obtained as:

$$|H(s)|^2 = H(s)H(-s) = \frac{1}{1 + (-s^2/\omega_c^2)^n} = \frac{\omega_c^{2n}}{\omega_c^{2n} + (-s^2)^n} = \frac{(-1)^n \omega_c^{2n}}{s^{2n} + (-1)^n \omega_c^{2n}}$$

The poles of $|H(s)|^2$ are:

$$|s_k| = \omega_c, \quad \angle s_k = \frac{(2k+1)\pi}{2n} + \frac{\pi}{2}, \quad k = 0, 1, \dots, 2n-1$$

This yields:

$$|H(s)|^2 = H(s)H(-s) = \frac{(-1)^n \omega_c^{2n}}{\prod_{k=0}^{2n-1} (s - s_k)} \quad (46)$$

## E. Stability of control and feedback systems

Consider the transfer function:

$$H(s) = \frac{Y(s)}{X(s)} = \frac{b_0 s^m + b_1 s^{m-1} + \dots + b_{m-1} s + b_m}{a_0 s^n + a_1 s^{n-1} + \dots + a_{n-1} s + a_n} \quad (47)$$

with $m \le n$.
Routh stability criterion is related to the denominator of $H(s)$ such its poles i.e. zeroes of $X(s) = \sum_{k=0}^{n} a_k s^k$ should be in the left hand plane.

As an example consider the quadratic characteristic polynomial:

$$a_0s^2 + a_1s + a_2 = 0 \qquad (48)$$

The array of Routh coefficients is given by:

$$ \begin{array}{l} \text{Row 2 } s^2: a_0 \ a_2 \ 0 \\ \text{Row 1 } s^1: a_1 \ 0 \ 0 \\ \text{Row 0 } s^0: a_2 \ 0 \ 0 \end{array} \qquad (49) $$

From the first column, we conclude that the system is stable when all $a_i$ are positive.

A quartic example is given by the transfer function:

$$H(s) = \frac{A}{s(s^2 + s + 1)(s + 2) + A} \qquad (50)$$

The denominator is: $s^4 + 3s^3 + 3s^2 + 2s + A = 0$.

The array of Routh coefficients is given by:

$$ \begin{array}{l} \text{Row 4 } s^4: \quad 1 \quad 3 \quad A \quad 0 \\ \text{Row 3 } s^3: \quad 3 \quad 2 \quad 0 \quad 0 \\ \text{Row 2 } s^2: \quad 7/3 \quad A \quad 0 \quad 0 \\ \text{Row 1 } s^1: \quad 2 - 9A/7 \quad 0 \quad 0 \quad 0 \\ \text{Row 0 } s^0: \quad A \quad 0 \quad 0 \quad 0 \end{array} \qquad (51) $$

From the first column we get the stability conditions: $A > 0$ and $2 - 9A/7 > 0$ i.e. $0 < A < 14/9$.

There exists another procedure to test stability based on the Routh-Hurwitz [6] criterion.

Given a characteristic polynomial of the form:

$$s^n + p_1 s^{n-1} + p_2 s^{n-2} + \dots + p_n = 0 \qquad (52)$$

Form the determinants $\Delta_1, \Delta_2, \Delta_3, \dots, \Delta_n$ with $\Delta_k$ of size $k \times k$ with first element on the main diagonal $p_1$ and last element $p_k$. Each row starts by an odd element followed by succeeding elements until last element is $p_1$ after which the elements are 1 then a succession of zeroes to complete the $k$ row elements:

$$ \Delta_1 = a_1, \Delta_2 = \begin{vmatrix} p_1 & 1 \\ p_3 & p_2 \end{vmatrix}, \Delta_3 = \begin{vmatrix} p_1 & 1 & 0 \\ p_3 & p_2 & p_1 \\ p_5 & p_4 & p_3 \end{vmatrix}, \dots, \Delta_k = \begin{vmatrix} p_1 & 1 & 0 & 0 & 0 & 0 & \cdots \\ p_3 & p_2 & p_1 & 1 & 0 & 0 & 0 & \cdots \\ p_5 & p_4 & p_3 & p_2 & p_1 & 1 & 0 & \cdots \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ p_{2k-1} & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & p_k \end{vmatrix}, \dots \qquad (53) $$

System stability is when all roots are in the left hand plane with $\Delta_1 > 0, \Delta_2 > 0, \Delta_3 > \cdots \Delta_n > 0$.

[1] J. R. Rice, *Numerical Methods, Software and Analysis*, Second Edition, Academic Press, New-York (1993).

[2] R. Wang, *Introduction to Orthogonal Transforms*, Cambridge University Press, New-York (2011).

[3] S. D. Conte, C. de Boor, *Elementary Numerical Analysis*, McGraw-Hill, New York, (1980).

[4] H. G. Schuster and W. Just, *Deterministic Chaos: An Introduction*, Fourth Edition, Wiley-VCH (2005).

[5] C. D. Meyer, *Matrix Analysis and Applied Linear Algebra*, SIAM, Philadelphia (2000).

[6] G. F. Franklin, J. D. Powell and A. Emami-Naeini, *Feedback control of dynamic systems*, 6th edition, Pearson Education Ltd, London (2010)

[7] E. Bright Wilson Jr., J.C. Decius and P. C. Cross, *Molecular Vibrations*, McGraw-Hill, New-York (1955).

[8] A. Nussbaum, Am. J. Phys., **36**, 529 (1968).

[9] M. A. Jenkins and J. F. Traub, *A three-stage variable-shift iteration for polynomial zeros and its relation to generalized Rayleigh iteration*, Numer. Math., **14**, 252, (1969/70), M. A. Jenkins and J. F. Traub, *A three-stage algorithm for real polynomials using quadratic iteration*, SIAM J. Numer. Analysis, **7**, 545 (1970), M. A. Jenkins, *Algorithm 493: Zeros of a Real Polynomial*, ACM TOMS, **1**, 178. (1975), M. A. Jenkins and J. F. Traub, *Algorithm 419: Zeros of a Complex Polynomial*, Comm. ACM, **15**, 97 (1972)
[10] W. H. Press, W. T. Vetterling, S. A. Teukolsky and B. P. Flannery, *Numerical Recipes in C: The Art of Scientific Computing*, Third Edition, Cambridge University Press, New-York (2007).

[11] D.E. Potter, *Computational physics*, John Wiley & Sons, London, New York [1973].

[12] M. R. Garey and D. S Johnson, *Computers and Intractability: A guide to the theory of NP completeness*, Bell Telephone Laboratories (1979).

[13] J. Stoer and R. Bulirsch, *Introduction to Numerical Analysis*, Second Edition, Springer-Verlag, (New York, 1992).

[14] D.E. Knuth, *The Art of Computer Programming: Semi-numerical Algorithms*, Vol. 2, Addison-Wesley (1997).

[15] T. H. O'Dell, *Circuits for Electronic Instrumentation*, Cambridge University Press, New-York (2005).

[16] A. B. Carlson and P. B. Crilly, *Communication systems: An Introduction to Signals and Noise in Electrical Communication*, 5th Edition, McGraw-Hill, New York (2010).

[17] G. F. Arfken and H. J. Weber, *Mathematical methods for physicists*, Sixth Edition, Elsevier Academic Press (London, 2005).