Erzeugen von Höhenlinien zu zweidimensionalen Funktionen mit contour()

Einordnung des Artikels

• Einführung in die Informatik
  • Einführung in die Programmiersprache R
    • Graphik
      • Erzeugen von dreidimensionalen Graphiken mit scatterplot3d
      • Erzeugen von Höhenlinien zu zweidimensionalen Funktionen mit contour()

Hier werden Kenntnisse vorausgesetzt, wie man mit Hilfe von outer(X, Y, FUN = "*", ...) Matrizen erzeugt. Die Funktion outer() wird ausführlich besprochen in Die Familie der apply-Funktionen in R Teil 2: Die Verarbeitung mehrerer Listen mit mapply(), Map() und outer().

Einführung

Um reellwertige Funktionen darzustellen, die auf einem zweidimensionalen Gebiet definiert sind, kann man entweder persp() verwenden oder scatterplot3d(). Erstere befindet sich im Paket The graphics package, Letztere wurde in Erzeugen von dreidimensionalen Graphiken mit scatterplot3d ausführlich erklärt. Hier werden derartige "Gebirge" mit Hilfe von persp() erzeugt (siehe Abbildung 1, 4, 7).

Manchmal ist es aber schwierig diese perspektivischen Darstellungen zu interpretieren. In diesen Fällen bietet es sich an:

• geeignete Höhenlinien auszuwählen und
• diese im Definitionsbereich darzustellen.

Die Bilder der Höhenlinien lassen sich oft leichter und vor allem eindeutig interpretieren:

• Entlang einer Höhenlinie ist die Funktion konstant. Die Richtung des steilsten Anstiegs ist immer senkrecht zur Höhenlinie.
• Dort wo die Höhenlinien sehr dicht liegen, steigt die Funktion sehr stark an.
• Liegen die Höhenlinien weit auseinander, verändern sich die Funktionswerte nur sehr wenig.

Im Paket The graphics package gibt es die Funktion contour(), mit der man Diagramme von Höhenlinien erstellen kann. Im Folgenden werden an einigen Beispielen die Eingabewerte von contour() erklärt.

Ein erstes Beispiel

Das folgende Beispiel soll die wichtigsten Eingabewerte der Funktion contour() vorstellen.

Betrachtet wird zunächst die Funktion

S(x, y) = ln x + 1.5 ln y,

wobei x und y nur positive Werte annehmen dürfen. Dabei ist ln der natürliche Logarithmus. Die Bezeichnung S wurde gewählt, weil diese Funktion (bis auf physikalische Konstanten) die Entropie S eines idealen einatomigen Gases in Abhängigkeit von Volumen V (hier x) und Temperatur T (hier y) beschreibt.

Eine mögliche Implementierung von S(x, y) zeigt das folgende Skript:

S <- function(x, y){     # x: V, y: T
    stopifnot(x > 0, y > 0)
    return( 1.5 * log(y) + log(x) )
}

Die graphische Darstellung der Funktion S(x, y) kann zum Beispiel mit Hilfe von persp() erfolgen, wobei hier die x- und y-Werte zwischen 0 und 2 gewählt werden:

N = 50
xs <- seq(0.01, 2, length.out = N)

persp(x = xs,
      y = xs,
      z = outer(X = xs, Y = xs, FUN = S),
      col = "red",
      xlab = "x", ylab = "y", zlab = "S(x, y)", 
      phi = 30, theta = 40,
      main = "S (x, y)",
      shade = 0.15, ticktype = "detailed")

Man definiert dazu ein Gitter aus N·N Punkten auf dem Definitionsbereich mit Hilfe des Vektors xs (Zeile 2).

Für jeden Gitterpunkt wird mit outer() der Funktionswert von S berechnet (Zeile 6). Der Rückgabewert des Aufrufs von outer() ist somit eine N×N-Matrix, deren Werte von persp() verwendet werden, um das Diagramm zu erstellen.

Die weiteren Anweisungen sollten sich von selbst erklären, da sie lediglich die Darstellung des Diagramms betreffen. Das Ergebnis ist in Abbildung 1 zu sehen.

Abbildung 1: Darstellung der Funktion S(x, y) als "Gebirge" über dem Definitionsbereich, wobei die x- und y-Werte zwischen 0 und 2 liegen.

Verfolgt man die Funktionswerte entlang der Koordinatenachsen, so erkennt man deutlich den logarithmischen Verlauf.

Aufschlussreicher ist manchmal eine explizite Darstellung der Höhenlinien. In der folgenden Abbildung 2 werden zu den Funktionswerten -5, -4,..., 0, 1 die Höhenlinien in den Definitionsbereich eingezeichnet.

Abbildung 2: Darstellung der Höhenlinien der Funktion S(x, y) zu den Funktionswerten S(x, y) = -5, -4,..., 0, 1.

Da die Funktion nicht symmetrisch in den Variablen x und y ist, ist die Winkelhalbierende des 1. Quadranten (grün eingetragen) keine Symmetrieachse für die Höhenlinien. Auch an den Höhenlinien ist es immer noch schwer zu erkennen, dass der Anstieg in y-Richtung steiler ist als in x-Richtung.

Das Skript zum Erzeugen von Abbildung 2 lautet:

xs <- seq(0.0001, 2, length.out = N)

contour(x = xs, y = xs,
        z = outer(X = xs, Y = xs, FUN = "S"),
        levels = (-5:1),
        col = heatColors(7),
        main = "Höhenlinien von S(x,y)",
        lwd = 2,
        asp = 1)

Wie bei der graphischen Darstellung des "Gebirges" mit persp() wird ein Gitter auf dem Definitionsbereich verwendet, um die mit outer() erzeugte Matrix auszuwerten:

• Dazu definiert Zeile 1 den Vektor xs, der die Gitterpunkte hervorbringen wird.
• In Zeile 3 wird xs an die Eingabewerte x und y von contour() übergeben.
• In Zeile 4 werden die Funktionswerte von S(x, y) auf dem Gitter mit Hilfe von outer() berechnet.

Mit levels = (-5:1) werden die Niveaus der Höhenlinien festgelegt (Zeile 5).

Zeile 9: Mit dem Argument asp wird das Verhältnis der Skalierungen der Achsen festgelegt, genauer das Verhältnis der Skalierung in y-Richtung zu der in x-Richtung. Die Erklärung von asp findet man im Paket The graphics package unter plot.window.

Die anderen Eingabewerte beziehen sich auf die graphische Darstellung und werden später erklärt beziehungsweise erklären sich selbst. Die Funktion heatColors() stammt nicht aus den R-Standardpaketen, auch sie wird weiter unten erklärt.

Die Eingabewerte von contour()

Das folgende Skript zeigt die Eingabewerte der Funktion contour():

contour(x, ...)

## Default S3 method:
contour(x = seq(0, 1, length.out = nrow(z)),
            y = seq(0, 1, length.out = ncol(z)),
            z,
            nlevels = 10, levels = pretty(zlim, nlevels),
            labels = NULL,
            xlim = range(x, finite = TRUE),
            ylim = range(y, finite = TRUE),
            zlim = range(z, finite = TRUE),
            labcex = 0.6, drawlabels = TRUE, method = "flattest",
            vfont, axes = TRUE, frame.plot = axes,
            col = par("fg"), lty = par("lty"), lwd = par("lwd"),
            add = FALSE, ...)

Dabei ist Zeile 1 so zu verstehen, dass man an x eine Liste übergibt, welche die Komponenten x$x , x$y und x$z hat. In der zweiten (S3) Version werden die drei Komponenten explizit gesetzt.

Im ersten Beispiel oben wurde contour() verwendet, um ein Diagramm zu erzeugen, in das dann die Höhenlinien eingetragen werden. Da contour() nur eine beschränkte Kontrolle über das Aussehen des Diagramms erlaubt, empfiehlt es sich:

• zuerst ein Diagramm mit plot() zu erzeugen,
• in das nachträglich die Höhenlinien mit contour() eingetragen werden.

Dazu muss das Argument add = TRUE gesetzt werden. Beispiele dazu werden weiter unten gezeigt.

Die Eingabewerte x, y, z

Mit den Eingabewerten x und y wird ein Gitter im Definitionsbereich erzeugt. Bei der Anwendung von persp() (zum Erzeugen des Graphen der Funktion S oben) hat man bereits gesehen, wie das Gitter eingesetzt wird, um die Funktionswerte zu berechnen und darzustellen.

Ganz ähnlich ist der Mechanismus wenn Höhenlinien erzeugt werden sollen. Dazu wird nochmals das Skript gezeigt, mit dem Abbildung 2 erzeugt wurde:

xs <- seq(0.0001, 2, length.out = N)

contour(x = xs, y = xs,
        z = outer(X = xs, Y = xs, FUN = "S"),
        levels = (-5:1),
        col = heatColors(7),
        main = "Höhenlinien von S(x,y)",
        lwd = 2,
        asp = 1)

Zunächst könnte man annehmen, dass die Eingabewerte x und y nicht nötig sind, da man die Gitterpunkte später an die Funktion outer() übergibt. Dazu muss man wissen, dass man die Eingabewerte x und y tatsächlich weglassen kann. Dann wird aber für den Definitionsbereich das Rechteck [0; 1]×[0; 1] gewählt. Die Eingabewerte x und y definieren somit nicht nur das Gitter, sondern legen zugleich den Definitionsbereich fest.

Dass das Rechteck [0; 1]×[0; 1] gewählt wird, erkennt man an den default-Werten für x und y: x = seq(0, 1, length.out = nrow(z)) und y = seq(0, 1, length.out = ncol(z))

Die Eingabewerte x und y erwarten jeweils einen Vektor, aus denen das Gitter aufgebaut wird (Zeile 3); der Eingabewert z erwartet eine Matrix, die hier mit Hilfe der Funktion outer() berechnet wird (Zeile 4). Auf allen Gitterpunkten wird der Funktionswert von S berechnet und diese Werte werden verwendet, um die Höhenlinien zu zeichnen.

Der Vektor xs unterteilt das Intervall [0; 2] in 49 (es war oben N = 50 gesetzt) Teilintervalle identischer Länge. Als Start kann aber nicht 0 gewählt werden, da die Funktion S für x = 0 beziehungsweise y = 0 nicht definiert ist.

Die folgende Abbildung zeigt das im letzten Skript erzeugte Gitter sowie die Höhenlinien. Dazu werden zuerst mit plot() das Koordinatensystem und das Gitter gezeichnet, anschließend werden mit contour() die Höhenlinien von S(x, y) hinzugefügt.

Abbildung 3: Darstellung des Gitters ('''blau''') und der Höhenlinien der Funktion S(x, y) zu den Funktionswerten -5, -4,..., 0, 1.

Man erkennt, dass die Höhenlinien einen sehr glatten Verlauf haben; nur dort, wo sie stark gekrümmt sind, sieht man wie sie aus Strecken zusammengesetzt werden. Um die Darstellung zu verbessern, muss man im folgenden Skript lediglich die Zahl N erhöhen:

N <- 50
v <- seq(0.0001, 2, length.out = N)
xs <- rep(v, times = N)
ys <- rep(v, each = N)

plot(x = xs, y = ys,
     col = "blue", 
     pch = 20,
     lwd = 1,
     xlab = "x", ylab = "y",
     main = "Gitter, Höhenlinien von S(x, y)",
     frame.plot = TRUE)
grid()
# Koord.achsen:
abline(h = 0, col = "darkgray")
abline(v = 0, col = "darkgray")

contour(x = v, y = v,
        z = outer(X = v, Y = v, FUN = "S"),
        levels = (-5:1),
        col = heatColors(7),
        lwd = 3,
        add = TRUE)

Zur Erklärung:

Zeile 1: Das Gitter ist hier quadratisch. Die Anzahl N der Gitterpunkte in jeder Richtung wird als Variable definiert. Damit kann das Skript leicht abgeändert werden, wenn ein feineres Gitter verwendet werden soll.

Zeile 2: Im Intervall von 0.0001 bis 2 werden 50 Punkte im gleichen Abstand festgelegt. Zeile 2 und 3: Zum Zeichnen aller Gitterpunkte benötigt man ihre x- und ihre y-Koordinaten. Sie können leicht mit rep() erzeugt werden, indem man einmal times und einmal each für die Wiederholungen verwendet.

Zeile 6 bis 16: Die Funktion plot() wird verwendet, um das Koordinatensystem und die Gitterpunkte zu zeichnen. Man beachte, dass xs und ys nur zur Darstellung des Gitters benötigt wird, nicht für die Höhenlinien.

Zeile 18 bis 23: An die Funktion contour() wird der Vektor v aus Zeile 2 an die Argumente x und y übergeben. Zusätzlich wird v verwendet, um mit outer() die Matrix der Funktionswerte von S auf den Gitterpunkten zu berechnen (Eingabewert z von contour()).

