Παραγωγή της μεθόδου FDTD

1. Εισαγωγή

Αυτή η σελίδα παρουσιάζει την λεπτομερή μαθηματική εξαγωγή των κύριων εξισώσεων ενημέρωσης της μεθόδου πεπερασμένων διαφορών στο πεδίο του χρόνου (FDTD), οι οποίες χρησιμοποιούνται για την επίλυση των εξισώσεων Maxwell στο πεδίο του χρόνου. Ο σκοπός αυτής της σελίδας διαφέρει από την κύρια σελίδα θεωρίας. Στόχος εδώ δεν είναι να παρέχεται μια συνοπτική επισκόπηση, αλλά να παρουσιαστεί ρητά πώς οι συνεχείς εξισώσεις στροβιλισμού μετατρέπονται σε διακριτές αλγεβρικές εκφράσεις που χρησιμοποιούνται σε έναν υπολογιστή. Αυτό είναι επίσης το τμήμα της μεθόδου FDTD που είναι πιο εύκολο να γίνει λάθος στον κώδικα: σφάλματα προσήμου, λανθασμένα τοποθετημένοι δείκτες ή ασυνεπής χρονική διαδοχή θα οδηγήσουν άμεσα σε ασταθή ή μη ρεαλιστικά αποτελέσματα.

Ως εκ τούτου, ξεκινάμε από τον νόμο του Ampère και τον νόμο του Faraday σε διαφορική μορφή, τους αναλύουμε σε καρτεσιανές συνιστώσες και στη συνέχεια τους εξειδικεύουμε σε ένα δισδιάστατο οπτικό πρόβλημα. Οι χωρικές παραγώγους αντικαθίστανται από κεντρικές πεπερασμένες διαφορές, οι χρονικές παραγώγους διακριτοποιούνται χρησιμοποιώντας το τυπικό σχήμα «leapfrog», και οι προκύπτουσες εκφράσεις αναδιατάσσονται για να δώσουν ρητές τύπους ενημέρωσης για τις συνιστώσες του ηλεκτρικού και του μαγνητικού πεδίου.

Για να διατηρήσουμε την παραγωγή ευανάγνωστη, απλοποιούμε τις εξισώσεις εδώ σε μια δισδιάστατη περίπτωση, όπου η \(\partial/\partial y = 0\) και η μαγνητική αγωγιμότητα ορίζονται σε μηδέν, \(\sigma_m = 0\). Αυτό καθιστά την άλγεβρα πολύ πιο εύκολη στην παρακολούθηση και δείχνει επίσης ρητά πώς πραγματοποιείται η μείωση των διαστάσεων στην πράξη. Στον πλήρη επιλυτή OghmaNano FDTD, ωστόσο, οι ηλεκτρομαγνητικές εξισώσεις επιλύονται σε τρεις διαστάσεις. Επομένως, η δισδιάστατη παραγωγή που παρουσιάζεται εδώ προορίζεται ως ένα σαφές παράδειγμα της μεθόδου, και όχι ως δήλωση ότι ο υποκείμενος επιλυτής περιορίζεται σε 2D.

Σε όλη αυτή τη σελίδα, οι συνιστώσες του ηλεκτρικού πεδίου υπολογίζονται σε ακέραια χρονικά βήματα \(n\Delta t\), ενώ οι συνιστώσες του μαγνητικού πεδίου υπολογίζονται σε ημι-ακέραια χρονικά βήματα \(\left(n+\tfrac{1}{2}\right)\Delta t\). Αυτή είναι η τυπική διάταξη Yee-grid leapfrog που χρησιμοποιείται στην μέθοδο FDTD. Σε αντίθεση με τη συντομότερη σελίδα επισκόπησης, η αρίθμηση δηλώνεται πιο ρητά εδώ, ώστε η εξαγωγή να μπορεί να διαβαστεί ως οδηγός υλοποίησης, καθώς και ως μαθηματική σημείωση.

2. Αρχικές εξισώσεις

