ΦΥΣ-36. Υπολογιστική Φυσική Ι
Διάλεξη 4

ΘΕΜΑΤΑ

• Αριθμητική Λύση Συνήθων Διαφορικών Εξισώσεων. Μέθοδοι Runge-Kutta.
• Σύγκριση μεθόδων Euler και Runge-Kutta 4ης Τάξης

Βιβλιογραφία

• Σημειώσεις διδάσκοντα: Παράγραφοι 3.3-3.4. Επανάληψη τις 3.1-3.2.

Επί πλέον Στοιχεία:

• Το τέταρτο κεφάλαιο των παλαιών σημειώσεων και η αντίστοιχη 5η Άσκηση .
• Κεφάλαια 0,1 από "Integration of Ordinary Differential Equations" από τις Numerical Recipes (στην παλιά έκδοση Fortran είναι κεφ. 16, στην καινούργια C++ κεφ. 17).
• Weisstein, Eric W. "Runge-Kutta Method." From MathWorld--A Wolfram Web Resource. mathworld.wolfram.com/Runge-KuttaMethod.html
• Runge-Kutta methods, Wikipedia
• J. H. E. Cartwright & O. Piro, The dynamics of Runge-Kutta methods, Int. J. Bifurcation and Chaos 2, 427-449, 1992.
• J. Mathews and K. Fink, κεφ. 9.5.
• Module for Runge-Kutta Method. From numerical analysis project of John Matthews.
• Adaptive Runge-Kutta-Fehlberg Method.

Εξάσκηση

• Γράψτε τον κώδικα για τη μέθοδο Runge-Kutta για δύο άγνωστες συναρτήσεις. Σκοπός είναι να λύσουμε το πρόβλημα κίνησης σωματιδίου υπό την επίδραση γνωστής δύναμης που του προσδίδει επιτάχυνση a(t,x,v). Ο κώδικας μπορεί να οργανωθεί έτσι ώστε το κυρίως πρόγραμμα να καλεί την υπορουτίνα RK(T,X1,X2,T0,TF,X10,X20,STEPS,PSIZE) όπου στην είσοδο θα της δίνει τον αρχικό και τελικό χρόνο T0, TF, τις αρχικές θέσεις X10, X20 καθώς και τον αριθμό των βημάτων STEPS. Στην έξοδο η ρουτίνα θα τοποθετεί τα αποτελέσματα στα arrays T(PSIZE), X1(PSIZE), X2(PSIZE) τα οποία και θα τυπώνονται σε αρχείο. Η ρουτίνα RK θα είναι οδηγός της κυρίας υπορουτίνας RKSTEP(t,x1,x2,dt) η οποία θα την καλεί να κάνει STEPS βήματα χρόνου dt από τη θέση (x1(t),x2(t)) -> (x1(t+dt),x2(t+dt)). Ο κύριος υπολογισμός θα γίνεται στην RK. Λύση
  Δοκιμάστε να πάρετε αποτελέσματα και να τα συγκρίνετε με την αναμενόμενη λύση. Πάρτε x₀=1, v₀=0 και δείτε πώς μεταβάλλεται η ακρίβεια των αποτελεσμάτων που παίρνετε:
  >gfortran rk.f90 -o rk;./rk
Runge-Kutta Method for 2-ODEs Integration
Enter STEPS,T0,TF,X10,X20:
10000 0 2 1 0
  εντολές που τρέχουν το πρόγραμμα και παράγουν αποτελέσματα στο αρχείο rk.dat. Ξεκινάμε το gnuplot για να δούμε τα αποτελέσματα και κάνουμε τη γραφική παράσταση της x(t). Τη συγκρίνουμε με τη γνωστή λύση x(t) = cos(ω t):
  >gnuplot
gnuplot>o2 = 100;o=sqrt(o2)
gnuplot>plot "rk.dat" u 1:2 w l,cos(o*x)
gnuplot>plot "rk.dat" u 1:(abs($2-cos(o*$1))) w l
  Η ενέργεια ανά μονάδα μάζας (1/2)(v²+ω²x²) διατηρείται. Κάνουμε τη γραφική παράσταση και τη συγκρίνουμε με την αναμενόμενη τιμή ( = 50):
  gnuplot>plot "rk.dat" u 1:(0.5*($3*$3+o2*$2*$2)) w l
gnuplot>plot "rk.dat" u 1:(abs(0.5*($3*$3+o2*$2*$2)-50)) w l
  Επαναλαμβάνουμε τον υπολογισμό για την ταχύτητα v(t) και τη συγκρίνουμε με την αναλυτική λύση v(t) = - ω sin(ω t):
  gnuplot>plot "rk.dat" u 1:3 w l,-o*sin(o*x)
gnuplot>plot "rk.dat" u 1:(abs($3+o*sin(o*$1))) w l
  Στη συνέχεια μεταβάλλετε τον αριθμό βημάτων από STEPS=100-100000 και μελετήστε τα σφάλματα σα συνάρτηση του δt= h = 2/STEPS. Υπολογίστε το μέγιστο σφάλμα στο χρονικό διάστημα [0,0.5] και στο [0,2]. Κάντε τη γραφική παράσταση σε log-log κλίμακα και από την κλίση υπολογίστε τον εκθέτη p, όπου σφάλμα ~ h^p.
  Αλλάξτε όλες τις μεταβλητές και συναρτήσεις από REAL σε REAL(8) στο πρόγραμμα rk.f90 και επαναλάβατε τα παραπάνω.
• Μεταβάλλετε τους κώδικές σας rk.f90, euler.f90 έτσι ώστε να ολοκληρώνουν τον αρμονικό ταλαντωτή με μεταβλητές διπλής ακρίβειας. Τυπώστε σε κάθε βήμα την ολική ενέργεια του ταλαντωτή, έτσι ώστε να εξετάζεται αν αυτή πράγματι διατηρήται. (Λύση: hoc_euler.f90,hoc_rk.f90).
• Υπολογίστε για διαφορετικές τιμές των βημάτων Δt την απόκλιση της αριθμητικής λύσης από τη θεωρητική, καθώς και αυτή της ενέργειας από την αρχική της τιμή. Εξετάστε αν η βάθμωση των αποκλίσεων ταυτίζεται με τη θεωρητική πρόβλεψη.
• Η παραπάνω διαδικασία είναι χρονοβόρα (δοκιμάστε τη...). Δοκιμάστε να χρησιμοποιήσετε το σενάριο φλοιού (shell script) hoc.csh για να την αυτοματοποιήσετε. Δοκιμάστε τις εντολές:
  hoc.csh -h
hoc.csh -x 0.3 -v 0 -t 6 100 500 1000 5000 10000 50000
  Δοκιμάστε να καταλάβετε πώς δουλεύει (ξεφεύγει από τα όρια του μαθήματος για την ώρα...).