Τετάρτη, 28 Ιανουαρίου 2015

Η μέθοδος Newton-Raphson (αριθμητική επίλυση εξισώσεων μέρος 1)

Θα εξετάσουμε δύο μεθόδους επίλυσης μη γραμμικών εξισώσεων. Ήδη στην ανάρτηση Επίλυση μη γραμμικών εξισώσεων με την μέθοδο των συναρτησιακών επαναλήψεων έχουμε εξετάσει μία τέτοια μέθοδο. Οι δύο μέθοδοι που θα παρουσιάσουμε σήμερα είναι η μέθοδος Newton-Raphson και η μέθοδος της διχοτόμησης.

Θα ξεκινήσουμε με την μέθοδο Newton-Raphson. Έστω ότι θέλουμε να επιλύσουμε την εξίσωση 

Για να το πετύχουμε αυτό θα χρησιμοποιήσουμε μία επαναληπτική διαδικασία στην οποία επιλέγουμε μία αρχική τιμή κοντά στην ρίζα που θέλουμε να βρούμε και την συμβολίζουμε με x0. Έπειτα ακολουθούμε τον τύπο
Δηλαδή θέτοντας i=0 στον παραπάνω τύπο και κάνοντας τις πράξεις (αφού έχουμε επιλέξει το x0) υπολογίζουμε το x1. Έπειτα θέτοντας i=1 στον ίδιο τύπο γνωρίζοντας πλέον το x1 βρίσκουμε το x2 κτλ. Συνεχίζουμε αυτή την διαδικασία μέχρι να έχουμε σύγκλιση, δηλαδή αν μετά από i επαναλήψεις για το xi ισχύει προσεγγιστικά
(όπου το ε ένας μικρός αριθμός) θεωρούμε το xi ως ρίζα. Αν δεν έχουμε σύγκλιση μετά από έναν αριθμό n επαναλήψεων σταματάμε την διαδικασία.

Ένας κώδικας που υλοποιεί την παραπάνω διαδικασία (αλγόριθμο) είναι ο ακόλουθος:

Private Sub Command1_Click()

'δήλωση των μεταβλητών
Dim x0 As Double, x1 As Double, n As Integer, i As Integer, flag As Boolean

'εισαγωγή από την φόρμα της αρχικής τιμής
x0 = Text1.Text

'επιλογή του μέγιστου αριθμού επαναλήψεων
n = Text2.Text

'θέτουμε το i=1
i = 1
'το flag είναι η μεταβλητή που σηματοδοτεί αν έχουμε βρει ρίζα.
'Όταν είναι false σημαίνει ότι δεν έχει βρεθεί η ρίζα
flag = False

'Ανοίγουμε ένα αρχείο για εξαγωγή των αποτελεσμάτων της κάθε επανάληψης
Open "results.txt" For Output As #1

'τυπώνουμε στο αρχείο το x0
Print #1, "0", x0

'θα γίνονται επαναλήψεις μέχρι να βρούμε ρίζα (flag=true) ή μέχρι να υπερβούμε το n
Do Until i > n Or flag = True

'υπολογίζουμε το xi+1
x1 = x0 - f(x0) / df(x0)

'θέτουμε ως xi το xi+1
x0 = x1

'τυπώνουμε στο αρχείο το i και το xi
Print #1, i, x1

'αν βρούμε ρίζα...
If Abs(f(x1)) < 0.0001 Then
'... θέτουμε το flag=true (ε=0.0001)
flag = True
End If

'αύξησε το i κατά 1
i = i + 1
Loop

'δώσε στην φόρμα το αποτέλεσμα
Text3.Text = Format(x1, "0.0000")

'κλείσε το αρχείο εξαγωγής δεδομένων
Close

End Sub

'ορισμός της f(x) και της παραγώγου της
Function f(x As Double) As Double
f = x ^ 2 - 5 * x + 6
End Function
Function df(x As Double) As Double
df = 2 * x - 5
End Function

Ο παραπάνω κώδικας έχει δημιουργηθεί για να επιλύει την εξίσωση
Γνωρίζουμε ότι οι ρίζες αυτής της εξίσωσης είναι το 2 και το 3. Με αυτό τον τρόπο θα μπορέσουμε να ελέγξουμε τον κώδικα μας. Στην εικόνα 1 φαίνεται η φόρμα εισαγωγής δεδομένων στον κώδικα.

εικόνα 1, η φόρμα εισαγωγής δεδομένων στον κώδικα.
Έπειτα θα εκτελέσουμε τον κώδικα χρησιμοποιώντας ως αρχική τιμή μία τιμή κοντά στις ρίζες. Ας επιλέξουμε ως αρχική τιμή το 0 και ας βάλουμε ως όριο τις 20 επαναλήψεις (n=20). Έτσι τρέχοντας τον κώδικα παίρνουμε ως ρίζα το 2 (εικόνα 2).

εικόνα 2, η φόρμα μετά την πρώτη εκτέλεση του κώδικα.

Για να δούμε την πορεία της επαναληπτικής διαδικασίας ανατρέχουμε στο αρχείο που δημιούργησε ο κώδικας. Αυτό είναι το παρακάτω:

i              xi
0             0 
1             1.2 
2             1.75384615384615 
3             1.9593973037272 
4             1.99847523980551 
5             1.9999976821746 

Βλέπουμε ότι έχουμε σύγκλιση πολύ γρήγορα (μόνο 5 επαναλήψεις). Όμως τι γίνεται με την άλλη ρίζα της εξίσωσης; Πολύ απλά το ποια ρίζα θα μας δώσει η διαδικασία εξαρτάται από την αρχική τιμή που θα επιλέξουμε. Αν επιλέξουμε ως αρχική τιμή το 5 και διατηρώντας το ίδιο όριο επαναλήψεων έχουμε ως ρίζα το 3 (εικόνα 3).

εικόνα 3, η φόρμα μετά την δεύτερη εκτέλεση του κώδικα.

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

i              xi
0             5 
1             3.8 
2             3.24615384615385 
3             3.0406026962728 
4             3.00152476019449 
5             3.0000023178254 

Επίσης έχουμε σύγκλιση μετά από 5 επαναλήψεις. Συνεπώς η αρχική τιμή καθορίζει σε ποια ρίζα θα συγκλίνει ο κώδικας. Ας δούμε τώρα τι θα γίνει αν βάλουμε μία αρχική τιμή που δεν είναι κοντά σε κάποια ρίζα πχ το 1000 και διατηρώντας πάλι το όριο των 20 επαναλήψεων. Πάλι παίρνουμε ως ρίζα το 3. Βλέποντας το αρχείο που δημιούργησε ο κώδικας έχουμε:

i              xi
0             1000
1             501.250125313283
2             251.875313283145
3             127.188157894076
4             64.8450814480139
5             33.674545693719
6             18.0912825283871
7             10.3036585648725
8             6.41784741021065
9             4.4908289803527
10            3.55820240477555
11            3.14722605207238
12            3.01674493041451
13            3.0002713066719
14            3.00000007356739

Παρατηρούμε ότι ο αλγόριθμος "δυσκολεύτηκε" αυτή την φορά καθώς έκανε 14 επαναλήψεις αλλά πάλι μας έδωσε αποτέλεσμα.

Η συνέχεια στην ανάρτηση Η μέθοδος της διχοτόμησης

Δεν υπάρχουν σχόλια:

Δημοσίευση σχολίου