Kapitel

3

Funktionen Die Darstellung von Funktionen und Daten ist eine wichtige Anwendung von MatLab. Sie wird in diesem Kapitel einen breiten Rahmen einnehmen, da die zugeh¨origen Fertigkeiten auch bei der Darstellung von Messdaten, z.B. im Praktikum, hilfreich sind. Da MatLab kein algebraisches System ist, ist eine Differentiation oder Integration von Funktionen nur numerisch m¨ oglich. Die zugeh¨ origen Verfahren werden im zweiten Teil dieses Kapitels beschrieben. Sie geben gleichzeitig eine Einf¨ uhrung in die Verwendung der if und do Strukturen in MatLab, wie sie in komplexen Problemen zur Steuerung ben¨otigt werden.

3.1

Darstellung von Funktionen

In diesem einf¨ uhrenden Teil werden wir uns auf Funktionen einer Variablen beschr¨anken, d.h. Funktionen der Form x(t) oder y = f (x). Die Darstellung von Funktionen von zwei Variablen, f (x, y), erfolgt ¨ ahnlich; sie wird in Abschn. 3.2 genauer dargestellt.

3.1.1

Einfache Plots

Die Erzeugung funktionaler Zusammenh¨ange und eine sehr einfache Form der Darstellung haben wir bereits im Zusammenhang mit Abb. 1.2 abgesprochen. Wir k¨onnen daher kommentarlos mit der Sequenz >> x=[0:0.01:5]; ←Testbild 1 1

1

berechnete Werte, VORLÄUFIG! 0.8

0.8

0.6

y=exp(x)*sin(10*x)

0.6 0.4 0.2 0 −0.2

0.4

0.2

0

−0.2

−0.4

−0.4 −0.6

−0.6 −0.8

−0.8 0

0.5

1

1.5

2

2.5

3

3.5

4

4.5

0

0.5

1

1.5

2

2.5

3

3.5

4

x−Achse

5

Abbildung 3.1: Einfacher (links) und beschrifteter (rechts) Funktionsplot 22

4.5

5

3.1. DARSTELLUNG VON FUNKTIONEN

23

Abbildung 3.2: Links: Axes-Property Editor zur interaktiven Erstellung einer Bildbeschriftung; Rechts: einfaches Pull-Down Men¨ u zur Ver¨anderung der Abbildung >> y=exp(-x).*cos(10.*x); ←>> plot(x,y) ←den im linken Teil von Abb. 3.1 gezeigten Plot der Funktion y = e−x cos(10x) erzeugen. Dieser Abbildung fehlt allerdings ein wichtiger Bestandteil: die Achsenbeschriftung. Diese kann interaktiv erg¨ anzt werden, auch k¨onnen Beschriftungsgr¨oßen, Abst¨ande der Ticks usw. ver¨andert werden. Dazu kann im Men¨ u des Figure-Fensters unter Edit → Axes Properties ein Men¨ u ge¨ offnet werden, entsprechend dem im linken Teil von Abb. 3.2 gezeigten. In diesem Men¨ u lassen sich die Achsenbeschriftungen f¨ ur alle Achsen sowie die Eigenschaften der Achsen, inklusive des dargestellten Wertebereichs, recht komfortabel einstellen. Eine noch einfachere Alternative ist im rechten Teil von Abb. 3.2 angedeutet: unter dem Men¨ upunkt insert erscheint eine Auswahl von einzuf¨ ugenden Objekten. W¨ahlt man X Label aus, so erscheint unterhalb der x-Achse ein kleines weißes K¨astchen, in das der Text eingetragen werden kann. Klickt man dieses K¨astchen sp¨ater mit der linken Maustaste an, so erscheint ein weiteres Men¨ u, mit dessen Hilfe die Eigenschaften der Achsenbeschriftung ver¨andert werden k¨ onnen. ¨ Uber die entsprechenden Punkte in diesen Men¨ us lassen sich auch weitere Eigenschaften der Abbildung ver¨ andern, wie z.B. Farben der Kurven, Symbole, Strichst¨arke usw. Das Speichern der Abbildung (inklusive der daran vor genommenen Manipulationen) erfolgt u u. Dort kann auch das Datenformat ausw¨ahlt ¨ber den Unterpunkte Export im File-Men¨ werden, in das die Abbildung exportiert werden soll, z.B. pdf, jpeg, eps oder bmp.

3.1.2

Fine-Tuning

Das interaktive Sch¨ onen der Abbildungen erscheint anfangs zwar angenehm, da die einzelnen Schritte leicht zu kontrollieren sind, erfordert aber umgekehrt viele Einzelschritte und wird daher leicht un¨ ubersichtlich und fehleranf¨allig, insbesondere, wenn wir f¨ ur ein Protokoll oder eine Arbeit viele Abbildungen in einheitlichem Stil erzeugen wollen. F¨ ur eine einfache Beschriftung k¨onnen wir die folgende Befehlssequenz in ein MatLabSkript einbauen >> x=[0:0.01:5]; ←c M.-B. Kallenrode

12. Januar 2005

24

KAPITEL 3. FUNKTIONEN

-: -. none

durchgezogene Linie (default) gestrichelte Linie gepunktete Linie strichpunktierte Linie keine Linie, nur die Marker an den Punkten

Tabelle 3.1: LineStyles beim Plotten in MatLab >> >> >> >> >> >> ←>> plot beschrift

y=exp(-x).*cos(10.*x); ←testbild = plot(x,y,’k’); ←xlabel(’x-Achse’,’Fontsize’,16); ←ylabel(’y=exp(x)*sin(10*x)’,’Fontsize’,16); ←title(’Testbild 1’,’Fontsize’,20); ←text(0.5,0.92,’berechnete Werte, VORL¨ AUFIG!’,’Color’,’red’,’FontWeight’,’bold’,’Fontsize’,14); saveas(testbild,’Testbild1.pdf’) ←-

und erhalten als Ergebnis den Plot im rechten Teil von Abb. 3.1. Dieses Beispiel findet sich im m-File plot beschrift. Kurvenparameter: Farben und Symbole

plot

Die ersten beiden Zeilen des obigen MatLab-Fragments sind bekannt, die dritte Zeile enth¨alt den normalen Plotbefehl, zus¨ atzlich wird zwischen den Hochkommata ein k als Parameter u bergeben. Dieser Parameter weist MatLab an, die Kurve nicht in blau sondern in schwarz ¨ zu plotten. Beim Aufruf von plot k¨ onnen verschiedene Parameter u ¨bergeben werden, die die Darstellung der Kurve spezifizieren. Dazu geh¨oren auch Strichst¨arke und Linienart. Der allgemeine Aufruf von plot beinhaltet daher auch die Gr¨oße PropertyName, die zusammen mit einem sie spezifizierenden Wert angegeben werden muss: plot(...,’PropertyName’,PropertyValue,...) Die wichtigsten Eigenschaften sind:

• LineWidth wird als Skalar u ¨bergeben, der die Linienst¨arke in pt angibt; Standardeinstellung ist 0.5pt. LineStyle • LineStyle spezifiziert die Art der Linie, als Eingabe werden die in Tabelle 3.1 gegebenen Zeichen, jeweils in Hochkommata eingeschlossen, akzeptiert: • zur Markierung der einzelnen Datenpunkte in einem Plot k¨onnen Marker verwendet werden. Diese Marker k¨ onnen einfach mit der weiter unten gegebene Kurzform Linespec Linespec definiert werden, die in einem String die drei Parameter Linientyp, Marker und Farbe u ¨bergibt. Bei der Markierung kann u ¨ber den PropertyName MarkerSize die Gr¨oße des MarkerSize Markers eingestellt werden. Wie bei der Linienweite wird sie in Punkten pt angegeben, die Voreinstellung entspricht einer Gr¨oße von 6 pt. Außerdem k¨onnen die Marker farbig MarkerEdgeColor dargestellt werden, wobei u ur die Um¨ber den PropertyName MarkerEdgeColor die Farbe f¨ randung des Markers gew¨ ahlt werden kann, u ¨ber den PropertyName MarkerFaceColor die MarkerFaceColor Farbe innerhalb des Markers. F¨ ur die Marker sind verschiedene Symbole vereinbart, die jeweils durch ein einzelnes Zeichen spezifiziert werden. Die Symbole f¨ ur die Marker k¨onnen u ¨ber den Shortcut Linespec Linespec u urzungen sind in Tabelle 3.2 zusammen gefasst. ¨bergeben werden, die Abk¨ • Farben werden u urzungen ¨ber ihre Darstellung im RGB-Code gew¨ahlt oder u ¨ber die Abk¨ in der untersten Zeile von Tabelle 3.3. LineWidth

Eine verk¨ uzte Form gibt nicht PropertyName und den dazugeh¨origen Wert an sondern enth¨alt alle Zusatzinformationen in einer Zeichenkette: >> plot(x,y,’-.or’) ←12. Januar 2005

c M.-B. Kallenrode

3.1. DARSTELLUNG VON FUNKTIONEN

