Das Differenzenverfahren

Als erster Schritt wird der Funktionsraum diskretisiert, indem man über das Grundgebiet ein zweidimensionales quadratisches Punktgitter (Gitterkonstante h) legt.

Nun werden alle Gitterpunkte, die im Inneren des Grundgebietes liegen (innere Gitterpunkte oder Netzpunkte) geeignet durchnumeriert. Der i -te von insgesamt i_max Gitterpunkten habe die Koordinaten x_i;y_i. Nun werden an jedem inneren Gitterpunkt die Differential-Quotienten in der Laplace-Gleichung durch entsprechende Differenzen-Quotienten ersetzt. Im gegebenen Programm werden dafür die einfachen "Dreipunkt-Formelnbeginequation &part#partial;^2 u&part#partial;x^2 &mid#mid;_x_i;y_i &ap#approx;u_i,r-2u_i+u_i,lh^2 beziehungsweise

$\displaystyle {\frac{{\partial^{2} u}}{{\partial y^{2}}}}$$ $\displaystyle \mid_{{x_{i};y_{i}}}^{}$$ $\displaystyle \approx$$ $\displaystyle {\frac{{u_{i,o}-2u_{i}+u_{i,u}}}{{h^{2}}}}$$ .

(3)

Setzt man (

) und (

) in die Dgl. (

) ein, so erhält man für den i-ten Punkt die Differenzengleichung

4u_i - u_{i, l} - u_{i, r} - u_{i, o} - u_{i, u} = 0 (i = 1,..., i_max) .

(4)

Die Gesamtheit aller Differenzengleichungen bildet ein homogenes, lineares Gleichungssystem. Inhomogene Elemente kommen über jene Gitterpunkte ins Spiel, die nicht nach allen 4 Seiten einen Gitterpunkt als Nachbarn haben. Solche fehlenden Nachbarpunkteliegen entweder auf der Randkurve, wenn das Gitter genau in das Grundgebiet hineinpaßt, wie das entlang der geradlinigen Begrenzungen (Punkte 5 $\rightarrow$$ 12 $\rightarrow$$ 1) der Fall ist, oder sie liegen außerhalb des Grundgebietes, z. B. außerhalb der Kreislinie (Punkte 1 $\rightarrow$$ 5). In letzterem Fall wird der Funktionswert an den fehlenden Nachbarpunktenünter Verwendung der nächstliegenden Randwerte linear extrapoliert.

In jedem Fall führen fehlende Nachbarpunkte zu inhomogenen Differenzengleichungen. Insgesamt erhält man also für die u_i ein inhomogenes, lineares Gleichungssystem mit einer für Differenzenverfahren charakteristischen schwach-besetzten Koeffizientenmatrix mit bandförmiger Struktur. Zur numerischen Auswertung solcher Systeme eignet sich ganz besonders das Gauss-Seidel- Verfahren.

Das Testbeispiel

$\displaystyle {\frac{{\partial^{2} u}}{{\partial x^{2}}}}$$ + $\displaystyle {\frac{{\partial^{2} u}}{{\partial y^{2}}}}$$ = 0

u(x, y) $\displaystyle \mid_{{\Gamma: A \to B}}^{}$$ = x²y²

u(x, y) $\displaystyle \mid_{{\Gamma: B \to C}}^{}$$ = 32 - $\displaystyle {\frac{{x^{4}}}{{8}}}$$

u(x, y) $\displaystyle \mid_{{\Gamma: C \to D}}^{}$$ = 32 - $\displaystyle {\frac{{y^{4}}}{{8}}}$$

u(x, y) $\displaystyle \mid_{{\Gamma: D \to E}}^{}$$ = 12x² - $\displaystyle {\frac{{x^{4}}}{{8}}}$$

u(x, y) $\displaystyle \mid_{{\Gamma: E \to A}}^{}$$ = 12y² - $\displaystyle {\frac{{y^{4}}}{{8}}}$$

u_exakt(x, y) = x²y² + $\displaystyle {\frac{{1}}{{8}}}$$ $\displaystyle \left\lbrack\vphantom{ 256-\left(x^{2}+y^{2} \right)^{2} }\right.$$ 256 - $\displaystyle \left(\vphantom{x^{2}+y^{2} }\right.$$ x² + y² $\displaystyle \left.\vphantom{x^{2}+y^{2} }\right)^{{2}}_{}$$ $\displaystyle \left.\vphantom{ 256-\left(x^{2}+y^{2} \right)^{2} }\right\rbrack$$