Verfasst am: 03.11.2012, 14:57
Titel: Zeitabhängige DGL mit ODE lösen (mit Daten in Matrizen)
Hallo,
ich wollte Euch bitte um Hilfe bitten, da ich nicht weiterkomme. Ich habe bereits ODE45 Bsp 3 angeschaut, nur hilft mir dass nicht so ganz für mein Problem bzw. ich weiß nicht wie ich das auf mein Problem ummünzen kann.
Beispiel Formel:
Und zwar habe ich in etwa sowas gegeben:
dy/dt = M(t) + y(t) * C
Erklärung:
- C ist eine Konstante (Skalar)
- M(t) habe ich als [Nx3] (N Zeilen, 3 Spalten) vorberechnet. N ca 4000 Werte. Es sind sozusagen X, Y, Z Koordinaten darin für jeden Zeitpunkt t.
- Als y(t) soll immer das zuvor berechnete y(t-1) genommen werden, wobei y(0) eine [1x3] Vektor mit Startwerten ist (X, Y, Z Komponente).
Meine Frage:
Das Problem ist nun, dass ich nicht weiß wie ich in der Funktion die ich für ODE schreiben muss, die Werte von M für jeden Zeitpunkt extra bekomme.
Ich könnte ja schreiben M(1, 2) * M(1,3), dass bedeutet aber ich bekommen für jeden Zeitpunkt immer nur die erste Zeile. Ich bräuchte aber für jeden Schritt einen Index (von t = 1 bis N): M(1,2) * M(1,3), nächster Schritt: M(2,2) * M(2,3)
Also kurz gesagt ich will (glaub ich) sowas in der Art:
Soviel es Zeilen in M gibt (sind ca 4000). Wie kann ich den ersten Index in M immer ändern?
Hier eine vereinfachte Version meiner Funktion, die ich mit ODE lösen will. Ich will sozusagen jedes gelöste dydt in einer Liste speichern und es für die neue Berechnung zum neuen Zeitpunkt verwenden.
Der erste Index von M und y_init soll bei jedem Schritt um 1 erhöt werden. Wie mache ich das? Oder macht das Matlab schon irgendwie? Wenn ja, wie soll ich M anschreiben? Muss ich noch irgendwas mitgeben dieser Funktion damit für jeden Aufruf der Index 1, 2, 3, ... mitgegeben und automatisch erhöht wird?
Ich hoffe Ihr könnt mir helfen. Ich habe leider noch nie mit ODE gearbeitet :-/.
Vielen Danke schon einmal!
MfG
[EDIT]
Oder ist dies nur möglich, indem ich in einer Schleife immer eine Zeile von M (M(counter,: )) und y(counter) = ode45(fun ... y(counter-1)) mitgebe?
was das Arbeiten mit M angeht, würde ich mir durch Interpolation die richtigen Werte holen, wie es ja auch in dem von dir refereenzierten Beispiel gemacht wird.
Wenn das möglich ist, fände ich es aber sinnvoller, M in jedem Aufruf separat zu berechnen.
Zitat:
Als y(t) soll immer das zuvor berechnete y(t-1) genommen werden
Da ode45 mit variabler Schrittweite rechnet, halte ich die Idee mathematisch nicht für sinnvoll. y(t) steht aber ohnehin in der DGL-Funktion zur Verfügung. Der Aufwand ist also nicht notwendig.
Ja, ich verstehe denn Sinn der Interpolation hier nicht ganz.
Grundlegend geht es darum, die Erdrotation mit äußeren Einflüssen über die Eulersche Gleichung zu lösen und diese ist für äußere Einflüsse nur numerisch lösbar.
In Wiki-Artikel, die erste Formel ist jene, welche ich lösen will.
I sind einfach konstanten und M ist von mir berechneter Einfluss zum Zeitpunkt t mit gegebene X,Y,Z Koordinaten von anderen Planeten. Vereinfacht berechnet sich M wie folgt:
ode45 ist hier gut geeignet. Du kannst aber doch M wunderbar zu jedem gewünschten Zeitpunkt berechnen? Falls du M nur zu bestimmten Zeitpunkten berechnen kannst, ist das der Moment, in dem du Interpolation brauchst.
Sind x_t, y_t und z_t denn als Funktionen oder Wertetabellen gegeben?
Ja stimmt, ich könnte M immer berechnen, jedoch ist es einfacher, sofort mit Matrizen die kompletten Werte für M vorzuberechnen und sie dann einfach in ode45 zu verwenden (denke ich).
x, y und z sind Positionen welche zu jeweils 10 min über ein Jahr gesampelt wurden (stehen in einer .txt Datei welche ich einlese - sind sehr viele Werte). Diese sind also vorgegeben.
Ich muss nur mit diesen Werten und jeweils den vorgegangen y_t (bzw. w(t), ein neues y_t+1 (bzw. w(t+1)) berechnen und immer speichern.
Das sollte, wie du schon sagtest, sehr gut mit ode45 Funktion zu lösen sein.
Nur ich weiß nicht wie ich das in Matlabschreibweise schreiben soll . Sry für meine umständliche Erklärung für mein Problem .
Ich muss dann alle y_t plotten (wobei dies dann 3 Dimensional ist, da y_t ein 3D Vektor ist).
Meine Idee ist nun Runge-Kutte 4ter Ordnung mit ode45 zu implementieren um die jeweiligen nächsten y_t+1 (bzw. w(t+1)) zu berechnen. Als Zeitintervallen nehme ich 10min (bzw. 10*60s), welche genau zu x,y,z passen, meiner Meinung nach.
x, y und z sind Positionen welche zu jeweils 10 min über ein Jahr gesampelt wurden (stehen in einer .txt Datei welche ich einlese - sind sehr viele Werte). Diese sind also vorgegeben.
Ok. Dann ist es sinnvoll, die M-Werte auf einmal zu berechnen und zu interpolieren. Interpolation ist entscheidend, da ode45 Schrittweitensteuerung hat (Schritte können deutlich größer oder kleiner als 10 Min. sein) und es wichtig ist, dass man die richtigen Werte aus M extrahiert.
Zitat:
Ich muss nur mit diesen Werten und jeweils den vorgegangen y_t (bzw. w(t), ein neues y_t+1 (bzw. w(t+1)) berechnen und immer speichern.
Dieser Ansatz ist nur sinnvoll, wenn du mit einem RK-Verfahren mit fester Schrittweite arbeitest. Das ist aber bei ode45 nicht der Fall. Bitte denke daran, dass du an ode45 eine Funktion der Form
ok das klingt logisch. So, ich habe jetzt ein wenig probiert, ich weiß aber nicht wie ich M interpolieren kann und welches Zeitintervall ich bei ode45 dann angeben muss.
Ich habs jetzt vereinfacht, so dass M immer 0 ist. (Das würde bedeuten, Erdrotation ohne äußere Einflüsse).
Hier einmal der Code:
Code:
function dw = getWdot(t, w, A, B, C, M)
dw = [M(1)/A - (C-B)/A* w(3) * w(2); M(2)/B - (A-C)/B * w(3) * w(1); M(3)/C];
So, die einzige und wichtigste Frage ist nun eben, wie kann ich meine M Matrix in die Funktion getWdot bringen und so dass sie Zeitinterpoliert wird.
M-Matrix mit ~50000 3-Dimensionalen Werten (aufgenommen in 10Minuten Intervall, das entspricht ungefähr einem Jahr).
(Dazu hätte ich für jeden aufgenommen X,Y,Z Wert die "Uhrzeit" (immer im 10Min Takt), falls ich das für das Zeitintervall wichtig sein sollte?).
Wie interpoliere ich 3 Dimensionale Werte in Matlab und wie welchen Zeitintervall soll ich ode mitgeben? Zurzeit habe ich ein Jahr (365*24*60*60) in Sekunden mitgeben.
Vielen Danke für die ganze Hilfe. Echt super dieses Forum !
ja 365 stimmt (wenn ich länger vor dem PC sitze, vertuh ich mich schnell) . Ok, ich dachte interp1 nimmt nur 2D Vektoren. Dann hab ich es wohl falsch verwendet weil es nie geklappt hatte. Dachte man müsste mit bei 3D Koordinaten mit Splines evtl. arbeiten aufgrund besserer Ergebnisse, etc.
Vielen Dank Harald, werde es sofort einmal probieren!
ich gehe mal davon aus, dass die Datenpunkte verhältnismäßig nah beieinanderliegen (d.h. keine großen Schwankungen und Oszillationen oder so; ein glatter Verlauf, wenn du die Datenpunkte selbst einfach mal plottest). Dann dürfte es relativ egal sein, wie man interpoliert. interp1 hat aber auch eine Option 'spline'.
Wenn die Punkte nicht im obigen Sinne nah beieinander liegen, ist es für mich fraglich, ob man überhaupt vernünftige Ergebnisse erwarten kann.
Ok, habs nun probiert. Es sollte theoretisch funktionieren . Habs mit spline und ohne Parameter probiert.
Ich erhalte zwar nicht mal ansatzweise mein gewünschtes Ergebnis, aber dass ist glaube ich ein anderes Problem welches ich noch nicht verstanden/gefunden habe
ohne weitere Informationen werde ich da nicht weiterhelfen können.
Zitat:
Ich erhalte zwar nicht mal ansatzweise mein gewünschtes Ergebnis,
Was ist das "gewünschte Ergebnis", und was bekommst du stattdessen?
Am besten nachvollziehbar wäre das anhand von Beispielcode (inkl. Beispieldaten, es müssen ja nicht die Daten für das ganze Jahr sein).
eingeschränkt und jetzt sieht das Ergebnis schon viel besser aus. Jetzt kenn ich mich langsam mit ODE aus, war aber ein sehr schwerer Anfang . Das mit den Interpolationen war super! Vielen Dank für Deine Hilfe!
statt 'MaxStep' würde ich eher 'RelTol' und/oder 'AbsTol' herunterzusetzen. Die berücksichtigen allerdings noch nicht die Interpolationsfehler.
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.
. Sry für meine umständliche Erklärung für mein Problem 
.
!!
. Das mit den Interpolationen war super! Vielen Dank für Deine Hilfe! 
RSS