Du bist nicht angemeldet.

Stilllegung des Forums
Das Forum wurde am 05.06.2023 nach über 20 Jahren stillgelegt (weitere Informationen und ein kleiner Rückblick).
Registrierungen, Anmeldungen und Postings sind nicht mehr möglich. Öffentliche Inhalte sind weiterhin zugänglich.
Das Team von spieleprogrammierer.de bedankt sich bei der Community für die vielen schönen Jahre.
Wenn du eine deutschsprachige Spieleentwickler-Community suchst, schau doch mal im Discord und auf ZFX vorbei!

implizites Runge Kutta??

Werbeanzeige

Wally

Frischling

  • »Wally« ist der Autor dieses Themas

Beiträge: 5

  • Private Nachricht senden

1

29.01.2008, 18:57

implizites Runge Kutta??

Abend erstmal,

also ich habe ein kleines Problem mit dem Runge-Kutta Verfahren. Ich frage mich ob das nicht auch implizit oder zumindest semi-implizit geht..

C-/C++-Quelltext

 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22

const double dt_by_2 = 0.5 * dt;
  const double dt_by_6 = (1./6.) * dt;

 
  double v_0     = v;                
  double accel_0 = a( x, v);  //dt wird nicht benötigt, da a() zeitdiskret ist

 
  double x_1     = x + dt_by_2*v_0;
  double v_1     = v + dt_by_2*accel_0;
  double accel_1 = a( x_1, v_1);
 
  double x_2     = x + dt_by_2*v_1;
  double v_2     = v + dt_by_2*accel_1;
  double accel_2 = a( x_2, v_2);
 
  // and at the end of the time interval

  double x_3     = x + dt*v_2;
  double v_3     = v + dt*accel_2;
  double accel_3 = a( x_3, v_3);
 
  v_neu = v + dt_by_6*(accel_0+accel_3+2.*(accel_1+accel_2));
  x_neu = x + dt_by_6*(v_0+v_3+2.*(v_1+v_2));


Sollte eigentlich richtig implementiert sein, oder bin ich zu blöd dafür? :confused:

Wenn a() jetzt eine gefederte Masse unter Erdschwerkraft beschreibt:

C-/C++-Quelltext

 
1
2
3
4
5
6

double a(double x, double v)
{
    double federkraft = 0.0;
    if (x < 1.0) federkraft = (1.0-federkraft)*FEDERHAERTE;
    return federkraft-9.81*MASSE;
}


dann wächst die Energie im System fast wie beim expliziten Euler-Verfahren.

Wie kann ich das nun vermeiden?

eine Idee wäre ja die letzte Zeile im Verfahren durch:

C-/C++-Quelltext

 
1

    x_neu = x + dt*v;



