%%% -------------------------------------------------------------------- 
%%% -- Folien zum Octavekurs im Oktober 1997 und 1999                 -- 
%%% -- Letzte Aenderung: 19-10-1999                                   -- 
%%% -- Autor: Hubert Selhofer (hubert@na.uni-tuebingen.de)            -- 
%%% -- Neu formatiert mit kleinen Aenderungen: 2001/10/16-2001/10/29  --
%%% -- Marcel Oliver (oliver@member.ams.org)                          -- 
%%% -------------------------------------------------------------------- 
\documentclass[a4paper]{article}
\usepackage{german}
\usepackage[latin1]{inputenc}
\usepackage{graphicx,amsmath}

\newcommand{\Octave}{{\textsf{Octave}}}
\newcommand{\Matlab}{{\textsf{Matlab}}}

\begin{document}

\title{Einführung in GNU Octave}
\date{2001/11/09}
\author{Hubert Selhofer, mit Änderungen von Marcel Oliver}
\maketitle

\begin{center}
  \includegraphics[height=10cm,angle=270]{sombr.ps}
\end{center}

\tableofcontents


\section{Grundlagen}

\subsection{Was ist \Octave?}

\Octave{} ist eine interaktive Scriptsprache, die speziell für
vektorisierbare Berechnungen optimiert ist und dabei Standardroutinen
der numerischen Mathematik (z.B. EISPACK oder LAPACK) auf einfache
Weise zugänglich macht.

Der Syntax von \Octave{} ist dem proprietären \Matlab{} sehr ähnlich,
d.h.  ein \Octave{} Programm kann meist auch von \Matlab{} ausgeführt
werden.  Die Rückwärtskompatibilität ist nicht immer gegeben, da
\Matlab{} besonders im Grafikbereich um einen erheblich größeren
Funktionsumfang verfügt.

\subsection{Hilfe!}

\begin{itemize}
\item \texttt{help} listet alle Befehle und internen Variablen auf.
\item \texttt{help name} gibt Hilfetext zur Variable oder Funktion
  "`name"' aus.
\end{itemize}

\paragraph*{Beispiel:}
\begin{verbatim}
  octave:1> help eig
\end{verbatim}


\subsection{Eingabekonventionen}

\begin{itemize}
\item Alle Befehle können sowohl interaktiv, als auch über
  Skriptdateien eingegeben werden.
\item Scripts sind Textdateien mit dem Suffix \texttt{.m}.  Sie werden
  durch Aufruf des Dateinamens ohne Suffix eingelesen, und verhalten
  sich so, als ob ihr Inhalt Zeile für Zeile am Prompt eingegeben
  würde.
\item \texttt{;} trennt mehrere Befehle in einer Zeile, und
  unterdrückt die Zusgabe von Werten.
\item \texttt{,} trennt Befehle, gibt aber Werte aus.
\item \texttt{...} am Zeilenende bedeutet, dass der Ausdruck in der
  nächsten Zeile weitergeht.
\item Kommentare werden mit \texttt{\%} eingeleitet.
\item \Octave{} unterscheidet zwischen Groß- und Kleinbuchstaben.
\end{itemize}


\subsection{Variablen und Standardoperationen}

\begin{itemize}
\item \texttt{varname = expression} weist \texttt{varname} den Wert
  \texttt{expression} zu.
\item Es gibt die "`üblichen"' Funktionen \texttt{+}, \texttt{-},
  \texttt{*}, \texttt{/}, \texttt{\^{}}, \texttt{sin}, \texttt{cos},
  \texttt{exp}, \texttt{acos}, \texttt{abs}, usw.
\item Die elementaren logischen Operatoren sind:

\begin{tabular}{clclcl}
  \texttt{<} & kleiner \qquad &
  \texttt{<=} & kleiner oder gleich \qquad &
  \texttt{\&} & und \\
  \texttt{>} & größer &
  \texttt{>=} & größer oder gleich &
  \texttt{|} & oder \\
  \texttt{==} & gleich &
  \texttt{<>} & ungleich &
  \verb+~+ & nicht 
