Verfasst am: 03.07.2011, 11:17
Titel: Zeitpropagation der Dichtematrix unter Verwendung von ode45
Hallo allerseits!
Im Rahmen meiner Bachelorarbeit muss ich mich nun zum ersten mal mit Matlab auseinandersetzen und habe schon zahlreiche Tutorials und Skripte gelesen um die Grundfunktionen zu erlernen.
Es geht darum, dass ich eine lineare DGL erster Ordnung für eine 16x16 Dichte-Matrix (Von-Neumann-Gleichung) mittels der ode45-Funktion lösen muss.
Da die Funktion ja nur mit Vektoren, nicht aber mit Matrizen als Eingabe arbeiten kann, übergebe ich die Matrix in Form eines 256 langen Spaltenvektors. Dazu habe ich zunächst die reshape-Funktion verwendet, danach jedoch eigene Methoden zur Multiplikation zweier Matrizen - eine in dieser 256-Vektorrepräsentation und die andere als normale 16x16-Matrix - geschrieben, die bei meinen Tests schneller als die reshape-Funktion liefen.
Die Lösung für die Dichtematrix soll auf einem Zeitgitter der Größe GRID berechnet werden. Die Matrix rho enthält dazu GRID Zeilen in denen jeweils die komplette Dichtematrix als 256-Vektor zum der Zeile entsprechenden Zeitpunkt steht.
Das Problem ist, dass das Programm viel zu langsam läuft. Tests mit einer 4x4 statt der 16x16 Matrix zeigten, dass zumindest die Berechnung funktioniert. Dazu habe ich im if-Zweig in der for-Schleife in der Methode FProp jedesmal das i ausgeben lassen. Da dauerte jeder Berechnungschritt allerdings ca. 15 Sekunden, das ganze mal 100 (Anzahl der Zeitgitterpunkte) macht 150 Sekunden. Meine tatsächliche 16x16-Matrix ist 16 mal größer, daher vermutlich auch 16-fache Laufzeit, d.h. 40 Minuten für eine einzige Propagation (d.h. einen kompletten Durchlauf der Methode FProp für alle 100 Zeitgitterpunkte). Das ist absolut inakzeptabel und sollte deutlich unter einer Minuten, am besten bei wenigen Sekunden liegen.
Im 16-dimensionalen Fall habe ich probeweise bei der Rückgabe von drhoVdt in der Funktion seq das "+ liouville(rhoV)" entfernt und das Programm läuft wesentlich schneller durch, so dass die lange Laufzeit wohl auf die Methode liouville zurückzuführen ist. Aber wie soll ich dort was optimieren?
Ich habe nun schon wirklich alles ausprobiert und auch mein Betreuer weiß nicht mehr weiter, was man da noch tun könnte. Daher wollte ich hier mal fragen, ob ihr eine Idee habt wie ich das Programm wesentlich schneller mache:
% Vektor mit den Schrittpunkten von 0 bis TFINAL grid = TSTEP*[0:GRID-1];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Festlegung der physikalischen Konstanten
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
global gamma1
gamma1 = 1;
global gamma2
gamma2 = 1;
global gamma3
gamma3 = 1;
global RabiRot
global RabiBlau
global E0
global E1
global Delta1
global Delta2
global U
RabiRot=4.5594895*10^-8;
RabiBlau=3.7995746*10^-9;
E0=0;
E1=1033*10^-9;
Delta1=9.1189791*10^-8;
Delta2=0;
U=7.5991492*10^-9;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Hamilton-Operator und Dissipations-Matrizen fuer Liouville-Superoperator
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
gam=1.801*10^9;
EPSILON = 0.1*gam;
epsilon0 = zeros(GRID,1);
for i = 0:GRID-1
epsilon0(i+1) = EPSILON*exp(-(i-50)*(i-50)/(25*25));
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Ausfuehrung der Vorwaertspropagation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% rho ist eine Matrix bei der die Zeile den Gitterpunkt angibt (1. Zeile: % Werte der Matrix am 1. Gitterpunkt (t=0) usw. % Die eigentliche Dichtematrix steht als 256 Elemente langer Zeilenvektor in rho
function ret = S1(t) % Shape-Funktion des ersten Lasers global TFINAL
ret = sin(pi/TFINAL*t)*sin(pi/TFINAL*t);
function ret = S2(t) % Shape-Funktion des zweiten Lasers global TFINAL
ret = sin(pi/TFINAL*t)*sin(pi/TFINAL*t);
function ret = rproduct(A,b) % Produkt zwischen Matrix A und Repraesentant der Matrix B als Spaltenvektor b % Ergebnis ist auch Spaltenvektorrepraesentant der eigentlichen % Ergebnismatrix
matDim = size(A);
bDim = size(b);
rowSize = matDim(1);
colSize = matDim(2);
bSize = bDim(1);
for m = 1:bSize/colSize
k = colSize*(m-1);
for i = 1:rowSize
ret((m-1)*rowSize+i) = 0;
for j = 1:colSize
ret((m-1)*rowSize+i) = ret(k+i) + A(i,j)*b(j+k);
end end end
ret = ret';
function ret = lproduct(b,A) % Produkt zwischen Repraesentant der Matrix B als Spaltenvektor b und % Matrix A % Ergebnis ist auch Spaltenvektorrepraesentant der eigentlichen % Ergebnismatrix
matDim = size(A);
bDim = size(b);
rowSize = matDim(1);
colSize = matDim(2);
bSize = bDim(1);
for j = 0:colSize-1 for m = 1:bSize/rowSize
k = colSize*(m-1);
ret(m+j*bSize/rowSize) = 0;
for i = 0:rowSize-1
ret(m+j*bSize/rowSize) = ret(m+j*bSize/rowSize) + b(bSize/rowSize*i + m) * A(i+1,j+1);
end end end
ret = ret';
um aufzuspüren, warum Code lange braucht, kannst du den Profiler verwenden. Unnötige for-Schleifen oder dynamisch wachsende Variablen (orange unterringelt) bremsen Code oft unnötig aus.
Konkret: rproduct und lproduct können mit ziemlicher Sicherheit mit Matrix-Vektor-Operationen umschrieben werden. Leider werde ich aus dem Code dort nicht schlau, deswegen habe ich nur mal ein wenig herumprobiert:
Für (3x1)-Vektoren b und (3x3)-Matrizen A scheint
rproduct(A,b) --> A*b
lproduct(b,A) --> (b'*A)'
Bei dir scheinen b 256x1-Vektoren und A 16x16-Matrizen zu sein. Dann scheint folgendes zu funktionieren:
rproduct(A,b) --> reshape(A*reshape(b, 16, 16), 256, 1);
lproduct(b,A) --> reshape(reshape(b, 16, 16)*A, 256, 1);
Bei meinen Versuchen auf die Schnelle bringt das Faktor 10.
Zudem ist das Arbeiten mit globalen Variablen generell nicht empfohlen, da es den Code unübersichtlich macht.
Desweiteren ist es äußerst gefährlich, i als Laufvariable zu verwenden UND mit komplexen Zahlen zu rechnen.
Ist es tatsächlich gewünscht, dass Zeile 132 nicht von i abhängt? Mir scheint, da wird sehr oft dasselbe gemacht?
Die Differentialgleichung scheint steif/starr zu sein. Versuche also auch ode23s oder ähnliche statt ode45.
Schau doch mal, wie dir das weiterhilft.
Wenn du weitere Unterstützung bei der Codeoptimierung brauchst, poste bitte eine abgespeckte Version, die in ein paar Minuten durchläuft.
Danke für den Hinweis auf den Profiler, den kannte ich noch nicht. Den Code in rproduct und lproduct habe ich so aus dem Bauch heraus überlegt, die Methoden arbeiten aber auf jeden Fall richtig, habe sie ausgiebig getestet.
Aber wenn ich statt dessen, wie von dir vorgeschlagen doch die reshape-Funktionen verwende, ist die Laufzeit doch tatsächlich deutlich kürzer. Ich meine mich zwar zu erinnern, dass das bei vorigen Tests von mir (ohne den Profiler, sondern mit ner Stoppuhr in der Hand ) anders war, aber da täusche ich mich ja nun offensichtlich.
Mir ist schon klar, dass der Programmierstil absolut grauenhaft ist, aber zunächst wollte ich sicherstellen, dass das Programm überhaupt funktioniert, bevor ich da was dran mache. Ich hatte diese Propagationsmethode von jemand anderem übernommen, der an einem ähnlichen Problem gearbeitet hat und das für meinen speziellen Hamilton-Operator und Dichte-Matrix anzupassen versucht. Von ihm habe ich auch diese ganzen globalen Variablen übernommen.
Was ich jedoch nicht verstehe, ist Folgendes: Lasse ich die Ausführung der liouville-Funktion in Zeile 148 weg, läuft das Programm im Profiler in weniger als zwei Sekunden durch, seq wird dabei 3050 mal aufgerufen.
Ich habe in der Hauptmethode die liouville-Methode 100mal direkt hintereinander zu Testzwecken aufgerufen und bekomme da eine Gesamtzeit von 0,027 s für diese Methode angezeigt.
Füge ich die Methode nun wieder an die ursprüngliche Stelle in Zeile 148 wie im obigen Code ein, so sollte sie ja 3050 mal aufgerufen werden, da seq eben so oft aufgerufen wird. Die erwartete zusätzliche Zeit wären doch dann 3050 * 0,027 s / 100 = 0,82 s.
Füge ich die Zeile aber wieder ein und lasse starte das Programm, so wird es wieder nicht fertig / läuft nicht durch, nichtmal für ein einziges i in der Schleife in Zeile 129.
Dass in Zeile 132 kein i auftaucht, ist übrigens gewollt. Jedes Mal wenn FProp aufgerufen wird, werden alle Zeilen aus dem OldRho in das NewRho kopiert ohne was zu verändern, außer bei einer speziellen Zeile (die wo i == (t0/TSTEP) + 3 gilt). Dort wird die entsprechende Zeile in NewRho dann eben mit der ode45-Funktion berechnet, was dann nur von der beim Aufruf von FProp übergebenen Zeitmarke time abhängt.
Die anderen ode's habe ich auch schon durchprobiert, leider mit dem selben Ergebnis.
Füge ich die Methode nun wieder an die ursprüngliche Stelle in Zeile 148 wie im obigen Code ein, so sollte sie ja 3050 mal aufgerufen werden, da seq eben so oft aufgerufen wird.
Nein. ode45 und ähnliche verwenden eine dynamische Schrittweitensteuerung, um den Simulationsfehler zu kontrollieren. Es kann durchaus sein, dass das Hinzufügen eines zusätzlichen Terms das dynamische Verhalten des Systems und damit auch seine Simulation, insbesondere Schrittweitensteuerung, gravierend beeinflusst.
Wenn du nur die Zeit für etwas stoppen willst, kannst du auch TIC, TOC verwenden.
Zitat:
Mir ist schon klar, dass der Programmierstil absolut grauenhaft ist, aber zunächst wollte ich sicherstellen, dass das Programm überhaupt funktioniert, bevor ich da was dran mache.
Im allgemeinen ist es kaum zusätzlicher Aufwand, ein Programm beim ersten Mal gleich "richtig" zu schreiben. Jedenfalls weniger Aufwand, als das im Nachhinein zu machen.
Zitat:
Die anderen ode's habe ich auch schon durchprobiert, leider mit dem selben Ergebnis.
Welches Ergebnis? Ich habe mir in den Iterationen mal t anzeigen lassen, und ode23s / ode15s scheinen deutlich schneller zu simulieren als ode45. Auch wenn es immer noch weit weg davon ist, bald fertig zu werden.
richtig klasse. da hat sieht man graphisch aufbereitet, wo evtl. Flaschenhälse sind.
noch was zum generellen: Matlab bedeutet ja ursprünglich Matrix Laboratory, der Vorteil ist dass man auch ohne for's komfortabel mit Feldern arbeiten konnte. Und die entsprechenden Funktionen sind sogut wie IMMER schneller wie ein for-Konstrukt. Die Vermeidung dessen nennt sich Vektorisierung (engl.: vectorisation).
_________________
Ich hasse es wenn die Leute Fragen stellen, man dann versucht sich Mühe zu geben, und diejenigen ihren Thread nie wieder besuchen...
Welches Ergebnis? Ich habe mir in den Iterationen mal t anzeigen lassen, und ode23s / ode15s scheinen deutlich schneller zu simulieren als ode45. Auch wenn es immer noch weit weg davon ist, bald fertig zu werden.
Mit Ergebnis meinte ich, dass es auch mit diesen Funktionen nicht durchläuft. Aber die sind bei mir im Test auch beide langsamer als der ode45. ode15s ruft seq zB 4,5-mal häufiger auf als ode45 und braucht dabei etwa 5mal so lang insgesamt.
Ich bin mir relativ sicher, dass die lange Laufzeit an den vielen Matrix-Multiplikationen in liouville liegt. Habe aber langsam echt keine Idee mehr, was ich da noch machen soll.
Habe aber langsam echt keine Idee mehr, was ich da noch machen soll.
Wie bereits vorgeschlagen: Eine reduzierte Version des Programms zur Verfügung stellen (kleineres Gitter, weniger Zeitschritte), die in erträglicher Zeit (= wenige Minuten) durchläuft, damit man das analysieren und Verbesserungsvorschläge machen kann.
Grüße,
Harald
Einstellungen und Berechtigungen
Du kannst Beiträge in dieses Forum schreiben. Du kannst auf Beiträge in diesem Forum antworten. Du kannst deine Beiträge in diesem Forum nicht bearbeiten. Du kannst deine Beiträge in diesem Forum nicht löschen. Du kannst an Umfragen in diesem Forum nicht mitmachen. Du kannst Dateien in diesem Forum posten Du kannst Dateien in diesem Forum herunterladen
MATLAB, Simulink, Stateflow, Handle Graphics, Real-Time Workshop, SimBiology, SimHydraulics, SimEvents, and xPC TargetBox are registered trademarks and The MathWorks, the L-shaped membrane logo, and Embedded MATLAB are trademarks of The MathWorks, Inc.
) anders war, aber da täusche ich mich ja nun offensichtlich.
RSS