HAT

Ein zweischaliges Hyperboloid. Die gefärbten Flächen sind eine Hyperfläche zweiter Ordnung (im dreidimensionalen Raum).

Die Hauptachsentransformation (HAT) ist ein Verfahren aus der linearen Algebra, um Gleichungen für sogenannte Hyperflächen zweiter Ordnung in einer Normalform darstellen und so klassifizieren zu können. In der theoretischen Physik werden Hauptachsentransformationen hauptsächlich zur Diagonalisierung einer Matrix verwendet, wie im Abschnitt „Anwendung“ näher beschrieben.

Als Hyperfläche werden allgemein geometrische Objekte bezeichnet, deren Dimension um genau eins geringer ist als die Dimension des Raumes, in dem sie dargestellt werden. So ist eine Ebene, deren mathematische Dimension 2 ist, im dreidimensionalen („euklidischen“) Raum unserer Anschauung eine Hyperfläche, ebenso wie eine Kurve mit Dimension 1 in einem zweidimensionalen Raum (beispielsweise auf einem Blatt Papier). Eine Hyperfläche zweiter Ordnung (oder zweiten Grades) ist eine Hyperfläche, in deren Bestimmungsgleichung die Variablen höchstens in der zweiten Potenz auftreten. Im zweidimensionalen Raum entsprechen die Hyperflächen zweiter Ordnung den Kegelschnitten.

Vorgehen

Hyperfläche 2. Ordnung im zweidimensionalen Raum, vor Drehung und Verschiebung
(f(x,y):x² − 4xy + y² + 10x + y + 12 = 0)

Die allgemeine Gleichung einer Hyperfläche lautet in euklidischen Räumen:

im $\R^2 \quad \alpha x^2 + \beta xy + \gamma y^2 + \delta x + \varepsilon y + \mu = 0$,

im $\R^3 \quad \alpha x^2 + \beta y^2 + \gamma z^2 + \delta xy + \varepsilon xz + \zeta yz + \eta x + \theta y + \iota z + \mu = 0$

sowie im $\R^n \quad \alpha_1 x_1^2 + \dots + \alpha_n x_n^2 + \beta_1 x_1 x_2 + \beta_2 x_1 x_3 + \dots + \gamma_1 x_1 + \dots + \gamma_n x_n + \mu = 0$

Theoretisch können Hauptachsentransformationen auch in unitären Räumen mit hermiteschen Matrizen über dem Körper der komplexen Zahlen durchgeführt werden, doch zur Vereinfachung und aus Gründen der Anschaulichkeit wird im Folgenden auf die analogen Vorgänge nicht eingegangen werden.

Um die Hyperflächen, die durch die obigen Gleichungen dargestellt werden, alleine anhand ihrer Bestimmungsgleichungen klassifizieren zu können (eine graphische Darstellung ist nur bis in die dritte Dimension sinnvoll möglich), müssen die Gleichungen in eine der drei folgenden Normalformen überführt werden:

$\frac {z_1^2} {\alpha_1^2} \pm \frac {z_2^2} {\alpha_2^2} \pm \dots \pm \frac {z_m^2} {\alpha_m^2} + \beta = 0 \qquad (1)$

$\frac {z_1^2} {\alpha_1^2} \pm \frac {z_2^2} {\alpha_2^2} \pm \dots \pm \frac {z_{r}^2} {\alpha_{r}^2} + \beta = 0 \qquad (2)$

$\frac {z_1^2} {\alpha_1^2} \pm \frac {z_2^2} {\alpha_2^2} \pm \dots \pm \frac {z_r^2} {\alpha_r^2} + \beta_{r+1}z_{r+1} + \dots + \beta_m z_m = 0 \qquad (3)$, jeweils mit $r < m \le n; \ r,m,n \in \mathbb N$

z_i sind dabei die Koordinaten des gewählten Koordinatensystems, β_i allgemeine reelle Koeffizienten und α_i reelle Koeffizienten ungleich Null. n ist die Dimension des Raums, in dem die Fläche dargestellt wird.

Die Umrechnung in Normalform erfolgt durch Verschieben und Drehen der Koordinatenachsen, also eine Transformation des Koordinatensystems.

Zur Vereinfachung der allgemeinen Gleichungen werden diese häufig zunächst mit Hilfe der symmetrischen Matrix $A \in \R^{n \times n}$ in eine Matrizengleichung umgeschrieben:

