Übungen zur Vorlesung Scientific Computing
Scilab Tutorial
Scilab ist ein freier MATLAB-clone entwickelt vom französischen Forschungszentrum INRIA.
Tippen Sie die Zeilen in Scilab ein und versuchen Sie mit den Kommentaren und der Ausga-be in Scilab, die Anweisungen zu verstehen. Bei Unklarheiten in der eingebauten Hilfe (helpBefehlsname) nachschlagen oder den Übungsleiter fragen.
zz = 4 + 3*%i // eine komplexe Zahlconj(zz)sqrt(zz*conj(zz)) == abs(zz)imag(zz)real(zz)
exp(-%i * %pi) // vordefinierte Konstanten beginnen mit %// siehe "help percent" für mehr Info
// Vektor- und Matrixoperationen ///////////////////////////////////////////
v = [0.1, 0.3, 0.8, -1.2] // Zeilenvektor (Kommas können weggelassen werden)
// zwei Arten Zeilenvektoren mit linearem Abstand zu definieren:x = 1:0.5:3 // Werte aus dem Intervall [1,3] mit Abstand 0.5y = linspace(0,1.0,5) // 5 Werten linear auf das Intervall [0,1] verteilt
z = x/3
x(1) // Vektorindizes beginnen mit 1!x(3)x($) // $ bezeichnet Index des letzten Elementes
sqrt(v) // Aufruf einer Funktion mit einem Vektorargument
a = rand(2,2) // generiere eine konstante 2x2 Matrix aus Pseudo-Zufallszahlenb = a >= 0.5 // der Vergleich erzeugt eine bool’sche Matrix (wahr/falsch)
m = [x; z] // baue eine Matrix aus zwei Zeilenvektorenm(2,3) // Matrixelement in Zeile 2 und Spalte 3u = [2; 5; 7; 11] // ein Spaltenvektor
clear m;m = [ 2, 4, 8; 0, x(2:3) ] // x(2:3) ergibt einen Teilvektor des Vektors xm’ // die transponierte Matrixz = m(1,:) // die erste Zeile
Sommersemester 2009 1
Übungen zur Vorlesung Scientific Computing
m(:,$) // die letzte Spalte
[rowidx, colidx] = find(m > 2) // liefert Indizes, für die m(rowidx, colidx) > 2m(rowidx, colidx) // die Matrixelemente > 2
max(m)min(m)
expm(a) // Matrix-Exponentialfunktion expm(a) = I + a + a^2 / 2 + ...expm(a) ~= exp(a) // Test auf Ungleichheit; ~= ist äquivalent zu <>
[rows, cols] = size(m)
length(x) == length(y) // Test ob die beiden Vektoren gleiche Länge haben
//x * y // nicht sinnvoll definiert!x .* y // elementweise Multiplikationx * y’ // Skalarprodukt (Zeile * Spalte)x’ * y // erzeugt eine Matrix (äußeres Produkt, Tensorprodukt)
ones(2,2)zeros(1,3)eye(3,3) // ~= eye(3)!diag([1,3,5,7]) // oder diag([1 3 5 7])diag(m)
// Erzeuge eine tridiagonale Matrix der Größe 2*n+1n=5; T = diag(-n:n) + diag(ones(2*n,1),1) + diag(ones(2*n,1),-1)
A=[1,2;3,4]; b=[5;6];x = A \ b // Löse lineares Gleichungssystem A*x = bnorm(A*x - b)
a ~= A // Scilab ist ’case-sensitive’!
det(A) // Determinanteinv(A) // Inversetrace(A) // Spurcond(A) // Konditionszahl - ein Maß dafür wie "singulär" die Matrix ist
[X, lambda] = spec(A) // bestimme Eigenwerte lambda und Eigenvektoren X von A
[L,U,P] = lu(A) // LU (LR)-Faktorisierung von A: PA = LUP*A == L*U
Sommersemester 2009 2
Übungen zur Vorlesung Scientific Computing
A=[1 2 3; 2 4 6; 7 3 5]; // eine singuläre Matrix
ker_A = kernel(A) // Nullraum, Kern von A[U,dim_im_A] = range(A,1) // Bildraum von A, im(A)
// die ersten dim Spalten von U’ spannen im(A) aufrank(A) // Rang von Adim_im_A == rank(A) // dim im(A) == Rang(A)size(ker_A,’c’) + rank(A) == 3 // Dimensionsformel
// Polynomeclear x;p = poly([1 3], ’x’) // definiert durch Nullstellenq = poly([1 2], ’x’, ’c’) // definiert durch Koeffizientenroots(q) // Nullstellen von q
p+q+1[p, 1-x; 1, p*q + x] // polynomiale Matrixr = 1 + 3*x + 2*x^2r/q // rationaler Ausdruck
// Funktionen //////////////////////////////////////////////////////////////
// Entweder direkt in Scilab eingeben oder in "Skript-Files" mit// Dateiendung ".sce" bzw. ".sci" (letztere für reine Funktionsdefinitionen)// Laden von Skript-Files mit "exec" bzw. "getf"
// Beispiel für zwei rekursive Funktionen die die Kontrollstrukturen// "if" bzw. "select" verwenden.
function y=fact(x)if x <= 1 theny=x;
elsey=x*fact(x-1);
endendfunction
function f = fibonacci(n)select ncase 0 thenf = 1
case 1 thenf = 1
else
Sommersemester 2009 3
Übungen zur Vorlesung Scientific Computing
f = fibonacci(n-1) + fibonacci(n-2)end
endfunction
// Vektorisierung ist meist schneller als übliche "for"/"while"-Schleifen:
// a) initialisiere einen Vektor über for-Schleife
n=10000;
timer(); // messe cpu-timeu = [];
for i=1:nu(i) = sin(i * %pi / n);
end
t1 = timer();
// b) selbe Initialisierung vektorisiert
timer();x = linspace(0,%pi,n);v = sin(x);t2 = timer();
// Vergleich der cpu-time:disp("cpu-times: "); // unformatierte Ausgabe[t1 t2] // wird automatisch angezeigt, weil kein ";" am Endemprintf("speedup = %.2f\n", t1/t2); // formatierte Ausgabe: syntax wie printf in C
// Regel: immer vektorisieren, wenn möglich;// insbesondere bei aufwendigen Berechnungen die oft ausgeführt werden
// Bsp. für "while" - Schleifex=1;while exp(x) <> %inf;x=x+1;
end
mprintf("exp(x) == inf (overflow!) for x = %d\n", x);// tabellenartige Ausgabemprintf("exp(%d) \t exp(%d) \n %e \t %e \n", x-1, x, exp(x-1), exp(x));
// Mit "break" kann man Schleifen vorzeitig verlassen// zur Illustration:
Sommersemester 2009 4
Übungen zur Vorlesung Scientific Computing
while %t // Endlosschleife// einige Anweisungen ...if Bedingung thenbreak;
endend
// Logische Operatoren// Mit & bzw | können Bedingungen mit logischem und bzw. oder verknüpft werden.// ~ steht für logische Verneinung
[a,b,c] = (rand(), rand(), rand());if (a > b & b > c) | (a < b & b < c) thendisp(’Tripel (a,b,c) ist geordnet’);
end~ (1 > 2)
// Graphik //////////////////////////////////////////////////////////////////
x = linspace(0,%pi,20)’; // Spaltenvektorplot2d(x,sin(x));help plot2d // konsultiere eingebaute Hilfeplot2d() // democlf(); // Grafikfenster löschen
clear f;function y=f(t); y=sin(t).*cos(t); endfunction // elementweise Operation!plot2d(x, [sin(x), f(x)], style = [-1, 2]);
// parametrischer Plot einer Kurveclf();t = linspace(0, 20*%pi, 1000);param3d1(sin(t), cos(t), exp(t/20).*t/50, 60, 170, ’x@y@z’)
// Plot einer Fläche im R^3clf();x = linspace(-%pi, %pi, 40);y = linspace(-%pi, %pi, 40);plot3d(x, y, sinh(x’)*cos(y)); // beachte: 3. Parameter = Matrix!
// Vektorfeld im R^2 definiert durch ein System von gewöhnlichen DGLn// mathematisches Pendel \ddot phi(t) = - sin phi(t)// als System 1. Ordnung:// Variablen: theta = [theta(1); theta(2)]// theta(1) = Auslenkungswinkel, phi(t)
Sommersemester 2009 5
Übungen zur Vorlesung Scientific Computing
// theta(2) = Zeitableitung von Auslenkungswinkel, \dot phi(t)//function [thetadot] = math_pendel(t,theta);thetadot = [theta(2); - sin(theta(1))];
endfunctiontheta1 = linspace(-%pi,%pi,15);theta2 = linspace(-%pi,%pi,15);clf();fchamp(math_pendel, 0, theta1, theta2); // = Phasenraumportrait!
// Histogramm von normalverteiltem Vektor und Normalverteilungclf();v = rand(1,2000,’n’);histplot([-4.5:0.25:4.5], v, style=2);function [y]=f2(x); y=exp(-x.*x/2) / sqrt(2*%pi); endfunctionx=-4.5:0.125:4.5; x=x’;plot2d(x, f2(x), 26, "000")xtitle(’normalized histogram plot’, ’classes’,.. // ".." = Zeilenfortsetzung
’number of data components in class (normalized)’);legends([’gaussian random sample histogram’,..
’exact gaussian density’],..[2 26], ’ur’)
// lineare Ausgleichsprobleme ///////////////////////////////////////////////clf();x = (1:10)’;y1 = 3*x + 4.5 + 3*rand(x,’normal’); // baue "fehlerbehaftete" Dateny2 = 1.8*x + 0.5 + 2*rand(x,’normal’);plot2d(x, [y1,y2], [-2,-3]);
A = [x, ones(x)]; b = [y1, y2];X = lsq(A,b); // lineare Regression
y1e = X(1,1)*x + X(2,1);y2e = X(1,2)*x + X(2,2);plot2d(x, [y1e, y2e], [2,3]);
// Nichtlineare Gleichungen /////////////////////////////////////////////////function y=g(x); y=sin(x**2) - 0.2; endfunctionx0 = 0; // Startwertroot = fsolve(x0, g) // suche Nullstelle von g(x)
Sommersemester 2009 6
Übungen zur Vorlesung Scientific Computing
// Interpolation ////////////////////////////////////////////////////////////x = [1 10 20 30 40];y = [1 30 -10 20 40];clf();plot2d(x’, y’, [-3], "011", " ", [-10, -40, 50, 50]);yi = interpln([x;y], -4:45); // lineare Interpolationplot2d((-4:45)’, yi’, [3], "000");yy = interp1(x, y, -4:45, ’spline’); // stückweise kubische Interpolationplot2d((-4:45)’, yy, [4]);
// Diskrete Fouriertransformation ///////////////////////////////////////////// [Motivation: kann durch Abschneiden einer Fourierreihe erzeugt werden// (dh. nur endlich viele Glieder); die Koeffizienten erhält man durch// Approximation des Integrals durch eine Summe.]
// Frequenzkomponenten eines Signals// Erzeuge ein verrauschtes Signal gesampled bei 1000hz, das// Sinusschwingungen bei 50 und 70 Hz enthält//// In der Praxis können Signale ja nur mit einer endlichen Abtastrate 1/\Delta t// aufgezeichnet werden. DFT ist diskrete Version der Fourier Transformation.
clear f;sample_rate = 1000;t = 0:1/sample_rate:0.6;N = size(t,’*’); // Anzahl der Sampless = sin(2*%pi*50*t) + sin(2*%pi*70*t + %pi/4) + grand(1,N,’nor’,0,1);clf()plot2d(t, s);y = fft(s); // schnelle Fouriertransformation// Nachdem die fft symmetrisch ist brauchen wir nur die ersten N/2 Punktef = sample_rate * (0:(N/2)) / N; // Frequenzvektorn = size(f,’*’)clf()plot2d(f, abs(y(1:n)));
// Numerische Integration ///////////////////////////////////////////////////clear f;function y=f(x), y = x*sin(30*x) / sqrt(1 - ((x/(2*%pi))^2)), endfunctionI_exact = -2.5432596188;I = intg(0, 2*%pi, f) // bestimmtes Integralabs(I_exact - I)
t = 0:0.1:%pi;
Sommersemester 2009 7
Übungen zur Vorlesung Scientific Computing
inttrap(t, sin(t)) // bestimmtes Integral über Trapezregel
// Numerische Differentiation ///////////////////////////////////////////////
// Berechne die Ableitung von sin(x) näherungsweisedx = 0.001;x = 0:dx:10; // Intervall [0,10]y = sin(x);dy = diff(y) / dx; // diff(y) liefert den Vektor y(2:$) - y(1:$-1)// Wir vergleichen dy mit dem bekannten Ausdruck für die Ableitung cos(x) und// berechnen die Unendlich-Norm (Maximumsnorm) von der Differenz:norm(dy - cos(x(1:$-1)), %inf)
function [y]=h(x); y=[x(2); -x(1)**2]; endfunction // Vektorwertige Funktionnumdiff(h, [3 2]) // Jacobimatrix von g ausgewertet bei x = [3 2]
numdiff(sin, 0) // Ableitung von sin(x) bei x=0
// Gewöhnliche DGLn ///////////////////////////////////////////////////////////// dy/dt = y^2 - y sin(t) + cos(t) = 0 // skalare DGL// Anfangsdaten y(0) = 0clear f;function ydot=f(t,y); ydot=y^2 - y*sin(t) + cos(t); endfunctiony0 = 0; t0 = 0; t = t0:0.1:%pi;y = ode(y0, t0, t, f)clf();plot(t, y)
// Probiere auch einige Demos (über das Menü aufrufbar) aus und studiere die// Dokumentation um einen weiteren Eindruck von der Funktionalität von Scilab// zu bekommen.
Literatur:
• Scilab Group, Introduction to Scilab, http://www.scilab.org/doc/intro, INRIA 1998, (ziem-lich veraltet, aber trotzdem eine gute Einführung)
• Bruno Pinçon, Eine Einführung in Scilab, http://www.scilab.org/publications/JARAUSCH/PinconD.pdf,(auch alt, aber sehr umfangreich)
• Weitere Dokumentation unter http://www.scilab.org/product
Sommersemester 2009 8
