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

Lösen einer DGL mit zeitabhängiger Einwirkung

 

eisklecks
Forum-Newbie

Forum-Newbie


Beiträge: 4
Anmeldedatum: 22.11.15
Wohnort: ---
Version: ---
     Beitrag Verfasst am: 22.11.2015, 22:47     Titel: Lösen einer DGL mit zeitabhängiger Einwirkung
  Antworten mit Zitat      
Hallo zusammen.

Im Rahmen meiner Bachelorarbeit muss ich folgende Bewegungsgleichung mit Matlab lösen:

M*x''(t)+D*x'(t)+K*x(t) = F(t)

Es handelt sich hierbei um die Berechnung eines 3 stöckigen Hauses z.B. unter Erdbebeneinwirkung. M, D, K (je 3x3 Matrizen) und F(t) (3x1 Matrix) sind bekannt.

Bisher bin ich zum Lösen einer Differentialgleichung nur auf den ODE-Befehl gestoßen, jedoch funktioniert dieser soweit ich weiß nicht, wenn in der DGL kein zeitabhängiger Parameter vorhanden ist. (bzw. ich habe noch keine Lösung zu meinem Problem gefunden).
Mein Betreuer meinte, man könnte die Bewegungsgleichung auch über das Zustandsraum-Modell lösen, das wollen wir aber vermeiden, wenn es eine Alternative gibt.

Kann mir einer von euch weiterhelfen? Da ich Matlab im Zuge meiner BA zum ersten Mal benutzt habe, kenne ich mich noch nicht soo gut aus.

Vielen Dank schonmal für jeden Tipp! Smile
Private Nachricht senden Benutzer-Profile anzeigen


Harald
Forum-Meister

Forum-Meister


Beiträge: 24.448
Anmeldedatum: 26.03.09
Wohnort: Nähe München
Version: ab 2017b
     Beitrag Verfasst am: 22.11.2015, 23:34     Titel:
  Antworten mit Zitat      
Hallo,

sollte an sich kein Problem sein.
Erster Schritt wird sein, das System dreier DGLen zweiter Ordnung in ein System von sechs DGLen erster Ordnung überzuführen. Dann sollte es z.B. mit ode45 funktionieren.

Grüße,
Harald
Private Nachricht senden Benutzer-Profile anzeigen
 
Jan S
Moderator

Moderator


Beiträge: 11.057
Anmeldedatum: 08.07.10
Wohnort: Heidelberg
Version: 2009a, 2016b
     Beitrag Verfasst am: 23.11.2015, 10:26     Titel: Re: Lösen einer DGL mit zeitabhängiger Einwirkung
  Antworten mit Zitat      
Hallo eisklecks,

Zitat:
Bisher bin ich zum Lösen einer Differentialgleichung nur auf den ODE-Befehl gestoßen, jedoch funktioniert dieser soweit ich weiß nicht, wenn in der DGL kein zeitabhängiger Parameter vorhanden ist.

Das "t" in Deiner Gleichung sieht doch sehr nach "Zeit" aus.
In der Dokumentation von ode45 bzw. ode15i wird genau erklärt, wie die Befehle funktionieren. Hilft Dir das bereits weiter?

Gruß, Jan
Private Nachricht senden Benutzer-Profile anzeigen
 
eisklecks
Themenstarter

Forum-Newbie

Forum-Newbie


Beiträge: 4
Anmeldedatum: 22.11.15
Wohnort: ---
Version: ---
     Beitrag Verfasst am: 23.11.2015, 12:43     Titel:
  Antworten mit Zitat      
Hallo, danke schonmal für eure Antworten.

Zitat:
Das "t" in Deiner Gleichung sieht doch sehr nach "Zeit" aus.
In der Dokumentation von ode45 bzw. ode15i wird genau erklärt, wie die Befehle funktionieren. Hilft Dir das bereits weiter?


Wie der ode-Befehl funktioniert habe ich soweit verstanden. Allerdings weiß ich momentan noch nicht wie ich F(t) einbinde, da dieses ja nicht konstant ist, sondern nach jeder Sekunde einen anderen Wert annimmt:
F(1)= -13.3980
F(2)= -103.3560
F(3)= -96.6570
F(4)= -84.2160
F(5)= -90.9150
.... (die entsprechenden Werte ziehe ich mir aus einer anderen mat-Datei!)

Wenn ich versuche mein F(t) (in der Form) in den ode-Befehl einzufügen, bekomme ich immer eine Fehlermeldung.


Zitat:
Erster Schritt wird sein, das System dreier DGLen zweiter Ordnung in ein System von sechs DGLen erster Ordnung überzuführen. Dann sollte es z.B. mit ode45 funktionieren.

Wie genau meinst du das mit der Überführung?
Private Nachricht senden Benutzer-Profile anzeigen
 
Jan S
Moderator

Moderator


Beiträge: 11.057
Anmeldedatum: 08.07.10
Wohnort: Heidelberg
Version: 2009a, 2016b
     Beitrag Verfasst am: 23.11.2015, 14:23     Titel:
  Antworten mit Zitat      
Hallo eisklecks,

Es gibt keinen "ode-Befehl" in Matlab. Es kommt darauf an, welchen Integrator Du für das System verwenden möchtest und dies wird zu verschiedenen Ergebnissen führen. Also meinst Du ode45 , ode15i , ode23 , ...? Weitere Integratoren findest Du unter http://www.mathworks.com/help/matlab/ref/ode45.html

Wenn F unstetig ist, können Matlabs Integratoren gar nicht damit umgehen. Die einzige numerisch zuverlässige Behandlung ist dann den Integrator nur von 0 bis 1, von 1 bis 2 etc laufen zu lassen.
Code:
F = [-13.3980,  -103.3560, -96.6570, -84.2160, -90.9150];
t0 = 0;
x0 = <startwert>
tC = cell(1, 5);
xC = cell(1, 5);
for k = 1:5
  fcn = @(t, x) myFcn(t, x, F(k));
  [t, x] = ode45(fcn, [t0, k], x0);
  tC{k} = t;
  xC{k} = x;
  x0 = x(end);
  t0 = t(end);
end

x = cat(1, xC{:});
t = cat(1, tC{:});
 


Zitat:
Wenn ich versuche mein F(t) (in der Form) in den ode-Befehl einzufügen, bekomme ich immer eine Fehlermeldung.

Bitte poste für eine Diskussion im Forum dann immer den Code und die vollständige Fehlermeldung. Wir wollen wirklich gerne helfen, können aber die Details nicht erraten.

Zitat:
Wie genau meinst du das mit der Überführung?

ODE45 kann wie die meisten anderen ODE-Integratoren nur Systeme erster Ordnung integrieren. Du kannst im Netz oder Deinem Vorlesungsskript danach suchen, wie man ein Sytsem n.ter Ordnung in n Systeme 1. Ordnung umwandelt.

Gruß, Jan
Private Nachricht senden Benutzer-Profile anzeigen
 
eisklecks
Themenstarter

Forum-Newbie

Forum-Newbie


Beiträge: 4
Anmeldedatum: 22.11.15
Wohnort: ---
Version: ---
     Beitrag Verfasst am: 25.11.2015, 17:49     Titel:
  Antworten mit Zitat      
Zitat:
Wenn F unstetig ist, können Matlabs Integratoren gar nicht damit umgehen. Die einzige numerisch zuverlässige Behandlung ist dann den Integrator nur von 0 bis 1, von 1 bis 2 etc laufen zu lassen.


So etwas in der Art hatte ich auch schonmal überlegt. Danke!
Ich habe allerdings mit meinem Betreuer besprochen, dass unser F doch erstmal wie folgt definiert werden soll: F(t) = a*sin(omega*t)

Wie ich die Differentialgleichung dann mit ode23 löse, wird ja sehr gut in Beispiel 3 dargestellt.

Mein nächstes Problem ist jetzt, dass ich mir nicht sicher bin, wie ich die Differentialgleichung mit den Matrizen gelöst bekomme.
Ich habe meine Funktion erstmal wie folgt geschrieben: (Die jeweiligen Werte aus der Matrix habe ich per Hand eingegeben und nicht zB mit K(1,1) abgefragt!)

Code:
function dzdt = rode2(z,t,F)  % M*x''(t) + D*x'(t) + K*x(t) = F

K = [ 1052702, -526351, 0;
      -526351, 1052702, -526351;
            0, -526351, 526351    ]; %Steifigkeitsmatrix
       
M = [ 957,   0,   0;
        0, 957,   0;
        0,   0, 1040  ]; %Massenmatrix
   
M_inv = [0.0010,      0, 0;
              0, 0.0010, 0;
              0,      0, 0.0010 ]; %Inverse der Massenmatrix
   
D = [ -2339, -5459, -5932;
       -798, -1206, -8534;
       3527,  7055, 14890 ]; %Dämpfungsmatrix
 
F = 120*sin(t*3);
 
dzdt = zeros(6,1);  
dzdt(1) = z(2);
dzdt(2) = z(4);
dzdt(3) = z(6);
dzdt(4) = [F-0.0010*(-2339*z(2)-5459*z(4)-5932*z(6)+1052702*z(1)-526351*z(3))];
dzdt(5) = [F-0.0010*(-798*z(2)-1206*z(4)-8534*z(6)-526351*z(1)+1052702*z(3)-526351*z(5))];
dzdt(6) = [F-0.0010*(3527*z(2)+7055*z(4)+14890*z(6)-526351*z(3)+526351*z(5))];


und diese dann mit dem gleichen Aufruf wie in Beispiel 3 versucht zu lösen:


Code:
Tspan = [1 100];
IC = [1, 1, 1, 0, 0, 0];
[T Z] = ode45(@(t,z) rode2(t,z,Ft,F,K,M,D,M_inv),Tspan,IC);


Ich erhalte nun folgende Fehlermeldung:
Code:
Attempted to access z(2); index out of bounds because numel(z)=1.

Error in rode2 (line 22)
dzdt(1) = z(2);

Error in @(t,z)rode2(t,z,F)


Error in odearguments (line 87)
f0 = feval(ode,t0,y0,args{:});   % ODE15I sets args{1} to yp0.

Error in ode23 (line 112)
[neq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0, odeArgs, odeFcn, ...


Kann mir einer helfen, wo der Fehler ist! Ist meine Herangehensweise bezüglich der Funktionsdefinition überhaupt richtig?
Private Nachricht senden Benutzer-Profile anzeigen
 
Harald
Forum-Meister

Forum-Meister


Beiträge: 24.448
Anmeldedatum: 26.03.09
Wohnort: Nähe München
Version: ab 2017b
     Beitrag Verfasst am: 25.11.2015, 21:15     Titel:
  Antworten mit Zitat      
Hallo,

ein Problem ist auf jeden Fall: dein anonymous function handle ruft rode2 mit rode2(t,z,Ft,F,K,M,D,M_inv) auf. rode2 nimmt aber nur drei Argumente entgegen, und diese auch noch in einer anderen Reihenfolge.
Ich weiß nicht, ob das auch noch ein Problem macht, aber ich würde konsequent bei Zeilenvektoren (momentan IC) oder Spaltenvektoren (momentan dzdt) bleiben.

Solchen Problemen kann man übrigens auch gut mit dem Debugger auf den Grund gehen.

Grüße,
Harald
Private Nachricht senden Benutzer-Profile anzeigen
 
eisklecks
Themenstarter

Forum-Newbie

Forum-Newbie


Beiträge: 4
Anmeldedatum: 22.11.15
Wohnort: ---
Version: ---
     Beitrag Verfasst am: 27.11.2015, 17:07     Titel:
  Antworten mit Zitat      
Vielen Dank!

Leider stehe ich gerade wieder vor einer neuen Fehlermeldung, da ich doch unseren ursprünglichen Plan verfolgen soll, dass bei jedem Zeitschritt eine andere äußere Einwirkung F(t) in die DGL eingelesen werden soll.
Die verschiedenen Einwirkungen wurden in einem Zeitschritt von 0.01s gemessen und importiere ich aus einer Workspace-Datei. Meine DGL soll in diesen Zeitschritten jeweils einmal ausgewertet werden. Mein Befehl sieht wie folgt aus:

Code:
Tspan = [0:0.01:10];
IC = [0 0 0 0 0 0];
[T Z] = ode45(@(t,z) rode3(t,z,j),Tspan,IC);
plot(T,Z(:,2),T,Z(:,4),T,Z(:,6))
legend('x1','x2','x3')


Meine function rode3 ist folgendermaßen definiert:

Code:
function dzdt = rode3(t,z,j)          %DGL: M*x''(t) + D*x'(t) + K*x(t) = F(t)

K = [ 1052702, -526351, 0;              %Steifigkeitsmatrix
      -526351, 1052702, -526351;
            0, -526351, 526351    ];    
       
M = [ 957,   0,   0;                    %Massenmatrix
        0, 957,   0;
        0,   0, 1040  ];                
   
M_inv = [0.0010,      0, 0;             %Inverse der Massenmatrix
              0, 0.0010, 0;
              0,      0, 0.0010 ];      
   
D = [ -2339, -5459, -5932;              %Dämpfungsmatrix
       -798, -1206, -8534;
       3527,  7055, 14890 ];            
   
load ('matlab')                         %Einwirkung
F = zeros(3000,1);
for j=1:3000
    F(j) = a_g_elcentro(j);
end
 
dzdt = zeros(6,1);  
dzdt(1) = z(2);
dzdt(2) = [F(t*100+1)-M_inv(1,1)*(D(1,1)*z(2)+D(1,2)*z(4)+D(1,3)*z(6)+K(1,1)*z(1)+K(1,2)*z(3))];
dzdt(3) = z(4);
dzdt(4) = [-M_inv(2,2)*(D(2,1)*z(2)+D(2,2)*z(4)+D(2,3)*z(6)+K(2,1)*z(1)+K(2,2)*z(3)+K(2,3)*z(5))];
dzdt(5) = z(6);
dzdt(6) = [-M_inv(3,3)*(D(3,1)*z(2)+D(3,2)*z(4)+D(3,3)*z(6)+K(3,2)*z(3)+K(3,3)*z(5))];


Nun erhalte ich folgende Fehlermeldung:
Code:
Attempted to access F(1.2); index must be a positive integer or logical.

Error in rode3 (line 27)
dzdt(2) = [F(t*100+1)-M_inv(1,1)*(D(1,1)*z(2)+D(1,2)*z(4)+D(1,3)*z(6)+K(1,1)*z(1)+K(1,2)*z(3))];


Ich weiß nicht wo der Fehler ist bzw wieso überhaupt versucht wird die Zeile 1.2 meines Vektors F einzulesen. Meine Schrittlänge bei der Integration hatte ich ja mit 0.01 definiert und es wird jeder F-Wert an der Stelle "t*100+1" abgefragt. Durch die Umrechnung "t*100+1" entstehen doch immer nur ganze Zahlen.

Kann mir nochmal jemand helfen??
Private Nachricht senden Benutzer-Profile anzeigen
 
Harald
Forum-Meister

Forum-Meister


Beiträge: 24.448
Anmeldedatum: 26.03.09
Wohnort: Nähe München
Version: ab 2017b
     Beitrag Verfasst am: 27.11.2015, 19:19     Titel:
  Antworten mit Zitat      
Hallo,

verwende doch den Debugger.

Zitat:
Durch die Umrechnung "t*100+1" entstehen doch immer nur ganze Zahlen.

Wenn du den Debugger verwendest, wirst du feststellen, dass dem eben nicht so ist, da ode45 Zwischenschritte macht um die Lösung mit der erforderlichen Genauigkeit zu berechnen. Wie wäre es mit Interpolation?

Grüße,
Harald
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 - 2024 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.