Beispiel: Funktion mit komplizierter Beschreibung

Ein dünner Stab der Länge l mit konstantem Querschnitt ist an einem Ende reibungsfrei gelagert. Er wird aus der vertikalen Lage um den Winkel φ0 ausgelenkt und ohne Anfangsgeschwindigkeit freigegeben (Luftwiderstand kann vernachlässigt werden).

Gegeben:  Erdbeschleunigung g = 9,81 m/s2 ,  φ0 = 3π/4.

Es soll die Abhängigheit der für eine volle Schwingung benötigten Zeit T (Schwingungsdauer, "einmal hin und zurück") von der Pendellänge l als Funktion T(l) für den Bereich 1 cm ≤ l ≤ 1 m ermittelt werden.

Das Problem besteht darin, dass selbst das Bewegungsgesetz φ(t) nicht mit einem analytischen Ausdruck formuliert werden kann. Im Kapitel "Kinetik starrer Körper" des Lehrbuchs "Dankert/Dankert: Technische Mechanik" wird gezeigt, dass für das Bewegungsgesetz das folgende Anfangswertproblem gilt:

Anfangswertproblem 2. Ordnung

Dieses nichtlineare Anfangswertproblem kann nur numerisch integriert werden. Man findet die Lösung dieses Problems als Beispiel 1 für das Programm "Anfangswertproblem" (unter "TM-interaktiv"), außerdem ist hier die Lösung mit Matlab beschrieben. Mit den dort demonstrierten Lösungen ist die Berechnung nur jeweils eines Punktes der Funktion T(l) möglich.

Lösung mit Matlab

Das oben angegebene Anfangswertproblem wird mit dem Matlab-Standardsolver ode45 gelöst. Es genügt, die Zeit für eine halbe Schwingung (φ = 3π/4 ... −3π/4) zu ermitteln, die Schwingungsdauer T ist das Doppelte dieses Wertes. Dies wird mit der Matlab-Event-Strategie realisiert, als "Ereignis" wird "Winkelgeschwindigkeit gleich Null" (am Umkehrpunkt der Bewegung) definiert. Dann ist ein Wert Ti(li) bekannt. Um diesen Algorithmus wird eine Schleife gelegt, die die Berechnung für alle gewünschten li-Werte durchführt.

Das nachfolgend gelistete Matlab-Script realisiert diese Rechnung:

function pendel ()

global pl g
          g = 9.81 ;

plbereich = [0.01 1] ;
nsteps    = 30 ;
plstep    = (plbereich(2) - plbereich(1))/nsteps ;
plen = zeros (nsteps+1 , 1) ;
tend = zeros (nsteps+1 , 1) ;


phi0  = [3*pi/4 ; 0] ;      % Anfangswerte
tspan = [0 1000]     ;      % Zeitintervall (groß genug)

% Optionen für Integration des Anfangswertproblems:
options = odeset ('MaxStep' , 0.01 , 'Events' , @Ereignis) ;

pl = plbereich(1) ;
for i=1:nsteps+1
   [t xv te xve ereig]  = ode45 (@BewDglPendel , tspan , phi0 , options) ;
   plen(i) = pl          ;
   tend(i) = t(end) * 2   ;     % Zeit für eine volle Schwingung
   pl      = pl + plstep ;
end


plot (plen , tend) , grid on , title ('pl-T-Diagramm:')

[plen tend]        %  ... schreibt Wertetabelle in das Command Window

% ===================================================================
function [value,isterminal,direction] = Ereignis (t,xvt)
global pl g

omega = xvt(2) ;
value(1)      = omega ; % Ereignis 1: "Winkelgeschwindigkeit Null"
isterminal(1) = 1 ;     % Abbruch, ...
direction(1)  = 0 ;     % ... in jedem Fall

% ===================================================================
% Funktion, die die "rechte Seite" des Dgl.-Systems definiert:

function xvp = BewDglPendel (t , xv)
global pl g                    % ... sind in oben definiert.
phi    = xv(1) ;               % Input-Werte kommen in ...
omega  = xv(2) ;               % ... einem Vektor an.
phip   = omega ;               % Die beiden Differenzialgleichungen
omegap = -1.5*g/pl*sin(phi) ;  % werden ausgewertet, ...
xvp    = [phip ; omegap] ;     % ... die Ergebnisse werden im Vektor xvp abgeliefert.

Schwingungsdauer in Abhängigkeit von der Pendellänge Nebenstehend sieht man die Wertetabelle Ti(li), nachfolgend die grafische Darstellung der Funktion T(l):

Schwingungsdauer in Abhängigkeit von der Pendellänge

Das oben gelistete Matlab-Script steht zum Download zur Verfügung.