Ξεκινάμε από τις εξισώσεις στροβιλισμού του Maxwell σε ένα γραμμικό μέσο με ηλεκτρική αγωγιμότητα \(\sigma\), ηλεκτρική διηλεκτρική σταθερά \(\epsilon\), μαγνητική διαπερατότητα \(\mu\) και, για λόγους γενικότητας, μαγνητική αγωγιμότητα \(\sigma_m\):

\[ \nabla \times \boldsymbol{H} = \sigma \boldsymbol{E} + \epsilon \frac{\partial \boldsymbol{E}}{\partial t}, \qquad \nabla \times \boldsymbol{E} = - \sigma_m \boldsymbol{H} - \mu \frac{\partial \boldsymbol{H}}{\partial t}. \]

Σε μορφή ορισμένου, οι μπούκλες είναι

\[ \nabla \times \boldsymbol{H} = \begin{vmatrix} \hat{\boldsymbol{x}} & \hat{\boldsymbol{y}} & \hat{\boldsymbol{z}} \\ \dfrac{\partial}{\partial x} & \dfrac{\partial}{\partial y} & \dfrac{\partial}{\partial z} \\ H_x & H_y & H_z \end{vmatrix}, \qquad \nabla \times \boldsymbol{E} = \begin{vmatrix} \hat{\boldsymbol{x}} & \hat{\boldsymbol{y}} & \hat{\boldsymbol{z}} \\ \dfrac{\partial}{\partial x} & \dfrac{\partial}{\partial y} & \dfrac{\partial}{\partial z} \\ E_x & E_y & E_z \end{vmatrix}. \]

Αυτές είναι οι εξισώσεις συνεχούς χρόνου και συνεχούς χώρου. Σκοπός της υπόλοιπης σελίδας είναι να τις μετατρέψουμε σε μια ρητή μέθοδο πεπερασμένων διαφορών, η οποία μπορεί να χρησιμοποιηθεί για την προσομοίωση της εξέλιξής τους στον χρόνο.

3. Η μορφή των εξισώσεων της στροβιλικότητας σε συντεταγμένες.

Αναπτύσσοντας τον νόμο του Ampère σε καρτεσιανές συνιστώσες, παίρνουμε:

\[ \sigma E_x + \epsilon \frac{\partial E_x}{\partial t} = \frac{\partial H_z}{\partial y} - \frac{\partial H_y}{\partial z}, \]

\[ \sigma E_y + \epsilon \frac{\partial E_y}{\partial t} = \frac{\partial H_x}{\partial z} - \frac{\partial H_z}{\partial x}, \]

\[ \sigma E_z + \epsilon \frac{\partial E_z}{\partial t} = \frac{\partial H_y}{\partial x} - \frac{\partial H_x}{\partial y}. \]

Αναπτύσσοντας τον νόμο του Faraday, προκύπτει:

\[ -\sigma_m H_x - \mu \frac{\partial H_x}{\partial t} = \frac{\partial E_z}{\partial y} - \frac{\partial E_y}{\partial z}, \]

\[ -\sigma_m H_y - \mu \frac{\partial H_y}{\partial t} = \frac{\partial E_x}{\partial z} - \frac{\partial E_z}{\partial x}, \]

\[ -\sigma_m H_z - \mu \frac{\partial H_z}{\partial t} = \frac{\partial E_y}{\partial x} - \frac{\partial E_x}{\partial y}. \]

Σε αυτό το στάδιο, δεν έχουμε κάνει ακόμη καμία υπόθεση σχετικά με τη διάσταση. Οι παραπάνω εκφράσεις είναι απλώς η μορφή των συνιστωσών των πλήρων τρισδιάστατων εξισώσεων στροβιλισμού.

4. Διδιάστατη μείωση.

Για να απλοποιήσουμε την εξαγωγή, περιορίζουμε τώρα το πρόβλημα σε δύο διαστάσεις, υποθέτοντας ότι δεν υπάρχει μεταβολή στην κατεύθυνση \(y\), επομένως \(\partial / \partial y = 0\). Αυτό το βήμα γίνεται αποκλειστικά για να απλοποιηθεί η άλγεβρα και να καταστεί ευκολότερη η κατανόηση της δομής των εξισώσεων ενημέρωσης. Η ίδια διαδικασία μπορεί να επεκταθεί άμεσα σε τρεις διαστάσεις και, στην πράξη, το OghmaNano επιλύει την πλήρη τρισδιάστατη μορφή των εξισώσεων του Maxwell.

