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.
or
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) **
Was ist eine eindimensionale Diffusionsgleichung?
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.
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.
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.
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.
Von den oben gezeigten Ausbreitungsmethoden ist die ** explizite Methode unter Verwendung des zentralen Unterschieds **
Speziell,
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,
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.
** * Dies ist ein Beispiel, um die Diffusionsatmosphäre zu erfassen **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)
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.
Ü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.
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)
Wenn Sie dies transformieren und den Wert der Zeit n + 1 auf die linke Seite bringen
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.
Die Typen sind wie folgt. Ich werde erklären, wann es gebraucht wird.
Ziehen Sie in Betracht, eine eindimensionale Simultangleichung zu lösen.
Dies
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.
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)}
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)}
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.
Details der früheren drei Punkte werden später beschrieben.
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.
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.
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.
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.
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.
$ \ 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))
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.
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.
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.
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)
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.
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.
Es gibt keine besondere Änderung in der Grafik.
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)
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.
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.
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)
~~ 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.
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.
Recommended Posts