zu ersetzen, nur glaub ich da nicht dran, dass das die Lösung sein soll. :([/b]

Black-Panther

Alter Hase

Beiträge: 1 443

Wohnort: Innsbruck

  • Private Nachricht senden

2

29.01.2008, 19:30

Hallo!

Also ich hab mit dem Runge-Kutta verfahren keinerlei Probleme! Funktioniert alles so wie es soll... Hier mal mein Code:

C-/C++-Quelltext

 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66

const bool ogpRigidBody::Update(float dt)
{
    const float fHalfdt     = 0.5f * dt;
    const float fSixthdt    = dt / 6.f;
    og3DVector x(0.f), v(0.f), L(0.f), w(0.f), p(0.f);
    ogQuaternion q;
    ogMatrix mRot = ogMatrixIdentity();

    //InternalForces werden nur im ersten Schritt verwendet, da es sich um RepulsiveForces handelt!


    //A1 = G(t, S0), B1 = S0 + (dt / 2) * A1        => S. 238

    const og3DVector    A1dxdt = m_v;   //A1 dx/dt = v

    const ogQuaternion  A1dqdt = 0.5f * ogQuaternion(m_w) * m_q;
    const og3DVector    A1dpdt = m_vInternalForce + m_vExternalForce;
    const og3DVector    A1dLdt = m_vInternalTorque + m_vExternalTorque;
    x = m_x + fHalfdt * A1dxdt;
    q = m_q + fHalfdt * A1dqdt;
    p = m_p + fHalfdt * A1dpdt;
    L = m_L + fHalfdt * A1dLdt;
    ComputeDerivations(q, p, L, &mRot, &v, &w);

    //A2 = G(t + dt / 2, B1), B2 = S0 + (dt / 2) * A2

    const og3DVector    A2dxdt = v;
    const ogQuaternion  A2dqdt = 0.5f * ogQuaternion(w) * q;
    const og3DVector    A2dpdt = m_vExternalForce;
    const og3DVector    A2dLdt = m_vExternalTorque;
    x = m_x + fHalfdt * A2dxdt;
    q = m_q + fHalfdt * A2dqdt;
    p = m_p + fHalfdt * A2dpdt;
    L = m_L + fHalfdt * A2dLdt;
    ComputeDerivations(q, p, L, &mRot, &v, &w);

    //A3 = G(t + dt / 2, B2), B3 = S0 + dt * A3

    const og3DVector    A3dxdt = v;
    const ogQuaternion  A3dqdt = 0.5f * ogQuaternion(w) * q;
    const og3DVector    A3dpdt = m_vExternalForce;
    const og3DVector    A3dLdt = m_vExternalTorque;
    x = m_x + dt * A3dxdt;
    q = m_q + dt * A3dqdt;
    p = m_p + dt * A3dpdt;
    L = m_L + dt * A3dLdt;
    ComputeDerivations(q, p, L, &mRot, &v, &w);

    //A4 = G(t + dt, B3), S1 = S0 + (dt / 6) * (A1 + 2 * A2 + 2 * A3 + A4)

    const og3DVector    A4dxdt = v;
    const ogQuaternion  A4dqdt = 0.5f * ogQuaternion(w) * q;
    const og3DVector    A4dpdt = m_vExternalForce;
    const og3DVector    A4dLdt = m_vExternalTorque;
    m_x = m_x + fSixthdt * (A1dxdt + 2.f * (A2dxdt + A3dxdt) + A4dxdt);
    m_q = m_q + fSixthdt * (A1dqdt + 2.f * (A2dqdt + A3dqdt) + A4dqdt);
    m_p = m_p + fSixthdt * (A1dpdt + 2.f * (A2dpdt + A3dpdt) + A4dpdt);
    m_L = m_L + fSixthdt * (A1dLdt + 2.f * (A2dLdt + A3dLdt) + A4dLdt);
    //Rotationsquaternion normalisieren!

    m_q = m_q.Normalize(); //Wegen numerischen Fehler bei der Integration!

    ComputeDerivations();

    //Transformationsmatrizen berechnen und anwenden, dann AABB aktualisieren

    ComputeTransformationMatrices();
    UpdateAABBFast();

    //Nach diesem Frame Externe und Interne Kräfte wieder auf Null setzen

    m_vInternalForce = m_vInternalTorque = og3DVector(0.f);
    m_vExternalForce = m_vExternalTorque = og3DVector(0.f);

    return true;
}
stillalive studios
Drone Swarm (32.000 Dronen gleichzeitig steuern!)
facebook, twitter

Wally

Frischling

  • »Wally« ist der Autor dieses Themas

Beiträge: 5

  • Private Nachricht senden

3

29.01.2008, 19:58

Wäre gut zu wissen, was ComputeDerivations bei dir macht. (Das die Ableitungen berechnet werden ist mir klar, aber wie?).

p nehme ich an ist die Geschwindigkeit?!

Wäre nett, wenn du auf meinen Code mal kurz eingehen könntest, thx. :D

Schonmal danke für die Antwort.

p0llux

Treue Seele

Beiträge: 101

Wohnort: Aachen

Beruf: HiWi (theo. Inf.)

  • Private Nachricht senden

4

29.01.2008, 20:43

Eigentlich ist (zumindest in der Physik) p der Impuls. Wenn mich mein Abiwissen nicht verlassen hat ist p=m*v.
http://www.scienceandcode.org http://www.mathestudium.org

Black-Panther

Alter Hase

Beiträge: 1 443

Wohnort: Innsbruck

  • Private Nachricht senden

5

29.01.2008, 22:26

Ja, p ist der Impuls (linearer Impuls!)

ComputeDerivations macht folgendes:

C-/C++-Quelltext

 
1
2
3
4
5
6
7
8
9
10
11
12
13
14

inline void ComputeDerivations()
    {
        m_mRotation = m_q.RotationMatrix();
        m_v         = m_fInvMass * m_p;
        m_w         = ogUtilsTransformCoords(m_L, m_mRotation * m_mInvInertiaTensor * ogMatrixTranspose(m_mRotation));
        m_fSpeed    = og3DVectorLength(m_v);
    }
    inline void ComputeDerivations(ogQuaternion& q, const og3DVector& p, const og3DVector& L,
                                   ogMatrix* mRot, og3DVector* v, og3DVector* w)
    {
        *mRot       = q.RotationMatrix();
        *v          = m_fInvMass * p;
        *w          = ogUtilsTransformCoords(L, *mRot * m_mInvInertiaTensor * ogMatrixTranspose(*mRot));
    }


Hab mir jetzt deinen Code mal angesehn. Mir scheint du berechnest v und x richtig. Bei der Beschleunigung bin ich mir nicht sicher... Könntest mal testhalber versuchen es anstatt über a über p zu integriern!
stillalive studios
Drone Swarm (32.000 Dronen gleichzeitig steuern!)
facebook, twitter

Werbeanzeige