Με αυτή την υπόθεση, οι εξισώσεις του ηλεκτρικού πεδίου απλοποιούνται σε:

\[ \sigma E_x + \epsilon \frac{\partial E_x}{\partial t} = - \frac{\partial H_y}{\partial z}, \]

\[ \sigma E_y + \epsilon \frac{\partial E_y}{\partial t} = - \frac{\partial H_z}{\partial x} + \frac{\partial H_x}{\partial z}, \]

\[ \sigma E_z + \epsilon \frac{\partial E_z}{\partial t} = \frac{\partial H_y}{\partial x}. \]

Εάν θέσουμε επίσης τη μαγνητική αγωγιμότητα στο μηδέν, \(\sigma_m = 0\), οι εξισώσεις του μαγνητικού πεδίου γίνονται:

\[ \mu \frac{\partial H_x}{\partial t} = \frac{\partial E_y}{\partial z}, \]

\[ \mu \frac{\partial H_y}{\partial t} = \frac{\partial E_z}{\partial x} - \frac{\partial E_x}{\partial z}, \]

\[ \mu \frac{\partial H_z}{\partial t} = - \frac{\partial E_y}{\partial x}. \]

Αυτές είναι οι συνεχείς εξισώσεις που θα διακριτοποιηθούν τώρα. Αξίζει να τονιστεί ρητά ότι πρόκειται για μια δισδιάστατη απλοποίηση του αρχικού τρισδιάστατου προβλήματος. Ο σκοπός δεν είναι να υποστηριχθεί ότι όλη η μέθοδος FDTD είναι εγγενώς δισδιάστατη, αλλά να εξαχθούν καθαρές εξισώσεις ανανέωσης για την κοινή περίπτωση όπου μία διάσταση παραμελείται.

5. Πλέγμα Yee, ευρετηρίαση και συμβολισμός.

Για να επιλύσουμε αριθμητικά τις παραπάνω εξισώσεις, τόσο ο χώρος όσο και ο χρόνος διακριτοποιούνται. Έστω ότι τα βήματα του υπολογιστικού πλέγματος στις κατευθύνσεις \(x\) και \(z\) είναι \(\Delta x\) και \(\Delta z\), αντίστοιχα, και έστω ότι το χρονικό βήμα είναι \(\Delta t\). Χρησιμοποιούμε τις ακέραιες τιμές \(i\) και \(k\) για τις κατευθύνσεις \(x\) και \(z\):

\[ x_i = i \Delta x, \qquad z_k = k \Delta z. \]

Ο χρόνος υποδηλώνεται με τον ακέραιο αριθμό \(n\), έτσι ώστε

\[ t^n = n \Delta t. \]

Στο σχήμα Yee, τα ηλεκτρικά και μαγνητικά πεδία δεν αποθηκεύονται στις ίδιες θέσεις, είτε στον χώρο είτε στον χρόνο. Αντίθετα, είναι διατεταγμένα με τρόπο που τα στοιχεία τους εναλλάσσονται. Οι συνιστώσες του ηλεκτρικού πεδίου υπολογίζονται σε ακέραιες χρονικές στάθμες \(n\Delta t\), ενώ οι συνιστώσες του μαγνητικού πεδίου υπολογίζονται σε ημι-ακέραιες χρονικές στάθμες \(\left(n+\tfrac{1}{2}\right)\Delta t\). Αυτό είναι η προέλευση της δομής "άλματος" (leapfrog) της μεθόδου FDTD.