$x^T A x + b^T x + \mu = 0\$ mit $\mu \in \R, \ b = (b_1, \dots , b_n)^T, \ x = (x_1, \dots , x_n )^T$

Ausmultipliziert in Summenschreibweise lautet die vorhergehende Gleichung:

$\sum_{i,j=1 \atop i=j}^n a_{ij}x_i^2 + 2 \sum_{i=1}^{n-1} \sum_{j=i+1}^n a_{ij}x_i x_j + \sum_{i=1}^n b_i x_i + \mu = 0 \quad (A=(a_{ij}))\,$

Je nach den Werten der Elemente a_ij der Matrix A treten die Koordinaten also in quadratischer, einfacher sowie in untereinander gemischter Form auf. Da A per definitionem eine symmetrische Matrix ist, gilt: a_ij = a_ji, jedes Koordinatenpaar tritt in der Summe der gemischten Terme also nur einmal auf.

Anhand der Elemente von A können an dieser Stelle bereits gewisse Aussagen über das Ergebnis der Klassifizierung getroffen werden: Sind diese so beschaffen, dass A nicht den vollen Rang n hat, so kann die Gleichung nur auf die Normalformen (2) oder (3) zurückgeführt werden, ist rg(A) jedoch gleich n, so ist die vorliegende Gleichung eine Variation von Gleichung (1).

Drehung

Graph von f(x,y) nach Elimination der gemischten Terme
(g(x,y):3x² − y² − 9x + 11y + 12 = 0)

Zunächst muss die Drehung des Objekts zum gewöhnlichen (Ziel-)Koordinatensystem, die durch die gemischten Terme a_ijx_ix_j beschrieben wird, eliminiert werden, also alle a_ij für $i\ne j$ gleich 0 sein. Die Matrix A muss dazu mit einer geeigneten Transformation in eine Diagonalmatrix umgewandelt werden, das heißt, alle Werte, die nicht auf der Hauptdiagonale der neuen Matrix liegen, müssen gleich 0 sein.

Da A eine symmetrische Matrix ist, existiert ein voller Satz von Eigenvektoren und damit eine orthogonale Matrix Q, deren Spalten die normierten Eigenvektoren von A sind. Es gilt mit Q ^{− 1} = Q^T aus dem Spektralsatz:

$Q^TAQ=D \Leftrightarrow A=QDQ^T$

Die Matrix $D=\operatorname{diag}(\lambda_i)_{1 \leq i \leq n}$ ist dabei eine Diagonalmatrix mit den Eigenwerten von A auf der Hauptdiagonale.

Zunächst müssen also die Eigenwerte λ_i sowie die zugehörigen Eigenvektoren v_i von A bestimmt und aus diesen die Matrizen D und Q zusammengesetzt werden. Wenn in der Bestimmungsgleichung keine einfachen Potenzen von x_i existieren, also die letzte Summe in der ausmultiplizierten Matrixgleichung (und damit die Verschiebung des Objekts) gleich 0 ist, wird die Matrix Q nicht benötigt und die Berechnung der Eigenvektoren kann ausgelassen werden.

Setzt man nun in obiger Matrixgleichung A = QDQ^T ein und definiert y: = Q^Tx und c^T: = b^TQ, erhält man

$x^TQD Q^T x + b^TQ Q^T x + \mu = y^TD y + c^T y + \mu = 0 \ (QQ^T = 1)$

Damit kommen die (neuen) Koordinaten y_i nur noch in einfacher und quadratischer, aber nicht mehr in gemischter Form vor, die Hyperfläche wurde gedreht.

Verschiebung

Graph von f(x,y) nach Drehung und Parallelverschiebung
(h(x,y): − 3 / 25x² + 1 / 25y² − 1 = 0)

Nun muss noch eine mögliche Parallelverschiebung des Objekts zum Ursprung des Koordinatensystems beseitigt werden. Dazu wird zunächst die Gleichung der gedrehten Hyperfläche ausmultipliziert:

$y^TD y + c^T y + \mu = \sum_{i=1}^n \lambda_i y_i^2 + \sum_{i=1}^n c_i y_i + \mu = 0$

Sind einer oder mehrere Eigenwerte λ_i gleich null, fallen die entsprechenden Terme von $y_i^2$ aus der Gleichung heraus. Im Folgenden ist also $\lambda \neq 0$.

