[PYTHON] Analyse numérique du champ électrique: méthode des différences finies (FDM) -pratique-
Préface
Poursuite de la précédente "Analyse numérique du champ électrique: méthode des différences finies (FDM) -Basics-".
Cette fois, j'écrirai enfin non seulement des formules mathématiques mais aussi une méthode aux différences finies avec python. Même ainsi, je suis un peu gêné d'écrire un programme sur le net car il semble que je ne suis pas bon dans ce domaine. Surtout en python, la propreté du script conduit directement à l'accélération et à l'optimisation, il semble donc que diverses choses seront précipitées.
En ce sens, C, C ++, Fortran, etc. sont plus rapides même si vous transformez la boucle for etc. dans le muscle cérébral, et récemment, j'apprécie celui-ci. ...... L'histoire a mal tourné. Commençons l'analyse numérique.
Le programme que j'ai réellement créé a github dans la section récapitulative, donc si vous voulez l'exécuter une fois, veuillez d'abord le télécharger et l'exécuter.
Problème de réglage
Il est difficile d'écrire un programme capable de résoudre tous les problèmes depuis le début, définissons donc le problème en premier.
Équation de Poisson
La zone de calcul est définie comme le trapèze équilatéral suivant, avec les côtés supérieur et inférieur comme condition aux limites de Diricle et les côtés gauche et droit comme condition aux limites de Neumann, et définis comme suit.
** Figure 1 **: Zone de calcul
Écrivons tout ensemble.
\begin{align}
\nabla \cdot \nabla u &= 0 \ \ \text{(Dans la zone)}\\
u &= 1 \ \ \text{(Face supérieure)} \\
u &= -1 \ \ \text{(La partie au fond)} \\
\frac{\partial u}{\partial n} &= 0 \ \ \text{(Côtés gauche et droit)}
\end{align}
Écrivons un programme qui résout ce problème. Cependant, il est ennuyeux d'écrire un programme qui ne résout que cela, alors implémentons-le tout en considérant une certaine généralité.
Les seules bibliothèques utilisées cette fois sont numpy et matplotlib.
la mise en oeuvre
Avant de dire le programme ..., écrivons d'abord un bref diagramme de classe. Même si vous dites un diagramme de classes, il n'est pas partagé par tout le monde, il est écrit dans le but d'organiser votre propre cerveau, donc il n'est pas écrit de manière bâclée, il semble qu'il existe une telle classe ~ Vous pouvez l'écrire avec un diagramme approprié Je pense (ligne de prévention).
Commençons par vérifier le flux FDM.
--Étape 1: Réglage du problème - Réglage de la zone / Réglage de l'élément source --Étape 2: Division du réseau
- Étape 3: Résoudre avec FDM - Générer des équations simultanées et les résoudre --Étape 4: Résultat
J'ai honnêtement divisé ces quatre classes en grandes classes.
Figure 2: Diagramme de classe
Laissez-moi vous expliquer brièvement. La classe PoissonEquation est la classe dite principale. Ensuite, le problème a été défini comme un problème et il a été divisé en deux parties: région Classe de domaine et terme source Classe source. Les conditions aux limites peuvent être séparées, mais j'ai décidé de les inclure dans la classe Domain.
Ensuite, dans la classe de division de grille, je l'ai mis dans la classe de grille. La grille est composée de classes (nœuds) qui gèrent les points de la grille. J'ai volontairement créé la classe Node parce que je pensais traiter le segment de ligne (Edge) et la zone (Cell) reliant chaque point de la grille, qui sera utilisé dans la méthode des éléments finis, ce que je prévois de faire la prochaine fois. ..
Ensuite, la classe Solver place également FDM sous elle, car elle place également d'autres méthodes telles que FEM en dessous. Bien qu'il ne soit pas spécifié par la ligne, doit-il s'agir d'une classe enfant? (Honnêtement, je ne sais pas parce que je ne l'ai pas encore fait)
Enfin, le résultat a été divisé en un potentiel (potentiel) et un champ électrique (densité de flux).
D'accord, nous avons une direction et écrivons un programme! Cependant, il y a de fortes chances que cela change petit à petit en le faisant.
Au fait, je suis du genre à écrire en détail depuis la fin du cours, mais est-ce courant ...? Je pense qu'il vaut mieux lire un livre comme "La manière royale d'écrire un programme". Eh bien, cette fois, je vais procéder selon ma propre méthode, donc je suis désolé si c'est une méthode étrange.
Étape 1: Configuration du problème
Commençons par définir le problème par programme. La première chose à laquelle il faut penser à ce stade est de savoir comment les enregistrer en tant que données.
Réglage de la zone (classe de domaine)
Par exemple, si l'on considère qu'elle correspond uniquement au problème trapézoïdal équilatéral mentionné ci-dessus, la forme peut être déterminée de manière unique s'il n'y a que «longueur du côté inférieur», «longueur du côté supérieur» et «hauteur», ce qui est suffisant comme données. .. Cependant, avec de telles données, il est possible de réaliser des trapèzes équilatéraux, des rectangles (côté supérieur = côté inférieur) et des carrés (côté supérieur = côté inférieur = hauteur), mais ce n'est pas compatible avec toutes les autres formes et n'est pas général. Pas intéressant. Au contraire, il est difficile de trouver trop de généralité, et il est presque impossible d'exprimer tous les sommets et courbes avec un nombre fini de données. Vous pouvez le faire en passant une fonction, mais il est évident que le programme devient assez compliqué. De cette façon, rechercher la généralité a tendance à compliquer le programme, et inversement, essayer de simplifier le programme a tendance à altérer la généralité (juste ma règle empirique, j'admets mon désaccord).
Cette fois, gardons les "coordonnées des sommets du polygone" sous forme de tableau. En faisant cela, il peut gérer tous les polygones, maintenir un certain degré de généralité et cela semble être simple comme des données. Je ne peux pas gérer les cercles, mais ... eh bien, si vous en faites un carré régulier à 20 côtés, vous pouvez l'approcher dans une certaine mesure, et vous devriez l'élargir si nécessaire. La source est la suivante.
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) #Si l'argument vertexes est une liste, np.Remplacez dans le tableau.
La raison pour laquelle j'ai créé la classe Domain et en ai fait la classe Polygon en tant que classe enfant est qu'il est possible que je crée une classe Circle afin qu'elle puisse gérer les cercles à l'avenir. Au contraire, il est également possible de créer une classe de rectangle et de la traiter à grande vitesse. Quoi qu'il en soit, je l'ai fait pour l'extensibilité. L'auteur n'est pas un professionnel, mais un amateur, donc je ne dis pas que c'est correct (ligne de prévention). Je pense que c'est une bonne idée d'écrire votre propre code.
Comme vous pouvez le voir dans le code ci-dessus, les données que nous avons ont shape, bc dans la classe parent et nVertexes, vertexes dans la petite classe.
Comme vous pouvez le voir à partir du nom et du code de la variable, elle est définie comme suit.
--shape: Le nom de la forme de la zone à utiliser. shape = "Polygone"
--bc: condition aux limites
--nVertexes: nombre de sommets de polygone
--vertices: coordonnées du sommet. Si vous les connectez dans l'ordre à partir de vertexes [0], vous pouvez créer un polygone.
Au fait, comment la condition aux limites devrait-elle être exprimée? Cette fois, je vais l'exprimer dans le type de dictionnaire suivant.
bc = {"bc":[{"bctype":"Dirichlet"Ou"Neumann" , "constant":constant}, ...], "priority":[...]}
Spécifiez la constante sur le côté droit de la condition aux limites dans «constant» et la condition aux limites Diricre / Neuman dans «bctype» ».
" bc " est constitué d'une liste de " constant " et de " bc type ". `bc [" bctype "] [0]" signifie la condition aux limites de la limite reliant les "vertexes [0]" et les "vertexes [1]".
Comme vous pouvez le voir, il ne prend en charge que les constantes ... c'est juste un problème. Je pense que d'autres que les constantes peuvent être facilement réalisées en passant le type de fonction comme argument.
Quelle condition aux limites la «priorité» donne-t-elle la priorité? Je l'ai attaché dans le sens où. Par exemple, laquelle des conditions aux limites sur les sommets de vertexes [1] est bc [" bctype "] [0] ʻou bc [" bctype "] [1]`? C'est pour préciser cela.
La "priorité" a décidé que les priorités devaient être stockées dans une liste dans l'ordre. Je n'explique pas en détail car je pense que les règles à mettre en œuvre diffèrent selon les personnes.
De plus, nous avons ajouté «calcRange», qui calcule les formes à angle droit, y compris les polygones, et «plot», qui trace pour le débogage. Il n'est pas nécessaire de l'écrire, je vais donc l'omettre.
En conséquence, le chiffre suivant a été obtenu.
Figure 3: tracé de classe de domaine
Très bien, la zone de calcul est maintenant définie.
Spécification du terme source (classe de problème)
Ensuite, créons la classe Source. Ceci est simple et peut prendre une expression (avec un tableau à deux dimensions comme argument) comme argument. La source de la classe Source elle-même est la suivante.
Source.py
import numpy as np
class Source(object):
def __init__(self, source):
self.source=source
…… Ça, ça tient dans une seule ligne. Je ne pense pas que ce genre de classement soit bon ... Faisons attention à ne pas faire ça! C'est un peu hors sujet, mais python peut changer et obtenir les variables membres de la classe directement de l'extérieur, donc je ne pense pas qu'il soit nécessaire d'avoir un getter et un setter. J'ai l'impression que python est facile à utiliser, donc je pense que c'est merveilleux. J'utilise personnellement python pour le prototype et C ++ pour la production.
Au fait, si cela est laissé tel quel, par exemple, même si le terme source est défini sur $ 0 $, la fonction anonyme lambda doit être utilisée et function = lambda x: 0 doit être donné comme argument.
C'est un peu ennuyeux, donc si le terme source est un terme constant, modifiez-le pour pouvoir passer une constante.
Le code est modifié comme suit.
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
Maintenant, lorsque l'argument source est de type int ou de type float, il semble fonctionner comme une fonction constante.
Domaine d'appel, source (classe de problème)
Écrivons brièvement la classe Problem qui rassemble les classes Domain et Source, qui sert également de résumé.
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)
Permettez-moi d'expliquer brièvement les arguments. Un exemple de l'argument est écrit sur la page suivante.
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]}}
)
#la source est une fonction constante f=Lors du réglage sur 0
setSource(
source = 0
)
#péché source(x)*sin(y)En faisant comme
setSource(
source = lambda x: np.sin(x[0])*np.sin(x[1])
)
Vous pouvez appeler ces fonctions à partir de la classe PoissonEquation.
Étape 2: Génération de grille
Ensuite, la zone est divisée en grilles. Pour être honnête, c'était le plus difficile ... Cette section n'aborde pas les algorithmes détaillés, mais je me concentrerai sur le type d'intention et le type de spécifications utilisées. Parce que [Computational Geometry](https://ja.wikipedia.org/wiki/%E8%A8%88%E7%AE%97%E5%B9%BE%E4%BD%95%E5%AD Est-ce un peu hors sujet car il est trop proche de% A6)? Parce que j'ai pensé. Bien sûr, la méthode de division de la grille est l'un des points profonds de FDM, mais elle est trop profonde, je vais donc la laisser un instant et en faire une histoire minimale. Si vous êtes intéressé, étudiez à partir de divers livres et articles.
Je l'ai fait dans les 4 étapes suivantes.
- Disposition des points de grille (nœuds) horizontalement à l’axe $ x $ et à l’axe $ y $
- Placez un nœud sur la frontière
- Effacer les nœuds en dehors de la zone
- Définition du nœud adjacent
Manipulation des points de grille (classe Node)
Étape 1. Placez le nœud
Au départ, je pensais le séparer en une classe Node et une classe Node Management, mais si vous y réfléchissez bien, python a un type de dictionnaire. Stockez les paramètres de chaque point de grille (nœud) dans un type de dictionnaire et faites de la classe Node une classe qui manipule l'ensemble.
Les paramètres du nœud sont les suivants.
node={
"point":np.array(point),
"position":"nd",
"nextnode":[
{"no":None, "position":"nd"},
{"no":None, "position":"nd"},
...
]
}
Chaque explication est résumée ci-dessous
- "point": coordonnées du nœud
- "position": Relation positionnelle des nœuds ("nd": indéfini "in": dans la zone "b (nombre n)": sur la frontière n "v (nombre n)": sur le sommet n)
- "nextnode": Liste des informations pour le nœud suivant
- "non": Index du nœud suivant
- "position": Relation de position des nœuds adjacents ("n": indéfini "l": gauche "r": droite "d": bas "u": haut)
L'argument __init __ (constructeur) cartésien nécessite une instance de domaine et le nombre de divisions (nommé le nom de la variable div).
Par exemple, si vous passez [10,10], définissons-le de manière à ce qu'il soit divisé en 10 parties égales verticalement et horizontalement. Ensuite, ça devient comme ça.
Figure 4: Grille équidistante
Vous voudrez peut-être aussi essayer quelque chose qui n'est pas également divisé (hein? Non?).
Par exemple, [[-10, -7, -3, -1,0,1,3,7,10], [-10, -7, -3, -1,0,1,3,7,10 ]] Si vous le passez comme , il le gérera bien,
Figure 5: Grille non équidistante
J'ai essayé de diviser la grille comme ça. Vous ne pouvez pas l'utiliser dans la pratique, mais c'est bon pour la recherche. Maintenant, ceci peut être réalisé comme suit.
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))]
C'est un problème, donc je ne prendrai pas la peine de l'expliquer. Les lecteurs pourront créer les leurs à leur manière. Pour le moment, pensez simplement que vous avez créé quelque chose qui crée un nœud comme celui illustré ci-dessus.
Étape 2. Placez un nœud sur la frontière
Ensuite, placez les nœuds sur la limite.
Peu importe s'il chevauche un nœud existant. Quoi qu'il en soit, placez-le là où il coupe l'horizon des axes $ x $, $ y $ du nœud existant.
J'ai fait cela en tant que getBorderPoint dans la classe Domain. Cette fois, c'est Polygon, mais je pensais que l'algorithme changeait en fonction de la forme.
Après avoir placé les nœuds sur la limite, si vous supprimez les nœuds qui se chevauchent, le résultat sera comme indiqué dans la figure suivante.
Figure 5: Bordures placées
... J'ai fait une erreur en définissant la couleur. Cela peut être difficile à voir car la bordure rouge et le nœud magenta se chevauchent, mais soyez patient. Vous pouvez voir que le deuxième nœud à partir de la gauche et le deuxième nœud à partir du haut sont très proches de la limite. Je l'ai remarqué la dernière fois que je l'ai testé, mais cela provoque beaucoup d'erreurs, il est donc raisonnable de le supprimer. À quelle distance voulez-vous effacer? Essayons de donner quelque chose comme argument.
Je l'ai écrit légèrement, mais c'était étonnamment gênant ... Si vous voulez voir la source, veuillez la voir sur github.
Étape 3. Effacez le nœud externe
Ensuite, un traitement de jugement interne est effectué pour supprimer des nœuds en dehors de la zone de calcul.
Tout d'abord, il faut juger de l'intérieur et de l'extérieur de chacun ..., mais divers algorithmes ont déjà été proposés.
Au fait, cette fois, le nœud sur la frontière est sauvegardé sur les données en tant que "position": "b" + (nombre n) dans le processus ci-dessus, donc ne modifiez que la" position ": nd. Simplement fais-le.
Ceci est très utile car le processus de jugement interne / externe à un nœud sur la frontière est assez gênant.
Il existe deux méthodes célèbres. La première consiste à tracer une demi-ligne dans n'importe quelle direction à partir du point à juger et à compter le nombre de croisements avec la limite (algorithme des nombres de croisement). Si le nombre d'intersections est pair, il est à l'extérieur, et s'il est impair, il est à l'intérieur.
Figure 6: Algorithme du nombre de croisements
Cette méthode est rapide car elle n'a que quatre règles, mais il y a certaines choses à faire attention. Comme vous pouvez le voir dans l'exemple de ③ dans la figure ci-dessus, que devez-vous faire lorsque le sommet touche le segment de ligne? Il faut y penser. En fait, il serait bon d'imposer la condition «Dans quelle direction est le vecteur de la ligne frontière qui est entrée en collision?», Mais cela semble gênant.
Alors, adoptons une autre méthode appelée Algorithme des nombres d'enroulement cette fois. Du point de vue du jugement, cela signifie que si la somme des angles entre chaque sommet est de 2π $, elle est à l'intérieur, et si elle est de 0 $, elle est à l'extérieur.
Figure 7: Algorithme du nombre d'enroulements
Puisque cette méthode utilise une fonction triangulaire inverse, cela prend un peu de temps, mais il n'est pas nécessaire de considérer les conditions ci-dessus.
Maintenant, j'écris un programme, mais il est difficile d'exécuter une boucle for sur chaque nœud.
Puisque j'utilise numpy, je parviendrai à réduire le nombre de lignes et à réduire le coût de calcul.
Le programme était le suivant.
class Polygon(Domain):
#réduction
def deleteOutDomain(self, node, thetaerr=0.1):
pts = np.array([node[i]["point"] for i in range(len(node))]) #["point"]Remplacé dans numpy
pos = np.array([node[i]["position"] for i in range(len(node))]) #["position"]Remplacé dans numpy
thetas = 0
for i in range(len(self.vertexes)):
p1 = self.vertexes[i] #Extrémité de ligne 1 p1
p2 = self.vertexes[(i + 1) % self.nVertexes] #Extrémité de ligne 2 p2
v1 = (pts - p1) #Vecteur du point que vous voulez juger à p1
v2 = (pts - p2) #Vecteur du point que vous voulez juger à p2
dot = v1[:,0] * v2[:,0] + v1[:,1] * v2[:,1] #produit intérieur(dot,Je ne visais pas l'intérieur ...)
cross = v1[:,0] * v2[:,1] - v1[:,1] * v2[:,0] #Produit extérieur(Comme ci-dessus)
vn1 = np.linalg.norm(v1, axis=1) #Obtenir la distance v1
vn2 = np.linalg.norm(v2, axis=1) #Obtenir la distance v2
theta = np.arccos(np.clip(dot / (vn1 * vn2), -1, 1)) #Clip car il peut être supérieur à 1 en raison d'une erreur numérique
thetas += np.sign(cross)*np.array(theta) #Ajoutez l'angle à chaque fois pour(la croix est un jugement plus / moins)
inx = np.where( ((pos == "nd") & ~(thetas < thetaerr)))[0] # ["position"]Est nd et thetas ne sont pas presque 0, mais obtenez l'index
#réduction
Si c'est bien fait, il peut être possible d'éliminer même la boucle vertexes, mais je pense que si vous en faites autant, cela vous pardonnera.
Le résultat est le suivant.
Figure 8: L'extérieur est effacé
4. Définissons ce qu'est le prochain nœud
Enfin, définissons node [next node]. La définition de ceci simplifie la génération de matrice pour l'étape suivante.
J'étais inquiet pour diverses choses, mais j'ai décidé de juger en haut, en bas, à gauche et à droite par des coordonnées.
La virgule flottante == (en fait np.is close) sort, donc ce n'est peut-être pas vraiment bon, mais je ne pouvais penser à rien d'autre.
Peut-être que j'aurais dû avoir un index entier comme données en premier, mais pardonnez-moi s'il vous plaît.
Et avant cela, trions la liste des nœuds par coordonnées. Je pense que c'est plus facile.
Vous n'avez pas à vous soucier de le lever, veuillez donc vous référer à github.
(En fait, quand je crée un programme, je dis "Ne pas utiliser pour!" Et c'est embarrassant car il est dupliqué ...)
Maintenant que nous avons fait cela, sortons le résultat de la liste node.
Même si le nombre de données est trop grand, il ne peut pas être visualisé de manière exhaustive, alors mettons-le autour de div = [4,4].
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'}]}
Ouais, ça fait du bien ... oh! J'ai oublié. Je voulais représenter les nœuds sur la frontière comme «f» (avant) et «b» (arrière). Puisqu'il peut être diagonal sur la frontière, il ne peut pas être exprimé uniquement par le haut, le bas, la gauche et la droite. J'ai donc revérifié le programme.
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'}]}
D'accord, c'est bon. Changeons un peu la forme et essayons-la.
Figure 9: Lorsque la forme est modifiée
Bonnes vibrations.
Étape 3: solveur FDM
Créer une matrice (générateur de matrice)
Je suis enfin là. Tout ce que j'ai à faire est de faire la queue ... mais il y a deux problèmes ici.
Expression discrète dans la zone de calcul
La première est l'équation discrète obtenue la dernière fois,
Comme le montre la figure ci-dessous, définissez l'espacement des nœuds gauche, droit, supérieur et inférieur sur $ \ Delta_l, \ Delta_r, \ Delta_u. \ Delta_d $. A ce moment, si la valeur de $ u $ à chaque point est $ u_l, u_r, u_u. U_d $, et le centre $ u_0 $ est utilisé pour l'expansion de Taylor, le résultat est le suivant.
\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}
Après cela, si vous ajoutez et soustrayez bien et effacez le terme de différenciation de premier ordre, vous pouvez le trouver comme suit.
\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}
Malheureusement, la précision est devenue la précision principale, mais elle ne peut pas être aidée. L'équation suivante peut être obtenue en la substituant dans l'équation de Poisson.
\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}
Cela semble gênant, mais un ordinateur peut le calculer immédiatement.
Une chose à noter est la taille du terme $ \ mathcal {O} (\ Delta) $. Comme vous pouvez le voir intuitivement, plus la différence entre les intervalles gauche et droit (supérieur et inférieur) est grande, plus l'erreur est grande. Lorsque vous le calculez, vous pouvez voir que si la différence entre $ \ Delta_l $ et $ \ Delta_r $ est grande, alors $ \ mathcal {O} (\ Delta_x) $ augmentera en conséquence. Ce n'est pas si difficile, et si vous êtes libre, essayez de le calculer (bien que vous puissiez le comprendre par un calcul secret).
Pour dire les choses dans l'autre sens, il est préférable que l'espacement de la grille soit le plus possible. Des choses comme la figure 5 ne sont pas vraiment bonnes. En fait, il vaut mieux augmenter le nombre de nœuds ou les diminuer.
Expression discrète sur la frontière
L'autre problème est plus gênant. Formule approximative de la condition aux limites de Neumann obtenue la dernière fois
\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}
Est juste une approximation lorsque la frontière est perpendiculaire à l'axe $ x, y $ (c'est-à-dire la différenciation partielle de $ x $ ou $ y $), et ne peut pas être utilisée dans un cas légèrement oblique comme cette fois.
Figure 10: Conditions aux limites diagonales
Ce serait bien s'il y avait un nœud dans le sens de la flèche dans la figure ci-dessus, mais ce n'est pas le cas. Revenons donc à l'essentiel et pensons au développement de Taylor.
Et avant ça
En d'autres termes, d'une manière ou d'une autre, la décentralisation suivante devrait être effectuée.
\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}
Les deuxième et troisième termes devraient être les mêmes que l'équation (1). Cela peut être fait en combinant deux formules développées par Taylor en deux points. Si vous osez l'écrire de manière redondante,
\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}
Il suffit d'ajouter ou de soustraire 3 éléments et 4 éléments à la formule (1). En d'autres termes, l'équation suivante doit être résolue pour $ c_1 $ et $ c_2 $.
\begin{align}
r_1c_1 + r_2c_2 = a \\
s_1c_1 + s_2c_2 = b
\end{align}
Si nous résolvons pour $ c_1 $ et $ c_2 $, nous pouvons l'approcher comme suit.
$ \ mathcal {E} $ représente l'erreur. Je ne savais pas comment l'écrire correctement, alors je l'ai écrit correctement.
Eh bien, c'est facile si vous pouvez le faire. Pour $ u_ {r_1, s_1} $, $ u_ {r_2, s_2} $, utilisez les deux près de la limite. Vous pouvez créer une matrice à ce sujet et la résoudre avec numpy.linalg.solve.
Cependant, il y a 3 nœuds près de la frontière, et c'est un gaspillage à résoudre avec 2 nœuds. Bien sûr, si vous pouvez résoudre avec deux, vous pouvez résoudre avec plus de précision si vous en avez trois.
En d'autres termes
Est établi. Si $ g $ est une constante, le côté droit est $ 0 $. Comme $ g = 0 $
Est. Si cette formule est faite avec trois points, une précision secondaire peut être obtenue. Pour le moment, notons les trois formules ci-dessous.
\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}
L'équation (2) a 5 termes, mais vous vous demandez peut-être si elle peut être résolue par 3 équations, mais $ k_1 $ et $ k_2 $ sont des valeurs arbitraires, donc elles peuvent être gérées. Écrivons-le comme une équation simultanée. Tout d'abord, les termes suivants sont créés en multipliant chaque terme par un coefficient et en les ajoutant.
\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}
Si vous regardez ces $ k_1 $ et $ k_2 $ comme des variables, vous aurez 5 variables et 5 équations. En le réécrivant dans une matrice, nous pouvons faire la formule suivante:
\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}
Cela semble être résolu d'une manière ou d'une autre.
Étape 4: sortie des résultats
Maintenant, enfin, le résultat est sorti.
C'est la chose la plus importante pour montrer les résultats aux gens, et c'est très amusant même si les développeurs le savent. Cependant, si de mauvais résultats sont obtenus ici, le sentiment de désespoir est grand. Non, il est peu probable que vous réussissiez la première fois ...
Potentiel
Avant de supposer le résultat, imaginons la réponse un instant. Le problème était que les côtés supérieur et inférieur du trapèze étaient respectivement 1 et -1 sous la condition aux limites de Diricre. Et si c'était un carré? C'est en fait exactement la même chose qu'un condensateur dans l'état idéal, c'est-à-dire que le champ électrique est constant de haut en bas et le potentiel est proportionnel. Ensuite, si la face supérieure devient plus petite ... Alors, j'ai écrit un diagramme en image ci-dessous.
Figure 10: Image de la solution
La flèche verte est le champ électrique et les lignes claires de couleur arc-en-ciel représentent les lignes de contour du potentiel. Si cela ressemble à ceci, j'ai pu obtenir la bonne réponse! Pour le moment, je l'ai essayé avec 100x100.
Je l'ai écrit comme si cela avait été fait en une seule fois, mais c'est le résultat d'écraser des bugs à plusieurs reprises et de les corriger avec une certaine ingéniosité.
Champ électrique (densité de flux)
Sortons également le champ électrique. Cette fois, nous allons calculer le champ électrique sur les points de la grille à partir des quatre potentiels environnants. En d'autres termes
\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}
Calculez comme suit. Le résultat est le suivant.
Ça fait du bien. Je voulais aussi le comparer avec la solution théorique, mais cette fois c'est devenu long, je vais donc l'arrêter. Si j'en ai envie un jour ...
Classe d'équation de Poisson
J'ai créé une classe PoissonEquation qui résume tout ce qui précède.
#réduction
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")
Résumé
Cette fois, nous avons conçu et développé le solveur FDM pour l'équation de Poisson dans les quatre catégories suivantes.
--Problème --Grid classe (Grille)
- Classe de solveur
- Classe de résultat
Le programme que j'ai créé cette fois est sur l'URL suivante.
[https://github.com/atily17/PoissonEquation]
Je pense que cela fonctionne si vous avez numpy et matplotlib. Je n'ai pas fait beaucoup de tests, donc je ne garantis pas que cela fonctionnera correctement. Peut-être que si vous modifiez un peu les conditions du problème, il s'effondrera immédiatement. Si vous souhaitez corriger le bogue, merci de me le faire savoir. Je vais le réparer quand j'en ai envie.
TODO Le prochain article parlera de la méthode des éléments finis (FEM). Cependant, dans certains domaines, FDM ne suffit pas, écrivons donc uniquement les mots-clés ici. J'écrirai un article quand j'en aurai envie un jour.
- Systèmes de coordonnées autres que le système de coordonnées orthogonales
-Système de coordonnées polaires
-Système de coordonnées courbé
--Comment analyser la condition aux limites à l'infini
- Boundary Relaxation Method
- FDM de haute précision --Formulation de la précision quaternaire
- Pochoir compact Précision du quatrième ordre
- Treillis non structural FDM --Sous-réseau
Galerie
Exemple en forme de U
étroit
(Cela semble un peu instable ...)
Charge au milieu
