Theoretische und methodische Grundlagen des Visual Computing – Partielle Differentialgleichungen

partielle Differentialgleichung: Eine partielle Differentialgleichung (PDE) ist eine Gleichung, die eine multivariate Funktion (z. B. \(u(x, t)\)) in Beziehung mit den partiellen Ableitungen (z. B. \(\partial _x u\), \(\partial _t u\)) und den Variablen (z. B. \(x\), \(t\)) setzt. Die Ordnung einer PDE ist die Ordnung der höchsten partiellen Ableitung. Eine PDE heißt linear, falls \(u\) und die partiellen Ableitungen nur linear auftauchen (wobei die Koeffizienten aber durchaus von den Variablen abhängen dürfen). Analytische Lösungen kann man explizit oft nur für sehr einfache Spezialfälle angeben.

1D-Advektionsgleichung

Advektionsgleichung: Die Gleichung \(\partial _t u + c \partial _x u = 0\) mit \(c \in \real \) zusammen mit Anfangs- und/oder Randbedingungen heißt 1D-Advektionsgleichung. Gesucht ist nun eine Lösung, d. h. eine Funktion \(u(x, t)\), die diese Gleichung erfüllt.

Lösungen konstant entlang Charakteristiken: Sei \(x(t) = ct + x_0\) eine Gerade im
\(t\)-\(x\)-Diagramm mit Steigung \(c\). Dann gilt \(\frac {\d }{\dt } u(x(t), t) = \partial _x u(x(t), t) x’(t) + \partial _t u(x(t), t)\)
\(= (\partial _t u + c \partial _x u)(x(t), t) = 0\) für jede Lösung \(u(x, t)\), d. h. Lösungen der Advektionsgleichung sind konstant entlang diesen sog. Charakteristiken.

Lösung der Advektionsgleichung: Sei \(u(x, t) := v(x - ct)\) für eine Funktion \(v\). Dann gilt \(\partial _t u + c \partial _x u = v’(x - ct) (-c) + c v’(x - ct) = 0\), d. h. \(u\) ist eine Lösung der Advektionsgleichung.

Klassifikation linearer PDEs zweiter Ordnung

lineare PDE zweiter Ordnung: Eine lineare PDE zweiter Ordnung auf einem Gebiet \(D \subset \real ^2\) kann man schreiben als \(A \partial _x^2 u + 2B \partial _x \partial _y u + C \partial _y^2 u + D \partial _x u + E \partial _y u + Fu = G\) mit \(A, \dotsc , G\) stückweise stetigen Funktionen von \(x, y\) auf \(D\).

Klassifikation: Wenn \(A^2 + B^2 + C^2 \not = 0\) in \(D\) ist, dann heißt die PDE

elliptisch, falls \(AC - B^2 > 0\) in \(D\),
parabolisch, falls \(AC - B^2 = 0\) in \(D\), und
hyperbolisch, falls \(AC - B^2 < 0\) in \(D\).

Diese Klassifikation ist analog zu der von Kegelschnitten \(Ax^2 + 2Bxy + Cy^2 + Dx + Ey + F = 0\).

Beispiele:

elliptische PDEs: Laplace-Gleichung (\(\Delta u = 0\)), Poisson-Gleichung (\(\Delta u = f\)),
Helmholtz-Gleichung (\(\Delta u + \varrho u = f\))
parabolische PDEs: 1D-Wärmeleitungsgleichung (\(k \partial _x^2 u - \partial _t u = 0\)),
Diffusionsgleichung (\(\partial _t u(\vec {x}, t) = \div (D(u, \vec {x}) \nabla u(\vec {x}, t))\))
hyperbolische PDE: Wellengleichung (\(\partial _t^2 u - c^2 \Delta u = 0\))

Randbedingungen: Damit die Lösung von elliptischen PDEs eindeutig bestimmt ist, muss man Randbedingungen (RBen) festlegen. Sei dazu \(D\) stückweise \(\C ^1\)-berandet mit den Stücken \(\Gamma _i\), \(\partial D = \Gamma = \bigcup _i \Gamma _i\) und \(\vec {n}\) der nach außen zeigenden Einheitsnormalen. Dann spricht man bei \(u = \varphi \) auf \(\Gamma _i\) von Dirichlet-RBen, bei \(\partial _{\vec {n}} u = \gamma \) auf \(\Gamma _i\) von Neumann-RBen und bei \(\partial _{\vec {n}} u + \alpha u = \beta \) auf \(\Gamma _i\) von Cauchy-RBen (\(\varphi , \gamma , \alpha , \beta \) feste Funktionen auf \(\Gamma _i\)).

Laplace-Gleichung in Polarkoordinaten

Laplace-Gleichung in Polarkoordinaten: Seien \(\Omega := B_1(0) \subset \real ^2\) und \(f \in \C ^0(\partial \Omega )\).
Gesucht ist eine Lösung \(u \in \C ^2(\Omega ) \cap \C ^0(\overline {\Omega })\) des RWPs \(\Delta u = 0\) in \(\Omega \) und \(u = f\) auf \(\partial \Omega \).

Mit Transformation in Polarkoordinaten \(x = r\cos \varphi \), \(y = r\sin \varphi \) gilt \(\partial _x = \cos \varphi \cdot \partial _r - \frac {\sin \varphi }{r} \partial _\varphi \) und \(\partial _y = \sin \varphi \cdot \partial _r + \frac {\cos \varphi }{r} \partial _\varphi \). Durch nochmalige Ableitung erhält man \(\Delta u = \partial _x^2 u + \partial _y^2 u = 0 \iff \partial _r^2 U + \frac {1}{r} \partial _r U + \frac {1}{r^2} \partial _\varphi ^2 U = 0\) mit \(U(r, \varphi ) := u(r\cos \varphi , r\sin \varphi )\). Die RB \(u = f\) auf \(\partial \Omega \) transformiert sich zu \(U(1, \varphi ) = f(\cos \varphi , \sin \varphi )\). Dabei ist \(r \in (0, 1)\) und \(\varphi \in (-\pi , \pi )\).

Damit \(u \in \C ^2(\Omega ) \cap \C ^0(\overline {\Omega })\) gilt, sollte \(U\) zusätzlich \(U(r, \pi ) = U(r, -\pi )\) und
\(\partial _\varphi U(r, \pi ) = \partial _\varphi U(r, -\pi )\) erfüllen und der Grenzwert \(\lim _{r \to 0} U(r, \varphi )\) sollte für jedes \(\varphi \) existieren und unabhängig von \(\varphi \) sein.

Mit dem Produktansatz \(U(r, \varphi ) := w(r) v(\varphi )\) erhält man \(v \partial _r^2 w + \frac {1}{r} v \partial _r w + \frac {1}{r^2} w \partial _\varphi ^2 v = 0\). Wenn man annimmt, dass \(w(r) \not = 0 \not = v(\varphi )\), dann kann man das umschreiben zu
\(\frac {r^2}{w} \partial _r^2 w + \frac {r}{w} \partial _r w = -\frac {1}{v} \partial _\varphi ^2 v =: \lambda \). \(\lambda \) ist unabhängig von \(r\) und \(\varphi \), weil die linke Seite nur von \(r\) abhängt und die rechte nur von \(\varphi \). \(w\) und \(v\) müssen also die ODEs

\(v’’ + \lambda v = 0\) mit \(v(\pi ) = v(-\pi )\) und \(v’(\pi ) = v’(-\pi )\) (harmonischer Oszillator)
\(r^2 w’’ + rw’ - \lambda w = 0\) (Besselsche DGL) und

erfüllen. Aus der ersten DGL folgt, dass \(v\) \(2\pi \)-periodisch sein muss, also \(\lambda = k^2\) für \(k \in \natural _0\) und \(v_k(\varphi ) = a_k \cos (k\varphi ) + b_k \sin (k\varphi )\). Lösungen der zweiten DGL sind \(w(r) = \log r\) für \(k = 0\) und \(w(r) = r^{\pm k}\) für \(k \ge 1\). Allerdings ist nur \(w(r) = r^k\) in \(r = 0\) stetig.

Somit ist \(U(r, \varphi ) = \frac {1}{2} a_0 + \sum _{k=1}^\infty U_k(r, \varphi )\) mit \(U_k(r, \varphi ) := (a_k \cos (k \varphi ) + b_k \sin (k \varphi )) r^k\). Um die Koeffizienten \(a_k\), \(b_k\) zu bestimmen, nutzt man die Orthogonalitätsrelationen der trigonometrischen Funktionen aus (Multiplikation z. B. mit \(\sin (n\varphi )\) und Integration über \(\varphi \) für \(r = 1\)).

1D-Diffusionsgleichung

Diffusionsgleichung: Die Gleichung \(\partial _t u = D \partial _x^2 u\) für \(x \in [0, L]\) mit \(D \in \real \), der Anfangsbedingung \(u(\cdot , 0) = f(\cdot )\) in \([0, L]\) und Dirichlet-Null-RBen \(u(0, t) = u(L, t) = 0\) für \(t > 0\) heißt 1D-Diffusionsgleichung.

analytische Lösung mit Trennung der Veränderlichen:
Mit dem Ansatz \(u(x, t) := X(x) \cdot T(t)\) erhält man \(\frac {1}{DT} \partial _t T = \frac {1}{X} \partial _x^2 X =: -\lambda \). Weil die linke Seite nur von \(t\) und die rechte nur von \(x\) abhängt, ist \(\lambda \) konstant und man bekommt zwei ODEs

\(X’’ + \lambda X = 0\) und
\(T’ + \lambda DT = 0\).

Die allgemeine Lösung der ersten Gleichung ist \(X(x) = c_1 \cos (\sqrt {\lambda } x) + c_2 \sin (\sqrt {\lambda } x)\) für \(\lambda > 0\). Mit den RBen erhält man \(c_1 = 0\) und \(\lambda = \lambda _n = (\frac {n \pi }{L})^2\) für \(n \in \natural \), also \(X(x) = C_n \sin (\frac {n \pi }{L} x)\). Aus der zweiten Gleichung bekommt man damit \(T(t) = B_n e^{-D \lambda _n t}\). Insgesamt erhält man also als allgemeine Lösung der Diffusionsgleichung \(u(x, t) = \sum _{n=1}^\infty A_n \sin (\frac {n \pi }{L} x) \exp (-D (\frac {n \pi }{L})^2 t)\).

Bestimmung der Koeffizienten: Für die Bestimmung der \(A_n\) benutzt man die Identität
\(\int _0^L \sin (\frac {n \pi }{L} \xi ) \sin (\frac {m \pi }{L} \xi ) \dxi = \frac {L}{2} \delta _{n,m}\) für \(n \not = 0 \not = m\). Damit bekommt man
\(\int _0^L f(\xi ) \sin (\frac {n \pi }{L} \xi ) \dxi = \int _0^L u(\xi , 0) \sin (\frac {n \pi }{L} \xi ) \dxi = \sum _{m=1}^\infty A_m \int _0^L \sin (\frac {n \pi }{L} \xi ) \sin (\frac {m \pi }{L} \xi ) \dxi = \frac {L}{2} A_n\)
\(\iff A_n = \frac {2}{L} \int _0^L f(\xi ) \sin (\frac {n \pi }{L} \xi ) \dxi \).

Finite-Differenzen-Methode