Χωρικά, τα στοιχεία του πεδίου είναι επίσης μετατοπισμένα, έτσι ώστε κάθε κεντραρισμένη πεπερασμένη διαφορά της στροβιλικότητας να ορίζεται φυσικά μεταξύ γειτονικών θέσεων του πλέγματος. Σε μια πλήρη υλοποίηση του κώδικα, κάθε στοιχείο του πεδίου φέρει το δικό του ζεύγος χωρικών δεικτών, για παράδειγμα \(E_x^n(i,k)\), \(H_y^{n+\frac{1}{2}}(i+\tfrac{1}{2},k+\tfrac{1}{2})\) και ούτω καθεξής. Για να διατηρηθεί η αλγεβρική έκφραση ευανάγνωστη, ορισμένοι από τους επαναλαμβανόμενους «παρατηρητικούς» δείκτες δεν γράφονται πάντα ρητά στα παρακάτω βήματα. Ωστόσο, η ερμηνεία είναι πάντα η ίδια:

Το \(i\) αντιστοιχεί στην κατεύθυνση \(x\).
Το \(k\) αντιστοιχεί στην κατεύθυνση \(z\).
\(n\) υποδηλώνει τις διακριτές χρονικές στάθμες του ηλεκτρικού πεδίου, και
Το \(n+\tfrac{1}{2}\) καταγράφει τις χρονικές διακυμάνσεις του μαγνητικού πεδίου σε διαστήματα μισού τόνου.

Για παράδειγμα, η κεντρική προσέγγιση πεπερασμένων διαφορών για το \(\partial H_y / \partial z\) είναι η εξής:

\[ \frac{\partial H_y}{\partial z} \approx \frac{ H_y^{n+\frac{1}{2}}\!\left(i,k+\frac{1}{2}\right) - H_y^{n+\frac{1}{2}}\!\left(i,k-\frac{1}{2}\right) }{\Delta z}, \]

και παρομοίως, η κεντρική προσέγγιση του \(\partial H_z / \partial x\) είναι

\[ \frac{\partial H_z}{\partial x} \approx \frac{ H_z^{n+\frac{1}{2}}\!\left(i+\frac{1}{2},k\right) - H_z^{n+\frac{1}{2}}\!\left(i-\frac{1}{2},k\right) }{\Delta x}. \]

Με αυτή την σαφή διατύπωση, η παρακάτω διαδικασία μπορεί να ερμηνευτεί χωρίς αμφιβολία ως μια ορθή μέθοδος ενημέρωσης του πλέγματος Yee.

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

Το ηλεκτρικό πεδίο υπολογίζεται σε διακριτά χρονικά βήματα \(n\Delta t\), ενώ το μαγνητικό πεδίο υπολογίζεται σε ημι-διακριτά χρονικά βήματα \(\left(n+\tfrac{1}{2}\right)\Delta t\). Για να εξαχθούν οι εξισώσεις ενημέρωσης, η χρονική παράγωγος προσεγγίζεται με μια κεντρική πεπερασμένη διαφορά, και ο αγώγιμος όρος υπολογίζεται χρησιμοποιώντας τον μέσο όρο των παλιών και νέων ηλεκτρικών πεδίων. Αυτό παρέχει μια αριθμητικά σταθερή τραπεζοειδή προσέγγιση της αγωγιμότητας.

Αυτό το βήμα της μέσης τιμής είναι σημαντικό. Εάν ο όρος της αγωγιμότητας αντιμετωπιζόταν με ασυνεπή τρόπο, οι συντελεστές ενημέρωσης δεν θα είχαν πλέον την τυπική μορφή με αποσβεστικό παράγοντα και θα χανόταν η αντιστοιχία με τις συνήθεις εξισώσεις FDTD.

6.1 Ενημέρωση της εξίσωσης για το \(E_x\).

Ξεκινήστε από

\[ \sigma E_x + \epsilon \frac{\partial E_x}{\partial t} = - \frac{\partial H_y}{\partial z}. \]

Υπολογίστε τον όρο της αγωγιμότητας χρησιμοποιώντας τη μέση τιμή των \(E_x^{n+1}\) και \(E_x^n\) και υπολογίστε την παράγωγο του χρόνου χρησιμοποιώντας μια κεντρική πεπερασμένη διαφορά:

\[ \sigma \frac{E_x^{n+1}(i,k) + E_x^n(i,k)}{2} + \epsilon \frac{E_x^{n+1}(i,k) - E_x^n(i,k)}{\Delta t} = - \frac{ H_y^{n+\frac{1}{2}}\!\left(i,k+\frac{1}{2}\right) - H_y^{n+\frac{1}{2}}\!\left(i,k-\frac{1}{2}\right) }{\Delta z}. \]

Συγκεντρώστε τους όρους \(E_x^{n+1}\) στην αριστερή πλευρά:

\[ \left( \frac{\sigma}{2} + \frac{\epsilon}{\Delta t} \right) E_x^{n+1}(i,k) = - \frac{ H_y^{n+\frac{1}{2}}\!\left(i,k+\frac{1}{2}\right) - H_y^{n+\frac{1}{2}}\!\left(i,k-\frac{1}{2}\right) }{\Delta z} + \left( \frac{\epsilon}{\Delta t} - \frac{\sigma}{2} \right) E_x^n(i,k). \]

Επομένως

\[ E_x^{n+1}(i,k) = \frac{ \left( \dfrac{\epsilon}{\Delta t} - \dfrac{\sigma}{2} \right) E_x^n(i,k) - \dfrac{ H_y^{n+\frac{1}{2}}\!\left(i,k+\frac{1}{2}\right) - H_y^{n+\frac{1}{2}}\!\left(i,k-\frac{1}{2}\right) }{\Delta z} }{ \dfrac{\epsilon}{\Delta t} + \dfrac{\sigma}{2} }. \]

Πολλαπλασιάζοντας τον αριθμητή και τον παρονομαστή με το \(2\Delta t\), λαμβάνουμε την γνωστή μορφή συντελεστών:

\[ E_x^{n+1}(i,k) = \frac{ 2\epsilon - \sigma \Delta t }{ 2\epsilon + \sigma \Delta t } E_x^n(i,k) - \frac{ 2 \Delta t }{ 2\epsilon + \sigma \Delta t } \frac{ H_y^{n+\frac{1}{2}}\!\left(i,k+\frac{1}{2}\right) - H_y^{n+\frac{1}{2}}\!\left(i,k-\frac{1}{2}\right) }{\Delta z}. \]

6.2 Ενημέρωση της εξίσωσης για το \(E_y\).

Ξεκινήστε από

\[ \sigma E_y + \epsilon \frac{\partial E_y}{\partial t} = - \frac{\partial H_z}{\partial x} + \frac{\partial H_x}{\partial z}. \]

Διακριτοποιήστε και τους δύο όρους της παραγώγου στο πλέγμα Yee:

\[ \sigma \frac{E_y^{n+1}(i,k) + E_y^n(i,k)}{2} + \epsilon \frac{E_y^{n+1}(i,k) - E_y^n(i,k)}{\Delta t} = - \frac{ H_z^{n+\frac{1}{2}}\!\left(i+\frac{1}{2},k\right) - H_z^{n+\frac{1}{2}}\!\left(i-\frac{1}{2},k\right) }{\Delta x} + \frac{ H_x^{n+\frac{1}{2}}\!\left(i,k+\frac{1}{2}\right) - H_x^{n+\frac{1}{2}}\!\left(i,k-\frac{1}{2}\right) }{\Delta z}. \]

Η αναδιάταξη προσφέρει.

\[ E_y^{n+1}(i,k) = \frac{ 2\epsilon - \sigma \Delta t }{ 2\epsilon + \sigma \Delta t } E_y^n(i,k) + \frac{ 2 \Delta t }{ 2\epsilon + \sigma \Delta t } \left[ - \frac{ H_z^{n+\frac{1}{2}}\!\left(i+\frac{1}{2},k\right) - H_z^{n+\frac{1}{2}}\!\left(i-\frac{1}{2},k\right) }{\Delta x} + \frac{ H_x^{n+\frac{1}{2}}\!\left(i,k+\frac{1}{2}\right) - H_x^{n+\frac{1}{2}}\!\left(i,k-\frac{1}{2}\right) }{\Delta z} \right]. \]

Αυτή είναι η μόνη συνιστώσα του ηλεκτρικού πεδίου στην δισδιάστατη αναγωγή που περιέχει συνεισφορές από δύο ξεχωριστές παραγώγους του μαγνητικού πεδίου, και ως εκ τούτου είναι η πιο εύκολη να εφαρμοστεί λανθασμένα στον κώδικα. Συγκεκριμένα, η σχετική σχέση μεταξύ των όρων παραγώγων \(x\) και \(z\) προκύπτει άμεσα από την επέκταση της στροβιλικότητας και δεν θα πρέπει να τροποποιείται.

6.3 Ενημέρωση της εξίσωσης για το \(E_z\).

Ξεκινήστε από

\[ \sigma E_z + \epsilon \frac{\partial E_z}{\partial t} = \frac{\partial H_y}{\partial x}. \]

Η διακριτοποίηση δίνει

\[ \sigma \frac{E_z^{n+1}(i,k) + E_z^n(i,k)}{2} + \epsilon \frac{E_z^{n+1}(i,k) - E_z^n(i,k)}{\Delta t} = \frac{ H_y^{n+\frac{1}{2}}\!\left(i+\frac{1}{2},k\right) - H_y^{n+\frac{1}{2}}\!\left(i-\frac{1}{2},k\right) }{\Delta x}. \]

Επομένως

\[ E_z^{n+1}(i,k) = \frac{ 2\epsilon - \sigma \Delta t }{ 2\epsilon + \sigma \Delta t } E_z^n(i,k) + \frac{ 2 \Delta t }{ 2\epsilon + \sigma \Delta t } \frac{ H_y^{n+\frac{1}{2}}\!\left(i+\frac{1}{2},k\right) - H_y^{n+\frac{1}{2}}\!\left(i-\frac{1}{2},k\right) }{\Delta x}. \]

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

7. Εξαγωγή των εξισώσεων ενημέρωσης του μαγνητικού πεδίου.

Το μαγνητικό πεδίο ενημερώνεται σε διαστήματα χρόνου που είναι μισές ακέραιες τιμές. Με τη χρήση της \(\sigma_m = 0\), οι εξισώσεις του μαγνητικού πεδίου δεν περιέχουν όρο απόσβεσης, επομένως οι τύποι ενημέρωσης παίρνουν την απλούστερη μορφή πρόβλεψης. Αυτό το τμήμα της απόδειξης είναι αλγεβρικά συντομότερο από την περίπτωση του ηλεκτρικού πεδίου, επειδή δεν υπάρχει όρος τραπεζοειδούς αγωγιμότητας που να πρέπει να υπολογιστεί.

7.1 Ενημέρωση της εξίσωσης για το \(H_x\).

Ξεκινήστε από

\[ \mu \frac{\partial H_x}{\partial t} = \frac{\partial E_y}{\partial z}. \]

Η διακριτοποίηση στον χρόνο μεταξύ των ενδιάμεσων βημάτων \(n-\tfrac{1}{2}\) και \(n+\tfrac{1}{2}\) και στον χώρο με τη χρήση μιας κεντρικής διαφοράς στην κατεύθυνση \(z\) δίνει:

\[ \mu \frac{ H_x^{n+\frac{1}{2}}(i,k+\frac{1}{2}) - H_x^{n-\frac{1}{2}}(i,k+\frac{1}{2}) }{\Delta t} = \frac{ E_y^n(i,k+1) - E_y^n(i,k) }{\Delta z}. \]

Στη συνοπτική μορφή σημειογραφίας με μισό δείκτη, η οποία χρησιμοποιείται συχνά στις μαθηματικές αποδείξεις, το ίδιο αποτέλεσμα μπορεί να γραφτεί ως:

\[ \mu \frac{ H_x^{n+\frac{1}{2}} - H_x^{n-\frac{1}{2}} }{\Delta t} = \frac{ E_y^n\!\left(k+\frac{1}{2}\right) - E_y^n\!\left(k-\frac{1}{2}\right) }{\Delta z}. \]

Επομένως

\[ H_x^{n+\frac{1}{2}} = H_x^{n-\frac{1}{2}} + \frac{\Delta t}{\mu} \frac{ E_y^n\!\left(k+\frac{1}{2}\right) - E_y^n\!\left(k-\frac{1}{2}\right) }{\Delta z}. \]

7.2 Ενημέρωση της εξίσωσης για το \(H_y\).

Ξεκινήστε από

\[ \mu \frac{\partial H_y}{\partial t} = \frac{\partial E_z}{\partial x} - \frac{\partial E_x}{\partial z}. \]

Η διακριτοποίηση δίνει

\[ \mu \frac{ H_y^{n+\frac{1}{2}}(i+\frac{1}{2},k+\frac{1}{2}) - H_y^{n-\frac{1}{2}}(i+\frac{1}{2},k+\frac{1}{2}) }{\Delta t} = \frac{ E_z^n(i+1,k+\frac{1}{2}) - E_z^n(i,k+\frac{1}{2}) }{\Delta x} - \frac{ E_x^n(i+\frac{1}{2},k+1) - E_x^n(i+\frac{1}{2},k) }{\Delta z}. \]

Επιπλέον, χρησιμοποιώντας τη συντομευμένη μορφή της μερικής σημειογραφίας,

\[ \mu \frac{ H_y^{n+\frac{1}{2}} - H_y^{n-\frac{1}{2}} }{\Delta t} = \frac{ E_z^n\!\left(i+\frac{1}{2}\right) - E_z^n\!\left(i-\frac{1}{2}\right) }{\Delta x} - \frac{ E_x^n\!\left(k+\frac{1}{2}\right) - E_x^n\!\left(k-\frac{1}{2}\right) }{\Delta z}. \]

Επομένως

\[ H_y^{n+\frac{1}{2}} = H_y^{n-\frac{1}{2}} + \frac{\Delta t}{\mu} \left[ \frac{ E_z^n\!\left(i+\frac{1}{2}\right) - E_z^n\!\left(i-\frac{1}{2}\right) }{\Delta x} - \frac{ E_x^n\!\left(k+\frac{1}{2}\right) - E_x^n\!\left(k-\frac{1}{2}\right) }{\Delta z} \right]. \]

7.3 Ενημέρωση της εξίσωσης για το \(H_z\).

Ξεκινήστε από

\[ \mu \frac{\partial H_z}{\partial t} = - \frac{\partial E_y}{\partial x}. \]

Η διακριτοποίηση δίνει

\[ \mu \frac{ H_z^{n+\frac{1}{2}}(i+\frac{1}{2},k) - H_z^{n-\frac{1}{2}}(i+\frac{1}{2},k) }{\Delta t} = - \frac{ E_y^n(i+1,k) - E_y^n(i,k) }{\Delta x}. \]

Εναλλακτικά, σε συμπυκνωμένη μορφή,

\[ \mu \frac{ H_z^{n+\frac{1}{2}} - H_z^{n-\frac{1}{2}} }{\Delta t} = - \frac{ E_y^n\!\left(i+\frac{1}{2}\right) - E_y^n\!\left(i-\frac{1}{2}\right) }{\Delta x}. \]

Επομένως

\[ H_z^{n+\frac{1}{2}} = H_z^{n-\frac{1}{2}} - \frac{\Delta t}{\mu} \frac{ E_y^n\!\left(i+\frac{1}{2}\right) - E_y^n\!\left(i-\frac{1}{2}\right) }{\Delta x}. \]

Αφού ολοκληρωθούν αυτές οι αντικαταστάσεις, οι τυπικές εξισώσεις ενημέρωσης FDTD για την μέθοδο "leapfrog" προκύπτουν άμεσα.

8. Τελικές εξισώσεις ενημέρωσης.

Αφού συγκεντρώσουμε τα αποτελέσματα, οι δισδιάστατες εξισώσεις ανανέωσης της μεθόδου FDTD είναι οι εξής:

\[ H_x^{n+\frac{1}{2}}(i,k+\tfrac{1}{2}) = H_x^{n-\frac{1}{2}}(i,k+\tfrac{1}{2}) + \frac{\Delta t}{\mu} \frac{ E_y^n(i,k+1) - E_y^n(i,k) }{\Delta z}, \]

