Verfasst am: 26.05.2010, 22:46
Titel: Spline aus einer Wertetabelle mit ODE 45 integrieren ?
Guten Abend!
Ich kämpfe nun schon einige Wochen mit meiner Matlab Aufgabe und komme keinen Schritt mehr weiter.
Was möchte ich?
Ich habe eine Wertetabelle mit Funktionswerten eines Drehmomentes über die Zeit gegeben. Über diese Momente möchte ich anhand des Drallsatzes der Mechanik die Beschleunigung über die Zeit erhalten.
Danach möchte ich eine Splinefunktion erzeugen, die mir für jeden möglichen Zeitpunkt den zugehörigen Beschleunigunswert zurückgibt.
Integriere ich nun diesen Beschleunigungs-Spline erhalte ich die Winkelgeschwindigkeit.
Integriere ich dann die Winkelgeschwindigkeit bekomm ich als Ziel meinen zurückgelegten Drehwinkel.
Was klappt/ was klappt nicht?
Aus meiner Momentenwertetabelle erhalte ich ganz einfach über den Zusammenhang PHI_dt_dt(i)=M(i)/I meine Winkelbeschleunigunswertetabelle.
Nun habe ich mit der Funktion spline einen kubischen Spline zu dieser Winkelbeschleunigunswertetabelle erhalten. Meine Funktion gibt mir für jeden gewünschetn Zeitpunt den zugehörigen Wert der Winkelbeschleunigung zurück, klappt sehr gut.
Nun möchte ich jedoch diese Winkelbeschleunigung mit dem Solver ode45 über die Zeit integrieren. Leider habe ich aber keine Funktion die ich integriere wie in den Matlab Beispielen sondern nur Funtionswerte.
Die ode45 hat ja eine automatische Schrittweitensteuerung. Das wäre kein Problem da ich ja dank des Spline zu jeden Zeitpunkt einen Beschleunigunswert erhalte, aber wie setze ich das genau in die Tat um?
Ich erhalte bei all meinen Versuchen leider nur Fehlermelungen.
Weiß möglicherweise jemand bescheid, wie ich die Integration mit der Ode45 zum Laufen bringe?
Ich hatte das gesamte Problem schon mit cumtrapz gelöst, jedoch war dies meinem Prof. zu ungenau.
Es macht schon einen großen Unterschied ob ich mit der Trapezmethode integriere oder mit RungeKutta 4ter Ordnung, deshalb verstehe ich auch warum mein Prof. nicht zufrieden war.
Meine gegebene Momentenfunktion habe ich diesem post als PLOT hinzugefügt.
Dieses "erregende" Moment kann positiv oder negativ sein.
Ich habe mal die folgende Funktion erzeugt,da die ODE ja mit einer Funktion versorgt werden muss:
>>
function a_Rad_s2= fun(m_kNm,Iz_kKgm2)
a_Rad_s2=m_kNm/Iz_kKgm2
<<
hierbei ist Iz_kKgm2=2 die konstante Massenträgheit und m_kNm ein Spaltenvektor mit 51 Elementen.
Diese Funtion möchte ich integrieren um auf die Winkelgeschwindigkeit zu kommen.
>>
timerange=[0 5];
initialacceleration=0
%fun(m_kNm,Iz_kKgm2)
[T,Y] = ode45(@fun,timerange,initialacceleration)
<<
Das Ergebnis Y erhale ich jedoch,
Y(1) = 0
alle anderen Y= NaN also keine Zahl
Außerdem wollte ich na nicht mit der Funktion fun in der ODE45 arbeiten, sondern die Werte für die Winkelbeschleunigung aus dem Spline (yi-Wert) holen.
Hier habe ich mit yi gleich einen Vektor gefüllt, was ich ja eigentlich nicht benötige. Ich kann zu jedem geforderten Zeitpunkt ti_s einen Funktionswert haben. Zum Integreieren der DG müsse ich auch noch jeden wert yi durch Iz_kKgm2 dividieren um die Beschleunigung vorliegen zu haben.
for i=1 : length(ti_s)
yi(i)=interpolation(t_s,m_kNm,ti_s(i))
end
<<
Ich hoffe ich hab mich in dieser Schilderung nicht zu unverständlich ausgedrückt. Ansonsten versuch ich gerne das Problem anders zu formulieren.
PS: die im vorhergehneden Beitrag erwähnte Fehlermeldungen hatte ich da ich nicht fun integrieren wollte, sondern versucht habe yi(i) aus dem Spline zu integrieren.
Das Problem unter dem Link habe ich noch nicht ganz verstanden, aber es geht dabei ja darum Parameter zu bekommen. Ich möchte eigentlich nur integrieren und brauche dann in der Integration noch events. Das ist notwendig, da zu meinem Erregermoment noch ein Bremsmoment hinzukommt, das jeweils gegen die momentane Winkelgeschwindigkeitsrichtung wirkt.
Das Problem bei der Variante mit Splines und ode45 ist, dass du von den Splines einen Interpolationsfehler bekommst und damit die Genauigkeit von ode45 knicken kannst. Ja, es ist ein Verfahren 5. (!) Ordnung - aber unter der Annahme von exakten Daten.
Wenn du es allerdings so angehen willst:
In dem anderen Thread ist die Lösung der Differentialgleichung (bzw. bei dir wäre es das Integrieren) ein Teilproblem, das darin gelöst wird, und du brauchst dir eigentlich nur den entsprechenden Teil herauszuziehen.
Dazu müsste man in jedem Fall auch wissen, zu welchen Zeitpunkten die Daten aufgenommen worden sind. Dann sollte es kein Problem sein, das umzusetzen - ob es vernünftig ist, ist die andere Frage.
Wenn es darum geht, die Lösung wirklich fundiert zu berechnen, und du die Lösung in gleichmäßigen Intervallen gegeben hast, würde ich nach einem Löser hoher Ordnung mit konstanter Schrittweite suchen - musst du evtl. selber nachprogrammieren.
In MATLAB gibt es übrigens interp1 für Interpolation...
Hallo!
Da hast du wohl recht, meine Lösung ist nicht wirklich optimal, aber ich bin gern für Anregungen offen.
Du hast ja auf die Funktion interp1 hingewiesen, diese hatte ich auch schon probiert. Dachte nur ob ich jetzt die Spline-Funktion verwende oder interp1-Funktion mit spline'als Methode anstatt default linear macht wenig Unterschied.
Ich habe im anderen Thread nach den Teilproblem (Lösung der DG) gesucht, aber kann diese Lösung leider nicht auf mein System, ohne dabei einen Fehler zu begehen, anwenden.
Ich habe erstens versucht das Bsp aus dem Thread auf meinem PC zum laufen zu bekommen, um es zu verstehen, aber bin daran gescheitert.
Der Interessante Teil für mich ist (habe den hoffentlich richtig erkannt):
>>
function f = diffeq(t, Y1in, x, p)
f = (x(2)*Y2eval(t,p) - (x(3)+2/15*(Y1in-Y3eval(t,p)))*(Y1in-Y3eval(t,p)))/(2*x(1));
end
[tsim, ysim] = ode23s(@(t,y) diffeq(t,y,x,p), p.t, p.Y1(1), opts);
<<
Sein.
Auf mein System übertragen müsste es dann heißen:
>>
function a_Rad_s2= fun( ti_s, yi, mBr_kNm, Iz_kKgm2 )
a_Rad_s2 = (yi - mBr_kNm)/Iz_kKgm2
end
%Die Funktion der Winkelgeschwindigkeit ist abhängig von:
%ti_s =Vektor mit Zeitpunkten, zu denen ich interpoliertes Erregermoment
%yi aus der ursprünglichen Erregermomenwertetabelle m_kNm ermittelt habe
%yi =Wertetbelle der Erregermoment zu Zeitpunkten ti_s
%mBr_kNm =Bremsmoment gegen die Bewegungsrichtung (ab Zeitpunkt t=0)
%aktiv, wirkt also auch schon bei Stillstand und wirkt gegen die
%Beschleunignug.. PS: Ich brauche noch Fallunterscheideung das Bremsmoment nur
%verzögern und nicht beschleunigen kann
%Iz_kKGm2 =Massenträgheitsmoment des Reibschwingers um die Drehachse
%Was natürlich schön wäre ist die jeweiligen, von der ODE45 genommene
%Integrtionszeitpunkte aus der ODE45 zu erhalten und diese in meinen ti_s
%Zeitvektor zu laden. (Outputarguments TE ,IE der ODE??)
<<
Diese Function habe ich als eigens m-file gespeichert.
Im Hauptskript steht:
>>
a_Rad_s2=fun(ti_s, yi, mBr_kNm, Iz_kKgm2)
timerange=[0 5];
initialacceleration=0
[T,Y] = ode45(@fun,timerange,initialacceleration)
<<
Wenn ich dises aber laufen lasse kommen folgende Fehlermeldungen:
??? Input argument "mBr_kNm" is undefined.
Error in ==> fun at 3
a_Rad_s2 = (yi - mBr_kNm)/Iz_kKgm2
Error in ==> odearguments at 110
f0 = feval(ode,t0,y0,args{:}); % ODE15I sets args{1} to
yp0.
Error in ==> ode45 at 173
[neq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0,
odeArgs, odeFcn, ...
Error in ==> file_V6 at 98
[T,Y] = ode45(@fun,timerange,initialacceleration)
Wie kann das sein?
Für >> a_Rad_s2=fun(ti_s, yi, mBr_kNm, Iz_kKgm2) <<brauche ich ja auch mBr_kNm und da regt Matlab sich nicht auf. mBr_kNm als globale Variable zu definieren hilft auch nicht.
Lasse ich die Zeile [T,Y] = ode45(@fun,timerange,initialacceleration) weg erhalte ich die Beschleunigung a_Rad_s2 ohne Fehlermeldungen.
Habe ich diesen Fehler, da ich der ODE gleich die Funtionswerte yi für die ODE eigenen Zeitpunkte übergeben muss?
Ich habe versucht micht besser an der Funktion im empfohlenen Thread zu halten. Danke Harald!
Beim Plot habe ich jedoch gemerkt, dass das Integrationsergebnis nicht ganz stimmen kann.
Das Moment wirkt die ersten ~1,7Sekunden in die negative Richtung > Die Geschwindigkeit ist negativ und wird immer Größer
das passt ja
Ab ~1,7 Sekunden hab ich aber eine positive Beschleunigung, dementsprechend müsste die Geschwindigkeit gegen 0 hin ja abnehmen und nicht weiter in negative Richtung zunehmen.
ein Problem sehe ich auf jeden Fall darin, dass deine fun - Funktion die Variable ti_s gar nicht verwendet. Da das Moment zeitabhängig ist, kann das nicht ja eigtl. nicht hinkommen.
Zudem wirst du interpolieren müssen, damit du das Moment zu dem von ode45 angeforderten Zeitpunkt bekommst - siehe der andere Thread (Funktionen wie Y2eval).
Du hast ja gesagt, du hast mit cumtrapz schon eine Lösung gefunden. Das sollte dir zumindest einen qualitativen Anhaltspunkt geben.
Bitte auch für Code die Code-Tags verwenden. Das macht das ganze zumindest etwas übersichtlicher.
Guten Abend!
Erstens, vielen Dank Harald, für deinen Hinweis zu ti_s und all die anderen wertvollen Informationen.
Ich habe jetzt mein Skript auf das wesentliche gekürzt um die Übersicht zu kriegen, da ich leider immer noch nicht verstehe wie mein Programm richtig funktioniert.
Hier das neue Skrikpt:
Code:
%Reibschwinger
clearall
importfile1('Momentenkurve.txt')
%Zeitvektor ##############################################################
p.t_s = [1 : length(data)]
p.t_s =p.t_s'
%Momentenvektor Erregung #################################################
p.m_kNm=[1 : length(data)]
p.m_kNm=p.m_kNm'
%Massentraegheit des Schwingers um die Z-Achse############################
Iz_kKgm2=2
%Bremsmomente#############################################################
mBrH_kNm=1 %für Versuch bis ODE läuft aus DG weggelassen
%Füllen des Momentenvektors und des Zeitvektors #########################
for i=1 : length(data)
p.t_s(i)=data(i, 1)
m_kNm(i)=data(i, 2) end
%Interpolationsfuntion für das Moment zu Zeitpunkten T welche die Ode
%braucht###########als m_kNmeval.m gespeichert###########################
%function m_T = m_kNmeval(T,p)
%m_T=interp1(p.t_s,p.m_kNm,T,'spline')
%end
dieses @(t,y) diffeq (t,y,x,p),.. steht.
Was ist dieses t? Tsim wären ja die Zeitpunkte an denen die ODE integrieren möchte, was ist dann die Aufgabe des t? Dieses p.t brauch ich nur für die Interpolation und das hat mit t ja auch nix zu tun.
Bei mir müsste ich anstelle von y ein ysim schreiben, aber was schreibe ich für das t hin?
Laut Matlab Hilfe schreibe ich als erstes nach ode45( ja den functionhandle @diffeq, aber auch bei dieser Variante kommen nur Fehlermeldungen.
Für das von mir oben gepostete Skript erhalte ich folgende Fehlermeldungen:
??? Error: File: diffeq.m Line: 2 Column: 47
Unbalanced or unexpected parenthesis or bracket.
Error in ==> @(T,Y)diffeq(T,Y,p.t_s,p.m_kNm,Iz_kKgm2)
Error in ==> odearguments at 110
f0 = feval(ode,t0,y0,args{:}); % ODE15I sets args{1} to
yp0.
Error in ==> ode45 at 173
[neq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0,
odeArgs, odeFcn, ...
Error in ==> Forum at 24
[T,Y]=ode45(@(T,Y)
diffeq(T,Y,p.t_s,p.m_kNm,Iz_kKgm2),timerange,p.m_kNm(1),opts)
>>
Weiß jemand wo der Fehler liegt?
Mein Ergebnis mit CUMTRAPZ hab ich als Plot hinzugefügt.
Darin ist die schwarze Kurve meine Momentenfuktion, die grüne Kurve das effektive Drehmoment unter Berücksichtigung der Bremse.
Rot ist die Winkelbeschleunigung. (nicht ganz ok).
Die dunkelblaue Kurve ist die Winkelgeschwindigkeit.
Hellblau der Betrag des zurückgelegten Drehwinkels.
Da möchte ich auch mit der ode45 hin.
Schönen Gruß,
Gerhard
Momentenkurve.txt
Beschreibung:
Wenn jemand das Skript probieren möchte, hier die Erregermomentenwertetabelle
ich fürchte, da gibt es ein grundlegendes Verständnisproblem. ode45 wertet wie die verwandten Methoden die DGL-Funktion nicht an den Stellen aus, die man vorgibt, sondern an durch die Schrittweite bestimmten Zeitpunkten. Man muss also eine Funktion schreiben, die t und y als Eingabeargumente annimmt und y' zurückgibt.
In deinem Skript ist der Fehler schlicht eine fehlende Klammer beim ode45-Aufruf. Das wird sicher auch im Editor rot unterringelt.
Ich habe hier mal ein Minimalskript zur Integration von Beschleunigung auf Geschwindigkeit.
Du hast recht, es war ein Verständnisproblem.
Vielen Dank für dein Skript.
Jetzt funktioniert´s.
Schönen Gruß,
Gerhard
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.
RSS