Richardson-Verfahren

Das Richardson-Verfahren ist in der numerischen Mathematik ein iterativer Algorithmus zur näherungsweisen Lösung von linearen Gleichungssystemen. Es wurde 1910 vom britischen Meteorologen Lewis Fry Richardson entwickelt und zählt wie das Gauß-Seidel-Verfahren zur Klasse der Splitting-Verfahren.

Richardson-Verfahren

Zur Lösung des linearen Gleichungssystems A x = b {\displaystyle Ax=b} wird der Teil des Systems, der durch Abspaltung der Einheitsmatrix I {\displaystyle I} von der Systemmatrix A {\displaystyle A} entsteht (Splitting-Verfahren mit B = I {\displaystyle B=I} ) auf der linken Seite belassen, alles andere auf die rechte Seite verschoben, so dass eine Fixpunktgleichung

x = Φ ( x ) {\displaystyle x=\Phi (x)}    mit    Φ ( x ) = b ( A I ) x {\displaystyle \Phi (x)=b-(A-I)x}

entsteht, deren Lösung auf die Iteration

x k + 1 = ( I A ) x k + b {\displaystyle x_{k+1}=(I-A)x_{k}+b}

führt. Dieses Verfahren konvergiert wie jede Fixpunktiteration dieser Art, falls der Spektralradius der Iterationsmatrix G := I A {\displaystyle G:=I-A} echt kleiner eins ist.

Relaxiertes Richardsonverfahren

Die Iterationsformel des relaxierten Richardson-Verfahrens lautet

x k + 1 = ( I ω A ) x k + ω b {\displaystyle x_{k+1}=\left(I-\omega A\right)x_{k}+\omega b} .

Dabei wird in jedem Schritt das Residuum mit einem Faktor ω {\displaystyle \omega } gewichtet. Falls A {\displaystyle A} eine symmetrisch positiv definite Matrix ist, so gilt für den optimalen Relaxionsparameter

ω o p t = 2 λ m a x ( A ) + λ m i n ( A ) {\displaystyle \omega _{\mathrm {opt} }={\frac {2}{\lambda _{\mathrm {max} }(A)+\lambda _{\mathrm {min} }(A)}}} .

Dabei bezeichnen λ m a x ( A ) {\displaystyle \lambda _{\mathrm {max} }(A)} und λ m i n ( A ) {\displaystyle \lambda _{\mathrm {min} }(A)} den maximalen und minimalen Eigenwert von A {\displaystyle A} . Für den Spektralradius der Iterationsmatrix G ω o p t {\displaystyle G_{\omega _{\mathrm {opt} }}} gilt

ρ ( G ω o p t ) = κ 2 ( A ) 1 κ 2 ( A ) + 1 < 1 {\displaystyle \rho (G_{\omega _{\mathrm {opt} }})={\frac {\kappa _{2}(A)-1}{\kappa _{2}(A)+1}}<1} ,

wobei κ 2 ( A ) {\displaystyle \kappa _{2}(A)} die spektrale Kondition der Matrix A {\displaystyle A} bezeichnet. Das relaxierte Richardson-Verfahren konvergiert dann genauso „schnell“ wie das Gradientenverfahren bei symmetrischen Matrizen, wofür man jedoch keinen Relaxionsparameter berechnen muss. Dafür kann man mit dem Richardsonverfahren auch bei unsymmetrischen Matrizen mit komplexen Eigenwerten noch Konvergenz erzwingen, solange deren Realteile alle positiv sind.

Das Verfahren ist als Glätter in Mehrgitterverfahren geeignet.

Zyklisches Richardsonverfahren

Die Konvergenz lässt sich erheblich verbessern, wenn man mehrere Schritte der Iteration mit unterschiedlichen Parametern ω k {\displaystyle \omega _{k}} betrachtet. Man führt dazu jeweils m {\displaystyle m} Schritte

x k = x k 1 + ω k ( b A x k 1 ) , k = 1 , , m {\displaystyle x_{k}=x_{k-1}+\omega _{k}(b-Ax_{k-1}),\quad k=1,\ldots ,m}

zyklisch durch. Das Verfahren konvergiert, wenn der Spektralradius des Matrixpolynoms

p m ( A ) = ( I ω 1 A ) ( I ω 2 A ) ( I ω m A ) {\displaystyle p_{m}(A)=(I-\omega _{1}A)(I-\omega _{2}A)\cdots (I-\omega _{m}A)}

kleiner als eins ist, und umso besser je kleiner er ist. Für eine Matrix A {\displaystyle A} mit reellen und positiven Eigenwerten kann der Spektralradius durch das Maximum des reellen Polynoms p m ( t ) = ( 1 ω 1 t ) ( 1 ω 2 t ) ( 1 ω m t ) {\displaystyle p_{m}(t)=(1-\omega _{1}t)(1-\omega _{2}t)\cdots (1-\omega _{m}t)} im Intervall [ λ m i n , λ m a x ] {\displaystyle [\lambda _{\mathrm {min} },\lambda _{\mathrm {max} }]} abgeschätzt werden. Besonders klein wird dieses Maximum, wenn man die Relaxationsparameter so wählt, dass ihre Kehrwerte gerade die Nullstellen des geeignet verschobenen Tschebyschow-Polynoms sind,

2 ω k = λ m i n + λ m a x + ( λ m a x λ m i n ) cos ( 2 k 1 2 m π ) ,   k = 1 , , m . {\displaystyle {\frac {2}{\omega _{k}}}=\lambda _{\mathrm {min} }+\lambda _{\mathrm {max} }+\left(\lambda _{\mathrm {max} }-\lambda _{\mathrm {min} }\right)\cos \left({\frac {2k-1}{2m}}\pi \right),\ k=1,\ldots ,m.}

Dann verbessert sich die Konvergenzaussage für symmetrische Matrizen und einen Zyklus der Länge m {\displaystyle m} zu

ρ ( p m ( A ) ) 2 ( κ 2 ( A ) 1 κ 2 ( A ) + 1 ) m . {\displaystyle \rho (p_{m}(A))\leq 2\left({\frac {{\sqrt {\kappa _{2}(A)}}-1}{{\sqrt {\kappa _{2}(A)}}+1}}\right)^{m}.}

Für realistische Probleme mit κ 2 ( A ) 1 {\displaystyle \kappa _{2}(A)\gg 1} stellt dies eine große Verbesserung gegenüber dem einfachen relaxierten Verfahren dar, da nur noch die Wurzel der Konditionszahl eingeht.

Für symmetrisch-definite Matrizen bietet dieses Verfahren kaum Vorteile gegenüber dem Verfahren der konjugierten Gradienten, da es die Schätzung der Eigenwerte λ m i n , λ m a x {\displaystyle \lambda _{\mathrm {min} },\lambda _{\mathrm {max} }} erfordert. Im unsymmetrischen Fall können aber die Parameter auch für komplexe Eigenwerte gut angepasst werden, vgl. Literatur. In den meisten Fällen ist aber die Tschebyschow-Iteration vorzuziehen, da sie die gleiche Fehlerschranke für jeden Iterationsschritt und nicht nur für Vielfache der Zykluslänge k = m , 2 m , 3 m , {\displaystyle k=m,2m,3m,\ldots } erreicht.

Literatur

  • Andreas Meister: Numerik linearer Gleichungssysteme. 2. Auflage. Vieweg 2005, ISBN 3-528-13135-7
  • Bernd Fischer, Lothar Reichel: A stable Richardson iteration method for complex linear systems. In: Numer. Math. 54. 1988, 225–242.