Zeile 20: Der Eingabewert levels legt fest, für welche Funktionswerte von S die Höhenlinien gezeichnet werden sollen; hier sind es die ganzen Zahlen von -5 bis +1.

Zeile 21: Um die Farbe der Höhenlinien festzulegen, wird hier eine eigene Funktion heatColors() definiert. Denn die Funktion heat.colors() hat den Nachteil, dass sie Farben erzeugt, die sehr wenig Kontrast zum weißen Hintergrund besitzen. (Die Dokumentation zu heat.colors() findet sich im Paket The grDevices package unter Palettes.) Mit heatColors(n) werden zunächst die Farben der Palette heat.colors(n+3) erzeugt, dann aber die letzten 3 Komponenten verworfen. Das folgende Skript zeigt die Implementierung:

heatColors <- function(n){
  stopifnot(n > 0)
  if(n < 5){
    return(heat.colors(n))
  } else {
    return( head(heat.colors(n+3), n = n ) )
  }
}

Zeile 22: Mit dem Eingabewert lwd wird die Breite der Höhenlinien festgelegt; sie wird hier gleich 3 gesetzt, um sie besser von den Gitterpunkten abzuheben.

Zeile 23: Das Argument add = TRUE sorgt dafür, dass kein eigener Plot erzeugt wird, sondern die Höhenlinien dem vorangegangenen Plot hinzugefügt werden.

Aufgabe: Setzen Sie in dem Skript, mit dem Abbildung 3 erzeugt wurde N <- 10 und geben Sie die Matrix outer(X = v, Y = v, FUN = "S") aus.

Die Eingabewerte levels und nlevels

Oben wurde in den Beispielen ausdrücklich festgelegt, welche Höhenlinien gezeichnet werden sollen (Eingabewert levels). Kennt man die Eigenschaften der Funktion noch nicht, so kann es schwer fallen geeignete Funktionswerte für levels festzulegen.

Einen Ausweg bietet hier der Eingabewert nlevels. Allerdings ist seine Verwendung nicht so einfach wie man wohl vermutet. Denn die Anzahl nlevels wird nicht wörtlich übernommen, sondern an den Verlauf der Funktion angepasst.

Das folgende Beispiel soll dies demonstrieren; wie nlevels "an den Verlauf der Funktion angepasst" wird, wird im Anschluss an das Beispiel erklärt.

Verwendet wird die Funktion

f(x, y) = x² + y², mit x ∈ [-2; 2], y∈ [-2; 2],

die in Abbildung 4 gezeigt wird. Mit dem gewählten Definitionsbereich lautet der Wertbereich

W_f = [0; 8].

Abbildung 4: Darstellung der Funktion f(x, y) mit dem Rechteck [-2; 2]×[-2; 2] als Definitionsbereich.

Im folgenden Skript werden die Höhenlinien gezeichnet, wobei nlevels = 10 gesetzt wird. Mit plot() wird lediglich das Koordinatensystem gezeichnet; alle relevanten Anweisungen für die Höhenlinien sind in contour() enthalten:

N <- 50
v <- seq(-2, 2, length.out = N)
w <- v*v

plot(x = c(-2, 2), y = c(-2, 2), 
     col = "white", 
     lwd = 1,
     xlab = "x", ylab = "y",
     main = "Höhenlinien von f(x, y) = x^2 + y^2",
     frame.plot = TRUE,
     asp = 1)
grid()
# Koord.achsen:
abline(h = 0, col = "darkgray")
abline(v = 0, col = "darkgray")

contour(x = v, y = v,    # Höhenlinien: v verwenden
        z = outer(w, w, FUN = "+"),
        nlevels = 10,
        col = heatColors(10),
        lwd = 2,
        add = TRUE)

Abbildung 5 zeigt das – unerwartete – Ergebnis.

Abbildung 5: Darstellung der Höhenlinien der Funktion f(x, y) über dem Definitionsbereich.

Die Höhenlinien sind Kreise wie man es bei der rotationssymmetrischen Funktion f(x, y) erwartet. Aber offensichtlich wird die Eingabe nlevels = 10 ignoriert, da die Höhenlinien zu den Funktionswerten 1, 2, ..., 7 dargestellt werden.

Erklären kann man dieses Verhalten mit den default-Werten nlevels = 10, levels = pretty(zlim, nlevels) :

