kramann.info
© Guido Kramann

Login: Passwort:










kramann.info
© Guido Kramann

Login: Passwort:




Integration mit Hilfe von Scilab

1:1 Umsetzung des C++ Programms in Scilab

  • Scilab ist ein Open-Source-Programm mit ähnlicher Syntax und Funktionalität wie das kommerzielle Programm für wissenschaftliches Rechnen Matlab
  • Ein ähnliches Programm, auch Open-Source wäre noch Octave
  • Nach Start->scilab-4.0->scilab-4.0 erscheint eine Konsole, in der direkt eingegebene Befehle ausgeführt werden:
  • Enden sie mit einem Semikolon, so wird das Ergebnis einer Operation nicht angezeigt.
  • Die Stärke von Scilab/Octave/Matlab liegt in der schnellen Verarbeitung von Matrizen und einer riesigen Bibliothek an mathematischen Funktionen und zur Ergebnisdarstellung.
  • Schleifen laufen dagegen eher langsam ab.
  • ACHTUNG: Im Gegensatz zu C,C++, oder Java, beginnen hier die Indizes von Arrays mit 1 und nicht mit 0.
  • Mit eckigen Klammern angelegte Vektoren liegen.
  • Mit einem nachfolgenden Hochkomma werden sie transponiert.
Beispielbefehle in Scilab

Bild 0-1: Beispielbefehle in Scilab

  • Statt einzelner Befehle, kann auch ein Skript in einem File mit Endung .m zusammengestellt werden.
  • Im folgenden ist das vorhergehende C++ - Programm fast 1:1 als Scilab-Makro umgesetzt.
  • Es wird z.B. in Notepad geschrieben und unter einem Namen mit .sce oder .m - Endung abgespeichert.
  • Im Scilab-Fenster kann man dann in das entsprechende Verzeichnis mit cd DATEIPFAD gehen und das Makro hier mit exec euler.m starten.
  • Die Angaben 1:ii-1 sind übrigens nötig, weil etwas mehr Speicher von Scicos automatisch allokiert wird, als nötig.
A=20
dt=60
k=1/872542

sekunden_pro_tag = 86400
tage = 30
schritte = (tage*sekunden_pro_tag)/dt

schritte_pro_tag = round(sekunden_pro_tag/dt)

t=0
ii=1

for i=0:1:schritte
    if(modulo(i,schritte_pro_tag)==0)
        merke_t(ii) = t;
        merke_A(ii) = A;
        ii=ii+1;
    end
    t=t+dt;
    A = A*(k*dt + 1);
end

plot(merke_t(1:ii-1),merke_A(1:ii-1),'*blu-');

//siehe help axes_properties
a = gca();
a.x_label.text = 't/s';
a.y_label.text = 'A/m2';
a.title.text = 'Flächenwachstum eines Seerosenteiches';
 

Code 0-1: Scilab-Makro euler.m zur Euler-Integration des Seerosenteich-Modells

Beispiel in Scilab etwas übersichtlicher ohne Grafik.

Bild 0-2: Beispiel in Scilab etwas übersichtlicher ohne Grafik.

seerose_scilab.zip, enthält euler.m

Wie man sieht, kommt man hier sehr leicht an eine ansprechende Visualisierung der Ergebnisse. Die folgende Scilab-plot-Grafik wurde mit dem Befehl File->export als gif-Bild exportiert. Leider gibt es bei der Y-Achsen-Beschriftung beim Exportieren ein Problem.

Ergebnis der Scilab-Berechnung mit plot(merke_t(1:ii-1),merke_A(1:ii-1),'*blu-'); dargestellt.

Bild 0-2: Ergebnis der Scilab-Berechnung mit plot(merke_t(1:ii-1),merke_A(1:ii-1),'*blu-'); dargestellt.

Vergleich mit dem analytischen Ergebnis

  • Nach Ablauf des Makros werden nun noch die folgenden Zeilen eingegeben, um auch die analytische Berechnung durchzuführen:
Aneu = 20*exp(merke_t(1:ii-1)/872542);

plot(merke_t(1:ii-1),merke_A(1:ii-1),merke_t(1:ii-1),Aneu);
 

Code 0-2: Vergleich zwischen analytischer und numerischer Berechnung

Übereinander liegende analytische und numerische Ergebniskurve.

Bild 0-3: Übereinander liegende analytische und numerische Ergebniskurve.

  • Die Beiden Kurven liegen augenscheinlich weitgehend deckungsgleich übereinander. Mehr Informationen bringt die Darstellung der absoluten Differenz beider Ergebnis, also quasi der absolute Fehler der numerischen Integration:
  • Folgende Befehle liefern die gewünschte Kurve:
    Fehler = abs( merke_A(1:ii-1) - Aneu );

    plot(merke_t(1:ii-1),Fehler);
 

Code 0-3: Absolute Differenz zwischen analytischer und numerischer Berechnung

Verlauf des absoluten Fehlers der numerischen Integration.

Bild 0-4: Verlauf des absoluten Fehlers der numerischen Integration.

Hier wird schon deutlich, dass sich der numerische Fehler bei der numerischen Integration mit dem Euler-Verfahren Schritt für Schritt aufsummiert und somit exponentiell anwächst.