% 28: Kinetik des Massenpunktes
% Beispiel 2: Berechnung eines Gleitsteines S.501
% Author:  Eugen Lass
% Revision: Prof. Dr. Thomas Frischgesell
% Datum:   03.04.2003

% Lösung der Differentialgleichung

clear all;

options = odeset('MaxStep',0.1); % Ohne Schrittweitensteuerung erhalten wir ein falsches Ergebnis

% setze Anfangswerte, Zeitspanne
x0 = [-4.493; 0]; % Anfangswerte [xi ; xipunkt]
tspan = [0 20];  % Zeitspanne

% Integration der DGL
[t x] = ode45('S501_xpunkt', tspan, x0,options);

% grafische Ausgabe von xi und xipunkt
xi = [1 0]*x';
xipunkt = [0 1]*x';

plot(t,xi,t,xipunkt)

% legt die Achsen fest
grid on
xlabel('Zeit [s]');title('Vertikale Führung')
legend('xi [rad]','xipunkt [rad/s]',4)

% Gleitstein_b
% Ableitungsberechnung des Gleitsteines

function xpunkt = xpunkt_Gleitstein(t,x)

% Gegeben
A = 1;             % A = c * a / (m * g)
B = 4;             % B = b / a
C = 9.81;          % C = g / a    in 1/s²

% Berechnung von xpunkt
xpunkt = [x(2); C * (1 - A * (sqrt(1 + x(1)^2) - B) * x(1) / sqrt(1 + x(1)^2))];

Nach dem Starten von Matlab wird File|Open gewählt und im Fenster "Open" (unter "Dateityp" sollte "All MATLAB Files" eingestellt sein) wird in das Verzeichnis gewechselt, in dem sich die Dateien S501.m und S501_xpunkt.m befindem. Doppelklick auf S501.m öffnet die Datei in einem neuen Fenster.

In diesem Fenster wird Debug|Run gewählt. Möglicherweise öffnet sich ein Fenster "MATLAB Editor", in dem man das Angebot "Change MATLAB current directory" annehmen sollte. Danach wird die Rechnung ausgeführt.

Die mit plot ... eingeleitete Zeile am Ende des Scripts führt zur grafischen Ausgabe des Bewegungsgesetzes (nebenstehende Abbildung).

Bei einem Fall aus Höhe x = −4,493 gelingt es der Masse gerade noch, sich an dem Lager "vorbeizudrängeln". Die Geschwindigkeit wird vorübergehend sehr klein, bis der Punkt passiert ist, an dem die Feder am stärksten zusammengedrückt ist.

Bei einem Fall aus der nur unwesentlich geringeren Höhe x = −4,492 (in der Datei S501.m wird nur eine Zeile geändert: x0 = [-4.492;0];) gelingt es der Masse nicht mehr, sich an dem Lager "vorbeizudrängeln". Sie führt eine Schwingung oberhalb des Punktes der maximalen Zusammendrückung der Feder (ausschließlich im Bereich negativer x-Werte) aus (nebenstehende Abbildung).

In der einfachsten Aufrufvariante der Funktion ode45 brauchen diese Optionen nicht angegeben zu werden (sie müssen dann auch nicht mit odeset gesetzt werden). Wenn in der Datei S501.m die entsprechende Zeile in der einfacheren Form

[t x] = ode45('S501_xpunkt', tspan, x0);

verwendet wird, arbeitet Matlab mit den Default-Einstellungen. Für die Anfangswerte

x0 = [-4.493; 0]; % Anfangswerte [xi ; xipunkt]

führt dies zu dem nebenstehend zu sehenden falschen Ergebnis (die numerischen Fehler führen dazu, dass sich die Masse nicht an der Feder "vorbeidrängeln" kann).

Dass die Rechnung fehlerhaft ist, kann man natürlich kaum erkennen, wenn man das richtige Ergbnis nicht kennt. Man sollte aber unbedingt gewarnt sein.

Manchmal sind fehlerhafte Ergebnisse allerdings sofort zu erkennen. Wenn man die gleiche Rechnung (mit den Matlab-Default-Werten) mit den Anfangswerten

x0 = [-4.5; 0]; % Anfangswerte [xi ; xipunkt]

ausführt, ergeben sich nebenstehend zu sehende Kurven: Die Masse schafft es (das ist richtig), sich bei ihrer Abwärtsbewegung an der Feder "vorbeizudrängeln", auf dem Rückweg allerdings nicht mehr. Weil das mathematische Modell aber keine Reibung berücksichtigt (und somit keine Energieverluste entstehen), müsste sie stets die Anfangshöhe wieder erreichen.