Finite-Differenzen-Methode

Grundidee und Einsatz

Die Finite-Differenzen-Methode (FDM), auch Differenzenverfahren, ist eine Klasse numerischer Verfahren zur Lösung partieller Differentialgleichungen (DGL). Sie ersetzt Orts- oder Zeitableitungen an Gitterpunkten eines vorgegebenen Intervalls durch Differenzenquotienten. Dadurch wird aus der Differentialgleichung ein Gleichungssystem, dessen Lösung die gesuchten Funktionswerte näherungsweise an den Gitterpunkten liefert.

Die Methode wird unter anderem für fluiddynamische Simulationen, in der Meteorologie, Astrophysik und Baustatik verwendet. In der Reaktorphysik dominierte sie von 1950 bis 1980 bei Programmen zur Berechnung des Neutronenflusses. Ein Beispiel ist ein Zylinderproblem in Zylinderkoordinaten: Wegen Unabhängigkeit von der polaren Koordinate wird ein reales 3D-Problem zu einem 2D-Problem in (r,z). Symmetrie kann Gitterpunkte einsparen; zudem können Randbedingungen verlangen, dass die Lösung an Außenrändern verschwindet.

Eindimensionale Randwertprobleme

Für die gewöhnliche Differentialgleichung zweiter Ordnung

u''(x)=g(x)

soll u im Intervall [a,b] bestimmt werden. Die gegebene Funktion g heißt Inhomogenität. Erst zwei Anfangs- oder Randbedingungen machen die Lösung eindeutig; für Randbedingungen werden finite Differenzen besonders verwendet.

Bei n inneren äquidistanten Stützstellen gibt es n+1 Intervalle und n+2 Stützstellen insgesamt. Die Gitterweite ist

h=(b-a)/(n+1), x_i=a+i·h für i=0,…,n+1.

An einem inneren Gitterpunkt wird die zweite Ableitung durch den zentralen Drei-Punkte-Differenzenquotienten angenähert:

u''(x) ≈ (u(x-h)-2u(x)+u(x+h))/h².

Daraus folgen für i=1,…,n die Differenzengleichungen

(1/h²)(u_{i-1}-2u_i+u_{i+1})=g_i.

Bei vorgegebenen Randwerten u(a) und u(b) werden diese Werte in die erste beziehungsweise letzte Gleichung eingesetzt. Es entsteht ein lineares Gleichungssystem für u_1,…,u_n. Seine Koeffizientenmatrix hat auf der Hauptdiagonalen 2 und auf den beiden Nebendiagonalen -1, nachdem mit -h² multipliziert wurde. Sie ist regulär, dünnbesetzt, eine Bandmatrix und Tridiagonalmatrix; bei konstanten Diagonalelementen auch eine Tridiagonal-Toeplitz-Matrix. Statt n² Speicherplätzen für eine voll besetzte Matrix benötigt eine Tridiagonalmatrix nur 3n.

Gleichmäßige und ungleichmäßige Gitter

Nichtäquidistante Stützstellen sind sinnvoll, wenn die Lösung Besonderheiten wie Singularitäten oder Grenzschichten hat. Für

a=x_0<x_1<…<x_{n+1}=b,

setzt man h_i=x_i-x_{i-1} und h=max_i h_i. Am Punkt x_i bezeichnet h_i die Schrittweite davor und h_{i+1} diejenige danach. Eine zentrale Drei-Punkte-Diskretisierung lautet

(-2/(h_i+h_{i+1}))·((u_i-u_{i-1})/h_i+(u_i-u_{i+1})/h_{i+1})=g_i.

Die daraus entstehende Matrix ist ebenfalls tridiagonal, aber bei unterschiedlichen Schrittweiten keine Tridiagonal-Toeplitz-Matrix. Werden alle h_i gleich h gesetzt, erhält man wieder

-u_{i-1}+2u_i-u_{i+1}=-g_i h².

Der Konsistenzfehler der nichtäquidistanten Formel ist

