Numerical Methods (Implicit, 2nd order ODE) - Part 4

We are learning how and when to use implicit numerical methods to solve second order differential equations.

Keywords: C, numerical methods, math, physical simulations

By Carmen Cincotti

April 11, 2022

Let’s review: Phase Space
Second Order Ordinary Differential Equations (ODE)
Linearization
- An aside: dimensions
Linear form
Next time
Resources

I am continuing my work through the Pixar Physics Simulation Animation Course, in particular the article on implicit numerical methods. Last week, we introduced first-order implicit methods. This week, we’ll be studying second-order implicit methods.

This lesson will open the door to understanding and creating cloth simulations. Of course, cloth simulations (woven and knitted) are all about using mass-spring systems. Unlike the simulation of the particle system we created a few weeks ago where we used an explicit solver, we’ll be using a second-order implicit Euler solver for our cloth simulation!

Let’s review: Phase Space

As I said, a cloth simulation is like a particle system simulation. Recall that the particles are represented in space in a configuration called phase space.

Each particle is represented in space by two generalized coordinate vectors (x, v) where x is a position vector and v is a velocity vector. The phase space vector $p*n$ builds itself from vectors x and v one after the other, thus it will have a size of 2*n elements (where _n* is the amount of dimensions represented in our system).

p_n = \begin{pmatrix} x_x \\ x_y \\ x_z \\ v_x \\ v_y \\ v_z \\ \end{pmatrix}

⚠️ If you are already lost, see the discussion on particles where I introduced this topic. It’s a good idea that you understand this topic before proceeding!

Second Order Ordinary Differential Equations (ODE)

The Pixar resource states:

Most dynamic problems are written in terms of a second-order differential equation:
$\bf \ddot x \textnormal{(t)} = f(x\textnormal{(t)}, \dot x\textnormal{(t)})$

where $\ddot x(t)$ is the second-order derivative of position x(t). We already know that this value is the acceleration. $\dot x(t)$ is the derivative of position… known as velocity.

Solving a second-order differential equation is not ideal because it is expensive and complicated. We’ll see that we can avoid this calculation by approximating the second-order term using the Taylor series in a process called linearization. By introducing and setting the variable $v = \dot x$ , we can rewrite these equations to first order:

\frac{d}{dt}\begin{pmatrix} \bf x\textnormal{(t)} \\ \bf v\textnormal{(t)} \end{pmatrix} = \begin{pmatrix} \bf v\textnormal{(t)} \\ \bf f(x\textnormal{(t)}, \bf v\textnormal{(t)}) \end{pmatrix}

Let us simplify this system by recalling that the implicit Euler step is defined as follows:

x_{new} = x_0 + hf(x_{new}) \\ \because \Delta x = x_{new} - x_0 \\ \therefore \Delta x = h f(x_{new})

and that:

\Delta x = x(t_0 + h) - x(t_0) = \frac{d}{dt} x(t)

…we can introduce variables $x_0 = x(t_0)$ and $v_0 = v(t_0)$ to finally have an equation to solve an implicit Euler step:

\begin{pmatrix} \bf \Delta x \\ \bf \Delta v \end{pmatrix} = h \begin{pmatrix} \bf v_0 + \Delta v \\ \bf f(x_0 + \Delta x, v_0 + \Delta v) \end{pmatrix}

Linearization

As I said, we can approximate $f(x_0 + \Delta x, v_0 + \Delta v)$ thanks to Taylor series expansion:

f(x_n + \Delta x, v_n + \Delta v) \approx f_n + \frac{\partial f}{\partial x} \Delta x + \frac{\partial f}{\partial v} \Delta v

💡 We can approximate a two-variable real function $f(x, y)$ using the Taylor series as follows:

f(x + \Delta x, y + \Delta y) \approx f(x,y) + \frac{\partial f(x,y)}{\partial x} \Delta x + \frac{\partial f(x,y)}{\partial y} \Delta y + ...

See this resource for more on the Taylor Series and this exact derivation.

So the equation becomes:

\begin{pmatrix} \bf \Delta x \\ \bf \Delta v \end{pmatrix} = h \begin{pmatrix} \bf v_0 + \Delta v \\ \bf f_n + \frac{\partial f}{\partial x} \Delta x + \frac{\partial f}{\partial v} \Delta v \end{pmatrix}

An aside: dimensions

It must be said that the dimensions of $\Delta x$ and $\Delta v$ are decided by us. In this example, we can assume that our particle is defined by a 3D position $\bf x \in \R^3$ and a velocity $\bf v \in \R^3$ .

This distinction helps us better understand the dimensions of this calculation. Let’s assume in this case that $f$ is 3D. A resource I recommend for creating cloth simulations tells us that $f$ is the force in the system, but I’m not so sure at this point because I haven’t studied it yet. In any case, we can assume that $\frac{\partial f}{\partial x} \in \R^{3x3}$

💡Thanks to vector calculus, we see how to calculate the derivation of a vector in 3D:

\begin{bmatrix} f_x(x,y,z) \\ f_y(x,y,z) \\ f_z(x,y,z) \end{bmatrix} = \begin{bmatrix} \dfrac{\partial f_x}{\partial x} & \dfrac{\partial f_x}{\partial y} & \dfrac{\partial f_x}{\partial z} \\[1em] \dfrac{\partial f_y}{\partial x} & \dfrac{\partial f_y}{\partial y} & \dfrac{\partial f_y}{\partial z} \\[1em] \dfrac{\partial f_z}{\partial x} & \dfrac{\partial f_z}{\partial y} & \dfrac{\partial f_z}{\partial z} \\[1em] \end{bmatrix}

The derivation is called the Jacobian matrix. We will come back to this with more details in a later part

Substituting $\Delta x$ in the equation for $\Delta v$ , we arrive at the following:

\bf \Delta v = h (\bf f_0 + \frac{\partial f}{\partial x} \textnormal{h}(v_0 + \Delta v) + \frac{\partial f}{\partial v} \Delta v)

Linear form

We have finally arrived at the heart of an implicit method. We can easily see that $\Delta v$ goes on either side of the equals sign. It is not enough to resolve the value of the variable $\Delta v$ by placing it on one side only. So we have to solve a linear system using the linear form $Ax = b$ .

💡 A linear system can be solved using the form:

Ax = b

where A is the matrix of coefficients, x is the solution vector and b is the constant vector. Consider this example with a linear system as follows:

3x + 2 = 2 \\ x - y = 4 \\ 5y + z = -1

Using the form $Ax = b$ , we arrive at:

\begin{bmatrix} 3 & 2 & 0 \\ 1 & -1 & 0 \\ 0 & 5 & 1 \\ \end{bmatrix} \begin{bmatrix} x \\ y \\ z \\ \end{bmatrix} = \begin{bmatrix} 2 \\ 4 \\ -1 \\ \end{bmatrix}

After arranging our terms in the form $Ax = b$ , we arrive at:

\bf(I - \textnormal{h}\frac{\partial f}{\partial v} - \textnormal{h}^2\frac{\partial f}{\partial x})\Delta v = \textnormal{h}(f_0 + \textnormal{h}\frac{\partial f}{\partial x}v_0)

where Iis the identity matrix ( $\bf I \in \R^{3x3}$ ).

The next step is of course to solve $\Delta v$ . I will do that in the next series - when we create a cloth simulation!

Next time

I continue to study physics by creating a cloth simulation!

Carmen's Graphics Blog