I/T5 - WS 04/05 CR – Übung 5 Termin: 03.-07.12.04 Modellierung eines Trampolins, mit Hilfe eines dünn besetzten Gleichungssystems Der Bewegungsablauf einer Kugel, die auf eine h =

(1)

I/T5 - WS 04/05 CR – Übung 5 Termin: 03.-07.12.04

Modellierung eines Trampolins, mit Hilfe eines dünn besetzten Gleichungssystems Der Bewegungsablauf einer Kugel, die auf eine h = 4 quadratische Trampolin-Plane fällt, soll gra-

phisch dargestellt werden. Die Kugel wird von der Plane aufgefangen und wieder nach oben geschleudert. Das Rechenproblem dabei ist die graphische Darstellung der Plane, die von der Kugel ausgebeult, aber nicht durchdrungen werden soll.

Plane h = 0

Wir verwenden die Einheitskugel (Kugel um den Nullpunkt mit Radius 1), die um die Höhe h in Richtung der z-Achse verschoben wird.

Kugelschale: x²+y²+ −(z h)² =1, Kugelinneres: x²+y²+ −(z h)² <1.

h = -2

Die Plane befindet sich in der Höhe h = 0 und belegt das Rechteck

( )

^{x y}^, ^{∈ −}

[

^{2, 2}

] [

^{× −}^2,^{2 .}

]

Die Kugel soll aus der Höhe h = 4 herab fallen und bei Erreichen der Höheh= −2 zum Stillstand kommen. Dabei wird die Plane stark ausgebeult (s. Skizze).

Kugel und Plane werden für die Graphik als räumliche Punktraster mit Matrizen der (x,y,z)- Koordinaten dargestellt:

• Matrizen zur Beschreibung der Einheitskugel: U,V,W; zur Beschreibung der Plane: X,Y,Z.

• Kugel in Höhe h: U,V,W+h, wobei h zu jeder Komponente von W addiert wird.

Die Form der Plane soll folgenden Bedingungen unterliegen:

• Der Planenrand bleibt fest. Die inneren Planenpunkte verschieben sich evtl. in z-Richtung, wenn die Kugel aufliegt. Die Koordinate z_mn jedes inneren Planenpunkts ergibt sich als Mit- telwert der z-Koordinaten seiner 4 Nachbarpunkte: ¹

(

^1, ^1, ^, ¹ ^, ¹

)

mn 4 m n m n m n m n

z = z ₊ +z ₋ +z ₊ +z ₋ .

• Falls das so berechnete z_mn im Kugelinneren liegen würde, wird sich die Plane an dieser Stelle anders verhalten und der Kugelschale anschmiegen. D.h. für z_mn wird statt des berechneten Mittelwerts der z-Wert eingesetzt, der sich aus der Gleichung der Kugelschale (s.o.) ergeben würde: z_mn = −h 1−x_mn² −y_mn² .

Sind die Planenpunkte mit m=1,...,M,n=1,...N durchnummeriert, entsteht so formal ein großes lineares Gleichungssystem für alle Unbekannten z_mn,m=2,...,N−1,n=2,...,M −1. Komplikation dabei ist die Lage einer variablen Anzahl von Planen-Punkten auf der Kugelschale. Wir wollen dieses Gleichungssystem iterativ mit dem Jacobi-Verfahren lösen und bei jedem Iterationsschritt überprüfen, welche Rasterpunkte z_mn auf der Kugelschale liegen.

Entwickeln Sie ein MatLab-Programm, das folgendes leistet (mit nur einer for-Schleife):

(2)

ERSTER SCHRITT: Startposition mit straffer Plane

Wichtig für die Aufgabe ist ein fixiertes Plot-Fenster in passender Größe. Verwenden Sie fi- gure, um ein Plot-Fenster zu erzeugen, das in etwa die linke Hälfte des Bildschirms einnimmt und komplett sichtbar ist. Das Matlab-Fenster stellen Sie auf die rechte Hälfte des Bild- schirms.

Erzeugen Sie mit Hilfe des Befehls sphere ein 3D-Punktraster für die Einheitskugel, bestehend aus 3 Matrizen U,V,W.

Erzeugen Sie mit Hilfe des Befehls meshgrid ein 2D-Punktraster für die Plane, bestehend aus 2 Matrizen X,Y für x und y-Werte der Plane. x und y sollen von -2 bis 2 laufen, in Schritten von 0,2. Erzeugen Sie für die z-Werte der Plane eine Null-Matrix Z gleicher Größe. Speichern Sie mit size das Format der Matrix Z in 2 Zahlen M,N. – Unsere Unbekannten sind damit die Zahlen Z(2:M-1,2:N-1).

Zeichen Sie mit surf die Kugel in Höhe h = 4 (!) sowie die Plane (hold on beachten). Wichtige Zusätze mit axis-Befehl: identische xy-Achsenskalierung, Plotbereich fixiert auf 2− ≤ ≤x 2,

, . pause-Befehl zum Betrachten des Bildes, dann hold off, damit das Bild anschließend für ein anderes h erneuert werden kann.

2 y 2

− ≤ ≤ − ≤ ≤3 z 5

ZWEITER SCHRITT: Zwischenposition mit gebogener Plane

Nächstes Bild: Kugel in Position h = 0, Plane hängt durch. Die Form der Plane muss zuerst mit den oben beschriebenen Formeln berechnet werden. In einer Iterationsschleife (while …) wird Z so lange verändert, bis die Planenform stimmt. D.h. es wird bei jedem Durchlauf

1

(

(2 : 1, 2 : 1) (.,.) (.,.) (.,.) (.,.)

Z M − N− = 4 Z +Z +Z +Z

)

gesetzt, mit entsprechend verschobe- nen Teilmatrizen Z(.,.). Die Indizes in den Klammern sollen Sie passend einsetzen.

In jedem Iterationsschritt wird mit Hilfe des Funktionsaufrufs K = find (..) festgestellt, welche Punkte aus Z innerhalb der Kugel liegen würden. Mit den Rasterpunkt-Koordinaten (x,y,z) bedeutet das (z−h)² < −1 x²−y² =b. In K sind nun alle Indizes dieser Punkte gespeichert, so dass Z(K) neu zu belegen ist durch Punkte auf der Kugelschale (s.o.): z= −h b.

In jedem Schleifendurchlauf wird die maximale Punktverschiebung d ermittelt. Dazu wird zuerst Z in einen Zwischenspeicher A kopiert, dann mit den obigen Formeln neu berechnet.

Die maximale Punktverschiebung ist d =max Z - A . Die Iteration endet, wenn d<0, 001 (Plot-Genauigkeit). Geben Sie die Anzahl der dazu benötigten Iterationsschritte aus.

Wie oben geschehen, werden Kugel und Plane mit surf, hold on und hold off in einem neuen Bild dargestellt, das Sie drehen und aus mehreren Richtungen betrachten sollen um festzustel- len, ob die Plane richtig geformt ist.

DRITTER SCHRITT: Animation

Die Berechnung der Planenform ist jetzt getestet und für jede beliebige Höhe h möglich. Sie wird nun in eine for-Schleife eingebettet, in der h mit Schrittweite -0,1 von 4 bis -2 läuft, an- schließend die Richtung wechselt und mit Schrittweite 0,1 bis 4 ansteigt. In jedem Schleifen- durchlauf wird ein Bild geplottet wie oben beschrieben und mit Film(k) = getframe in einem Vektor Film von Graphik-Strukturen gespeichert.

Die Funktion movie erlaubt, den Film beliebig oft hintereinander abzuspielen. Lassen Sie die Kugel 10 mal auf die Plane springen. – Viel Spaß!

Frage: Warum wird bei der Iteration nicht die relative Abweichung zwischen Z und A kontrol- liert?