Implementierung der Differentialoperatoren auf diskreten Skalarfeldern

Hallo,

ich programmiere eine diskrete Simulation, wofuer ich die Differentialoperatoren (grad, div, curl, laplace, ...) brauche. Mir ist klar, wie ich diese auf gewoehnliche differenzierbare mehrdimensionale Funktionen anwenden kann. Nun aber habe ich es mit einem diskreten Skalarfeld zu tun, zum Beispiel so (1D der Einfachheit halber):

+-+-+-+-+
|2|5|0|3|
+-+-+-+-+

Wie saehe es nun aus, wuerde ich den Gradienten dieses Skalarfeldes berechnen wollen? So?

+---+----+---+
|(3)|(-5)|(3)|
+---+----+---+

Ist das korrekt? Es kommt mir ein wenig spanisch vor, da die Anzahl der Elemente pro Differenzierung um 1 reduziert wird.
Wenn ich nun eine DGL a la \frac{\partial}{\partial t} f(\vec{p}, t) = - \operatorname{grad}(f) numerisch simulieren will, wie soll das dann gehen, wenn das rechte Feld 1 Element weniger besitzt als das linke?

LG

torenfroll420 schrieb:

Wenn ich nun eine DGL a la \frac{\partial}{\partial t} f(\vec{p}, t) = - \operatorname{grad}(f) numerisch simulieren will,

Upsi, da ist mir ein Fehler unterlaufen. Hier die Gleichung: \frac{\partial}{\partial t} f(\vec{p}, t) = - \left | \operatorname{grad}(f) \right |

scrontch

Naja, du hast eben das Stützgitter verschoben.
Das ist in gewisser Weise sinnvoll.
Aber die nächste Differenzierung würde wieder auf den ersten Gitterpunkten leben.
Es wird also nicht die Anzahl der Elemente immer um 1 reduziert, sondern nur ungeradzahlige Ableitungen leben auf dem versetzten Gitter mit (n-1) Punkten und ganzzahlige auf dem ursprünglichen Gitter mit n Punkten.
Das Problem kommt letztlich daher dass dein Gitter endlich ist.
Hier brauchst du also Randbedingungen - was ja bei Diff-Gls auch immer angegeben werden muss damit eine Lösung Sinn macht!
Du kannst z.b. die Werte an den Enden für höherwertige Ableitungen fest vorschreiben.

krümelkacker

Wenn Du den "Abtastpunkt" nicht verschieben willst, bietet sich die "zentrale Differenz" an, also

Eingabe: s[k] für diskrete k
Ausgabe a[k] mit a[k] = 0.5 * (s[k+1] - s[k-1])

...und das ist im Prinzip schon ein digitales Filter, welches die Amplituden und Phasen beeinflusst. Den Ableitungsoperator kann man nämlich auch als Filter verstehen, er verschiebt die Phase um 90° und die Amplitude verändert sich in Abhängigkeit der Frequenz.

d sin(f * t) / dt = f * cos(f * t)

Die recht einfache "zentrale Differenz" approximiert den idealen Differentiator aber nicht besonders gut. Für höhere Frequenzen ist dieser Ansatz nicht so gut. Es gibt auch Differentiatoren "höherer Ordnung". Das nächst bessere wäre

a[k] = (s[k-2] - 8 s[k-1] + 8 s[k+1] - s[k+2]) / 12.0

Dafür braucht man aber offensichtlich noch mehr Werte an den Rändern. Das ist ein grundsätzliches Problem.

Was ich aber auch empfehlen kann ist das Ableiten anhand einer eine kubischen B-spline Interpolation. Der vereinfacht sich zu einem "Vorfilter" (bidirektional angewandter IIR-Filter 1. Ordnung) und anschließender Zentraldifferenz. Das Randproblem wird dadurch nicht einfacher, aber die Güte ist super und das ist auch relativ CPU-schonend.

Siehe
"Image Interpolation and Resampling"
http://bigwww.epfl.ch/publications/thevenaz9901.pdf

Hier ein bisschen Matlab-Code für die B-Spline Variante:

eingabe = [0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0]; % Einheitsimpuls als Test
a = [1 2-sqrt(3)]; b = sum(a); % Filter-Parameter für Vorfilter
temp = filtfilt(b,a,eingabe); % vorwärts und rückwärts angewendet
ableitung = conv2(temp,[0.5 0 -0.5],'same'); % Zentrale Differenz

Der vektor "temp" speichert hier die Skalierungskoeffizienten für die kubischen Basis-Splines. Die kann man jetzt dazu benutzen an beliebigen Stellen die Funktion auszuwerten oder ihre Ableitungen, welche bis einschließlich der zweiten immer noch überall stetig ist. Wenn man aber an genau den Stützstellen die Ableitung haben will, kann man das einfach mit [0.5 0 -0.5] falten und ist fertig. Will man die zweite Ableitung, faltet man mit [1 -2 1] statt [0.5 0 -0.5] im zweiten Schritt, wobei derselbe Vorfilter benutzt wird.

Der obige Code sollte auch unter GNU/Octave laufen.

Hier eine Abbildung, die zeigt, was die verschienenen Filter mit verschiedenen Frequenzen anstellen. Der "ideale Differentiator" ist das, was man bekommen würde, wenn man eine perfekte sinc-Interpolation durchführen und die dann ableiten würde... Das ist aber auch selten das, was man machen will.

All das funktioniert auch in höhreren Dimensionen. Man kann die einzelnen Achsen separat handhaben, z.B. erst alle Spalten für die 2. Ableitung nach y filtern und dann alle Zeilen filtern für die 1. Ableitung nach x oder sowas.