Finite-Differenzen-Methode: Bei der Finite-Differenzen-Methode diskretisiert man das Gebiet \(D\) zunächst als gleichmäßiges Gitter mit Gitterweite \(\Delta t\), \(\Delta x\) und \(\Delta y\). Anschließend betrachtet man Approximationen \(u_{i,j} \approx u(x_i, y_j)\) und ersetzt alle Ableitungen durch Differenzenquotienten. Im Folgenden seien \(u_P := u_{i,j}\), \(u_N := u_{i,j+1}\), \(u_E := u_{i+1,j}\), \(u_S := u_{i,j-1}\) und \(u_W := u_{i-1,j}\).

einfache Ableitungen: Sei \(h := \Delta x = \Delta y\).

\(\partial _x u(x_i, y_j) \approx \frac {u_E - u_P}{h}\) (Vorwärts-Diff.quot.), \(\partial _x u(x_i, y_j) \approx \frac {u_P - u_W}{h}\) (Rückwärts-Diff.quot.) oder \(\partial _x u(x_i, y_j) \approx \frac {u_E - u_W}{2h}\) (zentraler Diff.quot.)
\(\partial _y u(x_i, y_j) \approx \frac {u_N - u_P}{h}\), \(\partial _y u(x_i, y_j) \approx \frac {u_P - u_S}{h}\) oder \(\partial _y u(x_i, y_j) \approx \frac {u_N - u_S}{2h}\)

Eine Approximation höherer Ordnung ist \(\partial _x u(x_i, y_j) \approx \frac {u_{WW} - 8u_W + 8u_E - 8u_{EE}}{12h}\).

zweifache Ableitungen: \(\partial _x^2 u \approx \frac {u_E - 2u_P + u_W}{h^2}\), \(\partial _y^2 u \approx \frac {u_N - 2u_P + u_S}{h^2}\), \(\partial _x \partial _y u \approx \frac {u_{NW} - u_{NE} - u_{SW} + u_{SE}}{4h^2}\)
Eine Approximation höherer Ordnung ist \(\partial _x^2 u(x_i, y_j) \approx \frac {-u_{WW} + 16u_W - 30u_P + 16u_E - u_{EE}}{12h^2}\).

FDM für die Poisson-Gleichung: Die Poisson-Gleichung ist gegeben durch \(-\Delta u = f\) mit einem Quellterm \(f\). Die diskretisierte Version ist gegeben durch \(4u_P - u_E - u_W - u_N - u_S = h^2 f_P\). Um Dirichlet-RBen zu erzwingen, setzt man die jeweiligen Terme auf den RB-Wert.
Für Neumann-RBen mit z. B. \(\vec {n} = (1, 0)^\tp \) und \(\partial _{\vec {n}} u(x_i, y_j) = 0\) erhält man die Approximation \(0 = \partial _{\vec {n}} u(x_i, y_j) \approx \frac {u_E - u_W}{2h}\), also \(u_E := u_W\). Daher erzwingt man Neumann-RBen, indem man in der Diskretisierung für den Randpunkt \(P\) den Term \(u_E\) durch \(u_W\) ersetzt, also
\(2u_P - u_W - \frac {1}{2} u_N - \frac {1}{2} u_S = \frac {1}{2} h^2 f_P\) (Teilen durch \(2\), damit LGS-Matrix nachher symmetrisch). Durch Lösen des entstehenden LGS kann man die Approximation ausrechnen.

FDM für Advektionsgleichung: Für die Advektionsgleichung \(\partial _t u + c \partial _x u = 0\) gibt es mehrere Möglichkeiten zur Diskretisierung, je nachdem, welche Differenzenquotienten man wählt. Dazu setzt man \(u_i^n \approx u(x_i, t_n)\) mit Gitterweiten \(\Delta x\) und \(\Delta t\).

Upwind-Methode: \(u_i^{n+1} := u_i^n - c \frac {\Delta t}{\Delta x} (u_i^n - u_{i-1}^n)\) (Zeit vorwärts, Ort rückwärts)
Downwind-Methode: \(u_i^{n+1} := u_i^n - c \frac {\Delta t}{\Delta x} (u_{i+1}^n - u_i^n)\) (Zeit vorwärts, Ort vorwärts)
zentrierte Methode: \(u_i^{n+1} := u_i^n - c \frac {\Delta t}{2 \Delta x} (u_{i+1}^n - u_{i-1}^n)\) (Zeit vorwärts, Ort zentriert)
Leap-Frog-Methode: \(u_i^{n+1} := u_i^{n-1} - c \frac {\Delta t}{\Delta x} (u_{i+1}^n - u_{i-1}^n)\) (Zeit zentriert, Ort zentriert)
Lax-Wendroff-Methode: \(u_i^{n+1} := u_i^n - c \frac {\Delta t}{2 \Delta x} (u_{i+1}^n - u_{i-1}^n) + \frac {1}{2} (c \frac {\Delta t}{\Delta x})^2 (u_{i+1}^n - 2u_i^n + u_{i-1}^n)\)

Für Stabilität sollte die CFL-Zahl (Courant-Friedrichs-Lewy) \(\sigma := c \frac {\Delta t}{\Delta x}\) kleiner als \(1\) sein.

FDM für Diffusionsgleichung: Als Beispiel betrachtet man die Wärmeleitung in einem Stab der Länge \(L = 1\) mit \(D = 1\). Anfangsbedingung soll \(u(x, 0) = 0\) für \(x \in [0, 1]\) und RBen sollen \(u(0, t) = \sin (\pi t)\) und \(\partial _x u(1, t) = 0\) für \(t > 0\). Diskretisiere \(D = [0, 1] \times [0, \infty )\) mit \(x_i := i \Delta x\) und \(t_n := n \Delta t\), wobei \(\Delta x := \frac {1}{N}\), \(i = 0, \dotsc , N\) und \(n \in \natural _0\). Durch Verwendung von Vorwärts-Diff.quot. für die Zeit und zentralen Diff.quot. für den Ort bekommt man mit \(u_i^n \approx u(x_i, t_n)\) die FTCS-Methode \(u_i^{n+1} := u_i^n + \frac {\Delta t}{\Delta x^2} (u_{i+1}^n - 2u_i^n + u_{i-1}^n)\) mit RBen \(u_0^n := \sin (\pi n \Delta t)\) und \(u_N^{n+1} := u_N^n + \frac {2\Delta t}{\Delta x^2} (u_{N-1}^n - u_N^n)\). Das Verfahren ist stabil \(\iff \Delta t < \frac {1}{2} (\Delta x)^2\).

Crank-Nicolson-Methode

Crank-Nicolson-Methode: Die Crank-Nicolson-Methode verwendet die Trapezregel für die Zeit und zentrale Diff.quot. für den Ort. Ist die Gleichung \(\partial _t u = F(x, t, u, \partial _x u, \partial _x^2 u)\) gegeben und bezeichnen \(\frac {u_i^{n+1} - u_i^n}{\Delta t} = F_i^n(x, t, u, \partial _x u, \partial _x^2 u)\) und \(\frac {u_i^{n+1} - u_i^n}{\Delta t} = F_i^{n+1}(x, t, u, \partial _x u, \partial _x^2 u)\) einen expliziten bzw. impliziten Euler-Schritt, so erhält man die Crank-Nicolson-Methode durch Mittelung \(\frac {u_i^{n+1} - u_i^n}{\Delta t} = \frac {1}{2} (F_i^n + F_i^{n+1})(x, t, u, \partial _x u, \partial _x^2 u)\).