\end{tabular}
\end{itemize}
Zur Diagnose von benutzerdefinierten Objekten sind folgende Befehle
hilfreich: 
\begin{itemize}
\item \texttt{whos} zeigt alle definierten Variablen und Funktionen
  (siehe Abschnitt~\ref{funktionen}) an.
\item \texttt{clear name} löscht \texttt{name} aus dem Speicher; falls
  "`name"' nicht angegeben wird, werden {\it alle} Variablen und
  Funktionen gelöscht.
\item \texttt{type name} gibt die Variable oder Funktion \texttt{name}
  am Bildschirm aus.
\end{itemize}

\paragraph*{Beispiele:}
\begin{verbatim}
  octave:1> x12 = 1/8, long_name = 'Ein String'
  x12 = 0.12500
  long_name = Ein String
  octave:2> sqrt(-1)-i
  ans = 0
  octave:3> x = sqrt(2); sin(x)/x
  ans = 0.69846
\end{verbatim} 
Und hier ein Script \texttt{doless}, der in der Datei
\texttt{doless.m} gespeichert ist:
\begin{verbatim}
  eins = 1;
  zwei = 2;
  drei = eins + zwei;
\end{verbatim}
Aufruf des Scripts:
\begin{verbatim}
  octave:1> doless
  octave:2> whos
  *** local user variables:
  prot  type               rows   cols  name
  ====  ====               ====   ====  ====
   wd   real scalar           1      1  drei
   wd   real scalar           1      1  eins
   wd   real scalar           1      1  zwei
\end{verbatim} 


\section{Vektor- und Matrixoperationen}

Matrizen und Vektoren sind die wichtigsten Grundbausteine zur
Programmierung in \Octave .

\subsection{Vektoren}
\label{vektor}

\begin{itemize}

\item Zeilenvektor $v = (1,2,3)$:

  \texttt{v = [ 1 2 3 ]}

\item Spaltenvektor: $v = (1,2,3)^T$:

  \texttt{v = [ 1; 2; 3 ]}
  
\item Automatische Erzeugung von Vektoren mit konstantem Inkrement:

  \texttt{Anfang[:Inkrement]:Ende}
\end{itemize}

\paragraph*{Beispiele:}
\begin{verbatim}
  octave:1> x = 3:6
  x =
    3  4  5  6
  octave:2> y = 0:.15:.7
  y =
  0.00000  0.15000  0.30000  0.45000  0.60000
  octave:3> z = pi:-pi/4:0
  z =
  3.14159  2.35619  1.57080  0.78540  0.00000
\end{verbatim}


\subsection{Matrizen}

Eine Matrix $\displaystyle A = \begin{pmatrix} 1 & 2 \\ 3 & 4
\end{pmatrix}$ wird erzeugt durch:
\begin{verbatim}
  octave:1> A = [ 1 2; 3 4]
  A =
    1  2
    3  4
\end{verbatim}
Matrizen können folgenderweise aus Teilmatrizen zusammengesetzt
werden:
\begin{verbatim}
  octave:2> b = [5; 6];
  octave:3> M = [A b]
  M =

    1  2  5
    3  4  6
\end{verbatim}

Für einige wichtige $(m \times n)$-Matrizen existieren eigene Befehle.
Falls $m$ gleich $n$ ist, muß nur ein Argument angegeben werden.
\begin{itemize}
\item \texttt{eye(m,n)} erzeugt eine Matrix mit $1$-ern auf der
  Hauptdiagonale. Im Spezialfall $m=n$ ergibt sich die Einheitsmatrix.
\item \texttt{zeros(m,n)} erzeugt eine Nullmatrix der Dimension $(m
  \times n)$.
\item \texttt{ones(m,n)} erzeugt eine Matrix der Dimension $(m \times
  n)$ mit $a_{ij} = 1$ $\forall i,j$.
\item \texttt{rand(m,n)} erzeugt eine Zufallsmatrix der Dimension $(m
  \times n)$ mit gleichverteilten Einträgen. Andere Verteilungen siehe
  \texttt{help randn}.
\end{itemize}


\subsection{Grundlegende Operationen}