Anschließend gruppiert man nach den Basisvektoren y_i und bringt die einzelnen Klammern durch affine Substitution (z. B. quadratische Ergänzung) in quadratische Form. Die Klammern sind die Basisvektoren z_i des neuen, verschobenen Koordinatensystems:

$\begin{align} \sum_{i=1}^n \lambda_i(y_i^2 + \frac {c_i}{\lambda_i}y_i) + \mu &= \sum_{i=1}^n \lambda_i(y_i+\frac{c_i}{2 \lambda_i})^2 + \nu \\ &= \sum_{i=1}^n \operatorname{sgn}(\lambda_i) \frac{z_i^2}{\left (\frac{1}{\sqrt{|\lambda_i|}}\right )^2} + \nu = 0 \quad (\forall \lambda _i \neq 0) \end{align}$

wobei wegen der quadratischen Ergänzung $\nu = \mu - \sum_{i=1}^n \frac{c_i^2}{4 \lambda_i}$ ist.

$\operatorname{sgn}(x)$ ist die sogenannte Signumfunktion, die das Vorzeichen von x ausgibt.

Für den Fall, dass der Summand ν ungleich 0 ist, wird dieser nun auf 1 normiert, das heißt, die ganze Gleichung wird durch den Betrag von ν geteilt. Diese ist nun – je nach Rang der Matrix A – auf eine der drei oben genannten Normalformen zurückgeführt.

$\left(\sum_{i=1}^n \operatorname{sgn}(\lambda_i) \frac{z_i^2}{\left (\frac{1}{\sqrt{|\lambda_i|}}\right )^2}\right ) \cdot \frac{1}{|\nu |} \pm 1 = 0, \qquad \nu \neq 0$

Klassifizierung

Anhand der Normalform kann nun eine Klassifizierung der Hyperfläche vorgenommen werden.

Im R²

Im zweidimensionalen Raum gibt es zehn verschiedene Hyperflächen zweiter Ordnung (beziehungsweise Kegelschnitte). Zwei von ihnen sind leere Mengen, die zugehörigen Gleichungen können im Bereich der reellen Zahlen nicht gelöst werden. Eine Lösung in den komplexen Zahlen ist möglich, aber nicht bildlich darstellbar.

Normalform	Kurventyp	rg(A)	Bild		Normalform	Kurventyp	rg(A)	Bild
$\frac{x^2}{\alpha^2}+\frac{y^2}{\beta^2}-1=0$	Ellipse (mit Halbachsen α und β)	2			$x^2-\alpha^2=0\,$	zwei parallele Geraden $(x=\pm \alpha)\,$	1
$\frac{x^2}{\alpha^2}-\frac{y^2}{\beta^2}-1=0$	Hyperbel (mit Halbachsen α und β)	2			$x^2+\alpha^2=0\,$	leere Menge	1
$\frac{x^2}{\alpha^2}+\frac{y^2}{\beta^2}+1=0$	leere Menge	2			$x^2-\delta y=0\,$	Parabel	1		Punkt (0,0)	2			$x^2=0\,$	Doppelgerade (y-Achse)	1
$\frac{x^2}{\alpha^2}-\frac{y^2}{\beta^2}=0$	zwei sich schneidende Geraden $\left( y=\pm \frac {\beta}{\alpha}x \right)$	2			$\kappa x + \lambda y + \nu = 0 \,$	allgemeine Gerade	0
$\alpha > 0; \beta > 0; \delta \neq 0; \kappa, \lambda, \nu \in \R$

Im R³

Im dreidimensionalen Raum gibt es 18 Hyperflächen zweiten Grades, wobei den Flächen mit rg(A)<3 deutlich die „Verwandtschaft“ mit den entsprechenden Flächen im zweidimensionalen Raum angesehen werden kann. Hier gibt es drei Formen der Gleichung, die nur in den komplexen Zahlen gelöst werden können, darunter die beiden Gleichungen, die bereits im R² leere Mengen darstellen.

