This involves solving ODEs (ordinary differential equations) that are considered to be ‘stiff’. We will see why we qualify such ODEs as ‘stiff’ and how we should approach solving for these types of ODEs by using implicit numerical methods such as the Backwards Euler implicit method.
⚠️ Pixar’s resource advises us to avoid an stiff ODE by reformulating our dynamic problem. We’ll be focusing on the cases in which this is not possible.
After having scoured the Internet for a definition of the term ‘stiff ordinary differential equation’, I ended up on Wikipedia’s entry which, admittedly for me, is not very sufficient but I will give it to you anyway:
A stiff differential equation is a differential equation whose sensitivity to parameters will make it difficult to solve by explicit numerical methods.
Definition aside, I find that Pixar Physics Simulation Animation Course gives us a fun example that we will work on to better visualize the Wikipedia definition… Let’s imagine that we want to simulate the dynamics of a particle whose change of position is described by the following differential equation:
where k is a positive and large constant. If k is large enough, the particle will never stray too far from y(t)=0 since −ky(t) will constrain y(t) back to the origin. We also see that the x position of the particle is freer to move around, as described by the absence of this large and positive constant k.
In this example, this ODE would therefore be ‘stiff’ if the variable k is so large that the particle stays close to the line y=0. As we will see in the next part, explicit methods will have a hard time simulating this expected behavior.
The problem with explicit methods
Recall that the equation of the explicit Euler method is described as follows:
and in this example, X0=(x0y0) and X˙(t)=dtd(x(t0)y(t0))=(−x0−ky0). Then, after the substitution, we have:
Which produces the following:
We learned in Part 1 of Numerical Methods that in order to find a solution, we must take small steps in time to approach the solution as massive steps could very well make the system unstable. A quick rearrangement of terms will show us that these conditions must be met in order for the position of the particle described in our example to converge to (0,0):
🚨 If we take k equal to 1000: the maximum permitted step is calculated to be 1−1000h=−1=>h=0.002. If we surpass this limit, the particle’s y position will never be able to converge towards 0!
Now, let’s take larger steps in time, where h > 0.002 (h = 0.003 and k = 1000.0). Notice how the solution diverges from 0:
OK, so we could use a very small step like h = 0.0015, right? Maybe… however, we would quickly realize that the x position converges to 0 very slowly, which will make our simulation slow, and expensive!
Can we do better? Could we use a more stable method to converge faster? Yes - with the Backwards Euler implicit method!
Backwards Euler (implicit)
Backwards Euler is the easiest of the implicit methods to understand and that’s why I’m going to focus on it.
Here is the equation:
If you had been an observant person, you would have seen that the difference between the implicit Euler equation and the explicit Euler equation is that the change in X is defined as hX˙∗new in the implicit Euler equation. We must therefore solve this change by calculating X˙∗new… to finally calculate X_new.
In fact, it’s like walking backwards.
In other words, it’s as if we were at position X_new and wanted to move backwards to get to position X0… Hence the name: ‘Backwards Euler’.
It’s a little awkward to solve X∗new, since we need to find X˙∗new… which is also unknown. Let’s see if we can rewrite this into something we can solve.
We start by providing the equation of the Backward Euler method:
where ΔX=X_new−X0. We can write it as follows:
💡 f(X∗0+ΔX) is unknown, and also ΔX (because ΔX=X∗new−X∗0 contains X∗new which is unknown). So we have to approximate f(X0+ΔX) using the Taylor series as follows:
See this resource to learn more about the Taylor Series and this exact derivation.
Using the last equation, we can rewrite our equation for ΔX
💡 To continue, remember that f(X0) is a vector and f′(X0) is a matrix. This is the case thanks to vector calculus: