Ruprecht-Karls-Universität Heidelberg

This is an archive of our old website and not updated anymore

The current website can be found at conan.iwr.uni-heidelberg.de

Octave

In den Übungen werden wir hin und wieder 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 (es gibt auch eine Windows-Variante). Das Programm kann in der Linux-Konsole mit

octave
gestartet werden.

Erste 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 der wichtigsten 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:
    
      help .*
      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 (beide Komponenten in einem Diagramm) geplottet werden kann

      
        % Define System:
        %     x' =  2 * y * (1 + x*x + y*y)
        %     y' = -4 * x * (1 + x*x + y*y)

        % set an inline function representing the system
        g = inline('[ 2 * y(2) * (1 + y(1)*y(1) + y(2)*y(2)); -4*y(1) * (1 + y(1)*y(1) + y(2)*y(2)) ]', 't', 'y');

        % plot solution (the classical way)
        %[t, y] = ode45(g, [0, 20], [1; 0]);
        %plot(t, y)

        % solve with Runge-Kutta (ode45),
        %     second parameter: time interval
        %     third           : initial conditions
        [t, y] = ode45(g, [0, 10], [0.2; 0]);

        % plot a phase portrait
        plot(y(:,1), y(:,2))
        grid;
        xlabel('y1')
        ylabel('y2')
        title("Phasen-Portrait eines einfachen Systems")
        print("portrait1.eps", "-color", "-deps")
      
    

Phasen-Portraits mit Nullklinen und Trajektorien

Der folgende Code zeigt, wie ein zweidimensionales dynamisches System mit Nulllkinen, Richtungsfeld und einigen Trajektorien, berechnet mit einem Forward-Euler-Verfahren, geplottet werden kann. Zunächst wird das dynamische System mit einer Funktion definiert, diese wird dann als Kernel an die Phasenportrait-Funktion übergeben. Der Kernel:

      
        % sample_field1.m
        % 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
      
    
Die Plot-Funktion:
      
        % phase_portrait.m
        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;
      
    
Aufruf in Octave z.B. mit:
      
        phase_portrait('sample_field1',-2:0.1:2,-2:0.1:2,200,0.02,30)
      
    


Revision: 1305     Letzte Änderung: 2012-08-09 10