Octave
In den Übungen werden wir das Mathematik-Algebrasystem Octave verwenden, das ein Open Source Clone von Matlab ist. Es ist im Pool installiert und kann von der Octave-Homepage heruntergeladen werden. Das Programm kann in der Linux-Konsole mit
octaveErste Schritte
In der Kommandozeile von Octave kann man einfach einmal folgendes ausprobieren (nachdem man die Octave-Kommandozeile wie oben angegeben gestartet hat):
x = 3.1415
format long
x = 3.1415
x = 3.1415;
A = [ 1 2; 3 4; ]
b = [ 1; 2; ]
A \ b
A^2
A.^2
A.' * A
Dies sind einige grundlegende Kommandos, mit denen man schon recht
viel arbeiten kann. In den ersten drei Zeilen wird eine Variable x angelegt
und auf dem Bildschirm ausgegeben, zunächst mit niedriger, dann mit hoher
Genauigkeit. Zeile vier legt die Variable nur an, ohne sie auf dem Bildschirm
auzugeben, die Ausgabe wird durch das Semikolon; am Zeilenende
unterdrückt. In der nächsten Zeile wird eine Matrix und ein Spaltenvektor
angelegt. Ohne die Semikolon in den eckigen Klammern wäre b ein
Zeilenvektor. Die Anweisung A \ b löst das Gleichungssystem A x =
b mit Ausgabe am Bildschirm. Die nächste Zeile berechnet A^2 = A *
A. Im Vergleich dazu arbeitet der gepunktete Potenz-Operator .^
der folgenden Zeile Element-Weise, d.h. jeder Matrix Eintrag wird separat
quadriert. Damit ist die letzte Zeile die Berechnung des Produkts aus
transponierter Matrix und Matrix.
Zu allen Kommandos gibt es übrigens Hilfe in der
Octave-Kommandozeile, z.B.:
help .*
help save
help format
Plotten von Funktionen und Export in Dateien
Funtktionen können recht einfach geplottet werden:
x = [ 0 : 0.1 : 3 ];
y = [ 0 : 0.4 : 2 ];
[ xx, yy ] = meshgrid(x, y);
z = sin(xx.^2 - yy.^2);
mesh(x,y,z)
print ("sin.png", "-color", "-dpng");
Die ersten beiden Zeilen legen die Plotbereiche und die jeweilige Auflösung
fest. Zeile drei erzeugt ein Gitter und in Zeile vier wird für jeden
Gitterpunkt (x,y) der Funktionswert z = sin(x*x - y*y) berechnet
und gespeichert. mesh plottet die Funktion und mit
print() wird anschließend in eine png-Datei geschrieben. Die
plot-Funktion erhält dazu den Dateinamen, ein Flag und die Treiberwahl für die
Ausgabe als Parameter.
Skripten
Alle Befehle können hintereindander in ein Skriptfile geschrieben werden, das die Endung .m haben und im gleichen Verzeichnis, in dem Octave gestartet wurde, liegen sollte. Dann kann das Skript in Octave einfach durch Aufruf des Skriptnamens ohne Endung aufgerufen werden. Stehen die Plot-Befehle des vorigen Abschnitts beispielsweise in plotsin.m,
x = [ 0 : 0.1 : 3 ];
y = [ 0 : 0.4 : 2 ];
[ xx, yy ] = meshgrid(x, y);
z = sin(xx.^2 - yy.^2);
mesh(x,y,z)
print ("sin.png", "-color", "-dpng");
, dann kann man das Skript in der Octave-Kommandozeile
octave:1>plotsin
ausgeführt werden.
Zeitmessung
Die verbrauchte CPU-Zeit kann in Octave mittels
t1 = cputime;
... tue was;
t2 = cputime - t1;
gemessen werden.
Phasenraum-Portraits
Der folgende Code zeigt, wie ein einfaches zweidimensionales dynamisches System mit einem Runge-Kutta-Verfahren gelöst und ein Phasenportrait (das Vektorfeld bestehend aus den Steigungen x'/y' an den Punkten (x,y) wird geplottet) erstellt werden kann.
% Define System:
% x' = 2 * y * (1 + x*x + y*y)
% y' = -4 * x * (1 + x*x + y*y)
[x1, x2] = meshgrid(-.5:.05:.5, -.5:.05:.5);
x1dot = 2 * x2 .* (1 + x1.^2 + x2.^2);
x2dot = -4 * x1 .* (1 + x1.^2 + x2.^2);
% plot vector field
quiver(x1, x2, x1dot, x2dot)
grid
xlabel('x1');
ylabel('x2');
% some output
title("Phasen-Portrait eines einfachen Systems")
print("phaseportrait.pdf", "-color", "-dpdf")
Phasen-Portrait mit Nullklinen und Trajektorien
Der folgende Code zeigt, wie ein zweidimensionales dynamisches System mit Nullklinen, Richtungsfeld und einigen Trajektorien, berechnet mit dem Forward-Euler-Verfahren, geplottet werden kann. Zunächst wird das dynamische System mit einer Funktion definiert, das als Kernel an die Plotfunktion übergeben wird:
% define system
Function [U,V] = sample_field1(X, Y)
U = X.^3 - 2 * X.^2 - 3 * X + 1 - Y.^2;
V = Y.^3 - 2 * Y.^2 - 3 * Y + 1 - X.^2;
endfunction
% plotting function
function r = phase_portrait(F, X, Y, traj, tau, time)
% (author: F. Piekniewski)
% Plots a vector field, nullclines and sample trajectories.
%
% Parameters:
% F function handle that computes the vector fields
% X, Y ranges of solutions U and V
% traj number of cmoputed trajectories
% tau timestep
% time time interval to compute trajectories
%
% sample call:
% phase_portrait('sample_field1',-2:0.1:2,-2:0.1:2,200,0.02,30)
% plot bounding box
xl = min(X(:));
xu = max(X(:));
yl = min(Y(:));
yu = max(Y(:));
% evaluate function at mesh points
[X,Y] = meshgrid(X,Y);
[U,V] = feval(F,X,Y);
% make contour plot (nullclines)
figure;
hold on;
contour(X, Y, U, [0 0], 'g')
contour(X, Y, V, [0 0], 'r')
% plot vector field
quiver(X, Y, U, V);
% compute and plot trajectories
for i = 1 : traj
% set initial values randomly
x0 = ((xu-xl) * rand(1,1)) + xl;
y0 = ((yu-yl) * rand(1,1)) + yl;
tr=[];
tr(1,:) = [x0, y0];
% compute current trajectory
for t = 0 : tau : time
% evaluate eqaution, FE step
[dx,dy] = feval(F, tr(end,1), tr(end,2));
tr(end+1,:) = tr(end,:) + tau * [dx,dy];
% limit reached, reset tr
if (tr(end,1)>xu) || (tr(end,1)<xl) || (tr(end,2)>yu) || (tr(end,2)<yl)
tr=tr(1:end-1,:);
break;
end;
end;
% plot trajectory
plot(tr(:,1), tr(:,2), 'k');
end;
% axis limits for plot
axis([xl xu yl yu]);
endfunction;
% run
phase_portrait('sample_field1',-2:0.1:2,-2:0.1:2,200,0.02,30)
