# 2D Rigidbody Physics Engine Bug

• ### Question

• [OT]

Thorsten, Tommy,

I've been working on another version of a library for gaming and general 2D graphics in windows forms.  It uses a 2D rigid body physics system based on the work (articles and C++ samples) located here.

I've almost got it working but there are still some bugs either in the translation and/or my implementation.  The physics works, but there can be undesired inertia build up resulting in too much force being applied to a body.

Here's a screen shot:

If I put this out on GitHub would you guys be interested in contributing and helping to resolve the bugs?  I could certainly use help on some of the math!

Reed Kimble - "When you do things right, people won't be sure you've done anything at all"

• Edited by Sunday, November 8, 2015 11:10 PM forgot link
• Split by Wednesday, November 11, 2015 3:40 PM side topic evolved into thread
Sunday, November 8, 2015 11:09 PM

• Hi,

My question is: Are you sure about the plus sign in (computing the relative velocity), method ApplyImpulse

rv = (B.velocity + MathF.Cross(B.angularVelocity, rb)) - (A.velocity - MathF.Cross(A.angularVelocity, ra))

Changing it to minus (= subtracting the "angular components" from the velocity vector for both shapes, since ra and rb will already contain the signs) at the two code lines where it appears makes the system stable immediately...

And both minus (b.vel - ... and a.vel - ...) (not both plus), because of viewing from A to B (standing on the edge of A looking at the mass-point of B moving by -> standing on earth watching a shootingstar)

(rv = (B.velocity - MathF.Cross(B.angularVelocity, rb)) - (A.velocity - MathF.Cross(A.angularVelocity, ra)))

[I came to this because some rotating objects will get increased angular velocity when colliding with other objects, but by "viewing" the collision, it should get a decreased angular velocity (or break into pieces)]

In the picture below the rotating shape should lower its av when colliding, but the opposite is the case...

Regards,

Thorsten

Wednesday, November 11, 2015 10:10 AM

### All replies

• Reed.

I will give it a try, as long as I am not expected to understand it. :)

Sunday, November 8, 2015 11:48 PM
• If I put this out on GitHub would you guys be interested in contributing and helping to resolve the bugs?  I could certainly use help on some of the math!

Hi,

yes, certainly I'll have a look, if I find the time...

The physics works, but there can be undesired inertia build up resulting in too much force being applied to a body

Do you use the Euler-method for solving the differential equations maybe? It generally tends to gain energy for a system. Doing a pendulum that way will probably break down all surrounding walls after some time... There's a trick for this, using a symplectic Euler method...

https://en.wikipedia.org/wiki/Semi-implicit_Euler_method

But thats just a wild guess here...

Regards,

Thorsten

Monday, November 9, 2015 1:40 AM
• If I put this out on GitHub would you guys be interested in contributing and helping to resolve the bugs?  I could certainly use help on some of the math!

Hi,

yes, certainly I'll have a look, if I find the time...

The physics works, but there can be undesired inertia build up resulting in too much force being applied to a body

Do you use the Euler-method for solving the differential equations maybe? It generally tends to gain energy for a system. Doing a pendulum that way will probably break down all surrounding walls after some time... There's a trick for this, using a symplectic Euler method...

https://en.wikipedia.org/wiki/Semi-implicit_Euler_method

But thats just a wild guess here...

Regards,

Thorsten

Thorsten,

That's amazing that you know that!

They describe Eulers here in an branch of the link Reed gave above.

Monday, November 9, 2015 3:14 AM
• ..

Thorsten,

That's amazing that you know that!

...

Thanks Tommy :-)

I had some good online professors when learning the basics.

[in german]

http://www.j3l7h.de/videos.html

[basics in english I & II]

and of course the courses of Leonard Susskind "The theoretical Minimum"

Regards,

Thorsten

Monday, November 9, 2015 4:51 AM
• https://en.wikipedia.org/wiki/Semi-implicit_Euler_method

But thats just a wild guess here...

Sounds like a wildly accurate guess!  Maybe that's not the whole issue but it sounds like it may indeed be part of it.  I'll check that link and compare against the code.

I was trying to get the physics portion working before incorporating all of the other game engine objects into the project but I may go ahead and get everything together (with the currently broken physics) just to have a mostly-complete project when making the initial GitHub check-in.  I'll post back here with a link when the repo is ready.

@Tommy: you are welcome to contribute in whatever way you can.  I think you sell yourself short sometimes!  :)  Also, on another note, you should really check out Babylon.js as it could be the perfect stepping stone to a hardware accelerated version of CadRail.  I saw a session on it at Summit and was very impressed by what you can do with it and how easy it was to use.

Reed Kimble - "When you do things right, people won't be sure you've done anything at all"

Monday, November 9, 2015 4:57 PM
• https://en.wikipedia.org/wiki/Semi-implicit_Euler_method

But thats just a wild guess here...

Sounds like a wildly accurate guess!  Maybe that's not the whole issue but it sounds like it may indeed be part of it.  I'll check that link and compare against the code.

I was trying to get the physics portion working before incorporating all of the other game engine objects into the project but I may go ahead and get everything together (with the currently broken physics) just to have a mostly-complete project when making the initial GitHub check-in.  I'll post back here with a link when the repo is ready.

@Tommy: you are welcome to contribute in whatever way you can.  I think you sell yourself short sometimes!  :)  Also, on another note, you should really check out Babylon.js as it could be the perfect stepping stone to a hardware accelerated version of CadRail.  I saw a session on it at Summit and was very impressed by what you can do with it and how easy it was to use.

Reed Kimble - "When you do things right, people won't be sure you've done anything at all"

Reed,

I will check out Babylon.js. If it is anything like Babylon 5 then I should understand it. :)

The write up on the basic engine in the link you gave describes what Thorsten said and apparently it is using the Symplectic but hard to tell for sure from the article:

This will work, and is commonly used as a starting point. However, it has numerical inaccuracies that we can get rid of without any extra effort. Here is what is known as Symplectic Euler:

 1 2 3 `// Symplectic Euler` `v += (1/m * F) * dt` `x += v * dt`

http://gamedevelopment.tutsplus.com/tutorials/how-to-create-a-custom-2d-physics-engine-the-core-engine--gamedev-7493

Monday, November 9, 2015 9:07 PM
• Reed,

I will check out Babylon.js. If it is anything like Babylon 5 then I should understand it. :)

The write up on the basic engine in the link you gave describes what Thorsten said and apparently it is using the Symplectic but hard to tell for sure from the article:

This will work, and is commonly used as a starting point. However, it has numerical inaccuracies that we can get rid of without any extra effort. Here is what is known as Symplectic Euler:

 1 2 3 `// Symplectic Euler` `v += (1/m * F) * dt` `x += v * dt`

http://gamedevelopment.tutsplus.com/tutorials/how-to-create-a-custom-2d-physics-engine-the-core-engine--gamedev-7493

The creator chose the name because he loved Babylon 5, so you're in luck.  :)  (that's a true story too!)

Thanks for finding that in the article.  The biggest issue with this article series is that the GitHub repo with the code is modeled quite differently than the code samples throughout the article - the author renamed, reorganized and refactored the code so the project and article are not in sync.  But now that I know what to look for I may be able to spot a fix in the differences between the code in the project and the code in the article.

Reed Kimble - "When you do things right, people won't be sure you've done anything at all"

Tuesday, November 10, 2015 5:42 PM
• The creator chose the name because he loved Babylon 5, so you're in luck.  :)  (that's a true story too!)

Thanks for finding that in the article.  The biggest issue with this article series is that the GitHub repo with the code is modeled quite differently than the code samples throughout the article - the author renamed, reorganized and refactored the code so the project and article are not in sync.  But now that I know what to look for I may be able to spot a fix in the differences between the code in the project and the code in the article.

Reed Kimble - "When you do things right, people won't be sure you've done anything at all"

The article also mentions that some dont save the mass as an object property, they save 1/m for the object property. One reason is thats what appears in the equations. So I got thinking if you calc 1/m in a loop 100,000 times and it has been rounded up by some amount mr then you have an extra 100000(mr) of mass in the system?

I dont know how to aviod it but just something I was thinking about, since you brought it up. If you start with 1/m then it does not get rounded. But then I realized the saved value 1/m would already be rounded as well so no difference?

Tuesday, November 10, 2015 9:59 PM
• The creator chose the name because he loved Babylon 5, so you're in luck.  :)  (that's a true story too!)

Thanks for finding that in the article.  The biggest issue with this article series is that the GitHub repo with the code is modeled quite differently than the code samples throughout the article - the author renamed, reorganized and refactored the code so the project and article are not in sync.  But now that I know what to look for I may be able to spot a fix in the differences between the code in the project and the code in the article.

Reed Kimble - "When you do things right, people won't be sure you've done anything at all"

The article also mentions that some dont save the mass as an object property, they save 1/m for the object property. One reason is thats what appears in the equations. So I got thinking if you calc 1/m in a loop 100,000 times and it has been rounded up by some amount mr then you have an extra 100000(mr) of mass in the system?

I dont know how to aviod it but just something I was thinking about, since you brought it up. If you start with 1/m then it does not get rounded. But then I realized the saved value 1/m would already be rounded as well so no difference?

I checked this out a little and the difference seems insignificant.

```Public Class Form3

Dim total1 As Double
Dim m As Double = 2.3
Dim loops As Double = 100000

For i = 1 To loops
total1 += (1 / m)
Next

Dim total2 As Double
total2 = loops * (1 / m)

Dim dt As Double = total1 - total2

MsgBox(total1.ToString & vbLf & total2.ToString & vbLf & dt.ToString)

End Sub
End Class```

Wednesday, November 11, 2015 12:22 AM

• ```Public Class Form3

Dim total1 As Double
Dim m As Double = 2.3
Dim loops As Double = 100000

For i = 1 To loops
total1 += (1 / m)
Next

Dim total2 As Double
total2 = loops * (1 / m)

Dim dt As Double = total1 - total2

MsgBox(total1.ToString & vbLf & total2.ToString & vbLf & dt.ToString)

End Sub
End Class```

Hi,

when seeing the minus sign here, I remember another possible issue, when computing numbers: [catastrophic] Cancellation.

http://en.wikipedia.org/wiki/Loss_of_significance

I mentioned it once here:

https://social.msdn.microsoft.com/Forums/vstudio/en-US/d27f9f72-dda9-482d-b55c-25229c92ccf6/the-big-failure-of-the-floating-point-computation?forum=vbgeneral

So this also could be an issue when values tend to go away from the ideal result... I have to refresh my Numerical maths, but I think we could estimate things via the Matrix-Norms...

Regards,

Thorsten

Wednesday, November 11, 2015 2:56 AM
• In the actual code implementation, the value for mass is stored as 1/m so it is only calculated once when the property is set. The inverse value is then used throughout the calculations.

I've put a copy of the project in a zip file on OneDrive if you want to check it out.  I haven't had a chance to do the whole GitHub thing yet, but hopefully I'll get to it in the next couple of days.

Reed Kimble - "When you do things right, people won't be sure you've done anything at all"

Wednesday, November 11, 2015 4:14 AM
• In the actual code implementation, the value for mass is stored as 1/m so it is only calculated once when the property is set. The inverse value is then used throughout the calculations.

I've put a copy of the project in a zip file on OneDrive if you want to check it out.  I haven't had a chance to do the whole GitHub thing yet, but hopefully I'll get to it in the next couple of days.

Reed Kimble - "When you do things right, people won't be sure you've done anything at all"

Ok, thanks.

I downloaded it and will have a look today or tomorrow. From a first view, I think, this is really good!!

I'll do some testing in the next days to lookup, why, when a lot of polygons or circels are added, the energy will increase for single shapes (getting kicked out of the "pot" when not expecting it)

Regards,

Thorsten

Wednesday, November 11, 2015 4:49 AM
• Hi,

My question is: Are you sure about the plus sign in (computing the relative velocity), method ApplyImpulse

rv = (B.velocity + MathF.Cross(B.angularVelocity, rb)) - (A.velocity - MathF.Cross(A.angularVelocity, ra))

Changing it to minus (= subtracting the "angular components" from the velocity vector for both shapes, since ra and rb will already contain the signs) at the two code lines where it appears makes the system stable immediately...

And both minus (b.vel - ... and a.vel - ...) (not both plus), because of viewing from A to B (standing on the edge of A looking at the mass-point of B moving by -> standing on earth watching a shootingstar)

(rv = (B.velocity - MathF.Cross(B.angularVelocity, rb)) - (A.velocity - MathF.Cross(A.angularVelocity, ra)))

[I came to this because some rotating objects will get increased angular velocity when colliding with other objects, but by "viewing" the collision, it should get a decreased angular velocity (or break into pieces)]

In the picture below the rotating shape should lower its av when colliding, but the opposite is the case...

Regards,

Thorsten

Wednesday, November 11, 2015 10:10 AM
• In the actual code implementation, the value for mass is stored as 1/m so it is only calculated once when the property is set. The inverse value is then used throughout the calculations.

I've put a copy of the project in a zip file on OneDrive if you want to check it out.  I haven't had a chance to do the whole GitHub thing yet, but hopefully I'll get to it in the next couple of days.

Reed Kimble - "When you do things right, people won't be sure you've done anything at all"

I have VS 2013 which does not include .net 4.6 and it appears 4.6 requires VS 2015? I tried 4.5.1 but it has a few statements not recognized. Do I need VS 2015 to run this project? I am also on a Win 7 system although I have win 10 on another system and could install VS 2015 there.

Wednesday, November 11, 2015 3:26 PM
• Hi,

I just experimented a bit with it, by changing some signs and values. My question is: Are you sure about the plus sign in (computing the relative velocity), method ApplyImpulse

rv = (B.velocity + MathF.Cross(B.angularVelocity, rb)) - (A.velocity - MathF.Cross(A.angularVelocity, ra))

Changing it to minus (= subtracting the angular momentum from the velocity vector for both shapes, since ra and rb will already contain the signs) at the two code lines where it appears makes the system stable immediately...

I'm sure that's how it was written in the C code I translated... but I'm not at all sure it was correct to begin with!

You are right in that changing the plus to minus seems to correct the issue.  Maybe that was the whole problem - a typo in the source.

Reed Kimble - "When you do things right, people won't be sure you've done anything at all"

Wednesday, November 11, 2015 3:26 PM

• I'm sure that's how it was written in the C code I translated... but I'm not at all sure it was correct to begin with!

You are right in that changing the plus to minus seems to correct the issue.  Maybe that was the whole problem - a typo in the source.

Reed Kimble - "When you do things right, people won't be sure you've done anything at all"

Hi,

I think so, since just putting in two shapes, both with initial angular velocity of 1, placing the first almost at the top of the form, the second somewhat below so that the first will hit the second one during falling, and looking at the values of ra and rb will result IMHO in a correct relative velocity...

And the gain of angular velocity - good to see when balls are surrounded by polygons - will just vanish...

Regards,

Thorsten

Wednesday, November 11, 2015 3:31 PM
• I just changed the project to target 4.5.1 and it compiled and ran OK.  What line(s) did you get errors on?  You would need 2015 to run 4.6, but shouldn't need Windows 10 (though I highly recommend both VS2015 and Win10).

Reed Kimble - "When you do things right, people won't be sure you've done anything at all"

Wednesday, November 11, 2015 3:37 PM
• I have VS 2013 which does not include .net 4.6 and it appears 4.6 requires VS 2015? I tried 4.5.1 but it has a few statements not recognized. Do I need VS 2015 to run this project? I am also on a Win 7 system although I have win 10 on another system and could install VS 2015 there.

Hi, I've changed the target framework to 4.5 for both projects. It runs without complaining anything...