• Der default-Wert für nlevels ist zwar gleich 10, aber um die Funktionswerte der tatsächlich dargestellten Höhenlinien zu berechnen, wird levels = pretty(zlim, nlevels) verwendet.
• Die Funktion pretty() sorgt dafür, dass die Funktionswerte so aus dem Wertbereich ausgewählt werden, dass sie ganzzahlig, halbzahlig und so weiter sind.

Im folgenden Skript werden für nlevels die Werte c(5, 10, 15, 30) gesetzt; Abbildung 6 zeigt das Ergebnis.

N <- 50
v <- seq(-2, 2, length.out = N)
w <- v*v

par(mfrow = c(2, 2))

for(nlvls in c(5, 10, 15, 30)){
  plot(x = c(-2, 2), y = c(-2, 2),
       col = "white", 
       lwd = 1,
       xlab = "x", ylab = "y",
       main = "Höhenlinien von f(x, y) = x^2 + y^2",
       frame.plot = TRUE,
       asp = 1)
  grid()
  # Koord.achsen:
  abline(h = 0, col = "darkgray")
  abline(v = 0, col = "darkgray")
  
  contour(x = v, y = v,
          z = outer(w, w, FUN = "+"),
          # levels = (-5:1),
          nlevels = nlvls,
          col = heatColors(nlvls),
          lwd = 2,
          add = TRUE)
}

Abbildung 6: Darstellung der Höhenlinien der Funktion f(x, y) über dem Definitionsbereich zu den Werten nlevels = 5, 10, 15, 30.

Aufgabe:

Informieren Sie sich über die Funktion pretty() und erklären Sie, welche Höhenlinien in Abbildung 6 ausgewählt wurden.

Die Eingabewerte labcex, drawlabel und method

In den bisherigen Abbildungen war zu erkennen, dass die Höhenlinien mit den zugehörigen Funktionswerten beschriftet wurden. Aber keiner der verwendeten Eingabewerte scheint für diese Beschriftung verantwortlich zu sein. Offensichtlich haben die default-Werte der entsprechenden Eingabewerte für die Beschriftung gesorgt.

Verantwortlich sind die drei Eingabewerte, die in der folgenden Tabelle kurz beschrieben werden:

Bei den bisherigen Beispielen wurde somit die Schriftgröße labcex = 0.6 sowie die Methode method = "flattest" verwendet. Die drei Methoden beschreibt die folgende Tabelle:

Bei der zweiten und dritten Methode kann es vorkommen, dass einige Höhenlinien nicht beschriftet werden.

Für das folgende Beispiel wird die Funktion verwendet:

g(x, y) = x² - y² mit x ∈ [-2; 2], y ∈ [-2; 2].

Ihre Darstellung mit persp() zeigt Abbildung 7. Man erkennt, dass durch den gewählten Definitionsbereich der Wertbereich das Intervall [-4; 4] ist.

Abbildung 7: Darstellung der Funktion g(x, y) mit Hilfe von persp().

In Abbildung 8 werden 7 Höhenlinien (-3, -2, ..., 3) von g(x, y) dargestellt, wobei die drei Methoden verwendet werden (siehe Titel des jeweiligen Diagramms). Die Schriftgröße wurde jetzt mit labcex = 1.0 festgelegt.

Abbildung 8: Darstellung der Höhenlinien der Funktion g(x, y) zu den Funktionswerten -3, -2,..., 3 mit den drei verschiedenen Werten des Argumentes method.

Das Skript, mit dem Abbildung 8 erzeugt wurde, lautet:

ms <- c("simple", "edge", "flattest")
N <- 201
v <- seq(-2, 2, length.out = N)
w <- v * v

par(mfrow = c(3, 1))
par(bg = "#F8FCFF")

for(i in (1:3)){
  plot(x = c(-2, 2), y = c(-2, 2), 
       col = "white", 
       lwd = 1,
       xlab = "x", ylab = "y",
       main = paste0("Höhenlinien von x^2 - y^2, method = ", ms[i]),
       frame.plot = TRUE)
  grid()
  abline(h = 0, col = "darkgray")
  abline(v = 0, col = "darkgray")

  contour(x = v, y = v,
          z = outer(w, w, FUN = "-"),
          levels = (-3:3),
          col = heatColors(7),
          lwd = 2,
          labcex = 1.0,
          method = ms[i],
          add = TRUE)
}