WICHTIG: Der Betrieb von goMatlab.de wird privat finanziert fortgesetzt. - Mehr Infos...

Mein MATLAB Forum - goMatlab.de

Mein MATLAB Forum

 
Gast > Registrieren       Autologin?   

Partner:




Forum
      Option
[Erweitert]
  • Diese Seite per Mail weiterempfehlen
     


Gehe zu:  
Neues Thema eröffnen Neue Antwort erstellen

Mit ode45 DGL 1.Ordnung in Matrixform lösen

 

Monika
Forum-Anfänger

Forum-Anfänger


Beiträge: 17
Anmeldedatum: 15.01.16
Wohnort: ---
Version: ---
     Beitrag Verfasst am: 18.02.2016, 12:29     Titel: Mit ode45 DGL 1.Ordnung in Matrixform lösen
  Antworten mit Zitat      
Hallo nochmal!

Ich probiere es noch einmal in einem neuen Thread, da ich kurz vorm Verzweifeln bin Sad

Ich möchte gerne eine DGL 1. Ordnung, die Matritzen unterschiedlicher Dimensionen enthält, mithilfe des Befehls ode45 lösen. Die Differentialgleichung lautet wie folgt:

dy/dt=(inv(L)*A*E(t))' + inv(L)*D*(y(t))'

Informationen zu den Matritzen:

1. dy/dt ist ein 4x1 Vektor, wobei für jeden dieser 4 Einträge ydot1, ydot2, ydot3 und ydot4 die DGL gelöst werden soll. Als Ergebnis soll dann eine nt x4 Matrix rauskommen, sodass die erste Zeile beispielsweise die Ergebnisse aller ydot(1-4) zum Zeitounkt=0 enthält.
2. L ist eine 4x4 Matrix und konstant
3. A ist eine 4x6 Matrix und konstant
4. D ist eine 4x4 Matrix und konstant
5. E(t) ist eine 6x nt Matrix, wobei nt für die Zeitschritte steht


Ich habe nun mehrere Fragen:

1. Die Matrix E ist in sofern zeitabhängig, dass jede Spalte dieses Vektors einem Zeitpunkt zugeordnet ist. Deswegen hat die Matrix auch nt Spalten. Wie baue ich diese Information in die ode45 ein? Es müsste also für jeden Zeitschritt die entsprechende Spalte verrechnet werden.

2. Ich habe zudem gelesen, dass für Y ein Spaltenvektor rauskommen muss, wobei jede Zeile einem Zeitschritt zugeordnet ist. Wenn ich die entsprechenden Matritzen transponiere, wie oben, dann bekomme ich für die Ergebnis-Matrix die Form nt x4 (ohne Transponieren, wäre die Matrix um 90° gedreht und dann würde nicht jede Zeile einem Zeitschritt entsprechen, wie es gefordert ist). Wie kann ich diese Ergebnis-Matrix vorher so als Array definieren, dass einzelne Vektoren abgespeichert werden, aber eben vier Stück für ydot(1-4)?

3: Überhaupt habe ich das Problem mit den Dimensionen. Welche Dimension wird für y(t) angenommen und kann ich das steuern? Sie müsste in irgendeiner Form die Dimensionen nt x4 bzw. nach dem Transponieren 4x nt haben.

4. Welche Dimension hat dann die Anfangsbedingung y0? wäre dann y0=[0,0,0,0] für die 4 zu berechnenden Variablen?

Hier nocheinmal die Codes. Zunächst das Hauptskript:

Code:

A=rand(4,6);
D=rand(4,4);
L=rand(4,4);

dt      = 1e-5;             %time step size
tstart  = 0;                %starting time
tend    = 1e-1;             %total time
tr      = 1e-3;             %loading time (ramp)
te       =[tstart:dt:tend]';
nt      = length(te);    %number time steps    

E=zeros(6,nt);
E(1,:) = te.*(0.01/tr); % Generate f(t)
E(1,tr/dt+1:end)=0.01;

yo = [0,0,0,0]; % y(t=1) = 1
[T Y] = ode45(@(t,y) myode2(t,y,E,A,D,L),te,yo); % Solve ODE
 


Und diem.-Datei:
Code:

function dydt = myode2(t,y,E,A,D,L)
dydt = (inv(L)*A*E)'+inv(L)*D*(y)'; % Evaluate ODE at time t
 


Ich bin für jede Hilfe sehr dankbar, weil ich wirklich nicht weiter komme.

Herzlichen Dank und viele Grüße,
Monika
Private Nachricht senden Benutzer-Profile anzeigen


Winkow
Moderator

Moderator



Beiträge: 3.842
Anmeldedatum: 04.11.11
Wohnort: Dresden
Version: R2014a 2015a
     Beitrag Verfasst am: 18.02.2016, 12:41     Titel:
  Antworten mit Zitat      
