_____________________________________________________________________________________
Inhalt....: Bewegung eines Kugelpendels im Gravitationsfeld
Kategorie.: Unterrichtsmaterial
Mathematik: Geometrie R^3, Physik
MuPAD.....: 3.1.1
Datum.....: 2006-02-23
Autoren...: Günter vom Stein <gvomstein@gmx.net>
Funktionen: linalg::crossProduct, linalg::ScalarProduct, norm
_____________________________________________________________________________________
Bewegung einer Kugel im Gravitationsfeld (z.B. der Erde). Die Kugel wird mit
einer (masselosen) Stange an einem ruhenden Punkt aufgehängt.
(Luft-)Reibung kann zugeschaltet werden.
Zur Idee:
Ein Massepunkt bewegt sich auf einer Kugelschale. Er wird zu Beginn auf die Kugeloberfläche
gesetzt und bekommt einen Geschwindigkeitsvektor verliehen, der in der Tangentialebene zum
Startpunkt liegt.
In der folgenden Iteration wird bei jedem Schritt der Geschwindigkeitsvektor unter dem Einfluss
der Schwerkraft neu orientiert. Anschließend legt der Massepunkt mit dieser Geschwindigkeit
ein Weg-Intervall auf demjenigen Großkreis zurück, den die Ursprungsebene, gebildet aus
momentanem Geschwindigkeitsvektor und Ortsvektor (= Fahrstrahl), aus der Kugel schneidet.
In der Simulation "kennt" der Massepunkt nicht mehr als in der Realität auch, nämlich die jeweilige
Tangentialkomponente der Gravitationsbeschleunigung und die Zwangsbedingung, an der Stange
festgemacht zu sein. Letztere Bedingung kommt in der Simulation zum Ausdruck, wenn der Weg
entlang des Großkreises als Drehung des Ortsvektors in der entsprechenden Ursprungsebene erfolgt.
Als mathematische Werkzeuge werden nur das Skalarprodukt und das Kreuzprodukt eingesetzt.
Insbesondere kann die Drehung elegant mit dem Kreuzprodukt realisiert werden.
___________________________________________________________________________
Als Anschauungshilfen stehen drei Graphiken zur Verfügung, auf die im Text der Iterations-Prozedur
Bezug genommen wird.
Graphik 1 und 2 sind aktive Graphiken (Doppelklick oder rechte-Maus-Taste-Graphik-Objekt-Open).
Graphik 1 stellt mit der blauen Gerade die Richtung der Rückstellkraft bzw. Rückstellbeschleunigung
dar. Die Beschleunigung zeigt auf der unteren Halbkugel vom Ort des Massenpunkts zur z-Achse hin,
auf der oberen Halbkugel von der z-Achse weg.
Graphik 2 zeigt die Ursprungsebene, in der der Ortsvektor des Massepunkts (= Fahrstrahl) gedreht wird,
sowie den momentanen Geschwindigkeitsvektor.
Graphik 3 beschreibt die Drehung eines Ortsvektors um eine vorgegebene Normale.
Wegen 
sind die drei Vektoren
und 
alle gleich lang. Dann aber gilt
(aus Kroll u.a.,Analytische Geometrie/Lineare Algebra, Dümmler-Verlag, S.70)
-------------------------------------------------------------------------------------------------------
Beginn des Programms
Bei 3000 Iterationsschritten werden für den Prozeduraufruf und die beiden Graphiken auf einem
Notebook 1,4 GHz jeweils ca. 10 Sekunden benötigt.
reset():
DIGITS:=10:
export(linalg):
alias(KreuzP = crossProduct):
alias(scP = scalarProduct):
Festlegung der Startwerte
x0 und y0 können unter Einhaltung von 
frei gewählt werden, z0 wird vom Programm
berechnet.
Der Startvektor der Geschwindigkeit muss in derjenigen Tangentialebene an die Kugel liegen,
die vom Startpunkt definiert wird.
Wenn vx0 und vy0 entgegengesetzt gleiche Werte haben, muss vz0 von Hand auf Null gesetzt werden,
denn solve(scP(s0,v0),vz0) führt dann nicht zu einem Ergebnis. Andernfalls wird vz0
ebenfalls vom Programm berechnet. solve(scP(s0,v0),vz0) macht nur einen Vorschlag
unter vielen möglichen.
sx0:=5:sy0:=5:vx0:=-15: vy0:=25: // Anfangswerte
//sx0:=1:sy0:=1:vx0:=0: vy0:=5: // Anfangswerte
g:=9.81: // Gravitationskonstante
r:=10:
//sz0:=float(sqrt(r^2-sx0^2-sy0^2));
sz0:=-float(sqrt(r^2-sx0^2-sy0^2));
s0:=matrix([sx0,sy0,sz0]):
v0:=matrix([vx0,vy0,vz0]):
L:=solve(scP(s0,v0),vz0):
vz0:=op(op(L,1),2);
//vz0:=0:
scP(s0,v0); // Kontrolle: v in der Tagentialebene ?
Error: expecting an arithmetical expression [normal]
Die Prozedur zur Iteration:
Es ist weiter unten möglich vor Aufruf der Prozedur einen Reibungsfaktor festzulegen, der zwar
hier nur willkürlich gewählt werden kann, aber einen Reibungsterm 
definiert, wie er bei
Luftreibung auftritt.
Bahn:=proc(s0,v0,n)
local a0,n0,KP,v,alpha,nz;
begin
dt:=delta_t/1000; // dt Zeitintervall in Sekunden
s[0]:=s0;
v:=v0;
nz:=matrix([0,0,1]);
for i from 0 to n do
// Rückstellbeschleunigung: Hierzu studiere man Graphik 1 (s.o)
KP:=float(KreuzP(s[i],KreuzP(s[i],nz)));
a0:=KP/float(norm(KP,2));
//a0:=-KP/float(sqrt(scP(KP,KP)));
// neue Geschwindigkeit:
v:=v+g*float(sin(arccos(abs(s[i][3]/r))))*dt*a0;
v:=v*(1-C*scP(v,v)*dt); // zusätzlich (Luft)Reibung (s.u.)
// Normalenvektor für Drehung bestimmen: s. Graphik 2
KP:=float(KreuzP(s[i],v));
//n0:=KP/float(sqrt(scP(KP,KP)));
n0:=KP/float(norm(KP,2));
// Drehwinkel im Bogenmaß bestimmen:
//alpha:=float(norm(v,2)*dt/r); // "norm" verlangsamt hier ungemein!?
alpha:=float(sqrt(scP(v,v))*dt/r);
// Drehung von s und v: s. Graphik 3
s[i+1]:=cos(alpha)*s[i]+sin(alpha)*KreuzP(n0,s[i]);
v:=cos(alpha)*v+sin(alpha)*KreuzP(n0,v);
end_for
end_proc:
C:=0:C:=0.001: // (Luft)Reibungsparameter
delta_t:=10: // Dauer eines Zeitintervalls in ms
n:=3000: // Anzahl der Iterationsschritte
Bahn(s0,v0,n):
k:=5:k_max:=(n div k): // Nur jeder k-te Punkt wird gezeichnet
Punkte:= plot::Point3d([s[k*i][1],s[k*i][2],s[k*i][3]],
PointSize=2,
PointStyle=FilledCircles,
Color=RGB::Red,
VisibleFromTo = i..i+1
)
$ i=0..k_max:
z_Achse:=plot::Line3d([0,0,0],[0,0,-r],
LineColor = RGB::Black
):
Seil:=plot::Line3d([s[k*i][1],s[k*i][2],s[k*i][3]],[0,0,0],
LineColor = RGB::Red,
VisibleFromTo = i..i+1
)
$ i=0..k_max:
Bahn:=plot::Line3d([s[k*i][1],s[k*i][2],s[k*i][3]],[s[k*(i+1)][1],s[k*(i+1)][2],s[k*(i+1)][3]],
LineColor = RGB::Red,
VisibleFromTo = i+1..i+1,
VisibleAfterEnd
)
$ i=1..k_max-1:
Suedpol:=plot::Point3d([0,0,-r] ,PointSize=1,
PointStyle=FilledCircles,Color=RGB::Black):
Nordpol:=plot::Point3d([0,0,r] ,PointSize=1,
PointStyle=FilledCircles,Color=RGB::Black):
b:=1.4:
Sc3:=plot::Scene3d(Punkte,Suedpol,Seil,Bahn,z_Achse,
ViewingBox = [-b*r..b*r, -b*r..b*r,-r..r]):
plot(Sc3, GridVisible = TRUE):
Zur Veranschaulichung mit Kugeloberfläche:
Kugel:=plot::Sphere(r,[0,0,0],FillColor=RGB::LightYellow.[0.85]):
Sc3:=plot::Scene3d(Punkte,z_Achse,Seil,Bahn,Suedpol,Nordpol,Kugel,
ViewingBox = [-r..r, -r..r,-r..r]):
plot(Sc3, GridVisible = TRUE):
Die Güte der Simulation kann weder an theoretischen Gleichungen und Ergebnissen noch am
Experiment gemessen werden, wohl aber an sich selbst.
Beispiel für den Vergleich zweier Durchgänge:
Anzahl der Iterationsschritte Dauer eines Zeit-Intervalls Endpunkt der Bahn
3000 10 ms [ 2.9726 | 1.1808 | -9.4747 ]
6000 5 ms [ 2.9500 | 1.2285 | -9.4757 ]
_______________________________________________________________________________
Anmerkungen:
1. Weitere Anregungen zum Einsatz von MuPAD in der Lehre finden Sie auf unserem WebPortal
MuPAD in Schule und Studium unter: http://schule.mupad.de bzw. http://studium.mupad.de.
_______________________________________________________________________________