\begin{itemize}
\item \texttt{+}, \texttt{-} und \texttt{*} bezeichnen Matrixaddition,
  -subtraktion und -multiplikation
\item \texttt{A'} transponiert und konjugiert \texttt{A}.
\item \texttt{A.'} transponiert \texttt{A}.
\end{itemize}

\paragraph*{Beispiele:}
\begin{verbatim}
  octave:1> A = [1 2; 3 4]; B = 2*ones(2,2);
  octave:2> A+B, A-B, A*B
  ans =
    3  4
    5  6
  ans =
    -1   0
     1   2
  ans =
     6   6
    14  14
\end{verbatim}


\subsection{Elementweise Operationen}
\label{element}

Während z.B. \texttt{*} die standard Matrixmultiplikation bezeichnet,
gibt es alle Operatoren mit vorangestelltem Punkt für elementweise
Verknüpfung: \texttt{.*}, \texttt{./}, \texttt{.\^{}}, u.s.w.

\paragraph*{Beispiele:}
\begin{verbatim}
  octave:1> A = [1 2; 3 4]; A.^2  % elementweise Potenz
  ans =

     1   4
     9  16

  octave:2> A^2  % Zum Vergleich die matrixweise Potenz
  ans =

     7  10
    15  22
\end{verbatim} 


\subsection{Indizierung und Slicing}

\begin{itemize}
\item \texttt{v(k)} ist das $k$-te Element des Zeilen- oder
  Spaltenvektors $v$.
\item \texttt{A(k,l)} ist das Matrixelement $a_{kl}$.
\item \texttt{v(m:n)} ist das ``Slice'' des Vektors $v$ vom $m$-ten
  bis zum $n$-ten Eintrag.
\item \texttt{A(k,:)} ist die $k$-te Zeile der Matrix $A$.
\item \texttt{A(:,l)} ist die $l$-te Spalte der Matrix $A$.
\end{itemize}
Bestimmung der Größe:
\begin{itemize}
\item \texttt{length(v)} gibt die Länge des Vektors $v$ zurück.
\item \texttt{[Zeilen,Spalten] = size(A)} gibt die Anzahl der Zeilen
  und Spalten von $A$ zurück.
\end{itemize}
Änderung der Gestalt:
\begin{itemize}
\item \texttt{reshape(A,m,n)} formt \texttt{A} in eine $(m \times
  n)$-Matrix um.
\item \texttt{diag(A)} liefert die Diagonale $a_{jj}$ der Matrix
  \texttt{A} in einem Spaltenvektor.
\item \texttt{diag(v)} erzeugt eine Matrix mit dem Vektor \texttt{v}
  auf der Hauptdiagonalen.
\item \texttt{diag(v,k)} erzeugt eine Matrix mit dem Vektor \texttt{v}
  in der \texttt{k}-ten Nebendiagonalen.
\item \texttt{A(k,:) = zv} weist der \texttt{k}-ten Zeile von
  \texttt{A} den Zeilenvektor \texttt{zv} zu.
\item
  \texttt{A(k,:) = []} löscht die \texttt{k}-te Zeile aus \texttt{A}.
\item \texttt{A(:,k) = sv} weist der \texttt{k}-ten Spalte von
  \texttt{A} den Spaltenvektor \texttt{sv} zu.
\item \texttt{A(:,k) = []} löscht die \texttt{k}-te Spalte aus
  \texttt{A}.
\end{itemize}

\paragraph*{Beispiele:}
\begin{verbatim}
  octave:1> A = [1 2 3; 4 5 6]; v = [7; 8];
  octave:2> A(2,3) = v(2)
  A =

    1  2  3
    4  5  8

  octave:3> A(:,2) = v
  A =

    1  7  3
    4  8  8

  octave:4> A(1,1:2) = v'
  A =

    7  8  3
    4  8  8
\end{verbatim} 


\subsection{Lösung linearer Gleichungssysteme}

\begin{itemize}
\item \verb+A\b+ löst die Gleichung $Ax = b$.
\end{itemize}

Falls $A$ eine $(n \times n)$-Matrix ist, wird Gauß\-elimination
verwendet, sonst wird via QR-Zerlegung eine Lösung im Sinne der
kleinsten Fehlerquadrate berechnet. Ist $A$ schlecht konditioniert
oder singulär, wird eine Warnung ausgegeben.


\subsection{Inverse, Zerlegungen, Eigenwerte, etc.}

\begin{itemize}
\item \texttt{B = inv(A)} berechnet die Inverse von $A$.
\item \texttt{[L,U,P] = lu(A)} berechnet eine LU-Zerlegung $LU=PA$.
\item \texttt{[Q,R] = qr(A)} berechnet die QR-Zer\-legung $QR=A$.
\item \texttt{R = chol(A)} berechnet die Choleskyzerlegung von $A$.
\item \texttt{S = svd(A)} berechnet die Singulärwerte von $A$.
\item \texttt{H = hess(A)} bringt $A$ in Hessenbergform.
\item \texttt{E = eig(A)} berechnet alle Eigenwerte von $A$.
\item \texttt{[V,D] = eig(A)} berechnet eine Diagonalmatrix $D$, die die
  Eigenwerte enthält und eine Matrix $V$ mit den dazugehörigen
  Eigenvektoren.  Es gilt also $AV=VD$.
\item \texttt{norm(X,p)} berechnet die $p$-Norm des Vektors $X$. Falls
  $X$ eine Matrix ist, kann $p$ die Werte \texttt{1}, \texttt{2} oder
  \texttt{inf} haben.  Der Default ist $p=2$.
\item \texttt{cond(A)} berechnet die Konditionszahl von $A$ in $2$-Norm.
\end{itemize}
Viele dieser Befehle unterstützen weitere Optionen. Diese können Sie
mit \texttt{help funcname} erfahren.


\subsection{Test auf Nulleinträge}

\begin{itemize}
\item \texttt{[i,j,v] = find(A)} findet die Indizes aller Elemente
  $\neq 0$ in \texttt{A} Es gilt: $A_{i(l),j(l)} = v_l \neq 0$.
\item \texttt{any(A)} ist $1$, falls \texttt{A} Elemente $\neq 0$
  enthält. Bei Matrizen operiert \texttt{any} auf den Spalten.
\end{itemize}


\section{Kontrollstrukturen}

\subsection{Funktionen}
\label{funktionen}

Traditionell sind Funktionen ebenfalls Textdateien mit Suffix
\texttt{.m}.  Im Unterschied zu Skripts können Funktionen mit
Argumenten aufgerufen werden, und alle innerhalb der Funk\-tion
verwendeten Variablen sind \emph{lokal} (sie beeinflussen die vorher
definierten Variablen nicht).

\paragraph*{Beispiel:}

Eine Funktion \texttt{f}, die in der Datei \texttt{f.m} gespeichert
ist.
\begin{verbatim}
  function y = f (x)
    y = cos(x/2)+x;
  end 
\end{verbatim} 

\paragraph*{Bemerkung:}

In \Octave\ kann man---sinnvollerweise---beliebig viele Funktionen in
einer Skriptdatei definieren.  \Matlab\ hingegen lässt strikt nur eine
Funktion pro \texttt{.m} Datei zu, und der Funktionsname muss mit dem
Dateinamen übereinstimmen.  Falls Kompatibilität mit \Matlab\ wichtig
ist, sollte diese Einschränkung auch auf \Octave-Programme angewandt
werden.

\paragraph*{Beispiel mit zwei Funktionswerten:}

Eine Funktion \texttt{dolittle}, die in der Datei \texttt{dolittle.m}
gespeichert ist.
\begin{verbatim}
  function [out1,out2] = dolittle (x)
    out1 = x^2;
    out2 = out1*x;
  end 
\end{verbatim} 
Aufruf der Funktion:
\begin{verbatim}
  octave:1> [x1,x2]=dolittle(2)
  x1 = 4
  x2 = 8
  octave:2> whos
  *** currently compiled functions:
  prot  type               rows   cols  name
  ====  ====               ====   ====  ====
   wd   user function         -      -  dolittle
  *** local user variables:
  prot  type               rows   cols  name
  ====  ====               ====   ====  ====
   wd   real scalar           1      1  x1
   wd   real scalar           1      1  x2
\end{verbatim} 
Offensichtlich waren die beiden Variablen \texttt{out1, out2} in
\texttt{dolittle} lokal.  Eine vorher definierte Variable
\texttt{out1} oder \texttt{out2} wird durch den Aufruf der Funktion
nicht beeinflußt.


\subsection{Globale Variablen}

\texttt{global name} deklariert \texttt{name} als globale Variable.

\paragraph*{Beispiel:}

Eine Funktion \texttt{foo} in der Datei \texttt{foo.m}:
\begin{verbatim}
  function out = foo(arg1,arg2)
  global N % verwende N als globale Variable
  <Berechnung>
  end
\end{verbatim} 
Wird $N$ in der Funktion verändert, ändert das auch den ursprünglichen
Wert von $N$ im "`workplace"'.


\subsection{Schleifen}

Die Syntax von \emph{for}- und \emph{while}-Schleifen wird an
folgenden Beipielen klar:
\begin{verbatim}
  for n = 1:10
    [x(n),y(n)]=dolittle(n);
  end

  while t<T
    t = t+h;
  end
\end{verbatim}
\emph{For}-Schleife rückwärts:
\begin{verbatim}
  for n = 10:-1:1
  ...
\end{verbatim}

\subsection{Verzweigungen}

Bedingte Verzweigungen lassen sich so programmieren:
\begin{verbatim}
  if x==0
    error('x ist 0!');
  else
    y = 1/x;
  end

  switch pnorm
    case 1; 
      sum(abs(v))
    case inf; 
      max(abs(v))
    otherwise 
      sqrt(v'*v)
  end
\end{verbatim} 


\subsection{Funktionen von Funktionen}

\begin{itemize}
\item \texttt{eval(string)} wertet \texttt{string} als \Octave-Code
  aus.
\item \texttt{feval(funcname,arg1,$\dots$)} ist äquivalent zum Aufruf
  von \texttt{funcname} mit Argumenten \texttt{arg1,$\dots$}.
\end{itemize}

\paragraph*{Beispiel}

Berechne mit der Mittelpunktsregel:
\[
  \int_a^b f(x) \, dx \approx \frac{b-a}{N} \sum_{j=0}^{N-1} 
  f \left(a + (j+\frac{1}{2}) \frac{b-a}{N}\right) 
\]
Die Funktionen \texttt{gauss.m} und \texttt{mpr.m} seien
folgendermaßen definiert:
\begin{verbatim}
  function y = gauss(x)
    y = exp(-x.^2/2); 
  end

  function S = mpr(fun,a,b,N)
    h = (b-a)/N;
    S = h*sum(feval(fun,[a+h/2:h:b]));
  end
\end{verbatim}
Jetzt kann man die Funktion \texttt{gauss} durch folgenden Aufruf
integrieren:
\begin{verbatim}
octave:1> mpr('gauss',0,5,500)
\end{verbatim}


\subsection{Effizienzüberlegungen}

Schleifen und Funktionsaufrufe, insbesondere über den Umweg
\texttt{feval}, sind sehr teuer. Deswegen möglichst alle Operationen
vektorisieren.

\paragraph*{Beispiel:}

Wir programmieren die Mittelpunktsregel aus dem vorherigen Abschnitt
mit einer \emph{for}-Schleife (Datei \verb+mpr_langsam.m+):
\begin{verbatim}
  function S = mpr_langsam(fun,a,b,N)
    h = (b-a)/N; S = 0;
    for k = 0:(N-1), 
      S = S + feval(fun,a+h*(k+1/2)); 
    end
    S = h*S;
  end
\end{verbatim}
Man vergleiche die Ausführungszeiten:
\begin{verbatim}
  octave:1> t = cputime;
  > Int1=mpr('gauss',0,5,500); t1=cputime-t;
  octave:2> t = cputime;
  > Int2=mpr_langsam('gauss',0,5,500); t2=cputime-t;
  octave:3> Int1-Int2, t2/t1
  ans = 0
  ans = 45.250
\end{verbatim}


\subsection{Ein- und Ausgabe}

\begin{itemize}
\item \texttt{save datei var1 [var2 $\dots$]} speichert die Variablen
  \texttt{var1} usw.\ in der Datei mit dem Namen \texttt{datei} ab.
\item \texttt{load datei} liest die Datei \texttt{datei}.
\item \texttt{fprintf(string[,var1,$\dots$])} lehnt sich stark an den
  Syntax in C an, siehe \texttt{man fprintf} unter Unix.
\item \texttt{format [long|short]} vergrößert oder verkleinert die
  Anzahl der ausgegebenen Kommastellen. Der Aufruf ohne Argument
  stellt das Standardverhalten wieder her.
\end{itemize}

\paragraph*{Beispiel:}
\begin{verbatim}
  octave:1> for k = .1:.2:.5, 
  > fprintf('1/%g = %10.2e\n',k,1/k); end
  1/0.1 =   1.00e+01
  1/0.3 =   3.33e+00
  1/0.5 =   2.00e+00
\end{verbatim} 


\section{Grafik}

\subsection{2D-Grafiken}

\begin{itemize}
\item \texttt{plot(x,y[,fmt])} zeichnet eine Linie durch die Punkte
  $(x_j,y_j)$. Mit dem String \texttt{fmt} können Linienart und Farbe
  gewählt werden (siehe \texttt{help plot}).
\item \texttt{semilogx(x,y[,fmt])} Syntax wie \texttt{plot}, die
  $x$-Achse wird logarithmisch gezeichnet.
\item \texttt{semilogy(x,y[,fmt])} Syntax wie \texttt{plot}, die
  $y$-Achse wird logarithmisch gezeichnet.
\item \texttt{loglog(x,y[,fmt])} Syntax wie \texttt{plot}, beide
  Achsen werden logarithmisch gezeichnet.
\end{itemize}

Soll eine Funktion $y=f(x)$ geplottet werden, geht man
praktischerweise folgendermassen vor:
\begin{enumerate}
\item Man erzeugt einen Vektor mit den $x$-Koordinaten aller
  zu plottenden Punkte:
\begin{verbatim}
  x = x_min:schrittweite:x_max;
\end{verbatim}
  (Siehe auch Abschnitt~\ref{vektor}.)
\item Man erhält einen Vektor mit den dazugehörigen $y$-Werten, indem
  man die funktion $f$ elementweise auf den $x$-Vektor wirken lässt:
\begin{verbatim}
  y = f(x);
\end{verbatim}
  \emph{Wichtig:} Damit $f$ elementweise operiert, müssen die
  Operatoren \texttt{.+}, \texttt{.-}, \texttt{.\^{}} usw.\ anstelle
  der üblichen \texttt{+}, \texttt{-} und \texttt{\^{}} verwendent
  werden!  (Siehe Abschnitt~\ref{element}.)
\item Jetzt kann man den Plot-Befehl aufrufen:
\begin{verbatim}
  plot(x,y); grid; replot
\end{verbatim}
\end{enumerate}

\paragraph*{Beispiel:}

\begin{verbatim}
  octave:1> x = -10:.1:10; 
  octave:2> y = sin(x).*exp(-abs(x));
  octave:3> plot(x,y); grid; replot
\end{verbatim}
\begin{center}
  \includegraphics[height=7cm,angle=270]{2d-plot1.ps}
\end{center}

\subsection{3D-Grafiken:}
\begin{itemize}
\item \texttt{[xx,yy] = meshgrid(x,y)} erzeugt aus den Vektoren $x$
  und $y$ Gitterdaten.
\item \texttt{mesh(x,y,z)} zeichnet eine Fläche in 3D.
\end{itemize}

\paragraph*{Beispiel:}
\begin{verbatim}
  octave:1> x = -2:0.1:2;
  octave:2> [xx,yy] = meshgrid(x,x);
  octave:3> z = sin(xx.^2-yy.^2);
  octave:4> mesh(x,x,z);
\end{verbatim} 
\begin{center}
  \includegraphics[height=8cm,angle=270]{3d-plot1.ps}
\end{center}


\subsection{Befehle für 2D- und 3D-Grafiken:}

\begin{itemize}
\item \texttt{title(string)} schreibt \texttt{string} als Titelzeile
  in die Grafik.
\item \texttt{xlabel(string)} beschriftet die $x$-Achse mit
  \texttt{string}.
\item \texttt{ylabel(string)} beschriftet die $y$-Achse mit
  \texttt{string}.
\item \texttt{zlabel(string)} beschriftet die $z$-Achse mit
  \texttt{string}.
\item \texttt{axis(v)} legt den Ausschnitt fest; \texttt{v} ist ein
  Vektor der Form \texttt{v = (xmin, xmax, ymin, ymax[, zmin zmax])}.
\item \texttt{hold [on|off]} legt fest, ob die nächste Grafikausgabe
  die alte Grafik löschen soll.
\item \texttt{clg} löscht das Grafikfenster.
\item \texttt{replot} erzwingt eine erneute Ausgabe der aktuellen
  Grafik.  Dies ist in \Octave{} bei Änderungen einer bestehenden
  Grafik nötig.  (Dieser Befehl existiert nur in \Octave{} und ist
  auch nur da notwendig.)
\end{itemize}


\section{Übungen}

\subsection{Lineare Algebra}

Erzeugen Sie eine Matrix $A$ und einen Vektor $b$ mit
\[
A = \begin{pmatrix}
  1 & 2 \\
  3 & 4 
\end{pmatrix} \quad \text{und} \quad
b = \begin{pmatrix}
  36 \\
  88
\end{pmatrix}.
\]
Lösen Sie nun das Gleichungssystem $Ax = b$ durch Eingabe von $4$
Zeichen. Berechnen Sie danach die LU- und QR-Zerlegung sowie die
Eigenwerte mit Eigenvektoren von $A$. Rechnen Sie die
Choleskyzerlegung von $A^T A$ aus, und überzeugen Sie sich, dass
$\mbox{cond}(A^T A) = \mbox{cond}(A)^2$ ist.

\paragraph*{Lösungsvorschlag:}

\begin{verbatim}
  A = reshape(1:4,2,2).'; b = [36; 88];
  A\b
  [L,U,P] = lu(A)
  [Q,R] = qr(A)
  [V,D] = eig(A)
  A2 = A.'*A;
  R = chol(A2)
  cond(A)^2 - cond(A2)
\end{verbatim} 

\subsection{Timing}

Berechnen Sie ein Matrix-Vektor-Produkt \mbox{einer} $(100 \times
100)$-Zufallsmatrix mit einem Zufallsvektor, indem Sie einerseits die
eingebaute Funktion \texttt{*}, andererseits \emph{for}-Schleifen
verwenden.  Vergleichen Sie die Rechenzeiten beider Lösungen.

\paragraph*{Lösungsvorschlag:}

\begin{verbatim}
  A = rand(100); b = rand(100,1);
  t = cputime;  
  v = A*b; t1 = cputime-t;
  w = zeros(100,1);
  t = cputime;
  for n = 1:100,
    for m = 1:100
      w(n) = w(n)+A(n,m)*b(m);
    end
  end
  t2 = cputime-t;

  norm(v-w), t2/t1
  ans = 0
  ans = 577.00
\end{verbatim}


\subsection{Stabilitätsfunktionen von BDF-Integratoren}

Berechnen Sie alle Nullstellen des Polynoms
\[
  \frac{147}{60}\zeta^6 -6\zeta^5 + \frac{15}{2}\zeta^4 
  - \frac{20}{3}\zeta^3 + \frac{15}{4}\zeta^2 - \frac{6}{5}\zeta +
  \frac{1}{6}
\]
Hinweis: Verwenden Sie den Befehl \texttt{compan}.

Zeichnen Sie diese Nullstellen in der kom\-plexen Ebene als Punkte
ein, und fügen Sie einen Einheitskreis hinzu (Hinweis: \texttt{hold,
real, imag}).
\paragraph*{Lösungsvorschlag:}

\begin{verbatim}
  bdf6 = [147/60 -6 15/2 -20/3 15/4 -6/5 1/6];
  R = eig(compan(bdf6));
  plot(R,'+'); hold on
  plot(exp(pi*i*[0:.01:2]));
  if any(find(abs(R)>1))
    fprintf('BDF6 ist instabil\n');
  else 
    fprintf('BDF6 ist stabil\n'); 
  end
\end{verbatim} 

\begin{center}
  \includegraphics[height=8cm,angle=270]{bdf6.ps}
\end{center}


\subsection{3D-Plot}

Zeichnen Sie den Graphen der Funktion 
\[
  f(x,y) = \exp(-x^2-y^2) .
\]

\paragraph*{Lösungsvorschlag:}

\begin{verbatim}
  x = -3:0.1:3;
  [xx,yy] = meshgrid(x,x);
  z = exp(-xx.^2-yy.^2);
  mesh(x,x,z);
  title('exp(-x^2-y^2)');
  replot;
\end{verbatim}


\subsection{Hilbertmatrix}
  
Berechnen sie für die $(n \times n)$-Hilbertmatrix $H$
($n=1,\dots,15$) die Lösung des linearen Gleichungssystems $Hx = b$,
$b = $ \texttt{ones(n,1)}. Berechnen Sie den Fehler und die Kondition
der Matrix und zeichnen Sie die beiden in einem Fenster
halblogarithmisch auf. (Hinweis: \texttt{hilb, invhilb}.)

\paragraph*{Lösungsvorschlag:}

\begin{verbatim}
  err = zeros(15,1); co = zeros(15,1);
  for k = 1:15
    H = hilb(k);
    b = ones(k,1);
    err(k) = norm(H\b-invhilb(k)*b);
    co(k) = cond(H);
  end
  semilogy(1:15,err,'r',1:15,co,'x');
\end{verbatim}


\subsection{Ausgleichsgerade}

Berechnen Sie zu beliebig vorgegebenen Punkten $(x_j,y_j)$, die durch
zwei Vektoren $x$ und $y$ gegeben sind, die Ausgleichsgerade (im Sinne
der kleinsten Fehlerquadrate). Zeichnen Sie die Punkte und die Gerade.

\paragraph*{Lösungsvorschlag:}

\begin{verbatim}
  function coeff = ausgl(x,y)
    n = length(x);
    A = [x ones(n,1)];
    coeff = A\y;
    plot(x,y,'x');
    hold on
    interv = [min(x) max(x)];
    plot(interv,coeff(1)*interv+coeff(2));
  end
\end{verbatim}


\subsection{Trapezregel}

Schreiben Sie ein Programm, das beliebige Funktionen $f$ in einer
Variablen auf einem Intervall $[a,b]$ numerisch integriert. Verwenden
sie die Trapezregel ($h = (b-a)/N$):
\[
  \int_a^b f(x) \, dx \approx h\left( \frac{f(a)+f(b)}{2} + 
  \sum_{j=1}^{N-1} f(a+jh) \right).
\]
Verifizieren Sie in einem doppellogarith\-mischen Bild anhand einer
beliebigen Funktion $f$, dass die Trapezregel Ordnung $2$ hat.

\paragraph*{Lösungsvorschlag:}

\begin{verbatim}
  function S = trapez(fun,a,b,N)
    h = (b-a)/N;
    % fy = feval(fun,[a:h:b]); besser:
    fy = feval(fun,linspace(a,b,N+1));
    fy(1) = fy(1)/2;
    fy(N+1) = fy(N+1)/2;
    S = h*sum(fy);
  end

  function y = f(x)
    y = exp(x);
  end  

  for k=1:15; 
    err(k) = abs(exp(1)-1-trapez('f',0,1,2^k));
  end
  loglog(1./2.^[1:15],err);
  hold on;
  loglog(1./2.^[1:15],err,'x');
  title('Trapezregel, f(x) = exp(x)');
  xlabel('Schrittweite');
  ylabel('Fehler');
  replot;
\end{verbatim} 
  
\end{document}
