Notebook zu: Mathematik für Informatiker
Band 1

©2001-2018 Gerald Teschl (http://www.mat.univie.ac.at/~gerald/) und Susanne Teschl  (http://staff.technikum-wien.at/~teschl/)

Dieses Notebook ist als Ergänzung zu unserem Buch Mathematik für Informatiker, Band 1, 4. Aufl., Springer 2013, gedacht. Dort finden Sie auch weitere Erklärungen.

Logik und Mengen

Mengen

Mengen können in Mathematica als Listen eingegeben werden:

A={a,b,c}; B={c,d,e};

Vereinigung

Union[A,B]

Durchschnitt

Intersection[A,B]

Zur Vereinfachung (Sortieren und Entfernen von Mehrfachelementen) kann man einfach folgendes verwenden:

Union[{a,c,b,c}]

A={a,c,b,c}; B={c,b,a};

Achtung: Unsere Mengen sind Listen und werden auch also solche behandelt. Um z.B auf Gleichheit zu testen muss man sich daher eines kleinen Tricks bedienen:

{A==B,Union[A]==Union[B]}

Schaltalgebra

Mathematica verwendet False, True für 0, 1  und kennt eine Reihe logischer Verknüpfungen: Negation Not[a] (oder ! a), Und And[a,b]} (oder a && b), Oder Or[a,b] (oder a || b).

a=True;b=False;
{Not[a],And[a,b],Or[a,b]}

{!a,a&&b,a || b}

Clear[a,b]

Mit folgendem Programm können wir leicht Wertetabellen erstellen:

LogicTable[f_,v_List]:=Module[{n=Length[v],tabl,values,rules},
tabl={Flatten[{v,f}]};
Do[
values= IntegerDigits[i,2,n] /. {0->False,1->True};
rules=Table[Rule[v[[i]],values[[i]] ],{i,n}];
tabl=Append[tabl,Flatten[{values,f /. rules}]];
,{i,0,(2^n)-1}];
TableForm[tabl]
]

Grübeln Sie nicht darüber, wie dieses Programm funktioniert, sondern rufen Sie es einfach mit einem logischem Ausdruck (oder einer Liste von logischen Ausdrücken) und einer Liste der Variablen auf:

LogicTable[{!a,!(a||b),!a&&!b},{a,b}]

Mathematica kann übrigens auch logische Ausdrücke vereinfachen :

LogicalExpand[!a&&!b||!a&&b||a&&b]

LogicalExpand[Implies[a,b]]

Zahlenmengen und Zahlensysteme

Approximation von (Sqrt[2])

Das auf im Buch beschriebene Programm zur Annäherung der Wurzel aus 2 kann mit Mathematica wie folgt  implementiert werden:

d[q_]:=Module[{p=q},
While[((((p/q)))^2)<2,p=p+1];
{((p-1)/q),(p/q)}]

Der Befehl Module fasst mehrere Befehle zusammen. Das erste Argument ist dabei eine Liste von lokalen Variablen.

d[100]

Clear[d]

Ungleichungen

Mathematica kann auch mit Ungleichungen umgehen. Der Simplify-Befehl kann zum Überprüfen von Ungleichungen verwendet werden :

Simplify[(x/((x^2)+(y^2)))<(1/y),x>0&&y>0]

Mit Reduce können Ungleichungen sogar aufgelöst werden:

Reduce[1-(x^2)>0,x,Reals]

Komplexe Zahlen

Mit komplexen Zahlen rechnet man folgendermaßen :

z1=1+2I;z2=3- I;(z1/z2)

Die imaginäre Einheit kann entweder über die Tastatur (als großes I) oder über die Palette (als i) eingegeben werden. Real - bzw. Imaginärteil, komplexe Konjugation und Absolutbetrag erhält man mit

{Re[z1],Im[z1],Conjugate[z1],Abs[z1]}

Manchmal muss man mit dem Befehl ComplexExpand noch nachhelfen, damit das Ergebnis in Real - und Imaginärteil aufgespalten wird :

(Sqrt[1+I (Sqrt[3])])

ComplexExpand[%]

Mehr noch, Mathematica geht bei allen Variablen standardmäßig davon aus, dass sie komplexwertig sind. Deshalb wird zum Beispiel der Ausdruck

Simplify[((Sqrt[a b])/(Sqrt[a]))]

nicht zu (Sqrt[b]) vereinfacht, denn das stimmt im Allgemeinen nur für a> 0! Abhilfe schafft in so einem Fall die Möglichkeit, im Simplify-Befehl die Zusatzinformation a > 0 zu geben :

Simplify[((Sqrt[a b])/(Sqrt[a])), a>0]

Achtung : Nicht alle Wurzelgesetzte bleiben im Komplexen richtig :

(Sqrt[-I])(Sqrt[-I])==(Sqrt[(-I)(-I)])

Clear[z1,z2]

Summen- und Produktzeichen

Das Summenzeichen kann entweder direkt über die Palette eingegeben werden,

(∑_(k=1)^n)k

oder auch als Sum[k, {k, 1, n}]. Wie Sie sehen, wertet Mathematica (falls möglich) Summen sofort aus. Analog für Produkte mit Product[k, {k, 1, n}] oder über die Palette:

(∏_(k=1)^n)k

Umwandlung zwischen Zahlensystemen

Der Mathematica-Befehl BaseForm[x, b] wandelt die Dezimalzahl x in eine Zahlendarstellung mit Basis b um. Zum Beispiel wird die Zahl (Subscript[(0.1), 10]) mit

BaseForm[0.1,2]

vom Dezimalsystem ins Dualsystem umgewandelt. Die Umwandlung einer Zahl x von einem System mit Basis b ins Dezimalsystem erhält man mit bˆˆx:

16^^FAD

wandelt die Hexadezimalzahl (Subscript[(FAD), 16]) ins Dezimalsystem um oder

8^^67

wandelt die Oktalzahl (Subscript[(67), 8]) ins Dezimalsystem um.

Teilbarkeit und Primzahlen

