Πρόβλημα με προσομοίωση fdtd σε matlab

[kostas123pao]

Προσπαθώ να προσομοιώσω 4ης τάξης στο χρόνο εξισώσεις αλλά δεν δουλεύει, ενώ για 2ης τάξης δούλευε. Ξέρει κανείς τι μπορεί να φταίει?

len=100;
i=zeros(1,len+2);
v=zeros(1,len+2);

dz=0.01;
c=0.1;
l=250;
dt=dz*5;


for n=1:2000
t=n*dt;


for k=3:len+1
v(k)=v(k)-(dt/(c*24*dz))*(i(k-2)-(27*i(k-1))+(27*i(k))-i(k+1));
end


v(1)=sin(2*pi*t);
v(len)=0;


for k=2:len
i(k)=i(k)-(dt/(l*24*dz))*(v(k-1)-(27*v(k))+(27*v(k+1))-v(k+2));
end


plot(v);
axis([0 len -5 5]);

frame = getframe;

end

[Technology fan]

τι πρόβλημα σου βγάζει?

[kostas123pao]

τι πρόβλημα σου βγάζει?

Koίτα το πρόγραμμα τρέχει αλλά το διαγραμμα δεν είναι σωστό και όταν τρέχει λαγγαρει και το pc (νομίζω βγάζει πολύ μεγάλες τιμές).

[Technology fan]

δεν κολλάει από τις μεγάλες τιμές αλλά από τα 3 for loop που έχεις. Αν το διάγραμμα δεν είναι σωστό αλλά το πρόγραμμά σου τρέχει τότε το πρόβλημα σου είναι "λογικής φύσης" και πρέπει να βεβαιωθείς οτι οι πράξεις που κάνεις δεν έχουν κανα λαθάκι...

 

Επίσης ποια παράμετρος σε ενδιαφέρει? εννοώ το διάγραμμα από τις τιμές του 'i' το βγάζεις?

 

 εδιτ: ακυρο, το v κοιτάμε...

[kostas123pao]

Κοίτα το v πέρνει τιμές πχ 1.0e+023 και όσο περναει η ώρα μεγαλώνει και για αυτό λαγκάρει.

Και γω αυτό σκέφτομαι ότι έχω κάνει λάθος στις εξισώσεις αλλά τις έχω ελένξει πολλές φορές και τις βρίσκω σωστά.

[V.I.Smirnov]

Δεν είδα τον κώδικα - κι ούτε σκοπεύω - αλλά
- το κριτήριο CFL για το βήμα ικανοποιείται ; (αναγκαία συνθήκη για τη σύγκλιση !)
- συνοριακές συνθήκες έβαλες (mur ή pml) ;

 

Δεν έχω μελετήσει τι γίνεται σε 4ης τάξης εξισώσεις αλλά τα παραπάνω είναι κρίσιμα θέματα.
Π.χ. αν για τις εξισώσεις 4ης τάξης χρησιμοποιείς το ίδιο βήμα όπως σε εκείνες της 2ης,
το πιθανότερο είναι ότι θα αποκλίνουν διότι δεν θα ικανοποιείται το CFL...

 

-

Επεξ/σία 12 Ιουλίου 2013 από V.I.Smirnov

[kostas123pao]

Aυτό που λες για το CFL δεν το χω κοιτάξει, μάλλον κάτι τέτοιο παίζει να προκαλεί το πρόβλημα.

[kostas123pao]

Ικανοποίησα το κριτήριο CFL αλλά τώρα έχω άλλο πρόβλημα, το πλάτος του v είναι αρκετά μικρότερο απο αυτό που ήταν με 2ης τάξης εξισώσεις.

[V.I.Smirnov]

Έχε υπόψην σου το εξής.
Πολλές φορές η άμεση ημιτονική διέγερση Α*sin(2πt) προξενεί μια παρασιτική DC συνιστώσα
η οποία οφείλεται σε αριθμητικά σφάλματα που επιφέρει η διακριτοποίηση (διασπορά, διάχυση)
και προφανώς είναι λαθος.
Δηλ. αν τεθεί E = Α*sin(2πt) η διακύμανση δεν θα είναι γύρω από το μηδέν όπως απλοϊκά αναμένεται
αλλά θα μετατοπιστεί γύρω από μιαν άλλη τιμή.
Η λύση είναι να μην τεθεί η διέγερση όπως είναι κατευθείαν αλλά να αυξηθεί γραμμικά για κάποιον χρόνο.
Π.χ.
Ε = Α*sin(2πt) * min(1, t*ramp_duration/Τ)
όπου Τ είναι η περίοδος και
ramp_duration είναι ο χρόνος (σε πολλαπλάσια της Τ) για να αναπτυχθεί πλήρως η κυματομορφή.

Το φαινόμενο αυτό υπάρχει και στις 2ης τάξης εξισώσεις.
Επειδή εφαρμόζεις 4ης τάξης, πιθανόν είναι ακόμη εντονότερο.
Έλενγξε αν η μέγιστη και η ελάχιστη τιμή που λαμβάνονται στην διάρκεια μιας περιόδου είναι απολύτως ίσες.
Αν όχι, (δηλ. αν δεν υπάρχει "κεντράρισμα" της κυματομορφής ως προς το μηδέν) εμφανίζεται DC συνιστώσα και
πρέπει να εφαρμόσεις το παραπάνω τέχνασμα....

-