Μαθηματική Μοντελοποίηση στον υπολογιστή
Σχεδιάστε καμπύλες και λύστε εξισώσεις χρησιμοποιώντας
το πακέτο matlab.
(Δείτε το Documentation.)
Παράδειγμα 1.
Μάζα m=1, κρέμεται από ράβρο μήκους l=1
(επίσης, η σταθερά βαρύτητας θεωρείται g=1).
Το δυναμικό αυτού του συστήματος είναι V = -cosθ
(θα το γράφουμε cos(x) γιά απλότητα στον συμβολισμό).
(α) Σχεδιάστε το δυναμικό.
Μπορούμε να χρησιμοποιήσουμε τις εντολές
x=-pi:0.02:pi; --> εκλέγουμε μία σειρά σημείων στον x-άξονα
y=-cos(x); --> υπολογίζουμε τα αντίστοιχα y(x)
plot(x,y) --> σχεδιάζουμε το γράφημα
axis([-pi pi -1 1]) --> ορίζουμε το εύρος των αξόνων
xlabel('angle') --> όνομα του άξονα x
ylabel('potential') --> όνομα του άξονα y
(β)
Έστω ότι η ενέργεια του σώματος είναι E0
(π.χ., E0=-0.9, E0=0.0).
Βρείτε το εύρος της κίνησης που θα κάνει το εκκρεμές.
Το εκκρεμές θα κινηθεί στην περιοχή όπου
E0-V(x) >= 0.
Βρίσκουμε τα δύο σημεία x1, x2 στα οποία
ισχύει η ισότητα E0-V(x) = 0
και έχουμε ότι το σώμα κινείται στην περιοχή
x1 < x < x2.
Γιά να βρούμε την λύση της εξίσωσης στο matlab γράφουμε:
x1 = fzero('-0.9+cos(x)',-0.1) --> λύση την εξ. -0.9+cos(x)=0, με αρχική εκτίμηση x=0.1
x2 = fzero('-0.9+cos(x)',0.1) --> βρίσκει την 2η λύση
(γ) Βρείτε και σχεδιάστε την λύση x=x(t).
Ισχύει η εξίσωση
1/2 (dx/dt)2 - cos(x) = E0.
(δ) Βρείτε την περίοδο ταλάντωσης T.
Συγκρίνετε με την περίοδο ταλάντωσης που βρίσκουμε με την αρμονική προσέγγιση
(T=2π).
Η περίοδος δίνεται από ένα ολοκλήρωμα:
T/2 = int{ dx/[2 Sqrt(E0+cos(x))] } (όρια ολοκήρωσης x=x1, x=x2).
Στο matlab αυτό μπορεί να
υπολογισθεί ως εξής:
f= inline('1./sqrt(2*(-0.9+cos(x)))') --> ορισμός της ολοκληρωτέας συνάρτησης
period = 2*quad(f,x1,x2) --> αριθμητικός υπολογισμός ολοκληρώματος
Επαναλάβετε γιά Ε0=0.0 και γιά άλλες τιμές.
Παράδειγμα 2.
Θεωρήστε δύο δίνες με φορτιά Q1, Q2 > 0.
(a) Λύστε αριθμητικά τις εξισώσεις κίνησης.
(β) Σχεδιάστε τις τροχιές των δινών γιά (μη-τετριμένη)
αρχική συνθήκη της επιλογής σας.
(a) Γιά λύση διαφορικών εξισώσεων και σχεδίαση της λύσης με matlab δείτε
εδώ.
Σημειώστε ότι πρέπει πρώτα να γράψουμε ένα αρχείο func.m
το οποίο περιέχει εξισώσεις κίνησης (δηλ. τις χρονικές παραγώγους
των (x1,y1), (x2,y2)).
Γιά να βρει το matlab την συνάρτησή μας πρέπει να πούμε σε ποιό
directory βρίσκεται. Π.χ., αν βρίσκεται στο newpath:
path(path,'newpath')
Περιγράφουμε ένα σύντομο πρόγραμμα:
1. Πρώτα πρέπει να γράψουμε ένα αρχείο
vortices.m
το οποίο περιέχει εξισώσεις κίνησης (δηλ. τις χρονικές παραγώγους
των (x1,y1), (x2,y2)).
Γιά να βρει το matlab την συνάρτησή μας πρέπει να πούμε σε ποιό
directory βρίσκεται. Π.χ., αν βρίσκεται στο newpath:
path(path,'newpath')
2. Μετά δίνουμε τις εντολές
t0=0;
tfinal=20;
y0 = [-1 0 1 0]';
[t,y] = ode45('vortices',[t0 tfinal],y0);
3. Γιά να κάνουμε γραφική παράσταση:
plot(y(:,1),y(:,2),'-',y(:,3),y(:,4),'-');
Πληροφορίες γιά ορισμένες εντολές που χρησιμοποιούμε:
function
path