Numerische Simulation

[0x00000001]

Hallo

ich habe Massepunkte gegeben, und jeder dieser Punkte hat eine Masse, Position, Geschwindkigkeit und anliegende Kraft.
Die Integration will ich nun mit dem Beeman Algorithmus machen, der wie folgt aussieht

$X_{t+h} \approx X\_t + h \cdot V\_t + h^2 \left( \frac{2}{3} \frac{F\_t}{m} + \frac{1}{6} \frac{F\_{t-h}}{m} \right)$

$V_{t+h} \approx V\_t + h \left( \frac{5}{12} \frac{F\_{t+h}}{m} + \frac{2}{3} \frac{F_{t}}{m} + \frac{1}{12} \frac{F_{t-h}}{m} \right)$

X ist Position, V Geschwindigkeit, F Kraft und h das Zeitintervall (1/60)

Mein Problem ist nun wie ich das konkret ausrechnen kann. Vor allem bei den Kräften im vorherigen und nächsten Zeitschritt bin ich mir unsicher.
Die Kraft aus dem letzten Zeitschritt habe ich mir einfach mal zusätzlich gespeichert.
Da F/m = a, hab ich die Beschleunigung im zukünftigen Zeitschritt über die Geschwindigkeit berechnet:

VEC3 V_next = (PosNew - PosOld) / h;
VEC3 A_next = V_next / h;

Weiß aber nicht ob das so OK ist. In meinem Skript steht nur dass man die Position aus dem nächsten Zeitschritt verwenden soll, und das habe ich ja gemacht (mit PosNew)

Irgendwas muss auf jeden Fall falsch sein, denn die Objekte fallen zunächst ganz normal nach unten (wegen Gravitation), werden dann aber plötzlich nach unten gezogen und zerschellen in der Unendlichkeit.

Kann jemand sagen was ich falsch gemacht hab?

[Werner Salomon]

Ich meine, dass Deine Gleichungen schon falsch sind. Wenn Du Vt=0 und Xt=0 und die Kraft als konstant ansetzt, so sollten die bekannten Bewegungsgleichungen herauskommen. Statt dessen erhält man
$X_{t+h} = \frac{5}{6} \frac{F}{m} h^2$
korrekt wäre
$X_{t+h} = \frac{1}{2} \frac{F}{m} h^2$

und bei der Geschwindigkeit
$V_{t+h} = \frac{7}{6} \frac{F}{m} h$
korrekt wäre
$V_{t+h} = \frac{F}{m} h$

das heißt in beiden Fällen zu viel, was vielleicht das von Dir beschriebene Verhalten erklärt.
Ich habe hier (Suche in der Datei nach 'Beeman') was gefunden, die Gleichungen sehen da ganz anders aus.

Gruß
Werner

[0x00000001]

Danke für Deine Antwort.
Ich habe jetzt mal die Variante aus dem PDF probiert, aber das Ergebnis ist ungefähr das gleiche. Das entspricht übrigens auch dem was man bei Wikipedia findet: [url]http://en.wikipedia.org/wiki/Beeman's_algorithm[/url]

Man muss ja die Beschleunigung bzw. die Kraft im nächsten Zeitschritt berechnen, und ich habe die Vermutung dass ich da was falsch gemacht habe. Wenns um Physik geht komm ich leider schnell an meine Grenzen.
Hier mal der komplette Code. Vielleicht ist ja ein billiger Anfängerfehler drin?

void BounceMesh::solveBeeman2()
{
	// Konstanten
	const float mass = m_Mass;
	const size_t n = m_iNumPoints;
	const float h = BOUNCE_TIMESTEP;

	// Für alle Punkte...
	for( size_t i=0 ; i<n ; ++i )
	{
		VEC3 F_prev = m_pMasspoints[i].force_old;
		VEC3 F_now  = m_pMasspoints[i].force;

		// Alte position merken
		VEC3 PosOld = m_pPositions[i];

		// Position berechnen:
		// x(t+h)
		m_pPositions[i] =	m_pPositions[i] +
							m_pMasspoints[i].velocity * h +
							((h*h) / (6*mass)) * (4*F_now - F_prev);

		// Geschwindigkeit, Beschleunigung und Kraft im nächsten Zeitschritt vorhersagen
		// Das kommt von mir selbst, stimmt das so?
		VEC3 V_next = (m_pPositions[i] - PosOld) / h;
		VEC3 A_next = V_next / h;
		VEC3 F_next = A_next * mass;

		// Geschwindigkeit berechnen:
		// v(t+h)
		m_pMasspoints[i].velocity = m_pMasspoints[i].velocity +
									(h/(6*mass)) * ( 2*F_next + 5*F_now - F_prev );

		// Kraft für nächsten Durchgang merken
		m_pMasspoints[i].force_old = m_pMasspoints[i].force;

		// Kraft zurücksetzen
		m_pMasspoints[i].force = VEC3(0,0,0);
	}
}

[lustig]

Wie findest du denn die Kraft für den aktuellen Zustand heraus?
Hast du nicht eine Funktion um die Kraft in Abhängigkeit vom Ort zu berechnen oder etwas in der Art?

Ich kenne den Algorithmus nicht, aber es scheint mir etwas seltsam, zuerst einen Wert für v_next ab zu schätzen und dann aus diesem den "richtigen" Wert für v_next zu berechnen. Kannst du nicht direkt die Kraft am neuen Ort herausfinden?
Wenn nicht solltest du dir trotzdem überlegen, ob deine Berechnung für f_next richtig ist. Du berücksichtigst nämlich keine "äusseren Kräfte", wenn du also z.B. eine langsam fliegende Masse im Gravitataionsfeld hast, wird das m*a wahrscheinlich kleiner sein als das m*g, oder schlimmstenfalls in eine ganz andere Richtung zeigen.

[0x00000001]

lustig schrieb:

Wie findest du denn die Kraft für den aktuellen Zustand heraus?
Hast du nicht eine Funktion um die Kraft in Abhängigkeit vom Ort zu berechnen oder etwas in der Art?

Die Kraft setzt sich aus mehrern Komponenten zusammen. Die Massepunkte sind untereinander mit Federn verbunden die die Punkte zu einem Körper verbinden. Außerdem wirkt noch Gravitation, und es werden Reibekräfte von Kollisionen mit anderen Objekten usw. dazukommen.
Aber das ist OK so wie es ist. Die Simulation funktioniert gut wenn ich zum Integrieren den Euler Algo nehme. Problem ist nur dass der explodiert wenn die Kräfte zu groß werden. Deshalb wollte ich mal Beeman austesten.

Da ich jetzt garnicht mehr weiterkomme werde ich wohl warten müssen und nach den Ferien den Prof fragen, denke mal der weiß wie's geht

[doppelmuffe]

hi,

0x00000001 schrieb:

Da F/m = a, hab ich die Beschleunigung im zukünftigen Zeitschritt über die Geschwindigkeit berechnet:

VEC3 V_next = (PosNew - PosOld) / h;
VEC3 A_next = V_next / h;

also für mich sieht das ziemlich falsch aus... V_next ist ok, aber die zweite zeile ist fragwürdig :xmas1:. mit

A_next = (V_next - V)/h

wäre ich schon eher zufrieden.

du kannst ja bem algorithmus probieren, F_next = F zu setzen. wenn's dann besser klappt, ist die extrapolation nicht so gut.

es wäre vermutlich auch angebracht, die differenzierbarkeit der kraftfunktion auszunutzen, um f_n+1 vorherzusagen. dazu müsstest du aber mehrere werte für F auf vorrat halten.

[0x00000001]

Danke doppelmuffe dass Du Dir das durchgelesen hast.
Ich war irgenwie der Meinung dass v=a*t, aber egal

Ich hab es jetzt mit Hilfe eines fremden Quelltextes zum Laufen bekommen, und ich glaube die machen es so wie Du vorgeschlagen hast. Es werden für jeden Punkt 3 Beschleunigungswerte mitgeführt. Der Trick ist vor allem, die Kräfte nicht vor der Simulation zu berechnen sondern mittendrin.

Letzenendes ist es aber eine Enttäuschung. Stabilität ist so gut wie garnicht dazugekommen, lediglich die Energieerhaltung funktioniert etwas besser. Ich denke ich bleibe bei Euler-Cromer. Der ist schneller und dafür kann ich dann den Zeitschritt etwas kleiner wählen was extrem viel bringt.

Danke an alle die geantwortet haben.

Zur Vollständigkeit noch die Lösung:

struct MASSPOINT
{
	// Positionen sind fürs Rendern gesondert gespeichert
	VEC3 velocity;
	VEC3 force;
	VEC3 acc;
	VEC3 acc1;
	VEC3 acc2;
};

void BounceMesh::solveBeeman()
{
	const float mass = m_Mass;
	const size_t n = m_iNumPoints;

	// Positionen
	for( size_t i=0 ; i<n ; ++i )
	{		
		// x(t+h)
		m_pPositions[i] =	m_pPositions[i] +
							BOUNCE_TIMESTEP * m_pMasspoints[i].velocity + 
							BOUNCE_TIMESTEP*BOUNCE_TIMESTEP * ( (2.0f/3.0f)*m_pMasspoints[i].acc - (1.0f/6.0f)*m_pMasspoints[i].acc1 );
	}

	// Kräfte zurücksetzen
	for( size_t i=0 ; i<n ; ++i )
	{
		m_pMasspoints[i].force = VEC3(0,0,0);
	}

	// Kräfte neu berechnen
	computeInternalForces(); // federn
	addGravity(); // schwerkraft
	collision(); // kollisionen

	// Bechleunigungen durchschieben
	for( size_t i=0 ; i<n ; ++i )
	{
		m_pMasspoints[i].acc2 = m_pMasspoints[i].acc1;
		m_pMasspoints[i].acc1 = m_pMasspoints[i].acc;
		m_pMasspoints[i].acc = m_pMasspoints[i].force / mass;
	}

	// Geschwindigkeiten
	for( size_t i=0 ; i<n ; ++i )
	{
		// v(t+h)
		m_pMasspoints[i].velocity = m_pMasspoints[i].velocity +
									BOUNCE_TIMESTEP * ( (5.0f/12.0f) * m_pMasspoints[i].acc +
														(2.0f/3.0f)  * m_pMasspoints[i].acc1 -
														(1.0f/12.0f) * m_pMasspoints[i].acc2 );
	}
}