Crank-Nicolson-Methode für Diffusionsgleichung:
Man setzt \(\partial _x^2 u(x_i, t_n) \approx \frac {1}{2(\Delta x)^2} ((u_{i+1}^n - 2u_i^n + u_{i-1}^n) + (u_{i+1}^{n+1} - 2u_i^{n+1} + u_{i-1}^{n+1}))\) und erhält damit
\(-r u_{i-1}^{n+1} + (2+2r)u_i^{n+1} - ru_{i+1}^{n+1} = ru_{i-1}^n + (2-2r)u_i^n + ru_{i+1}^n\) mit \(r := \frac {\Delta t}{(\Delta x)^2}\).
Die RBen lauten dabei \((2+2r)u_1^{n+1} - ru_2^{n+1} = (2-2r)u_1^n + ru_2^n + r(\sin (\pi n \Delta t) + \sin (\pi (n+1) \Delta t)\) und \(-2ru_{N-1}^{n+1} + (2+2r)u_N^{n+1} = 2ru_{N-1}^n + (2-2r)u_N^n\). Für jede Zeit \(t_n\) erhält man also ein tridiagonales LGS mit \(N\) Gleichungen für die \(N\) Unbekannten \(u_1^{n+1}, \dotsc , u_N^{n+1}\).

Anisotrope 1D-Diffusionsgleichung

anisotrope 1D-Diffusionsgleichung:
Die Gleichung \(\partial _t u = \partial _x (D(x) \partial _x u)\) heißt anisotrope 1D-Diffusionsgleichung.

FDM für anisotrope 1D-Diffusionsgleichung: Diskretisiere zunächst die rechte Seite
\(\partial _x g(x)\) mit \(g(x) := D(x) \partial _x u\) und halber Schrittweite durch \(\partial _x g(x) \approx \frac {g(x + h/2) - g(x - h/2)}{h}\). Approximiere nun \(g(x) \approx D(x) \frac {u(x + h/2) - u(x - h/2)}{h}\). Durch Einsetzen erhält man
\(\frac {u_i^{n+1} - u_i^n}{\Delta t} = \frac {1}{\Delta x^2} (D_{i+1/2} (u_{i+1}^n - u_i^n) - D_{i-1/2} (u_i^n - u_{i-1}^n))\) mit \(D_{i\pm 1/2} := D(x_i \pm \Delta x/2)\).

Crank-Nicolson für anisotrope 1D-Diffusionsgleichung:
Durch Mittelung für die Schritte \(n\) und \(n + 1\) erhält man
\(\frac {u_i^{n+1} - u_i^n}{\Delta t} = \frac {1}{2\Delta x^2} (D_{i+1/2} (u_{i+1}^{n+1} + u_{i+1}^n - u_i^{n+1} - u_i^n) - D_{i-1/2} (u_i^{n+1} + u_i^n - u_{i-1}^{n+1} - u_{i-1}^n))\) bzw.
\(-D_{i+1/2} u_{i+1}^{n+1} + (\frac {2\Delta x^2}{\Delta t} + D_{i+1/2} + D_{i-1/2}) u_i^{n+1} - D_{i-1/2} u_{i-1}^{n+1}\)
\(= D_{i+1/2} u_{i+1}^n (\frac {2\Delta x^2}{\Delta t} - D_{i+1/2} - D_{i-1/2}) u_i^n + D_{i-1/2} u_{i-1}^n\).
Dies kann man in der Form \(A\vec {u}^{n+1} = B\vec {u}^n\) schreiben mit \(\vec {u}^n := (u_1^n, \dotsc , u_N^n)^\tp \),
\(A_{j,j-1} := -D_{j-1/2}\), \(A_{j,j} := \frac {2\Delta x^2}{\Delta t} + D_{j+1/2} + D_{j-1/2}\), \(A_{j,j+1} := -D_{j+1/2}\) und
\(B_{j,j-1} := D_{j-1/2}\), \(B_{j,j} := \frac {2\Delta x^2}{\Delta t} - D_{j+1/2} - D_{j-1/2}\), \(B_{j,j+1} := D_{j+1/2}\) (sonst \(A_{i,j} := 0 =: B_{i,j}\)).
Nach Hinzufügung der RBen zu \(A\) und \(B\) erhält man dann \(\vec {u}^{n+1} = A^{-1} B\vec {u}^n\).

Perona-Malik-Diffusion

Perona-Malik-Diffusion: Die Perona-Malik-Diffusionsgleichung ist gegeben durch
\(\partial _t u = \div (g(|\vec {\nabla } u|) \vec {\nabla } u)\) mit AB \(u(\cdot , 0) = f\) und homogenen Neumann-RBen \(\partial _{\vec {n}} u = 0\).
Die Diffusivitätsfunktion \(g\) sollte \(g(0) = 1\), \(g \ge 0\) und \(\lim _{s \to \infty } g(s) = 0\) erfüllen. PM schlagen z. B. \(g(s) := \frac {1}{1 + s^2/\lambda ^2}\) oder \(g(s) := e^{-s^2/\lambda ^2}\) mit dem Kontrastparameter \(\lambda \) vor.

1D-PM-Diffusion: In einer Dimension ist die PM-Diffusion mit der Flussfunktion
\(\Phi (s) := s \cdot g(s)\) gegeben durch \(\partial _t u = \partial _x (g(|\partial _x u|) \partial _x u) = \Phi ’(|\partial _x u|) \partial _x^2 u\). Für die beiden von PM vorgeschlagenen Diffusivitätsfunktionen gilt \(\Phi ’(|\partial _x u|) > 0\) für \(|\partial _x u| < \lambda \) (Vorwärts-Diffusion) und \(\Phi ’(|\partial _x u|) < 0\) für \(|\partial _x u| > \lambda \) (Rückwärts-Diffusion). Obwohl also die Diffusivität \(g\) nicht-negativ ist, gibt es Vor- und Rückwärts-Diffusion.

2D-PM-Diffusion: In zwei Dimensionen ist die PM-Diffusion gegeben durch
\(\partial _t u = \partial _x (g(|\vec {\nabla } u|) \partial _x u) + \partial _y (g(|\vec {\nabla } u|) \partial _y u)\). Die Norm des Gradienten kann man mittels zentraler Diff.quot.en approximieren: \(|\vec {\nabla } u_{i,j}| \approx \sqrt {(\frac {u_{i+1,j} - u_{i-1,j}}{2})^2 + (\frac {u_{i,j+1} - u_{i,j-1}}{2})^2}\) mit \(g_{i,j} := g(|\vec {\nabla } u_{i,j}|)\).
Damit kann man ein Finites-Differenzen-Verfahren herleiten.
Eine Anwendung ist die Rauschunterdrückung von Bildern, ohne dass Kanten unscharf werden (im Gegensatz zur linearen Diffusion).

Dilatation und Erosion

Dilatation und Erosion: Bei der Dilatation und Erosion ist eine Teilmenge \(A \subset \real ^2\) (z. B. ein Quadrat) und eine Maske \(M \subset \real ^2\) mit \(0 \in M\) (z. B. ein Kreis). Die Dilatation ist das Ergebnis von \(\bigcup _{x \in A} (M + x)\) und die Erosion ist das Ergebnis von \(\{x \in A \;|\; M + x \subset A\}\).

Dilatationsgleichung: Dilatation und Erosion können mithilfe der Dilatationsgleichung
\(\partial _t u = \pm |\vec {\nabla } u|\) simuliert werden. Bei der Rouy-Tourin-Methode setzt man
\(u_{i,j}^{n+1} := u_{i,j}^n \pm \Delta t \sqrt {(\partial _x u)^2 + (\partial _y u)^2}\) mit \((\partial _x u)^2 \approx \frac {1}{\Delta x^2} (\max \{0, u_{i+1,j} - u_{i,j}, -(u_{i,j} - u_{i-1,j})\})^2\), analog \((\partial _y u)^2\).
Bei der Osher-Sethian-Upwind-Methode verwendet man stattdessen
\((\partial _x u)^2 \approx \frac {1}{\Delta x^2} ((\max (0, u_{i+1,j} - u_{i,j}))^2 + (\max (0, u_{i-1,j} - u_{i,j}))^2)\), analog \((\partial _y u)^2\).
Anwendungen von Dilatation und Erosion umfassen das Unscharfmachen von Kanten oder die Kontrasterhöhung von Fingerabdrucks-Bildern.