Iterative Verfahren

Fixpunktiteration

Ein Beispiel aus der Rohrhydraulik:

Zur Bestimmung der Rohrreibungszahl $\lambda$ kann bei glatten, turbulent durchströmten Rohren die implizite Formel von Prandtl verwendet werden:

$$\frac{1}{\sqrt\lambda} = 2\cdot\log\left(\text{Re}\cdot\sqrt\lambda\right) - 0,8$$

Mathe-Nerds lösen die Gleichung mithilfe der Lambertschen W-Funktion, Ingenieure jedoch meist mit einem numerischen Verfahren.

Die Lösung für die Prandtl-Formel lässt sich grafisch als Schnittpunkt der linken und rechten Seite finden.

Wir definieren für die linke und rechte Seite separat:

$$LHS\left(\lambda\right) = \frac{1}{\sqrt\lambda}$$$$RHS\left(\lambda\right) = 2\cdot\log\left(\text{Re}\cdot\sqrt\lambda\right) - 0,8$$

Für eine Reynoldszahl $\text{Re} = 6,4\cdot 10^6$ ergeben sich dann folgende Kurvenverläufe:


In [1]:
import matplotlib.pyplot as plt
import numpy as np
import math
%config InlineBackend.figure_format = 'svg'
%matplotlib inline

# linke Seite der Gleichung (left hand side)
def LHS(lamb):
    return 1/np.sqrt(lamb)

# rechte Seite der Gleichung (right hand side)
def RHS(lamb, Re):
    return 2.0 * np.log10(Re * np.sqrt(lamb)) - 0.8

# Array mit äquidistanten Werten für lambda und Festlegen der Re-Zahl:
lamb = np.arange(0.0001, 0.05, 0.0001);
Re = 6.4e6

plt.plot(lamb, LHS(lamb), label="LHS")
plt.plot(lamb, RHS(lamb, Re), label="RHS")
plt.axis([0, 0.02, 5, 15])
plt.ylabel('LHS, RHS')
plt.xlabel('Rohrreibungszahl $\lambda$')
plt.grid()
plt.legend()
plt.show();


Als Lösung ("Fixpunkt"), bei dem linke und rechte Seite denselben Wert annehmen lässt sich für $\lambda$ ein Wert von etwa 0,0085 ablesen.

Mithilfe einer Fixpunktiteration lässt sich dieser auch iterativ bestimmen. Hierzu formen wir die Prandtl-Formel so um, dass auf der linken Seite nur noch $\lambda$ steht:

$$\lambda = \frac{1}{\left[RHS\left(\lambda\right)\right]^2}$$

Damit lässt sich nun eine Iterationsvorschrift formulieren:

$$\lambda_{i+1} = \frac{1}{\left[RHS\left(\lambda_i\right)\right]^2}$$

In [2]:
# Startwert für lambda:
lamb_alt = 100

# Liste, um Zwischenergebnisse zu speichern
lambda_i = []
lambda_i.append(lamb_alt)

# Fixpunkt-Algorithmus
for iteration in range(0, 5):
    lamb_neu = 1 / (RHS(lamb_alt, Re)**2)
    lambda_i.append(lamb_neu)
    lamb_alt = lamb_neu
    
fehler = (RHS(lamb_neu, Re)-LHS(lamb_neu)) / RHS(lamb_neu, Re)

plt.plot(lamb, LHS(lamb))
plt.plot(lamb, RHS(lamb, Re))
plt.plot(lambda_i, LHS(lambda_i), 'o')
for i, txt in enumerate(lambda_i):
    plt.annotate(i, (lambda_i[i], LHS(lambda_i[i])))
plt.axis([0, 0.02, 5, 15])
plt.ylabel('LHS, RHS')
plt.xlabel('Rohrreibungszahl $\lambda$')
plt.show();

print ("lambda = ", lamb_neu, ", Fehler in %: ", fehler*100)


lambda =  0.008653767448978819 , Fehler in %:  -0.00149506315410304

Newton-Verfahren

Das Newton-Verfahren ist eine weitere Möglichkeit, um Gleichungen iterativ zu lösen. Hierzu wird die Gleichung so umgestellt, dass sich das Problem in eine Nullstellensuche konvertiert. Die oben behandelte Rohrreibungsgleichung ergibt dann:

$$f(\lambda) = \frac{1}{\sqrt\lambda} - 2\cdot\log\left(\text{Re}\cdot\sqrt\lambda\right) + 0,8 = 0$$

Ausgehend von einem geschätzten Startwert für die Nullstelle $\lambda_i$ wird die Steigung $f'(\lambda_i)$ der Funktion berechnet. Die Tangente im Punkt $(\lambda_i,f(\lambda_i))$ ist dann:

$$t(\lambda) = f(\lambda_i) + f'(\lambda_i)\cdot (\lambda - \lambda_i)$$

Der Schnittpunkt dieser Tangente mit der $\lambda$-Achse ergibt den neuen Näherungswert für die Nullstelle und damit die Iterationsvorschrift:

$$\lambda_{i+1} = \lambda_i - \frac{f(\lambda_i)}{f'(\lambda_i)}$$

Im Beispiel mit der Rohrreibungsgleichung ist die Ableitung:

$$f'(\lambda) = -\frac{1}{2\cdot\lambda^{3/2}} - \frac{1}{\lambda}$$

In [8]:
# Startwert für lambda:
lamb_alt = 0.01

# Liste, um Zwischenergebnisse zu speichern
lambda_newton_i = []
lambda_newton_i.append(lamb_alt)

# die Funktion f
def f(lamb, Re):
    return 1/np.sqrt(lamb) - 2.0 * np.log10(Re * np.sqrt(lamb)) + 0.8

# die Ableitung von f
def f_strich(lamb):
    return -1/(2*lamb**1.5) - 1/(lamb*math.log(10))

# die Tangente an f (nur zur Visualisierung, wird eigentlich nicht benötigt)
def tangente_f(x, lamb, Re):
    return f(lamb,Re)+f_strich(lamb)*(x-lamb)

# Newton-Verfahren:
for iteration in range(0, 15):
    lamb_neu = lamb_alt - f(lamb_alt, Re)/f_strich(lamb_alt)
    lambda_newton_i.append(lamb_neu)
    lamb_alt = lamb_neu

fehler = (RHS(lamb_neu, Re)-LHS(lamb_neu)) / RHS(lamb_neu, Re)

# Ergebnisse im Diagramm darstellen:
plt.plot(lamb, f(lamb, Re))
plt.plot(lambda_newton_i, f(lambda_newton_i, Re), 'o')
for i, txt in enumerate(lambda_newton_i):
    plt.annotate(i, (lambda_newton_i[i], f(lambda_newton_i[i], Re)))
    plt.plot(lamb, tangente_f(lamb, lambda_newton_i[i], Re))
plt.plot([0,0.02],[0,0],'k', linewidth=1)
plt.axis([0, 0.02, -4, 6])
plt.ylabel('f = LHS - RHS')
plt.xlabel('Rohrreibungszahl $\lambda$')
plt.show();

print ("lambda = ", lamb_neu, ", Fehler in %: ", fehler*100)


lambda =  0.008654006863945804 , Fehler in %:  0.0

Verfahren zur Lösung von Gleichungssystemen

Weitere Verfahren, zur Lösung von ganzen Gleichungssystemen werden in Kapitel 4 vorgestellt. Diese Verfahren werden z.B. verwendet, um die riesigen Gleichungssysteme zu lösen, die bei der Diskretisierung von Transportgleichungen mithilfe der Finite-Differenzen- (FDM), Finite-Elemente- (FEM) oder Finite-Volumen-Methode (FVM) entstehen.

Prominente Vertreter sind das Gauß-Verfahren (Gauß-Seidel-Verfahren) und der Thomas-Algorithmus.

Hier geht's weiter oder hier zurück zur Übersicht.


Der folgende Python-Code darf ignoriert werden. Er dient nur dazu, die richtige Formatvorlage für die Jupyter-Notebooks zu laden.


In [1]:
from IPython.core.display import HTML
def css_styling():
    styles = open('TFDStyle.css', 'r').read()
    return HTML(styles)
css_styling()


Out[1]:
/* */

In [ ]: