# 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

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!

## 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