[PYTHON] Numerische Analyse des elektrischen Feldes: Finite-Differenzen-Methode (FDM) -Praxis-
Vorwort
Fortsetzung der vorherigen "Numerische Analyse des elektrischen Feldes: Finite-Differenzen-Methode (FDM) -Basics-".
Dieses Mal werde ich endlich nicht nur mathematische Formeln schreiben, sondern auch Finite-Differenzen-Methoden mit Python. Trotzdem ist es mir ein wenig peinlich, dass ich nicht gut darin bin, Programme im Internet zu schreiben. Insbesondere in Python führt die Sauberkeit des Skripts direkt zu einer Beschleunigung und Optimierung, sodass es den Anschein hat, als würden verschiedene Dinge schnell erledigt.
In diesem Sinne sind C, C ++, Fortran usw. schneller, selbst wenn Sie die for-Schleife usw. in den Gehirnmuskel verwandeln, und in letzter Zeit hat mir dieser gefallen. ...... Die Geschichte ging schief. Beginnen wir mit der numerischen Analyse.
Das Programm, das ich tatsächlich erstellt habe, hat Github im Abschnitt "Zusammenfassung". Wenn Sie es also einmal ausführen möchten, laden Sie es bitte zuerst herunter und führen Sie es aus.
Problemstellung
Es ist schwierig, ein Programm zu schreiben, das alle Probleme von Anfang an lösen kann. Stellen wir also zuerst das Problem ein.
Poisson-Gleichung
Der Berechnungsbereich wird als das folgende gleichseitige Trapez definiert, wobei die obere und untere Seite als Richtungsgrenzbedingung und die linke und rechte als Neumann-Randbedingung definiert und wie folgt festgelegt werden.
** Abbildung 1 **: Berechnungsbereich
Lassen Sie uns alles zusammen schreiben.
\begin{align}
\nabla \cdot \nabla u &= 0 \ \ \text{(In der Gegend)}\\
u &= 1 \ \ \text{(Oberseite)} \\
u &= -1 \ \ \text{(Unterseite)} \\
\frac{\partial u}{\partial n} &= 0 \ \ \text{(Linke und rechte Seite)}
\end{align}
Schreiben wir ein Programm, das dieses Problem löst. Es ist jedoch langweilig, ein Programm zu schreiben, das nur dieses Problem löst. Lassen Sie uns es also implementieren, während wir einige allgemeine Aspekte berücksichtigen.
Die einzigen Bibliotheken, die dieses Mal verwendet werden, sind numpy und matplotlib.
Implementierung
Bevor wir das Programm sagen ... zeichnen wir zunächst ein kurzes Klassendiagramm. Selbst wenn Sie ein Klassendiagramm sagen, wird es nicht von allen geteilt, sondern zum Zweck der Organisation Ihres eigenen Gehirns geschrieben. Es ist also nicht schlampig geschrieben, aber es scheint, dass es eine solche Klasse gibt. Ich denke (Präventionslinie).
Lassen Sie uns zunächst den FDM-Fluss überprüfen.
--Schritt 1: Problemeinstellung --Bereichseinstellung / Einstellung des Quellelements --Schritt 2: Gitterteilung --Schritt 3: Lösen mit FDM - Generieren Sie simultane Gleichungen und lösen Sie sie --Schritt 4: Ergebnis
Ich habe diese vier ehrlich in große Klassen unterteilt.
Abbildung 2: Klassendiagramm
Lassen Sie mich kurz erklären. Die PoissonEquation-Klasse ist die sogenannte Hauptklasse. Dann wurde Problem als Problem festgelegt und in zwei Bereiche unterteilt, die Region Domain-Klasse und den Quellbegriff Source-Klasse. Randbedingungen können getrennt werden, aber ich habe beschlossen, sie in die Domain-Klasse aufzunehmen.
Als nächstes habe ich es in der Grid-Division-Klasse in die Grid-Klasse eingefügt. Das Raster besteht aus Klassen (Knoten), die Rasterpunkte verarbeiten. Ich habe die Node-Klasse absichtlich erstellt, weil ich daran gedacht habe, das Liniensegment (Kante) und den Bereich (Zelle) zu behandeln, die jeden Gitterpunkt verbinden, was in der Finite-Elemente-Methode verwendet wird, die ich beim nächsten Mal planen werde. ..
Als nächstes setzt die Solver-Klasse auch FDM darunter, weil sie auch andere Methoden wie FEM darunter legt. Sollte es sich um eine untergeordnete Klasse handeln, obwohl dies nicht in der Zeile angegeben ist? (Ehrlich gesagt, ich weiß es nicht, weil ich es noch nicht geschafft habe)
Schließlich wurde das Ergebnis in ein Potential (Potential) und ein elektrisches Feld (Flussdichte) unterteilt.
Okay, wir haben eine Richtung und lasst uns ein Programm schreiben! Es besteht jedoch eine hohe Wahrscheinlichkeit, dass es sich während der Herstellung nach und nach ändert.
Übrigens ist der Autor der Typ, der ab dem Ende der Klasse ausführlich schreibt, aber ist das üblich ...? Ich denke, es ist besser, ein Buch wie "Die königliche Art, ein Programm zu schreiben" zu lesen. Nun, diesmal werde ich nach meiner eigenen Methode vorgehen, also tut es mir leid, wenn es eine seltsame Methode ist.
Schritt 1: Problemstellung
Lassen Sie uns das Problem zunächst programmgesteuert festlegen. Das erste, worüber Sie zu diesem Zeitpunkt nachdenken sollten, ist, wie Sie es als Daten speichern können.
Bereichseinstellung (Domänenklasse)
Wenn beispielsweise berücksichtigt wird, dass es nur dem oben erwähnten gleichseitigen Trapezproblem entspricht, kann die Form eindeutig bestimmt werden, wenn nur "untere Seitenlänge", "obere Seitenlänge" und "Höhe" vorhanden sind, was als Daten ausreicht. .. Mit solchen Daten ist es jedoch möglich, gleichseitige Trapezoide, Rechtecke (Oberseite = Unterseite) und Quadrate (Oberseite = Unterseite = Höhe) zu realisieren, aber es ist nicht mit allen anderen Formen kompatibel und nicht allgemein. Uninteressant. Im Gegenteil, es ist schwierig, zu viel Allgemeinheit zu finden, und es ist fast unmöglich, alle Eckpunkte und Kurven mit einer endlichen Anzahl von Daten auszudrücken. Es kann möglich sein, eine Funktion zu übergeben, aber es ist offensichtlich, dass das Programm ziemlich kompliziert wird. Auf diese Weise kann das Streben nach Allgemeinheit das Programm komplizieren, und umgekehrt kann der Versuch, das Programm zu vereinfachen, die Allgemeinheit beeinträchtigen (nur meine empirische Regel, ich gebe zu, dass ich nicht einverstanden bin).
Lassen Sie uns diesmal die "Koordinaten der Polygonscheitelpunkte" als Array beibehalten. Auf diese Weise kann es alle Polygone verarbeiten, einen gewissen Grad an Allgemeinheit beibehalten und es scheint so einfach wie Daten zu sein. Ich kann mich nicht mit Kreisen befassen, aber ... nun, wenn Sie es zu einem regulären 20-seitigen Quadrat machen, können Sie es bis zu einem gewissen Grad approximieren, und Sie sollten es bei Bedarf erweitern. Die Quelle ist wie folgt.
domain.py
import numpy as np
import matplotlib.pyplot as plt
class Domain(object):
def __init__(self, shape, bc):
self.shape=shape
self.bc = bc
class Polygon(Domain):
def __init__(self, vertexes, bc, priority=None):
super().__init__("Polygon", bc)
self.nVertexes=len(vertexes)
self.vertexes = np.array(vertexes) #Wenn das Argument Scheitelpunkte eine Liste ist, np.Im Array ersetzen.
Ich habe absichtlich eine Domänenklasse erstellt und sie als untergeordnete Klasse zu einer Polygonklasse gemacht, da die Möglichkeit besteht, dass ich eine Kreisklasse erstelle, damit sie in Zukunft Kreisen entsprechen kann. Im Gegenteil, es ist auch möglich, eine Rechteckklasse zu erstellen und mit hoher Geschwindigkeit zu verarbeiten. Jedenfalls habe ich dies aus Gründen der Erweiterbarkeit getan. Der Autor ist kein Profi, sondern ein Amateur, daher sage ich nicht, dass dies korrekt ist (Präventionslinie). Ich denke, es ist eine gute Idee, eigenen Code zu schreiben.
Wie Sie dem obigen Code entnehmen können, haben die Daten, die wir haben, "Form", "bc" in der übergeordneten Klasse und "nVertexes", "Scheitelpunkte" in der kleinen Klasse. Wie Sie dem Variablennamen und dem Code entnehmen können, ist dieser wie folgt definiert.
--shape: Der Name der Form des zu verwendenden Bereichs. shape = "Polygon" --bc: Randbedingung --nVertexes: Anzahl der Polygonscheitelpunkte --vertices: Koordinaten der Vertices. Wenn Sie sie in der Reihenfolge von "Scheitelpunkten [0]" verbinden, können Sie ein Polygon erstellen.
Wie soll übrigens die Randbedingung ausgedrückt werden? Dieses Mal werde ich es im folgenden Wörterbuchtyp ausdrücken.
bc = {"bc":[{"bctype":"Dirichlet"Oder"Neumann" , "constant":Konstante}, ...], "priority":[...]}
Geben Sie die Konstante auf der rechten Seite der Randbedingung in "Konstante" und die Diricre / Neumann-Randbedingung in "bctype" an. "bc" wird als Liste von "Konstanten" und "bc-Typ" erstellt. "bc [" bctype "] [0]" bedeutet die Randbedingung der Grenze, die "Scheitelpunkte [0]" und "Scheitelpunkte [1]" verbindet. Wie Sie sehen können, werden nur Konstanten unterstützt. Dies ist nur ein Aufwand. Ich denke, dass andere als Konstanten leicht realisiert werden können, indem der Funktionstyp als Argument übergeben wird.
Welche Randbedingung priorisiert "Priorität"? Ich habe es in dem Sinne angehängt, dass. Welche der Randbedingungen an den Eckpunkten von "Scheitelpunkten [1]" ist beispielsweise "bc [" bctype "] [0]" oder "bc [" bctype "] [1]"? Es ist für die Angabe, dass. Die "Priorität" entschied, dass die Prioritäten in einer Liste in der richtigen Reihenfolge gespeichert werden sollten. Ich erkläre nicht im Detail, weil ich denke, dass die Regeln, die implementiert werden sollten, je nach Person unterschiedlich sind.
Zusätzlich haben wir "calcRange" hinzugefügt, das rechtwinklige Formen einschließlich Polygonen berechnet, und "plot", das zum Debuggen plottet. Es ist nicht notwendig, es zu schreiben, deshalb werde ich es weglassen.
Als Ergebnis wurde die folgende Figur erhalten.
Abbildung 3: Diagramm der Domänenklasse
Okay, der Berechnungsbereich ist jetzt eingestellt.
Quelltermspezifikation (Problemklasse)
Als nächstes erstellen wir die Quellklasse. Dies ist unkompliziert und kann einen Ausdruck (mit einem zweidimensionalen Array als Argument) als Argument verwenden. Die Quelle der Source-Klasse selbst lautet wie folgt.
Source.py
import numpy as np
class Source(object):
def __init__(self, source):
self.source=source
…… Das passt in eine Zeile. Ich denke nicht, dass diese Art der Klassifizierung gut ist ... Seien wir vorsichtig, um dies nicht zu tun! Es ist ein wenig abseits des Themas, aber Python kann die Mitgliedsvariablen der Klasse direkt von außen ändern und abrufen, sodass ich Getter und Setter nicht benötige. Ich habe den Eindruck, dass Python einfach zu bedienen ist, deshalb finde ich das wunderbar. Ich persönlich benutze Python für den Prototyp und C ++ für die Produktion.
Übrigens, wenn dies so bleibt, wie es ist, zum Beispiel, selbst wenn der Quellterm auf $ 0 $ gesetzt ist, muss die anonyme Funktion Lambda verwendet werden und function = lambda x: 0 muss als Argument angegeben werden.
Es ist ein bisschen nervig. Wenn der Quellterm ein konstanter Term ist, bearbeiten Sie ihn so, dass Sie eine Konstante übergeben können.
Der Code wird wie folgt geändert.
Source.py
import numpy as np
class Source(object):
def __init__(self, source):
print("-----Source-----")
if isinstance(source, (int, float)):
self.source = lambda x: source
else:
self.source = source
Wenn das Argument source vom Typ int oder float ist, scheint es als konstante Funktion zu funktionieren.
Domäne, Quelle aufrufen (Problemklasse)
Lassen Sie uns kurz über die Problemklasse schreiben, die die Domänen- und Quellklassen zusammenstellt, die auch als Zusammenfassung dient.
Problem.py
from . import Domain
from . import Source
class Problem(object):
def __init__(self, problem):
self.setDomain(**problem["domain"])
self.setSource(problem["source"])
def setDomain(self, **domain):
if domain["shape"]=="Polygon":
self.domain=Domain.Polygon(domain["vertexes"], domain["bc"])
self.domain.calcRange()
def setSource(self, source):
self.source=Source.Source(source)
Lassen Sie mich die Argumente kurz erläutern. Ein Beispiel für das Argument finden Sie auf der nächsten Seite.
setDomain(
shape = "Polygon"
vertexes = np.array([[-1,-1], [1,-1], [0.4,1],[-0.4,1]])
bc = {"bc": [{"bctype":"Dirichlet", "constant": 1},
{"bctype":"Neumann" , "constant": 0},
{"bctype":"Dirichlet", "constant":-1},
{"bctype":"Neumann" , "constant": 0}]
"priority":[0,2,1,3]}}
)
#Quelle ist eine konstante Funktion f=Bei Einstellung auf 0
setSource(
source = 0
)
#Quelle Sünde(x)*sin(y)Wenn du magst
setSource(
source = lambda x: np.sin(x[0])*np.sin(x[1])
)
Sie können diese Funktionen über die PoissonEquation-Klasse aufrufen.
Schritt 2: Gittererzeugung
Als nächstes wird der Bereich in Gitter unterteilt. Um ehrlich zu sein, war dies das schwierigste ... In diesem Abschnitt werden keine detaillierten Algorithmen behandelt, aber ich werde mich darauf konzentrieren, welche Art von Absicht und welche Art von Spezifikationen verwendet wurden. Weil [Computational Geometry](https://ja.wikipedia.org/wiki/%E8%A8%88%E7%AE%97%E5%B9%BE%E4%BD%95%E5%AD Ist es ein wenig vom Thema entfernt, weil es zu nahe an% A6) liegt? Weil ich dachte. Natürlich ist die Methode der Gitterteilung auch einer der Tiefenpunkte von FDM, aber sie ist zu tief, also werde ich sie für einen Moment verlassen und sie zu einer minimalen Geschichte machen. Wenn Sie interessiert sind, studieren Sie aus verschiedenen Büchern und Papieren.
Ich habe es in den folgenden 4 Schritten geschafft.
- Anordnung der Gitterpunkte (Knoten) horizontal zur $ x $ -Achse und $ y $ -Achse
- Setzen Sie einen Knoten an die Grenze
- Löschen Sie Knoten außerhalb des Bereichs
- Definition des benachbarten Knotens
Bearbeiten von Gitterpunkten (Knotenklasse)
Schritt 1. Platzieren Sie den Knoten
Anfangs dachte ich darüber nach, es in eine Node-Klasse und eine Node Management-Klasse zu unterteilen, aber wenn Sie sorgfältig darüber nachdenken, hat Python einen Wörterbuchtyp. Speichern Sie die Parameter jedes Gitterpunkts (Knotens) in einem Wörterbuchtyp und machen Sie die Node-Klasse zu einer Klasse, die die Menge manipuliert.
Die Parameter des Knotens sind wie folgt.
node={
"point":np.array(point),
"position":"nd",
"nextnode":[
{"no":None, "position":"nd"},
{"no":None, "position":"nd"},
...
]
}
Jede Erklärung ist unten zusammengefasst
- "Punkt": Knotenkoordinaten
- "Position": Positionsbeziehung der Knoten ("nd": undefiniert "in": innerhalb des Bereichs "b (Nummer n)": an der Grenze n "v (Nummer n)": am Scheitelpunkt n)
- "nextnode": Liste der Informationen für den nächsten Knoten
- "no": Index des nächsten Knotens
- "Position": Positionsbeziehung benachbarter Knoten ("n": undefiniert "l": links "r": rechts "d": unten "u": oben)
Cartesians Argument "__init __" (Konstruktor) erfordert eine Domäneninstanz und die Anzahl der Unterteilungen (mit dem Variablennamen "div"). Wenn Sie beispielsweise "[10,10]" übergeben, stellen Sie es so ein, dass es vertikal und horizontal in 10 gleiche Teile unterteilt ist. Dann wird es so.
Abbildung 4: Gleichmäßiges Gitter
Vielleicht möchten Sie auch etwas ausprobieren, das nicht gleichmäßig verteilt ist (wie? Nein?). Zum Beispiel: [[-10, -7, -3, -1,0,1,3,7,10], [-10, -7, -3, -1,0,1,3,7,10 ]] Wenn Sie es wie `übergeben, wird es gut damit umgehen,
Abbildung 5: Nicht äquidistantes Gitter
Ich habe versucht, das Gitter so zu teilen. Sie können es in der Praxis nicht verwenden, aber es ist gut für die Forschung. Dies kann nun wie folgt erreicht werden.
class Cartesian(Node):
def __init__(self, domain, div):
np.array(div)
if isinstance(div[0], (int)):
self.nDivX = div[0]
self.xs = np.linspace(domain.left, domain.right, self.nDivX)
else:
self.nDivX = len(div[0])
self.xs = np.array(div[0])
self.xs = np.sort(self.xs)
d = (domain.right - domain.left) / (self.xs[-1] - self.xs[0])
self.xs = d * self.xs
if isinstance(div[1], (int)):
self.nDivY = div[1]
self.ys = np.linspace(domain.down, domain.up , self.nDivY)
else:
self.nDivY = len(div[1])
self.ys = div[1]
self.ys = np.sort(self.ys)
d = (domain.up - domain.down) / (self.ys[-1] - self.ys[0])
self.ys = d * self.ys
self.nDiv = self.xs * self.ys
self.nodes = [{"point":np.array([self.xs[i],self.ys[j]]), "position":"nd", "nextnode":{"no":None, "position":"nd"} }
for j,i in itertools.product(range(self.nDivX), range(self.nDivY))]
Es ist ein Ärger, also werde ich mich nicht darum kümmern, es zu erklären. Die Leser werden in der Lage sein, ihre eigenen auf ihre eigene Weise zu machen. Denken Sie vorerst nur daran, dass Sie etwas erstellt haben, das einen Knoten wie den oben gezeigten erstellt.
Schritt 2. Setzen Sie einen Knoten an die Grenze
Als nächstes platzieren Sie die Knoten an der Grenze. Es spielt keine Rolle, ob es sich mit einem vorhandenen Knoten überschneidet. Platzieren Sie es auf jeden Fall dort, wo es den Horizont der $ x $, $ y $ -Achsen des vorhandenen Knotens schneidet. Ich habe dies als "getBorderPoint" in der Domain-Klasse gemacht. Diesmal ist es Polygon, aber ich dachte, der Algorithmus hat sich je nach Form geändert.
Wenn Sie nach dem Platzieren der Knoten an der Grenze die überlappenden Knoten löschen, ist das Ergebnis wie in der folgenden Abbildung dargestellt.
Abbildung 5: Grenzen gesetzt
... Ich habe einen Fehler beim Einstellen der Farbe gemacht. Es ist möglicherweise schwer zu erkennen, da sich der rote Rand und der magentafarbene Knoten überlappen. Bitte haben Sie etwas Geduld. Sie können sehen, dass der zweite Knoten von links und der zweite Knoten von oben sehr nahe an der Grenze liegen. Ich habe dies beim letzten Testen bemerkt, aber es verursacht viele Fehler, daher ist es sinnvoll, es zu entfernen. Wie nah möchten Sie löschen? Lassen Sie uns versuchen, etwas als Argument zu geben.
Ich habe es leicht geschrieben, aber es war unerwartet mühsam ... Wenn Sie die Quelle sehen möchten, sehen Sie sie bitte auf Github.
Schritt 3. Löschen Sie den äußeren Knoten
Als nächstes wird eine interne Beurteilungsverarbeitung durchgeführt, um Knoten außerhalb des Berechnungsbereichs zu löschen. Zunächst müssen wir das Innere und Äußere jedes einzelnen beurteilen, aber es wurden bereits verschiedene Algorithmen vorgeschlagen. Übrigens wird diesmal der Knoten an der Grenze in den Daten als "Position": "b" + (Nummer n) "im obigen Prozess gespeichert, also bearbeiten Sie nur die" Position ": nd". TU es einfach. Dies ist sehr hilfreich, da der Prozess der internen / externen Beurteilung an einem Knoten an der Grenze ziemlich mühsam ist.
Es gibt zwei berühmte Methoden. Eine besteht darin, eine halbe Linie in eine beliebige Richtung von dem zu beurteilenden Punkt zu zeichnen und die Anzahl der Kreuzungen mit der Grenze zu zählen (Crossing Number Algorithm). Wenn die Anzahl der Kreuzungen gerade ist, ist sie außerhalb und wenn sie ungerade ist, ist sie innen.
Abbildung 6: Crossing Number-Algorithmus
Diese Methode ist schnell, da sie nur vier Regeln enthält, aber es gibt einige Dinge, auf die Sie achten müssen. Wie Sie aus dem Beispiel von ③ in der obigen Abbildung sehen können, was sollten Sie tun, wenn der Scheitelpunkt das Liniensegment berührt? Darüber muss man nachdenken. Eigentlich wäre es gut, die Bedingung "Welche Richtung ist der Vektor der Grenzlinie, die kollidierte?" Auferlegen, aber es scheint problematisch zu sein.
Nehmen wir dieses Mal eine andere Methode namens Winding Number Algorithm an. Wenn aus Sicht des Urteils die Summe der Winkel zwischen jedem Scheitelpunkt $ 2π $ beträgt, befindet sie sich innerhalb und wenn sie $ 0 $ beträgt, befindet sie sich außerhalb.
Abbildung 7: Wicklungszahlalgorithmus
Da diese Methode eine inverse Dreiecksfunktion verwendet, dauert es etwas, aber es ist nicht erforderlich, die obigen Bedingungen zu berücksichtigen.
Nun, ich schreibe ein Programm, aber es gibt keinen Trick, auf jedem Knoten eine for-Schleife auszuführen. Da ich numpy verwende, werde ich es schaffen, die Anzahl der Zeilen und die Berechnungskosten zu reduzieren.
Das Programm war wie folgt.
class Polygon(Domain):
#Kürzung
def deleteOutDomain(self, node, thetaerr=0.1):
pts = np.array([node[i]["point"] for i in range(len(node))]) #["point"]In numpy ersetzt
pos = np.array([node[i]["position"] for i in range(len(node))]) #["position"]In numpy ersetzt
thetas = 0
for i in range(len(self.vertexes)):
p1 = self.vertexes[i] #Linienendpunkt 1 p1
p2 = self.vertexes[(i + 1) % self.nVertexes] #Linienendpunkt 2 p2
v1 = (pts - p1) #Vektor von dem Punkt, den Sie beurteilen möchten, bis p1
v2 = (pts - p2) #Vektor von dem Punkt, den Sie beurteilen möchten, bis p2
dot = v1[:,0] * v2[:,0] + v1[:,1] * v2[:,1] #Innenprodukt(dot,Ich habe nicht auf innere ...)
cross = v1[:,0] * v2[:,1] - v1[:,1] * v2[:,0] #Äußeres Produkt(Das gleiche wie oben)
vn1 = np.linalg.norm(v1, axis=1) #Holen Sie sich v1 Abstand
vn2 = np.linalg.norm(v2, axis=1) #Holen Sie sich v2 Abstand
theta = np.arccos(np.clip(dot / (vn1 * vn2), -1, 1)) #Clip, da er aufgrund eines numerischen Fehlers mehr als 1 sein kann
thetas += np.sign(cross)*np.array(theta) #Fügen Sie den Winkel jedes Mal für hinzu(Kreuz ist ein Plus / Minus-Urteil)
inx = np.where( ((pos == "nd") & ~(thetas < thetaerr)))[0] # ["position"]Is nd und thetas sind nicht fast 0, sondern erhalten den Index
#Kürzung
Wenn es gut gemacht wird, kann es möglich sein, sogar die "Vertexes" -Schleife zu eliminieren, aber ich denke, wenn Sie so viel tun, wird es Ihnen vergeben.
Das Ergebnis ist wie folgt.
Abbildung 8: Die Außenseite wird gelöscht
4. Definieren wir den nächsten Knoten
Zum Schluss definieren wir node [next node]. Die Definition vereinfacht die Matrixgenerierung für den nächsten Schritt.
Ich machte mir Sorgen um verschiedene Dinge, beschloss aber, nach Koordinaten nach oben, unten, links und rechts zu urteilen.
Gleitkomma == (eigentlich np.is close) kommt heraus, daher ist es vielleicht nicht wirklich gut, aber mir fiel nichts anderes ein.
Vielleicht hätte ich zuerst einen ganzzahligen Index als Daten haben sollen, aber bitte verzeihen Sie mir.
Und vorher sortieren wir die Knotenliste nach Koordinaten. Ich denke es ist einfacher.
Sie müssen sich nicht die Mühe machen, es zu erhöhen, also beziehen Sie sich bitte auf Github. (Wenn ich ein Programm mache, sage ich "Nicht verwenden" für "!" Und es ist peinlich, weil es dupliziert ist ...)
Nachdem wir dies getan haben, geben wir das Ergebnis der "Knoten" -Liste aus. Selbst wenn die Anzahl der Daten zu groß ist, können sie nicht umfassend angezeigt werden. Lassen Sie es uns also um "div = [4,4]" setzen.
node0 {'point': array([-1., -1.]), 'position': 'c0', 'nextnode': [{'no': 1, 'position': 'r'}]}
node1 {'point': array([-0.333, -1. ]), 'position': 'b0', 'nextnode': [{'no': 5, 'position': 'u'}, {'no': 0, 'position': 'l'}, {'no': 2, 'position': 'r'}]}
node2 {'point': array([ 0.333, -1. ]), 'position': 'b0', 'nextnode': [{'no': 6, 'position': 'u'}, {'no': 1, 'position': 'l'}, {'no': 3, 'position': 'r'}]}
node3 {'point': array([ 1., -1.]), 'position': 'c1', 'nextnode': [{'no': 2, 'position': 'l'}]}
node4 {'point': array([-0.8 , -0.333]), 'position': 'b3', 'nextnode': [{'no': 5, 'position': 'r'}]}
node5 {'point': array([-0.333, -0.333]), 'position': 'in', 'nextnode': [{'no': 1, 'position': 'd'}, {'no': 9, 'position': 'u'}, {'no': 4, 'position': 'l'}, {'no': 6, 'position': 'r'}]}
node6 {'point': array([ 0.333, -0.333]), 'position': 'in', 'nextnode': [{'no': 2, 'position': 'd'}, {'no': 10, 'position': 'u'}, {'no': 5, 'position': 'l'}, {'no': 7, 'position': 'r'}]}
node7 {'point': array([ 0.8 , -0.333]), 'position': 'b1', 'nextnode': [{'no': 6, 'position': 'l'}]}
node8 {'point': array([-0.6 , 0.333]), 'position': 'b3', 'nextnode': [{'no': 9, 'position': 'r'}]}
node9 {'point': array([-0.333, 0.333]), 'position': 'in', 'nextnode': [{'no': 5, 'position': 'd'}, {'no': 13, 'position': 'u'}, {'no': 8, 'position': 'l'}, {'no': 10, 'position': 'r'}]}
node10 {'point': array([0.333, 0.333]), 'position': 'in', 'nextnode': [{'no': 6, 'position': 'd'}, {'no': 14, 'position': 'u'}, {'no': 9, 'position': 'l'}, {'no': 11, 'position': 'r'}]}
node11 {'point': array([0.6 , 0.333]), 'position': 'b1', 'nextnode': [{'no': 10, 'position': 'l'}]}
node12 {'point': array([-0.4, 1. ]), 'position': 'c3', 'nextnode': [{'no': 13, 'position': 'r'}]}
node13 {'point': array([-0.333, 1. ]), 'position': 'b2', 'nextnode': [{'no': 9, 'position': 'd'}, {'no': 12, 'position': 'l'}, {'no': 14, 'position': 'r'}]}
node14 {'point': array([0.333, 1. ]), 'position': 'b2', 'nextnode': [{'no': 10, 'position': 'd'}, {'no': 13, 'position': 'l'}, {'no': 15, 'position': 'r'}]}
Ja, es fühlt sich gut an ... oh! Ich habe vergessen. Ich wollte die Knoten an der Grenze als "f" (vorwärts) und "b" (zurück) darstellen. Da es an der Grenze diagonal ist, kann es nicht nur durch Auf, Ab, Links und Rechts ausgedrückt werden. Also habe ich das Programm nochmal überprüft.
node0 {'point': array([-1., -1.]), 'position': 'c0', 'nextnode': [{'no': 1, 'position': 'f'}, {'no': 4, 'position': 'b'}, {'no': 1, 'position': 'r'}]}
node1 {'point': array([-0.333, -1. ]), 'position': 'b0', 'nextnode': [{'no': 0, 'position': 'b'}, {'no': 2, 'position': 'f'}, {'no': 0, 'position': 'l'}, {'no': 2, 'position': 'r'}]}
node2 {'point': array([ 0.333, -1. ]), 'position': 'b0', 'nextnode': [{'no': 1, 'position': 'b'}, {'no': 3, 'position': 'f'}, {'no': 1, 'position': 'l'}, {'no': 3, 'position': 'r'}]}
node3 {'point': array([ 1., -1.]), 'position': 'c1', 'nextnode': [{'no': 2, 'position': 'b'}, {'no': 7, 'position': 'f'}, {'no': 2, 'position': 'l'}]}
node4 {'point': array([-0.8 , -0.333]), 'position': 'b3', 'nextnode': [{'no': 0, 'position': 'f'}, {'no': 8, 'position': 'b'}, {'no': 5, 'position': 'r'}]}
node5 {'point': array([-0.333, -0.333]), 'position': 'in', 'nextnode': [{'no': 4, 'position': 'l'}, {'no': 6, 'position': 'r'}]}
node6 {'point': array([ 0.333, -0.333]), 'position': 'in', 'nextnode': [{'no': 5, 'position': 'l'}, {'no': 7, 'position': 'r'}]}
node7 {'point': array([ 0.8 , -0.333]), 'position': 'b1', 'nextnode': [{'no': 3, 'position': 'b'}, {'no': 11, 'position': 'f'}, {'no': 6, 'position': 'l'}]}
node8 {'point': array([-0.6 , 0.333]), 'position': 'b3', 'nextnode': [{'no': 4, 'position': 'f'}, {'no': 12, 'position': 'b'}, {'no': 9, 'position': 'r'}]}
node9 {'point': array([-0.333, 0.333]), 'position': 'in', 'nextnode': [{'no': 8, 'position': 'l'}, {'no': 10, 'position': 'r'}]}
node10 {'point': array([0.333, 0.333]), 'position': 'in', 'nextnode': [{'no': 9, 'position': 'l'}, {'no': 11, 'position': 'r'}]}
node11 {'point': array([0.6 , 0.333]), 'position': 'b1', 'nextnode': [{'no': 7, 'position': 'b'}, {'no': 15, 'position': 'f'}, {'no': 10, 'position': 'l'}]}
node12 {'point': array([-0.4, 1. ]), 'position': 'c3', 'nextnode': [{'no': 13, 'position': 'b'}, {'no': 8, 'position': 'f'}, {'no': 13, 'position': 'r'}]}
node13 {'point': array([-0.333, 1. ]), 'position': 'b2', 'nextnode': [{'no': 12, 'position': 'f'}, {'no': 14, 'position': 'b'}, {'no': 12, 'position': 'l'}, {'no': 14, 'position': 'r'}]}
node14 {'point': array([0.333, 1. ]), 'position': 'b2', 'nextnode': [{'no': 13, 'position': 'f'}, {'no': 15, 'position': 'b'}, {'no': 13, 'position': 'l'}, {'no': 15, 'position': 'r'}]}
node15 {'point': array([0.4, 1. ]), 'position': 'c2', 'nextnode': [{'no': 11, 'position': 'b'}, {'no': 14, 'position': 'f'}, {'no': 14, 'position': 'l'}]}
Okay, es ist in Ordnung. Lassen Sie uns die Form ein wenig ändern und es versuchen.
Abbildung 9: Wenn die Form geändert wird
Gute Stimmung.
Schritt 3: FDM-Löser
Erstellen Sie eine Matrix (Matrix Generator)
Ich bin endlich da. Ich muss nur eine Warteschlange erstellen ... aber hier gibt es zwei Probleme.
Diskreter Ausdruck innerhalb des Berechnungsbereichs
Die erste ist die diskrete Gleichung, die beim letzten Mal abgeleitet wurde.
Wie in der folgenden Abbildung gezeigt, werden die linken, rechten, oberen und unteren Knotenabstände auf $ \ Delta_l, \ Delta_r, \ Delta_u. \ Delta_d $ festgelegt. Wenn zu diesem Zeitpunkt der Wert von $ u $ an jedem Punkt $ u_l, u_r, u_u. U_d $ ist und das Zentrum $ u_0 $ für die Taylor-Erweiterung verwendet wird, ist das Ergebnis wie folgt.
\begin{align}
u_l &= u_0 - \Delta_l \frac{\partial u_0}{\partial x} + \frac{\Delta_l^2}{2} \frac{\partial^2 u_0}{\partial x^2} -\frac{\Delta_l^3}{6} \frac{\partial^3 u_0}{\partial x^3} + \mathcal{O}(\Delta_x^4) \\
u_r &= u_0 + \Delta_r \frac{\partial u_0}{\partial x} + \frac{\Delta_r^2}{2} \frac{\partial^2 u_0}{\partial x^2} +\frac{\Delta_r^3}{6} \frac{\partial^3 u_0}{\partial x^3} + \mathcal{O}(\Delta_x^4) \\
u_d &= u_0 - \Delta_u \frac{\partial u_0}{\partial y} + \frac{\Delta_u^2}{2} \frac{\partial^2 u_0}{\partial y^2} -\frac{\Delta_u^3}{6} \frac{\partial^3 u_0}{\partial y^3} + \mathcal{O}(\Delta_y^4) \\
u_u &= u_0 + \Delta_d \frac{\partial u_0}{\partial y} + \frac{\Delta_d^2}{2} \frac{\partial^2 u_0}{\partial y^2} +\frac{\Delta_d^3}{6} \frac{\partial^3 u_0}{\partial y^3} + \mathcal{O}(\Delta_y^4)
\end{align}
Wenn Sie danach den Begriff der Differenzierung erster Ordnung gut addieren und subtrahieren und löschen, können Sie ihn wie folgt finden.
\begin{align}
\frac{\partial^2 u}{\partial x^2} &= \frac{2}{\Delta_l \Delta_r(\Delta_l + \Delta_r)}\left( \Delta_r(u_l - u_0) + \Delta_l (u_r - u_0) \right) + \mathcal{O}(\Delta_x)\\
\frac{\partial^2 u}{\partial y^2} &= \frac{2}{\Delta_d \Delta_u(\Delta_d + \Delta_u)}\left( \Delta_u(u_d - u_0) + \Delta_d (u_u - u_0) \right) + \mathcal{O}(\Delta_y)
\end{align}
Leider ist die Genauigkeit zur primären Genauigkeit geworden, aber es kann nicht geholfen werden. Die folgende Gleichung kann erhalten werden, indem sie in die Poisson-Gleichung eingesetzt wird.
\begin{multline}
\frac{2}{\Delta_l \Delta_r(\Delta_l + \Delta_r)}\left( \Delta_r(u_l - u_0) + \Delta_l (u_r - u_0) \right) \\+ \frac{2}{\Delta_d \Delta_u(\Delta_d + \Delta_u)}\left( \Delta_u(u_d - u_0) + \Delta_d (u_u - u_0) \right) = f + \mathcal{O}(\Delta_{x,y})
\end{multline}
Es sieht unangenehm aus, aber ein Computer kann es sofort berechnen.
Zu beachten ist die Größe des Ausdrucks $ \ mathcal {O} (\ Delta) $. Wie Sie intuitiv sehen können, ist der Fehler umso größer, je größer der Unterschied zwischen dem linken und rechten (oberen und unteren) Intervall ist. Die tatsächliche Berechnung zeigt, dass $ \ mathcal {O} (\ Delta_x) $ entsprechend zunimmt, wenn der Unterschied zwischen $ \ Delta_l $ und $ \ Delta_r $ groß ist. Es ist nicht so schwierig, und wenn Sie frei sind, versuchen Sie es zu berechnen (diejenigen, die es durch geheime Berechnung verstehen können).
Um es anders herum auszudrücken, ist es besser, wenn der Rasterabstand so groß wie möglich ist. Dinge wie Abbildung 5 sind nicht wirklich gut. Eigentlich ist es besser, die Anzahl der Knoten zu erhöhen oder zu verringern.
Diskreter Ausdruck an der Grenze
Das andere Problem ist problematischer. Ungefähre Formel der zuletzt abgeleiteten Neumann-Randbedingung
\begin{align}
\frac{u_{m+1,n} - u_{m,n}}{\Delta_x} = g \\
\frac{u_{m,n+1} - u_{m,n}}{\Delta_y} = g
\end{align}
Ist nur eine Annäherung, wenn die Grenze senkrecht zur $ x, y $ -Achse ist (dh teilweise Differenzierung von $ x $ oder $ y $) und kann in einem leicht schrägen Fall wie diesem nicht verwendet werden.
Abbildung 10: Diagonale Randbedingungen
Es wäre schön, wenn es in der obigen Abbildung einen Knoten in Pfeilrichtung gäbe, aber das ist nicht der Fall. Kehren wir also zu den Grundlagen zurück und denken wir an die Taylor-Entwicklung.
Und davor
Mit anderen Worten, irgendwie sollte die folgende Dezentralisierung durchgeführt werden.
\begin{align}
u_{r, s}&=\sum_{p=0}^\infty \frac{1}{p!}\left( r \frac{\partial}{\partial x} + s \frac{\partial}{\partial y} \right)^p u_0 \\
&=u_{0,0} + r \frac{\partial u_0}{\partial x} + s \frac{\partial u_0}{\partial y} + \cdots
\end{align}
Der zweite und dritte Term sollten mit Gleichung (1) identisch sein. Dies kann erreicht werden, indem zwei von Taylor entwickelte Formeln an zwei Punkten kombiniert werden. Wenn Sie es wagen, es redundant zu schreiben,
\begin{align}
u_{r_1, s_1}=u_{0,0} + r_1 \frac{\partial u_0}{\partial x} + s_1 \frac{\partial u_0}{\partial y} + \cdots \\
u_{r_2, s_2}=u_{0,0} + r_2 \frac{\partial u_0}{\partial x} + s_2 \frac{\partial u_0}{\partial y} + \cdots
\end{align}
Sie können der Formel (1) 3 Elemente und 4 Elemente hinzufügen oder davon abziehen. Mit anderen Worten sollte die folgende Gleichung für $ c_1 $ und $ c_2 $ gelöst werden.
\begin{align}
r_1c_1 + r_2c_2 = a \\
s_1c_1 + s_2c_2 = b
\end{align}
Wenn wir nach $ c_1 $ und $ c_2 $ auflösen, können wir dies wie folgt approximieren.
$ \ mathcal {E} $ repräsentiert den Fehler. Ich wusste nicht, wie ich es richtig schreiben sollte, also habe ich es richtig geschrieben.
Nun, es ist einfach, wenn Sie dies tun können. Verwenden Sie für $ u_ {r_1, s_1} $, $ u_ {r_2, s_2} $ die beiden in der Nähe der Grenze. Sie können eine Matrix dazu erstellen und mit numpy.linalg.solve lösen.
Es gibt jedoch 3 Knoten in der Nähe der Grenze, und es ist eine Verschwendung, mit 2 Knoten zu lösen. Wenn Sie mit zwei lösen können, können Sie natürlich genauer lösen, wenn Sie drei haben.
Mit anderen Worten
Ist festgelegt. Wenn $ g $ eine Konstante ist, ist die rechte Seite $ 0 $. Als $ g = 0 $
Ist. Wenn diese Formel mit drei Punkten erstellt wird, kann eine sekundäre Genauigkeit erreicht werden. Schreiben wir vorerst die drei folgenden Formeln auf.
\begin{align}
u_{r_1, s_1}=u_{0,0} + r_1 \frac{\partial u_0}{\partial x} + s_1 \frac{\partial u_0}{\partial y} + \frac{r_1^2}{2} \frac{\partial^2 u_0}{\partial x^2} + r_1s_1 \frac{\partial^2 u_0}{\partial x \partial y} + \frac{s_1^2}{2} \frac{\partial^2 u_0}{\partial y^2} + \cdots \\
u_{r_2, s_2}=u_{0,0} + r_2 \frac{\partial u_0}{\partial x} + s_2 \frac{\partial u_0}{\partial y} + \frac{r_2^2}{2} \frac{\partial^2 u_0}{\partial x^2} + r_2s_2 \frac{\partial^2 u_0}{\partial x \partial y} + \frac{s_2^2}{2} \frac{\partial^2 u_0}{\partial y^2} + \cdots \\
u_{r_3, s_3}=u_{0,0} + r_3 \frac{\partial u_0}{\partial x} + s_3 \frac{\partial u_0}{\partial y} + \frac{r_3^2}{2} \frac{\partial^2 u_0}{\partial x^2} + r_3s_3 \frac{\partial^2 u_0}{\partial x \partial y} + \frac{s_3^2}{2} \frac{\partial^3 u_0}{\partial y^2} + \cdots
\end{align}
Sie mögen sich fragen, ob Gleichung (2) durch Gleichung 3 gelöst werden kann, während es sich um 5 Terme handelt, aber $ k_1 $ und $ k_2 $ sind willkürliche Werte, sodass sie verwaltet werden können. Schreiben wir es als simultane Gleichung auf. Zunächst werden die folgenden Begriffe erstellt, indem jeder Begriff mit einem Koeffizienten multipliziert und addiert wird.
\begin{equation}
\left\{
\begin{alignedat}{5}
r_1 &c_1& + r_2 &c_2 &+ r_3 c_3&=& a \\
s_1 &c_1& + s_2 &c_2 &+ s_3 c_3&=& b \\
r_1^2&c_1& + r_2^2 &c_2 &+ r_3^2 c_3&=& 2 k_1a \\
r_1 s_1 &c_1& + r_2 s_2 &c_2 &+ r_3 s_3 c_3&=& k_1 b + k_2 a \\
s_1^2&c_1& + s_2^2 &c_2 &+ s_3^2 c_3&=& 2 k_2b
\end{alignedat}
\right.
\end{equation}
Wenn Sie diese $ k_1 $ und $ k_2 $ als Variablen betrachten, haben Sie 5 Variablen und 5 Gleichungen. Durch Umschreiben in eine Matrix können wir die folgende Formel erstellen:
\begin{equation}
\left[
\begin{matrix}
r_1 & r_2 & r_3 & 0 & 0 \\
s_1 & s_2 & s_3 & 0 & 0 \\
r_1^2 & r_2^2 & r_3^2 & -2a & 0 \\
r_1s_1 & r_2s_2& r_3s_3& -b & -a \\
s_1^2 & s_2^2 & s_3^2 & 0 & -2b \\
\end{matrix}
\right]\left[
\begin{matrix}
c_1 \\
c_2 \\
c_3 \\
k_1 \\
k_2
\end{matrix}
\right] =
\left[
\begin{matrix}
a \\
b \\
0 \\
0 \\
0
\end{matrix}
\right]
\end{equation}
Dies scheint irgendwie gelöst zu sein.
Schritt 4: Ergebnisausgabe
Nun wird endlich das Ergebnis ausgegeben.
Es ist das Wichtigste, um den Menschen Ergebnisse zu zeigen, und ich freue mich sehr darauf, auch wenn die Entwickler es wissen. Wenn hier jedoch schlechte Ergebnisse erzielt werden, ist das Gefühl der Verzweiflung groß. Nein, es ist unwahrscheinlich, dass Sie das erste Mal Erfolg haben ...
Potenzial
Bevor wir das Ergebnis annehmen, stellen wir uns die Antwort für einen Moment vor. Das Problem war, dass die Ober- und Unterseite des Trapezes unter der Diricre-Randbedingung 1 bzw. -1 waren. Was wäre, wenn es ein Quadrat wäre? Dies ist genau das gleiche wie bei einem Kondensator im Idealzustand, dh das elektrische Feld ist von oben nach unten konstant und das Potential ist proportional. Wenn dann die Oberseite kleiner wird ... Also habe ich unten ein Bilddiagramm geschrieben.
Abbildung 10: Bild der Lösung
Der grüne Pfeil ist das elektrische Feld, und die hellen regenbogenfarbenen Linien repräsentieren die Konturlinien des Potentials. Wenn es so aussieht, konnte ich die richtige Antwort bekommen! Vorerst habe ich es mit 100x100 versucht.
Ich habe es so geschrieben, als ob es sofort gemacht worden wäre, aber es ist das Ergebnis der mehrmaligen Zerstörung von Fehlern und der Behebung mit einigem Einfallsreichtum.
Elektrisches Feld (Flussdichte)
Lassen Sie uns auch das elektrische Feld ausgeben. Dieses Mal berechnen wir das elektrische Feld an den Gitterpunkten aus den vier umgebenden Potentialen. Mit anderen Worten
\begin{align}
\left. \frac{\partial u}{\partial x} \right|_{m,n} = \frac{1}{2\Delta_x}(u_{m+1, n} - u_{m-1, n}) \\
\left. \frac{\partial u}{\partial y} \right|_{m,n} = \frac{1}{2\Delta_y}(u_{m, n+1} - u_{m, n-1})
\end{align}
Berechnen Sie wie folgt. Das Ergebnis ist wie folgt.
Es fühlt sich sehr gut an. Ich wollte es auch mit der theoretischen Lösung vergleichen, aber dieses Mal ist es lang geworden, also werde ich es aufhalten. Wenn ich eines Tages Lust dazu habe ...
Poisson-Gleichungsklasse
Ich habe eine PoissonEquation-Klasse erstellt, die alle oben genannten Punkte zusammenfasst.
#Kürzung
if __name__ == "__main__":
#Step 1: Problem
domain = {"shape":"Polygon",
"vertexes":[[-1,-1],[1,-1], [0.4,1],[-0.4,1]],
"bc":{
"bc": [
{"bctype":"Dirichlet", "constant":-1},
{"bctype":"Neumann", "constant":0},
{"bctype":"Dirichlet", "constant":1},
{"bctype":"Neumann", "constant":0}
],
"priority":[0,2,1,3]
}
}
source = {"source": 0}
poisson = PoissonEquation(domain,source)
# Step 2: generateGrid
grid = {"type":"Cartesian", "div":[25,25]}
poisson.generateGrid(grid, "Node")
# Step 3: Solve
method = "FDM"
poisson.solve(method)
# Step 4: Result
poisson.result()
poisson.plot("Potential")
poisson.plot("DensityFlux")
Zusammenfassung
Dieses Mal haben wir den FDM-Löser für die Poisson-Gleichung in den folgenden vier Kategorien entworfen und entwickelt.
--Problem --Grid-Klasse (Gitter)
- Löserklasse
- Ergebnisklasse
Das Programm, das ich dieses Mal gemacht habe, ist unter der folgenden URL.
[https://github.com/atily17/PoissonEquation]
Ich denke, es funktioniert, wenn Sie Numpy und Matplotlib haben. Ich habe nicht viel getestet, daher kann ich nicht garantieren, dass es richtig funktioniert. Wenn Sie die Problembedingungen ein wenig ändern, wird es möglicherweise sofort zusammenbrechen. Wenn Sie den Fehler beheben möchten, lassen Sie es mich bitte wissen. Ich werde es reparieren, wenn ich Lust dazu habe.
TODO Der nächste Artikel befasst sich mit der Finite-Elemente-Methode (FEM). Es gibt jedoch noch einige Bereiche, in denen FDM nicht ausreicht. Schreiben wir daher hier nur die Schlüsselwörter. Ich werde einen Artikel schreiben, wenn ich eines Tages Lust dazu habe.
- Andere Koordinatensysteme als das orthogonale Koordinatensystem
- Polarkoordinatensystem
- Gebogenes Koordinatensystem
- Wie man die Randbedingung im Unendlichen analysiert
- Boundary Relaxation Method
- Präzisions-FDM hoher Ordnung
- Formulierung der quaternären Genauigkeit --Kompakte Schablone Genauigkeit vierter Ordnung
- Nichtstrukturelles Gitter-FDM
- Subgrid
Galerie
U-förmiges Beispiel
eng
(Es fühlt sich ein wenig instabil an ...)
Ladung in der Mitte
