Εύρεση ριζών

6. Εύρεση ριζών#

6.1. Διατύπωση του προβλήματος#

Ρίζες μιας συνάρτηση $f (x)$ ονομάζονται οι τιμές του $x$ που μηδενίζουν την συνάρτηση:

f (x) = 0

Μια απλή περίπτωση είναι το τριώνυμο $f (x) = a x^{2} + b x + c$ , που έχει την γνωστή αναλυτική λύση:

x_{1, 2} = \frac{- b \pm \sqrt{b^{2} - 4 a c}}{2 a}

6.2. Επαναληπτικές μέθοδοι#

Στην πράξη συχνά οι εξισώσεις είναι μη γραμμικές και δεν υπάρχει διαθέσιμη αναλυτική λύση για τις ρίζες. Οι αριθμητικές μέθοδοι που έχουν αναπτυχθεί για αυτές τις περιπτώσεις είναι επαναληπτικές, δηλαδή ξεκινάνε από μια εκτίμηση για την ρίζα και σε κάθε επανάληψη την βελτιώνουν προσεγγίζοντας σταδιακά την πραγματική τιμή.

Η επαναληπτική διαδικασία σταματά, όταν ικανοποιηθεί το κριτήριο ακρίβειας. To κριτήριο ακρίβειας μπορεί να αφορά την απόλυτη ή σχετική απόκλιση μεταξύ δύο διαδοχικών εκτιμήσεων της ρίζας ή το υπόλοιπο της συνάρτησης $f (x)$ . Για να τερματιστούν οι επαναλήψεις πρέπει το επιλεγμένο μέγεθος να γίνει μικρότερο από την ανοχή $T O L$ . Ο ορισμός του κριτηρίου ακρίβειας και η τιμή της ανοχής μπορεί να διαφέρουν κατά περίπτωση.

6.3. Μέθοδος διχοτόμησης#

Σύμφωνα με το θεώρημα ενδιάμεσης τιμής, μια συνεχής συνάρτηση $f (x)$ με πεδίο ορισμού που περιλαμβάνει το διάστημα $[x_{a}, x_{b}]$ παίρνει οποιαδήποτε τιμή ανάμεσα στο $f (x_{a})$ και στο $f (x_{b})$ . Αυτό οδηγεί σε ένα χρήσιμο πόρισμα, γνωστό κι ως θεώρημα Boltzano:

Αν τα $f (x_{a})$ και στο $f (x_{b})$ είναι μη μηδενικά και έχουν διαφορετικό πρόσημο:

f (x_{a}) f (x_{b}) < 0

τότε πρέπει να υπάρχει τουλάχιστον ένα σημείο $x_{m}$ ανάμεσά τους ( $x_{a} < x_{m} < x_{b}$ ) για το οποίο ισχύει $f (x_{m}) = 0$ .

Η μέθοδος της διχοτόμησης ξεκινά με ένα διάστημα $[x_{a}, x_{b}]$ στον οποίο είναι γνωστό ότι υπάρχει ρίζα και περιορίζει σταδιακά το διάστημα για να βρεθεί η ρίζα.

Με δεδομένα τα $x_{a}$ , $x_{b}$ και $f (x_{a}) f (x_{b}) < 0$ :

Υπολόγισε το $x_{m} = \frac{x_{a} + x_{b}}{2}$
Υπολόγισε το $f (x_{m})$
Αν $f (x_{m})$ και $f (x_{a})$ έχουν το ίδιο πρόσημο (μπορεί να ελεγχθεί με διάφορους τρόπους, π.χ. $f (x_{m}) f (x_{a}) > 0$ )
- τότε $x_{a} = x_{m}$ και $f (x_{a}) = f (x_{m})$
- αλλιώς $x_{b} = x_{m}$ και $f (x_{b}) = f (x_{m})$
Αν έχει επιτευχθεί το κριτήριο σύγκλισης,
- τότε επίστρεψε την τιμή $x_{m}$ ,
- αλλιώς επανάλαβε τα βήματα με το ενημερωμένο $x_{a}$ ή $x_{b}$ .

Exercise 6.1

Να βρεθεί μία ρίζα της συνάρτησης

f (x) = e^{x} + x - 5

στο διάστημα $(0, 10)$ με την μέθοδο της διχοτόμησης και όριο ανοχής $10^{- 5}$ στην απόλυτη διαφορά μεταξύ δύο διαδοχικών εκτιμήσεων της ρίζας.

import numpy as np

def f(x: float) -> float:
    return np.exp(x)+x-5


xa = 0
xb = 10

