[c++] berekening klopt niet, op papier wel

[elgringo]

Ik heb een functie ala: H=2x^2+3x^2y^2-10py+1 (1)

De daadwerkelijke functie, volledig uitgeschreven:

code:

0.5*px*px+0.5*py*py+0.5*pz*pz+0.5*pa*pa+0.5*pb*pb+0.5*pc*pc+0.5*x*x*x*x-2*x*x*x*a+1*x*x*a*a+1*x*x*y*y-2*x*x*y*b+1*x*x*b*b+1*x*x*z*z-2*x*x*z*c+1*x*x*c*c-1*x*x+2*x*a*x*a-2*x*a*a*a-2*x*a*y*y+4*x*a*y*b-2*x*a*b*b-2*x*a*z*z+4*x*a*z*c-2*x*a*c*c+2*x*a+0.5*a*a*a*a+1*a*a*y*y-2*a*a*y*b+1*a*a*b*b+1*a*a*z*z-2*a*a*z*c+1*a*a*c*c-1*a*a+0.5*y*y*y*y-2*y*y*y*b+1*y*y*b*b+1*y*y*z*z-2*y*y*z*c+1*y*y*c*c-1*y*y+2*y*b*y*b-2*y*b*b*b-2*y*b*z*z+4*y*b*z*c-2*y*b*c*c+2*y*b+0.5*b*b*b*b+1*b*b*z*z-2*b*b*z*c+1*b*b*c*c-1*b*b+0.5*z*z*z*z-2*z*z*z*c+1*z*z*c*c-1*z*z+2*z*c*z*c-2*z*c*c*c+2*z*c+0.5*c*c*c*c-1*c*c+0.5

Deze functie sla ik op in 3 arrays:
Array 1 geeft de coefficient aan: [2,3,-10,1] (met fucntie boven genoemd bij (1) ) (double array)
Array 2 hoeveel vars er zijn: [2,4,1,0] (int array)
Array 3 welke vars [0,0,0,0,1,1,4] (0=x en 1=y, 3=px, 4=py) (deze array is de som van array 2 groot) (int array)

Verder heb ik een state vector, bij dit voorbeeld [statex,statey, statepx, statepy] (statex & statey zijn dus getallen, statepx & statepy zijn de impulsen, doubles)
En een afgeleide: [statexdot, stateydot, statepxdot, statepydot].

Dus elk punt heeft een coordinaat en een impuls. Stel dat je 2 dimensies hebt en 1 punt, dan heb je zoals boven staat. Het aantal states is dus het aantal dimensies*2.

Om de afgeleide te berkenen loop is alsvolgt door de 3 arrays heen:

C++:

memset(dStatesDot, 0x00,iNMax*sizeof(double));
            while (iWhichCoeff<iMaxTerms){     
                  if (dCoefficient[iWhichCoeff]==0) break;
                  for (unsigned int iLoop=0;iLoop<iNumberOfVariabeles[iWhichCoeff];iLoop++){  
                      dTemp=dCoefficient[iWhichCoeff];
                      for (unsigned int iLoop2=0;iLoop2<iNumberOfVariabeles[iWhichCoeff];iLoop2++){
                          if (iLoop!=iLoop2){
                             dTemp=dTemp*dStates[iVariabeles[iWhichVar+iLoop2]];

                          }
                      }
                      dStatesDot[iVariabeles[iWhichVar+iLoop]]+=dTemp;
                     
                  }
                  iWhichVar+=iNumberOfVariabeles[iWhichCoeff];
                  iWhichCoeff++;
            }

(dCoefficient is array 1, iNumberOfVariabeles is array 2 en iVariabeles is array 3, dStates is toestandsvector en dStatesDot de afgeleide hiervan)

Dit deel maakt van 2x^2 de afgeleide 2x+2x=4x
Dit wordt on the fly ingevuld. Nadat deze code klaar is staat in dStateDot de uitkomst van de afgeleide van H naar dState met de waarden alvast ingevuld.

dState is van het type (x1,y1, px1, py1, x2, y2, px2, py2). dus eerst de lokatie dan de impuls vervolgens het volgende punt.
dState wordt berekend met:

C++:

for (int iLoop=0;iLoop<iN;iLoop++){  
                for (int iLoop2=0;iLoop2<iDimension;iLoop2++){

                    dStates[iLoop*2*iDimension+iLoop2+iDimension]+=dDeltaT*(-dStatesDot[iLoop*2*iDimension+iLoop2]-dFriction*dStatesDot[iLoop*2*iDimension+iLoop2+iDimension]);                          
                    dStates[iLoop*2*iDimension+iLoop2]+=dDeltaT*(dStatesDot[iLoop*2*iDimension+iLoop2+iDimension]);      
                }
            }

Oftewel in het geval zoals bij functie (1) (irl is deze veel complexer):
dState[0]=dState[0]+dDeltaT*dStateDot[2] (x=x+dT*px)
dState[1]=dState[1]+dDeltaT*dStateDot[3] (y=y+dT*py)
dState[2]=dState[2]+dDeltaT*(-dStateDot[0] -dFriction*dStateDot[2]) (px=px+dT(-x-wrijving*px)
dState[3]=dState[3]+dDeltaT*(-dStateDot[1] -dFriction*dStateDot[3]) (py=py+dT(-y-wrijving*py)

Nu werkt dit proces prima als dFriction 0 is. (0<=dFriction<1). dDeltaT is de stapgrote in seconde, meestal 0.1

Als dFriction bijv 0.5 is zal dStateDot een afwijking gaat tonen. Deze zal in een situatie dat ie 0 moet bijven (zie de daadwerkelijke functie bovenin, met 2 nodes, 3 dimensies) langzaam groter dan 0 worden. Bij een stap wordt hij van 0 naar 1,4E-14, terwijl hij 0 moet blijven.

Waarom blijft ie niet nul? is dit een afrondingsfout oid?
Ik hoop dat alles duidelijk is

[ Voor 9% gewijzigd door elgringo op 18-06-2006 21:32 ]

[Robtimus]

Grote kans dat dit idd komt door afrondingsfouten.

Probeer het volgende stuk code maar eens:

C:

float f = 0.0;
while (f < 5) {
    printf("f is now %f\n", f);
    f = f + 0.1;
}

Je zult zien dat je niet 0.0, 0.1, 0.2, 0.3, ... zult zien maar na een tijdje bv 0.8999999 of 0.90000001. Dit komt omdat floating point getallen vaak niet precies kunnen worden opgeslagen, dat kan alleen met getallen in de vorm van a*(1/2) + b*(1/4) + c*(1/8) + .....

Edit: even getest, hier komt het eerst voor met 2.799999 en later ook 3.899998.

[ Voor 17% gewijzigd door Robtimus op 18-06-2006 21:34 ]

[Daos]

Dat is inderdaad een afronding fout. In plaats van vergelijken met nul moet je kijken of de absolute waarde kleiner is dan bv 1e-5.

[elgringo]

Daos schreef op zondag 18 juni 2006 @ 21:31:
Dat is inderdaad een afronding fout. In plaats van vergelijken met nul moet je kijken of de absolute waarde kleiner is dan bv 1e-5.

Maar ik vergelijk niets met 0. Hoe kan ik dit gaan afronden?

[Daos]

elgringo schreef op zondag 18 juni 2006 @ 21:36:
[...]


Maar ik vergelijk niets met 0. Hoe kan ik dit gaan afronden?

Er komen vaak topics langs over berekeningen die helemaal fout gaan door het vergelijken van floatingpointgetallen met == en !=. Ik maakte uit je startpost op dat zoiets ook hier gebeurt. Of heeft alleen eindresultaat een afwijkinkje?

[elgringo]

Daos schreef op zondag 18 juni 2006 @ 21:42:
[...]


Er komen vaak topics langs over berekeningen die helemaal fout gaan door het vergelijken van floatingpointgetallen met == en !=. Ik maakte uit je startpost op dat zoiets ook hier gebeurt. Of heeft alleen eindresultaat een afwijkinkje?

wordt niets vergelijken, alleen berekend. De uitkomst heeft dus een miniscule afwijking. Maar omdat hiermee de volgende stap berekend wordt gaat het mis.

Lees dat wordt 1 bijvoorbeeld verkeerd afgerond waarna hij van bijna nul naar oneindig gaat....

[ Voor 9% gewijzigd door elgringo op 18-06-2006 21:56 ]

[Daos]

elgringo schreef op zondag 18 juni 2006 @ 21:45:
[...]


wordt niets vergelijken, alleen berekend. De uitkomst heeft dus een miniscule afwijking. Maar omdat hiermee de volgende stap berekend wordt gaat het mis.

Heb je ooit wel eens van significante cijfers gehoord? Je heb hier altijd mee te maken als je iets numeriek op gaat lossen of benaderen.

Als je exact wil rekenen, dan kan je proberen een symbolisch rekenprogramma (zoals Maple) te gebruiken.

[elgringo]

Daos schreef op zondag 18 juni 2006 @ 21:58:
[...]

Heb je ooit wel eens van significante cijfers gehoord? Je heb hier altijd mee te maken als je iets numeriek op gaat lossen of benaderen.

Als je exact wil rekenen, dan kan je proberen een symbolisch rekenprogramma (zoals Maple) te gebruiken.

Ja ik heb vansignificante cijfers gehoord. Ik verwacht niet dat iets wat naar 0 gaat naar oneindig gaat door een afrondingsfout.

Maple is hiet niet voor geschikt. De functie H bevat meer dan 1000 coefficienten. Maple zou dat dagen moeten rekenen voordat er een uitkomst komt

[windancer]

Het probleem hier is de machine preciesie. Reele getallen kunnen niet altijd in een eindige hoeveelheid geheugen worden opgeslagen. In het binaire stelsel dat door de computer wordt gebruikt kan bijvoorbeel 3/4 opgeslagen worden als (1/2)^(-1)+(1/2)^(-2) maar 1/10 kan bijvoorbeeld niet precies worden opgeslagen.

Precies hetzelfde probleem heb je in het decimale stelsel. Als je 1/3 uit wil schrijven krijg je 0.33333333333333 en zo voort. Het maakt niet uit hoe ver je door gaat, je maakt altijd een afbreekfout. In dit geval is die ook ongeveer 10^(-14), net als in jouw geval.

[elgringo]

windancer schreef op zondag 18 juni 2006 @ 22:15:
Het probleem hier is de machine preciesie. Reele getallen kunnen niet altijd in een eindige hoeveelheid geheugen niet worden opgeslagen. In het binaire stelsel dat door de computer wordt gebruikt kan bijvoorbeel 3/4 opgeslagen worden als (1/2)^(-1)+(1/2)^(-2) maar 1/10 kan bijvoorbeel niet precies worden opgeslagen.

Precies hetzelfde probleem heb je in het decimale stelsel. Als je 1/3 uit wil schrijven krijg je 0.33333333333333 en zo voort. Het maakt niet uit hoe ver je door gaat, je maakt altijd een afbreekfout. In dit geval is die ook ongeveer 10^(-14), net als in jouw geval.

Dat stap ik, maar ik ben zoiets aan het uitrekenen als y=1/x. deze moet naar 0 gaan en bij mij schiet hij opeens naar oneindig.

Hoe kan ik dit oplossen?

[D-Raven]

Ik weet niet of het in C++ mogelijk is om bij een float het aantal significante cijfers op te geven. (Lijkt mij van niet, maar zovaak speel ik daar niet mee dus ik kan er naast zitten).
Maar anders zou je zelf een klasse kunnen schrijven die een float representeerd met een in te stellen aantal significante cijfers..

[MisterData]

Gebruik anders eens een double in plaats van een float, als je probleem dan minder wordt dan weet je in ieder geval dat je het in de precisie moet gaan zoeken

[GarBaGe]

elgringo schreef op zondag 18 juni 2006 @ 22:18:
Dat stap ik, maar ik ben zoiets aan het uitrekenen als y=1/x. deze moet naar 0 gaan en bij mij schiet hij opeens naar oneindig.

Hoe kom je nu aan "oneindig" ??
In je eerste post geef je aan "1,4E-14" ipv nul
Nu is 1,4E-14 niet oneindig, maar BIJNA nul, namelijk:
0,000000000000014
Lijkt me verre van oneindig en toch echt wel bijna nul....

[Michali]

elgringo schreef op zondag 18 juni 2006 @ 22:18:
Dat stap ik, maar ik ben zoiets aan het uitrekenen als y=1/x. deze moet naar 0 gaan en bij mij schiet hij opeens naar oneindig.

Naast wat GarBaGe zegt, ligt dit ook maar net aan wat voor waarde je aan X geeft. Als X oneindig klein wordt, dan wordt Y oneindig groot. Maar dat was niet het probleem hier.

[H!GHGuY]

je kan natuurlijk ook proberen om met breuken te gaan werken...

class Breuk en dan alle operatoren overloaden.

[Daos]

D-Raven schreef op maandag 19 juni 2006 @ 00:33:
Ik weet niet of het in C++ mogelijk is om bij een float het aantal significante cijfers op te geven. (Lijkt mij van niet, maar zovaak speel ik daar niet mee dus ik kan er naast zitten).
Maar anders zou je zelf een klasse kunnen schrijven die een float representeerd met een in te stellen aantal significante cijfers..

Hij gebruikt al een double. In een float past een getal met ongeveer 6-7 significante cijfers en in een double past een getal met ongeveer 14-15 significante cijfers.

[elgringo]

MisterData schreef op maandag 19 juni 2006 @ 08:15:
Gebruik anders eens een double in plaats van een float, als je probleem dan minder wordt dan weet je in ieder geval dat je het in de precisie moet gaan zoeken

zoals in ts staat gebruik ik doubles (long doubles wordt het beter, +/- 70 stappen later pas oneindig (met doubles na 830 stappen, long double ongeveer 900))

De afronding van de double die geen 0 is maar net iets meer heeft uiteindelijk het gevolg dat hij oneindig wordt, aangezien hier mee doorgerekend wordt en elke node van een andere afhankelijk is

breuken zijn ook geen oplossing. Ik wil gewoon dat hij een significante heeft < de double significantie, zodat ik deze fouten er uit kan halen. Of zijn er andere manieren om deze nummerieke fouten er uit te halen

[ Voor 20% gewijzigd door elgringo op 21-06-2006 12:42 ]

[elgringo]

Ik heb gevonden waar het fout gaat.

Na 749 stappen gerekend te hebben staat 1 impuls op 0 en de andere nog niet. Deze andere had ook op nul moeten staan, maar waarschijnlijk door wat afrondingsfouten is dat niet het geval.
doordat de ene wel op 0 staat en de andere nog niet wordt de eerste nagatief en dan gaat het mis.
Hoe kan dit opgelost worden?

[MSalters]

Beter numeriek algoritme. Als jij weet dat beide impulsen tegelijk 0 worden, zorg ervoor dat je software die randvoorwaarde dan ook gebruikt. En het allerbelangrijkste: Singulariteiten zijn punten als je met oneindige precisie rekent, maar intervallen als je met beperkte precisie rekent. Zorg er dus voor dat je uit de buurt blijft.

[elgringo]

MSalters schreef op zondag 02 juli 2006 @ 15:59:
Beter numeriek algoritme. Als jij weet dat beide impulsen tegelijk 0 worden, zorg ervoor dat je software die randvoorwaarde dan ook gebruikt. En het allerbelangrijkste: Singulariteiten zijn punten als je met oneindige precisie rekent, maar intervallen als je met beperkte precisie rekent. Zorg er dus voor dat je uit de buurt blijft.

Maar in snap niet waar het verschil van beide vandaan komt
beide worden met de zelfde formule uitgerekend allen is de volgorde net iets anders (dus bij de ene A+B en bij de andere B+A (maar dan met een x aantal termen))

[.oisyn]

Het probleem is dat floating point berekeningen niet associatief zijn. Dus A+(B+C) hoeft niet dezelfde uitkomst te hebben als (A+B)+C. Dit is logisch als je bedenkt dat de getallen worden opgeslagen dmv een aantal significante cijfers en een exponent. Stel er wordt in het decimale talstelsel gewerkt en er zijn 4 significante cijfers, en je wilt 12340+3+4 berekenen

12340 wordt opgeslagen als 1234e1, 3 als 3e0 en 4 als 4e0
Voor (12340+3)+4:
1234e1 + 3 = 12343 -> naar 4 significante cijfers: 1234e1
1234e1 + 4 = 12344 -> naar 4 significante cijfers: 1234e1

Nu voor 12340+(3+4)
3+4 = 7 -> 7e0
1234e1 + 7 = 12347 -> 1235e1

Verschillende uitkomsten, omdat je na verschillende berekeningen afrondt.

[ Voor 8% gewijzigd door .oisyn op 02-07-2006 18:16 ]