שיטת אברת'

**שִׁיטַת אָבֵּרְתּ**', או **שִׁיטַת אָבֵּרְתּ' אֶרְלִיךְ** (באנגלית: **Aberth-Ehrlich Method**) היא אלגוריתם איטרטיבי למציאת שורשים מרובים (ממשיים ומרוכבים) של פונקציה פולינומית בעלת משתנה אחד באופן סימולטני. כלומר, היא מוצאת את הפתרונות למשוואה f(x) =0, כאשר f(x) פונקציה פולינומית. היא פותחה בשנת 1967 על ידי המתמטיקאי אוליבר אברת' וכן בגרסה מעט שונה על ידי המתמטיקאי לואיס ו. ארליך ב-1973. השיטה מתחילה במציאה של קירובים לשורשי הפונקציה, אשר מתכנסים אליהם יותר ויותר בכל איטרציה. לשיטה מקדם התכנסות של 3, ולפיכך יותר יעילה מאלגוריתם דורנד-קרנר, שיטה נוספת למציאת שורשים של פולינום בבת אחת, אשר מתכנסת באופן ריבועי (מקדם של 2).

## רקע
בתחומים רבים, ובייחוד ובפיזיקה ובהנדסה, נעשה שימוש במתמטיקה על מנת לתאר מצבים מורכבים. לעיתים קרובות, על מנת לפתור בעיות בתחומים אלו יש לבטא אותן כמשוואה, או כמערכת משוואות, ולמצוא את הפתרונות המתאימים.  ניתן להפוך כל משוואה למשוואה מן הצורה f(x) = 0, ולכן נהוג להתייחס לפתרונות של משוואה כשורש של פונקציה זו, כלומר, החיתוך של הפונקציה f(x)עם ציר ה x. אפוא, מציאת שורשי פונקציה זו שקולה לפתירת המשוואה. אחד מסוגי המשוואות הנפוצים ביותר אשר נתקלים בהם הוא משוואה פולינומית. פולינום הוא ביטוי אשר יכול להכיל רק חיבור, חיסור, כפל, חילוק, והעלאה בחזקה כאשר המעריך הוא מספר לא שלילי. פולינום יכול להכיל מספר משתנים, או משתנה אחד בלבד.  לדוגמה, הביטוי  3x^2 + 4y + 1 הוא פולינום, עבור המשתנים x ו-y. פתירת משוואות פולינומיות מכמה משתנים דורשת בדרך כלל יצירה של מערכת משוואות, אך שיטת אברת' מתייחסת לפתירה של משוואה פולינומית בעלת משתנה אחד בלבד. ניתן לפתור משוואות פולינומיות מסוימות באמצעות נוסחאות ישירות. למשל, על מנת לפתור משוואה ממעלה 2, ניתן להשתמש בנוסחה הריבועית המוכרת, ובשביל משוואה ממעלה 3, ניתן להשתמש בנוסחת קרדאנו. עם זאת, כפי שהוכח במשפט אבל-רופיני במאה ה-19, לא ניתן לנסח נוסחה כללית וישירה לפתירה של משוואות ממעלה 5 או יותר. מכאן נובע הצורך לפתח שיטות אחרות לפתרון משוואות אלו. על מנת לענות על קושי זה, פותחו שיטות נומריות למציאת קירובים לפתרונות המשוואה, קרי, לשורשי הפונקציה. שיטות אלו נקראות שיטות למציאת שורשים ( באנגלית : **Root Finding Methods** ) ואחת מן השיטות המוכרות והמוקדמות ביותר לכך היא שיטת ניוטון-רפסון : x_{n+1} = x_n - \frac {f(x)}{f'(x)}.  בשיטה זו, בדומה לשיטות רבות אחרות למציאת שורשים ( כולל שיטת אברת'), בוחרים x התחלתי, ומשנים אותו שוב ושוב עד שהוא מתכנס לשורש של הפונקציה. שיטות מסוג זה הן איטרטיביות, שכן מבצעים חזרות ( איטרציות ) עד שהקירוב ההתחלתי מתכנס לשורש. עם זאת, בשיטת ניוטון-רפסון, ובאלגוריתמים רבים אחרים כגון שיטת היילי ( **Halley's Method** ) ושיטת החצייה ( **Bisection Method** ) ניתן למצוא רק פתרון אחד של המשוואה. מכיוון שפעמים רבות עולה הצורך למצוא את כל פתרונות המשוואה, פותחו שיטות למציאת שורשים מרובים. בין השיטות האיטרטיביות לכך נמנות שיטת דורנד-קרנר ( Durand-Kerner Method ) ,שיטת אברת' ארליך ( **Aberth-Ehrlich Method** ) ושיטת ניוטון מהיילי. בנוסף, ניתן למצוא שורשים מרובים באמצעות מטריצות מלוות ( Companion Matrices ). בשיטות האיטרטיביות למציאת כל שורשי הפונקציה, נהוג לחשב מספר קירובים התחלתיים לפתרונות, כאשר כל קירוב יתכנס לשורש המתאים לו בפונקציה. למרבית המשוואות הפולינומיות יש גם פתרונות מרוכבים, ולפיכך בשיטות אלה נהוג ליצור קירובים התחלתיים מרוכבים. למשל, אין פתרונות ממשיים למשוואה הפולינומית הידועה x^2 + 1 = 0, אלא פתרונות מרוכבים: \pm\sqrt{-1}, שמסומנים כ \pm i.

## היסטוריה
השיטה פותחה לראשונה על ידי המתמטיקאי אוליבר אברת' בשנת 1967  כשיפור לשיטה שהציגו המתמטיקאים וורסטראס (1891), דורנד (1960) וקרנר (1966). אברת' השתמש בשדות וקטורים ובאלקטרוסטטיקה על מנת להגיע לנוסחה. גרסה דומה לנוסחתו של אברת התגלתה מאוחר יותר ב-1973 על ידי המתמטיקאי לואיס ו. ארליך, אשר נעזר בשיטת ניוטון-רפסון על מנת להגיע אליה. 

בתקציר המאמר המקורי של אברת' נכתב:

כלומר, אברת' פיתח שיטה שבה הקירובים ההתחלתיים מתכנסים מהר יותר לשורשי הפונקציה, מאשר בשיטה שהציעו דורנד וקרנר.

## תיאור
ניתן ליישם את השיטה על פונקציות פולינומיות בעלות משתנה אחד, כלומר:

p(x)=p_nx^n+p_{n-1}x^{n-1}+\cdots+p_1x+p_0

או בקיצור: \sum_{i=0} ^ {n} a_ix^{i}

לדוגמה, f(x) = 2x^4 - 3x^3 + 2x^2 - 7x + 8

ראשית, יש לחשב את הקירובים ההתחלתיים של השורשים z_1,\dots,z_n\in\mathbb C. למשל, למשוואה ממעלה רביעית יש 4 שורשים (פתרונות), ולכן למשוואה מסוג זה ניצור 4 קירובים התחלתיים, כאשר כל אחד יתכנס לשורש אחר. למשוואה ממעלה חמישית ניצור 5 קירובים, למשוואה ממעלה שישית ניצור 6, וכן הלאה. קירובים התחלתיים אלו אינם מספרים מרוכבים רנדומליים, אלא מספרים מרוכבים אשר מפוזרים באופן שווה על מעגל במישור גאוס שמרכזו ב (0,0) ורדיוסו R, כאשר  

R = \sqrt[n]{\left \vert \frac{p_0}{p_n} \right \vert}, וכאשר n הוא החזקה הגבוהה ביותר, p_0 הוא המקדם של המספר החופשי של הפונקציה, ו p_n הוא המקדם של החזקה הגבוהה ביותר בפולינום. ניתן לחשב את R גם בשיטה נוספת, כאשר R = 1 + \max\big(|p_{n-1}|, |p_{n-2}|, ...,p_0). על מנת להשתמש בשיטה האחרונה, יש לוודא שהמקדם של החזקה הגבוהה ביותר בפולינום שווה ל-1.