fa = f(xa)
fb = f(xb)

print(fa*fb < 0)  # Check

xm = np.Inf
TOL = 1e-5
print("i\tx\t\tf(x)")

for i in range(100):    # 100 is max number of iterations
    xm_old = xm
    xm = 0.5*(xa+xb)
    fm = f(xm)

    if fa*fm > 0:
        xa = xm
        fa = fm
    else:
        xb = xm
        fb = fm

    print(f"{i}\t{xm:0.7f}\t{fm:0.7f}")

    e = abs(xm-xm_old)  # absolute error
    if e < TOL:
        break

True
i	x		f(x)
0	5.0000000	148.4131591
1	2.5000000	9.6824940
2	1.2500000	-0.2596570
3	1.8750000	3.3958191
4	1.5625000	1.3332332
5	1.4062500	0.4868743
6	1.3281250	0.1020856
7	1.2890625	-0.0815551
8	1.3085938	0.0095593
9	1.2988281	-0.0361726
10	1.3037109	-0.0133506
11	1.3061523	-0.0019066
12	1.3073730	0.0038236
13	1.3067627	0.0009578
14	1.3064575	-0.0004746
15	1.3066101	0.0002416
16	1.3065338	-0.0001165
17	1.3065720	0.0000625
18	1.3065529	-0.0000270
19	1.3065624	0.0000178

Για τον χαρακτηρισμό της σύγκλισης χρησιμοποιείται ο μαθηματικός ορισμός:

lim_{n \to \infty} \frac{| x_{n + 1} - x^{*} |}{| x_{n} - x^{*} |^{q}} = μ

Η μέθοδος της διχοτόμησης συγκλίνει στο $x^{*}$ με τάξη $q = 1$ (γραμμική σύγκλιση) και ταχύτητα $μ = 1 / 2$ , επομένως το σφάλμα προσεγγιστικά υποδιπλασιάζεται σε κάθε βήμα.

Η μονάδα κώδικα SciPy παρέχει μια έτοιμη υλοποίηση της μεθόδου (scipy.optimize.bisect).

6.4. Μέθοδος Newton-Raphson#

Η μέθοδος Newton-Raphson ή σκέτο μέθοδος Newton βασίζεται στην προσεγγιστική σχέση:

x_{n + 1} = x_{n} - \frac{f (x_{n})}{f^{'} (x_{n})}

όπου $x_{n + 1}$ η ρίζα που προκύπτει ξεκινώντας από μια εκτίμηση $x_{n}$ και προσεγγίζοντας την συνάρτηση $f (x)$ με ευθεία. Η συνάρτηση $f (x)$ πρέπει να είναι συνεχής. Η παραπάνω σχέση μπορεί να αποδειχθεί με αποκοπή των παραγώγων δεύτερης και μεγαλύτερης τάξης στην σειρά Taylor.

Η διαδικασία εφαρμόζεται ξεκινώντας από μια αρκετά καλή εκτίμηση και επαναλαμβάνοντας την ίδια διαδικασία μέχρι να ικανοποιηθεί το κριτήριο ακρίβειας. Το πληθος των σωστών ψηφίων διπλασιάζεται περίπου σε κάθε επανάληψη και η τάξη της σύγκλισης είναι τετραγωνική ( $q = 2$ ).

Πλεονέκτημα της μεθόδου είναι η ταχύτητα σε σχέση με την μέθοδο διχοτόμησης και μειονεκτήματα της μεθόδου είναι ότι:

πρέπει να υπολογιστεί η παράγωγος, που δεν είναι πάντα εύκολο,
η μέθοδος δεν συγκλίνει σε κάποιες περιπτώσεις, όπως
- ασυνέχεια της παραγώγου
- μέγιστα και ελάχιστα της συνάρτησης ( $f^{'} (x) = 0$ )
- κακή αρχική εκτίμηση (μπορεί να προσεγγίσει άλλη ρίζα)
- πολλαπλές ρίζες

Η μέθοδος Newton-Raphson έχει υλοποιηθεί στην συνάρτηση scipy.optimize.newton.

Exercise 6.2

Να επιλυθεί η Άσκηση 6.1 με χρήση της συνάρτησης scipy.optimize.newton.

Για την λύση της άσκησης πρέπει να υπολογιστεί η πρώτη παράγωγος:

\begin{array}{r} \begin{aligned} f (x) & = & e^{x} + x - 5 \\ f^{'} (x) & = & e^{x} + 1 \end{aligned} \end{array}

και να δηλωθούν οι δύο συναρτήσεις:

import numpy as np
from scipy.optimize import newton

def f(x: float) -> float:
    return np.exp(x)+x-5

def fp(x: float) -> float:
    return np.exp(x)+1

x = newton(f, x0=0, fprime=fp, tol=1e-5)
print(x, f(x))

1.3065586410393504 8.881784197001252e-16

6.5. Μέθοδος της τέμνουσας#

Για την εφαρμογή της μεθόδου Newton-Raphson είναι απαραίτητος ο υπολογισμός της πρώτης παραγώγου $f^{'} (x)$ . Όταν αυτό είναι δύσκολο ή ανέφικτο, η παράγωγος μπορεί να εκτιμηθεί εφαρμόζοντας προς τα πίσω πεπερασμένες διαφορές:

f^{'} (x_{n}) \approx \frac{f (x_{n}) - f (x_{n - 1})}{x_{n} - x_{n - 1}}

Η μέθοδος αυτή ονομάζεται μέθοδος της τέμνουσας. Σε σχέση με την κλασσική Newton-Raphson, η μέθοδος της τέμνουσας είναι πιο οικονομική από άποψη πράξεων ανά επανάληψη καθώς όλοι οι όροι είναι ήδη γνωστοί. Από την άλλη, η προσεγγιστική εκτίμηση της παραγώγου μπορεί να οδηγήσει σε περισσότερες επαναλήψεις (υποτετράγωνη τάξη σύγκλισης, $1 < q < 2$ ). Σημειώνεται ότι απαιτούνται δύο αρχικές εκτιμήσεις ( $x_{0}$ και $x_{1}$ ) για την εκκίνηση της μεθόδου.

Η μέθοδος της τέμνουσας θεωρείται παραλλαγή της μεθόδου Newton-Raphson και γιαυτό εφαρμόζεται με την ίδια συνάρτηση (scipy.optimize.newton) παραλείποντας απλώς το προαιρετικό όρισμα για την πρώτη παράγωγο και εισάγοντας ένα όρισμα $x_{1}$ για την δεύτερη αρχική εκτίμηση. Αν το $x_{1}$ παραληφθεί, η συνάρτηση κάνει αυτόματα μια εκτίμηση διαφοροποιώντας ελαφρώς το $x_{0}$ .

import numpy as np
from scipy.optimize import newton

def f(x: float) -> float:
    return np.exp(x)+x-5

x = newton(f, x0=0, x1=1e-4, tol=1e-5)
print(x, f(x))

1.3065586415194725 2.253425890330618e-09

6.6. Εύρεση πολλαπλών ριζών#

Για την εύρεση πολλαπλών ριζών μπορεί να χρησιμοποιηθεί η συνάρτηση scipy.optimize.fsolve, η οποία βασίζεται σε μια παραλλαγή της κατάβασης βαθμίδας (gradient descent) που έχει υλοποιηθεί στην βιβλιοθήκη MINPACK. Όπως φαίνεται κι από το όνομα της επιλεγμένης μονάδα κώδικα (scipy.optimize), υπάρχει επικάλυψη μεταξύ των τεχνικών που εφαρμόζονται στην εύρεση ριζών έχουν και στην βελτιστοποίηση. Ο ορισμός της αρχικής εκτίμησης είναι πολύ σημαντικός, καθώς επηρεάζει σε ποια ρίζα θα συγκλίνει η μέθοδος όσο και το πλήθος των επαναλήψεων.

Όταν η εξίσωση είναι πολυωνυμικής μορφής ενδείκνυται η χρήση της εξειδικευμένης συνάρτησης numpy.roots, η οποία υπολογίζει και τις μιγαδικές ρίζες.

Exercise 6.3

Να βρεθούν οι ρίζες του παρακάτω πολυωνύμου με χρήση της scipy.optimize.fsolve και της numpy.roots.

f (x) = 5 \cdot 10^{- 2} x^{4} - 0.2 x^{3} - 6 x - 50

Στην πρώτη περίπτωση να χρησιμοποιηθούν οι αρχικές εκτιμήσεις $- 100, - 10, 0, 10, 100$ .

import numpy as np
from scipy.optimize import fsolve


def f(x: float) -> float:
    return 5e-2*x**4-0.2*x**3-6*x-50

x = fsolve(f, [-100,-10,0,10,100])
print(x)

p=[5e-2,-0.2,0,-6.,-50.]
x=np.roots(p)
print(x)

[-4.01403447 -4.01403447  7.92266417  7.92266417  7.92266417]
[ 7.92266417+0.j          0.04568515+5.60737258j  0.04568515-5.60737258j
 -4.01403447+0.j        ]