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:
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?
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:
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
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
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:
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 :/
Dein Code funktioniert bei mir, aber leider verstehe ich nicht alle Schritte 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
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!
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.
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?
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:
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.