\[ H_y^{n+\frac{1}{2}}(i+\tfrac{1}{2},k+\tfrac{1}{2}) = H_y^{n-\frac{1}{2}}(i+\tfrac{1}{2},k+\tfrac{1}{2}) + \frac{\Delta t}{\mu} \left[ \frac{ E_z^n(i+1,k+\tfrac{1}{2}) - E_z^n(i,k+\tfrac{1}{2}) }{\Delta x} - \frac{ E_x^n(i+\tfrac{1}{2},k+1) - E_x^n(i+\tfrac{1}{2},k) }{\Delta z} \right], \]

\[ H_z^{n+\frac{1}{2}}(i+\tfrac{1}{2},k) = H_z^{n-\frac{1}{2}}(i+\tfrac{1}{2},k) - \frac{\Delta t}{\mu} \frac{ E_y^n(i+1,k) - E_y^n(i,k) }{\Delta x}. \]

Όπως διατυπώνεται, η έννοια του σχήματος «άλματος» είναι σαφής. Το ηλεκτρικό πεδίο στην χρονική στιγμή \(n+1\) προκύπτει από το μαγνητικό πεδίο στην χρονική στιγμή \(n+\tfrac{1}{2}\), ενώ το μαγνητικό πεδίο στην χρονική στιγμή \(n+\tfrac{1}{2}\) προκύπτει από το ηλεκτρικό πεδίο στην χρονική στιγμή \(n\). Επομένως, τα πεδία εξελίσσονται σε μια εναλλασσόμενη ακολουθία στον χρόνο.

9. Σχόλια σχετικά με την εφαρμογή.

Σε έναν πρακτικό κώδικα, οι συντελεστές που πολλαπλασιάζουν τις παλιές τιμές του πεδίου και τους όρους της στροβιλικότητας υπολογίζονται συχνά εκ των προτέρων και αποθηκεύονται, ειδικά όταν οι \(\epsilon\), \(\mu\) ή \(\sigma\) μεταβάλλονται από κελί σε κελί. Για παράδειγμα, η ενημέρωση του ηλεκτρικού πεδίου συνήθως έχει την εξής μορφή:

\[ E^{n+1} = C_a E^n + C_b (\nabla \times H), \]

όπου οι τοπικοί συντελεστές \(C_a\) και \(C_b\) κωδικοποιούν τις παραμέτρους του υλικού και το χρονικό βήμα. Αυτό καθιστά την υλοποίηση πιο απλή και αποφεύγει την επανάληψη του υπολογισμού των ίδιων παραγόντων σε κάθε χρονικό βήμα.

Επιπλέον, αξίζει να τονιστεί ότι η παραπάνω εξαγωγή παρέχει μόνο το βασικό σχήμα ενημέρωσης FDTD. Μια πρακτική προσομοίωση απαιτεί επιπλέον πηγές, οριακές συνθήκες, απορροφητικά στρώματα, όπως τα τέλεια προσαρμοσμένα στρώματα (PML), και ένα χρονικό βήμα που ικανοποιεί την συνθήκη σταθερότητας Courant. Καμία από αυτές τις προσθήκες δεν αλλάζει τη θεμελιώδη δομή άλματος που προκύπτει εδώ, αλλά καθεμία από αυτές τροποποιεί τον τρόπο με τον οποίο εφαρμόζονται οι εξισώσεις ενημέρωσης σε συγκεκριμένες θέσεις του πλέγματος.

Τέλος, παρόλο που η παραπάνω σημειογραφία έχει διατυπωθεί με αρκετή ακρίβεια ώστε να είναι χρήσιμη για την υλοποίηση, η ακριβής χωρική θέση κάθε στοιχείου εξαρτάται από τη διάταξη πλέγματος που χρησιμοποιείται στον κώδικα. Η κρίσιμη απαίτηση δεν είναι η ίδια η τυπογραφία, αλλά η συνέπεια: ο τελεστής curl, η σταδιακή μετατόπιση και τα πρόσημα πρέπει όλα να αντιστοιχούν στην ίδια διάταξη κελιών Yee σε όλο τον κώδικα.