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

Stückweise definierte Funktion in f.-Datei plotten

 

Monika
Forum-Anfänger

Forum-Anfänger


Beiträge: 17
Anmeldedatum: 15.01.16
Wohnort: ---
Version: ---
     Beitrag Verfasst am: 31.03.2016, 12:19     Titel: Stückweise definierte Funktion in f.-Datei plotten
  Antworten mit Zitat      
Hallo zusammen!

Ich möchte gerne eine stückweise definierte Funktion in einer f.-Datei plotten bzw. eigentlich gerne im Workspace hinterlegen. Hier sind die notwendigen Files stark vereinfacht:

m.-File:

Code:

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

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  
f_1     =0.01;  
tr=1.000000000000000e-03;

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


f.-file:

Code:

function dydt = myode3(t,y,A,D,L_inv,f_1,tr)
 dydt = L_inv*A*E(t)+L_inv*D*y;
function E_val = E(t)
if     (t ==0)
    E_val=zeros(6,1) ;
elseif (t<tr)                          %mit tr= Zeit für die Rampe
    E=[1;0;0;0;0;0];
    E_val = E.*(t*(f_1/tr)) ;     %Geradengleichung für die Rampe
else
    E=[1;0;0;0;0;0];
    E_val=E.*f_1; %Konstante Dehnung f_1=0.01
end
end
end

 


Dabei geht es mir um die Funktion E_val. Bis zum Zeitpunkt tr gleicht die Funktion einer Gerade. Ab tr ist E_val konstant und beträgt den Wert 0.01. Das ganze wird für einen Startvektor E (6x1) berechnet. Demnach müsste E_val eine Matrix mit 6 Spalten und nt Zeilen sein. Die hätte ich gerne im Workspace gespeichert und ausgegeben.

Hat jemand eine Idee, wie das funktioniert bzw. was ich dafür ändern oder eingeben muss?

Vielen herzlichen Dank im voraus!

Beste Grüße,
Moni
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: 31.03.2016, 14:16     Titel: Re: Stückweise definierte Funktion in f.-Datei plotten
  Antworten mit Zitat      
Hallo Monika,

Eine nicht-differenzierbare Funktion mit ODE45 zu integrieren, ist eine wirklich schlechte Idee. Die Genauigkeit des Ergebnisses wird durch die Schrittweiten-Steuerung im Integrator kontrolliert. An den Knicken im Parameter E_val muss diese Steuerung jedoch scheitern. Entweder wird die Schrittweite auf das Minimum heruntergefahren und die Integration stoppt mit einem Fehler. Oder die Schrittweite wird nur so winzig, dass der Knick wegen der Rundungsfehler nicht mehr feststellbar ist. Das beeinträchtigt aber sowohl die Rechenzeit als auch die akkumulierten Rundungsfehler des Ergebnisses. Oder aber der Integrator "bemerkt" den Sprung nicht und hudelt einfach drüber, was die Ergebnisse ebenfalls deutlich verfälschen kann.
Du bekommst also entweder ein ungenaues Ergebnis oder gar keins. Für eine wissenschaftliche Berechnung ist es jedenfalls ungeeignet.

Der einzige sinnvolle Weg ist es, die Integration 3 mal mit den entsprechenden Parametern für E_val zu starten und den jeweiligen End-Wert für y als Startwert für die nächste Integration zu nutzen.

Eine Mutliplikation mit der Inversen einer Matrix ist aus numerischer Sicht immer eine schlechtere Lösung als L\A.

Was ist ein "f.-file"? Fortran? Der Code sieht aus wie eine Funktion in einem "M-File".

Die eigentliche Frage, wie Du die Matrix für E_val erhältst, löst man so:
Code:

function dydt = myode3(t,y,A,D,L_inv,f_1,tr)

if length(t) == 1
  dydt = L_inv * A * E(t) + L_inv * D * y;
else
  dydt = E(t);  % Reply E if t is a vector
end

  function E_val = E(t)   % Vectorized: Logical indexing instead of IF
    E_val = zeros(6, numel(t));
    index = (t > 0 & t < tr);
    E_val(1, index) = t(index) * (f_1 / tr);
    index = (t >= tr);     % EDITED: ">"  ==>  ">="
    E_val(1, index) = f_1;
  end
end

Nun kannst Du zuerst integrieren (in 3 Schritten, siehe oben!) und bekommst Deine Lösung für T und Y heraus. Mit diesem T kannst Du dann myode3 wieder aufrufen, und diesmal wird nicht dydt ausgegeben, sondern E. Statt der Länge des Inputs "t" könntest Du auch einen weiteren Parameter als Flag verwenden. Oder aber E extern berechnen mit der gleichen Methode.

Gruß, Jan

Zuletzt bearbeitet von Jan S am 31.03.2016, 18:44, insgesamt einmal bearbeitet
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: 31.03.2016, 14:50     Titel:
  Antworten mit Zitat      
Hallo Jan!

Vielen herzlichen Dank für diese ausführliche Antwort! Die Gleichung in 3 Schritten zu integrieren ist einleuchtend. Ich werde es probieren und melde mich, sobald ich Probleme habe oder auch wenn es geklappt hat. Mir ist nur nicht klar, was es mit der ersten if-Schleife auf sich hat..ich werde es mir aber zunächst genauer anschauen Smile

Eine Bitte hätte ich noch. Ich habe nun versucht anstatt einer verschachtelten Funktion, zwei m.-Files mit den beiden Funktionen für dtdy und für E(t) anzulegen. Das sieht dann so aus:

Code:

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

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  
f_1     =0.01;  
tr=1.000000000000000e-03;

yo = [0,0,0,0];                % y(t=1) = 1
[T Y] = ode45(@(t,y) myode4(t,y,E,A,D,L_inv,f_1,tr),te,yo); % Solve ODE
 


Code:

function dydt = myode4(t,y,A_ai,D_ab,L_ab_inv)
 dydt = L_ab_inv*A_ai*E(t)+L_ab_inv*D_ab*y;
end
 


Code:

function E_val = E(t)
if     (t ==0)
    E_val=zeros(6,1) ;
elseif (t<1.000000000000000e-03) %mit tr= Zeit für die Rampe
    E=[1;0;0;0;0;0];
    E_val = E.*(t*(0.01/1.000000000000000e-03)) ; %Geradengleichung für die Rampe
else
    E=[1;0;0;0;0;0];
    E_val=E.*0.01; %Konstante Dehnung f_1=0.01
end
end
 


Wenn man das ganz durchrechnen lässt und dann plot(E(T)) anzeigen lässt, sieht die Kurve anders aus als gewollt (siehe Datei im anhand). Gewollt war eine Rampe, die von 0 bis 1e-3 linear ansteigt bis zum Wert 0.01 und diesen dann beibehält. Hast du eine Idee, woran das liegt? Auch die Ergebnisse für y kommen mir genau falsch herum vor, sprich, als würde man die Zeit rückwärts abfahren :/

Vielen, vielen Dank nochmal!

Beste Grüße,
Monika

Plot.fig
 Beschreibung:

Download
 Dateiname:  Plot.fig
 Dateigröße:  11.45 KB
 Heruntergeladen:  295 mal
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: 31.03.2016, 15:48     Titel:
  Antworten mit Zitat      
Hallo noch einmal Jan!

Dein Code funktioniert bei mir, aber leider verstehe ich nicht alle Schritte Sad Was hat es mit der Länge von t in der if-Schleife auf sich? Ich wäre Dir überaus dankbar, wenn Du mir in wenigen Zeilen erklären könntest, was da genau passiert, weil ich so nicht in der Lage bin die nächsten 2 Schritte für die Integration anzugehen Sad

Zudem verwirrt mich (wie oben erwähnt), dass die abschnittsweise definierte Funktion trotz korrekter Definition seltsam aussieht..also wäre sie an der x-Achse gespiegelt worden. Liegt das an der Zeit-Skallierung? Ich bekomme mit deinem Code nämlich dieselben Ergebnisse für Y wie mit meinem, sodass auch Deine Funktion E(t) ähnlich komisch aussehen würde (kann ich leider nicht plotten, weil Deine Funktion in die Haupt-Funktion integriert ist).

Ich hoffe Du weißt was ich meine. Vielen Dank nochmal!

Beste Grüße,
Monika
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: 31.03.2016, 18:55     Titel:
  Antworten mit Zitat      
Hallo Monika,

Wie wäre es, alles in ein M-File zu packen, um es hier mal zu vereinfachen:
Code:

function myIntegration
A = rand(4,6);
D =rand(4,4);
L = rand(4,4);
L_inv = inv(L);    % "böse"

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  
f_1     =0.01;  
tr=1.0e-3;

yo = [0,0,0,0];                % y(t=1) = 1
[T, Y] = ode45(@(t,y) myode4(t,y,E,A,D,L_inv,f_1,tr),te,yo);
E_val = E(T, f_1, tr);

plot(T, E_val);
end

function dydt = myode4(t,y,A_ai,D_ab,L_ab_inv, f_1, tr)
 dydt = L_ab_inv*A_ai*E(t, f_1, tr)+L_ab_inv*D_ab*y;
end

function E_val = E(t, f_1, tr)
    E_val = zeros(6, numel(t));
    index = (t > 0 & t < tr);
    E_val(1, index) = t(index) * (f_1 / tr);
    index = (t >= tr);     % EDITED: ">"  ==>  ">="
    E_val(1, index) = f_1;
end

Wenn man E(T) gegen T aufträgt, sieht das nach einer Rampe aus, ganz genau wie Du sie schilderst.

In meinem Code-Beispiel wird myode4 vonm Integrator mit einem skalaren Werte für t aufgerufen. Beim Aufruf mit dem Vektor T wird dann nicht dydt zurückgegeben, sondern der Wert von E(T). Mit einer eigenständigen Funktion ist das aber hübscher, so wie im Code hier.

Bemerkung: "IF" ergibt keine "Schleife". Nur WHILE und FOR erzeugen Schleifen.

Gruß, Jan
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: 01.04.2016, 09:03     Titel:
  Antworten mit Zitat      
Hallo Jan,

ich habe Deinen Code übernommen (an einigen Stellen mussten A_ai, D_ab und L_ab_inv in A,D und L_inv geändert werden) und die Datei als "myIntegration" gespeichert. Wenn ich alles anpasse, läuft der Code trotzdem nicht durch. Das hier ist eine der Fehlermeldungen:

Not enough input arguments.

Error in myIntegration>E (line 28 )
E_val = zeros(6, numel(t));

Liegt es daran, dass die erste Funktion "myIntegration" nicht zugeordnet ist?

Vielen Dank,
Monika
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: 01.04.2016, 18:38     Titel:
  Antworten mit Zitat      
Hallo Monika,

Das ist eine sehr seltsame Fehlermeldung!
Entweder "zeros" oder "numel" muss hier mit einer von Dir definierten Funkltion überschrieben worden sein, die mehr Inputs erwartet.

Setze einen Breakpoint in die Zeile, die scheitert, und untersuche die verwendeten Symbole:
Code:


Gruß, Jan
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.