> For the complete documentation index, see [llms.txt](https://jona-projects.gitbook.io/kiwicalc/llms.txt). Markdown versions of documentation pages are available by appending `.md` to page URLs; this page is available as [Markdown](https://jona-projects.gitbook.io/kiwicalc/numerical-methods/multi-root-algorithms.md).

# Multi-Root Algorithms

**Durand-Kerner Method**

The Durand-Kerner method, also known as the Weierstrass method is an iterative approach for finding all of the real and complex roots of a polynomial. It was first discovered by the German mathematician Karl Weierstrass in 1891, and was later discovered by Durand(1960) and Kerner (1966). This method requires the function and a collection of its coefficients.

Here are some examples:

```python
func = lambda x: x ** 4 - 16
coefficients = [1, 0, 0, 0, -16]
print(durand_kerner(func, coefficients))               
```

Output:

```python
{(2+0j), -2j, 2j, (-2+0j)}      
```

You can also use the `durand_kerner2()` method if you only have the coefficients. For instance:

```python
print(durand_kerner2([1, 0, 0, 0, -16]))           
```

You can also run it from a  polynomial

**Aberth-Ehrlich Method**

Aberth's method is another method for simultaneous calculations of multiple roots of a polynomial function. It was first developed in 1967 and it's named after the mathematicians Oliver Aberth and Louis W. Ehrlich. Overall, in most implementations it is considered faster than the Durand-Kerner method, as it converges faster to the roots, and thus less iterations are performed. This implementation of the method consists of several steps:

1. creating $$n$$ complex approximations for the n roots, each one will convert to a different root. These approximations will be placed evenly on a circle on the complex plane that its center will be $$(0,0)$$, and its radius will be determined by the following formula:

$$
R = \sqrt\[n]{\left| \frac{p\_0}{p\_n} \right|}
$$

, where $$n$$ is the highest power in the expression, $$p\_0$$ is the coefficient of the free number of the expression, and $$p\_n$$ is the coefficient of the highest power in the expression.

2. For each approximation to the root, we will compute an offset, by the following formula: 

$$
w\_k = \frac{\frac{p(z\_k)}{p'(z\_k)}}{1 - \frac{p(z\_k)}{p'(z\_k)} \sum\_{j \neq k} \left( \frac{1}{z\_k - z\_j} \right)}
$$

3. Subtract each offset from its corresponding root. This will bring the approximations closer to the roots.
4. Repeat steps 2 and 3 until all the approximations have converged to the solutions.

The method requires the original function, its derivative, and the coefficients of the function.It will return a set of the complex approximations of the roots. Here is the signature of the method:

```python
 def aberth_method(f_0: Callable, f_1: Callable, coefficients, epsilon: float = 0.000001, nmax: int = 100000) -> set:
```

Here are some examples of using this method:                   

```python
func = lambda x: 5 * x ** 4 - 1
derivative = lambda x: 20 * x ** 3
coefficients = [5, 0, 0, 0, -1]
print(aberth_method(func, derivative, coefficients))

# output:
# {(0.66874+0j), 0.66874j, -0.66874j, (-0.66874+0j)}         
```

#### Generalized Newton

The generalized newton's method is a numerical method for finding the solutions of systems of nonlinear equations. Here are the steps of the underlying algorithm:

1. Creating $$X$$ - a vector of initial guesses to the solutions
2. Computing $$J(X)$$ - a corresponding jacobian matrix from the system of equations
3. If all of the guesses have converged to the appropriate solutions, return them and exit
4. Otherwise, modify the next set of solutions: $$X\_{n+1} = X\_n - J(X\_n)^{-1}F(X\_n)$$
5. Repeat phases 3 and 4 until the guesses have converged to the solutions.

For example, given this system of equations: 

$$
\left{
\begin{array}{l}
x^2 + y^2 = 25 \\
2x + 3y = 18
\end{array}
\right.
$$

We can extract $$F(X)$$, a system of functions: 

$$
\left{
\begin{array}{l}
f\_1(x, y) = x^2 + y^2 - 25 \\
f\_2(x, y) = 2x + 3y - 18
\end{array}
\right.
$$

Hence the jacobian matrix $$J(X)$$ will be:  $$\left\[ \begin{array}{cc}  \frac{\partial f\_1}{\partial x} & \frac{\partial f\_1}{\partial y} \  \frac{\partial f\_2}{\partial x} & \frac{\partial f\_2}{\partial y}  \end{array} \right]$$

Lets compute the partial\_derivatives: $$\frac{\partial f\_1}{\partial x} = 2x$$, $$\frac{\partial f\_1}{\partial y} = 2y$$ , $$\frac{\partial f\_2}{\partial x} = 2$$ $$\frac{\partial f\_2}{\partial y} = 3$$\
And therefore the \`J(X)\` will be:

$$
J(X) = \begin{bmatrix}
2x & 2y \\
2 & 3
\end{bmatrix}
$$

Now we need to choose two guesses for the solutions for \`x\` and \`y\`:\
$$X\_0 = \begin{bmatrix}  1 \  2  \end{bmatrix}$$

We can check whether the current guesses are accurate by checking if $$F(X\_0) \approx \begin{bmatrix}  0 \  0  \end{bmatrix}$$.

However $$F(X\_0) = F\left( \begin{bmatrix}  1 \  2  \end{bmatrix} \right) = \begin{bmatrix}  -20 \  -10  \end{bmatrix}$$ and therefore we have to keep iterating.

In order to update the guesses and find $$X\_1$$ we need to apply $$X\_1 = X\_0 - J(X\_0)^{-1}F(X\_0)$$, And for that we need to compute $$J(X\_0)$$ and $$J(X\_0)^{-1}$$.

$$
J(X\_0) = \begin{bmatrix}
2 & 4 \\
2 & 3
\end{bmatrix}
$$

$$
J(X\_0)^{-1} = \begin{bmatrix}
-1.5 & 2 \\
1 & -1
\end{bmatrix}
$$

And therefore:

$$
X\_1 = \begin{bmatrix}
1 \\
2
\end{bmatrix} - \begin{bmatrix}
-1.5 & 2 \\
1 & -1
\end{bmatrix} \begin{bmatrix}
-20 \\
-10
\end{bmatrix} = \begin{bmatrix}
-9 \\
12
\end{bmatrix}
$$

After several iterations, we will reach the solutions $$x \approx 2.5384$$ and $$y \approx 4.3077$$. However, for different initial guesses, we might get different sets of solutions. For instance, for $$X\_0 = \begin{bmatrix}  2 \  1  \end{bmatrix}$$, we will actually get $$x = 3$$ and $$y = 4$$.

Currently, you can use this method only to solve systems of polynomial equations, via the `solve_poly_system()` method.

**Parameters**

* `equations`- a collection of equations(type str) or polynomials (type Poly)
* `initial_vals` - a dictionary that represents the initial guesses. For instance: `{'x':2, 'y':1}` will be interpreted as  $$X\_0 = \begin{bmatrix}  2 \  1  \end{bmatrix}$$, 
* `epsilon` Determines the negligible difference. Default is $$0.00001$$
* `nmax` - The maximum number of iterations. Default is $$10000$$

For instance, this is how we would solve the aforementioned example:

```python
# Solving systems of polynomial equations via KiwiCalc
solutions = solve_poly_system(["x^2 + y^2 = 25", "2x + 3y = 18"], {'x': 2, 'y': 1})
print(solutions)
# output: '{'x': 3.0000001628514434, 'y': 3.999999891432371}'
                    
```