|
Die Definition von “Events” wird in der (hellblauen) Options-Zeile angekündigt:
Das Schlüsselwort ‘Events’ verweist auf ein Handle einer Funktion Ereignis (Name kann beliebig gewählt werden), die die “Events” definiert. Hier werden die drei Ereignisse
"Geschwindigkeit Null", "Punkt x = a erreicht, weiter in Richtung x > a" und "Punkt x = a erreicht, weiter in Richtung x < a" definiert:
function [value,isterminal,direction] = Ereignis (t,xvt)
global m c my a g r myq
x = xvt(1) ;
v = xvt(2) ;
value(1) = v ; % Ereignis 1: "Geschwindigkeit Null"
isterminal(1) = 1 ; % Abbruch, ...
direction(1) = 0 ; % ... in jedem Fall
value(2) = x-a ; % Ereignis 2: "Punkt x=a erreicht"
isterminal(2) = 1 ; % Abbruch, ...
direction(2) = 1 ; % ... wenn Bewegung nach rechts
value(3) = x-a ; % Ereignis 3: "Punkt x=a erreicht"
isterminal(3) = 1 ; % Abbruch ...
direction(3) = -1 ; % ... wenn Bewegung nach links
Die Function bekommt als Input die gleichen Werte (Zeit und Vektor der Bewegungsgrößen) wie die Funktion für die Definition der Differenzialgleichungen (hier: RechteSeite). Sie
muss für jedes Ereignis 3 Werte abliefern: Eine “Nullfunktion” (hier auf Vektor value) definiert durch ihre Nullstelle das Ereignis, außerdem: Was soll beim Eintreten dieses Ereignisses geschehen (hier auf isterminal:
1 bedeutet Abbruch, 0 würde Weiterrechnung bedeuten), und soll dies in jedem Fall (die 0 auf direction) geschehen oder nur beim Nulldurchgang in positive Werte (1) oder nur beim Nulldurchgang in negative
Werte (-1) geschehen?
Alle drei Ereignisse führen also in diesem Fall zum Abbruch der Rechnung. Die ode-Funktion liefert bei gesetzter Event-Option drei zusätzliche Ergebnisse:
[t xv te xve ereig] = ode45 (@RechteSeite , tspan , x0 , options) ;
Auf te werden die Zeiten abgeliefert, zu denen die Ereignisse eingetreten sind, auf xve findet man die Zustandsgrößen der Bewegung zu diesen Zeitpunkten (hier: Wegkoordinaten und Geschwindigkeiten), der Vektor ereig enthält die Ereignisnummern (hier 1, 2 oder 3). Nach jedem Abbruch werden diese Informationen ausgewertet, um die Parameter (hier myq und r) zu aktualisieren:
lent = length (t) ;
if (t(lent) >= tend)
cont = 0 ; % Ende erreicht, Abbruch der while-Schleife
else % ... war ein "Ereignis" der Grund fuer den abbruch
telen = length (te) ;
tspan = [te(telen) tend] ; % ... weiter mit "Restintervall ...
x0 = [xve(telen,1) xve(telen,2)] ; % ... und neuen Anfangsbedingungen
if ereig(telen) == 1 % Ereignis "Geschwindigkeit Null"
if abs(c*xve(telen,1)) < my0*m*g ...
& xve(telen,1) > a % Federkraft kleiner Haftkraft, ...
cont = 0 ; % ... also Abbruch
tout = [tout ; tend] ;
xvout = [xvout ; [xv(lent,1) 0]] ;
elseif xve(telen,1) < 0 r=1 ; % Es geht nach rechts weiter ...
else r=-1 ; % ... bzw. nach links
end
elseif ereig(telen) == 2 myq=my ; % Ereignis "Ab sofort Reibung"
elseif ereig(telen) == 3 myq=0 ; % Ereignis "Ab sofort reibungsfrei"
end
end
Weil die Rechnung mehrfach abgebrochen wird, ist die ode45-Funktion in eine while-Schleife eingebettet. Die Ergebnisse, die nach jedem ode45-Aufruf auf dem Vektor t und der Matrix xv abgeliefert werden, werden deshalb auf tout und xvout kumuliert (grüne Zeilen):
tout = [] ; % Auf tout und xvout werden die Ergebnisse gesammelt
xvout = [] ;
...
while cont == 1
[t xv te xve ereig
] = ode45 (@FederMasseReibung_Dgl , tspan , x0 , options) ;
tout = [tout ; t] ; % Die gerade abgelieferten Ergebnisse werden an ...
xvout = [xvout ; xv] ; % ... die bereits gespeicherten Ergebnisse angehaengt
...
end
Außerdem wird nach jedem Abbruch die Rechnung für das Restintervall vorbereitet, indem tspan vom aktuellen Zeitpunkt bis zum Ende des Intervalls gesetzt und als Anfangsbedingungen die Endwerte des gerade abgearbeiteten Intervalls gesetzt werden:
telen = length (te) ;
tspan = [te(telen) tend] ; % ... weiter mit "Restintervall ...
x0 = [xve(telen,1) xve(telen,2)] ; % ... und neuen Anfangsbedingungen
Bei einem Aufruf der Function Aufg28_5Ev ohne Spezifizierung von Parametern werden die Werte der Aufgabenstellung genommen, und man erhält die Lösung, die auch mit der “nicht ganz sauberen Strategie” erzielt wurde. Bei einem Aufruf der Function
Aufg28_5Ev aus dem Command Window mit den Parametern μ = 0,4 (Gleitreibungskoeffizient) und μ 0 = 0,8 (Haftungskoeffizient), für die die einfache Realisierung der Parameteränderung versagte, erhält man nun das erwartete Ergebnis mit einer Endlage der Masse im Bereich x > a.
Nachfolgend sind Command Window und Graphik-Fenster dieser Berechnung zu sehen:
|
Bewegungsdifferenzialgleichung:
