Παρασκευή, 12 Ιουνίου 2015

Μέθοδος Runge-Kutta 4ης τάξης

Σε αυτή την ανάρτηση θα μελετήσουμε την μέθοδο Runge-Kutta 4ης τάξης που είναι μία αριθμητική μέθοδος επίλυσης διαφορικών εξισώσεων. Η διαφορά αριθμητικών και αναλυτικών μεθόδων επίλυσης διαφορικών εξισώσεων είναι η εξής: χρησιμοποιώντας αναλυτικές μεθόδους επίλυσης διαφορικών εξισώσεων προσδιορίζουμε απευθείας την συνάρτηση που αναζητούμε ενώ με την βοήθεια αριθμητικών μεθόδων επίλυσης διαφορικών εξισώσεων, όπως η Runge-Kutta, υπολογίζουμε σημεία που ανήκουν στην συνάρτηση που αναζητούμε. Έτσι με την μέθοδο Runge-Kutta προσπαθούμε να προσδιορίζουμε τα σημεία που ανήκουν στην λύση της διαφορικής εξίσωσης μας.

Το πρόβλημα που καλούμαστε να λύσουμε είναι ο προσδιορισμός, μέσα σε ένα διάστημα [a,b], της συνάρτησης y(x) για την οποία ισχύει η εξής διαφορική εξίσωση:
όπου f μία συνάρτηση των x και y. Για να εκτελέσουμε την μέθοδο Runge-Kutta πρέπει να γνωρίζουμε μία αρχική συνθήκη

Αν, λοιπόν, γνωρίζουμε την αρχική συνθήκη μπορούμε να υπολογίσουμε τα y που αντιστοιχούν στα σημεία a+h, a+2h+, a+3h,..., b-2h, b-h, b. Αυτό επιτυγχάνεται με την ακόλουθη διαδικασία:

Όπως είπαμε το y που αντιστοιχεί σε x=a είναι το y0. Έτσι το y που αντιστοιχεί σε x1=a+h (το συμβολίζουμε με y1) δίνεται από τον τύπο
όπου

Τώρα που γνωρίζουμε το y1 μπορούμε να υπολογίσουμε το y που αντιστοιχεί σε x2=a+2h=x1+h (το συμβολίζουμε με y2). Έτσι αυτό δίνεται από τον τύπο
όπου

Οπότε γενικεύοντας την διαδικασία αυτή, αν γνωρίζουμε το yn που αντιστοιχεί στο σημείο xn, μπορούμε να υπολογίσουμε το yn+1 που αντιστοιχεί στο σημείο xn+1=xn+h, με τον τύπο
όπου


Αυτή η διαδικασία συνεχίζεται μέχρι να υπολογιστούν τα y για κάθε x που ανήκει στο διάστημα [a,b].

Επειδή αυτή η διαδικασία απαιτεί πολλές πράξεις μπορούμε να γράψουμε έναν κώδικα ο οποίος να εκτελεί αυτή την διαδικασία. Αυτός ο κώδικας είναι ο

Private Sub Command1_Click()

Dim x_in As Double, y_in As Double, h As Double, b as Double

'εισαγωγή των a, y0, h και b από την φόρμα
x_in = Text1.Text
y_in = Text2.Text
h = Text3.Text
b = Text4.Text

'άνοιγμα του αρχείου results.txt για εξαγωγή των δεδομένων
Open "results.txt" For Output As #1

'όσο τα x είναι μικρότερα από το b ή και ίσα με αυτό
While x_in < = b
'υπολογισμός των k1, k2, k3 και k4
k1 = h * f(x_in, y_in)
k2 = h * f(x_in + h / 2, y_in + k1 / 2)
k3 = h * f(x_in + h / 2, y_in + k2 / 2)
k4 = h * f(x_in + h, y_in + k3)
'υπολογισμός του επόμενου y
y_next = y_in + (1 / 6) * (k1 + 2 * k2 + 2 * k3 + k4)
'εκτύπωση του σημείου στο αρχείο
Print #1, x_in, y_in
'προχωράμε στο επόμενο σημείο
y_in = y_next
x_in = x_in + h
Wend

Close

End Sub
'ορισμός της συνάρτησης f(x,y)
Public Function f(x As Double, y As Double) As Double
f = (y + 1) / x
End Function

