[Calcul scientifique / technique par Python] Comparaison des vitesses de convergence de la méthode SOR, de la méthode Gauss-Seidel et de la méthode Jacobi pour l'équation de Laplace, équation différentielle partielle, problème des valeurs aux limites
introduction
La solution de l'équation différentielle partielle elliptique à l'aide de Python est décrite.
Ensuite, l'équation de Laplace a été résolue en utilisant la méthode Jacobi. On dit que la méthode Jacobi est lente à converger et peu pratique. ** Dans cet article, nous résolvons l'équation de Laplace en utilisant la méthode de Gauss-Seidel et la méthode de sur-relaxation séquentielle (méthode SOR), qui sont dites plus rapides que la méthode Jacobi, et comparons les vitesses de convergence. ** **
** Conclusion: la méthode SOR est extrêmement rapide (dans des conditions optimales) </ font> **
La méthode est brièvement résumée dans l'addendum, veuillez donc vous y référer.
Contenu
Utilisation presque du même code de calcul utilisé dans article Qiita
Comparez les vitesses de convergence de la méthode Jacobi, de la méthode Gauss-Seidel et de la méthode SOR. Plus précisément, le nombre d'itérations pour le nombre de mailles N est examiné dans chaque schéma. N a choisi quatre points: N = 100, 1000, 10000, 100000.
La méthode SOR devient la méthode Gauss-Seidel dans des cas particuliers.
SOR
%%time
"""
Équation de Laplace:Solution numérique,Méthode SOR:
ω =Gauss à 1-Devenez la méthode Seidel
Zone rectangle
14 Aug. 2017
"""
#%matplotlib nbagg
import numpy as np
import matplotlib.pyplot as plt
#
delta_Lx=0.03
delta_Ly=0.03
LLx = 10 #Largeur du rectangle(direction x)
LLy= 10 #Largeur du rectangle(direction y)
Lx = int(LLx/delta_Lx)
Ly = int(LLy/delta_Ly)
V = 2.0 #Tension
convegence_criterion = 10**-5
phi = np.zeros([Lx+1,Ly+1])
phi_Bound= np.zeros([Lx+1,Ly+1])
phi_Bound[0,:] = V #condition limite
# for SOR method
aa_recta=0.5*(np.cos(np.pi/Lx)+np.cos(np.pi/Ly)) #
omega_SOR_recta = 2/(1+np.sqrt(1-aa_recta**2)) #Paramètres d'accélération optimaux pour les zones rectangulaires
print("omega_SOR_rect=",omega_SOR_recta)
#Principale
delta = 1.0
n_iter=0
conv_check=[]
while delta > convegence_criterion:
phi_in=phi.copy()
if n_iter % 40 ==0: #Suivi de l'état de la convergence
print("iteration No =", n_iter, "delta=",delta)
conv_check.append([n_iter, delta])
for i in range(Lx+1):
for j in range(Ly+1):
if i ==0 or i ==Lx or j==0 or j==Ly:
phi[i,j] = phi_Bound[i,j]
else:
phi[i,j] = phi[i,j]+omega_SOR_recta *((phi[i+1,j] + phi[i-1,j] + phi[i,j+1] + phi[i,j-1])/4-phi[i,j]) #Mise à jour par la méthode SOR
delta = np.max(abs(phi- phi_in))
n_iter+=1
print ("The number of total iteration =",n_iter)
print("data_points=",Lx*Ly)
#for plot
plt.xlabel('X')
plt.ylabel('Y')
plt.imshow(phi,cmap='hsv')
plt.show()
Résultat (1) Affichage de la solution
Tout d'abord, la solution de l'équation de Laplace elle-même est illustrée dans la figure ci-dessous. Ceci est obtenu à partir de Article Qiita. Calculé avec un maillage 100x100.
Résultat (2) Comparaison de la vitesse de convergence
Le tableau montre la dépendance du score des données ($ Nx \ times Ny $) du nombre d'itérations nécessaires pour converger.
| Nombre de grilles: |
Méthode SOR | Gauss-Méthode Seidel | Méthode Jacobi |
|---|---|---|---|
| 100 (10x10) | 21 | 89 | 166 |
| 1089 (33x33) | 68 | 721 | 1302 |
| 10000 (100x100) | 201 | 4416 | 7476 |
| 110889 (333x333) | 663 | 22334 | 30913 |
** Vous pouvez voir que la méthode SOR est extrêmement rapide </ font>. ** ** On constate également que, comme prévu, le nombre de répétitions de la méthode [Addendum] ** SOR est proportionnel à N. La vitesse de convergence de la méthode Gauss Seidel est environ deux fois plus rapide que celle de la méthode Jacobi. Cependant, il existe toujours une différence de vitesse remarquable par rapport à la méthode SOR, et ce n'est pas pratique.
Un graphique logarithmique double de ce tableau est présenté ci-dessous. La pente montre la dépendance du score de grille du nombre d'itérations.
L'axe horizontal indique le nombre de points de grille et l'axe vertical indique le nombre d'itérations nécessaires pour que l'axe vertical converge. Le changement linéaire indique que l'indice de puissance de la dépendance du score de grille du nombre d'itérations est constant. On constate que la pente de la méthode SOR est plus petite que celle de la méthode Gauss-Seidel-Jakobi et que l'ordre de convergence est grand.
Conclusion
** Si vous connaissez les paramètres d'accélération optimaux, vous devez utiliser la méthode SOR. </ font> Un calcul réaliste ($ Nx, Ny \ gt $ 100) est plus de 20 fois plus rapide que la méthode Gauss Seidel. ** **
Addenda
Addendum (1): Méthode Gauss-Seidel: -Légère amélioration de la méthode Jacobi-
[Méthode Gauss-Seidel](https://ja.wikipedia.org/wiki/%E3%82%AC%E3%82%A6%E3%82%B9%EF%BC%9D%E3%82%B6% E3% 82% A4% E3% 83% 87% E3% 83% AB% E6% B3% 95) est Jacobi dans l'addendum de Article Qiita Formule par loi (12)
$ \phi(x_i,y_j) = \frac{1}{4} [\phi(x_{i+1},y_{j})+\phi(x_{i-1},y_{j})+\phi(x_{i},y_{j+1})+\phi(x_{i},y_{j-1})]+\frac{\rho(x_{i},y_{j})}{4\epsilon_0}\Delta^2 $
$ (N + 1) $ Représentation par étapes
$ \phi(x_i,y_j)^{(n+1)} = \frac{1}{4} [\phi(x_{i+1},y_{j})^{(n)}+\phi(x_{i-1},y_{j})^{(n)}+\phi(x_{i},y_{j+1})^{(n)}+\phi(x_{i},y_{j-1})^{(n)}]+\frac{\rho(x_{i},y_{j})}{4\epsilon_0}\Delta^2 \tag{1} $
Dans, $ \ phi (x_ {i-1}, y_j) ^ {(n)} $ et $ \ phi (x_ {i}, y_ {j-1}) ^ {(n)} $ sont ** en fait Puisqu'il a été mis à jour à l'avance, la valeur de l'étape $ n $ n'est pas utilisée. Utilisez $ \ phi (x_ {i-1}, y_j) ^ {(n + 1)} $ et $ \ phi (x_ {i}, y_ {j-1}) ^ {(n + 1)} $ Méthode **. En d'autres termes
$ \phi(x_i,y_j)^{(n+1)} = \frac{1}{4} [\phi(x_{i+1},y_{j})^{(n)}+\phi(x_{i-1},y_{j})^{(n+1)}+\phi(x_{i},y_{j+1})^{(n)}+\phi(x_{i},y_{j-1})^{(n+1)}]+\frac{\rho(x_{i},y_{j})}{4\epsilon_0}\Delta^2 \tag{2} $ Et. C'est ** parce que $ \ phi (x_i, y_j) $ est mis à jour séquentiellement, de sorte que l'utilisation de la mémoire du code peut être sauvegardée par rapport à la méthode Jacobi ** On peut s'attendre à ce que la convergence soit légèrement plus rapide que la méthode Jakobi (le nombre d'itérations pour un grand maillage sera environ la moitié de celui de la méthode Jakobi [1]), mais cette méthode a également une convergence plus lente.
Addendum (2): Solution numérique: Méthode SOR: -Légère amélioration de la méthode Jacobi-
Un des algorithmes pratiques est la ** méthode de sur-relaxation séquentielle (méthode SOR [2]) **. Comme indiqué ci-dessous, cette méthode accélère la mise à jour de $ \ phi (x_i, y_j) ^ {(n)} $ par la méthode ** Gauss-Seidel en appliquant le paramètre d'accélération $ \ omega $. ** **
Le processus itératif de la méthode Gauss-Seidel peut être réécrit comme suit [1].
Ici, $ \ Delta \ phi (x, y) $ représente le montant de la mise à jour de $ \ phi (x, y) ^ {(n)} $ par la méthode Gauss-Seidel. Au fait, ** Dans la méthode Gauss-Seidel, $ \ phi (x, y) ^ {(n)} $ change souvent de façon monotone dans le processus itératif **. A titre d'exemple, la figure montre les changements dans le processus itératif de certaines valeurs $ \ phi $ lors de la résolution de l'équation de Laplace sur un maillage 10x10. Vous pouvez voir que cela change de façon monotone.
Changements dans $ \ phi $ pendant le processus itératif.
** La méthode SOR est une méthode pour utiliser activement ce changement monotone. Puisque $ \ phi $ change de manière monotone, le montant de la mise à jour $ \ Delta \ phi (x, y) $ (équation (4)) prédit par la méthode Gauss-Seidel peut être augmenté pour augmenter la vitesse de convergence. C'est vrai. ** Cette idée est l'essence de la méthode SOR.
Par conséquent, processus itératif de multiplication de $ \ Delta \ phi (x, y) $ dans l'équation (4) par $ \ omega_ {SOR} $
$ \phi(x_i,y_j)^{(n+1)} = \phi(x_i,y_j)^{(n)}+\omega_{SOR}[\frac{\phi(x_{i+1},y_{j})^{(n)}+\phi(x_{i-1},y_{j})^{(n)}+\phi(x_{i},y_{j+1})^{(n)}+\phi(x_{i},y_{j-1})^{(n)}}{4} +\frac{\rho(x_{i},y_{j})}{4\epsilon_0}\Delta^2-\phi(x_i,y_j)^{(n)}]\tag{5} $
penser à.
** Il s'agit de la méthode SOR. $ \ Omega_ {SOR} $ est appelé le paramètre de sur-relaxation. Lorsque $ \ omega_ {SOR} $ = 1, cela devient la méthode Gauss-Seidel. ** **
Il existe une option pour donner $ \ omega_ {SOR} $. Le choix de la méthode optimale (meilleur taux de convergence) peut entraîner une convergence beaucoup plus rapide que la méthode de Gauss-Seidel. Différentes façons de trouver le meilleur $ \ omega_ {SOR} $ ont été étudiées [2].
** Nous connaissons le $ \ omega_ {SOR} $ optimal pour l'équation de Laplace (Poisson) dans la région rectangulaire, **
Est. Pour une grille suffisamment grande $ Nx \ times Ny >> 1 $, cela devient $ \ omega_ {SOR} \ approx 2 $.
Disons ** $ Nx = Ny = N $. À ce stade, le nombre d'itérations nécessaires à la convergence du calcul à l'aide de la méthode SOR est proportionnel à $ N ^ 1 $. Celle de la méthode Jacobi et de la méthode Gauss Seidel est proportionnelle à $ N ^ 2 $. Lorsque N est suffisamment grand, la vitesse de convergence de la méthode SOR est nettement plus rapide que celle de la méthode Gauss-Seidel et de la méthode Jacobi [2] </ font> **. Les résultats examinés dans cet article le soutiennent fortement.
Notez que ** facile à mettre en œuvre **. Utilisez $ \ omega_ {SOR} $ pour garder l'équation (5) telle quelle
phi[i,j] = phi[i,j]+omega_SOR*((phi[i+1,j] + phi[i-1,j] + phi[i,j+1] + phi[i,j-1])/4-phi[i,j])
Tout ce que vous avez à faire est d'écrire [2].
Les références
[1] Mon article Qiita, [[Calcul scientifique / technique par Python] Solution numérique de l'équation de Laplace-Poisson bidimensionnelle par la méthode Jacobi pour la position électrostatique, équation différentielle partielle elliptique, problème de valeur aux limites](http: // qiita. com / sci_Haru / items / 6b80c7eb8d4754fb5c2d)
[2] ["Recette numérique en version japonaise de la mer-Recette pour le calcul numérique en langage C"](https://www.amazon.co.jp/%E3%83%8B%E3%83%A5 % E3% 83% BC% E3% 83% A1% E3% 83% AA% E3% 82% AB% E3% 83% AB% E3% 83% AC% E3% 82% B7% E3% 83% 94-% E3% 82% A4% E3% 83% B3-% E3% 82% B7% E3% 83% BC-% E6% 97% A5% E6% 9C% AC% E8% AA% 9E% E7% 89% 88- C% E8% A8% 80% E8% AA% 9E% E3% 81% AB% E3% 82% 88% E3% 82% 8B% E6% 95% B0% E5% 80% A4% E8% A8% 88% E7% AE% 97% E3% 81% AE% E3% 83% AC% E3% 82% B7% E3% 83% 94-Press-William-H / dp / 4874085601 / ref = dp_return_1? _ Encoding = UTF8 & n = 465392 & s = livres) Revue technique, 1993.
Recommended Posts