+ o * . x s d ∧

v > < p h

25

Plus-Zeichen Kreis Sternchen (Asterisk) Punkte Kreuzchen Quadrate (square) Rhomben (diamond) aufw¨arts gerichtetes Dreieck abw¨arts gerichtetes Dreieck nach rechts gerichtetes Dreieck nach links gerichtetes Dreieck Pentagram (5zackiger Stern) Hexagram (6zackiger Stern)

Tabelle 3.2: Marker und die zu ihrer Eingabe verwendeten Abk¨ urzungen FontAngle FontName Fontsize FontWeight Color

z.B. kursiv Schrifttype Gr¨ oße Schriftsatz in Punkten z.B. Fettdruck Textfarbe

normal, italic, oblique Helvetica 14, 16, ... light, normal, demi, bold y, m, c, r, g, b, w, k

Tabelle 3.3: Parameter zur Definition der Achsenbeschriftung und zul¨assige Werte. Bei der Farbe k¨ onnen statt der Abk¨ urzungen auch die langen Farbnamen verwendet werden: yellow, magenta, cyan, red, green, blue, white, black plottet die y-Werte gegen die x-Werte, wobei die einzelnen Werte mit einer strichpunktierten Linie (-.) verbunden sind, durch einen Kreis (o) markiert sind und wobei Linie und Marker in rot (r) dargestellt werden. Anschließend k¨onnen weitere Spezifikationen (z.B. Markergr¨oße) u ¨ber einen PropertyName mit seinem Wert angegeben werden. Auf diese Weise kann auch die Farbe der Marker gesetzt werden – dann haben Linie und Marker gegebenenfalls unterschiedliche Farben. Im Beispiel wurde der Befehl plot(x,y) nicht einfach aufgerufen sondern es wurde ihm durch phandle = plot(x,y) ein Handle zugewiesen, mit dem die Abbildung eindeutig identifiziert werden kann. Dieser Handle ist z.B. notwendig beim Speichern der Abbildung (letzte Zeile im Beispiel). Er kann auch verwendet werden, um Eigenschaften der Abbildung wie Kurven- oder Achsenparameter nachtr¨aglich zu ver¨andern. Achsenparameter: Griechisch-r¨ omischer Freistil Die folgenden drei Zeilen im zum rechten Teil von Abb. 3.1 geh¨origen Fragment legen die Achsen- und Bildbeschriftung fest. Die allgemeine Syntax f¨ ur den Befehl xlabel ist xlabel(...,’Textstring’,’PropertName’,PropertyValue,...) wobei Textstring den in Hochkomma gesetzten Text der Achsenbeschriftung enth¨alt und PropertyName den Namen f¨ ur eine bestimmte Eigenschaft der Achsenbeschriftung, im Beispiel f¨ ur die Schriftgr¨ oße. MatLab akzeptiert f¨ ur die Eigenschaften u.a. die in Tabelle 3.3 gegebenen PropertyNames und PropertyValues, weitere Details k¨onnen in der Hilfe unter Text ←- Property List gefunden werden. Die einfachste Variante von xlabel enth¨alt nur den in Hochkommata gesetzten Text der Achsenbeschriftung: xlabel(’x-Achse’). Die Achsenbeschrftung erfolgt in schwarz mit einer Schriftgr¨ oße von 10 pt. Bei der Darstellung physikalischer Gr¨oßen ben¨otigen wir h¨aufig Sonderzeichen, z.B. π, griechische Buchstaben, Verkn¨ upfungen wie ± oder Vergleichsoperationen wie ≤. Diese stellt c M.-B. Kallenrode

12. Januar 2005

xlabel

26

KAPITEL 3. FUNKTIONEN

\alpha \beta \gamma \delta \epsilon \zeta \eta \theta \vartheta \iota \kappa \lambda \mu \nu \xi \pi \rho \sigma \varsigma \tau \equiv \Im \otimes \cap \supset \int \rfloor \lfloor \perp \wedge \rceil \vee \langle

α β γ δ  ζ η θ ϑ ι κ λ µ ν ξ π ρ σ ς τ ≡ = ⊗ ∩ ⊃ R c b ⊥ ∧ e ∨ h

\upsilon \phi \chi \psi \omega \Gamma \Delta \Theta \Lambda \Xi \Pi \Sigma \Upsilon \Phi \Psi \Omega \ forall \exists \ni \cong \approx \Re \oplus \cup \subseteq \in \lceil \cdot \neg \times \surd \varpi \rangle

υ φ χ ψ ω Γ ∆ Θ Λ Ξ Π Σ Υ Φ Ψ Ω ∀ ∃ 3 ∼ = ≈ < ⊕ ∪ ⊆ ∈ d · ¬ × √ $ i

\sim \leq \infty \clubsuit \diamondsuit \heartsuit \spadesuit \leftrightarrow \leftarrow \uparrow \rightarrow \downarrow \circ \pm \geq \propto \partial \bullet \div \neq \aleph \wp \oslash \supseteq \subset \o \nabla \ldots \prime \O \mid \copyright

∼ ≤ ∞ ♣ ♦ ♥ ♠ ↔ ← ↑ → ↓ ◦ ± ≥ ∝ ∂ • ÷ 6 = ℵ ℘  ⊇ ⊂ ø ∇ ... 0 Ø | c

Tabelle 3.4: Eingabe von Sonderzeichen bei der Beschriftung von Abbildungen

MATLAB in der gleichen Form dar, wie sie bei der Eingabe in LATEX verwendet werden. Die Symbole sind in Tabelle 3.4 zusammen gefasst. Text innerhalb der Abbildung text

Auch innerhalb einer Abbildung kann mit Hilfe der Funktion text Text angegeben werden, z.B. um Kurven zu beschriften oder besondere Punkte zu markieren. Der Befehl dazu lautet text(x,y,’String’,’PropertyName’,PropertyValue,...)

saveas

Die Parameter x und y geben die Koordinaten der linken unteren Ecke der Textbox an, der String ’String’ enth¨ alt den Text und u ¨ber PropertyName kann dessen Darstellung ver¨andert werden. Die vorletzte Zeile des Beispiels enth¨alt diesen Befehl und erzeugt damit die warnende Beschriftung innerhalb der Abbildung. Diese Zeile liefert gleichzeitig auch ein Beispiel f¨ ur die Anwendung der in Tabelle 3.3 vorgestellten Eigenschaften. In der letzten Zeile des Beispiels dient der Befehl saveas der Speicherung der Abbildung in einer Datei: der erste Parameter in der Klammer ist der Handle des Bildes, d.h. dient der Identifikation der zu speichernden Abbildung, der zweite Parameter u ¨bergibt als String 12. Januar 2005

c M.-B. Kallenrode

3.1. DARSTELLUNG VON FUNKTIONEN

27

Winkelfunktionen ↓ cos(x)

1 0.8 0.6

y−Achse)

0.4 0.2 ← sin(x)

0 −0.2 −0.4 −0.6 −0.8 −1

−6

−4

−2

0

2

4

6

x−Achse

Abbildung 3.3: Zwei Kurven in einer Abbildung, erzeugt mit dem Skript winkelfunktionen den Namen der Datei, in der die Abbildung gespeichert werden soll.1 Die Dateiendung gibt gleichzeitig das Format, in dem die Abbildung gespeichert wird, in diesem Fall pdf.

3.1.3

Linear oder Logarithmisch?

Erstrecken sich die Funktionswerte u ¨ber mehrere Gr¨oßenordnungen, so ist es sinnvoller, sie logarithmisch darzustellen. Die MatLab-Funktion plot stellt beide Achsen linear dar, al¨ lerdings kann mit Hilfe der Men¨ us nachtr¨aglich interaktiv eine Anderung auf logarithmische Darstellung vorgenommen werden, z.B. durch Verwendung des Pop-Up Men¨ us unten im Axes-Property Editor (vgl. linkes Teilbild in Abb. 3.1). Eleganter ist nat¨ urlich die direkte Wahl einer logarithmischen Darstellung. Dazu wird der Befehl plot durch semilog ersetzt und wir erhalten eine halblogarithmische Darstellung mit linearer x und logarithmischer y-Achse. Mit dem Befehl loglog statt plot l¨asst sich eine Abbildung erzeugen, in der beide Achsen logarithmisch skaliert sind (doppeltlogarithmischee Darstellung). Die Syntax beider Befehle entspricht der von plot.

3.1.4

semilog loglog

Mehrere Kurven in ein Bild

Die Darstellung mehrerer Kurven erfordert die mehrfache Verwendung des Befehls plot. Da dieser jedoch jedes mal das Fenster erst l¨oscht und dann den Plot neu erstellt, lassen sich auf diese Weise nicht mehrere Kurven in einer Abbildung darstellen. Abhilfe schafft der Befehl hold: mit hold on h¨ alt MatLab die vorherige Abbildung fest. Jeder folgende Aufruf von plot zeichnet die Abbildung in den gleichen Rahmen. Auf diese Weise lassen sich mehrere Funktionen in einer Abbildung darstellen. Am Ende sollte allerdings auch ein hold off gegeben werden, da MatLab sonst auch alle weiteren folgenden Abbildungen in diese zeichnen w¨ urde. Abbildung 3.3 zeigt ein Beispiel f¨ ur eine Abbildung mit zwei Kurven, erzeugt mit dem Skript winkelfunktionen: clear x=[-2*pi:0.01:2*pi]; y1=cos(x);y2=sin(x); 1 Die Umkehrung, das Einlesen einer Abbildung aus einer Datei, erfolgt mit dem Befehl imread. Das Einlesen von Bildern in MatLab eher f¨ ur fortgeschrittene Bildbearbeitung von Interesse als f¨ ur unsere mathematischen Anwendungen, daher wird im Rahmen des Skripts nicht weiter darauf eingegangen. Alle n¨ aheren Informationen finden Sie in der MatLab-Hilfe.

c M.-B. Kallenrode

12. Januar 2005

hold on

hold off

winkelfunktionen

28

KAPITEL 3. FUNKTIONEN

20

8

0

6

−20

4

−40

2

−60

0 −10

−8

−6

−4

−2

0

2

4

6

8

y2=−x2*x2+9

y1=x1*x1/10

Plot zweier Funktion in einen Graphen 10

−80 10

x−Achse

Abbildung 3.4: Abbildung mit zwei y-Achsen, erzeugt mit dem MatLab-Skript plot 2yachs phandle = plot(x,y1,’-k’); hold on; plot(x,y2,’-.k’); axis([-2*pi 2*pi -1.2 1.2]); xlabel(’x-Achse’,’Fontsize’,16); ylabel(’y-Achse)’,’Fontsize’,16); title(’Winkelfunktionen’,’Fontsize’,20); text(0-0.1,cos(0)+0.04,’\downarrow cos(x)’); text(0,sin(0),’\leftarrow sin(x)’); saveas(phandle,’winkelfunktionen.pdf’); hold off

axis

Neu in diesem MatLab-Fragment sind die bereits besprochenen Befehle hold on und hold off in den Zeilen 5 und 14, die Verwendung von Funktionen zur Bestimmung der Orte des Textes in den Zeilen 11 und 12 sowie der Befehl axis, der eine Ver¨anderung der Achsen erlaubt – in diesem Fall eine Festlegung des dargestellten Bereichs in der Form axis([xmin xmax ymin ymax]).

3.1.5

plot 2yachs

Plots mit zwei y-Achsen

Stellt man zwei Funktionen in einer Abbildung dar, so haben nicht zwingend beide die gleiche y-Achse: zum einen, auf Grund unterschiedlicher Wertebereiche, zum anderen, da die Funktionen unterschiedliche Dinge darstellen, z.B. Weg und Geschwindigkeit oder Bruttosozialprodukt und Geburtenrate. Das MatLab-Skript plot 2yachs gibt ein Beispiel f¨ ur die Darstellung von zwei Funktionen mit getrennten y-Achsen in einem Bild; das Ergebnis ist in Abb. 3.4 gezeigt. Ohne Kommentierung ist die wesentliche Befehlssequenz des Skriptes x=[-10:0.02:10]; x2 = [-7:0.01:9]; ←y1 = x1.*x1/10; y2 = -x2.*x2 +9; ←[ax,h1,h2] = plotyy(x1,y1,x2,y2); ←title(’Plot zweier Funktion in einen Graphen’,’Fontsize’,20); ←xlabel(’x-Achse’,’Fontsize’,16); ←set(get(ax(1),’Ylabel’),’String’,’y1=x1*x1/10’,’Fontsize’,16); ←set(get(ax(2),’Ylabel’),’String’,’y2=-x2*x2+9’,’Fontsize’,16); ←set(h1,’LineStyle’,’--’); ←12. Januar 2005

c M.-B. Kallenrode

3.1. DARSTELLUNG VON FUNKTIONEN

29

set(h2,’LineStyle’,’-.’); ←saveas(h1,’plot2yachs.pdf’) ←In der ersten Zeile wird der x-Wertebereich f¨ ur die beiden Funktionen fest gelegt – trotz gemeinsamer x-Achse haben sie unterschiedliche Wertebereiche. Dies kann bei Messdaten, insbesondere Zeitreihen, h¨ aufiger vorkommen: Klimatologen vergleichen z.B. Temperaturen und den Kohlendioxid-Gehalt der Atmosph¨are, wobei erstere bereits seit ca. 1870 zuverl¨assig und kontinuierlich gemessen werden, letztere jedoch erst seit ca. 1950. In der zweiten Zeile werden die beiden Funktionen berechnet, in der dritten geplottet, wobei der Handle etwas aufwendig ist, da nicht nur f¨ ur die Abbildung sondern auch f¨ ur jede der y-Achsen ein Handle u ¨bergeben werden muss, damit diese getrennt ver¨andert werden ¨ k¨onnen. Diese Anderungen erfolgen jeweils mit den set-Befehlen: der erste bestimmt die Beschriftung der ersten y-Achse, der zweite die der zweiten. Die beiden anderen set Befehle bestimmen die Darstellung der Funktionsgraphen. Weitere Informationen finden sich in der MatLab-Hilfe unter plotyy.

3.1.6

set

plotty

Funktionen in Parameterdarstellung

Funktionen in Parameterdarstellung k¨onnen in MatLab einfach dargestellt werden. Dazu wird als erstes der Wertebereich des Parameters definiert und ein Vektor mit Werten dieses Parameters erzeugt. Zu diesem Vektor werden jeweils die x- und y-Werte bestimmt und anschließend gegeneinander geplottet. Beispiel 5 Eine Archimedische Spirale entsteht, wenn sich ein K¨orper mit konstanter Geschwindigkeit v entlang eines Strahls nach außen bewegt und dieser Strahl gleichzeitig mit konstanter Winkelgeschwindigkeit ω um den Ursprung rotiert. F¨ ur den radialen Abstand r(t) und den Drehwinkel ϕ(t) erhalten wir r(t) = v t

und

archimedespolar

φ(t) = ω t .

F¨ ur eine graphische Darstellung k¨ onnen wir zwei Formen w¨ ahlen, Polarkoordinaten oder kartesische Koordinaten. F¨ ur die Darstellung in Polarkoordinaten ben¨ otigen wir außer der oben gegebenen Parameterdarstellung nur den Befehl polar f¨ ur die graphische Darstellung. Die zugeh¨ orige Befehlssequenz ist

polar

archimedeskartesisch

>> v=1; ohm = 1.; t=[0:0.01:20]; ←>> r = v*t; phi = ohm*t; ←>> phandle=polar(phi,r,’k’); ←In der ersten Zeile werden die beiden Parameter v und ohm, d.h. die Geschwindigkeit v und die Winkelgeschwindigkeit ω vorgegeben. In der zweiten Zeile werden die zugeh¨ origen Abst¨ ande und Drehwinkel berechnet und in der dritten geplottet. Das Ergebnis ist im linken Teil von Abb. 3.5 gezeigt. F¨ ur die Darstellung in kartesischen Koordinaten m¨ ussen wir die Polarkoordinaten erst mit Hilfe von pol2cart umwandeln und dann als kartesische Koordinaten konventionell mit plot darstellen:

pol2cart

>> [x,y] = pol2cart(phi,r); ←>> plot(x,y,’k’); ←>> axis equal; ←Das Ergebnis ist im rechten Teil von Abb. 3.5 dargestellt. Der Befehl axis equal bewirkt, dass beide Achsen gleiche Skalierung haben – sonst wird die Spirale verzerrt dargestellt. 2

3.1.7

Im Vorgriff aus Praktikum: Darstellung von Messdaten

Die Darstellung von Messdaten folgt den gleichen Regeln wie die Darstellung von Funktionen. Die wesentlichen Unterschiede d¨ urften darin liegen, dass bei den Messdaten nicht unbedingt eine Linie dargestellt wird sondern die Betonung auf der Darstellung der einzelnen Messpunkte durch Marker liegt und dass die Messdaten m¨oglicherweise in einer externen Datei liegen, aus der sie erst eingelesen werden m¨ ussen. c M.-B. Kallenrode

12. Januar 2005

axis equal

30

KAPITEL 3. FUNKTIONEN

90

20

120

60

15

15 10

150

30

10

5

5

180

0

0

−5

330

210

−10

−15

300

240 270

−20

−15

−10

−5

0

5

10

15

20

Abbildung 3.5: Archimedische Spirale als Beispiel f¨ ur eine Funktion in Parameterdarstellung. Links Darstellung in Polarkoordinaten (Skript archimedespolar), rechts Darstellung in kartesischen Koordinaten (Skript archimedeskartesisch) 10

120

100

9.5

Ort

Temperatur

80

9

60

40

20

8.5

0

2

4

6

8

10

12

14

16

18

20

0

1

Messpunkt #

2

3

4

5

6

7

8

9

10

Zeit

Abbildung 3.6: Zwei Beispiele f¨ ur die Darstellung von Messwerten