Στην εικόνα 1 φαίνεται η φόρμα από την οποία ο χρήστης επιλέγει τα a, y0, h και b. Αφού ο χρήστης εισάγει αυτά τα στοιχεία, εκτελεί τον κώδικα και δημιουργείται ένα αρχείο με τα σημεία που προκύπτουν από την μέθοδο Runge-Kutta.

εικόνα 1, η φόρμα εισαγωγής δεδομένων στον κώδικα

Ο παραπάνω κώδικας είναι γραμμένος για να λύνει την διαφορική εξίσωση
και αυτό φαίνεται από το γεγονός πως στον κώδικα ισχύει f(x,y)=(y+1)/x. Η αναλυτική λύση αυτής της διαφορικής εξίσωσης (1) είναι η
όπου c μία σταθερά που υπολογίζεται από την αρχική συνθήκη. Αν η αρχική συνθήκη μας είναι η
προκύπτει πως c=2, επομένως η λύση της (1) είναι η

Ας προσπαθήσουμε να λύσουμε την διαφορική εξίσωση (1) με την μέθοδο Runge-Kutta 4ης τάξης. Έτσι στην φόρμα του κώδικα που παρουσιάστηκε προηγουμένως εισάγουμε a=1, y0=1, h=0.1 και b=5. Οι αριθμητικές τιμές των a και y0 προκύπτουν από την αρχική συνθήκη y(1)=1, ενώ οι τιμές του βήματος και του b σχετίζονται με το πόσο πυκνά θέλουμε να υπολογίζουμε σημεία και με το ποια περιοχή τιμών μας ενδιαφέρει. Έτσι εκτελούμε τον κώδικα και προκύπτουν τα εξής σημεία:

1             1 
 1.1           1.2 
 1.2           1.4 
 1.3           1.6 
 1.4           1.8 
 1.5           2 
 1.6           2.2 
 1.7           2.4 
   .               .
   .               .
   .               .
 4.4           7.8 
 4.5           8 
 4.6           8.2 
 4.7           8.4 
 4.8           8.6 
 4.9           8.8 
 5             9 

Παρατηρούμε πως κάθε σημείο επαληθεύει την εξίσωση (2) οπότε ο κώδικας μας έλυσε επιτυχώς την διαφορική εξίσωση.

Ας δοκιμάσουμε να λύσουμε την διαφορική εξίσωση
με αρχική συνθήκη την
Επιλέγουμε ως βήμα της διαδικασίας το 0.1 και μας ενδιαφέρει η περιοχή [0,10]. Έτσι εισάγουμε στον κώδικα τα εξής στοιχεία: a=0, y0=2, h=0.1 και b=10 και ορίζουμε ως συνάρτηση f την
Εκτελώντας τον κώδικα παίρνουμε τα εξής σημεία:

 0             2 
 0.1           2.01047260813742 
 0.2           2.04352648331579 
 0.3           2.10095245309671 
 0.4           2.18357556680841 
 0.5           2.29128861639933 
 0.6           2.42322200285055 
 0.7           2.57798490826328 
 0.8           2.75390765634508 
 0.9           2.94923861638538 
 1              3.16227912957838 
 1.1           3.39146134055386 
   .               .
   .               .
   .               .
 8.7           40.2959801383064 
 8.8           40.9475764278526 
 8.9           41.6026202939515 
 9.0           42.2610933996082 
 9.1           42.9229777081922 
 9.2           43.5882554749024 
 9.3           44.2569092385878 
 9.4           44.9289218139018 
 9.5           45.6042762837725 
 9.6           46.2829559921735 
 9.7           46.9649445371775 
 9.8           47.6502257642785 
 9.9           48.33878375997 
 10.0         49.0306028455642 

Συνεπώς αν πλοτάρουμε αυτά τα σημεία προκύπτει το γράφημα 1.

γράφημα 1

Η αναλυτική λύση της διαφορικής εξίσωσης (3), για την αρχική συνθήκη y(0)=2, είναι η
Η γραφική παράσταση της συνάρτησης αυτής στο διάστημα [0,10] φαίνεται στο γράφημα 2.

γράφημα 2
Παρατηρούμε πως τα δύο γραφήματα είναι πανομοιότυπα. Αυτό σημαίνει πως τα σημεία που προέκυψαν από την μέθοδο Rnge-Kutta είναι πράγματι σημεία που επαληθευθούν την συνάρτηση (4).

Πηγή: Αριθμητική ανάλυση με χρήση H/Y, Ε. Σιδηρόπουλος, Χ. Φωτιάδης

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

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