Błąd algorytmu czy ograniczenia operacji zmiennoprzecinkowych?

Witam. Wzdęło nie na napisanie symulatora ośrodka sprężystego (link: https://github.com/pirogronian/ElasticBodySimulator - proszę wybaczyć obecne tam pliki tymczasowe, robiłem repo pierwszy raz i na szybko). Rzecz wiąże się z pewnymi średnio zaawansowanymi obliczeniami na wektorach 3D. Rzecz tyczy się (najpewniej) symulatora, czyli klasy ElasticBodySimulator, zawartej w plikach https://github.com/pirogronian/ElasticBodySimulator/blob/master/Simulator/v1.0/ElasticBodySimulator.cpp oraz https://github.com/pirogronian/ElasticBodySimulator/blob/master/Simulator/v1.0/ElasticBodySimulator.h. Chodzi o to, że w miarę trwania symulacji energia układu rośnie, prawdopodobnie wykładniczo. Na obecnym etapie (po wyłapaniu kilku ewidentnych błędów) nie mam już żadnej pewności, czy jest to dalej błąd w moim algorytmie, czy też wynik niedoskonałości obliczeń zmiennoprzecinkowych.

Aby samodzielnie dokonać testu symulacji, należy po uruchomieniu programu kliknąć (nie ma tam nic innego) zakładkę Symulacja, następnie Nowa symulacja. Domyślnie układ składa się z sześciennej kostki. Rzutowanie na płaszczyznę odbywa się na stałe w poprzek osi X na zerowym indeksie, więc rezultatem będą cztery kropki na planie kwadratu. Z dostępnych obecnie opcji nowej symulacji istotne są Okres czasu - odcinek czasu na krok - jego zwiększenie bardzo przyspiesza przyrost energii, więc mam wrażenie, jakbym gdzieś w kodzie dał o jedno mnożenie za dużo; reszta powinna być jasna. Masa wiązania jest obecnie nieużywana.

Szczególnie podejrzane dla mnie są funkcje symulatora, odpowiedzialne za wyliczanie wektorów odległości oraz siły i położenia:

QVector3D inline distances(ElasticBodySimulator::Granule a, ElasticBodySimulator::Granule b) {
   QVector3D v;
    
   v.setX(b.position.x() - a.position.x());
   v.setY(b.position.y() - a.position.y());
   v.setZ(b.position.z() - a.position.z());
   
   return v;
}

QVector3D ElasticBodySimulator::forces(ElasticBodySimulator::Granule a, ElasticBodySimulator::Granule b) {
   QVector3D dis, fp;
    
   dis = distances(a, b);
    
   qreal f = _couple * (dis.length() - _eqspace);
    
   fp.setX(f * dis.x());
   fp.setY(f * dis.y());
   fp.setZ(f * dis.z());
    
   return fp;
}

void ElasticBodySimulator::updateForce(int i, int j, int k) {
   ElasticBodySimulator::Granule &g =  _gMatrix[i][j][k];

   g.force.setX(0);
   g.force.setY(0);
   g.force.setZ(0);
	
   if (i < _xSize - 1) g.force += forces(g, _gMatrix[i+1][j][k]);
   if (j < _ySize - 1) g.force += forces(g, _gMatrix[i][j+1][k]);
   if (k < _zSize - 1) g.force += forces(g, _gMatrix[i][j][k+1]);
	
   if (i > 0)  g.force += forces(g, _gMatrix[i-1][j][k]);
   if (j > 0)  g.force += forces(g, _gMatrix[i][j-1][k]);
   if (k > 0)  g.force += forces(g, _gMatrix[i][j][k-1]);
}

void ElasticBodySimulator::updatePosition(ElasticBodySimulator::Granule &g) {
   QVector3D accel, newspeed, movement;
	
   accel = g.force / granMass();
	
   newspeed = g.speed + accel * _timeSlice;
	
   movement = (g.speed + newspeed) / 2 * _timeSlice;
	
   g.position += movement;
   g.speed = newspeed;
}

Przeglądałem je ju jednak tyle razy, że jeśli faktycznie jest tam błąd, to po prostu jestem na niego ślepy. Będę wdzięczny za wszelkie sugestie.

Nie chce mi się analizować, ale z "efektu" wynika że zastosowałeś metodę Eulera, która zwykle wprowadza błąd systematyczny.
Poczytaj o metodzie Verleta albo LeapFrog.

Dzięki za info, jednakowoż po zapoznaniu się z rzeczonymi algorytmami stwierdziłem, że to raczej kosmetycznie poprawiony Euler. Jednak zaimplementowałem oba, z takim samym jak oryginał rezultatem (może oprócz tego, że przy tradycyjnym węzły odlatywały na koniec w jedną stronę, a przy leap frog - w drugą). Na wiki oraz w jednej pracy, którą znalazłem w necie piszą, że algorytm ten jest stabilny, ale pod paroma warunkami. Być może ich nie spełniłem. sprawdzę to w najbliższym czasie. Dzięki jeszcze raz - przynajmniej już wiem, gdzie szukać winnego.

polecam Verleta, jest naprawdę prosty w implementacji (wszystkie jego wady są dla ciebie nieistotne). Wystarczy ci taki wzór:

Zamiast pamiętać położenie i prędkość, masz pamiętać obecne i poprzednie położenie.

Z algorytmami Verleta jak i LeapFrog zapoznałem się na Wiki (stamtąd chyba nawet wziąłeś wzór), Jednak opisałem też, że sprawiły kłopoty. Dopiero po drugim twoim poście (o opatrzności!) jeszcze raz zaimplementowałem podstawową formę algorytmu i zauważyłem błąd w inicjacji przeszłego położenia. Zamiast inicjować obecnym położeniem, zerowałem je. Więc efektem było wyrzucenie węzłów w kosmos ;). Obecny kluczowy kod wygląda tak:

void ElasticBodySimulator::updatePosition(ElasticBodySimulator::Granule &g) {
	QVector3D accel, newspeed, movement, newpos;
	
//	accel = g.force / granMass();
	
//	newspeed = g.speed + accel * _timeSlice;
	
//	movement = (g.speed + newspeed) / 2 * _timeSlice;
	
	newpos = 2 * g.position - g.prevpos + g.force * timeSlice() * timeSlice() / granMass();
	
//	g.position += movement;
//	g.speed = newspeed;
	g.prevpos = g.position;
	g.position = newpos;
}

Teraz algorytm rzeczywiście działa zgodnie z oczekiwaniami i stabilnie. Jeszcze raz dzięki za poświęcony czas.