[You always can replace the async await part by implementing a control thread yourself, you know that I'm sure...]

Regards,

Thorsten

Wednesday, November 11, 2015 3:37 PM
• I have VS 2013 which does not include .net 4.6 and it appears 4.6 requires VS 2015? I tried 4.5.1 but it has a few statements not recognized. Do I need VS 2015 to run this project? I am also on a Win 7 system although I have win 10 on another system and could install VS 2015 there.

Hi, I've changed the target framework to 4.5 for both projects. It runs without complaining anything...

[You always can replace the async await part by implementing a control thread yourself, you know that I'm sure...]

Regards,

Thorsten

When I open the project it gives me the choice to switch to .net 4.5.0 which I do and it shows the errors. Then I change it to 4.5.1 and still the same errors.

Then I added the suggested end properties and a couple other suggestions and the clock errors went away but it still says duration is not defined with no suggestions.

These are the original errors:

So it may be a reverse conversion problem of some kind from 2015 to 2013? I am not sure how to fix it up but you guys can probably see the problem in the first image?

I have been wanting to try 2015 so will install it.

Sounds like you have found the problem. Atta boy Thorsten!

Wednesday, November 11, 2015 3:58 PM
• Guys, I split this out into its own thread since it has evolved beyond the scope of a comment in someone else's thread.  Thorsten also deserved an "answer" for finding that bug in the code!  :)

Reed Kimble - "When you do things right, people won't be sure you've done anything at all"

Wednesday, November 11, 2015 4:00 PM
• Hi,

something to think about/ play with:

Increase the gravity by a factor of 10 -> so the value is (0;980)

Give balls some more friction and some restitution:

```    Private Sub Form1_MouseClick(sender As Object, e As MouseEventArgs) Handles Me.MouseClick
If e.Button = MouseButtons.Left Then
Dim poly As New Polygon
'poly.SetBox(MathF.Random(6, 24), MathF.Random(6, 24))
Dim count = MathF.Random(3, Polygon.MaxPolyVertexCount) '
Dim vertices(count - 1) As Vector2
'vertices = MathF.GenerateEvenPolygon(MathF.Random(16, 32), count, New Vector2(e.X, e.Y), False)
Dim ee = MathF.Random(25, 55)
For i = 0 To count - 1
vertices(i).Set(MathF.Random(-ee, ee), MathF.Random(-ee, ee))
Next

poly.Set(vertices, count)
Dim b As New Rigidbody(poly, e.X, e.Y)
'b.density = 0.01
b.SetOrient(MathF.Random(-MathF.PI, MathF.PI))
b.restitution = 1.5F
b.dynamicFriction = 0.1F
b.staticFriction = 0.2F

b.ApplyInitialImpulse()

ElseIf e.Button = MouseButtons.Right
Dim ballShape As New Circle(MathF.Random(16, 32))
Dim ball As New Rigidbody(ballShape, e.X, e.Y)
ballShape.SetOrient(MathF.Random(-MathF.PI, MathF.PI))
'ball.density = 20
ball.restitution = 1.85F
ball.dynamicFriction = 0.2F
ball.staticFriction = 0.25F
End If
End Sub```

...

and, of course, make the bottom vibrate :-)

Regards,

Thorsten

Wednesday, November 11, 2015 4:11 PM
• So I added the gets to the top of the class but still get the errors shown on the durations.

Just FYI.

Wednesday, November 11, 2015 4:18 PM

• When I open the project it gives me the choice to switch to .net 4.5.0 which I do and it shows the errors. Then I change it to 4.5.1 and still the same errors.

Hi Tommy,

this looks to me like a problem with invisible/control chars in the file, or an error somewhere else in the file...

Those oneline property declarations should work without problems in fw 4.0/VS 2010...

Update: I tested with VS2010  the class clock. It gave me the same errors. When removing the word ReadOnly it works... (So, if you want to have the properties readonly, you need to full-implement them)

Regards,

Thorsten

Wednesday, November 11, 2015 4:20 PM
• So I added the gets to the top of the class but still get the errors shown on the durations.

...

The _Duration and _FixedDuration variables are the automatically created fields by the property declarations in that one-line way.

Regards,

Thorsten

Wednesday, November 11, 2015 4:24 PM
• ReadOnly auto-properties are new to VS2015; in earlier versions we only have Read/Write auto-properties.

LOL I changed the WFG library to 4.5 but not the example project.  :P  That's where the compatibility issues are at.

You'll have to fully implement those readonly properties; as Thorsten mentioned, you need to declare the backing fields (the property name with an underscore) and return them from the property accessor.

Reed Kimble - "When you do things right, people won't be sure you've done anything at all"

Wednesday, November 11, 2015 4:30 PM
• ReadOnly auto-properties are new to VS2015; in earlier versions we only have Read/Write auto-properties.

LOL I changed the WFG library to 4.5 but not the example project.  :P  That's where the compatibility issues are at.

You'll have to fully implement those readonly properties; as Thorsten mentioned, you need to declare the backing fields (the property name with an underscore) and return them from the property accessor.

Reed Kimble - "When you do things right, people won't be sure you've done anything at all"

Well I think I did what you describe? Now I get this which I cant figure what it means.

PS no errors are shown now with the adds.

```Public Class Clock
Private _FixedDuration As Double
Private _Duration As Double

Public ReadOnly Property Duration As Double
Get
Return _Duration
End Get
End Property
Public ReadOnly Property FixedDuration As Double
Get
Return _FixedDuration
End Get
End Property

Public ReadOnly Property FixedUpdateTime As Double
Get
Return lastFixed
End Get
End Property
Public ReadOnly Property UpdateTime As Double
Get
Return lastUpdate
End Get
End Property
Private lastUpdate As Double
Private lastFixed As Double
Private sw As New Stopwatch
Private swf As New Stopwatch
Public Sub New()
sw.Start()
swf.Start()
End Sub

Public Sub FixedUpdate()
swf.Stop()
lastFixed = swf.Elapsed.TotalSeconds
_FixedDuration += lastFixed
swf.Restart()
End Sub

Public Sub Update()
sw.Stop()
lastUpdate = sw.Elapsed.TotalSeconds
_Duration += lastUpdate
sw.Restart()
End Sub
End Class```

Wednesday, November 11, 2015 4:43 PM

• Well I think I did what you describe? Now I get this which I cant figure what it means.

PS no errors are shown now with the adds.

set the other project as StartupProject

Regards,

Thorsten

Wednesday, November 11, 2015 4:52 PM
• Ok now I am running.

PS It does look cool Reed.

Wednesday, November 11, 2015 5:37 PM
• Ok now I am running.

PS It does look cool Reed.

Now you are in for it!

In divideorzero:

```    Public Shared Function DivideOrZero(numerator As Single, denominator As Single) As Single
'If Equal(denominator, 0F) Then Return 0F
If denominator = 0F Then Return 0
Return numerator / denominator
End Function```

the function returns 0 for divide by zero. However isnt the value infinity? I know you cant return that (or maybe you can?) but it seems it should return a large value of some kind. Like say you are computing the slope of a vertical line by dy/dx and dx = 0 then the function returns zero for the slope of a vertical line, but that is a horz line.

What I have always done is done not knowing the "proper" how to do it is check

if x = 0 then x = .00001

Depends on what one is doing I am sure but so far its all I can say other than looks cool. Also I use this function heavily so just curious. In a CAD type thing it can make for interesting results I think.

I wonder what the graphic is? Just a test? Or is it a way to build stone walls?

I see the err you had as once the "wall" should have settled but something was buzzing around below and then eventually a rock poped up into the air and it was really spinning.

:)