Normalform	Flächentyp	rg(A)	Bild		Normalform	Flächentyp	rg(A)	Bild
$\frac{x^2}{\alpha^2}+\frac{y^2}{\beta^2}+\frac{z^2}{\gamma^2}-1=0$	Ellipsoid (mit Halbachsen α, β und γ)	3			$\frac{x^2}{\alpha^2}-\frac{y^2}{\beta^2}-1=0$	hyperbolischer Zylinder	2
$\frac{x^2}{\alpha^2}+\frac{y^2}{\beta^2}+\frac{z^2}{\gamma^2}+1=0$	leere Menge	3			$\frac{x^2}{\alpha^2}+\frac{y^2}{\beta^2}+1=0$	leere Menge	2
$\frac{x^2}{\alpha^2}+\frac{y^2}{\beta^2}-\frac{z^2}{\gamma^2}-1=0$	einschaliges Hyperboloid	3			$\frac{x^2}{\alpha^2}+\frac{y^2}{\beta^2}=0$	Gerade (z-Achse)	2
$\frac{x^2}{\alpha^2}+\frac{y^2}{\beta^2}-\frac{z^2}{\gamma^2}+1=0$	zweischaliges Hyperboloid	3			$\frac{x^2}{\alpha^2}-\frac{y^2}{\beta^2}=0$	Ebenenpaar mit Schnittgerade (z-Achse)	2
$\frac{x^2}{\alpha^2}+\frac{y^2}{\beta^2}+\frac{z^2}{\gamma^2}=0$	Punkt (0,0,0)	3			x2 − α2 = 0	paralleles Ebenenpaar	1
$\frac{x^2}{\alpha^2}+\frac{y^2}{\beta^2}-\frac{z^2}{\gamma^2}=0$	elliptischer Kegel	3			x2 + α2 = 0	leere Menge	1
$\frac{x^2}{\alpha^2}+\frac{y^2}{\beta^2}-\delta z=0$	elliptisches Paraboloid	2			x2 − δz = 0	parabolischer Zylinder	1
$\frac{x^2}{\alpha^2}-\frac{y^2}{\beta^2}-\delta z=0$	hyperbolisches Paraboloid	2			x2 = 0	Ebene (y-z-Ebene)	1
$\frac{x^2}{\alpha^2}+\frac{y^2}{\beta^2}-1=0$	elliptischer Zylinder	2			κx + λy + μz + ν = 0	allgemeine Ebene	0
$\alpha > 0; \beta > 0; \gamma > 0; \delta \neq 0; \kappa, \lambda, \mu ,\nu \in \R$

Siehe Hauptachsentransformation#Code_zur_Visualisierung zur Visualisierung.

Anwendung

Neben der rein mathematisch-geometrischen Bedeutung der Hauptachsentransformation zur Bestimmung des Typs von Hyperflächen zweiter Ordnung wird sie in zahlreichen Disziplinen der theoretischen Physik sowie in der Informatik und den Geowissenschaften eingesetzt.

In der theoretischen Physik in der klassischen Mechanik wird die Hauptachsentransformation zur Beschreibung der Kinematik starrer Körper verwendet: Hier können über eine Hauptachsentransformation des Trägheitstensors, der die Trägheiten des Körpers bezüglich Drehungen um verschiedene Achsen angibt, eventuell – zum Beispiel bei einem Kreisel – vorhandene Deviationsmomente zum Verschwinden gebracht werden.

Ein Deviationsmoment ist ein Maß für das Bestreben eines starren Körpers, seine Drehachse zu verändern. Deviationsmomente werden mit den Trägheitsmomenten in Trägheitstensoren zusammengefasst, wobei die Trägheitsmomente sich auf der Hauptdiagonalen des Tensors, die Deviationsmomente auf den Nebendiagonalen befinden. Wie oben gezeigt, kann der symmetrische Trägheitstensor auf eine Diagonalform gebracht werden. Die durch die Hauptachsentransformation festgelegten Achsen des neuen, angepassten Koordinatensystems bezeichnet man als Hauptträgheitsachsen, das neue Koordinatensystem als Hauptachsensystem. Die Diagonalelemente des transformierten Tensors werden konsequent Hauptträgheitsmomente genannt.

Anschaulich fallen nach einer Hauptachsentransformation bei symmetrisch gebauten Körpern die Koordinatenachsen des Hauptachsensystems mit dem Symmetrieachsen des starren Körpers zusammen, gemischte Terme verschwinden also aus der Gleichung.

Auch in weiteren Teilgebieten der klassischen Mechanik wird die Hauptachsentransformation eingesetzt, so zum Beispiel in der Festigkeitslehre zur Berechnung der Hauptspannungen, die auf einen Körper einwirken. Häufig angewandt werden Hauptachsentransformationen weiterhin in der relativistischen Mechanik zur Basisdarstellung der Raumzeit im vierdimensionalen Minkowski-Raum.

