Einführung in Mathematica  - Teil II  -
zur Vorlesung Mathematische Methoden der Physik im WiSe 2013/14

: :    Prof. Dr. Olaf Lechtenfeld und PD Dr. Michael Flohr    : :    13. 11. 2013    : :

Analysis

Mathematica kann exellent mit allem umgehen, was mit Analysis zu tun hat. Differenzieren haben wir bereits in Teil I kennen gelernt. Natürlich kann Mathematica auch integrieren, und zwar so gut, dass, wenn es denn überhaupt eine in elementaren Funktionen ausdrückbare Stammfunktion gibt, sie von Mathematica auch gefunden wird. Und Mathematica löst auch Differentialgleichungen - die fundamentalen Gleichungen in der Physik.

Differentialgleichungen

Analytische Integration (Man beachte die Integrationskonstante C[1].)

DSolve[y'[x] - x y[x] == 0,y[x],x]

Beim Zuweisen zu einer neuen Funktion z ist wiederum der Ersetzungsoperator notwendig; [[1]] “löst” sozusagen die verschachtelnde Liste auf.

z[x_] = y[x]/.DSolve[y'[x] - x y[x] == 0,y[x],x][[1]]

Nun kann, bis auf die verbleibende Konstante, mit der neuen Funktion numerisch gerechnet werden ...

z[1]

N[%]

Die Eleminierung der Konstanten erfolgt, wie gehabt, ü̈ber eine Ersetzung mit einer Transformationsregel (I ist die komplexe Zahl, siehe Teil I):

%/.C[1]-> I

Analytische Integration mit Randbedingung:

DSolve[{y''[x]+ω^2 y[x] == 0,y[0] == 1,y'[0] == 1/2},y[x],x]

Allgemeine Lösung einer Differentialgleichung in Form der sog. “reinen” Funktion (beachte, dass das Argument bei y nicht mit angegeben wird):

DSolve[y''[x]+ω^2 y[x] == 0,y,x]

Lösung für einen bestimmten Fall durch Substitution, vergleiche mit dem ersten Beispiel:

(y[x]/.%)[[1]]

Hinweis: Die Substitutionsanweisung schreibt man vorteilhafterweise in runde Klammern, bevor die Liste aufgelöst wird, denn Mathematica ist sehr logisch: Geschweifte Klammern erzeugten gleich wieder eine neue Liste, für die eine weitere Operation mit [[1]] nötig wäre.

{y[x]/.%%}[[1]][[1]]

Vektoranalysis

Die Vorlesung “Mathematische Methoden der Physik” macht Sie vor allem mit den Techniken der Vektoranalysis vertraut, die Mathematica natürlich auch beherrscht.

Clear[x,y,z,f,g,h]

Früher musste man hier erst ein Zusatzpaket laden, das aber mittlerweile in den Kern der grundlegenden Mathematicai Funktionen, die immer zur Verfügung stehen, integriert wurde. Überhaupt gibt es so viele spezielle Funktionen und Algorithmen in Mathematica, dass es keinen Sinn macht, alle diese bei jedem Programmstart vollständig zu laden. Statt dessen wir nur der als absolut essentiell erachtete Teil immer geladen, den Rest kann man bei Bedarf in Form von Paketen dazuladen, die sich jeweiligen speziellen Themen widmen.

<<VectorAnalysis`

Hinweis: Zuerst sollte man das Koordinatensystems und seine Variablen festlegen (Mathematica gibt zwar Cartesian, d.h. “kartesisch”, aber konsistent zu seinen Groß- und Kleinschreibungsregeln Xx,Yy und Zz vor. Ein anderes mögliches Koordinatensystem ist z. B. Spherical mit den Variablen Rr, Ttheta, Pphi.)

SetCoordinates[Cartesian[x,y,z]]

Definition einer (skalaren) Funktion:

f[x_,y_,z_] := (x^2)(3 - 2(y^2)) - 2 ((x^2) + (y^2))(z^2) + 3 ((y^2) + (z^2))

SIe sollen Äquipotentiallinien in der derzeitigen Hausübung zeichen. Hier plotten wir eine Äquipotentialfläche:

Manipulate[ContourPlot3D[f[x,y,z] == a,{x,-2,2},{y,-2,2},{z,-2,2}],{{a,3},1,5}]

Berechnung des Gradienten von f im oben vorgegebenen Koordinatensystem:

g = Grad[f[x,y,z]]

Achtung: g ist keine Funktion gemäß der fü̈r f gebrauchten Definition, sondern eine “vektorartige” Liste!

Darstellung eines dreidimensionalen Vektorfeldes: Unter den benutzten Darstellungsoptionen sei die universelle, auch in Plot verwendbare AxesLabel hervorgehoben.

VectorPlot3D[g,{x,0,2},{y,0,2},{z,0,2},AxesLabel -> {"x","y","z"},VectorStyle -> "Arrow3D",VectorScale -> {Small,Small,None},VectorColorFunction -> Hue]

Wir berechnen nun die zwei grundlegenden Differentationen für Vektorfelder, die Divergenz und die Rotation. Letztere sollte für ein Gradientenfeld natürlich verschwinden:

h = Div[g]
r = Curl[g]

FullSimplify[%%]

Achtung: h ist ebenfalls keine Funktion gemäß der fü̈r f gebrauchten Definition, sondern eine “skalarartige” (also einelementige) Liste!

Wie man nach der Vereinfachung mittels FullSimplify sieht, ist diese Funktion rotationssymmetrisch:

ContourPlot3D[h == 0,{x,-2,2},{y,-2,2},{z,-2,2}]

Listen

Mathematica arbeitet intern unglaublich stark mit dem Konzept von Listen. Das kann man daher ausnutzen, um sehr effizienten Code zu schreiben, oder manche mathematischen Probleme sehr einfach in Mathematica zu implementieren. Deshalb ist es wichtig, die grundlegenden Befehle zum Hantieren mit Listen zu kennen.

Listenerzeugung

Ein einfaches Beispiel zum Einstieg:

Table[n,{n,0,16}]

Table[n,{n,0,16,2}]

Tabellen sind im wesentlichen geordnete Listen. Und sie können ganz einfach, ineinander verschachtelt werden. Wir kennen das im Grunde schon ähnlich von den Array, die wir in Teil I zur Erzeugung von Vektoren und Matrizen genutzt haben. Table kann also auch eine Liste von Listen generieren, beispielsweise Wertetabellen.

Table[{n,2^n},{n,0,16}]

Listen dieser Art mit numerischen Daten lassen sich unkompliziert graphisch darstellen.

ListPlot[%]

Symolisch geht's mit der Listenerzeugung freilich auch:

Table[(a+b)^n,{n,0,4}]

liste=Expand[%]

Listenmanipulationen

Zunächst wird beispielhaft das erste Listenelement der vorletzten Programmausgabe entfernt ...

Drop[%%,1]

... und dann das erste "extrahiert".

First[%]

Den letzten Listeneintrag erhält man "natürlich" mit Last!

Last[liste]

Ausschneiden von Listenteilen erlaubt Part.

Part[liste,2;;3]

Oder alternativ auf diesem Weg (vergleiche mit dem Extrahieren von Array Elementen in Teil I):

liste[[2;;3]]

Der Zugriff auf ein einzelnes Element erfolgt in gewohnter Manier:

%[[1]]

Table[{n,(1/n)},{n,1,7}]

Ein extrem nützlicher Befehl: Flatten reduziert die beiden Ebenen zuvor generierter Liste auf eine.

Flatten[%]

Die Summe aller Listenelemente läßt sich mit Total berechnen.

Total[%]

Bedingungen

Mathematica ist nicht nur eine Ansammlung zahlloser Befehle und Datenstrukturen zum Umgang mit mathematischen Objekten, sondern eine vollwertige “High-Level” Programmiersprache. Wie jede vernüfntige Programmiersprache bietet es die Möglichkeit, Bedingungen zu formulieren, auszuwerten, und dann in Abhängigkeit vom Ergebnis unterschiedlich weiter zu verfahren.

Abschnittsweise definierte Funktionen

Eine reelle Funktion, die auf zwei Intervallen unterschiedlich definiert ist,

f[x_]:=If[x<0,x^3,x]

Plot[f[x],{x,-1,1}]

sowie deren bestimmtes Integral,

Integrate[f[x],{x,-1,1}]

deren "Stammfunktion"

Integrate[f[x],x]

und deren erste Ableitung

Clear[g]

g[x_]=D[f[x],x]

Wir wollen uns das in einem Plot ansehen. Wir hätten bei diesem Plot allerdings auch gerne eine schöne Legende, die die verschiedenen Funktionen, die wir plotten, auflistet. Es gibt zahllose solche “Sonderwünsche”, für die es in Mathematica bereits Lösungen in Form von Zusatzpaketen gibt. Damit das Programm beim Starten jedoch nicht eine riesige Menge solcher Zusatzfunktionen laden muss, wird nur ein kleiner Teil grundlegender Funktionen geladen, der Rest kann in Form dieser Zusatzpakete, falls benötigt, jederzeit nachgeladen werden. 
So musste das Paket PlotLegends zum Erzeugen einer Legende zuerst geladen werden. Das Laden eines Paketes geht mit folgender Syntax, wobei dieses spezielle Paket inzwischen zum grundlegenden Kern von Mathematica zählt und nicht mehr geladen werden müsste.

<<PlotLegends`

Plot[{f[x],g[x]},{x,-1,1},AxesLabel->{"x"},PlotLegend->{"f(x)","(\[DifferentialD]/(\[DifferentialD]x))f(x)"},LegendPosition->{1,0}]

Weitere Anwendungen

Hier kombinieren wir das bisher gelernte über Listen und Bedingungen in ein paar wenigen Beispielen, die sich um Primzahlen drehen.

primzahlen=Table[Prime[i],{i,20}]

j=Last[primzahlen]

Die Funktion MemberQ durchsucht eine Liste darauf hin, ob ein bestimmtes Element enthalten ist.

DiscretePlot[If[MemberQ[primzahlen,i],1,-1],{i,1,j}]

Durchsuchen ist aber eigentlich gar nicht nötig, da die Liste aufsteigend sortiert ist und aufsteigend "abgearbeitet" wird.

DiscretePlot[If[i==First[primzahlen],primzahlen=Drop[primzahlen,1];1,-1],{i,1,j}]

Achtung: Nach obiger "Abarbeitung" der Primzahlenliste ist diese leer, denn die Liste wurde tatsächlich während der Auswertung bei jedem Schritt um ein Element verkürzt.

primzahlen

Logische Ausdrücke

Mathematica kann in vielfältiger Weise logische Ausdrücke konstruieren und auswerten. Jeder auswertbare logische Ausdruck basiert auf Grundeinheiten, die meist als Vergleich formuliert werden können. Alle Vergleiche können nur zwei Werte annehmen: Wahr oder falsch.

(2^0)==1

2<0

Verknüpfung zweier Vergleiche mit "und"  bzw. "oder"

3>E&&3<Pi

Äquivalente Eingabe:

3 > E ∧3<Pi

3<E||3<Pi

Äquivalente Eingabe:

3>E∨3>Pi

Negation

!(a>b)

Verschachtelte If-Befehle

Integrate[If[x<0,If[x==-1,1,2],3],x]

Schleifenstrukturen und Iteration

Der nächste wichtige Typ von Strukturbefehlen, der Mathematica zu einer vollen und extrem leistungsfähigen Programmiersprache macht, sind Schleifen, mit denen die gleiche Sequenz von Befehlen mehrfach ausgeführt werdden können.

i=0;
While[i<3,Print[i];i=i+1]

Eine klassische "for"-Schleife, wie sie wohl nahezu jede (imperative) Programmiersprache kennt:

For[i=0,i<3,i=i+1,Print[i]]

In Mathematica gibt es auch das einfachere, zu For äquivalente Do,

Do[Print[i],{i,0,2}]

welches von der Struktur beispielsweise analog zu Table ist.

Table[i,{i,0,2}]

Anwendungen

Wozu ist das gut? Hier ein paar wirklich extrem simple Beispiele. Wie alles in Mathematica können auch diese Strukturen beliebig ineinander verschachtelt werden. Es hat auch seinen Grund, dass viele Dinge in Mathematica auf verschiedene Weisen erreicht werden können. So hat jede der verschiedenen Methoden einer For Schleife besondere Vorzüge für bestimmte Anwendungen, was zu schnellerem und effizienterem Code führt, wenn man sich damit auskennt.

For[i=0,i<3,i=i+1,k=i]

k

Eine (Summations)schleife in mathematisch etwas vertrauterer Darstellung

k=Sum[1,{i,1,10}]

Dieselbe sieht "programmiert" wie folgt aus:

k=0;
For[i=1,i≤10,i=i+1,k=k+1]
k

Ein weiteres Beispiel, die Gauß-Summe

k=(∑_(i=1)^10)i

k=0;
For[i=1,i≤10,i=i+1,k=k+i]
k

Am Rande sei erwähnt, dass Mathamtica die Gauß-Summe natürlich auch für eine unbestimmte obere Grenze symbolisch ausrechnen kann.

k=(∑_(i=1)^N)i

Fakultät, iterativ

Und hier noch ein Klassiker, die iterative Definition der Fakultät:

fi[n_]:=Product[m,{m,2,n}]

fi[137]

Und zur Kontrolle hier mit dem in Mathematica eingebauten Befehl

137!

Rekursion

Mathematica kann hervorragend mit rekursiv definierten Strukturen umgehen. Das ist zwar nicht immer die effizientere Lösung, aber oft elegant und in vielen Fällen einfacher zu programmieren

Fakultät, rekursiv

fr[n_]:=n fr[n-1]

fr[1]=1;

fr[137]

Die Funktion Timing gibt die Verarbeitungszeit in Sekunden an.

Timing[fr[10000];]

Um eine Systemvariable wie $RecursionLimit nur in einer Eingabezeile zu verändern, muß die beeinflußte Funktion in Block verschachtelt werden.

Timing[Block[{$RecursionLimit=10004},rf[10000];]]

Hier ein Geschwindigkeitsvergleich der drei Möglichkeiten:

Timing[Block[{$RecursionLimit=10004},rf[10000];]]

Timing[fi[10000];]

Timing[10000!;]

Aha! Die iterative Definition ist die schlechteste. Dies liegt daran, dass bei der iterativen Definition sehr viele Zwischenergebnisse erzeugt werden, was wiederum daran liegt, dass die Product Funktion im Grunde in einer Art For Schleife implementiert ist. Die rekursive Definition hingegen arbeitet “in place”,und erledigt daher die 10.000-fache Multiplikation effizienter. Und die eingebaute Funktion ist ohnehin wesentlich komplizierter programmiert, und hat einen sogenannten “overhead”, weil sie erst einmal prüft, ob das Argument eine einfache Berechnung mit natürlichen Zahlen zulässt, eine symbolisch exakte Berechnung, eine numerische, oder gar keine. 
Ergebnis des Geschwindigkeitsvergleiches. Die eingebauten Funktionen von Mathematica sind hochoptimiert und meist unschlagbar schnell. Dass die rekursive Definition hier gewinnt, liegt daran, dass Mathematica einmal berechnete Werte cached, und das die eingebaute Funktion für die Fakultät wesentlich mehr kann. Dahinter steckt nämlich die Γ-Funktion, die analytische Fortsetzung der Fakultät auf die ganzen komplexen Zahlen.

Und wieder eine Randbemerkung: Die interne Definition der Fakultät ist in der Tat in Wirklichkeit viel allgemeiner, denn sie kann auch auf beliebige reelle (und komplexe!) Zahlen angewandt werden. In Wirklichkeit ist dies die Gammafunktion:

5.5!

(11/2)!