δ=(1/3)u'''(x_i)(h_{i+1}-h_i)+O(h²).

Auf einem äquidistanten Gitter verschwindet der erste Term; der Konsistenzfehler ist dann von Ordnung 2. Auf beliebigen Gittern ist er von Ordnung Eins, der Fehler der Lösung in der Maximumnorm bleibt jedoch von Ordnung Zwei:

|u(x_i)-u_i|=O(h²).

Für das Randwertproblem u''(x)=2 auf 0≤x≤1 mit u(0)=u(1)=3 ist die exakte Lösung u(x)=3+x(x-1). Bei vier äquidistanten Intervallen mit h=0,25 ergeben sich u_1=2,8125, u_2=2,7500 und u_3=2,8125. Bei den nichtäquidistanten Schrittweiten h_1=1/6, h_2=1/3, h_3=1/3, h_4=1/6 ergeben sich u_1=2,86111, u_2=2,75000 und u_3=2,86111. In beiden Fällen stimmen die Werte an den Gitterpunkten exakt mit der quadratischen Lösungsfunktion überein; dies ist ein Ausnahmefall, weil der Konsistenzfehler hier Null ist.

Upwind-Verfahren und zeitabhängige Wärmeleitung

Beim Konvektions-Diffusionsproblem

-εu''+u'=0, u(0)=0, u(1)=1, 0<ε<<1,

dominiert wegen des kleinen ε der Konvektionsterm u'. Die zentrale Diskretisierung verwendet für u' den Quotienten (u_{i+1}-u_{i-1})/(2h). Bei h>2ε wird der Faktor r=-(2ε+h)/(h-2ε) negativ; die diskrete Lösung oszilliert dann und nähert die exakte Lösung sehr schlecht an.

Das Upwind-Verfahren ersetzt die zentrale Näherung von u' abhängig von der Stromrichtung durch einen einseitigen Differenzenquotienten. Hier lautet die Formel

-ε(u_{i+1}-2u_i+u_{i-1})/h²+(u_i-u_{i-1})/h=0.

Sie liefert auch bei h>2ε eine robuste Approximation. Im eindimensionalen Fall sorgt sie für positive Diagonalelemente und negative Außerdiagonalelemente der Koeffizientenmatrix, was Stabilität garantiert.

Für die Wärmeleitungsgleichung ∂_t u-Δu=f auf (0,T)×Ω mit u=0 auf ∂Ω und u(0,·)=u_0 wird in 1D Ω=(a,b) gewählt. Mit N+2 äquidistanten Punkten und h=(b-a)/(N+1) werden nur die N inneren Werte im Vektor U_h=(u_1(t),…,u_N(t))^T betrachtet. Die Ortsableitung wird durch

∂{xx}u(t,x_i) ≈ (u(t,x{i+1})-2u(t,x_i)+u(t,x_{i-1}))/h²

ersetzt. Es entsteht das System gewöhnlicher Differentialgleichungen erster Ordnung

Ů_h(t)=-A_hU_h+f_h,

wobei A_h=(1/h²)·tridiag(-1,2,-1) ist. Dieses System kann etwa mit dem Euler- oder Runge-Kutta-Verfahren gelöst werden.

Poissongleichung, Fehler und Gleichungsklassen

Für die zweidimensionale Poissongleichung -Δu=f in Ω=(0,1)² mit u|_{∂Ω}=0 wird ein Gitter mit x_i=ih, y_j=jh und h=1/(n+1) verwendet. An inneren Punkten gilt der Fünf-Punkte-Stern:

(1/h²)(4u_{ij}-u_{i-1,j}-u_{i+1,j}-u_{i,j-1}-u_{i,j+1})=f(x_i,y_j).

Das Gleichungssystem besitzt eindeutig eine Lösung. Ist die exakte Lösung ausreichend glatt, gilt |u(x_{ij})-u_{ij}|=O(h²); in Gebieten mit Ecken ist diese Voraussetzung nicht trivial.

Allgemein hat eine FDM die Form A_hu_h=b_h. U_i=u(x_i) beschreibt die exakte Lösung an den Gitterpunkten. Konsistenz der Ordnung n bedeutet

||A_hU-b_h||≤C_Kh^n.

Stabilität bedeutet für alle Gitterfunktionen v:

||A_hv||≥C_S||v||.

Dann folgt Konvergenz derselben Ordnung:

||U-u_h||≤(C_K/C_S)h^n.

Als Norm wird häufig die Maximumnorm ||v||=max_i|v(x_i)| verwendet. Stabilität vorausgesetzt, ist Konsistenz der Ordnung n für Konvergenz der Ordnung n hinreichend, aber nicht notwendig.

Partielle Differentialgleichungen zweiter Ordnung werden über Δ(x,y)=B²(x,y)-4A(x,y)C(x,y) klassifiziert: hyperbolisch bei Δ>0, parabolisch bei Δ=0 und elliptisch bei Δ<0. Elliptische Probleme erhalten typischerweise Randbedingungen, parabolische und hyperbolische Probleme Anfangs- und Randbedingungen. Beispiele sind Laplace-, Poisson-, Helmholtz-, Wellen- und Wärmeleitungsgleichung.