In the previous post, I showed how enforcing $\mathbf{Jv}=0$ imposes one or more constraints. If you have gone ahead and implemented this, you will notice that the system often collapses (especially when gravity is involved) - once the constraint is slightly violated, it never recovers, and the associated rigidbodies plummet down. However, consider what would happen if we were to add in a “bias” term (should be familiar to anyone who’s tried their hand at machine learning): $\mathbf{Jv+b}=0$. Roughly, this allows for less leeway in the positional drift, since we can feed the position error (which we had initially discounted, opting straight for velocity analysis) back into the velocity constraint. Additionally, a time-dependent bias allows us to incorporate objects with prescribed motion, like motors.
It is clear that $\boldsymbol\lambda$ will now be solved as
$$
\boldsymbol\lambda=(-\mathbf{Jv\color{red}{+b}})\left(\mathbf{JM}^{-1}\mathbf J^\top\right)^{-1}
$$
Let’s focus on the position error. This is represented by the violation of the original constraint $C({q_i})\ne0$, which we would like to indirectly enforce through the velocity. However, to account for the positional drift, we can set the bias to be $\mathbf b=\frac{C}{\Delta t}$, where $\Delta t$ is added into the denominator to make the expression dimensionally consistent (this is extremely important in physics engines, since computationally, missing factors of $\Delta t$ will cause $\mathcal O(10^{2})$ variations in values!). Why is this?
Look at the ordinary differential equation in terms of the constraint:
$$
\mathbf{Jv}+\frac C{\Delta t}=0
\\\frac{\mathrm d C}{\mathrm d t}+\frac C{\Delta t}=0
\\C=C_0\exp(-\frac{t}{\Delta t})
$$
In other words, if the constraint is initially violated by $C_0$, then the bias term will assist in damping this value to zero (the desired value), exponentially, rather than letting positional errors accumulate. This is known as Baumgarte stabilisation. The bias term is typically equipped with a control parameter $\beta\in[0,1]$ to yield $\mathbf b = \frac{\beta C}{\Delta t}$, and this determined through trial and error, and may even be dependent on the joint type.
An implementation along these lines will be extremely sturdy in general. But the devil lies in contact constraints, which enforce collisoin resolution. Let’s take a look.
Contact constraints are an example of inequality constraints, since clearly we want to objects to have a penetration depth of anything less than zero. In practice, determining the penetration depth for two generic convex colliders is very complicated, but we can blackbox that for now - or alternatively work with a simplified system in which we only have circles (to avoid cheating and not implementing rotation, imagine each circle has a line drawn from the radius to a point on the circumference).
The applied velocity must be such that the bodies are just separated, with zero relative velocity along the normal, upon separation (restitution, i.e. elastic collisions, will be implemented later with minimal changes). This can be implemented through the position constraint
$$
C=(\vec p_1 -\vec p_2)\cdot\vec n=0
$$
where $\vec p_1$ and $\vec p_2$ are the deepest penetrating points of the two rigibodies, and $\vec n$ is the normal vector separating the two bodies. In other words, a translation of either body through a displacement of $\lVert\vec p_1 -\vec p_2\rVert\vec n$ would just about separate the two - but of course we wish to implement this at the velocity level. In doing so, we can eyeball the Jacobian (this is the general procedure for any constraint).
$$
\frac{\mathrm dC}{\mathrm dt}=(\vec v_1|_{\vec p_1} -\vec v_2|_{\vec p_2})\cdot\vec n=0
$$
Here the product rule was used, but the value of $\frac{\mathrm d\vec n}{\mathrm dt}$ is $0$ since the translation will occur along the normal and hence it will remain constant. Importantly, this expression does not involve the motion of the centre of mass $\vec x_1$: it is instead the velocity of another point $\vec p_1$, typically on the boundary of that collider! This means we also have to account for rotational motion, and hence $\vec v_1\vert_{\vec p_1}=\vec v_1\vert_{\vec x_1}+\vec\omega_1\times(\vec p_1-\vec x_1)$, by definition. Thus the velocity constraint becomes
$$
\big(\vec v_1|_{\vec x_1}+\vec\omega_1\times(\vec p_1-\vec x_1)-\vec v_1|_{\vec x_2}+\vec\omega_2\times(\vec p_2-\vec x_2)\big)\cdot\vec n=0
$$
Finally, using the vector identity $(A\times B)\cdot C=(B\times C)\cdot A$, the previous equation can be converted to the form
$$
\big[\vec n\quad (\vec p_1-\vec x_1)\times\vec n\quad-\!\vec n\quad-\!(\vec p_2-\vec x_2)\times\vec n\big]\begin{bmatrix}\vec v_1\\\vec\omega_1\\\vec v_1\\\vec\omega_1\end{bmatrix}=0
\\\mathbf{Jv}=0
$$
The Jacobian becomes manifest, without having to calculate a vector of partial derivatives! Finally, we can introduce our newfound Baumgarte stabilisation to regulate the position error. Adding restitution is straightforward from here: the velocity along the normal should not be zero after the update (inelastic collision, as we have done here), but rather a fraction of the incoming relative velocity during collision. In other words
$$
\frac{\mathrm dC}{\mathrm dt}=\mathbf{Jv+b}\ge-\min(e_1,e_2)\ \bigg[(\vec v_1|_{\vec p_1} -\vec v_2|_{\vec p_2})\cdot\vec n\bigg]_\text{pre-impulse}
\\\frac{\mathrm dC}{\mathrm dt}=\mathbf{Jv+b}'\ge0
\\\mathbf b'=\frac{\beta C}{\Delta t}+\min(e_1,e_2)\ \bigg[(\vec v_1|_{\vec p_1} -\vec v_2|_{\vec p_2})\cdot\vec n\bigg]_\text{pre-impulse}
$$
The $\min(e_1, e_2)$ refers to the lower value of the restitutions of the two colliders. This is conventional, and other schemes like the average between the two can also be employed (but clearly $e\in[0,1]$).
Coupled with the solvers described in the previous post, we now have a reasonably stable, fast and accurate engine that can handle all sorts of constraints, including collisions. However, the collisions can be a bit shaky and unreliable, so I will describe some creative methods to deal with this in a future post. Stay tuned!
Flutter, developed by Google, is a (relatively) recent entrant to the mobile app development scene, where it competes with the likes of React Native to deliver a framework for frontend UI development and smooth backend integration.
Since I haven’t personally used React Native all that much, I will primarily focus on the elements of Flutter that I enjoyed, rather than stage a boxing match between the two. A bit of a backstory first: I originally envisioned myself as a “fundamentals” programmer, and I didn’t really see myself enjoying app development, since it seemed to be primarily design rather than algorithms, and this ambivalence was stoked by Android Studio’s sprawling clunkiness. However, I can honestly say that Flutter has converted me - I now enjoy front-end development, and am exploring back-end technologies and proto service architecture too!
The primary point of appeal, especially for someone who is aiming to go the whole hog and widely publish their app, is that Flutter is natively cross-platform - write once, run anywhere in its modern incarnation. Saves you from having to splurge on a Mac to get Xcode, and then rewrite your Android code for iOS. Oh and - this doesn’t mean the end result will be exactly the same for the two, which might lead to unnatural formatting and components. Upon writing the code once, the UI elements take the form of the native components on their respective system! As of late, it’s no longer limited to mobile devices, Flutter has burst onto the web app and IoT scene too!
Flutter also integrates seamlessly with Firebase (Google’s service platform for mobile and web applications), so there are solid backend options already in place. The community-developed SQLite package is also very sturdy and apt for local storage: together, they form a very effective duo.
The UI system is a pleasure to work with. Flutter’s motto is “Everything is a Widget”, so immediately after conceiving the layout of your app, you can jump straight into an extremely modular but compact design. It does use its own language, known as Dart, but this is hardly a downside since the language is extremely similar to Java/C#, and nobody who’s been programming for more than a couple of years will have any problem picking it up simultaneously. That said, I feel like they could have employed a slightly more modern, quirky language - Dart is a bit staid and unadventurous, but it gets the job done.
Features from http requests to JSON serialisation are a breeze, and traditional UI elements like animation, navigation and themes are extremely intuitive to use and look fantastic, build times are reasonable, and the “hot reload” feature is very handy. Once you have a handle on the numerous widgets and the core principles of their sizing, one can develop impeccable-looking apps very rapidly - I was pleasantly surprised, since the only UI I had done prior to this was Java’s Swing API (I’m old school, I learned from books written in the late ’90s). The error messages are not the best, but most problems are generic and solvable via the Stack Overflow corpus.
The Flutter community also seems to be very integrated, helpful, and proactive in package development, and the package management system is as efficient as cargo in Rust. However Flutter is still maturing, and a uncomfortably large reliance on external, only moderately stable packages is often inevitable for large apps. But the size of the community support and effort on these repositories is just indicative of how effectively Flutter is bringing in and retaining app developers - despite its fledgling status, it boasts an impressive array of widgets, solid documentation and toolchain integration.
I have personally used Flutter to develop the frontends of a wide array of applications, including a news aggregator and an app to connect students with propsective mentors - I’m loving the experience, and I wholeheartedly recommend that any app developer, whether a beginner or an Android Studio veteran, give Flutter a go.
For reasons of stability and accuracy, we typically work with velocity-based constraints to introduce interesting features into our physics engine, like joints and collisions. This approach outstrips its peers, since:
- Acceleration-based updates lead to extreme instabilities as the accumulated acceleration can go to infinity
- Position-based updates either look jittery or unrealistic (“hugging”), without a good balance in between. You can’t even implement friction at a velocity level.
These constraints are well-defined, and are essentially 100% physically accurate. That’s why the same story carries forth to rendered physics, realistic simulations and robotics. The key difference is that game physics engines are time-critical, with ever-increasing demands for “higher FPS” by gamers leading to ever-decreasing amounts of self-contentment in low-level game programmers. In other words, it’s not the constraints themselves that are the problem, but rather solving them - the time complexity of obtaining the exact solution is cubic in the number of constraints, which is a nightmare for any physics-based game after the 1990s. As the programmer’s saying goes: “It’s $\mathcal O(n)$ or nothing”.
So what is a constraint, mathematically? We have at our disposal the kinematic variables ${q_i}$ of all the rigid bodies in the scene. In practice, $q_i \equiv{\vec x_i, \vec v_i, \vec\theta_i, \vec\omega_i}$ (okay, the vectors on the angular variables are gratuitous for 2D. Whatever, it looks cooler). Then a constraint is simply $C({q_i})\overset!=0$ (equality constraint) or $C({q_i})\overset!\ge0$ (inequality constraint). For instance, a ball-and-socket or fixed-distance joint between $\mathbf 1$ and $\mathbf 2$ is given by $C(\vec x_1, \vec x_2)=\lVert\vec x_1-\vec x_2\rVert=0$. Nothing to it! Unfortunately, apropos the disclaimer I issued above, rigidly enforcing $C=0$ is a mistake, so we work with $\partial_t C=0$ instead. In words, we want the constraint to remaining unchanging over time, and correct the velocity of the rigidbodies in question in order to enforce this. Through the kindergarten multivariable chain rule, we have
$$
0=\frac{\partial}{\partial t}C(\{q_i\})=\sum_i\frac{\partial q_i}{\partial t}\frac{\partial}{\partial q_i} C(\{q_i\})\tag{1}
$$
Specialising to constraints between two rigidbodies (i.e. practically all of them), we can adapt this to a more linear algebraic form by defining the the stacked velocity vector $\mathbf v = [\vec v_1\ \vec\omega_1\ \vec v_2\ \vec\omega_2 ]$ - this is a column vector with length 12 (or 6 in 2D). Similarly, $\mathbf J\equiv\frac{\partial C({q_i})}{\partial q_i}$, which is a row vector or the same length, known as the Jacobian (a familiar term from multivariable calc). If you are aware of Jacobians acting roughly as coordinate transformations, then you can view it here as a projection map down from ordinary Euclidean space to constraint space. Then equation $(1)$ becomes $\mathbf{Jv}=0$. We can interpret this as meaning that the Jacobian vector is perpendicular to the velocity. But this is precisely the same property obeyed by the constraint force!
Why? Recall the principle of virtual work. Since constraint forces are nudging the system back into static equilibrium, they crucially cannot do any work. Consequently for the constraint forces, $\mathbf F\cdot\mathbf v=0$, which implies that $\mathbf F=\mathbf J^\top\lambda$, where the transpose is merely due to the distinction between row and column vectors, and $\lambda$ acts as a Lagrange multiplier - it is this value we need to solve to obtain the requisite update, since the value of the Jacobian can be easily computed using the current velocities and the type of constraint.
Now at any given time step, the constraint will typically be violated due to the presence of external forces acting on the body, say due to user input, powered motion, etc. A correctional velocity must be applied to force $\partial_tC$ back to $0$. The velocity update (this is shorthand for the updates to the velocity and angular velocity of both bodies) is by definition $\delta\mathbf v=\mathbf M^{-1}\mathbf F=\mathbf M^{-1}\mathbf J^\top\lambda$. Here $\mathbf M=\mathrm{diag}[m_1\mathbb{I}^3, \mathbf{I}_1, m_2\mathbb{I}^3, \mathbf{I}_2]$, where $m$ and $\mathbf I$ are the mass and inertia tensor respectively - this definition should be intuitively obvious from the definition of the generalised velocity. Thus $\lambda$ is determined by
$$
\mathbf {Jv}\ne0,\quad \mathbf J+(\mathbf v+\delta\mathbf v)=0
\\\Rightarrow\mathbf J(\mathbf v+\mathbf M^{-1}\mathbf J^\top\lambda)=0
\\\mathbf{Jv}=-\mathbf{JM}^{-1}\mathbf J^\top\lambda
\\\Rightarrow\lambda =-\mathbf{Jv}\left(\mathbf{JM}^{-1}\mathbf J^\top\right)^{-1}
$$
In other words, at each time step, for each constraint, we update $\mathbf v$ of the involved bodies as
$$
\mathbf v\rightarrow\mathbf v-\mathbf M^{-1}\mathbf J^\top\mathbf{Jv}\left(\mathbf{JM}^{-1}\mathbf J^\top\right)^{-1}
$$
The problem here of course is that updating one constraint at a time generically breaks all of the other ones, so naïvely it appears we cannot enforce all the constraints: are we doomed?
Of course not. The matrix notation above is very suggestive. Rather than having a length-12 row vector as a Jacobian, we can simultaneously incorporate multiple constraints by horizontally stacking the two individual Jacobians to form a $2\times12$ matrix. $\lambda$ is also promoted to a matrix, but $\delta\mathbf v$ remains a column vector, as you can check. Thus, $\mathbf{Jv}=-\mathbf{JM}^{-1}\mathbf J^\top\boldsymbol\lambda\sim\mathbf{b}=\mathbf A\boldsymbol\lambda$ is now a matrix equation in $\boldsymbol\lambda$. As with any other matrix equation, it is best to solve it using something like LU decomposition rather than inverting $\mathbf A$, which would be extremely inefficient.
This will essentially yield perfect results, provided a global solution does indeed exist for all the constraints. However, this method is downright terrible for game physics engines! This is because the size of the stacked Jacobian scales as $\mathcal O(N^2)$, where $N$ is the number of bodies, and so you will have to solve extremely large matrix equations to determine $\boldsymbol\lambda$, with a time complexity of $\mathcal O(N^{\sim2.8})$. Also, from experience, it can lead to single-frame collapse, where there is no global solution, so it just ends up destroying all the constraints in the absence of time-consuming mitigation tactics. That said, this method is perfect for pre-built rendered physics, where speed is not a major concern.
However this is just one method. The second, espoused by Box2D creator Erin Catto, is sequential impulse, where we solve the constraints sequentially and iteratively, in the hope that after several runs, the accumulated velocity correction converges to the global solution. The theoretical time complexity is $\mathcal O(N)$ space and time, though as with most Big-O notations, this hides a constant prefactor. The linear time is because there are no matrix equations to be solved, just fast, constant-size matrix multiplications. This works fantastically … provided you iterate enough times - in a naïve implementation, this may still be too slow. Nevertheless, plenty of work has been channeled into optimisation and stability schemes in this regime, which makes it the most popular framework today.
However, there are certain cases where you actually do want to solve multiple constraints simultaneously for stability reasons - the best example I can think of is a resting inequality constraint: there are two contact normals generated with the following constraint: $C=(\vec x_2-\vec x_1)\cdot\vec n=0$, but solving them sequentially will lead to the top rigidbody oscillating back and forth rather than remaining still. To counter this, the block-solving method identifies “blocks” of objects that you want to solve together as a matrix equation. Then, each block is solved iteratively as in sequential impulse above, bringing together the advantages off both previous methods.
All of these methods, in their raw forms, still raise questions of stability. In the next post, I will satiate the enthusiastic reader with concrete examples of the concepts and methods discussed above, and also describe inequality constraints and various methods to improve stability.
I’ve found myself exploring some really intriguing and novel (at least to me) topics recently. In particular, an underlying theme of many of these is that they are topics at the crossroads of advanced field theory and elements from other areas of physics:
- The instanton fluid model
- The thermodynamics and bulk propagation in the chiral condensate
- The XY model and the Kosterlitz-Thouless transition - it’s a statistical model which exhibits duality with magnetic fields. It also provides a beautiful example of topological defects, spin waves and a non-trivial phase transition, with vortex-antivortex unbinding.
- The Skyrme model of baryons. I’ve always loved chiral perturbation theory, i.e. the pion model, but I only recently learned about how baryons can be described as solitons, or topological twists of the pion field. Simply fantastic. I explored non-linear sigma models in more detail, symmetry breaking, ABJ and ‘t Hooft anomalies, and skyrmions even link with many-body nuclear systems.
- Lattice gauge theory. One day I’d love to implement a reduced-dimension toy model. It provides insight on glueballs, flux tubes, the deconfinement phase of QCD and what goes on in condensed matter theory
- The Ryu-Takayanagi formula and Page curves for black hole evaporation
- Asymptotically conserved surface charges, some more AdS3/CFT2, and extremal BTZ black holes
- Type II String compactifications - I am attending a few lectures in the Strings conference after all!
- Galois theory, just for fun!
Some stuff which I would like to get started with is:
- The SYK model, since it provides a solvable example of AdS/CFT
- JT gravity and matrix models
- Group cohomology in ‘t Hooft anomalies
- The D1/D5 system and the near-horizon geometry of extremal black holes
- K3 surfaces
- Thermal QFT and symmetry restoration