PS I guess you are not even using this for divide by zero like I was thinking. If so nevermind.
Wednesday, November 11, 2015 8:03 PM

• In divideorzero:

```    Public Shared Function DivideOrZero(numerator As Single, denominator As Single) As Single
'If Equal(denominator, 0F) Then Return 0F
If denominator = 0F Then Return 0
Return numerator / denominator
End Function```

the function returns 0 for divide by zero. However isnt the value infinity? I know you cant return that (or maybe you can?) but it seems it should return a large value of some kind. Like say you are computing the slope of a vertical line by dy/dx and dx = 0 then the function returns zero for the slope of a vertical line, but that is a horz line.

What I have always done is done not knowing the "proper" how to do it is check

if x = 0 then x = .00001

Tommy, I think its a convention here, because he checks for im being zero in the methods and either continues the loop or jumps out of the method-body.

Regards,

Thorsten

Wednesday, November 11, 2015 9:32 PM
• Looks like you guys figured out the DivideOrZero method; yes, zero is the intended result and is supposed to be taken into consideration when using the method.

Note that there is also an epsilon value for the "approximately zero" check that you mentioned.  Both are useful techniques.  :)

Reed Kimble - "When you do things right, people won't be sure you've done anything at all"

Wednesday, November 11, 2015 11:25 PM
• I haven't even begun to experiment with the physics variables yet.  I've been killing myself trying to squash that bug!

I plan to refactor a lot of this and one of the refactorings is to implement the Material type.  It will take some experimenting to find suitable values for different materials like metal, wood, concrete, rubber and so on.

Reed Kimble - "When you do things right, people won't be sure you've done anything at all"

Wednesday, November 11, 2015 11:29 PM