[Python] Fluidsimulation: Implementieren Sie die Diffusionsgleichung

Einführung

Ich möchte (in mehreren Artikeln) das Wissen zusammenfassen, das zum Erstellen eines Codes für die numerische Flüssigkeitsanalyse für Wasser erforderlich ist, und gleichzeitig die numerische Fluiddynamik (CFD) untersuchen, eine Studie zur Simulation von Flüssigkeiten wie Luft und Wasser.

Ich bin der Meinung, dass die numerische Fluiddynamik eine ziemlich schwierige Studie ist, daher möchte ich sie so schreiben, dass sie für Anfänger leicht verständlich ist. Es scheint, dass es viele Fehler gibt. Bitte kontaktieren Sie uns, wenn Sie diese finden. Ich würde mich auch freuen, wenn Sie kommentieren könnten, was schwer zu verstehen ist. Wir werden es von Zeit zu Zeit aktualisieren.

Zielgruppe

or

Serie

Grober Inhalt dieses Artikels

Als ersten Schritt zur Behandlung der für die Wassersimulation erforderlichen Grundgleichungen von Flüssigkeit werden wir auch die Diffusionsgleichung kurz zusammenfassen und implementieren. ** Ich habe es verständlich gemacht, ohne den vorherigen Artikel zu lesen (wahrscheinlich) **

  1. [Mittendifferenz](# 1-Mittendifferenz)
  2. [Crank Nicholsons implizite Methode](# 2-Crank Nicholsons implizite Methode)
  3. [Direkte Methode](# 2-1-Direkte Methode)
  4. [Wiederholungsmethode](# 2-2-Wiederholungsmethode)
  5. [Konstante iterative Methode](# 2-2-1-Konstante iterative Methode)
  6. [Nichtstationäre Wiederholungsmethode](# 2-2-2-Nichtstationäre Wiederholungsmethode Krylov-Subraummethode)
  7. [Implementierung der stationären Iterationsmethode](# 2-3-Implementierung der stationären Iterationsmethode)
  8. [Jacobi-Methode](# 2-3-1-Jacobi-Methode)
  9. [Einfaches Beispiel](# 2-3-1-1-Einfaches Beispiel)
  10. [Vorsicht](# 2-3-1-2-Jacobi-Methode Vorsichtsmaßnahme)
  11. [Implementierung](Implementierung der Methode # 2-3-1-3-jacobi)
  12. [Gauß-Seidel-Methode](# 2-3-2-Gauß-Seidel-Methode)
  13. [Achtung](# 2-3-2-1-Gauß-Seidel-Methode Vorsichtsmaßnahme)
  14. [Implementierung](Implementierung der # 2-3-2-2-Gauß-Seidel-Methode)
  15. [SOR-Methode](# 2-3-3-Sor-Methode)
  16. [Implementierung](Implementierung der Methode # 2-3-3-1-sor)

Diffusionsgleichung (Formel zur Vereinheitlichung des Feldes)

Was ist eine eindimensionale Diffusionsgleichung? $ \frac{\partial T}{\partial t} = \kappa \frac{\partial^2 T}{\partial x^2} $

Es wird durch die Formel ausgedrückt. In Bezug auf die physikalische Bedeutung bedeutet dies, dass physikalische Größen zufällig einheitlich sind. Es bedeutet Wärmeleitung und Diffusion von Substanzen sowie Wärme. Wenn beispielsweise Milch zu Kaffee gegeben wird, breitet sie sich langsam ohne Rühren aus, oder wenn das Ende eines Metallstabs in einen Kochtopf gelegt wird, wie in der folgenden Abbildung gezeigt, wird der gesamte Stab allmählich zu heißem Wasser. Das Phänomen der gleichen Temperatur entspricht der Diffusion.

diffusion_phenomenon.png

Verfahren zur Unterscheidung dieser eindimensionalen Diffusionsgleichung umfassen die explizite zentrale Differenzmethode und die implizite Methode Crank Nicholsons implizite Methode. Zur Erinnerung: Die explizite Methode sagt den zukünftigen Wert eine Stunde später basierend auf dem bekannten Wert der aktuellen Zeit voraus, und die implizite Methode sagt den zukünftigen Wert auch eine Stunde später voraus. Wir werden diese beiden implementieren. Wenn Sie jeweils eine Differenzformel aufschreiben, lautet diese wie folgt.

  1. Mitteldifferenz Sie kann abgeleitet werden, indem die im vorherigen Artikel erläuterte Genauigkeitsdifferenz erster Ordnung genau durchgeführt wird. Dies ist eine Art expliziter Methode, die den Wert des nächsten Males n + 1 vorhersagt, wobei nur die Daten zum aktuellen Zeitpunkt n verwendet werden. $ \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. Crank Nicholsons implizite Methode

Ein Verfahren zum Mitteln eines bekannten Wertes (Zeitpunkt n) und eines unbekannten Werts einen Schritt später (Zeitpunkt n + 1) für die partielle Differenzierung in Bezug auf die räumliche Differenzierung. Wie der Name schon sagt, handelt es sich um eine implizite Methode, die auch den Wert des nächsten Males n + 1 verwendet, um den Wert des nächsten Males n + 1 vorherzusagen. Details werden später beschrieben. $ \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) $

Bestimmung der Stabilität anhand der Anzahl der Spreads

Von den oben gezeigten Ausbreitungsmethoden ist die ** explizite Methode unter Verwendung des zentralen Unterschieds ** $ d = \kappa \frac{\Delta t}{\Delta x^2} $ Es ist notwendig, die Stabilität der numerischen Berechnung unter Verwendung der sogenannten ** Diffusionszahl d ** zu bestimmen. Andererseits ist die implizite Methode unbedingt stabil.

Speziell,

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

Muss im Bereich von liegen. Das Wichtigste bei dieser ** Diffusionszahlbedingung ist, dass Sie, wenn Sie die Schrittgröße $ \ Delta x $ reduzieren möchten, die Zeitschrittgröße $ \ Delta t $ um ihr Quadrat ** reduzieren müssen. Wenn Sie beispielsweise $ \ Delta x $ mit 1/2 multiplizieren, müssen Sie $ \ Delta t $ auf 1/4 reduzieren. Infolgedessen erhöht das Lösen der Diffusionsgleichung mit der expliziten Methode tendenziell die Berechnungslast, sodass im Allgemeinen häufig die implizite Methode verwendet wird.

Darüber hinaus ergibt sich diese Bedingung aus der Stabilitätsanalyse von Von Neumann, aber ich werde sie hier nicht erklären, sondern nur ein intuitives Bild.

Wenn Sie die Formel des zentralen Unterschieds neu schreiben,

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

Es wird sein. Diffusion bedeutet, dass es ungeordnet uniformiert ist. Wenn also der Koeffizient $ (1-2d) $ von $ T_j ^ n $ negativ wird, passiert er mehr als die physikalische Größe des j-ten Gitters daneben. Daher sollte es $ d \ leq0.5 $ sein.

Ich denke, es ist schwer zu verstehen, deshalb erkläre ich es anhand des folgenden Beispiels für Geld (** Es ist nur ein Beispiel, um ein Bild zu haben ). Angenommen, Sie haben eine Gruppe von N Personen, und irgendwann hat die j-1. Person 30 Yen, die j-te Person 100 Yen und die j + 1-Person 50 Yen. Das Verbreiten ist eine Operation, bei der diese N Personen den gleichen Geldbetrag haben ( Vereinheitlichung **). Tauschen Sie also einen bestimmten Prozentsatz des Geldes mit den Nachbarn aus, bis der gleiche Geldbetrag erreicht ist. Ich werde. Wenn Sie sich die Formel für die Mittendifferenz genau ansehen, ist dieser konstante Prozentsatz die Diffusionszahl d. Wenn Sie also d = 0,1 setzen und berechnen, hat die j-te Person eine Stunde später 88 nen bei n + 1 Mal. Wenn Sie so denken, können Sie leicht erkennen, dass die Diffusionszahl d nicht mehr als 1/2 haben darf.

money_diffusion_r2.png ** * Dies ist ein Beispiel, um die Diffusionsatmosphäre zu erfassen **

Berechnungsbedingung

problem.png

answer.png

Die Gitterbreite $ \ Delta x $, der Wärmeleitungskoeffizient $ \ kappa $ und die Zeitschrittbreite $ \ Delta t $ werden alle als Konstanten berechnet. Diesmal ist $ \ Delta x = 1 $, $ \ Delta t = 0,2 $, $ \ kappa = 0,5 $, und die Berechnung wird unter stabilen Bedingungen mit der Diffusionszahl $ d = 0,1 $ durchgeführt.

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)

Implementierung

1. Mitteldifferenz

\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}

Die folgende Abbildung zeigt das Ergebnis, wenn die Diffusionszahl d 0,1 beträgt. Ich habe versucht, das Diagramm auszugeben, indem ich es in mehrere Zeitschritte (Zeit) unterteilt habe. Nach 100 Stunden haben beide Enden fast 150 K erreicht, und die Temperaturverteilung wird wie eine Sinuskurve daran gezogen. Nach 5000 Stunden nähert sich die Temperatur des gesamten Stabes allmählich 150 K, was darauf hinweist, dass die Implementierung der Mittendifferenz erfolgreich ist. Da es einfach ist, dem Ausdruck zu folgen, werde ich ihn schreiben, ohne die numpy-Funktion so oft wie möglich zu verwenden.

explicit.png

Übrigens, wenn Sie $ \ kappa = 5 $ setzen und mit der Diffusionszahl $ d = 1 $ berechnen, können Sie sehen, dass die Berechnung divergiert und die Lösung nicht erhalten werden kann, wie in der folgenden Abbildung gezeigt.

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. Crank Nicholsons implizite Methode

\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)

Wenn Sie dies transformieren und den Wert der Zeit n + 1 auf die linke Seite bringen $ -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} $

Unter der Annahme, dass die Gitterpunkte im Bereich von 1 bis M liegen, wird der Grenzwert durch $ T_0, T_ {M + 1} $ dargestellt und in eine Matrixanzeige konvertiert.


\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)

Ziehen Sie in Betracht, diese simultane eindimensionale Gleichung zu lösen. Wenn die Anzahl der Gitterpunkte M klein ist, ist die unbekannte Anzahl klein, sodass eine direkte Methode verwendet wird, die direkt eine Lösung findet, z. B. eine Methode namens Gaußsche direkte Methode. Je größer M ist, desto schwieriger ist es jedoch, eine Berechnungslösung zu finden. Daher wird im Allgemeinen zur Lösung einer simultanen eindimensionalen Gleichung im großen Maßstab eine Methode verwendet, die als ** iterative Methode ** bezeichnet wird und die Näherungslösung zu einer exakten Lösung konvergiert.

2-1. Direkte Methode

Die Typen sind wie folgt. Ich werde erklären, wann es gebraucht wird.

2-2. Iterative Methode

Ax=b

Ziehen Sie in Betracht, eine eindimensionale Simultangleichung zu lösen. Dies $ r=b-Ax' $ Das Verfahren zum Finden einer ungefähren Lösung der exakten Lösung $ x ^ * = A ^ {-1} b $ wird wiederholt, indem der geschätzte Wert $ x '$ iterativ geändert wird, bis r ausreichend klein wird. Es ist das Gesetz. Die iterative Methode kann in die stationäre und die instationäre Methode (Krylov-Subraummethode) unterteilt werden.

2-2-1. Methode der konstanten Wiederholung

Die stationäre Methode (stationäre iterative Methode) ist eine Methode, die sich während der iterativen Berechnung mit Ausnahme des Lösungsvektors x nicht ändert. Da es im Allgemeinen langsam ist, verwenden Sie entweder die später beschriebene instationäre Methode, wenn Sie numerische Berechnungen durchführen, oder verwenden Sie sie zusammen als Vorverarbeitung (Verarbeitung, um eine ungefähre ungefähre Lösung zu finden) der instationären Methode. Überlegen.

Die folgenden drei Methoden können als typische Methoden genannt werden. Im Folgenden wird der geschätzte Wert der Lösung in m Iterationen als $ x ^ {(m)} $ ausgedrückt, und der geschätzte Wert $ x ^ {(m)} $ in der m-ten Iteration wird als m + 1-te Iteration bezeichnet. Hier ist ein Beispiel, wie der geschätzte Wert von $ x ^ {(m + 1)} $ ermittelt wird.

  1. Jacobi-Methode

Eine Methode zur Berechnung der geschätzten Lösung $ x_i ^ {(m + 1)} $ in der i-ten Zeile der m + 1-Iteration unter Verwendung nur der bekannten geschätzten Lösung $ 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. Gauß-Seidel-Methode

Wenn die geschätzte Lösung $ x_i ^ {(m + 1)} $ in der i-ten Zeile der m + 1-Iteration gefunden wird, wird die bekannte geschätzte Lösung $ x_i ^ {(m)} $ bereits als m + 1 berechnet Geschätzte Lösung für Iterative $ x_ {0} ^ {(m + 1)}, \ cdots, x_ {i-1} ^ {(m + 1)} . Grundsätzlich endet die Konvergenzberechnung schneller als die Jacobi-Methode. $ 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. SOR-Methode (Sukzessive Überentspannung) / Relaxationsmethode mit sequentieller Beschleunigung

Gauß-Seidel-Methode plus Relaxationsfaktor $ \ omega $. Wenn der Relaxationskoeffizient $ \ omega = 1 $ ist, wird die Gauß-Seidel-Methode verwendet. Wenn es 1 oder weniger ist, ist die Konvergenz langsamer als bei Gauß-Seidel, aber Probleme, die mit der Gauß-Seidel-Methode nicht gelöst werden können, können gelöst werden. Stellen Sie den Wert auf 1 oder mehr ein. Wenn er jedoch zu groß eingestellt ist, wird er abweichen. Daher ist es wichtig, einen geeigneten Wert auszuwählen. Es scheint, dass 1.9 oft ausgewählt wird. $ 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. Multi-Grid-Methode Es gibt die geometrische Multi-Grid-Methode und die algebrische Multi-Grid-Methode (AMG). Die letztere AMG-Methode wird häufig als Vorbehandlung für instationäre Methoden verwendet. Ich denke, dass viele Leute es zum ersten Mal wissen werden, also werde ich die Details in einem anderen Artikel schreiben.

Details der früheren drei Punkte werden später beschrieben.

2-2-2. Nichtstationäre Wiederholungsmethode (Krylov-Subraummethode)

Ist es wie ein guter Junge mit direkten und iterativen Methoden? Es wird als iterative Methode klassifiziert.

Die bekannteste instationäre Methode ist die Conjugate Gradient-Methode (CG-Methode). Eine ausführliche Erklärung finden Sie in diesem Artikel. Es gibt eine Figur und sie ist wahnsinnig leicht zu verstehen. Darüber hinaus erklärt dieser Artikel den Unterschied zu der Methode des steilsten Abstiegs, die häufig beim maschinellen Lernen verwendet wird. Das Bild der CG-Methode ist ebenfalls geschrieben und leicht zu verstehen.

  1. Symmetrische Matrix
  2. Gradientenmethode konjugieren (CG-Methode)
  3. Vorkonditionierte konjugierte Gradientenmethode (PCG-Methode)
  4. Unvollständige Cholesky-Konjugat-Gradientenmethode (ICCG-Methode)
  5. Unvollständige LU-Konjugat-Gradienten-Methode (ILUCG)
  6. Andere
  7. Asymmetrische Matrix
  8. Bi-Konjugat-Gradient (BiCG-Methode)
  9. Konjugieren Sie das verbleibende Quadrat (CGS-Methode)
  10. BiCG-Stabilisierung (BiCG STAB-Methode)
  11. Generalized Conjugate Residual (GCR-Methode)
  12. Generalized Minimal Residual (GMRES-Methode)

Diese Techniken werden grundsätzlich als Bibliothek bereitgestellt. Ich denke auch, dass es auf Github viele Berechnungscodes gibt. Wenn Sie genauer studieren möchten, lesen Sie bitte auch dies. In einem anderen Artikel möchte ich so viele Methoden wie möglich für die instationäre Methode erklären und implementieren.

In diesem Artikel möchte ich vorerst die Jacobi-Methode, die Gauß-Seidel-Methode und die SOR-Methode verwenden, die die stationären Iterationsmethoden unter den bisher erläuterten Algorithmen sind.

2-3. Implementierung der stationären Iterationsmethode

Die Implementierung wird unter Bezugnahme auf das Material von Professor Nakajima von der Universität Tokio durchgeführt. Wenn Sie die Erklärung schwer verständlich finden, würden wir uns freuen, wenn Sie sich darauf beziehen könnten.

Im folgenden Beispiel

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)

Definieren Sie die Matrix wie folgt. D ist nur die diagonale Komponente von A, L ist die Unterseite der Diagonalkomponente und U ist die Oberseite der Diagonalkomponente.

2-3-1. Jacobi-Methode

Wenn $ A = D + L + U $ unter Verwendung der Matrix D mit diagonalen Komponenten der Matrix A, der oberen Matrix U und der unteren Matrix L,

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

$ Ax = b $ kann so transformiert werden, und die Methode, wiederholt die ungefähre Lösung mit der Formel unten zu finden, wird als Jacobi-Methode bezeichnet.

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. Einfaches Beispiel

Sei A eine $ 3 \ times3 $ Matrix


\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)

Ziehen Sie in Betracht, eine eindimensionale Simultangleichung zu lösen.

Wenn Sie die Anzahl der Iterationen als (m) ausdrücken, ist der Wert von $ x ^ {(m + 1)} $ bei (m + 1) nach einer Iteration

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}

Es kann in Form von ausgedrückt werden. Wie Sie aus der Form dieser Gleichung sehen können, gilt $ x ^ {(m + 1)} = D ^ {-1} (b- (L + R) x ^ {(m)}) $. Ich werde.

Dann kann durch Wiederholen der obigen Berechnung, bis die unten gezeigte Konvergenzbedingung erfüllt ist, die genaue Lösung $ x ^ * $ erhalten werden.

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

$ \ Epsilon $ sollte jedoch eine ausreichend kleine positive Zahl sein.

Lassen Sie uns nun mit einem einfachen Problem mit der Jacobi-Methode üben.

Als einfaches Beispiel

\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)

Denken Sie über das Lösen nach. Die Antwort ist

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

Im Gegensatz zu früher werden wir es implementieren, indem wir numpy-Funktionen stark nutzen. Weil es überwiegend leicht zu sehen ist. Die Implementierung ist wie folgt. Wenn Sie tatsächlich rechnen

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

Wird angezeigt und es kann bestätigt werden, dass die Antwort gesucht wird.

def jacobi(a_matrix, b_array, target_residual):
    """
Jacobi-Methode
    Ax = b
    A = (D+L+U)Dann
    Dx = b - (L+U)x
    x = D^{-1} (b - (L+U)x)
D ist jedoch eine Diagonalmatrix L.+Sei U die Restmatrix.
    
    Parameters
    ----------
    a_matrix : numpy.float64
m × n Matrix
    b_array : numpy.float64
m-Zeilenmatrix
    target_residual : numpy.float64
Positive Zahl. Zielreste.

    Returns
    -------
    x : numpy.float64
m-Zeilenmatrix
    """
    x = b_array.copy()
    x_old = b_array.copy()
    
    diag_matrix= np.diag(a_matrix)  #Diagonale Matrix
    l_u_matrix = a_matrix - np.diagflat(diag_matrix)  #Restmatrix
    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. Vorsichtsmaßnahmen für die Jacobi-Methode

Die Jacobi-Methode (und die Gauß-Seidel-Methode) konvergieren nur, wenn die Matrix A den diagonalen Vorteil (im engeren Sinne) erfüllt. Um kurz den diagonalen Vorteil (im engeren Sinne) anzugeben, heißt es, dass die diagonale Komponente größer ist als die Summe der Absolutwerte der anderen Komponenten in derselben Zeile. Versuchen Sie in der obigen Übung, mit den nicht diagonalen Komponenten zu spielen, die größer sind.

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

2-3-1-3. Implementierung der Jacobi-Methode

Ich werde die zu lösende Formel erneut veröffentlichen.


\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)

Das Ergebnis selbst entspricht in etwa dem zentralen Unterschied. Die Diffusionszahl d beträgt 0,1.

jacobi.png

Die folgende Abbildung zeigt das Ergebnis der Berechnung, wobei die Diffusionszahl d wie im Fall der Mittendifferenz auf 1 gesetzt ist. ** Sie können sehen, dass die Lösung auch dann fest gefunden werden kann, wenn die Diffusionszahl d 0,5 oder mehr beträgt **. Dies ist der Vorteil der impliziten Methode. Aus physikalischer Sicht sind die Zeitschrittbreite und die Raumschrittbreite mit zunehmender Diffusionszahl d konstant, so dass die Wärmeübertragungsgeschwindigkeit schneller wird und schneller zu 150 K konvergiert als bei einer Diffusionszahl d = 0,1. Sie können sehen, dass es tut.

jacobi.png

Ein Implementierungsbeispiel ist unten gezeigt. Die Jacobi-Funktion in der Syntax wurde im vorherigen Abschnitt erstellt.


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. Gauß-Seidel-Methode (Gauß-Seidel-Methode)

Wenn $ A = D + L + U $ unter Verwendung einer Matrix D mit diagonalen Komponenten der Matrix A, einer oberen Dreiecksmatrix U und einer unteren Dreiecksmatrix 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)})

$ Ax = b $ kann so transformiert werden, und die Methode, wiederholt die ungefähre Lösung mit der Formel unten zu finden, wird als Gauß-Seidel-Methode bezeichnet.

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. Vorsichtsmaßnahmen für die Gauß-Seidel-Methode

Wie in der Jacobi-Methode erwähnt, konvergiert die Gauß-Seidel-Methode nur, wenn die Koeffizientenmatrix A den diagonalen Vorteil (im engeren Sinne) erfüllt.

2-3-2-2. Implementierung der Gauß-Seidel-Methode

Es gibt keine besondere Änderung in der Grafik.

gauss-seidel.png

def gauss_seidel(a_matrix, b_array, target_residual):
    """
    Gauss-Seidel-Methode
    Ax = b
    A = (D+L+U)Dann
    x^{(m+1)} = D^{-1}(b- L x^{(m+1)} - U x^{(m)})
    x^{(m+1)} = (D+L)^{-1}(b - U x^{(m)})
D ist jedoch eine Diagonalmatrix, L ist eine untere Dreiecksmatrix und U ist eine obere Dreiecksmatrix.
    
    Parameters
    ----------
    a_matrix : numpy.float64
n × m Matrix
    b_array : numpy.float64
n-reihige Matrix
    target_residual : numpy.float64
Positive Zahl. Zielreste.

    Returns
    -------
    x : numpy.float64
n-reihige Matrix
    """
    x_old = b_array.copy()
    lower_matrix = np.tril(a_matrix) #Untere dreieckige Matrix
    upper_matrix = a_matrix - lower_matrix  #Obere Dreiecksmatrix
    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. SOR-Methode

Gauß-Seidel-Methode mit hinzugefügtem Relaxationsfaktor $ \ omega $. Entspricht der Gauß-Seidel-Methode, wenn $ \ omega = 1 $. Wenn $ A = D + L + U $ unter Verwendung der Matrix D mit diagonalen Komponenten der Matrix A, der oberen Matrix U und der unteren Matrix 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 $ kann so transformiert werden, und die Methode zum wiederholten Finden der Näherungslösung mit der Formel unten wird als SOR-Methode (sequentielle Beschleunigungsrelaxationsmethode) bezeichnet.

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. Implementierung der SOR-Methode

Diesmal ist der Relaxationskoeffizient $ \ omega $ eine Konstante. Der Relaxationskoeffizient $ \ omega $ ist ebenfalls verfügbar wenn der optimale Wert bestimmt werden kann. Es gibt keine Änderung in der Grafik.

sor.png

def sor(a_matrix, b_array, target_residual):
    """
SOR-Methode
    Ax = b
    A = (D+L+U)Dann
    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)})
D ist jedoch eine Diagonalmatrix, L ist eine untere Dreiecksmatrix und U ist eine obere Dreiecksmatrix.
    
    Parameters
    ----------
    a_matrix : numpy.float64
n × m Matrix
    b_array : numpy.float64
n-reihige Matrix
    target_residual : numpy.float64
Positive Zahl. Zielreste.

    Returns
    -------
    x : numpy.float64
n-reihige Matrix
    """
    x_old = b_array.copy()
    lower_matrix = np.tril(a_matrix) #Untere dreieckige Matrix
    upper_matrix = a_matrix - lower_matrix  #Obere Dreiecksmatrix
    inv_lower_matrix = np.linalg.inv(lower_matrix)
    omega = 1.9  #In diesem Beispiel kann die Konvergenz langsam sein.
    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)

Zusammenfassung

~~ Das nächste Mal möchte ich die instationäre Iterationsmethode (Krylov-Subraummethode) zusammenfassen. ~~ Ich beschloss, eines Tages die instationäre Iterationsmethode zusammenzufassen. Die Verwendung der instationären Iterationsmethode in Python ist in [Python] -Artikel zusammengefasst, der die Hochgeschwindigkeitsberechnung mit geringer Dichte ermöglicht. Das nächste Mal die (lineare) Übertragungsdiffusionsgleichung, die eine Kombination aus der Übertragungsgleichung von vorheriger Artikel und der in diesem Artikel behandelten Diffusionsgleichung und ihrer Nichtlinearisierung ist Ich möchte mich mit der Burgers-Gleichung befassen.

Verweise

  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 ← Wahnsinnig einfach zu verstehen
  6. http://ri2t.kyushu-u.ac.jp/~watanabe/RESERCH/MANUSCRIPT/KOHO/GEPP/GEPP.pdf ← Umfassend und detailliert

Iterative Referenz

Die iterative Methode sind die Materialien und Videos, die von Professor Nakajima der Universität Tokio aufgelistet wurden, und diese [Site], die die numerische Analyse zusammenfasst. Wenn Sie sich (http://www.slis.tsukuba.ac.jp/~fujisawa.makoto.fu/cgi-bin/wiki/index.php?Numerical%20Calculation) ansehen, können Sie es ungefähr sehen. Herr Nakajima ist ein Gott.

Jacobi-Methodenreferenzen

Implementierungsreferenz für die Methode der konstanten Relaxation

Recommended Posts

[Python] Fluidsimulation: Implementieren Sie die Diffusionsgleichung
[Python] Flüssigkeitssimulation: Implementieren Sie die Übertragungsgleichung
[Python] Flüssigkeitssimulation: Inkompressible Navier-Stokes-Gleichung
Implementieren Sie das Singleton-Muster in Python
[Python] Fluidsimulation: Von linear zu nichtlinear
[Python] Beobachten wir den Kalman-Wirbel nach der Gitter-Boltzmann-Methode. [Flüssigkeitssimulation]
Finden Sie die Lösung der Gleichung n-ter Ordnung mit Python
Lösen von Bewegungsgleichungen in Python (odeint)
Implementieren Sie REPL
Implementieren Sie die Lösung der Riccati-Algebra in Python
Finden Sie das maximale Python
der Zen von Python
Implementieren Sie XENO mit Python
Implementieren Sie sum in Python
[Python] Teilen Sie das Datum
Implementieren Sie Traceroute in Python 3
[Python] LASSO-Regression mit Gleichungsbeschränkung unter Verwendung der Multiplikatormethode
Lassen Sie das Gleichungsdiagramm der linearen Funktion in Python zeichnen
Ich habe versucht, die Mail-Sendefunktion in Python zu implementieren