Außerdem ist die Hauptachsentransformation ein Teil der Karhunen-Loève-Transformation, die vor allem in der Bildverarbeitung auch als Hauptkomponentenanalyse bezeichnet wird. Manchmal werden die Begriffe synonym gebraucht, doch sind beide Transformationen nicht identisch^[1].

Praktisch wird die Hauptachsentransformation als Teil der Hauptkomponentenanalyse dazu verwendet, die Größe umfangreicher Datensätze ohne wesentlichen Datenverlust zu vermindern. Dabei werden vorhandene Beziehungen zwischen einzelnen statistischen Variablen durch Überführung in ein neues, linear unabhängiges problemangepasstes Koordinatensystem so weit wie möglich reduziert. Beispielsweise kann die Anzahl der benötigten Signalkanäle verringert werden, indem diese nach Varianz geordnet und die Kanäle geringster Varianz gegebenenfalls ohne relevanten Datenverlust aus dem Datensatz entfernt werden. Dadurch können Effizienz und Ergebnis einer späteren Analyse der Daten verbessert werden^[2].

In der elektronischen Bildverarbeitung wird die Reduktion der Datensatzgröße durch Hauptkomponentenanalysen besonders in der Fernerkundung durch Satellitenbilder sowie den zugehörigen naturwissenschaftlichen Disziplinen der Geodäsie, Geographie, Kartographie und Klimatologie eingesetzt. Hier kann die Qualität der Satellitenaufnahmen durch Unterdrückung des Rauschens mittels Hauptkomponentenanalyse deutlich verbessert werden^[3].

In der Informatik wird die Hauptachsentransformation vor allem bei der Mustererkennung und in der neuronalen Netztheorie, also der Schaffung künstlicher neuronaler Netze, einem Teilgebiet der künstlichen Intelligenz, zur Datenreduktion angewandt^[3][4].

Beispiel

f(x,y,z) = x² − 2xy + y² − 2yz + z² − 2xz + 2x + 2y + 2z − 3 = 0

Gegeben sei die Gleichung

$f(x,y,z)=x^2-2xy+y^2-2yz+z^2-2xz+2x+2y+2z-3=0\,$

Mittels Hauptachsentransformation soll diese Gleichung nun in eine Normalform überführt werden und der Typ der durch die Gleichung dargestellten Hyperfläche bestimmt werden.

In Matrizenschreibweise lautet die Gleichung:

$(x,y,z) \begin{pmatrix} 1 &amp;amp;amp; -1 &amp;amp;amp; -1 \\ -1 &amp;amp;amp; 1 &amp;amp;amp; -1 \\ -1 &amp;amp;amp; -1 &amp;amp;amp; 1 \end{pmatrix} \begin{pmatrix} x \\ y \\ z \end{pmatrix} + (2,2,2) \begin{pmatrix} x \\ y \\ z \end{pmatrix} -3 = 0$

Der Rang von A ist gleich 3, das heißt, die Matrix hat vollen Rang und die diskutierte Gleichung ist eine Variation von Gleichung (1). Die Eigenwerte von A sind λ₁ = − 1 mit algebraischer Vielfachheit 1 und λ_2,3 = 2 mit algebraischer Vielfachheit 2; die zugehörigen Eigenvektoren lauten:

$v_1=\begin{pmatrix} 1 \\ 1 \\ 1 \end{pmatrix},\quad v_2=\begin{pmatrix} -1 \\ 1 \\ 0 \end{pmatrix},\quad v_3=\begin{pmatrix} -1 \\ 0 \\ 1 \end{pmatrix}$

Noch stellen v₁, v₂ und v₃ aber keine Orthonormalbasis eines Vektorraums dar, dazu müssen die Vektoren mittels eines geeigneten Verfahrens – beispielsweise des Gram-Schmidtschen Orthogonalisierungsverfahrens – orthogonalisiert und anschließend ihre Beträge auf 1 normiert werden:

$v_1'=\frac{1}{\sqrt{3}}\begin{pmatrix} 1 \\ 1 \\ 1 \end{pmatrix},\quad v_2'=\frac{1}{\sqrt{2}}\begin{pmatrix} -1 \\ 1 \\ 0 \end{pmatrix},\quad v_3'=\frac{1}{\sqrt{6}}\begin{pmatrix} -1 \\ -1 \\ 2 \end{pmatrix}\quad (|v_1'|=|v_2'|=|v_3'|=1)$

