Einführung in Mathematica  - Teil 2 -
zur Vorlesung Mathematische Methoden der Physik im WiSe 2011/12

Norbert Dragon und Michael Flohr (mit Unterstützung von Martin Paech)                                                                                                       16. 12. 2011

Analysis

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]] “lö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, über eine Ersetzung mit einer Transformation- sregel (I ist die komplexe Zahl):

%/.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 Lö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]

Lösung für einen bestimmten Fall durch Substitution, vgl. erstes Beispiel:

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

Hinweis: Die Substitutionsanweisung schreibt man vorteilhafterweise in runde Klammern, bevor die Liste aufgelöst wird, denn Mathematica ist sehr logisch: Geschweifte Klammern erzeugten gleich wieder eine neue Liste, für die eine weitere Operation mit [[1]] nötig wäre.

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

Vektoranalysis

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

Laden eines Paketes für erweitere Funktionen - das kann einen Augenblick dauern ...

<<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 mö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 haben schon Äquipotentiallinien in einer Hausübung zeichen müssen. 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 gemäß der fü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 gemäß der fü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

Listenerzeugung

Ein einfaches Beispiel zum Einstieg:

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

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

Table kann 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[%]

Symbolisch 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]

Ein alternativer Weg geht mittels des uns schon bekannten [[...]]:

liste[[2;;3]]

Der Zugriff auf ein Element erfolgt in gewohnter Manier:

%[[1]]

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

Flatten reduziert die beiden Ebenen zuvor generierter Liste auf eine. Dies ist ein wichter Operator in Mathematica.

Flatten[%]

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

Total[%]

Bedingungen

Abschnittsweise definierte Funktionen

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

Clear[f,g]
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

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

Das Paket PlotLegends muß zum Erzeugen einer Legende geladen werden.

<<PlotLegends`

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

Weitere Anwendungen

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

j=Last[primzahlen]

MemberQ durchsucht die Liste, 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}]

Nach obiger "Abarbeitung" der Primzahlenliste ist diese leer.

primzahlen

Logische Ausdrücke

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

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 (und meist schnellere), 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

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

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

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

Fakultät, iterativ

Clear[fi,fr]

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

fi[137]

Kontrollrechung: Haben wir das richtig programmiert?

137!

Rekursion

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];]]

Geschwindigkeitsvergleich. 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.

Timing[fi[10000];]

Timing[10000!;]