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 Β 

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βˆ—np*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).

pn=(xxxyxzvxvyvz)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:

x¨(t)=f(x(t),x˙(t)) \bf \ddot x \textnormal{(t)} = f(x\textnormal{(t)}, \dot x\textnormal{(t)})

where xΒ¨(t)\ddot x(t) is the second-order derivative of position x(t). We already know that this value is the acceleration. xΛ™(t)\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=xΛ™v = \dot x, we can rewrite these equations to first order:

ddt(x(t)v(t))=(v(t)f(x(t),v(t))) \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:

xnew=x0+hf(xnew)βˆ΅Ξ”x=xnewβˆ’x0βˆ΄Ξ”x=hf(xnew) x_{new} = x_0 + hf(x_{new}) \\ \because \Delta x = x_{new} - x_0 \\ \therefore \Delta x = h f(x_{new})

and that:

Ξ”x=x(t0+h)βˆ’x(t0)=ddtx(t) \Delta x = x(t_0 + h) - x(t_0) = \frac{d}{dt} x(t)

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

(Ξ”xΞ”v)=h(v0+Ξ”vf(x0+Ξ”x,v0+Ξ”v)) \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}


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

f(xn+Ξ”x,vn+Ξ”v)β‰ˆfn+βˆ‚fβˆ‚xΞ”x+βˆ‚fβˆ‚vΞ”v 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)f(x, y) using the Taylor series as follows:

f(x+Ξ”x,y+Ξ”y)β‰ˆf(x,y)+βˆ‚f(x,y)βˆ‚xΞ”x+βˆ‚f(x,y)βˆ‚yΞ”y+... 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:

(Ξ”xΞ”v)=h(v0+Ξ”vfn+βˆ‚fβˆ‚xΞ”x+βˆ‚fβˆ‚vΞ”v) \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 Ξ”x\Delta x and Ξ”v\Delta v are decided by us. In this example, we can assume that our particle is defined by a 3D position x∈R3\bf x \in \R^3 and a velocity v∈R3\bf v \in \R^3.

This distinction helps us better understand the dimensions of this calculation. Let’s assume in this case that ff is 3D. A resource I recommend for creating cloth simulations tells us that ff 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 βˆ‚fβˆ‚x∈R3x3\frac{\partial f}{\partial x} \in \R^{3x3}

πŸ’‘Thanks to vector calculus, we see how to calculate the derivation of a vector in 3D:

[fx(x,y,z)fy(x,y,z)fz(x,y,z)]=[βˆ‚fxβˆ‚xβˆ‚fxβˆ‚yβˆ‚fxβˆ‚zβˆ‚fyβˆ‚xβˆ‚fyβˆ‚yβˆ‚fyβˆ‚zβˆ‚fzβˆ‚xβˆ‚fzβˆ‚yβˆ‚fzβˆ‚z] \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 Ξ”x\Delta x in the equation for Ξ”v\Delta v, we arrive at the following:

Ξ”v=h(f0+βˆ‚fβˆ‚xh(v0+Ξ”v)+βˆ‚fβˆ‚vΞ”v) \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 Ξ”v\Delta v goes on either side of the equals sign. It is not enough to resolve the value of the variable Ξ”v\Delta v by placing it on one side only. So we have to solve a linear system using the linear form Ax=bAx = b.

πŸ’‘ A linear system can be solved using the form:

Ax=b 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=2xβˆ’y=45y+z=βˆ’1 3x + 2 = 2 \\ x - y = 4 \\ 5y + z = -1

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

[3201βˆ’10051][xyz]=[24βˆ’1] \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=bAx = b, we arrive at:

(Iβˆ’hβˆ‚fβˆ‚vβˆ’h2βˆ‚fβˆ‚x)Ξ”v=h(f0+hβˆ‚fβˆ‚xv0)\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 (I∈R3x3\bf I \in \R^{3x3}).

The next step is of course to solve Ξ”v\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!


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

Written by Carmen Cincotti, computer graphics enthusiast, language learner, and improv actor currently living in San Francisco, CA. Β Follow @CarmenCincotti


Interested in contributing to Carmen's Graphics Blog? Click here for details!