אם כן, ניתן לבטא את הקירובים ההתחלתיים כz_k = R \cos{\frac{2\pi k}{n}} + iR\sin{\frac{2\pi k}{n}}, \text{for k=0....n-1}, או בקיצור לפי נוסחת אוילר: Re^{\frac{i2\pi k}{n}} \text{for k=0....n-1}.

לאחר שחושבו הקירובים ההתחלתיים, יש לשנות אותם בכל איטרציה כך שיתכנסו לשורשים המתאימים. לשם כך, בכל איטרציה נחסיר מכל קירוב הפרש מסוים המתאים לו, כך שיתכנס יותר לשורש המתאים לו.

לפיכך הפרשים אלה מוגדרים כ w_1,\dots,w_n\in\mathbb C, והם מחושבים כך:

w_k=\frac{\frac{p(z_k)}{p'(z_k)}}{1-\frac{p(z_k)}{p'(z_k)}\cdot \sum_{j\ne k}\frac1{z_k-z_j}},\text{for k=1....n}

כאשר p'(z_k) היא הנגזרת של הפונקציה עבור z_k.

אם כן, הפרשים אלו יוחסרו מהקירובים, כך שהסט הבא של הקירובים לשורשי הפונקציה יהיה z_1-w_1,\dots,z_n-w_n.

חישוב ההפרשים והחסרתם מהקירובים יתבצעו שוב ושוב עד שהקירובים יתכנסו לשורשי הפונקציה.

## מימושים תכנותיים של השיטה
מימוש השיטה בשפת התכנות פייתון:
```python
# Implementation of Aberth's method in Python, written by Jonathan.

from numpy import linspace
from math import sin,cos,pi
from typing import Callable

def random_approximations(coefficients):
 """ הפעולה מחזירה רשימה של קירובים התחלתיים, בהתאם לרשימת המקדמים אשר היא מקבלת"""
    n = len(coefficients) - 1
    radius = abs(coefficients[-1] / coefficients[0]) ** (1 / n)
    return [complex(radius * cos(angle), radius * sin(angle)) for angle in linspace(0, 2 * pi, n)]

def negligible_complex(expression: complex, epsilon:float)->bool:
    """ הפעולה מחזירה האם מספר מרוכב כה קטן כך שהוא זניח וניתן להחשיבו כ-0"""
    return abs(expression.real) < epsilon and abs(expression.imag) < epsilon

def aberth_method(f_0: Callable, f_1: Callable, coefficients, epsilon=0.000001, nmax=100000):
    """ הפעולה מקבלת את הפונקציה, את נגזרתה, ואת המקדמים של הפונקציה. בנוסף מקבלת אפסילון אשר קובע את ההפרש שלפיו מחשיבים מספר כזניח ואת מספר האיטרציות המקסימלי"""
    try:
        random_guesses = random_approximations(coefficients)
        for n in range(nmax):
            offsets = []
            for k, zk in enumerate(random_guesses):
                m = f_0(zk) / f_1(zk)  # save it as m, so it won't be calculated many times
                sigma = sum(1 / (zk - zj) for j, zj in enumerate(random_guesses) if k != j and zk != zj)
                denominator = 1 - m * sigma
                offsets.append(m / denominator)
            random_guesses = [approximation - offset for approximation, offset in zip(random_guesses, offsets)]
            if all(negligible_complex(f_0(guess), epsilon) for guess in random_guesses):
                break
        return set(random_guesses)
    except ValueError:
        return set()
```
מימוש של השיטה בשפת C++ :
```c++
// Implementation of Aberth-Ehrlich Method in C++, by Jonathan.
#include <iostream>
#include <vector>
#include <complex>
#include <cmath>

using std::vector;
using std::complex;

const int PI = 3.141592653589793;

vector<complex<float>> initial_approximations(vector<float> coefficients){
    // Generating the initial approximations that will converge to the roots.
    unsigned int n = coefficients.size();
    float radius =  std::sqrt(std::abs(coefficients[n-1] /  coefficients[0])); // the radius of the circle in the complex plane
    vector<complex<float>> initials;
    float angle;
    for(size_t k=0; k<n-1 ; ++k){
        angle = 2*PI*k /n;
        initials.push_back(complex<float>{radius*std::cos(angle),radius*std::sin(angle)});
    }
    return initials;
}

vector<complex<float>> aberth_method(complex<float>(*f_0)(complex<float>),complex<float> (*f_1)(complex<float>),
                                         vector<float>coefficients,float epsilon=0.00001,int nmax=100000){
    vector<complex<float>> initials = initial_approximations(coefficients);
    unsigned int n = coefficients.size(); // the number of coefficients
    complex<float> offsets[n-1]; // array of complex numbers that represents the differences between the approximations and the solutions
    complex<float> y,y_tag,m,sum;
    complex<float> one = complex<float>{1,0};
    bool all_converged;
    for(size_t iteration = 0; iteration < nmax ; ++iteration){ // Continue with this loop until convergence
        all_converged = true; // assuming that all of the roots have converged
        for(size_t i = 0; i<n-1 && all_converged ; ++i){ // checking if the roots have converged
                y = f_0(initials[i]); 
                if(y.real() > epsilon || y.imag() > epsilon){
                    all_converged = false;}
        }
        if(all_converged){return initials;} // if the roots have converged, return them
        // otherwise, continue with the algorithm.
        for(size_t k = 0 ; k < n-1 ; ++k){ // create n offsets, and subtract each offset from the corresponding approximation.
            y = f_0(initials[k]); // 
            y_tag = f_1(initials[k]);
            m = y/y_tag;
            sum = 0;
            for(size_t j = 0; j < n-1 ; ++j){ // computing the result of the sigma in the denominator of the formula 
            // for computing the offsets
                if(j == k || initials[k] == initials[j]){ // avoid dividing by 0 in later  
                    continue;}
                sum += one/(initials[k]-initials[j]);
            }
            offsets[k] = m / (one - m * sum);
            initials[k] -= offsets[k]; // subtract the current offset from the approximation
        }
        
    }
    printf("WARNING: it's highly likely the initial values haven't converged sufficiently to the solutions\n");
    return initials;

}
```
בדרך כלל, לאלגוריתם אברת' ארליך יש סיבוכיות זמן ריצה של O(m*n^2), כאשר n הוא מספר השורשים שיש למצוא, ו m הוא מספר האיטרציות המקסימלי ללולאה החיצונית, כפי שמתואר בדוגמה שלעיל. לפיכך, ככל שהשיטה נדרשת למצוא יותר שורשים, כך זמן הריצה שלה גדל.  בתחומים מדעיים שונים עולה לעיתים הצורך למצוא את השורשים של משוואה פולינומית בעלת מספר גדול של שורשים במהירות רבה. לכן, חוקרים שונים פיתחו אלגוריתמים מקביליים ואסינכרונים יעילים , אשר מצליחים לגרום להפחתה של זמן הריצה באמצעות שימוש יעיל בו זמני במספר משאבי חומרה ( מעבדים, מעבדים גרפיים, וכו' ) ו/או מספר תהליכים אשר מריצים מספר גדול של תהליכונים. 

## ראו גם
* אנליזה נומרית
* שורש של פונקציה
* שיטת ניוטון-רפסון

## לקריאה נוספת
*
Aberth, O. (1973, January 1)  *Iteration methods for finding all zeros of a polynomial simultaneously*. American Mathematical Society. Retrieved October 26, 2021, from https://www.ams.org/journals/mcom/1973-27-122/S0025-5718-1973-0329236-7.

## קישורים חיצוניים
* סרטון הסבר על שיטת Aberth Ehrlich

## הערות שוליים

קטגוריה:אנליזה נומרית