Σχεδίαση της συνθήκης επανάληψης
Ένα κύριο σημείο του αλγορίθμου είναι η συνθήκη επανάληψης της διαδικασίας υποδιαίρεσης του διαστήματος. Αυτό μπορεί να επιτευχθεί με δύο τρόπους.
- Στην Απόδειξη του Θεωρήματος των Ενδιαμέσων Τιμών είδαμε ότι οι δύο ακολουθίες \(a_n\) και \(b_n\) συγκλίνουν σε μία ρίζα της εξίσωσης \(f(x)=\eta\), έστω \(x_0\). Επειδή: $$a_n\leq x_0\leq b_n\ \&\ b_n-a_n=\frac{(\beta-\alpha)}{2^n}$$ τότε: $$\forall n\in\mathbb{Z}:\ 0\leq x_0 -a_n\leq \frac{(\beta-\alpha)}{2^n}$$ Αυτό σημαίνει ότι μπορούμε να προσεγγίσουμε τη ρίζα \(x_0\) από το \(a_n\) με οποιαδήποτε ακρίβεια θέλουμε αφού αυτή δηλωθεί στην αρχή. Δίνοντας έναν θετικό πραγματικό \(\varepsilon\) και ψάχνοντας τον ακέραιο \(n\) έτσι ώστε \(\displaystyle{\frac{(\beta-\alpha)}{2^n}\leq \varepsilon}\).
- Ας ξαναδούμε λίγο τον ορισμό της συνέχειας στη γλώσσα των δεκαδικών αριθμών. (\(\forall m\exists n\)) Για κάθε δεδομένο \(m\) υπάρχει ένα \(n\) έτσι ώστε για τον υπολογισμό της τιμής της \(f(x)\) με ακρίβεια της τάξης του \(10^{-m}\), να είναι αρκετό να γνωρίζουμε τον \(x\) μέχρι το \(n-\)οστό ψηφίο μετά την υποδιαστολή, δηλαδή να βρούμε τον \(x\) με \(n\) δεκαδικά ψηφία. Χρησιμοποιώντας αυτόν τον ορισμό, μπορούμε να γράψουμε το αλγόριθμο με την διχοτόμηση διαστήματος. Ας υποθέσουμε ότι θέλουμε να βρούμε το \(x\) έτσι ώστε \(f(x)=\eta\).
- Βρίσκουμε το μεγαλύτερο δυνατό \(x(1)\) για το οποίο \(f(x(1))<\eta\) και μετά το μεγαλύτερο \(x(2)\) έτσι ώστε \(f(x(2))<\eta\), κ.ο.κ
- Συνεχίζοντας έτσι τη διαδικασία κατoρθώσαμε να κατασκευάσουμε τον μεγαλύτερο \(n\) έτσι ώστε \(f(x(n))<\eta\), \(\forall n\). (Μπορεί να μην είναι το supremum του συνόλου των \(y\) έτσι ώστε \(f(y)<\eta\), αλλά με ελαφρώς διαφορετικό όμως επιχείρημα μπορούμε να έχουμε ένα τέτοιο \(x\)).
- Με ένα απλό επαγωγικό επιχείρημα έχουμε ότι:$$f(x(n))<\eta\wedge f(x(n)+10^{-n})\geq \eta$$ Αφού \(f\) είναι συνεχής, μπορούμε να επιλέξουμε \(n\) έτσι ώστε τα \(f(x(n))\) και \(f(x(n)+10^{-n})\) να απέχουν από το \(f(x)\) λιγότερο από \(10^{-m}\) τότε το \(f(x)\) απέχει το πολύ \(10^{-m}\) από τον \(\eta\).
- Αφού αυτό μπορούμε να το κάνουμε για κάθε \(m\), ο \(f(x)\) πρέπει να είναι ίσος με \(\eta\).
Εύρεση ριζών σε διάστημα, Python
Ας δούμε τον απλούστερο αλγόριθμο στην Python που κατασκεύασε μια τάξη του Λυκείου.
import mathdef fun(x): #return math.cos(x) #return 360*x*x*x*x*x-6001704/5*x*x*x*x+57005893/50*x*x*x-39668427/100*x*x+3000047/50*x-333333/100 #return x*x*x*x*x*x*x*x*x*x*x*x*x*x*x*x*x-2*(6*x-1)*(2*x-1) #return x*x*x*x*x-2*(3*x-1)*(3*x-1) #return x**59-2*(6*x-1)**2 return x*x*x-17*x-13def main(): a=float(input("Δώσε το αριστερό όριο του διαστήματος: ")) b=float(input("Δώσε το δεξί όριο του διαστήματος: ")) precision=float(input("Δώσε την επιθυμητή ακρίβεια της ρίζας: ")) x_0=(a+b)/2 while abs(fun(x_0)) >= precision: x_0=(a+b)/2 if fun(x_0)==0: break elif fun(a)*fun(x_0)<0: b=x_0 else: a=x_0 print(x_0) print("Approximation of solution is:",x_0)main()
Το πολυώνυμο στην είσοδο είναι το \(x^3-17x-13\). Το πολυώνυμο αυτό έχει 3 πραγματικές ρίζες:
$$-3.668252782,\ -0.7941698527,\ 4.462422635$$ Μπορείτε να δουλέψετε με τον αλγόριθμο στoν On line Python Compiler.
Για το \(x^{17}-2(6x-1)^2\) το οποίο όπως έχουμε δει, έχει 2 πραγματικές ρίζες πολύ κοντά,
$$0.1666666380,\ 0.1666666953,\ 1.305918236$$ Στο διάστημα \((0,0.2)\) με ανοχή λάθους της τάξης \(10^{-10}\) (στην Python γράψτε \(0.0000….n,\ n\in\mathbb{N})\) ο αλγόριθμος δίνει 1 ρίζα, ενώ το πολυώνυμο έχει 2 ρίζες πάρα πολύ κοντά η μια με την άλλη. Δείτε, όλες τις αρνητικές απαντήσεις του αλγορίθμου όταν προσπαθείστε να επιλέξτε το διάστημα για να υπολογίσει τις κοντινές ρίζες.
Ένας άλλος αλγόριθμος ικανότερος δίνει καλύτερα αποτελέσματα στο πολυώνυμο Mignotte.
import mathdef f(x): #return 360*x*x*x*x*x-6001704/5*x*x*x*x+57005893/50*x*x*x-39668427/100*x*x+3000047/50*x-333333/100 #return x**59-2*(6*x-1)**2 # #.16666666666666382426, .16666666666666950908, 1.12738466824392800662 #return math.sin(x) #return math.log(x)-x+2 ##0.01023855228144327807 η ρίζα της λογαριθμικής sto [0.1,0.3] #return math.e**x-1 #return x-x**(1/3)-2 return x**17-2*(6*x-1)**2 #.16666666666666666667, .16666666666666666667, 1.07154057527829735302def bisection(x0,x1,e): step = 1 print('\n\n*** H Μέθοδος Διχοτομίας Διαστήματος ***') condition = True while condition: x2 = (x0 + x1)/2 print('Iteration-%d, x2 = %0.6f and f(x2) = %0.6f' % (step, x2, f(x2))) if f(x0) * f(x2) < 0: x1 = x2 else: x0 = x2 step = step + 1 condition = abs(f(x2)) > e #abs(x1-x0)/2**step > e # print('\nRequired Root is : %0.8f' % x2)# Input Sectionx0 = input('Αριστερό Όριο: ')x1 = input('Δεξιό Όριο: ')e = input('Aνοχή λάθους: ')# Converting input to floatx0 = float(x0)x1 = float(x1)e = float(e)# Έλεγχος δεδομένων εισόδουif f(x0) * f(x1) > 0.0: print('To Θεώρημα Bolzano δεν είναι ικανο να εντοπισει ρίζα.') print('Δοκίμασε ξανά με νέα όρια')else: bisection(x0,x1,e)
(οι 2 παραπάνω αλγόριθμοι υλοποιήθηκαν το έτος 2022)
Αποτελεσματικότητα των αλγορίθμων
Ο έλεγχος της αποτελεσματικότητας θα γίνει με τη χρήση των ακραίων καταστάσεων σχετικά με την πυκνότητα των ριζών που δίνουν τα πολυώνυμα Mignotte. Εδώ θα δοκιμάσουμε με τις ρίζες του πολυωνύμου $$x^{59}-2(6x-1)^2$$ To πολυώνυμο έχει 3 πραγματικές ρίζες όπως μας πληροφορεί ο αλγόριθμος sturm
$$ \hbox{sturm}(x^{59}-2(6x-1)^2,x, -\infty,+\infty)$$
$$3$$ Στο Maple με στρογγυλοποίηση υπολογισμών στο 1000-οστό ψηφίο και στρογγυλοποίηση εμφάνισης 100 ψηφίων, οι τρείς πραγματικές ρίζες του πολυωνύμου είναι:
$$\begin{array}{lcl}
\rho_1&:=&0.1666666666666666666666536608810952185499…51968244214,
\\
\rho_2&:=&0.1666666666666666666666679724522381147833603…732685769,\\
\rho_3&:=&1.0715405752782973530218708513184…4612023643215342
\end{array}$$
Η συνάρτηση sturm του Maple δίνει στο διάστημα:
$$[0.1666666666666666666666653, 0.1666666666666666666666655]$$
μία ακριβώς ρίζα:
$$\hbox{sturm} (x^{59}-2(6x-1)^2,x, 0.1666666666666666666666653, 0.1666666666666666666666655)$$ $$1$$
Επίσης, στο διάστημα
$$[0.1666666666666666666666655, 0.1888]$$
δίνει μια ρίζα (τη δεύτερη ρίζα, \(\rho_2\), που βρήκαμε προηγουμένως).
$$ \hbox{sturm}(x^{59}-2(6x-1)^2,x, 0.1666666666666666666666655, 0.1888)$$ $$1$$
O αλγόριθμος 1 στη Python, με ανοχή λάθους της τάξης \(10^{-20}\), γυρνάει ρίζα του πολυωνύμου στο πρώτο διάστημα τον αριθμό: $$0.16666666666666666$$
ενώ στο δεύτερο διάστημα τον αριθμό: $$0.16666700439453125$$
Σε σύγκριση με την αποτελεσματικότητα του δευτέρου αλγορίθμου, στο πρώτο διάστημα δεν βρίσκει ρίζα και γυρνά το μήνυμα ότι:
Ο αλγόριθμος δεν είναι ικανός να εντοπίσει ρίζα
στο δεύτερο διάστημα δίνει ρίζα τον αριθμό: $$0.16666667$$
Επομένως, ο πρώτος αλγόριθμος φαίνεται ότι είναι αποτελεσματικότερος, και αυτό οφείλεται στην συνθήκη επανάληψης. Και οι δύο όμως αλγόριθμοι δεν είναι αξιόπιστοι ως προς τον εντοπισμό των ριζών. Δεν βρίσκουν ρίζες που μπορεί να θεωρηθούν ότι προσεγγίζουν τις πραγματικές ρίζες του πολυωνύμου. Αν τους συγκρίνουμε με την αποτελεσματικότητα του Maple, αυτό φαίνεται να ωφείλεται στο ότι δεν ελέγχουμε την επιθυμητή προσέγγιση στον πυρήνα του συστήματος ο οποίος εκτελεί τις αριθμητικές πράξεις. Εκεί φαίνεται να υπερισχύουν οι αλγόριθμοι των συστημάτων CAS.
Σε πολυώνυμα όπου οι ρίζες είναι διακριτές και όχι πυκνές, όπως συναντήσαμε στα πολυώνυμα Mignotte, οι αλγόριθμοι είναι αποτελεσματικοί.
Το άλλο πολυώνυμο που παρουσιάζει ενδιαφέρον λόγω του ότι έχει απομακρυσμένες μεταξύ τους ρίζες είναι το εξής:
το πολυώνυμο έχει 5 πραγματικες ρίζες, οι οποίες με προσεγγιση \(10^{-10}\) είναι οι εξής: $$0.1666666…66667,\ 0.200000…0000,\ 0.2500000….000,\ 0.3333333…333,\ 3333.330000$$ Βλέπουμε ότι μια ρίζα η \(3333.330000\) ξεχωρίζει. Για να αρχίσουμε να ζητάμε τις ρίζες αυτές με την μέθοδο της διχοτομίας διαστήματος, πρέπει να εντοπίσουμε το διάστημα μέσα στο οποίο βρίσκονται, Το διάστημα αυτό είναι το $$(-3335.3,3335.3)$$ όπως αυτό προέκυψε από το θεώρημα Cauchy. Με τον αλγόριθμο της ακολουθίας sturm στο Maple εντοπίζουμε τα διαστήματα μέσα στα οποία βρίσκονται οι ρίζες: $$\begin{array}{clc} \text{Επιλογή διαστήματος}& \text{Αλγόριθμος στο Maple}&\text{Αριθμός πραγματικών ριζών}\\ \hline (-3335.3, 3335.3) &\text{sturm(p,x,-3335.3, 3335.3)=}& 5\\ (0,3335.3)&\text{sturm(p,x,0,3335.3)=}&5\\ (0,1700)&\text{sturm(p,x,0,1700)=}&4\\ (1700,3335.3)&\text{sturm(p,x,1700,3335.3)=}&1\\ \dots&&\\ (0,0.5)&\text{sturm(p,x,0,0.5)=}&4\\ (0,0.25)&\text{sturm(p,x,0,0.25)=}&3\\(0.25,0.5)&\text{sturm(p,x,0.25,0.5)=}&1\\ (0.125,0.1875)&\text{sturm(p,x,0.125,0.1875)=}&1\\ (0.1875,0.21875)&\text{sturm(p,x,0.1875,0.21875)=}&1\\(0.2875,0.25)&\text{sturm(p,x,0.21875,0.25)=}&1\\ \end{array}$$
Θα βρουμε τώρα αριθμητικά τις ρίζες με ανοχή λάθους \(10^{-5}\) χρησιμοποιώντας τον πρώτο αλγόριθμο στην Python μέσα σε κάθε διάστημα στο οποίο υπάρχει ρίζα σύμφωνα με τον παραπάνω πίνακα. $$\begin{array}{cc} \text{ Διάστημα} & \text{Ρίζα σε δεκαδική μορφή}\\ \hline (1700,3335.3)& 3333.33\\ (0.25,0.5)& 0.3333333358168602 \\ (0.125,0.1875)&0.1666666567325592 \\ (0.1875,0.21875)&0.19999998807907104\\ (0.2875,0.25) & 0.2500000178813934 \end{array}$$ Καθόλου άσχημα!
O αλγόριθμος δίνει ρίζες συναρτήσεων σε διαστήματα, γενικά. Υπάρχουν ήδη στην είσοδο του αλγορίθμου 1 στην Python, συναρτήσεις όπως: $$\displaystyle{y=\sin(x)-\frac{x}{2},\quad y=x\log(2)-2\log(x),\quad y=(2x+1)e^{-2x},\quad y=\frac{x-\sin(x-3)}{2x+\sin(x)},\dots}$$