was soll denn dieses E sein? ist das eine zeitabhängigkeit? wie man die benutzt steht doch in der doc unter example 3 ich würd mich daran halten. oder soll das was anderes sein? E sollte also 6x1 sein und y0 4x1. irgendwie finde ich auch komisch das 5 einträge von E dann 0 sind. da kann man bestimmt irgendwo auch was einspaaren dann.
_________________

richtig Fragen
Private Nachricht senden Benutzer-Profile anzeigen
 
Monika
Themenstarter

Forum-Anfänger

Forum-Anfänger


Beiträge: 17
Anmeldedatum: 15.01.16
Wohnort: ---
Version: ---
     Beitrag Verfasst am: 18.02.2016, 12:53     Titel:
  Antworten mit Zitat      
Wenn man E plottet zeigt es eine zunächst linear ansteigende Rampe. Ab einem bestimmten Zeitpunkt ist E dann konstant. Dass die anderen Einträge gleich Null sind ist nicht immer der Fall, deswegen kann E auch in keinem Fall vereinfacht werden.

Und wie bereits gesagt: Die 6 Einträge pro Spalte zeigen einen Zustand zu einem bestimmten Zeitschritt nt.

Und ich habe mir schon unzähliges dazu durchgelsen, aber ich verstehe es nicht. Deswegen stelle ich die Frage hier ja.
Private Nachricht senden Benutzer-Profile anzeigen
 
Winkow
Moderator

Moderator



Beiträge: 3.842
Anmeldedatum: 04.11.11
Wohnort: Dresden
Version: R2014a 2015a
     Beitrag Verfasst am: 18.02.2016, 12:57     Titel:
  Antworten mit Zitat      
wie ich bereits sagte. einfach an die doc halten
Code:
A=rand(4,6);
D=rand(4,4);
L=rand(4,4);

dt      = 1e-5;             %time step size
tstart  = 0;                %starting time
tend    = 1e-1;             %total time
tr      = 1e-3;             %loading time (ramp)
te       =[tstart:dt:tend];
nt      = length(te);    %number time steps    

E=zeros(6,nt);
E(1,:) = te.*(0.01/tr); % Generate f(t)
E(1,tr/dt+1:end)=0.01;

yo = [0,0,0,0]'; % y(t=1) = 1
[T Y] = ode45(@(t,y) myode2(t,y,E,A,D,L,te),te,yo); % Solve ODE

Code:
function dydt = myode2(t,y,E,A,D,L,te)
E1=interp1(te,E',t)';
dydt = (inv(L)*A*E1)+inv(L)*D*(y); % Evaluate ODE at time t
end

_________________

richtig Fragen
Private Nachricht senden Benutzer-Profile anzeigen
 
Monika
Themenstarter

Forum-Anfänger

Forum-Anfänger


Beiträge: 17
Anmeldedatum: 15.01.16
Wohnort: ---
Version: ---
     Beitrag Verfasst am: 18.02.2016, 13:11     Titel:
  Antworten mit Zitat      
Ok, es ist schonmal super, dass es durchläuft, also vielen herzlichen Dank dafür!

Eine Frage habe ich aber noch: ist es richtig E zu interpolieren bzw. gibt das nicht ein falsches Ergebnis? E liefert nur bis zum Teitpunkt tr linear ansteigende Werte, ab dann beträgt der Wert immer nur 0.01. Im Prinzip sind das ja dann zwei Funktionen, wenn man so will. Kann E dann als eine einzige Funktion interpoliert werden?

Kann ich das irgendwie mit einer if-Schleife lösen?
Beispielsweise so in der Art:

Code:

function dydt = myode3(tr,t,y,E,A,D,L,te)
if t<tr
E=interp1(te,E',t)';
else
    E(1,:)=0.01;
end
dydt = (inv(L)*A*E)+inv(L)*D*(y); % Evaluate ODE at time t
end
 


Der Code funktioniert so leider nicht, aber ich hoffe ihr wisst, was ich meine.

Gruß,
Monika
Private Nachricht senden Benutzer-Profile anzeigen
 
Winkow
Moderator

Moderator



Beiträge: 3.842
Anmeldedatum: 04.11.11
Wohnort: Dresden
Version: R2014a 2015a
     Beitrag Verfasst am: 18.02.2016, 14:05     Titel:
  Antworten mit Zitat      
Zitat:
: ist es richtig E zu interpolieren bzw. gibt das nicht ein falsches Ergebnis? E liefert nur bis zum Teitpunkt tr linear ansteigende

solange der anstieg nicht unendlich ist und ist es doch egal für die interpolation.
du kannt auch eine funktion da hin schreiben für E
_________________

richtig Fragen
Private Nachricht senden Benutzer-Profile anzeigen
 
Neues Thema eröffnen Neue Antwort erstellen



Einstellungen und Berechtigungen
Beiträge der letzten Zeit anzeigen:

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
.





 Impressum  | Nutzungsbedingungen  | Datenschutz | FAQ | goMatlab RSS Button RSS

Hosted by:


Copyright © 2007 - 2025 goMatlab.de | Dies ist keine offizielle Website der Firma The Mathworks

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.