tempplot

Messdaten liegen in der Regel als eine Matrix vor. Im einfachsten Fall handelt es sich dabei um einen Spaltenvektor, der eine Reihe von Zahlen enth¨alt. Diese Art von Messdaten entsteht u.a. bei der Aufnahme einer Zeitserie, d.h. wenn alle 5 min die Temperatur gemessen und weg geschrieben wird. Als Beispiel dient der Datenfile temp.dat. Seine Darstellung im linken Teil von Abb. 3.6 erfolgte mit Hilfe des Skripts tempplot: T=load(’temp.dat’); tplot=plot(T,’--sr’); xlabel(’Messpunkt #’,’FontSize’,14); ylabel(’Temperatur’,’FontSize’,14); saveas(tplot,’tplot.pdf’)

plot

speedplot

Die wesentliche Neuerung ist die Verwendung von plot in der zweiten Zeile: hier wird nur ein Vektor als Parameter u ¨bergeben; die Elemente des Vektors werden gegen ihre Ordnungsnummer aufgetragen. Daten k¨ onnen auch als mehrspaltige Matrix vorliegen, z.B. bei der Beobachtung einer Bewegung in einer Datei, die in der ersten Spalte die Zeitpunkte enth¨alt, zu denen gemessen wurde, und in den folgenden Spalten jeweils die Messungen des Ortes zu diesen Zeiten f¨ ur unterschiedliche Messreihen. Ein Beispiel ist in der Datei speed.dat gegeben; die im rechten Teil von Abb. 3.6 gezeigte Abbildung wurde daraus mit Hilfe des Skripts speedplot erzeugt: v=load(’speed.dat’); tplot=plot(v(:,1),v(:,2),’--sr’); hold on; 12. Januar 2005

c M.-B. Kallenrode

3.1. DARSTELLUNG VON FUNKTIONEN

31

plot(v(:,1),v(:,3),’--og’); plot(v(:,1),v(:,4),’-->b’); for i=1:length(v) vm(i) = sum(v(i,2:4))/3; end plot(v(:,1),vm,’k’,’LineWidth’,2); xlabel(’Zeit’,’FontSize’,14); ylabel(’Ort’,’FontSize’,14); hold off saveas(tplot,’splot.pdf’) In der ersten Zeile wird die Tabelle der Messwerte aus der Datei speed.dat in die Matrix v eingelesen. In den folgenden Zeilen werden die Spalten 2 bis 4 der Matrix jeweils gegen die erste Spalte geplottet. Die for-Schleife mittelt u ¨ber die Spalten 2 bis 4 der Matrix, d.h. es wird zu jedem Zeitpunkt der Mittelwert der Messungen des Ortes gebildet. Dieses Information wird anschließend als dicke schwarze Linie in den Plot eingetragen. Der Rest ist Achsenbeschriftung und Speichern der Abbildung. Gerade bei der Darstellung von Messwerten kann es zur Orientierung hilfreich sein, ein Gitternetz aus waagerechten und senkrechten Linien in der Zeichnung zu haben. Dieses l¨asst sich, nach dem Aufruf der Funktion plot, mit Hilfe von grid on ein- und mit Hilfe des Befehls grid off wieder ausschalten.

3.1.8

grid on grid off

Erg¨ anzung: fplot statt plot

Zur Darstellung eines gegebenen funktionalen Zusammenhangs kann statt der MatLabFunktion plot die Funktion fplot verwendet werden. W¨ahrend bei plot die x- und y-Werte als Vektoren u ¨bergegeben werden, akzeptiert fplot als Eingabe die Funktion als String: >> fplot(’x + sin(x).*cos(x)’,[-2*pi 2*pi],’k’) ←Die weiteren Modifikationen an der Abbildung erfolgen wie bei plot. Die darzustellende Funktion wird an fplot entweder, wie im Beispiel, als in Hochkommata gesetzter String u ¨bergeben, als Name eines m-Files, in dem die Funktion definiert ist, oder als Handle einer an anderer Stelle im m-File definierten Funktion. Eine etwas l¨ angere Formulierung des obigen Beispiels weist erst die Funktion und ihre Grenzen jeweils einer Variablen zu und ruft anschließend fplot auf: >> fun=’x + sin(x).*cos(x)’; lims=[-2*pi 2*pi]; fplot(fun,lims,’k’) ←Eine andere M¨ oglichkeit, die Funktion an fplot zu u ¨bergeben, ist die Verwendung eines Handles: >> fhandle = @(x) x+sin(x).*cos(x); fplot(fhandle,[-2*pi 2*pi],’k’) ←Der Parameter hinter dem Klammeraffen @ gibt an, welches (bzw. bei einer Liste welche) Argumente an die im anschließenden Ausdruck dargestellte Funktion u ¨bergeben werden sollen. In diesem Fall ist es ein Parameter. Die Funktion wird bei der Definition des Handles nicht als String angegeben sondern ohne Hochkommata, so wie wir sie auch direkt in einem m-File eingeben w¨ urden. Genauere Informationen finden sich in der MatLab-Hilfe unter Anonymous Functions. In der oben gegebenen Form erzeugt fplot keine Vektoren der x und y-Werte – im Workspace treten keine zus¨ atzlichen Variablen auf. fplot kann jedoch auch verwendet werden, um anstelle der Abbildung zwei Vektoren X und Y mit den Ordinaten- und Abszissenwerten zu erstellen: >> [x y] = fplot(’x.*x + 2*x + 4’,[-5 5]); ←Diese k¨ onnen dann mit plot dargestellt werden. c M.-B. Kallenrode

12. Januar 2005

fplot

32

KAPITEL 3. FUNKTIONEN

eval und feval Die Funktion fplot greift auf die MatLab-Funktion eval zur¨ uck, die eine als String oder u ¨ber einen Handle gegebene Funktion an verschiedenen Stellen bestimmt: >> fun=’x + sin(x).*cos(x)’; x=2; eval(fun) ←y = 1.6216 liefert den Funktionswert von y = x + sin x cos x an der Stelle x = 2; >> fun=’x + sin(x).*cos(x)’; x=[1 2 3 4]; eval(fun) ←y = 1.4546 1.6216 2.8603 4.4947 feval

liefert die Funktionswerte an den im Vektor x spezifizierten Punkten.2 Eng mit eval verwandt ist die universeller verwendbare Funktion feval: >> fhandle = @(x) x+sin(x).*cos(x); x=2; y=feval(fhandle,x) ←y = 1.6216 Der Vorteil von feval liegt darin, dass feval ein Handle als Argument akzeptiert. Die M¨ oglichkeit, eine Funktion als String, u ¨ber einen Handle oder durch Aufruf einer Datei auszuwerten, hat immer dann Vorteile, wenn die Funktion nicht direkt in den m-File geschrieben werden kann (z.B. bei einem GUI oder einem Algorithmus, der mehrfach auf verschiedene Funktionen angewandt werden soll) oder wenn die Funktion innerhalb eines Programms mehrfach ausgewertet werden muss. Letzteres ist z.B. bei vielen numerischen Verfahren der Fall.

3.2

Funktionen mehrerer Variablen

Funktionen von mehreren Variablen lassen sich auf Grund der Beschr¨anktheit unserer Wahrnehmung nur in geringem Maße darstellen: f¨ ur Funktionen in Abh¨angigkeit von zwei Variablen kann der Funktionsgraph als eine Fl¨ ache im dreidimensionalen Raum dargestellt werden, f¨ ur die parametrische Darstellung einer Funktion im 3D kann diese als Bahn dargestellt werden.

3.2.1

meshgrid mesh surf

Funktionen als Fl¨ achen im 3D

Die Darstellung von Funktionen als Fl¨achen im 3D erfolgt analog zur der von Funktionen in Abh¨angigkeit von einer Variablen. Als erstes definieren wir ein Gitter in der xy-Ebene mit einem vorgegebenen Bereich und Gitterabstand mit Hilfe des Befehls meshgrid. Dies entspricht der Definition der x-Achse als Vektor bei der Darstellung von Funktionen einer Variablen. Anschließend wird die Funktion berechnet und mit mesh oder surf dargestellt. Als Beispiel verwenden wir wieder das Tal f (x, y) = x2 + y f¨ ur den Bereich x ∈ [−2, 2] und y ∈ [0, 4]. Als erstes definieren wir das Grundgitter: >> [x,y] = meshgrid([-2:0.2:2],[0:0.1:4]); ←-

meshgrid

Die beiden Parameter von meshgrid sind Vektoren, wie wir sie bei der Definition der x-Achse 2 Achten Sie auf das Multiplikationszeichen! fplot arbeitet auch dann vern¨ unftig, wenn bei der Definition der Funktion ein Multiplikationszeichen * statt der punktweisen Multiplikation .* verwendet wurde. Auch eval liefert beim ersten Beispiel einen vern¨ unftigen Wert zur¨ uck, da die Funktion nur f¨ ur einen Skalar ausgewertet wird. Beim zweiten Beispiel dagegen wird die Funktion f¨ ur einen Vektor ausgewertet, so dass es bei Verwendung von * Probleme bei der Ausf¨ uhrung des Produkts sin(x)*cos(x) gibt und MatLab sich mit der u ¨blichen Fehlermeldung ??? Error using ==> mtimes; Inner matrix dimensions must agree. verabschiedet.

12. Januar 2005

c M.-B. Kallenrode

eval

3.2. FUNKTIONEN MEHRERER VARIABLEN

8

8

6

6

4

4

2

2

0 4

0 4 3

2 1

2

0

1

−1 0

33

3

2 1

2

0

1

−2

−1 0

−2

Abbildung 3.7: mesh und surf zur Darstellung von Funktionen als Fl¨achen im 3D als Vektor verwendet haben. Sind x und y-Bereich identisch, so ist es ausreichend, einen Vektor anzugeben, z.B. meshgrid([-2:0.2:2]). Anschließend berechnen wir die Funktionswerte >> z=x.∧ 2+y; ←und stellen sie mit Hilfe von

mesh

>> mesh(x,y,z) als Gitternetz (linkes Teilbild in Abb. 3.7) oder mit Hilfe von

surf

>> surf(x,y,z) als Oberfl¨ ache (rechtes Teilbild in Abb. 3.7) im 3D dar. Achsenbeschriftung und Kurvenparameter werden wie bei normalen Funktionsplots im 2D manipuliert. Beispiel 6 Mit Hilfe des Skripts sombrero >> >> >> >> >> >> >> >> >> >> >> >>

sombrero

[x,y] = meshgrid([-15:1:15]);←r = sqrt(x.*x + y.*y) + eps; ←z = sin(r)./r;←subplot(2,2,1), mesh(x,y,z); ←title(’Mesh’)←subplot(2,2,2), surf(x,y,z); ←title(’Surf’) ←subplot(2,2,3),surf(x,y,z); ←shading interp ←title(’Surf + interpolated Shading’) ←subplot(2,2,4), contour(z,9); ←title(’Contour’) ←-

wurden die verschiedenen in Abb. 3.8 gezeigten Darstellungsformen f¨ ur eine Funktion von zwei Variablen erzeugt. Neben den bekannten Formen mesh (links oben) und surf (rechts oben) findet sich eine Variante von surf im linken unteren Teil. Hier wurde durch den Befehl shading interp auf eine interpolierende Schattierung umgeschaltet und das normalerweise bei surf schwarz dargestellte Gitter unterdr¨ uckt. Die Darstellung im rechten unteren Teilbild unterscheidet sich von den anderen, da hier die Funktion nicht im 3D sondern u ¨ber Isolinien auf die Ebene projiziert dargestellt wird. Der zugeh¨ orige Befehl ist contour. 2

Verwendet man an Stelle der Befehle mesh und surf die Befehle meshc und surfc, so werden zus¨ atzlich zu den Gittern bzw. Fl¨achen im 3D auf der xy-Ebene die Isolinien dargestellt. Außerdem haben wir in dem Beispiel die Funktion subplot kennen gelernt: sie erm¨oglicht es, mehrere Abbildungen auf einem Blatt darzustellen. Die ersten beiden Parameter in c M.-B. Kallenrode

12. Januar 2005

contour meshc surfc subplot

34

KAPITEL 3. FUNKTIONEN

Mesh

Surf

1

1

0

0

−1 20

−1 20

20

0 −20 −20

20

0

0

−20 −20

Surf + interpolated Shading

0

Contour 30

1 20

0

10

−1 20

20

0 −20 −20

0 10

20

30

Abbildung 3.8: Darstellung von Funktionen im 3D subplot spezifizieren die Zahl der Teilabbildungen pro Spalte bzw. Zeile, die dritte Zahl entspricht der laufenden Nummer des Feldes, in das der folgende plot-Befehl die Abbildung platzieren soll. GUI plottest popup

ezsurf ezsurfc ezmesh ezmesh contour setplot

Zur Darstellung von Funktionen zweier Variablen kann das GUI plottest popup verwendet werden, vgl. Abb. 3.9. Es erlaubt die Eingabe einer Funktion, die Eingabe eines (um Null symmetrischen) Wertebereichs f¨ ur die x- und y-Achse sowie u u die ¨ber ein Pop-Up Men¨ graphische Darstellung als Oberfl¨ ache (ezsurf), als Oberfl¨ache mit zugeh¨origer Darstellung der H¨ohenlinien in der xy-Ebene (ezsurfc), als Gitter (ezmesh), gegebenenfalls ebenfalls mit Konturen (ezmeshc), oder als einfacher Konturplot (H¨ohenlinien) in der xy-Ebene (contour). plottest popup ruft die Funktion setplot auf, die sich im gleichen Verzeichnis befinden muss. In diesem GUI werden die Funktionen ezmesh oder ezmeshc statt mesh oder meshc verwendet. Bei mesh und ¨ ahnlichen Plotfunktionen haben wir die Funktionswerte erst explizit ausgerechnet. Das bedeutet, das im MatLab-Skript an einer Stelle diese Funktion wirklich berechnet werden muss.3 Beim GUI wollen wir die Funktion jedoch nicht als Zeile eines MatLab-Skripts eingeben sondern als Zeichenkette, die auch f¨ ur nicht-MatLab-User verst¨andlich ist. Genau hier kommen die Plotroutinen mit dem Pr¨afix ez ins Spiel. Ihre Syntax akzeptiert die Eingabe einer Funktion (statt der bereits berechneten Funktionswerte) und plottet diese routinem¨ aßig im Bereich −2π ≤ x ≤ 2π und −2π ≤ y ≤ 2π. Werden die 3 Dies entspricht dem Unterschied zwischen den Funktionen plot und fplot bei der Darstellung von Funktionen.

12. Januar 2005

c M.-B. Kallenrode

3.2. FUNKTIONEN MEHRERER VARIABLEN

35

Abbildung 3.9: Screenshot des GUI plottest popup y/(1 + x2 + y2)

0.5

0

−0.5 5 5 0

0 −5

y

−5

x

Abbildung 3.10: Funktionsplot mit ezmeshc als Einzeiler: ezmeshc(’y/(1 + x∧ 2 + y∧ 2)’,[-5,5,-2*pi,2*pi]) x und y Bereiche explizit angegeben, wie im GUI m¨oglich, so stellt ezmesh die Funktion in diesem Bereich dar. Mit Hilfe der ez-Routinen lassen sich Funktionsdarstellungen als Einzeiler erzeugen. So liefert die Zeile >> ezmeshc(’y/(1 + x∧ 2 + y∧ 2)’,[-5,5,-2*pi,2*pi]) ←den in Abb. 3.10 gezeigten 3D-Plot mit den zugeh¨origen Konturen. Ver¨anderungen von Achsenbeschriftung und -parametern wieder wie bei plot, weitere Informationen in der MatLabHilfe unter ezmeshc.

3.2.2

Funktionen als Linien im 3D

Beispiel 7 Gyriert ein Elektron um eine Magnetfeldlinie, so entsteht eine Bahn im dreidimensionalen Raum. In Parameterform k¨ onnen wir diese Helix in Polarkoordinaten beschreiben als %(t) = % = const , c M.-B. Kallenrode

ϕ(t) = ω t

und

z = vz t . 12. Januar 2005

36

KAPITEL 3. FUNKTIONEN

4

3

2

1

0 1

2 1

0 −1

0 −1

Abbildung 3.11: Bahn eines geladenen Teilchens um eine Magnetfeldlinie, dargestellt mit plot3 Darin ist % der Abstand des Elektrons vom Mittelpunkt seines Gyratioorbits (und damit von der Feldlinie), ω die Winkelgeschwindigkeit der Gyration um die Feldlinie und vz die konstante Geschwindigkeit der Translation entlang der Feldlinie. Zur Darstellung der sich ergebenden Bewegung definieren wir zuerst die Konstanten und einen Vektor, der den Parameter Zeit enth¨ alt: >> ohm=2.; vz = 0.4; r = 2.; t=[0:0.02:10]; ←Damit berechnen wir zu jedem Zeitpunkt den Aufenthaltsort des Elektrons in Polarkoordinaten mit >> phi=ohm*t; z = vz*t; ←und wandeln anschließend in kartesische Koordinaten um >> [X,Y,Z] = pol2cart(phi,r,z); ←plot3

Die kartesischen Koordinaten lassen sich dann mit dem Befehl plot3 darstellen >> plot3(X,Y,Z,’k’); axis equal; ←der Befehl axis equal wird wieder ben¨ otigt, um Verzerrungen zu vermeiden. Das Ergebnis ist in Abb. 3.11 gezeigt. 2

3.3

Differentialrechnung

MatLab ist kein symbolisches oder algebraisches sondern ein numerisches System. Daher kann MatLab weder differenzieren noch integrieren in dem Sinne, dass es einen geschlossenen Ausdruck f¨ ur die Ableitung oder das Integral liefert. MatLab bedient sich jedoch verschiedener numerischer Verfahren mit deren Hilfe die Steigung einer Funktion in einem bestimmten Intervall als Wertetabelle bzw. in graphischer Form gegeben werden kann ebenso wie ein bestimmtes Integral.4

3.3.1

Nullstellen

Eines der einfachsten Merkmale zur Charakterisierung einer Funktion sind deren Nullstellen, d.h. die Punkte, an denen der Graph der Funktion die x-Achse schneidet. Analytisch werden diese durch Nullsetzen der Funktion und aufl¨osen nach x bestimmt. Da MatLab die 4 Es ist nicht ganz korrekt, dass MatLab uberhaupt keine algebraischen Rechnungen vornehmen kann. Mit ¨ der Symbolic Math Toolbox als Erweiterung k¨ onnte MatLab korrekt symbolisch rechnen. In der normalen MatLab-Version ist dies nicht allgemein m¨ oglich sondern auf Polynome beschr¨ ankt, vgl. Anhang ??.

12. Januar 2005

c M.-B. Kallenrode

3.3. DIFFERENTIALRECHNUNG

37

dazu erforderlichen Operationen nicht durchf¨ uhren kann, werden in MatLab die Nullstellen numerisch gesucht: in einer Wertetabelle der Funktion wird in der Umgebung eines vorgegebenen Wertes nach einem Vorzeichenwechsel gesucht und dieser zur genaueren Bestimmung der Nullstelle weiter eingeschachtelt. Diese Nulltstellensucher erfolgt mit Hilfe der Funktion fzero. Die Argumente sind die Funktion als String (oder als String in eine Variable geschrieben) sowie ein Wert in dessen Umgebung oder ein Intervall in dem die Nullstelle gesucht werden soll. Der Befehl >> x=fzero(’x*x-1’,0.9) ←x = 1 sucht die Nullstelle der Funktion f (x) = x2 − 1 in der N¨ahe von 0.9 und findet ganz korrekt die +1. Mit dem Befehl >> x=fzero(’x*x-1’,[-1.1 0.9]) ←x = - 1 sucht MATLAB im Intervall von -1.1 bis 0.9 und findet als Nullstelle die -1. Da die Nullstellensuche u ¨ber einen Vorzeichenwechsel der Funktionswerte erfolgt, wird f¨ ur die Funktion x2 keine Nullstelle gefunden: >> x=fzero(’x*x’,0.01) ←Exiting fzero: aborting search for an interval containing a sign change because NaN or Inf function value encountered during search (Function value at -1.553326e+154 is Inf) Check function or try again with a different starting value. x = NaN Ebenso scheitert MatLab bei der Nullstellensuche, wenn bei der uns bereits bekannten Funktion f (x) = x2 − 1 das Intervall so gross gew¨ahlt wird, dass es beide Nullstellen enth¨alt: >> fzero(’x*x-1’,[-1.1 1.1]) ←??? Error using ==> fzero The function values at the interval endpoints must differ in sign. Auch hier kann die Funktion nat¨ urlich vorher definiert und dann aus fzero aufgerufen werden: >> fun=’x*x-1’;lims=[-1.1 0.9];fzero(fun,lims) ←Ein derartiges Verfahren kann z.B. dann sinnvoll sein, wenn man nacheinander mehrere Funktion oder nacheinander mehrere Intervalle bei einer Funktion auf Nullstellen untersuchen m¨ochte und die verschiedenen Werte von fun bzw. lims mit Hilfe einer Schleife nacheinander in fzero f¨ uttert.

3.3.2

Numerisches Differenzieren

Die Ableitung einer Funktion ist definiert als der Quotient aus der Differenz der Funktionswerte an zwei um ∆x aus einander liegenden Stellen dividiert durch ∆x f¨ ur den Grenzwert ∆x → 0: f 0 (x) = lim

∆x→0

f (x) + f (x + ∆x) . ∆x

Numerische Differentiation verwendet den gleichen Quotienten, jedoch ohne den Grenzu ur die praktische Durchf¨ uhrung des numerischen Verfahrens gibt es ¨bergang ∆x → 0. F¨ verschiedene M¨ oglichkeiten: c M.-B. Kallenrode

12. Januar 2005

fzero

38

KAPITEL 3. FUNKTIONEN

• Die Vorw¨ arts-Differenz (forward finite difference) δf+ (x) ist der eigentlichen Definition der Ableitung am ¨ ahnlichsten: f (xo + ∆x) − f (xo ) . ∆x In ihrer Genauigkeit ist diese Differenz von erster Ordnung: die Genauigkeit h¨angt von der zweiten Ableitung der Funktion ab, d.h. von der ersten Ableitung der gesuchten Ableitung. F¨ ur stark ver¨ anderliche Funktionen ist diese aber relativ groß, so dass die Vorw¨artsDifferenz in der N¨ ahe starker Ver¨ anderungen ungenau wird. Die Ordnung der Genauigkeit l¨asst sich durch Taylor-Entwicklung zeigen. • Die R¨ uckw¨ arts-Differenz (backward finite difference) δf− (x) ist auf ¨ahnliche Weise definiert δf+ (xo ) =

f (xo ) − f (xo − ∆x) ∆x und ebenfalls von einer Genauigkeit erster Ordnung. • Die zentrale Differenz (centered finite difference) δf (x) ist definiert als δf− (xo ) =

f (xo + ∆x) − f (xo − ∆x) . 2∆x Sie hat eine Genauigkeit von zweiter Ordnung und ist daher wesentlich robuster als die beiden anderen Differenzen-Verfahren. δf (xo ) =

diff

Ein einfaches Verfahren zur numerischen Differentiation verwendet die MatLab-Funktion diff, die aus einem n-elementigen Vektor einen (n-1)-elementigen Vektor erzeugt, der die Differenzen benachbarter Werte enth¨alt: >> x=[1 2 4 7 11 16 22 29 37]; ←>> dx=diff(x) ←dx = 1 2 3 4 5 6 7 8 Beispiel 8 Die Ableitung der Funktion y = sin(x) ist y 0 = cos(x), die zweite Ableitung ist y 00 = − sin(x). Mit der folgenden Befehlssequenz soll diese Ableitung numerisch im Intervall ug bis og bestimmt und mit der analytischen L¨ osung verglichen werden. Da der Vektor, der die Werte der Ableitung enth¨ alt, um ein Element k¨ urzer ist als der Vektor der Funktionswerte, m¨ ussen wir f¨ ur die graphische Darstellung drei x-Achsen definieren: x f¨ ur die Funktion, xv f¨ ur die Ableitung und xvv f¨ ur die zweite Ableitung. >> >> >> >> >> >> >> >> >> >>

winkelableit

ug=0;og=10;deltx=0.1; ←x=[ug:deltx:og];xv=[ug+deltx/2:deltx:og-deltx/2];xvv=[ug+deltx:deltx:og-deltx];←dx=diff(x);dy=diff(sin(x));abl1=dy./dx; ←d2x=diff(xv);d2y=diff(abl1);abl2=d2y./d2x; ←plot(x,sin(x),’k’);hold on;plot(xv,cos(xv),’r’);plot(xvv,-sin(xvv),’g’);←plot(xv,abl1,’vr’);plot(xvv,abl2,’og’); ←xlabel(’x’,’Fontsize’,16); ylabel(’y’,’Fontsize’,16); ←title(’sin(x) und seine Ableitungen (Symbole)’,’Fontsize’,20); ←text(2.2,0.8,’ \leftarrow sin(x)’); text(3*pi/2,0,’ \leftarrow cos(x)’);←text(5.9,-0.8,’-sin(x) \rightarrow’); hold off ←-

Abbildung 3.12 zeigt die Funktion (schwarz) mit ihren beiden analytisch bestimmten Ableitungen (rote bzw. gr¨ une Kurve) sowie die numerisch mit Hilfe der obigen Befehlssequenz bestimmten Ableitungen (Symbole). Trotz des relativ groben Gitters (anschaulich durch den Abstand benachbarter Symbole), stimmt selbst die zweite numerische Ableitung sehr gut mit den analytischen Werten u ¨berein. Das Skript winkelableit basiert auf obiger Befehlssequenz. 2

Das von uns in MATLAB verwendete Verfahren ¨ahnelt dem einer zentralen Differenz, wobei allerdings die Ableitung nur an jedem zweiten Gitterpunkt bestimmt wurde. In der 12. Januar 2005

c M.-B. Kallenrode

3.4. NUMERISCHE INTEGRATION

39

sin(x) und seine Ableitungen (Symbole) 1 ← sin(x)

0.8 0.6 0.4

y

0.2 ← cos(x)

0 −0.2 −0.4 −0.6

−sin(x) →

−0.8 −1 0

1

2

3

4

5

6

7

8

9

10

x

Abbildung 3.12: Funktion y = sin(x) (schwarz) und analytisch bestimmte Ableitungen y 0 = cos(x) (rote Kurve) bzw. y 00 = − sin(x) (gr¨ une Kurve) sowie zugeh¨orige numerisch bestimmte Ableitungen (Symbole) echten zentralen Differenz w¨ urden wir Funtkionswerte auch an den Gitterpunkten ben¨otigen, in denen wir jetzt die Steigung angeben, und w¨ urden auch Steigungen in den Gitterpunkten erhalten, f¨ ur die wir die Funktionswerte angegeben haben.

3.3.3

Extrema

Neben der Bestimmung von Nullstellen und Ableitungen ist ein Aspekt von Kurvendiskussionen die Bestimmung von Extrema. Diese sind die Punkte mit waagerechter Tangente und damit die Nullstellen der Ableitung. MatLab stellt neben fzero f¨ ur die Nullstellensuche die Funktion fminbnd f¨ ur die Suche nach Extrema. Auch hier wird die Funktion als String eingegeben. Als weitere Argumente werden zwei x-Werte u ¨bergeben, die das Intervall definieren, in dem fminbnd das Minimum sucht: >> fminbnd(’x*x’,-3,3) ←ans = 1.1102e-016 Das Ergebnis weicht um die interne Rechnergenauigkeit von Null ab. Bei der Verwendung von fminbnd ist jedoch Vorsicht geboten, da die Funktion das Minimum innerhalb des betrachteten Intervalls bestimmt – was nicht zwingend ein Extremum der Funktion ist. So liefert auch der Befehl >> fminbnd(’x*x’,-3,-1) ←ans = -1.0000 ein Minimum (eben das in dem Intervall), aber keinen Extremwert.

3.4

Numerische Integration

Auch die Integration von Funktionen erfolgt numerisch. MatLab verf¨ ugt u ¨ber fest eingebaute Funktionen, wir k¨ onnen MatLab jedoch auch verwenden, um verschiedene einfachere ¨ Verfahren wie Mittelpunkts- oder Trapezformel zu benutzen. Eine Ubersicht u ¨ber die Codierung derartiger numerischer Verfahren in verschiedenen Programmiersprachen mit etlichen Beispielen und guten Erl¨ auterungen der Verfahren gibt [1]. c M.-B. Kallenrode

12. Januar 2005

fminbnd

40

KAPITEL 3. FUNKTIONEN

F¨ ur alle folgenden Beispiele ist die Aufgabe die Bestimmung des bestimmten Integrals der Funktion f (x) im Intervall von a bis b: Zb f (x) dx .

I= a

Das Interval [a, b] kann in M Schritte der L¨ange ∆x eingeteilt werden: ∆x = (b − a)/M . Beispiel 9 Als Beispiel werden wir die Funktion f (x) = 3x2 + 2 im Intervall von 2 bis 4 betrachten, d.h. es ist das Integral Z4

(3x2 + 2) dx

I= 2

zu bestimmen. Als Referenz f¨ ur unsere numerischen Verfahren erhalten wir die analytische L¨ osung Z4



4

(3x2 + 2) dx = x3 + 2x

I=

2

= 60 .

2

2

3.4.1

Mittelpunktsformel

Bei Verwendung der Mittelpunktsformel wird das Integral aus Rechtecken der Breite ∆x zusammengesetzt, die H¨ ohe der Rechtecke wird durch den Funktionswert in der Intervallmitte gegeben:   M M X X xk−1 + xk f f (xmk ) ∆x = ∆x . IMP = 2 k=1

k=1

In MatLab k¨ onnen wir dieses Integral auswerten, in dem wir einen Vektor der Intervallmitten erzeugen, den zugeh¨ origen Vektor der Funktionswerte ausrechnen, beide Vektoren punktweise multiplizieren und die Komponenten des sich dabei ergebenden Vektors addieren: >> a=2;b=4;M=100;deltx=(b-a)/M; ←>> x=[a+deltx/2:deltx:b-deltx/2]; y=(3*x.*x+2).*deltx; ←>> I = sum(y) ←I = 59.9998 ¨ Das Ergebnis ist in guter Ubereinstimmung mit dem analytischen Wert aus Beispiel 9. Die Abh¨angigkeit der Genauigkeit des numerischen Ergebnisses von der Zahl M der Schritte bzw. der Schrittweite ∆x ist in folgender Tabelle gegeben: M ∆x I

3.4.2

1 2 58

2 1 59.5

5 0.4 59.92

10 0.2 59.98

20 0.1 59.995

50 0.04 59.9992

100 0.02 59.9998

200 0.01 60.0000

500 0.004 60.0000

Trapezformel

Bei der Trapezformel wird das Integral durch Trapeze der Breite ∆x und der durch die Funktionswerte an den Intervallr¨ andern gegebenen H¨ohen zusammen gesetzt: IT =

M X 1 k=1

2

(f (xk ) + f (xk−1 )) ∆x .

Die Regeln f¨ ur die numerische Integration in MatLab entsprechen denen bei der Mittelpunktsformel. Dort ist nur die zweite Zeile (Vektor der Punkte, an denen der Funktionswert zu nehmen ist und Verfahren der Bestimmung der Fl¨ache) zu ersetzen durch 12. Januar 2005

c M.-B. Kallenrode

3.4. NUMERISCHE INTEGRATION

41

x=[a:deltx:b]; for k=1:M; y(k)=(3*x(k)*x(k) +2 + 3*x(k+1)*x(k+1) + 2)*deltx/2; end ¨ und wir erhalten als Ergebnis 60.0004, ebenfalls in guter Ubereinstimmung mit dem in Beispiel 9 gefundenen analytischen Wert. Das vollst¨andige Verfahren ist im Skript trapez zusammen gefasst. Die Abh¨ angigkeit des Ergebnis von der Schrittweite ist in folgender Tabelle gegeben: M ∆x I

1 2 64

2 1 61

5 0.4 60.16

10 0.2 60.04

20 0.1 60.01

50 0.04 60.0016

100 0.02 60.0004

200 0.01 60.0001

trapez

500 0.004 60.0000

Im oben gegebenen Beispiel f¨ ur die Trapezformel ist die Funktion zwei mal im m-File eingetragen. Hier ist eine alternative Codierung zu empfehlen, bei der die Funktion u ¨ber einen Handle definiert und mit feval ausgewertet wird:

feval

x=[a:deltx:b]; fhandle = @(x) 3*x.*x +2; for k=1:M; y(k) = (feval(fhandle,x(k))+feval(fhandle,x(k+1)))*deltx/2; end Mit dieser Vor¨ ubung k¨ onnen wir satt eines Skripts jetzt eine Funktion trapezfunk schreiben, die die Integration gem¨ aß Trapezformel ausf¨ uhrt: function I = trapezfunk(f,a,b,M) % Trapezformel % Eingabe: Funktion f, Grenzen a und b, Zahl der Schritte M % Ausgabe: Skalar I als Fl¨ ache deltx=(b-a)/M; x=[a:deltx:b]; for k=1:M; y(k) = (feval(f,x(k))+feval(f,x(k+1)))*deltx/2; end I=sum(y); Diese Funktion kann aus einem Skript oder dem MatLab Befehlsfenster heraus aufgerufen werden. In beiden F¨ allen m¨ ussen die zu u ¨bergebenden Parameter vorher vereinbart werden. F¨ ur den Aufruf aus dem Kommandofenster heraus k¨onnte dies wie folgt aussehen: >> a=2; b=4; M=100; fun = @(x) 3*x.*x +2; >> I=trapezfunk(fun,a,b,M) I = 60 Der Vorteil der Darstellung der Trapezformel als eine Funktion ist offensichtlich: jetzt m¨ ussen Integrationsgrenzen, Schrittzahlen oder die zu integrierende Funktion nicht mehr im m-File ge¨ andert werden sondern nur in der ersten Eingabezeile im Befehlsfenster bzw. falls die Funktion trapezfunk aus einem m-File aufgerufen wird in diesem. Auf Grund dieses offensichtlichen Vorteils werden die folgenden numerischen Verfahren stets als Funktionen und mit Hilfe eines Handles zur Beschreibung der zu integrierenden mathematischen Funktion eingef¨ uhrt.

3.4.3

Simpson-Regel

Bei der Simpson-Regel wird der Integrand durch ein Polynom zweiten Grades angen¨ahert. F¨ ur die Teilintegrale werden dazu in einem Intervall ∆x Funktionswerte an den Intervallgrenzen c M.-B. Kallenrode

12. Januar 2005

trapezfunk

42

KAPITEL 3. FUNKTIONEN

xk und xk−1 bestimmt sowie in der Intervallmitte xmk . F¨ ur das Integral ergibt sich dann M

ISimpson =

1X (f (xk−1 ) + 4f (xmk ) + f (xk ))∆x . 6 k=1

simpsonfunk

F¨ ur das Integral aus Beispiel 9 l¨ asst sich das Simpson-Verfahren in folgender MatLabFunktion simpsonfunk realisieren: function I = simpsonfunk(f,a,b,M) % Simpson-Regel % Eingabe: Funktion f, Grenzen a und b, Zahl der Schritte M % Ausgabe: Skalar I als Fl¨ ache deltx=(b-a)/M; x=[a:deltx:b]; xm=[a+deltx/2:deltx:b-deltx/2]; for k=1:M; y(k) = (feval(f,x(k))+4*feval(f,xm(k))+feval(f,x(k+1)))*deltx/6; end I=sum(y); Der Aufruf erfolgt wie bei trapezfunk; zum Vergleich k¨onnen auch beide Funktionen aus einem Skript heraus gestartet werden. Als Ergebnis f¨ ur die Standard-Eingabeparameter erhalten wir 60.000. Das Simpson-Verfahren liefert also bereits f¨ ur recht kleine Schrittzahlen M bzw. recht große Schrittweiten ∆x bessere Ergebnisse als die beiden anderen Verfahren, wie aus der folgenden Tabelle zu erkennen ist: M ∆x I

3.4.4

quad

1 2 60

2 1 60

5 0.4 60

10 0.2 60

20 0.1 60

50 0.04 60.0000

100 0.02 60.0000

200 0.01 60.0000

500 0.004 60.0000

Adaptives Simpson-Verfahren

Durch geeignete Wahl der Schrittweite kann man mit Hilfe des Simpson-Verfahrens sehr genaue Ergebnisse f¨ ur die numerische Integration erzielen. Diese werden, wie bei jedem numerischen Verfahren, mit abnehmender Schrittweite immer besser (sorry, in obigem Beispiel waren sie leider zu gut). Ein sehr feines Gitter bedeutet jedoch auch eine lange Rechenzeit. Dieses feine Gitter ist aber nicht erforderlich in den Bereichen, in denen sich eine Funktion nur schwach ¨ andert: im Extremfall konstanter Intensit¨at in einem Bereich der Funktion reicht ein einziges Rechteck oder Trapez. Allgemein muss die Gitterweite um so feiner sein, je st¨arker sich die Funktion ¨andert. Bei Funktionen, die sich in einigen Bereich stark, in anderen kaum ver¨andern, ist daher die Verwendung eines feinen Gitters u ¨ber den gesamten Integrationsbereich zwar eine sichere aber auch eine rechenintensive Methode. Besser w¨are die Verwendung eines feinen Gitters in einigen Bereichen und eines groben Gitters in anderen. Dann kann bei gleicher Rechengenauigkeit Rechenzeit gespart werden. Eine derartige Anpassung des Integrationsverfahrens an die Funktion ist in MatLab als adpatives Simpson-Verfahren realisiert und wird mit der Funktion quad aufgerufen. Das adaptive Simpson-Verfahren angewandt auf unsere Beispielfunktion ist >> I = quad(f,a,b) ←wobei f wieder die Funktion und a und b die Grenzen geben, in denen diese Funktion integriert werden soll. Beim Aufruf von quad l¨ asst sich eine Toleranz als weiterer Parameter angeben, die MatLab anweist, bis zu welcher Genauigkeit die Integration erfolgen soll (und damit auch, wie zeitintensiv diese sein darf). Statt der Zuweisung der Funktion zu einer Variablen kann die Funktion auch direkt in als Parameter beim Funktionsaufruf u ¨bergeben werden: >> I = quad(’3*x.∧ 2+2’,2,4) ←-

12. Januar 2005

c M.-B. Kallenrode

3.4. NUMERISCHE INTEGRATION

43

Abbildung 3.13: GUI numint zur numerischen Integration

3.4.5

Weitere MatLab Funktionen

F¨ ur die numerische Integration stehen in MatLab abgesehen von den oben gegebenen einfachen ‘Selbstbau-Skripten’ verschiedene Funktionen zur Verf¨ ugung: • quad verwendet ein adaptives Simpson-Verfahren f¨ ur die numerische Integration. • quadl verwendet ein adaptives Lobatto-Verfahren. • dblquad kann f¨ ur Doppelintegrale verwendet werden. Sie basiert auf einem adaptiven Simpson-Verfahren, ist also eine Anwendung von quad. • triplequad hilft bei der numerischen L¨osung von Dreifachintegralen, ebenfalls basierend auf einem adaptiven Simpson-Verfahren. Die MatLab-Funktionen folgen in ihrer Syntax alle dem bei quad besagten; n¨ahere Informationen bietet die MatLab-Hilfe.

3.4.6

GUI zur numerischen Integration

Die Datei numint enth¨ alt ein einfaches GUI zu numerischen Integration von Funktionen, vgl. Abb. 3.13. Die Eingabeparameter (rechts oben) beschr¨anken sich auf die Funktion, die untere ¨ und obere Grenze sowie die Schrittzahl. Uber das Pop-Up Men¨ u k¨onnen die verschiedenen bisher diskutierten Verfahren zur numerischen Integration ausgew¨ahlt werden. Als Ergebnis wird im linken Teil der Funktionsgraph im zu untersuchenden Bereich geplottet, der Wert des sich ergebenden Integrals ist ebenfalls angegeben. Das GUI besteht aus dem eigentlichen Skript numint, einer Routine setintegration, die f¨ ur den Plot und die Zuordnung zu den verschiedenen Integrationsverfahren erforderlich ist, sowie den Funktionen mittelgui, trapezgui und simpsongui, die jeweils die numerischen Verfahren ausf¨ uhren. Diese Funktionen unterscheiden sich von den weiter oben eingef¨ uhrten Funktionen nur durch die Art und Weise, in der die mathematische Funktion u bergeben wird ¨ (Handle vs. String) und damit wie sie ausgewertet wird (feval vs. eval).

3.4.7

Monte-Carlo Integration

Ein ganz anderes Verfahren zur Integration ist die Monte-Carlo Methode. Monte-Carlo Verfahren basieren auf der Verwendung von Zufallszahlen, d.h. sie unterscheiden sich in ihrem Ansatz von den bisher betrachteten Verfahren. Diese waren streng deterministisch und damit reproduzierbar, w¨ ahrend in Monte-Carlo Verfahren ein Zufallselement ins Spiel kommt. Monte-Carlo Verfahren finden in allen Bereichen Anwendung, in denen Prozesse formal eher c M.-B. Kallenrode

12. Januar 2005

numint

44

KAPITEL 3. FUNKTIONEN

auf der Basis von Wahrscheinlichkeitsverteilungen beschrieben werden als in Gleichungen mit genau definierten Parametern. Ein Beispiel f¨ ur ein Monte-Carlo Verfahren ist die Bestimmung von π = 3.14159265358979. Die Zahl π k¨ onnen wir weder als rationale Zahl noch u ¨ber einen funktionalen Zusammenhang einf¨ uhren – außer u ur die Bestimmung der Fl¨ache oder des ¨ber die Bedeutung von π f¨ Umfangs eines Kreises. So gibt π z.B. die Fl¨ache eines Einheitskreises. Wenn wir ein numerisches Verfahren zur Bestimmung der Fl¨ache des Einheitskreises anwenden, bestimmen wir gleichzeitig auch π. Malen wir also einen Einheitskreis auf ein quadratisches St¨ uck Pappe, das diesen Kreis genau umschreibt. Legen wir das Ganze auf den Boden und lassen aus großer H¨ohe Reisk¨ orner auf die Pappe rieseln. Diese treffen statistisch gleichm¨aßig verteilt auf. Daher sollte sich die Zahl der Reisk¨ orner im Kreis zur Gesamtzahl der Reisk¨orner so verhalten, wie die Fl¨ ache des Einheitskreises zur Gesamtfl¨ache: Reisk. in Einheitsk Einheitskreis π = = . Reisk. gesamt Quadrat 4

rand

Durch Ausz¨ ahlen der Reisk¨ orner l¨ asst sich also π bestimmen. Das Verfahren l¨ asst sich auch in MATLAB anwenden. Dazu ben¨otigen wir Zufallszahlen, die mit der Funktion rand erzeugen. Wird rand ohne Argument auf gerufen, so erzeugt es im Intervall (0,1) gleich verteilte Zufallszahlen. Werden Argumente an rand u ¨bergeben, so kann rand z.B. Matrizen aus Zufallszahlen erzeugen: >> rand(3) ←ans = 0.2028 0.2722 0.1987 0.1988 0.6038 0.0153

0.7468 0.4451 0.9318

oder >> rand(1,3) ←ans = 0.7095 0.4289 montecarlopi

0.3046

F¨ ur die Bestimmung der Fl¨ ache des Kreises ben¨otigen wir die folgende Sequenz aus dem Skript montecarlopi: n=1000; z=0; x=2*rand(1,n)-1; y=2*rand(1,n)-1; r= sqrt(x.*x + y.*y); for i=1:n if r(i) =1 zaehl = zaehl + 1; end end vol = (n-zaehl)/n

Literaturverzeichnis [1] Press, W.H., B.P. Flannery, S.A. Teukolsky and W.T. Vetterling, Numerical Recepies in C, Numerical Recepies in Fortran, Numerical Recepies in C++, Cambridge University Press; die C und Fortran Versionen sind zum Download unter http://www.nr.com/ zu finden. 39

12. Januar 2005

c M.-B. Kallenrode