[Python] Simulation de fluides: implémenter l'équation de diffusion

introduction

Je voudrais résumer (en plusieurs articles) les connaissances nécessaires pour construire un code d'analyse numérique des fluides pour l'eau, tout en étudiant également la dynamique numérique des fluides (CFD), qui est une étude liée à la simulation de fluides tels que l'air et l'eau.

Je pense que la dynamique des fluides numérique est une étude assez difficile, alors j'aimerais l'écrire d'une manière facile à comprendre pour les débutants. Il semble qu'il y ait beaucoup d'erreurs, alors veuillez nous contacter si vous le trouvez. Je vous serais également reconnaissant de bien vouloir commenter ce qui est difficile à comprendre. Nous le mettrons à jour de temps en temps.

Public cible

or

séries

Contenu approximatif de cet article

Comme étape préliminaire pour traiter les équations de base du fluide requises pour la simulation de l'eau, nous allons également résumer brièvement et mettre en œuvre l'équation de diffusion. ** Je l'ai rendu compréhensible sans lire l'article précédent (probablement) **

  1. [Différence centrale](différence # 1-Center)
  2. [Méthode implicite de Crank Nicholson](# 2-Méthode implicite de Crank Nicholson)
  3. [Méthode directe](# 2-1-Méthode directe)
  4. [Méthode de répétition](# 2-2-Méthode de répétition)
  5. [Méthode itérative constante](# 2-2-1-Méthode de répétition constante)
  6. [Méthode de répétition non stationnaire](# 2-2-2-Méthode de répétition non stationnaire Méthode du sous-espace krylov)
  7. [Implémentation de la méthode d'itération stationnaire](# 2-3-Implémentation de la méthode d'itération stationnaire)
  8. [Méthode Jacobi](méthode # 2-3-1-jacobi)
  9. [Exemple simple](# 2-3-1-1-Exemple simple)
  10. [Attention](précaution de la méthode # 2-3-1-2-jacobi)
  11. [Implémentation](implémentation de la méthode # 2-3-1-3-jacobi)
  12. [Méthode Gauss-Seidel](méthode # 2-3-2-gauss-seidel)
  13. [Attention](précaution de la méthode # 2-3-2-1-gauss-seidel)
  14. [Mise en œuvre](mise en œuvre de la méthode # 2-3-2-2-gauss-seidel)
  15. [Méthode SOR](méthode # 2-3-3-sor)
  16. [Mise en œuvre](mise en œuvre de la méthode # 2-3-3-1-sor)

Équation de diffusion (formule pour l'uniformisation du champ)

Qu'est-ce qu'une équation de diffusion unidimensionnelle? $ \frac{\partial T}{\partial t} = \kappa \frac{\partial^2 T}{\partial x^2} $

Il est exprimé par la formule. En termes de signification physique, cela signifie que les quantités physiques sont uniformisées de manière aléatoire. Cela signifie la conduction thermique et la diffusion de substances ainsi que la chaleur. Par exemple, lorsque du lait est ajouté au café, il se répand lentement sans remuer, ou lorsque l'extrémité d'une tige métallique est placée dans une marmite bouillante comme le montre la figure ci-dessous, la tige entière devient progressivement de l'eau chaude. Le phénomène de la même température correspond à la diffusion.

diffusion_phenomenon.png

Les méthodes pour discriminer cette équation de diffusion unidimensionnelle comprennent la méthode de différence centrale explicite et la méthode implicite de la méthode implicite de Crank Nicholson. Pour rappel, la méthode explicite prédit la valeur future une heure plus tard en fonction de la valeur connue de l'heure actuelle, et la méthode implicite prédit également la valeur future une heure plus tard. Nous allons mettre en œuvre ces deux. Si vous écrivez chacun avec une formule de différence, ce sera comme suit.

  1. Différence de centre Il peut être dérivé en effectuant deux fois la différence de précision au vent du premier ordre expliquée dans l'article précédent. Il s'agit d'un type de méthode explicite qui prédit la valeur de la prochaine fois n + 1 en utilisant uniquement les données à l'instant courant n. $ \frac{T_j^{n+1} - T_j^n}{\Delta t} = \kappa \frac{T_{j+1}^n - 2 T_j^n + T_{j-1}^n }{\Delta x^2} $

  2. Méthode implicite de Crank Nicholson

Méthode de moyennage d'une valeur connue (temps n) et d'une valeur inconnue un pas plus tard (temps n + 1) pour la différenciation partielle liée à la différenciation spatiale. Comme son nom l'indique, il s'agit d'une méthode implicite qui utilise également la valeur de la prochaine fois n + 1 pour prédire la valeur de la prochaine fois n + 1. Les détails seront décrits plus tard. $ \frac{T_j^{n+1} - T_j^n}{\Delta t} = \frac{\kappa}{2} \left(\frac{T_{j+1}^n - 2 T_j^n + T_{j-1}^n }{\Delta x^2} + \frac{T_{j+1}^{n+1} - 2 T_j^{n+1} + T_{j-1}^{n+1} }{\Delta x^2}\right) $

Détermination de la stabilité par le nombre de spreads

Parmi les méthodes de dispersion présentées ci-dessus, la ** méthode explicite utilisant la différence centrale ** $ d = \kappa \frac{\Delta t}{\Delta x^2} $ Il est nécessaire de déterminer la stabilité du calcul numérique en utilisant ce qu'on appelle le ** nombre de diffusion d **. En revanche, la méthode implicite est inconditionnellement stable.

En particulier,

d \leq \frac{1}{2} \quad or \quad \Delta t \leq \frac{(\Delta x)^2}{2 \kappa}

Doit être compris entre. L'important dans cette ** condition de nombre de diffusion est que si vous voulez réduire la taille du pas $ \ Delta x $, vous devez réduire la taille du pas de temps $ \ Delta t $ de son carré **. Par exemple, si vous multipliez $ \ Delta x $ par 1/2, vous devez réduire $ \ Delta t $ à 1/4. En conséquence, la résolution de l'équation de diffusion par la méthode explicite tend à augmenter la charge de calcul, donc en général, la méthode implicite est souvent utilisée.

De plus, cette condition est obtenue à partir de l'analyse de stabilité de Von Neumann, mais je ne vais pas l'expliquer ici, mais seulement une image intuitive.

Si vous réécrivez la formule de la différence centrale,

T_j^{n+1} = d T_{j+1}^n + (1 - 2d) T_j^n + d T_{j-1}^n

Ce sera. La diffusion signifie qu'il est uniformisé de manière désordonnée, donc lorsque le coefficient $ (1-2d) $ de $ T_j ^ n $ devient négatif, il passe plus que la quantité physique du jème réseau à côté. Par conséquent, ce devrait être $ d \ leq0.5 $.

Je pense que c'est difficile à comprendre, alors je vais l'expliquer avec l'exemple de l'argent ci-dessous (** C'est juste un exemple pour avoir une image ). Supposons que vous ayez un groupe de N personnes et qu'à un certain point n, la j-1ère personne a 30 yens, la jème personne 100 yens et la j + 1ère personne 50 yens. Le spread est une opération qui permet à ces N personnes d'avoir le même montant d'argent ( uniformisation **), alors continuez à échanger un certain pourcentage d'argent avec les voisins jusqu'à ce que le même montant soit atteint. Je vais. Si vous regardez attentivement la formule de la différence de centre, ce pourcentage constant est le nombre de diffusion d. Par conséquent, si vous définissez d = 0,1 et calculez, la jième personne aura 88 yens à n + 1 fois une heure plus tard. En pensant de cette façon, vous pouvez facilement voir que le nombre de diffusion d ne doit pas être supérieur à 1/2.

money_diffusion_r2.png ** * Ceci est un exemple pour saisir l'atmosphère de diffusion **

Condition de calcul

problem.png

answer.png

La largeur du réseau $ \ Delta x $, le coefficient de conduction thermique $ \ kappa $ et la largeur du pas de temps $ \ Delta t $ sont tous calculés comme des constantes. Cette fois, $ \ Delta x = 1 $, $ \ Delta t = 0,2 $, $ \ kappa = 0,5 $, et le calcul est effectué dans des conditions stables avec le nombre de diffusion $ d = 0,1 $.

import numpy as np
import matplotlib.pyplot as plt
Num_stencil_x = 101
x_array = np.float64(np.arange(Num_stencil_x))
temperature_array = x_array + 100
temperature_lower_boundary = 150
temperature_upper_boundary = 150
Time_step = 100
Delta_x = max(x_array) / (Num_stencil_x-1)
Delta_t = 0.2
kappa = 0.5
d = kappa * Delta_t / Delta_x**2
exact_temperature_array = (temperature_upper_boundary - temperature_lower_boundary) / (x_array[-1] - x_array[0]) * x_array + temperature_lower_boundary
plt.plot(x_array, temperature_array, label="Initial condition")
plt.plot(x_array, exact_temperature_array, label="Answer", c="black")
plt.legend(loc="lower right")
plt.xlabel("x")
plt.ylabel("Temperature [K]")
plt.xlim(0, max(x_array))
plt.ylim(100, 200)

la mise en oeuvre

1. Différence de centre

\frac{T_j^{n+1} - T_j^n}{\Delta t} = \kappa \frac{T_{j+1}^n - 2 T_j^n + T_{j-1}^n }{\Delta x^2}

La figure ci-dessous montre le résultat lorsque le nombre de diffusion d est de 0,1. J'ai essayé de sortir le graphique en le divisant en plusieurs pas de temps (temps). Après 100 heures, les deux extrémités ont atteint près de 150K, et la distribution de température est entraînée par elle comme une courbe sinusoïdale. Au bout de 5000 heures, la température de l'ensemble de la tige s'approche progressivement de 150K, indiquant que la mise en œuvre de la différence de centre est réussie. Puisque la mise en œuvre est facile à suivre l'expression, je l'écrirai sans utiliser autant que possible la fonction numpy.

explicit.png

Au fait, si vous définissez $ \ kappa = 5 $ et calculez avec le nombre de diffusion $ d = 1 $, vous pouvez voir que le calcul diverge et que la solution ne peut pas être obtenue comme le montre la figure ci-dessous.

explicit_diffuse.png


temperature_explicit = temperature_array.copy()
for n in range(Time_step):
    temperature_old = temperature_explicit.copy()
    temperature_explicit[0] += kappa * Delta_t / Delta_x**2 * \
        (temperature_explicit[1] - 2*temperature_old[0] + temperature_lower_boundary)
    temperature_explicit[-1] += kappa * Delta_t / Delta_x**2 * \
        (temperature_upper_boundary - 2*temperature_old[-1] + temperature_old[-2])
    for j in range(1, Num_stencil_x-1):
        temperature_explicit[j] += kappa * Delta_t / Delta_x**2 * \
            (temperature_old[j+1] - 2*temperature_old[j] + temperature_old[j-1])
plt.plot(x_array, exact_temperature_array, label="Answer", c="black")
plt.plot(x_array, temperature_explicit, label="Explicit(100 steps)")
plt.legend(loc="lower right")
plt.xlabel("x")
plt.ylabel("temperature")
plt.xlim(0, max(x_array))
plt.ylim(100, 200)

2. Méthode implicite de Crank Nicholson

\frac{T_j^{n+1} - T_j^n}{\Delta t} = \frac{\kappa}{2} \left(\frac{T_{j+1}^n - 2 T_j^n + T_{j-1}^n }{\Delta x^2} + \frac{T_{j+1}^{n+1} - 2 T_j^{n+1} + T_{j-1}^{n+1} }{\Delta x^2}\right)

Si vous transformez cela et amenez la valeur du temps n + 1 sur le côté gauche $ -T_{j+1}^{n+1} +2 \left(\frac{1}{d} + 1 \right) T_j^{n+1} - T_{j-1}^{n+1} = T_{j+1}^{n} +2 \left(\frac{1}{d} - 1 \right) T_j^{n} + T_{j-1}^{n} \\ d = \kappa \frac{\Delta t}{\Delta x^2} $

En supposant que les points de grille existent dans la plage de 1 à M, la valeur limite est représentée par $ T_0, T_ {M + 1} $, et elle est convertie en un affichage matriciel.


\left(
\begin{array}{cccc}
      2 \left(\frac{1}{d} + 1 \right) & -1 &  0 & \ldots  & \ldots & 0 \\
      -1 & 2 \left(\frac{1}{d} + 1 \right) & -1 & 0 & \ldots & 0 \\
      0  &-1 & 2 \left(\frac{1}{d} + 1 \right) & -1 & 0 & \ldots  \\
      \vdots & \ddots & \ddots & \ddots & \ddots & \ddots\\
      0  & 0 & \ddots & -1 & 2 \left(\frac{1}{d} + 1 \right) & -1 \\
      0  & 0 & \ldots & 0 & -1 & 2 \left(\frac{1}{d} + 1 \right)
    \end{array}
  \right)
  \left(
    \begin{array}{c}
      T_1^{n+1}  \\
      T_2^{n+1}  \\
      T_3^{n+1}    \\
      \vdots \\
      T_{M-1}^{n+1} \\
      T_M^{n+1} 
\end{array}
\right)


= \left( \begin{array}{c}
    T_2^{n} + 2 \left(\frac{1}{d} - 1 \right) T_1^{n} + \left(T_0^n + T_0^{n+1}\right) \\
      T_3^{n} + 2 \left(\frac{1}{d} - 1 \right) T_2^{n} + T_1^n  \\
      T_4^{n} + 2 \left(\frac{1}{d} - 1 \right) T_3^{n} + T_2^n  \\
      \vdots \\
      T_M^{n} + 2 \left(\frac{1}{d} - 1 \right) T_{M-1}^{n} + T_{M-2}^n  \\
      \left(T_{M+1}^n + T_{M+1}^{n+1}\right) + 2 \left(\frac{1}{d} - 1 \right) T_{M}^{n} + T_{M-1}^n
    \end{array} \right)

Envisagez de résoudre cette équation unidimensionnelle simultanée. Lorsque le nombre de points de grille M est petit, le nombre inconnu est petit, utilisez donc une méthode directe qui trouve directement une solution, comme une méthode appelée méthode directe de Gauss, qui calcule comme un calcul de pinceau. Cependant, plus M est grand, plus il est difficile de trouver une solution en termes de calcul. Par conséquent, en général, afin de résoudre une équation unidimensionnelle simultanée à grande échelle, une méthode appelée ** méthode itérative ** qui fait converger la solution approchée vers une solution exacte est utilisée.

2-1. Méthode directe

Les types sont les suivants. J'expliquerai quand c'est nécessaire.

2-2. Méthode itérative

Ax=b

Envisagez de résoudre une équation simultanée unidimensionnelle. cette $ r=b-Ax' $ La méthode pour trouver une solution approximative de la solution exacte $ x ^ * = A ^ {-1} b $ est répétée en modifiant la valeur estimée $ x '$ de manière itérative jusqu'à ce que r devienne suffisamment petit. C'est la loi. La méthode itérative peut être divisée en méthode stationnaire et méthode non stationnaire (méthode du sous-espace de Krylov).

2-2-1. Méthode de répétition constante

La méthode stationnaire (méthode itérative stationnaire) est une méthode qui ne change pas sauf pour le vecteur solution x lors du calcul itératif. Comme il est généralement lent, utilisez la méthode non stationnaire décrite plus loin lors de l'exécution de calculs numériques, ou utilisez-la ensemble comme prétraitement de la méthode non stationnaire (traitement effectué pour trouver une solution approximative approximative). pense.

Les trois méthodes suivantes peuvent être mentionnées comme méthodes typiques. Dans ce qui suit, la valeur estimée de la solution en m itérations est exprimée comme $ x ^ {(m)} $, et la valeur estimée $ x ^ {(m)} $ dans la mième itération est connue sous le nom de m + 1e itération. Voici un exemple de la façon de trouver la valeur estimée de $ x ^ {(m + 1)} $.

  1. Méthode Jacobi

Une méthode de calcul de la solution estimée $ x_i ^ {(m + 1)} $ sur la i-ième ligne de l'itération m + 1 en utilisant uniquement la solution estimée connue $ x_i ^ {(m)} . $ x_i^{(m+1)} = \frac{1}{a_{ii}}\left( b_i - \sum_{j=1,j\not=i}^{n} a_{ij} x_j^{(m)} \right) $$ 2. Méthode Gauss-Seidel

Lors de la recherche de la solution estimée $ x_i ^ {(m + 1)} $ sur la i-ième ligne de l'itération m + 1, la solution estimée connue $ x_i ^ {(m)} $ est déjà calculée comme m + 1 Solution estimée pour les itératifs $ x_ {0} ^ {(m + 1)}, \ cdots, x_ {i-1} ^ {(m + 1)} . Fondamentalement, le calcul de convergence se termine plus rapidement que la méthode Jacobi. $ x_i^{(m+1)} = \frac{1}{a_{ii}}\left( b_i - \sum_{j=1}^{i-1} a_{ij} x_j^{(m+1)} - \sum_{j=i+1}^{n} a_{ij} x_j^{(m)} \right) $$ 3. Méthode SOR (Successive Over-Relaxation) / Méthode de relaxation séquentielle d'accélération

Méthode de Gauss-Seidel plus facteur de relaxation $ \ omega $. Lorsque le coefficient de relaxation $ \ omega = 1 , la méthode de Gauss-Seidel est utilisée. De plus, si elle est égale ou inférieure à 1, la convergence est plus lente que Gauss-Seidel, mais les problèmes qui ne peuvent pas être résolus par la méthode de Gauss-Seidel peuvent être résolus. En gros, définissez-le sur une valeur de 1 ou plus, mais s'il est trop grand, il divergera, il est donc important de sélectionner une valeur appropriée. Il semble que 1,9 soit souvent sélectionné. $ x_i^{(m+1)} = x_i^{(m)} + \omega \frac{1}{a_{ii}}\left( b_i - \sum_{j=1}^{i-1} a_{ij} x_j^{(m+1)} - \sum_{j=i}^{n} a_{ij} x_j^{(m)} \right) $$

  1. Méthode multi-grille Il existe la méthode multi-grille géométrique et la méthode multi-grille algébrique (AMG). Cette dernière méthode AMG est souvent utilisée comme prétraitement pour les méthodes non stationnaires. Je pense que beaucoup de gens le sauront pour la première fois, alors j'écrirai les détails dans un autre article.

Les détails des trois premiers éléments seront décrits plus loin.

2-2-2. Méthode de répétition non stationnaire (méthode du sous-espace de Krylov)

Est-ce comme un bon garçon avec des méthodes directes et itératives? Elle est classée comme méthode itérative.

La méthode non stationnaire la plus connue est la méthode du gradient conjugué (méthode CG). Des explications détaillées sont données dans cet article. Il y a un chiffre et c'est incroyablement facile à comprendre. En outre, cet article explique la différence par rapport à la méthode de descente la plus raide souvent utilisée dans l'apprentissage automatique. L'image de la méthode CG est également écrite et elle est facile à comprendre.

  1. Matrice symétrique
  2. Méthode du gradient conjugué (méthode CG)
  3. Méthode du gradient conjugué préconditionné (méthode PCG)
  4. Méthode incomplète du gradient conjugué de Cholesky (méthode ICCG)
  5. Méthode de gradient conjugué LU incomplète (ILUCG)
  6. Autre
  7. Matrice asymétrique
  8. Gradient bi-conjugué (méthode BiCG)
  9. Conjugate Residual Squared (méthode CGS)
  10. Stabilisation BiCG (méthode BiCG STAB)
  11. Résidu conjugué généralisé (méthode GCR)
  12. Résiduel minimal généralisé (méthode GMRES)

Ces techniques sont essentiellement fournies sous forme de bibliothèque. De plus, je pense qu'il y a beaucoup de codes de calcul sur Github. Si vous souhaitez étudier plus en détail, veuillez vous y référer également. Dans un autre article, je voudrais expliquer et mettre en œuvre autant de méthodes que possible pour la méthode non stationnaire.

Pour le moment, dans cet article, j'aimerais utiliser la méthode Jacobi, la méthode Gauss-Seidel et la méthode SOR, qui sont les méthodes d'itération stationnaire, parmi les algorithmes expliqués jusqu'à présent.

2-3. Mise en œuvre de la méthode d'itération stationnaire

Nous l'implémenterons en référence au Material donné par le professeur Nakajima de l'Université de Tokyo. Si vous trouvez l'explication difficile à comprendre, nous vous serions reconnaissants de bien vouloir vous y référer.

Dans l'exemple ci-dessous

A = \left(
    \begin{array}{ccc}
      a_{1,1} & a_{1,2} & \ldots &a_{1,n}  \\
      a_{2,1} & a_{2,2} & \ldots &a_{2,n} \\
      \vdots & \vdots &  \ddots &\vdots \\
       a_{n,1} & \ldots &  \ldots &a_{n,n}
    \end{array}
  \right) \\
D = \left(
    \begin{array}{ccc}
      a_{1,1} & 0 & \ldots &0  \\
      0 & a_{2,2} & \ldots &0 \\
      \vdots & \vdots &  \ddots &\vdots \\
       0 & \ldots &  \ldots &a_{n,n}
    \end{array}
  \right) \\
L = \left(
    \begin{array}{ccc}
      0 & 0 & \ldots &0  \\
      a_{2,1} & 0 & \ldots &0 \\
      \vdots & \vdots &  \ddots &\vdots \\
       a_{n,1} & \ldots &  \ldots &0
    \end{array}
  \right) \\
U = \left(
    \begin{array}{ccc}
      0 & a_{1,2} & \ldots &a_{1,n}  \\
      0 & 0 & \ldots &a_{2,n} \\
      \vdots & \vdots &  \ddots &\vdots \\
       0 & \ldots &  \ldots &0
    \end{array}
  \right)

Définissez la matrice comme suit. D est uniquement le composant diagonal de A, L est le côté inférieur du composant diagonal et U est le côté supérieur du composant diagonal.

2-3-1. Méthode Jacobi

Lorsque $ A = D + L + U $ en utilisant la matrice D ayant des composantes diagonales de la matrice A, de la matrice supérieure U et de la matrice inférieure L,

　Ax = b\\
 (D+L+U) x = b\\
 Dx = b - (L+U)x \\
 x = D^{-1}(b-(L+U)x)

$ Ax = b $ peut être transformé comme ceci, et la méthode de recherche répétée de la solution approximative avec la formule en bas est appelée la méthode Jacobi.

x_i^{(m+1)} = \frac{1}{a_{ii}}\left( b_i - \sum_{j=1,j\not=i}^{n} a_{ij} x_j^{(m)} \right)

2-3-1-1. Exemple simple

Soit A une matrice $ 3 \ times3 $


\left(
    \begin{array}{ccc}
      a_{1,1} & a_{1,2} &  a_{1,3}  \\
      a_{2,1} & a_{2,2} &  a_{2,3} \\
      a_{3,1} & a_{3,2} &  a_{3,3}
    \end{array}
  \right)
  \left(
    \begin{array}{c}
      x_{1} \\
      x_{2} \\
      x_{3} 
    \end{array}
  \right)
  =
  \left(
    \begin{array}{c}
      b_{1} \\
      b_{2} \\
      b_{3} 
    \end{array}
  \right)

Envisagez de résoudre une équation simultanée unidimensionnelle.

En exprimant le nombre d'itérations comme (m), la valeur de $ x ^ {(m + 1)} $ à (m + 1) après une itération est

x_1^{(m+1)} = \left(b_1 - a_{1,2} x_2^{(m)} - a_{1,3} x_3^{(m)} \right) / a_{1,1} \\
x_2^{(m+1)} = \left(b_2 - a_{2,1} x_1^{(m)} - a_{2,3} x_3^{(m)} \right) / a_{2,2} \\
x_3^{(m+1)} = \left(b_3 - a_{3,1} x_1^{(m)} - a_{3,2} x_2^{(m)} \right) / a_{3,3}

Il peut être exprimé sous la forme de. Comme vous pouvez le voir sous la forme de cette équation, $ x ^ {(m + 1)} = D ^ {-1} (b- (L + R) x ^ {(m)}) $ tient. Je vais.

Ensuite, en répétant le calcul ci-dessus jusqu'à ce que la condition de convergence ci-dessous soit satisfaite, la solution exacte $ x ^ * $ peut être obtenue.

\sum_{j=1}^{n} \left| \frac{x_j^{(m+1)}-x_j^{(m)}}{x_j^{(m+1)}} \right|= \epsilon

Cependant, $ \ epsilon $ doit être un nombre positif suffisamment petit.

Maintenant, pratiquons la méthode Jacobi avec un problème simple.

À titre d'exemple simple

\left(
    \begin{array}{ccc}
      3 & 2 &  -0.5  \\
      1 & 4 &  1  \\
      -1 & 0 &  4  
    \end{array}
  \right)
  \left(
    \begin{array}{c}
      x_1   \\
      x_2  \\
      x_3
    \end{array}
  \right)
  =
  \left(
    \begin{array}{c}
      3   \\
      2  \\
      1
    \end{array}
  \right)

Pensez à résoudre. La réponse est

\left(
    \begin{array}{c}
      1   \\
      0.125  \\
      0.5
    \end{array}
  \right)

Contrairement à avant, nous allons l'implémenter en utilisant beaucoup de fonctions numpy. Parce qu'il est extrêmement facile à voir. La mise en œuvre est la suivante. Lorsque vous calculez réellement

Solution [1.00000063 0.12499957 0.50000029]
Answer [1.    0.125 0.5  ]

S'affiche et il peut être confirmé que la réponse est recherchée.

def jacobi(a_matrix, b_array, target_residual):
    """
Méthode Jacobi
    Ax = b
    A = (D+L+U)Puis
    Dx = b - (L+U)x
    x = D^{-1} (b - (L+U)x)
Cependant, D est une matrice diagonale, L+Soit U la matrice résiduelle.
    
    Parameters
    ----------
    a_matrix : numpy.float64
matrice m × n
    b_array : numpy.float64
matrice de m lignes
    target_residual : numpy.float64
Nombre positif. Cibler les résidus.

    Returns
    -------
    x : numpy.float64
matrice de m lignes
    """
    x = b_array.copy()
    x_old = b_array.copy()
    
    diag_matrix= np.diag(a_matrix)  #Matrice diagonale
    l_u_matrix = a_matrix - np.diagflat(diag_matrix)  #Matrice résiduelle
    count = 0
    while True:
        count += 1
        x = (b_array - np.dot(l_u_matrix, x_old))/diag_matrix
        residual = np.linalg.norm(x - x_old) / np.linalg.norm(x)
        x_old = x
        if residual <= target_residual:
            break
        elif count >= 10000:
            import sys
            print(residual)
            sys.exit()
    return x

A = np.array([[ 3.0,  2.0, -0.5],
              [ 1.0,  4.0,  1.0],
              [-1.0,  0.0,  4.0]])
b = np.array([3.0, 2.0, 1.0])

x = jacobi(A, b, 10**-6)
print("Solution", x)
print("Answer", np.dot(np.linalg.inv(A),b))

2-3-1-2. Précautions relatives à la loi Jacobi

Si la matrice A ne satisfait pas la dominance diagonale (au sens étroit), la méthode de Jacobi (et la méthode de Gauss-Seidel) ne convergeront pas. Pour énoncer brièvement ce qu'est l'avantage diagonal (au sens étroit), il dit que la composante diagonale est supérieure à la somme des valeurs absolues des autres composantes de la même ligne. Dans l'exercice ci-dessus, essayez de jouer avec les composants non diagonaux plus grands.

\left| a_{ii} \right| > \sum_{j=1, j\not=i}^{n} \left| a_{ij} \right|

2-3-1-3. Implémentation de la méthode Jacobi

Je republierai la formule à résoudre.


\left(
\begin{array}{cccc}
      2 \left(\frac{1}{d} + 1 \right) & -1 &  0 & \ldots  & \ldots & 0 \\
      -1 & 2 \left(\frac{1}{d} + 1 \right) & -1 & 0 & \ldots & 0 \\
      0  &-1 & 2 \left(\frac{1}{d} + 1 \right) & -1 & 0 & \ldots  \\
      \vdots & \ddots & \ddots & \ddots & \ddots & \ddots\\
      0  & 0 & \ddots & -1 & 2 \left(\frac{1}{d} + 1 \right) & -1 \\
      0  & 0 & \ldots & 0 & -1 & 2 \left(\frac{1}{d} + 1 \right)
    \end{array}
  \right)
  \left(
    \begin{array}{c}
      T_1^{n+1}  \\
      T_2^{n+1}  \\
      T_3^{n+1}    \\
      \vdots \\
      T_{M-1}^{n+1} \\
      T_M^{n+1} 
\end{array}
\right)


= \left( \begin{array}{c}
    T_2^{n} + 2 \left(\frac{1}{d} - 1 \right) T_1^{n} + \left(T_0^n + T_0^{n+1}\right) \\
      T_3^{n} + 2 \left(\frac{1}{d} - 1 \right) T_2^{n} + T_1^n  \\
      T_4^{n} + 2 \left(\frac{1}{d} - 1 \right) T_3^{n} + T_2^n  \\
      \vdots \\
      T_M^{n} + 2 \left(\frac{1}{d} - 1 \right) T_{M-1}^{n} + T_{M-2}^n  \\
      \left(T_{M+1}^n + T_{M+1}^{n+1}\right) + 2 \left(\frac{1}{d} - 1 \right) T_{M}^{n} + T_{M-1}^n
    \end{array} \right)

Le résultat lui-même est à peu près le même que la différence centrale. Le nombre de diffusion d est de 0,1.

jacobi.png

La figure ci-dessous montre le résultat du calcul avec le nombre de diffusion d fixé à 1 comme dans le cas de la différence de centre. ** Vous pouvez voir que la solution peut être trouvée fermement même si le nombre de diffusion d est égal ou supérieur à 0,5 **. C'est l'avantage d'utiliser la méthode implicite. De plus, d'un point de vue physique, à mesure que le nombre de diffusion d augmente, la largeur du pas de temps et la largeur du pas d'espace sont constantes, de sorte que la vitesse de transfert de chaleur devient plus rapide et converge vers 150 K plus rapidement que lorsque le nombre de diffusion d = 0,1. Vous pouvez voir que ça marche.

jacobi.png

Un exemple de mise en œuvre est présenté ci-dessous. La fonction jacobi dans la syntaxe a été construite dans la section précédente.


temperature_jacobi = temperature_array.copy()
for n in range(Time_step):
    a_matrix = np.identity(len(temperature_jacobi)) * 2 *(1/d+1) \
                - np.eye(len(temperature_jacobi), k=1) \
                - np.eye(len(temperature_jacobi), k=-1)
    temp_temperature_array = np.append(np.append(
                        temperature_lower_boundary, 
                        temperature_jacobi), temperature_upper_boundary)
    b_array = 2 * (1/d - 1) * temperature_jacobi + temp_temperature_array[2:] + temp_temperature_array[:-2]
    b_array[0] += temperature_lower_boundary
    b_array[-1] += temperature_upper_boundary
    temperature_jacobi = jacobi(a_matrix, b_array, 1e-8)
plt.plot(x_array, exact_temperature_array, label="Answer", c="black")
plt.plot(x_array, temperature_jacobi, label="Implicit Jacobi(100 steps)")
plt.legend(loc="lower right")
plt.xlabel("x")
plt.ylabel("temperature")
plt.xlim(0, max(x_array))
plt.ylim(100, 200)

2-3-2. Méthode Gauss-Seidel (méthode Gauss-Seidel)

Lorsque $ A = D + L + U $ en utilisant une matrice D ayant des composantes diagonales de la matrice A, une matrice triangulaire supérieure U et une matrice triangulaire inférieure L,

　Ax = b\\
 (D+L+U) x = b\\
 (D+L)x^{(m+1)} = b - Ux^{(m)} \\
 x^{(m+1)} = D^{-1}(b- L x^{(m+1)} - U x^{(m)})\\
x^{(m+1)} = (D+L)^{-1}(b - U x^{(m)})

La méthode Gauss-Seidel est une méthode qui peut transformer $ Ax = b $ et trouve une solution approximative de manière itérative avec la formule en bas.

x_i^{(m+1)} = \frac{1}{a_{ii}}\left( b_i - \sum_{j=1}^{i-1} a_{ij} x_j^{(m+1)} - \sum_{j=i+1}^{n} a_{ij} x_j^{(m)} \right)

2-3-2-1. Précautions relatives à la méthode Gauss-Seidel

Comme mentionné dans la méthode de Jacobi, la méthode de Gauss-Seidel ne converge que si la matrice de coefficients A satisfait l'avantage diagonal (au sens étroit).

2-3-2-2. Implémentation de la méthode Gauss-Seidel

Il n'y a pas de changement particulier dans le graphique.

gauss-seidel.png

def gauss_seidel(a_matrix, b_array, target_residual):
    """
    Gauss-Méthode Seidel
    Ax = b
    A = (D+L+U)Puis
    x^{(m+1)} = D^{-1}(b- L x^{(m+1)} - U x^{(m)})
    x^{(m+1)} = (D+L)^{-1}(b - U x^{(m)})
Cependant, D est une matrice diagonale, L est une matrice triangulaire inférieure et U est une matrice triangulaire supérieure.
    
    Parameters
    ----------
    a_matrix : numpy.float64
matrice n × m
    b_array : numpy.float64
matrice à n lignes
    target_residual : numpy.float64
Nombre positif. Cibler les résidus.

    Returns
    -------
    x : numpy.float64
matrice à n lignes
    """
    x_old = b_array.copy()
    lower_matrix = np.tril(a_matrix) #Matrice triangulaire inférieure
    upper_matrix = a_matrix - lower_matrix  #Matrice triangulaire supérieure
    inv_lower_matrix = np.linalg.inv(lower_matrix)
    count = 0
    while True:
        count += 1
        x = np.dot(inv_lower_matrix, (b_array - np.dot(upper_matrix, x_old)))
        residual = np.linalg.norm(x - x_old) / np.linalg.norm(x)
        x_old = x
        if residual <= target_residual:
            break
        elif count >= 10000:
            import sys
            print(residual)
            sys.exit()
    return x

temperature_gs = temperature_array.copy()
for n in range(Time_step):
    a_matrix = np.identity(len(temperature_gs)) * 2 *(1/d+1) \
                - np.eye(len(temperature_gs), k=1) \
                - np.eye(len(temperature_gs), k=-1)
    temp_temperature_array = np.append(np.append(
                        temperature_lower_boundary, 
                        temperature_gs), temperature_upper_boundary)
    b_array = 2 * (1/d - 1) * temperature_gs + temp_temperature_array[2:] + temp_temperature_array[:-2]
    b_array[0] += temperature_lower_boundary
    b_array[-1] += temperature_upper_boundary
    temperature_gs = gauss_seidel(a_matrix, b_array, 1e-8)
plt.plot(x_array, exact_temperature_array, label="Answer", c="black")
plt.plot(x_array, temperature_gs, label="Implicit Gauss-Seidel(100 steps)")
plt.legend(loc="lower right")
plt.xlabel("x")
plt.ylabel("temperature")
plt.xlim(0, max(x_array))
plt.ylim(100, 200)

2-3-3. Méthode SOR

Méthode de Gauss Seidel avec facteur de relaxation $ \ omega $ ajouté. Identique à la méthode Gauss-Seidel lorsque $ \ omega = 1 $. Lorsque $ A = D + L + U $ en utilisant la matrice D ayant des composantes diagonales de la matrice A, de la matrice supérieure U et de la matrice inférieure L,

　Ax = b\\
 (D+L+U) x = b\\
 \tilde{x}^{(m+1)} = D^{-1}(b- L x^{(m+1)} - U x^{(m)}) \quad : Gauss-Seidel\\
 x^{(m+1)} = x^{(m)} + \omega \left(\tilde{x}^{(m+1)} - x^{(m)}  \right)

$ Ax = b $ peut être transformé comme ceci, et la méthode de recherche répétée de la solution approximative avec la formule en bas est appelée la méthode SOR (méthode de relaxation d'accélération séquentielle).

x_i^{(m+1)} = x_i^{(m)} + \omega \frac{1}{a_{ii}}\left( b_i - \sum_{j=1}^{i-1} a_{ij} x_j^{(m+1)} - \sum_{j=i}^{n} a_{ij} x_j^{(m)} \right)

2-3-3-1. Mise en œuvre de la méthode SOR

Cette fois, le coefficient de relaxation $ \ omega $ est une constante. Le coefficient de relaxation $ \ omega $ est également disponible [lorsque la valeur optimale peut être déterminée](méthode https://ja.wikipedia.org/wiki/SOR). Il n'y a pas de changement dans le graphique.

sor.png

def sor(a_matrix, b_array, target_residual):
    """
Méthode SOR
    Ax = b
    A = (D+L+U)Puis
    x~^{(m+1)} = D^{-1}(b- L x^{(m+1)} - U x^{(m)}) : Gauss-Seidel
    x^{(m+1)} = x^{(m)} + omega (x~^{(m+1)} - x^{(m)})
Cependant, D est une matrice diagonale, L est une matrice triangulaire inférieure et U est une matrice triangulaire supérieure.
    
    Parameters
    ----------
    a_matrix : numpy.float64
matrice n × m
    b_array : numpy.float64
matrice à n lignes
    target_residual : numpy.float64
Nombre positif. Cibler les résidus.

    Returns
    -------
    x : numpy.float64
matrice à n lignes
    """
    x_old = b_array.copy()
    lower_matrix = np.tril(a_matrix) #Matrice triangulaire inférieure
    upper_matrix = a_matrix - lower_matrix  #Matrice triangulaire supérieure
    inv_lower_matrix = np.linalg.inv(lower_matrix)
    omega = 1.9  #Dans cet exemple, la convergence peut être lente.
    count = 0
    while True:
        count += 1
        x_tilde = np.dot(inv_lower_matrix, (b_array - np.dot(upper_matrix, x_old)))
        x = x_old + omega * (x_tilde - x_old)
        residual = np.linalg.norm(x - x_old) / np.linalg.norm(x)
        x_old = x
        if residual <= target_residual:
            break
        elif count >= 10000:
            import sys
            print(residual)
            sys.exit()
    return x

temperature_sor = temperature_array.copy()
for n in range(Time_step):
    a_matrix = np.identity(len(temperature_sor)) * 2 *(1/d+1) \
                - np.eye(len(temperature_sor), k=1) \
                - np.eye(len(temperature_sor), k=-1)
    temp_temperature_array = np.append(np.append(
                        temperature_lower_boundary, 
                        temperature_sor), temperature_upper_boundary)
    b_array = 2 * (1/d - 1) * temperature_sor + temp_temperature_array[2:] + temp_temperature_array[:-2]
    b_array[0] += temperature_lower_boundary
    b_array[-1] += temperature_upper_boundary
    temperature_sor = sor(a_matrix, b_array, 1e-8)
plt.plot(x_array, exact_temperature_array, label="Answer", c="black")
plt.plot(x_array, temperature_sor, label="Implicit SOR(100 steps)")
plt.legend(loc="lower right")
plt.xlabel("x")
plt.ylabel("temperature")
plt.xlim(0, max(x_array))
plt.ylim(100, 200)

Résumé

~~ La prochaine fois, je voudrais résumer la méthode d'itération non stationnaire (méthode du sous-espace de Krylov). ~~ J'ai décidé de résumer un jour la méthode d'itération non stationnaire. L'utilisation de la méthode d'itération non stationnaire en Python est résumée dans l'article [Python] qui permet le calcul de matrice éparse à grande vitesse. La prochaine fois, l'équation de diffusion de transfert (linéaire), qui est une combinaison de l'équation de transfert de article précédent et l'équation de diffusion traitée dans cet article, et sa non-linéarisation Je voudrais parler de l'équation de Burgers.

Les références

  1. https://qiita.com/sci_Haru/items/960687f13962d63b64a0
  2. https://qiita.com/IshitaTakeshi/items/cf106c139660ef138185
  3. http://nkl.cc.u-tokyo.ac.jp/14n/PDE.pdf
  4. https://home.hiroshima-u.ac.jp/nakakuki/other/simusemi-20070907.pdf
  5. https://www.sit.ac.jp/user/konishi/JPN/L_Support/SupportPDF/InplicitMethod.pdf ← Incroyablement facile à comprendre
  6. http://ri2t.kyushu-u.ac.jp/~watanabe/RESERCH/MANUSCRIPT/KOHO/GEPP/GEPP.pdf ← Complet et détaillé

Référence itérative

La méthode itérative est le matériel et les vidéos répertoriés par Professeur Nakajima de l'Université de Tokyo, et ce [site] qui résume l'analyse numérique. Si vous regardez (http://www.slis.tsukuba.ac.jp/~fujisawa.makoto.fu/cgi-bin/wiki/index.php?Numerical%20Calculation), vous pouvez le voir à peu près. M. Nakajima est un dieu.

Références de la méthode Jacobi

Référence d'implémentation de la méthode de relaxation constante

Recommended Posts

[Python] Simulation de fluides: implémenter l'équation de diffusion
[Python] Simulation de fluides: implémenter l'équation de transfert
[Python] Simulation de fluide: équation de Navier-Stokes incompressible
Implémenter le modèle Singleton en Python
[Python] Simulation de fluide: de linéaire à non linéaire
[Python] Observons le vortex de Kalman par la méthode de Boltzmann en treillis. [Simulation de fluide]
Trouvez la solution de l'équation d'ordre n avec python
Résolution d'équations de mouvement en Python (odeint)
Mettre en œuvre REPL
Implémenter la solution de l'algèbre de Riccati en Python
Trouvez le maximum de Python
Mettre en œuvre des recommandations en Python
le zen de Python
Implémenter XENO avec python
Implémenter sum en Python
[Python] Fractionner la date
Implémenter Traceroute dans Python 3
[Python] Régression LASSO avec contrainte d'équation utilisant la méthode du multiplicateur
Avoir le graphique d'équation de la fonction linéaire dessiné en Python
J'ai essayé d'implémenter la fonction d'envoi de courrier en Python