Mit dem Mathematica-Befehl PrimeQ kann man feststellen, ob eine Zahl eine Primzahl ist :

PrimeQ[4]

Das „Q "steht dabei für „question ".Mit einer Do-Schleife können wir zum Beispiel die Liste aller Primzahlen bis 5 ausgeben lassen :

Do[
If[PrimeQ[n],Print[n]],
{n,1,5}];

Der Befehl zur Primfaktorzerlegung heißt FactorInteger und liefert die Liste aller Primfaktoren, zusammen mit der zugehörigen Vielfachheit:

FactorInteger[180]

also 180=(2^2)·(3^2)·(5^1). Der größte gemeinsame Teiler kann mit dem Befehl GCD (“greatest common divisor”) berechnet werden:

GCD[75,38]

Die Zahlen 75 und 38 sind also teilerfremd.Der Rest der Division einer ganzen Zahl x durch die natürliche Zahl m wird mit Mod[x, m] erhalten :

Mod[22,5]

Der Quotient der Division wird mit

Quotient[22,5]

berechnet. Also ist 22 = 4·5 + 2.

Elementare Begriffe der Zahlentheorie

Rest modulo m

Der Rest von a modulo m wird mit Mod[a, m] berechnet :

Mod[17,5]

Multiplikatives Inverses

Das multiplikative Inverse 1 in (ℤ^m) kann mit PowerMod[e, −1, m] berechnet werden :

PowerMod[4,-1,9]

also 1 = 7 in (ℤ^9). Allgemein berechnet PowerMod[e, k, m] die Potenz (e^k) modulo m. In Mathematica kann man nicht einfach e − 1 schreiben, denn woher soll das arme Programm wissen, ob Sie in ℝ oder in (ℤ^m) rechnen wollen!

Euklidscher Algorithmus

Der Euklid' sche Algorithmus kann wie folgt implementiert werden :

Euklid[a_Integer,b_Integer]:= Module[{r=a,rr=b},
While[ rr!= 0, {r,rr}={rr,Mod[r,rr]}];
r];

Euklid[75,38]

(„ != "bedeutet „ungleich "). Der ggT (75, 38) ist also 1. Natürlich hätten wir auch gleich den internen Mathematica-Befehl GCD für den größten gemeinsamen Teiler verwenden können :

GCD[75,38]

Alternativ können wir den Algorithmus auch rekursiv implementieren:

Euklid[a_Integer,0]=a;Euklid[a_Integer,b_Integer]:=Euklid[b,Mod[a,b]];

Euklid[75,38]

Analog kann der erweiterte Euklid’sche Algorithmus so programmiert werden:

ExtendedEuklid[a_Integer,b_Integer]:= Module[{r=a,rr=b,x=1,xx=0,y=0,yy=1, Q},
While[ rr != 0, 
Q=Quotient[r,rr];
{r,rr,x,xx,y,yy}={rr,Mod[r,rr],xx,x-Q xx,yy,y-Q yy}
];
{r,x,y}];

ExtendedEuklid[75,38]

Ausgegeben werden also der ggT (75, 38) = 1 und ganzzahlige Lösungen x = −1 und y = 2 der diophantischen Gleichung 75 x + 38 y = ggT (75, 38). Wieder gibt es einen internen Mathematica-Befehl dazu : ExtendedGCD[a, b] gibt die Liste {g, {x, y}} aus, wobei g = ggT (a, b) und x, y ganzzahlige Lösungen von ax + by = g sind.

ExtendedGCD[75,38]

RSA

Die Verschlüsselung mittels RSA-Algorithmus ist natürlich zu aufwendig, um sie von Hand durchzuführen (aber auch für langsame Computer – Stichwort Chipkarten – kann die Geschwindigkeit bei RSA zu einem Problem werden). Wir wollen uns hier von Mathematica helfen lassen : Der Befehl ToCharacterCode wandelt ein Zeichen (Buchstabe, Ziffer, . . .) und sogar eine ganze Zeichenkette in eine Liste von Zahlen gemäß dem ASCII-Code um :

ToCharacterCode["KLEOPATRA"]

Da im ASCII-Code der Buchstabe A der Zahl 65, B der Zahl 66, usw. entspricht, müssen wir – um die gewünschte Zuordnung A = 0, B = 1 usw.zu erhalten – noch 65 subtrahieren :

x=%-65

Das ist nun der in Zahlen codierte Klartext, der verschlüsselt werden soll. Für die Verschlüsselung von x benötigen wir den öffentlichen Schlüssel

n= 1147; e=29;

Nun können wir mit der Vorschrift y = (x^e) (mod n) verschlüsseln :

y=PowerMod[x,e,n]

Hier ist PowerMod[x, e, n] eine effektivere Variante von Mod[(x^e), n]. Der Empfänger kann mit dem geheimen Schlüssel d und der Vorschrift x = (y^d) (mod n) entschlüsseln:

d=149;PowerMod[y,d,n]

Bei unserem kurzen Spielzeugschlüssel ist es natürlich für einen Angreifer kein Problem den Algorithmus zu knacken, d.h.n zu faktorisieren :

FactorInteger[n]

zerlegt den Modul n = 1147 in seine Primfaktoren p = 31 und q = 37. Damit können wir m berechnen :

m=(31-1)(37-1)

m kann auch alternativ mittels m = EulerPhi[n] berechnet werden:

EulerPhi[n]

Der geheime Schlüssel d ist nun die Lösung der Gleichung e d = 1 (mod m), wobei e = 29 der öffentliche Schlüssel ist und m = 1080 gerade vom Angreifer gefunden wurde. d kann mit dem Befehl PowerMod berechnet werden :

PowerMod[e,-1,m]

Clear[x,y,n,e,d,m]

Chinesischer Restsatz

Das System von Kongruenzen x = (Subscript[a, 1]) (mod (Subscript[m, 1])), ..., x = (Subscript[a, k]) (mod (Subscript[m, k])) kann mit dem Befehl ChineseRemainder[{(Subscript[a, 1]), . . . , (Subscript[a, k])}, {(Subscript[m, 1]), . . . , (Subscript[m, k])}] gelöst werden :

ChineseRemainder[{2,3,2},{3, 5, 7}]

Ausgegeben wird die kleinste nichtnegative Lösung x, hier x = 23.

Polynomringe und endliche Körper

Faktorisierung

Polynome können mit dem Befehl

Factor[x^3-1]

faktorisiert und mit dem Befehl

Expand[%]

ausmultipliziert werden. Eine Faktorisierung erfolgt nur, wenn die Nullstellen rationale Zahlen sind. Beispielsweise wird die Faktorisierung (x^2) − 2 = (x-(Sqrt[2]))(x-+(Sqrt[2]))  von Mathematica nicht durchgeführt:

Factor[(x^2)-2]

Factor[(x^2)-2,Extension->(Sqrt[2])]

Polynomdivision

Den Quotienten s(x) der Polynomdivision p(x) = s(x) q(x) + r(x) können wir mit PolynomialQuotient[p (x), q (x), x]

PolynomialQuotient[3x^4+x^3-2 x,x^2+1,x]

und den Rest r(x) mit PolynomialRemainder[p(x), q(x), x]

PolynomialRemainder[3x^4+x^3-2 x,x^2+1,x]

berechnen.

Euklid'scher Algorithmus

In Mathematica können wir den Euklid' schen Algorithmus analog wie für ganze Zahlen definieren :

PolynomialEuklid[p_,q_]:= Module[{r=p,rr=q},
While[ rr=!= 0, {r,rr}={rr,PolynomialRemainder[r,rr,x]}];
(r/(CoefficientList[r,x][[-1]]))];

(Im letzten Schritt wird dabei nur noch der höchste Koeffizient von r auf eins normiert.)

PolynomialEuklid[(x^3)-2x+1,(x^2)-1]

Wir hätten übrigens auch gleich den internen Befehl PolynomialGCD[a(x), b(x), x] verwenden können. Ich wollte aber zeigen, dass bei der Implementierung nur die Division ganzer Zahlen durch die Polynomdivision ersetzt werden muss. Außerdem normiert Mathematica den ggT nicht : So gibt es als ggT von 5 (x3 − 2 x + 1) und 5 (x2 − 1) das Polynom 5 (x − 1) anstelle des zugehörigen normierten Polynoms x − 1 aus.

Analog wird der erweiterte Euklid' sche Algorithmus implementiert :

PolynomialExtendedEuklid[p_,q_]:= Module[{r=p,rr=q,s=1,ss=0,t=0,tt=1,Q},
While[ rr =!= 0, 
Q=PolynomialQuotient[r,rr,x];
{r,rr,s,ss,t,tt}={rr,PolynomialRemainder[r,rr,x],ss,s-Q ss,tt,t-Q tt}
];
(({r,s,t})/(CoefficientList[r,x][[-1]]))];

Als Ergebnis werden also die drei Polynome r(x), s(x), t(x) ausgegeben, wobei r(x) der größte gemeinsame Teiler von p(x) und q(x) ist und s(x), t(x) die Beziehung s(x)p (x) + t(x) q(x) = r(x) erfüllen. Damit berechnen wir

PolynomialExtendedEuklid[(x^3)-2x+1,(x^2)-1]

Also sind r(x) = x − 1, s(x) = −1 und t(x) = x.

Die analogen internen Befehle lauten:

PolynomialGCD[(x^3)-2x+1,(x^2)-1]

PolynomialExtendedGCD[(x^3)-2x+1,(x^2)-1]

Multiplikatives Inverses

In Mathematica kann der Befehl PolynomialRemainder[a (((x))^(-1)), m(x), x, Modulus → p}] verwendet werden, um das multiplikative Inverse von a(x) in (Subscript[ℤ, p])(Subscript[[x], m(x)]) zu berechnen :

PolynomialRemainder[(((x+1))^(-1)),x^2+x+1,x,Modulus->2]

Faktorisierung

Mathematica kann natürlich auch Polynome aus (Subscript[ℤ, m])[x] faktorisieren. Ein Beispiel aus (Subscript[ℤ, 2])[x] :

Factor[(x^8)+(x^4)+(x^3)+x+1, Modulus->2]

Das Polynom (x^8) + (x^4) + (x^3) + x + 1 ist also über (Subscript[ℤ, 2]) irreduzibel. Über  (Subscript[ℤ, 3]) ist es reduzibel:

Factor[(x^8)+(x^4)+(x^3)+x+1, Modulus->3]

Folgen und Reihen

Folgen

Mit Mathematica können wir eine Folge einfach als eine Funktion definieren,

a[n_]:=(1/(2^n))

und damit leicht Folgenglieder berechnen:

Table[a[n],{n,4}]

Analog können rekursiv definierte Folgen eingegeben werden, wie hier die Heron' sche Folge mit Startwert (Subscript[a, 1]) = 2,

a[1]=2.;a[n_]:=(1/2)(a[n-1]+(2/(a[n-1])))

und die Folgenglieder berechnet werden:

Table[a[n],{n,5}]

Bei rekursiven Folgen muss aber unbedingt die Definition mit “:= “erfolgen, damit die rechte Seite der Definition nicht sofort, sondern erst bei Angabe einer konkreten Zahl n ausgewertet wird. Außerdem darf auch der Startwert a1 nicht vergessen werden! Wenn Mathematica (Subscript[a, n]) berechnet, so drückt es gemäß dem Bildungsgesetz (Subscript[a, n]) durch (Subscript[a, n-1]) aus. Das geschieht so lange, bis es auf den Startwert trifft. Wurde dieser nicht angegeben (oder wurde der Doppelpunkt vergessen), so schicken Sie Mathematica in eine Endlosschleife. Sie können diese durch den Menüpunkt “Kernel/Abort Evaluation” oder durch die Tastenkombination “ALT“+“.”abbrechen.

Die Folge für die Euler' sche Zahl lautet

a[n_]:=(((1+(1/n)))^n)

und die im Buch gesuchten Folgenglieder sind:

N[{a[2],a[3],a[10],a[100],a[1000],a[10000]}]

Erinnern Sie sich, dass der Befehl N[ ...] bewirkt, dass das Ergebnis numerisch angezeigt wird. Der Grenzwert einer Folge kann mit dem Befehl Limit berechnet werden :

Limit[a[n],n->∞]

Der Pfeil →, das Symbol ∞ und die Konstante e werden in Mathematica entweder über die Palette oder über die Tastatur als ->, Infinity bzw. E eingegeben.

Einen numerische Wert von e auf 30 Stellen genau erhalten wir mit dem Befehl N[E, 30].

N[E,30]

Clear[a]

Heron'sche Folge

Ein Programm zur Berechnung von Wurzeln mithilfe der Heron' schen Folge kann in Mathematica so implementiert werden:

MySqrt[x_]:=Module[{h=N[x],hh=N[x]+1},
While[hh-h>0,hh=h;h=(1/2)(h+(x/h))];
hh]

Dabei wird h mit dem numerischen Wert von x initialisiert (rechnen wir symbolisch, so würden wir endlos iterieren; -). Dann wird in jedem Schritt der letzte Wert als hh gespeichert und ein neuer Wert h berechnet, bis sich der Wert nicht mehr ändert.

MySqrt[2]

Reihen

Mit Mathematica berechnet man den Grenzwert einer Reihe mit:

(∑_(k=0)^∞)(1/(2^k))

Das kann auch ohne Palette als

Sum[1/2^k,{k,0,Infinity}]

eingegeben werden. Die n-te Teilsumme erhalten wir mit

s[n_]=(∑_(k=1)^n)(1/(k^2));

bzw.numerische Werte davon mit

N[{s[1],s[3],s[10], s[100], s[1000], s[10000]}]

Der Grenzwert dieser Reihe ist

(∑_(k=1)^∞)(1/(k^2))

bzw. als Kommazahl ausgegeben:

N[%]

Kombinatorik

Fakultät und Bimomialkoeffizient

In Mathematica wird die Fakultät mit

6!

und der Binomialkoeffizient mit

Binomial[45,6]

berechnet.

Rekursionen und Wachstum von Algorithmen

Iteration und Chaos

Um Folgenglieder einer rekursiv definierten Folge effektiv zu berechnen, bietet sich eine Do-Schleife an:

f[μ_,x_]:=x+μ(1-x)x;
F[x_,μ_,n_]:=Module[{xx=x},Do[xx=f[μ,xx],{n}];xx];

Hier wurde die Rekursion (Subscript[x, n]) = (Subscript[x, n]) − 1 + μ (1 − (Subscript[x, n-1])) (Subscript[x, n-1]) definiert.Den Buchstaben μ können wir (wie jeden griechischen Buchstaben) mithilfe der Palette oder über die Tastatur mit der Tastenkombination “ESC m ESC” eingeben.

Noch schneller, aber dafür auch noch kryptischer, würde es mit folgender Variante gehen:

f[μ_,x_]:=x+μ(1-x)x;
F[x_,μ_,n_]:=Nest[f[μ,#1]&,x,n];

Folgenglieder kann man sich nun bequem mit

Table[F[0.1,2.2,n], {n,1,10}]

ausgeben lassen. Ein “Spinnennetz” entsteht mit folgendem kleinen Programm :

ShowWeb[f_,xstart_,nmax_,opts___]:=Module[{x,xmin,xmax,delta,graph,web,lines},
lines=Type /. {opts,Type->Line};
x[0]:=xstart;
x[n_]:=x[n]=f[x[n-1]];
If[x[0]==x[1],Return[Print["Fixed Point!"]]];
web=Flatten[Table[{{x[n],x[n]},{x[n],x[n+1]}},{n,0,nmax}],1];
xmax=Max[web];xmin=Min[web];delta=0.1(xmax-xmin);
graph=Plot[{f[x],x},{x,xmin-delta,xmax+delta}];
Show[graph,Graphics[lines[web]]]
] /;(IntegerQ[nmax]&&Positive[nmax])

Gehen wir es kurz Schritt für Schritt durch. Es übernimmt als Argument die zu iterierende Funktion f, den Startwert xstart, die Anzahl der durchzuführenden Iterationen nmax und weitere Optionen opts (der dreifache Unterstrich steht für eine beliebige Anzahl (inkl. keiner) weiterer Argumente). Als Erstes wird in der Variable lines der Verbindungstyp festgelegt. Der Defaultwert ist Line und kann mit einer optionalen Ersetzungsregel Type → wert (als letztes Argument unseres Programms) geändert werden. Nun werden rekursiv die zu verbindenden Punkte ((Subscript[x, n]), (Subscript[x, n])) und ((Subscript[x, n]), (Subscript[x, n+1])) berechnet (der Flatten-Befehl wirft nur ein paar überschüssige Klammern weg). Über das Minimum bzw. Maximum der berechneten Werte (Subscript[x, n]) wird der darzustellende Bereich ermittelt (dieser Wert wird noch um 10 % vergrößert, damit nichts Wertvolles abgeschnitten wird). Nun wird noch die Funktion f(x) zusammen mit x gezeichnet und die Grafik mit den durch Linien verbundenen Punkten zusammengefügt.

Wenn wir die Linien durch Pfeile ersetzen möchten, dann können wir als Verbindungstyp folgende Funktion, die eine Liste von Punkten durch Pfeile verbindet, verwenden:

Arrows[l_List]:={Arrowheads[Small],Table[Arrow[{l[[i]],l[[i+1]]}],{i,Length[l]-1}]};

Setzen wir

f[x_]:=x+μ x(1-x);

dann zeichnet

μ=1;ShowWeb[f,0.1,6,Type->Arrows]

das Spinnennetz der ersten 12 Folgenglieder für μ = 1 mit dem Startwert (Subscript[x, 0]) = 0.1, das im Buch abgebildet ist.

Einige  weitere Grafiken:

μ=2.2;ShowWeb[f,0.7,6,Type->Arrows]

μ=3;ShowWeb[f,0.7,50]

Clear[f,F]

Lösung einer Rekursion

In Mathematica können Rekursionen mit dem Befehl RSolve gelöst werden.

RSolve[{a[n]==3a[n-1]-4,a[0]==5},a[n],n]

Analog für Rekursionen höherer Ordnung:

RSolve[{a[n]==a[n-1]+6a[n-2], a[0]==1, a[1]==8},a[n],n]

Landausymbol

Zuerst die schlechte Nachricht : Mathematica kann nicht mit Landausymbolen rechnen (es kennt zwar das Landausymbol und verwendet es bei der Darstellung von Potenzreihen, kann aber nur sehr eingeschränkt damit rechnen). Nun die gute Nachricht : Mathematica ist eine Programmiersprache, also kann uns niemand daran hindern, Mathematica etwas auf die Sprünge zu helfen, und dem Programm die Regeln für das Landausymbol beizubringen!

Wir definieren zunächst die Funktion LandauO für das Landausymbol:

LandauO[f_+g_]:=LandauO[If[tmp==0,g,f] ]/; FreeQ[tmp=Limit[Abs[(f/g)],n->∞],n];
LandauO[c_  f_]:= LandauO[f] /;FreeQ[c,n];

Die erste Definition LandauO[f_ + g_] versucht festzustellen, welcher der beiden Summanden schneller wächst. Dazu wird der Grenzwert (lim_(n→∞)) = ((| f | )/(|g|)) berechnet und getestet, ob es Mathematica gelungen ist, den Grenzwert zu berechnen (in diesem Fall darf das Ergebnis kein n mehr enthalten). Die bedingte Definition F[x_] := expr /; Bedingung wendet die Definition für F nur an, falls die Bedingung erfüllt ist. Die zweite Definition wirft Konstanten weg.

Jetzt noch eine handliche Funktion

LandauExpand[expr_]:= LandauO[Expand[expr]]/.Log[x_]/; PolynomialQ[x,n]->Log[n];

die den übergebenen Ausdruck expandiert (alle Produkte ausmultipliziert) und berücksichtigt, dass log((Subscript[p, k]) (n^k) + (Subscript[p, k-1]) (n^(k-1)) − 1 + . . .+(Subscript[p, 0])) = O (log (n)) ist. Die Ersetzungsregel Alt /; Bedingung → Neu wird nur ausgeführt, wenn die Bedingung erfüllt ist.

Nun unser erster Test:

LandauExpand[((n+1)/((n^2)+1))+Exp[-n+3]+11 n Log[3 n+1]]

Sieht doch ganz gut aus! Ganz ehrlich, hätten Sie das mit der Hand auch geschafft?

Vektorräume

Vektoren

In Mathematica können Vektoren einfach mit geschwungenen Klammern eingegeben werden :

a={1,2,0,-1}; b={3,0,1,2};

Wir erhalten dann zum Beispiel 3 a − b mit

3a-b

Die Länge eines Vektors wird mit

Norm[a]

berechnet.

Clear[a,b]

Matrizen und Lineare Abbildungen

Matrizen

Matrizen werden in Mathematica als verschachtelte Listen eingegeben:

A={{1,3,2},{-1,4,5}}; B={{-2,4,0},{6,2,-8}};

Manche Großbuchstaben (wie z.B.C, D oder N) sind in Mathematica bereits vordefiniert, und können daher nicht als Variablennamen verwendet werden.

Summe und Multiplikation mit einem Skalar werden wie erwartet gebildet:

2(A+B)

Das können wir noch in die vertraute Matrixschreibweise bringen:

MatrixForm[%]

Wir können sogar

$PrePrint = MatrixForm;

setzen, dann wird jedes Ergebnis automatisch in Matrixform dargestellt.

Matrizen können mithilfe einer Palette bereits in Matrixform eingegeben werden. Von der Palette wird dabei eine (2, 2)-Matrix als Schablone vorgegeben; mit den Tastenkombinationen “Ctrl”+“Enter” bzw. “Ctrl”+“, “ kann eine Zeile bzw.eine Spalte hinzugefügt werden. Alternativ können Sie eine Schablone mit der gewünschten Zeilen und Spaltenanzahl auch mit der rechten Maustaste über das Menü “Create Table/Matrix/Palette” eingeben (engl. : Spalte = column, Zeile = row). Mit der Tabulatortaste kommen Sie in der Schablone von einem Platzhalter zum nächsten.

Die Dimension einer Matrix kann mit dem Befehl Dimensions überprüft werden:

Dimensions[A]

Die transponierte Matrix erhalten wir mit

Transpose[(({{3, 2}, {5, -1}, {0, 7}}))]//MatrixForm

Die (n, n)-Einheitsmatrix wird mit IdentityMatrix[n] erhalten:

IdentityMatrix[2]//MatrixForm

Eine Diagonalmatrix diag((Subscript[a, 1]), . . . , (Subscript[a, n])) erhalten Sie mit DiagonalMatrix[{ (Subscript[a, 1]), . . . , (Subscript[a, n])}] und das Kronecker Delta (Subscript[δ, jk]) mit KroneckerDelta[j, k].

DiagonalMatrix[{1,2,3}] //MatrixForm

KroneckerDelta[1,2]

Multiplikation von Matrizen

Achtung : Das Matrixprodukt wird mit einem Punkt “.” (und nicht mit dem üblichen Mal-Symbol “∗ “ oder mit Leerzeichen) berechnet!

A=(({{4, 2, 0}, {-1, 3, 5}}));B=(({{2, 1}, {3, 7}, {1, 0}}));  A.B//MatrixForm

Die Eingabe A ∗ B (oder A B) würde als Resultat eine Matrix liefern, die durch elementweises Ausmultiplizieren von A und B entsteht. In diesem Sinn muss das Quadrat (M^2) der Matrix

M = (({{7, 1}, {2, 5}}));

mit Punkt

M.M//MatrixForm

eingegeben werden.Im Gegensatz dazu wird

M^2//MatrixForm

von Mathematica wieder als M ∗ M (elementweises Quadrieren) interpretiert. Höhere Potenzen können wir effizient mit MatrixPower[A, n] berechnen.

MatrixPower[M,2]//MatrixForm

Achtung : Die Eingabe A = {{1, 2}, {3, 4}} // MatrixForm kann zu unerwarteten Ergebnissen führen, da MatrixForm vor der Zuweisung ( = ) angewendet wird:

A={{1,2},{3,4}}//MatrixForm

A+M

Um das zu vermeiden, verwenden Sie Klammern : (A = {{1, 2}, {3, 4}}) // MatrixForm.

Clear[A,B,M]

Eingabe von Gleichungen mit Matrizen

Lineare Gleichungen können in der Form Ax = b eingegeben

A={{5,3},{2,-1}};x={x1,x2};b={2,3};
A.x==b

und gelöst werden:

Solve[%,{x1,x2}]

Clear[A,B,b,x]

Inverse Matrix

Die inverse Matrix wird mit dem Befehl Inverse berechnet

A={{2,4},{-1,3}};Inverse[A] //MatrixForm

Matrizen größerer Dimension sind für Mathematica kein Problem:

Inverse[(({{0.30, 0.59, 0.11}, {0.60, -0.28, -0.32}, {0.21, -0.52, 0.31}}))]//MatrixForm

Clear[A]

$PrePrint=.;

Lineare Gleichungen

Gauß-Algorithmus

Wir wissen bereits, dass in Mathematica der Befehl Solve zur Lösung eines Gleichungssystems zur Verfügung steht:

Solve[{y +3z== 2, x  +2y+5z==0, 2x +5y + 13z == 2}, {x, y, z}]

Alternativ kann Mathematica ein lineares Gleichungssystem auch auf die Endform des Gauß-Jordan-Algorithmus (reduzierte Zeilenstufenform) bringen. Wir geben dazu die erweiterte Koeffizientenmatrix ein, und verwenden dann die Anweisung RowReduce. Für das Beispiel aus dem Buch ergibt sich dann:

A = (({{1, -1, -1, 0}, {-1, 0, -3, -11}, {4, -1, 2, 15}}));

RowReduce[A]//MatrixForm

Einige weitere Beispiele aus dem Buch:

RowReduce[(({{1, -1, 7, 6}, {-1, 4, -13, 3}, {2, 1, 8, 17}}))]//MatrixForm

RowReduce[(({{0, 1, 3, 2}, {1, 2, 5, 0}, {2, 5, 13, 2}}))]//MatrixForm

RowReduce[(({{1, 2, -1, 3}, {2, 4, -2, 6}, {3, 6, -3, 9}}))]//MatrixForm

Clear[A]

Rang, Kern und Bild

In Mathematica kann eine Basis für den Kern mit dem Befehl NullSpace[A] berechnet werden.

A= {{1,1,2}, {0,1,1}, {1,0,1}};NullSpace[A]

Eine Basis für das Bild von A kann zwar nicht direkt berechnet werden, aber es gibt einen Trick : Um die lineare Hülle der Spaltenvektoren nicht zu verändern, müssten wir statt elementarer Zeilenumformungen elementare Spaltenumformungen machen.Der Befehl RowReduce führt aber nur elementare Zeilenumformungen durch. Da aber Transponieren Zeilen und Spalten vertauscht, ist eine elementare Spaltenumformung von A gleich einer elementaren Zeilenumformung von (A^T). Somit erhalten wir mit dem Befehl

Transpose[RowReduce[Transpose[A]]]//MatrixForm

Vektoren, die das Bild aufspannen.Der letzte Vektor ist der Nullvektor und kann natürlich weggeworfen werden. Die beiden ersten Vektoren (1, 0, 1) und (0, 1, −1) sind linear unabhängig und bilden somit eine Basis für das Bild. Daher ist insbesondere die Dimension des Bildes (der Rang der Matrix) gleich 2. Das kann auch direkt mit

MatrixRank[A]

berechnet werden.

Clear[A]

Determinante

Determinanten werden mit dem Befehl Det berechnet:

A = (({{1, -1, -1}, {-1, 0, -3}, {4, -1, 2}}));Det[A]

A = (({{0, 1, 3}, {1, 2, 5}, {2, 5, 13}}));Det[A]

Clear[A]

Lineare Optimierung

Lineare Ungleichungen

Mit dem Befehl RegionPlot können Ungleichungen dargestellt werden:

RegionPlot[x≥0&&y≥0&&x+y≤20&&x+y≥10&&y≤2 x,{x,0,20},{y,0,20}]

Lineare Optimierung

Mit Mathematica kann die Lösung eines linearen Optimierungsproblems mit dem Befehl LinearProgramming[c, A, b] gefunden werden. Im Argument des Befehls werden zuerst die Koeffizienten c = ((Subscript[c, 1]), ..., (Subscript[c, n])) der Zielfunktion f eingegeben (da die Konstante d für das Finden des Maximums unerheblich ist, wird sie weggelassen); danach die Koeffizientenmatrix A = ((Subscript[a, ij])) aller Variablen, und zuletzt die zugehörigen rechten Seiten b = ((Subscript[b, 1]), ..., (Subscript[b, m])) der Ungleichungen. Es wird dabei immer das Minimum von f gesucht. Benötigt man das Maximum, so ist einfach f durch − f zu ersetzen. Zum Beispiel:

LinearProgramming[{4,3},(({{-1, -1}, {1, 1}, {2, -1}})),{-20,10,0}]

Wir haben hier im Vergleich zur Angabe aus dem Buch zuvor alle Ungleichungszeichen auf die Richtung >= gebracht und das Vorzeichen der Zielfunktion geändert, da wir ja das Maximum suchen.

LinearProgramming[{-3,-4},(({{-1, -2}, {-2, -1}})),{-20,-28}]

Wenn man zu faul ist, das Problem auf die Standardform zu bringen, kann man auch die Befehle Maximize und Minimize verwenden:

Maximize[1.8-0.04 x1-0.03 x2,x1≥0&&x2≥0&&20-x1-x2≥ 0&&
x1+x2≥ 10&& x2≤ 2x1,{x1,x2}]

Dieser Befehl kann auch für nichtlineare Optimierungsprobleme verwendet werden. Allerdings ist dann meist nur noch eine numerische Lösung möglich.

Skalarprodukt und Orthogonalität

Skalarprodukt

Das Skalarprodukt wird mit einem Punkt berechnet:

a={1,2,3};b={2,-4,1};a.b

Die Länge eines Vektors können wir also auch mit

(Sqrt[a.a])

erhalten.

Kreuzprodukt

Das Kreuzprodukt kann mithilfe des Befehls Cross[a, b] berechnet werden:

Cross[{1,2,0}, {3,4,5}]

Ray Tracing

Der reflektiere Strahl erhält man durch lösen folgender Gleichung:

solr=Solve[(((e.n + kr ))^2)==(((e .n))^2)(1+2kr (e.n)+kr^2),kr]

Einsetzten der Lösung ergibt

er=Simplify[e + kr n /. solr[[2]] ]

und konkret für die Beispielvektoren:

% /. {n->{0,0,1}, e->{1,0,-1}/(Sqrt[2])} //Simplify

Analog erhält man den gebrochenen Strahl:

solb=Solve[1 - (((e.n+kb )/Sqrt[(1+2(e.n)kb  +kb ^2)]))^2==nb^2 (1-(e.n)^2),kb]//Simplify

kb =kb/.solb[[1]]

1+2kb (e.n)+kb^2 //Simplify

Clear[solr,solb,er,kb]

Gram-Schmidt-Verfahren

Das Gram-Schmidt-Verfahren kann mit folgendem Befehl durchgeführt werden :

Orthogonalize[{{1,0,1},{0,1,1},{1,2,4}}]

Diskrete Kosinustransformation

Die diskrete Kosinustransformation können wir wie folgt berechnen:

DCT[n_]:=Table[(Sqrt[((2-If[k==1,1,0])/n)])Cos[(((2j-1)(k-1)π)/(2n))],{j,n},{k,n}]

DCT[2]//MatrixForm //Simplify

DCT[3]//MatrixForm //Simplify

Eigenwerte und Eigenvektoren

Koordinatentransformationen

Die Drehmatrix um die z-Achse ist:

U[ϕ_]=(({{Cos[ϕ], -Sin[ϕ], 0}, {Sin[ϕ], Cos[ϕ], 0}, {0, 0, 1}}));

Konkret für eine Drehung um π/4

U[π/4] //MatrixForm

Ein gedrehter Vektor

U[((-π)/4)].{1,3,4} //Simplify//MatrixForm

und eine gedrehte Matrix:

A=(({{0, 1, 0}, {1, 0, 0}, {0, 0, 1}}));

Inverse[U[(π/4)]].A.U[(π/4)]//MatrixForm

Clear[A,U]

Eigenwerte und Eigenvektoren

In Mathematica können die Befehle Eigenvalues und Eigenvectors zur Berechnung von Eigenwerten und Eigenvektoren verwendet werden:

A=(({{0, -1}, {1, 0}}));Eigenvalues[A]

Eigenvectors[A]

Clear[A]

Implizite Kurven

Implizit gegebene Kurven, wie die Ellipse 3 (Subscript[x, 1]^2) + 2 (Subscript[x, 1])  (Subscript[x, 2])+ 3 (Subscript[x, 2]^2) = 4, kann man mit folgendem Befehl zeichnen:

ContourPlot[{x==y,x==-y,3 x^2-2 x y +3 y^2==4},{x,-1.5,1.5},{y,-1.5,1.5}]

Grundlagen der Graphentheorie

Darstellung von Graphen

Graphen können mit dem Befehl AdjacencyGraph[A] eingegeben werden. Dabei ist A die Adjazenzmatrix und P eine Liste der gewünschten Koordinaten der Knoten im (ℝ^2).

A={{0,1,0,1,1}, {1,0,1,1,0}, {0,1,0,1,0}, {1,1,1,0,0},{1,0,0,0,0}};
gr=AdjacencyGraph[A]

Mathematica wählt dabei die Koordinaten der Knoten automatisch. Man kann sie aber auch über die Option VertexCoordinates angeben. Die Knoten können mit der Option VertexLabels bezeichnet werden und die Option ImagePadding stellt sicher, dass die Namen am Rand des Bildes nicht abgeschnitten werden.

AdjacencyGraph[A,VertexCoordinates ->{{0,1}, {1,1}, {2,1}, {1,0},{0,0}},VertexLabels->"Name",ImagePadding->10]

Zusammenhang

Mit dem Befehl

ConnectedGraphQ[gr]

können wir auf Zusammenhang testen.

EulerZug

In Mathematica könnten wir das Beispiel aus dem Buch wie folgt lösen:

gr1=AdjacencyGraph[{ {0,1,0,0,0,0,0,0,0,0,1},{1,0,1,0,0,0,0,0,0,0,0},{0,1,0,1,1,0,0,0,0,0,1},{0,0,1,0,1,0,0,0,0,0,0},{0,0,1,1,0,1,0,1,0,0,0},{0,0,0,0,1,0,1,0,0,0,0},{0,0,0,0,0,1,0,1,0,0,0},{0,0,0,0,1,0,1,0,1,0,1},{0,0,0,0,0,0,0,1,0,1,0},{0,0,0,0,0,0,0,0,1,0,1},{1,0,1,0,0,0,0,1,0,1,0}}, VertexCoordinates ->{ {0,2},{0,1},{1,1},{1.5,0},{2,1},{3,1},{3,2},{2,2},{2,3},{1,3},{1,2}},VertexLabels->"Name",ImagePadding->10]

FindEulerianCycle[gr1]

Mathematica gibt uns einen Euler-Zug mit Start im Knoten 1 (= a) aus. (Die Knoten sind entsprechend der Adjazenzmatrix, über die der Graph in Mathematica eingegeben wurde, durchnummeriert.)

Wir können den Euler-Zug auch hervorheben und animieren:

Animate[HighlightGraph[gr1,Part[First[%],1;;i]],{i,1,Length[First[%]],1}]

Hamilton Kreis

Einen Hamilton-Kreis erhalten wir mit

FindHamiltonianCycle[gr1]

Wiederum können wir die Lösung veranschaulichen:

Animate[HighlightGraph[gr1,Part[First[%],1;;i]],{i,1,Length[First[%]],1}]

Bäume und kürzeste Wege

Gewichtete Graphen

Gewichtete Graphen können analog wie ungewichtete Graphen mit dem Befehl WeightedAdjacencyGraph eingegeben werden. Dabei müssen nicht-vorhandene Kanten durch ein Gewicht ∞ (Eingabe als Infinite oder über die Palette) gekennzeichnet werden.

A=(({{∞, 6, ∞, ∞, ∞, 10, 1}, {6, ∞, 5, ∞, ∞, ∞, 7}, {∞, 5, ∞, 6, ∞, ∞, 9}, {∞, ∞, 6, ∞, 3, ∞, 4}, {∞, ∞, ∞, 3, ∞, 3, 2}, {10, ∞, ∞, ∞, 3, ∞, 10}, {1, 7, 9, 4, 2, 10, ∞}}));gr2=WeightedAdjacencyGraph[A,VertexLabels->"Name",ImagePadding->5,EdgeLabels->"EdgeWeight"]

Mit der Option EdgeLabels wird bei der Darstellung das Gewicht neben jede Kante geschrieben.

Travelling Salesman Problem

Eine billigste Rundreise in einem Graphen erhalten wir mit

ts=FindShortestTour[Range[Length[A]],DistanceFunction->(A[[#1,#2]]&)]

Die Ausgabe enthält die Länge zusammen mit der Rundreise in Form einer Liste der durchlaufenen Knoten. Um das Ergebnis zu veranschaulichen, können wir den gefundenen Weg mit folgendem Befehl hervorheben:

HighlightGraph[gr2,PathGraph[ts[[2]]]]

Der Befehl PathGraph macht aus der Knotenliste den zugehörigen Graphen, der diese Knoten verbindet.

Den vollständige Graphen (Subscript[K, n]) bekommen Sie übrigens mit CompleteGraph[n].

CompleteGraph[3]

Clear[A,ts]

Minimal aufspannende Bäume

Mathematica implementiert keinen Algorithmus für minimal aufspannende Bäume, aber das können wir leicht selbst beheben:

Kruskal[g_Graph]:=Module[{el,i=2,n=VertexCount[g],mst,cmst,mstn,cmstn},
el=Sort[EdgeList[g],PropertyValue[{g,#1},EdgeWeight]<PropertyValue[{g,#2},EdgeWeight]&];
mst=Graph[VertexList[g],{el[[1]]}];cmst=n-1;
While[i≤ n&&cmst≠1,
mstn=GraphUnion[mst,Graph[{el[[i]]}]];
cmstn=Length[ConnectedComponents[mstn]];
If[cmstn<cmst,mst=mstn;cmst=cmstn];
i++];
mst
]

Zum Beispiel

mst=Kruskal[gr2];
HighlightGraph[gr2,mst]

Clear[mst]

Kürzeste Wege

Mit Mathematica kann der kürzeste Weg zwischen zwei Knoten i, j mit dem Befehl ShortestPath[graph, i, j] berechnet werden.

A=(({{∞, 4, 6, 7}, {4, ∞, ∞, 1}, {6, ∞, ∞, 2}, {7, 1, 2, ∞}}));gr3=WeightedAdjacencyGraph[A,VertexLabels->"Name",ImagePadding->5]

FindShortestPath[gr,1,4]

Das Ergebnis können wir natürlich wieder mit

HighlightGraph[gr3,PathGraph[%]]

veranschaulichen.

Clear[A,gr3]

Flüsse in Netzwerken und Matchings

Netzwerke

Definieren wir zunächst einen Graphen

A={{∞,10,∞,12,∞,∞},{∞,∞,9,∞,∞,∞},{∞,∞,∞,4,∞,8},{∞,∞,∞,∞,12,∞},{∞,∞,6,∞,∞,14},{∞,∞,∞,∞,∞,∞}};v={{-1,0},{-0.5,0.86},{0.5,0.86},{-0.5,-0.86},{0.5,-0.86},{1,0}};
gr=WeightedAdjacencyGraph[A,VertexCoordinates->v,ImagePadding->5,GraphStyle->"SmallNetwork",EdgeLabels->"EdgeWeight"]

Nun kann der maximale Fluss mittels

mf=FindMaximumFlow[gr,1,6,"OptimumFlowData",EdgeCapacity->EdgeWeight];

berechnet werden. Mit

mf["FlowValue"]

erhalten wir den maximalen Fluss und mit mf[" FlowGraph"] den zugehörigen Graphen. In diesem können wir die Kantenbezeichnungen auf die Angabe von “Kapazität, Fluss“ erweitern:

SetProperty[mf["FlowGraph"],{EdgeLabels->(#->Row[{PropertyValue[{gr,#},EdgeWeight],",",mf[#]}]&/@EdgeList[gr])}]

Matchings

Definieren wir zunächst einen Graphen

gr=Graph[{1->5,1->6,2->6,2->7,2->8,3->7,4->6},VertexCoordinates->{{0,3},{2,3},{2,1},{0,2},{2,2},{2,0},{0,1},{0,0}},VertexLabels->"Name",ImagePadding->10];

Hier haben wir zur Abwechslung den Graphen als eine Liste von gerichteten Kanten angegeben (die Knotenliste wird automatisch aus den Kanten generiert). Nun kann mit dem Befehl FindIndependentEdgeSet  ein maximales Matching gefunden werden:

FindIndependentEdgeSet[gr]

Veranschaulichen können wir es mit

HighlightGraph[gr,Graph[%/. {a_,b_}-> (a->b)]]

Clear[gr];