Nun werden die Matrizen Q und D aus den Eigenvektoren und Eigenwerten von A zusammengesetzt:

$Q=(v_1',v_2',v_3')=\begin{pmatrix} \frac{1}{\sqrt{3}} &amp;amp;amp; \frac{-1}{\sqrt{2}} &amp;amp;amp; \frac{-1}{\sqrt{6}} \\ \frac{1}{\sqrt{3}} &amp;amp;amp; \frac{1}{\sqrt{2}} &amp;amp;amp; \frac{-1}{\sqrt{6}} \\ \frac{1}{\sqrt{3}} &amp;amp;amp; 0 &amp;amp;amp; \frac{2}{\sqrt{6}} \end{pmatrix},\quad D=\begin{pmatrix} -1 &amp;amp;amp; 0 &amp;amp;amp; 0 \\ 0 &amp;amp;amp; 2 &amp;amp;amp; 0 \\ 0 &amp;amp;amp; 0 &amp;amp;amp; 2 \end{pmatrix}$

Es ergibt sich für die gedrehte Hyperfläche mit Gruppierung nach den Basisvektoren x', y' und z' sowie quadratischer Ergänzung:

$(x',y',z')\begin{pmatrix} -1 &amp;amp;amp; 0 &amp;amp;amp; 0 \\ 0 &amp;amp;amp; 2 &amp;amp;amp; 0 \\ 0 &amp;amp;amp; 0 &amp;amp;amp; 2 \end{pmatrix}\begin{pmatrix} x' \\ y' \\ z' \end{pmatrix} + (2\sqrt{3},0,0) \begin{pmatrix} x' \\ y' \\ z' \end{pmatrix} -3 = 0$

ein elliptischer Kegel (Kreiskegel) mit Gleichung in Normalform
(g(x,y,z): − x² + 2y² + 2z² = 0)

$-{x'}^2+2{y'}^2+2{z'}^2 + 2\sqrt{3}\ x'-3=-1\left({x'}^2-2\sqrt{3}\ x'\right)+2{y'}^2+2{z'}^2-3=$

$=-1\left(x'-\sqrt{3}\right)^2+2{y'}^2+2{z'}^2-3+3=0$

Oder transformiert in Basisvektoren x″, y″ und z″:

$-\frac{{x''}^2}{1^2}+ \frac{{y''}^2}{\left(\frac{1}{\sqrt{2}}\right)^2}+ \frac{{z''}^2}{\left(\frac{1}{\sqrt{2}}\right)^2}=0$

Die Gleichung beschreibt also einen elliptischen Kegel. In diesem speziellen Fall sind die beiden Halbachsen der Kegelgrundfläche (beschrieben durch die Vorfaktoren zu y''² und z''²) gleich groß, es handelt sich also um einen Kreiskegel.

Einzelnachweise

1. ↑ vgl. Skript Mustererkennung, Kap. 5.2: Karhunen-Loeve-Transformation. Laboratorium für Nachrichtenverarbeitung, Universität Hannover (pdf)
2. ↑ Gesellschaft für Datenanalyse und Fernerkundung Hannover
3. ↑ ^{a b} Probabilistische Hauptachsentransformation zur generischen Objekterkennung, Diplomarbeit im Fach Informatik, Friedrich-Alexander-Universitat Erlangen-Nürnberg (Postscript) (Textversion)
4. ↑ Neuronale Netze, Heinrich-Heine-Universität Düsseldorf (pdf)

Code zur Visualisierung

MATLAB

%p = impl3('x^2+y^2/2+z^2/3-1=0');
%p = impl3('x^2+y^2-z^2-1=0',20);
%p = impl3('x^2+y^2-z^2+1=0',20);
%p = impl3('x^2+y^2-z^2=0',20);
%p = impl3('x^2+y^2/2-z=0',20);
%p = impl3('x^2-y^2/2-z=0',20);
%p = impl3('x^2+y^2/2-1=0',30);
%p = impl3('x^2-y^2-1=0',20);
%p = impl3('x^2-y^2=0');
%p = impl3('x^2-1=0',20);
%p = impl3('x^2-z=0',30);
p = impl3('x-y+z+1=0',20);
 
[n,m] = size(get(p, 'Vertice'));
set(p, 'FaceColor', 'flat');
set(p, 'CData', linspace(0,255,n));
 
%print -dpng ellipsoid.png;
%print -dpng einsch_hyp.png;
%print -dpng zweisch_hyp.png;
%print -dpng ell_keg.png;
%print -dpng ell_par.png;
%print -dpng hyp_par.png;
%print -dpng ell_zyl.png;
%print -dpng hyp_zyl.png;
%print -dpng ebenenp_ger.png;
%print -dpng ebenenp_par.png;
%print -dpng par_zyl.png;
print -dpng allg_ebene.png;

%IMPL3   3-D plot of implicit defined functions.
%   IMPL3(FUN_STR) plots the function FUN_STR as patched isosurface.
%
%   IMPL3(FUN_STR, N) plots the function over the default domain using
%   an N-by-N grid. The default value is 60.
%
%   IMPL3(FUN) plots the function FUN(X,Y,Z)=0 as patched isosurface.
%
%   IMPL3(FUN, C) plots the function FUN(X,Y,Z)=C as patched isosurface.
%
%   An implicit function is defined this way:
%      FUN_STR(X,Y,Z) = C, where C is a constant
%      FUN(X,Y,Z) = C
%
%   IMPL3 returns the patch object (handle to surface plot object?).
%
%   Notes:
%   A way to improve speed is to use 'pcode impl3' once.
%   One programming technique is to vectorize the string expression, look
%   at EZSURF for more information. In EZSURF this is optional but here it
%   is needed, thus the operators are converted.
%
%   Example:
%      impl3('X.^2+Y.^2-Z.^2-5=0');
%      impl3('x^2+y^2-z^2-5=0');
%
%   See also quadr, patch, isosurface, isonormals, view, ndgrid, pcode.
%
%   (implemented by [e-mail] - 20080530 - with special
%   thanks to Michael!)
function [p] = impl3(fun_str, varargin)
ip = inputParser;   % Create an instance of the class.
if ischar(fun_str)
    ip.addRequired('fun_str', @ischar); % validated by 'ischar'
    ip.addOptional('n', 60);            % no validator
    ip.parse(fun_str, varargin{:});
    n = ip.Results.n;
 
    fun_str = regexprep(fun_str, 'x', 'X');
    fun_str = regexprep(fun_str, 'y', 'Y');
    fun_str = regexprep(fun_str, 'z', 'Z');
    fun_str = regexprep(fun_str, '\.\*', '\*');
    fun_str = regexprep(fun_str, '\.\/', '\/');
    fun_str = regexprep(fun_str, '\.\^', '\^');
    fun_name = fun_str;
    fun_str = regexprep(fun_str, '\*', '\.\*');
    fun_str = regexprep(fun_str, '\/', '\.\/');
    fun_str = regexprep(fun_str, '\^', '\.\^');
 
    i = strfind(fun_str,'=');
    c = eval(fun_str((i+1):end));
    fun = fun_str(1:(i-1));
    fun = inline(fun, 'X', 'Y', 'Z');
else
    ip.addRequired('fun');              % no validator
    ip.addOptional('c', 0);             % no validator
    ip.parse(fun_str, varargin{:});
    c = ip.Results.c;
 
    fun = fun_str;
    fun_name = func2str(fun);
    n = 60;
end
 
[Y,X,Z] = ndgrid(linspace(-5,5,n),linspace(-5,5,n),linspace(-5,5,n));
V = fun(X,Y,Z);
p = patch(isosurface(X,Y,Z,V,c));
isonormals(X,Y,Z,V,p);
set(p,'FaceColor','b','EdgeColor','k','FaceAlpha',0.5);
grid on;
view(3);
title(fun_name);
 
end

Octave

Eine mit MATLAB vergleichbare Qualität liefert die Verwendung von Octaviz. Durch die Einbindung von VTK ist eine professionelle Darstellung von Oberflächen möglich.

Literatur

• Nolting, W.: Grundkurs Theoretische Physik 1: Klassische Mechanik. 7. Auflage. Springer, Berlin 2004, ISBN 3-540-21474-7
• Fließbach, T.: Mechanik. 2. Auflage. Spektrum, Heidelberg 1996, ISBN 3-86025-686-6
• Greiner, W.: Theoretische Physik Band 2. Mechanik Teil 2. 5. Auflage. Verlag Harri Deutsch, Thun/Frankfurt am Main, 1989, ISBN 3-8171-1136-3