Dzhelil Rufat Джелил Руфат en-us Sun, 04 Dec 2016 00:00:00 -0800 <![CDATA[Propagation of a Wave over a Non-Flat 2D Space Emedded in 3D]]> Propagation of a Wave over a Non-Flat 2D Space Emedded in 3D

We are simulating the wave equation applied to a zero form subject to Neumann boundary condtions

\[\ddot{\phi} = \star \mathbf{d} \star \mathbf{d} \phi\]

We use an explicit time integrator.

First, we observe how the wave propagates on a flat surface.

Second, we observe how the wave propagates on a surface with a small ditch in the middle.

Observe that as the wave passes through the hole it has to transverse an extra distance and thus the wave front is delayed relative to the front outside of the hole.

Let us also observe what happens when the hole gets deeper.

Sun, 04 Dec 2016 00:00:00 -0800 <![CDATA[Orthonormal Bases]]> Orthonormal Bases

Given a standard orthornormal reference frame, consisting of the unit vectors along the x and y directions, the video below shows a naive mapping of that frame to the surface. When the surface is distorted, one can see that the frame is no longer orthonormal with respect to the metric of the space in which the surface is embedded.

In the video below we compute an orthonormal frame for each point using the square root method.

In the video below we domonstrate the texture unrolled onto a flat surface.

Mon, 28 Nov 2016 00:00:00 -0800 <![CDATA[Different Boundary Conditions under Spectral DEC]]> Different Boundary Conditions under Spectral DEC

To validate the code and ensure that it does what is expected of it, we start off with a simple example of a linear wave travelling in one direction.

The checkerboard that the surface is shaded with corresponds to the actual cellular complex.




Having validated the code with a simple example, now we can consider a more complicated wave form. In the examples, below we will start off with a Gaussian bump in the middle of the domain and see how it spreads out under different boundary conditions.




Thu, 27 Oct 2016 00:00:00 -0700 <![CDATA[Locally Orthonormal Frame]]> Locally Orthonormal Frame

At each point we will also store a matrix \(J\) such that \(J^{T}J=g\) that will map to a locally orthonormal basis. The hodge star operator will be implemented in that basis as

\[\begin{split}\left(\begin{array}{cc} 0 & -1\\ 1 & \phantom{+}0 \end{array}\right)\end{split}\]

Coordinate bases for vectors and forms

\[\frac{\partial}{\partial x^{i}}\qquad dx^{i}\]

Non-coordinate orthonormal bases for vectors and forms


Transformations to orthornormal basis

\[\hat{e}_{i}=J_{i}^{\phantom{j}j}\frac{\partial}{\partial x^{j}}\quad\hat{\theta}^{i}=J_{\phantom{i}j}^{i}\,dx^{j}\]

where \(J_{i}^{\phantom{j}j}\) and \(J_{\phantom{i}j}^{i}\) are transposed inverses of each other


The above property is required because the pairing of one forms with vector fields remains invariant

\[\hat{e}_{i}\hat{\theta}^{i}=\frac{\partial}{\partial x^{i}}dx^{i}\]

Relationship to the metric


Transforming one-forms


Transforming two-forms

\[\beta=\beta_{ij}\,dx^{i}\wedge dx^{j}=\hat{\beta}_{ij}\,\hat{\theta}^{i}\wedge\hat{\theta}^{j}\]
Fri, 01 Jul 2016 00:00:00 -0700 <![CDATA[Riemannian Metric and the Hodge-Star in 3D]]> Riemannian Metric and the Hodge-Star in 3D

Euclidean Space

The flat metric is

\[\begin{split}g=\left[\begin{array}{ccc} 1 & 0 & 0\\ 0 & 1 & 0\\ 0 & 0 & 1 \end{array}\right]\end{split}\]

Given one-forms \(\mathbf{U}=U_{x}dx+U_{y}dy+U_{z}dz\) and \(\mathbf{V}=V_{x}dx+V_{y}dy+V_{z}dz\), for brevity we will write them as column vectors

\[\begin{split}\mathbf{U}=\left[\begin{array}{c} U_{x}\\ U_{y}\\ U_{z} \end{array}\right]_{1}\qquad\mathbf{V}=\left[\begin{array}{c} V_{x}\\ V_{y}\\ V_{z} \end{array}\right]_{1}\end{split}\]

and their inner product is given by


Give two-forms \(\mathbf{U}=U_{yz}\,dy\wedge dz+U_{zx}\,dz\wedge dx+U_{xy}\,dx\wedge dy\) and \(\mathbf{V}=U_{yz}\,dy\wedge dz+U_{zx}\,dz\wedge dx+U_{xy}\,dx\wedge dy\), again we will express them as column vectors

\[\begin{split}\mathbf{U}=\left[\begin{array}{c} U_{yz}\\ U_{zx}\\ U_{xy} \end{array}\right]_{2}\qquad\mathbf{V}=\left[\begin{array}{c} V_{yz}\\ V_{zx}\\ V_{xy} \end{array}\right]_{2}\end{split}\]

and their inner product is given by


The Hodge star must be such that we must have the following properties satisfied for one-forms


Clearly, since \(\star dx=dy\wedge dz\), \(\star dy=dz\wedge dx\) and \(\star dz=dx\wedge dy\), then for a one-form \(\mathbf{U}\)

\[\begin{split}\mathbf{U}=\left[\begin{array}{c} U_{x}\\ U_{y}\\ U_{z} \end{array}\right]_{1}\qquad\star\mathbf{U}=\left[\begin{array}{c} U_{x}\\ U_{y}\\ U_{z} \end{array}\right]_{2}\end{split}\]

and for a two-form

\[\begin{split}\mathbf{U}=\left[\begin{array}{c} U_{yz}\\ U_{zx}\\ U_{xy} \end{array}\right]_{2}\qquad\star\mathbf{U}=\left[\begin{array}{c} U_{yz}\\ U_{zx}\\ U_{xy} \end{array}\right]_{1}\end{split}\]

General (Non-Euclidean) Space

In this case the metric (inner product between one-forms) is given by

\[\begin{split}g^{ij}=\left[\begin{array}{ccc} g^{xx} & g^{xy} & g^{zx}\\ g^{xy} & g^{yy} & g^{yz}\\ g^{zx} & g^{yz} & g^{zz} \end{array}\right]\end{split}\]

The inner product between two one-forms is then given by

\[\begin{split}\mathbf{U}\cdot\mathbf{V}=\mathbf{U}^{T}g\mathbf{V}=\left[\begin{array}{ccc} U_{x} & U_{y} & U_{z}\end{array}\right]_{1}\left[\begin{array}{ccc} g^{xx} & g^{xy} & g^{zx}\\ g^{xy} & g^{yy} & g^{yz}\\ g^{zx} & g^{yz} & g^{zz} \end{array}\right]\left[\begin{array}{c} V_{x}\\ V_{y}\\ V_{z} \end{array}\right]_{1}\end{split}\]

Expanding the matrix expressions

\[\begin{split}\begin{aligned} \mathbf{U}\cdot\mathbf{V} & =g^{xx}U_{x}V_{x}+g^{yy}U_{y}V_{y}+g^{zz}U_{z}V_{z}+\\ & \quad+g^{xy}\left(U_{x}V_{y}+U_{y}V_{x}\right)+g^{yz}\left(U_{y}V_{z}+U_{z}V_{y}\right)+g^{zx}\left(U_{x}V_{z}+U_{z}V_{x}\right)\end{aligned}\end{split}\]

To compute the inner product between two-forms we need to consider terms of the form \(\left(dx^{i}\wedge dx^{j}\right)\cdot\left(dx^{\alpha}\wedge dx^{\beta}\right)\). The wedge product is an anti-symmetric tensor product, it is given by

\[dx^{i}\wedge dx^{j}=dx^{i}\otimes dx^{j}-dx^{j}\otimes dx^{i}\]

The inner product between two forms is then given by

\[\begin{aligned} \left(dx^{i}\wedge dx^{j}\right)\cdot\left(dx^{\alpha}\wedge dx^{\beta}\right) & =\left(dx^{i}\otimes dx^{j}-dx^{j}\otimes dx^{i}\right)\cdot\left(dx^{\alpha}\otimes dx^{\beta}-dx^{\beta}\otimes dx^{\alpha}\right)\end{aligned}\]

where up to a normalization constant \(N\)

\[\begin{split}\begin{aligned} \left(dx^{i}\otimes dx^{j}\right)\cdot\left(dx^{\alpha}\otimes dx^{\beta}\right) & =N\left(dx^{i}\cdot dx^{\alpha}\right)\left(dx^{j}\cdot dx^{\beta}\right)\\ & =Ng^{i\alpha}g^{j\beta}\end{aligned}\end{split}\]

we can deduce that

\[\begin{split}\begin{aligned} \left(dx^{i}\wedge dx^{j}\right)\cdot\left(dx^{\alpha}\wedge dx^{\beta}\right) & =N\left(g^{i\alpha}g^{j\beta}-g^{i\beta}g^{j\alpha}-g^{j\alpha}g^{i\beta}+g^{j\beta}g^{i\alpha}\right)\\ & =2N\left(g^{i\alpha}g^{j\beta}-g^{i\beta}g^{j\alpha}\right)\end{aligned}\end{split}\]

The normalization constant will be chosen so that \(dx^{i}\wedge dx^{j}\) form an orthonormal basis, thefore \(N=1/2\), and

\[\left(dx^{i}\otimes dx^{j}\right)\cdot\left(dx^{\alpha}\otimes dx^{\beta}\right)=\frac{1}{2}g^{i\alpha}g^{j\beta}\]

The inner product between two two-forms then will be given by

\[\begin{split}\mathbf{U}\cdot\mathbf{V}=\left[\begin{array}{ccc} U_{yz} & U_{zx} & U_{xy}\end{array}\right]_{2}\left[\begin{array}{ccc} g^{yy}g^{zz}-g^{yz\,2} & g^{yz}g^{zx}-g^{xy}g^{zz} & g^{xy}g^{yz}-g^{zx}g^{yy}\\ g^{yz}g^{zx}-g^{xy}g^{zz} & g^{zz}g^{xx}-g^{zx\,2} & g^{zx}g^{xy}-g^{yz}g^{xx}\\ g^{xy}g^{yz}-g^{zx}g^{yy} & g^{zx}g^{xy}-g^{yz}g^{xx} & g^{xx}g^{yy}-g^{xy\,2} \end{array}\right]\left[\begin{array}{c} V_{yz}\\ V_{zx}\\ V_{xy} \end{array}\right]_{2}\end{split}\]

where we have used the fact that

\[\begin{split}\begin{aligned} \left(dx^{i}\wedge dx^{j}\right)\cdot\left(dx^{\alpha}\wedge dx^{\beta}\right) & =\left(dx^{i}\cdot dx^{\alpha}\right)\left(dx^{j}\cdot dx^{\beta}\right)-\left(dx^{i}\cdot dx^{\beta}\right)\left(dx^{j}\cdot dx^{\alpha}\right)\\ & =g^{i\alpha}g^{j\beta}-g^{i\beta}g^{j\alpha}\end{aligned}\end{split}\]

Expanding the matrix expression

\[\begin{split}\begin{aligned} \mathbf{U}\cdot\mathbf{V} & =g^{yy}g^{zz}U_{yz}V_{yz}+g^{zz}g^{xx}U_{zx}V_{zx}+g^{xx}g^{yy}U_{xy}V_{xy}+\\ & +g^{yz}g^{zx}\left(U_{yz}V_{zx}+U_{zx}V_{yz}\right)+g^{zx}g^{xy}\left(U_{zx}V_{xy}+U_{xy}V_{zx}\right)\\ & +g^{xy}g^{yz}\left(U_{xy}V_{yz}+U_{yz}V_{xy}\right)\end{aligned}\end{split}\]

The hodge-star must be such that


where \(\mu\) is

\[\mathbf{\mu}=\frac{1}{\sqrt{\det g^{ij}}}dx\wedge dy\wedge dz\]


\[\det g^{ij}=g^{xx}g^{yy}g^{zz}+2g^{xy}g^{yz}g^{zx}-g^{xx}g^{yz\,2}-g^{yy}g^{zx\,2}-g^{zz}g^{xy\,2}\]

The wedge product in 3D between two one-forms is

\[\begin{split}\begin{aligned} \mathbf{U} & =U_{x}\,dx+U_{y}\,dy+U_{z}\,dz\\ \mathbf{V} & =V_{x}\,dx+V_{y}\,dy+V_{z}\,dz\\ \mathbf{U}\wedge\mathbf{V} & =\left(U_{y}V_{z}-U_{z}V_{y}\right)\,dy\wedge dz+\\ & +\left(U_{z}V_{x}-U_{x}V_{z}\right)\,dz\wedge dx+\\ & +\left(U_{x}V_{y}-U_{y}V_{x}\right)\,dx\wedge dy\end{aligned}\end{split}\]

and between a one-form and a two form is

\[\begin{split}\begin{aligned} \mathbf{U} & =U_{x}\,dx+U_{y}\,dy+U_{z}\,dz\\ \mathbf{V} & =V_{yz}\,dy\wedge dz+V_{zx}\,dz\wedge dx+V_{xy}\,dx\wedge dy\\ \mathbf{U}\wedge\mathbf{V} & =\left(U_{x}V_{yz}+U_{y}V_{zx}+U_{z}V_{xy}\right)dx\wedge dy\wedge dz\end{aligned}\end{split}\]

Now we need to look for a \(\mathbf{V}_{2}=\star\mathbf{U}_{1}\) such that \(\mathbf{U}_{1}\wedge\mathbf{V}_{2}=(\mathbf{U}_{1}\cdot\mathbf{U}_{1})\mu_{3}\) (where the subscript in \(\mathbf{U}_{k}\) indicates its a k-form).

\[\begin{aligned} U_{x}V_{yz}+U_{y}V_{zx}+U_{z}V_{xy} & =\frac{1}{\sqrt{\det}}\left(g^{xx}U_{x}^{\phantom{x}2}+g^{yy}U_{y}^{\phantom{y}2}+g^{zz}U_{z}^{\phantom{y}2}+2g^{xy}U_{x}U_{y}+2g^{yz}U_{y}U_{z}+2g^{zx}U_{z}U_{x}\right)\end{aligned}\]

Clearly one way to satisfy that is

\[\begin{split}\begin{aligned} V_{yz} & =\frac{1}{\sqrt{\det}}\left(g^{xx}U_{x}+g^{xy}U_{y}+g^{zx}U_{z}\right)\\ V_{zx} & =\frac{1}{\sqrt{\det}}\left(g^{xy}U_{x}+g^{yy}U_{y}+g^{yz}U_{z}\right)\\ V_{xy} & =\frac{1}{\sqrt{\det}}\left(g^{zx}U_{x}+g^{yz}U_{y}+g^{zz}U_{z}\right)\end{aligned}\end{split}\]

Then the expression for the hodge-star is given by

\[\begin{split}\star\mathbf{U}_{1}=\frac{1}{\sqrt{\det}}\left[\begin{array}{ccc} g^{xx} & g^{xy} & g^{zx}\\ g^{xy} & g^{yy} & g^{yz}\\ g^{zx} & g^{yz} & g^{zz} \end{array}\right]\mathbf{U}_{1}\end{split}\]

Now to compute the hodge-star of a two-form, we need to look for \(\mathbf{V}_{1}=\star\mathbf{U}_{2}\) such that \(\mathbf{U}_{2}\wedge\mathbf{V}_{1}=(\mathbf{U}_{2}\cdot\mathbf{U}_{2})\mu_{3}\)


and the equality that needs to be satisfied becomes

\[\begin{split}\begin{aligned} g^{yy}g^{zz}U_{yz}^{\phantom{yz}2}+g^{zz}g^{xx}U_{zx}^{\phantom{yz}2}+g^{xx}g^{yy}U_{xy}^{\phantom{yz}2}+2g^{yz}g^{zx}U_{yz}U_{zx}+2g^{zx}g^{xy}U_{zx}U_{xy}+2g^{xy}g^{yz}U_{xy}U_{yz} & =\\ \frac{1}{\sqrt{\det}}\left(V_{x}U_{yz}+V_{y}U_{zx}+V_{z}U_{xy}\right)\end{aligned}\end{split}\]


\[\begin{split}\begin{aligned} V_{x} & =\frac{1}{\sqrt{\det}}\left(g^{yy}g^{zz}U_{yz}+g^{yz}g^{zx}U_{zx}+g^{xy}g^{yz}U_{xy}\right)\\ V_{y} & =\frac{1}{\sqrt{\det}}\left(g^{yz}g^{zx}U_{yz}+g^{zz}g^{xx}U_{zx}+g^{zx}g^{xy}U_{xy}\right)\\ V_{z} & =\frac{1}{\sqrt{\det}}\left(g^{xy}g^{yz}U_{yz}+g^{zx}g^{xy}U_{zx}+g^{xx}g^{yy}U_{xy}\right)\end{aligned}\end{split}\]

Clearly, the following satisfies the above equality

\[\begin{split}\star\mathbf{U}_{2}=\frac{1}{\sqrt{\det}}\left[\begin{array}{ccc} g^{yy}g^{zz}-g^{yz\,2} & g^{yz}g^{zx}-g^{xy}g^{zz} & g^{xy}g^{yz}-g^{zx}g^{yy}\\ g^{yz}g^{zx}-g^{xy}g^{zz} & g^{zz}g^{xx}-g^{zx\,2} & g^{zx}g^{xy}-g^{yz}g^{xx}\\ g^{xy}g^{yz}-g^{zx}g^{yy} & g^{zx}g^{xy}-g^{yz}g^{xx} & g^{xx}g^{yy}-g^{xy\,2} \end{array}\right]\mathbf{U}_{1}\end{split}\]


  • zero-form \(\mathbf{U}=U\,1\)

    \[\star\mathbf{U}=U\,\frac{1}{\sqrt{\det}}\,dx\wedge dy\wedge dz\]
  • one-form \(\mathbf{U}=U_{x}\,dx+U_{y}\,dy+U_{z}\,dz\)

    \[\begin{split}\star\mathbf{U}=\frac{1}{\sqrt{\det}}\left[\begin{array}{ccc} g^{xx} & g^{xy} & g^{zx}\\ g^{xy} & g^{yy} & g^{yz}\\ g^{zx} & g^{yz} & g^{zz} \end{array}\right]\mathbf{U}\end{split}\]
  • two-form \(\mathbf{U}=U_{yz}\,dy\wedge dz+U_{zx}\,dz\wedge dx+U_{xy}\,dx\wedge dy\)

    \[\begin{split}\star\mathbf{U}=\frac{1}{\sqrt{\det}}\,\left[\begin{array}{ccc} g^{yy}g^{zz}-g^{yz\,2} & g^{yz}g^{zx}-g^{xy}g^{zz} & g^{xy}g^{yz}-g^{zx}g^{yy}\\ g^{yz}g^{zx}-g^{xy}g^{zz} & g^{zz}g^{xx}-g^{zx\,2} & g^{zx}g^{xy}-g^{yz}g^{xx}\\ g^{xy}g^{yz}-g^{zx}g^{yy} & g^{zx}g^{xy}-g^{yz}g^{xx} & g^{xx}g^{yy}-g^{xy\,2} \end{array}\right]\mathbf{U}\end{split}\]
  • three-form \(\mathbf{U}=U_{xyz}\,dx\wedge dy\wedge dz\)

Mon, 27 Jul 2015 00:00:00 -0700 <![CDATA[Riemanian Metric and the Hodge-Star]]> Riemanian Metric and the Hodge-Star

In elementary geometry the inner product between two vectors \(\mathbf{U}=U^{i}\frac{\partial}{\partial x^{i}}\) and \(\mathbf{V}=V^{i}\frac{\partial}{\partial x^{i}}\) is given by


In the general case, the Riemanian metric is a type \((0,2)\) tensor field which at each point of the manifold \(m\in M\) satisfies

\[\begin{split}\begin{aligned} g_{m}(\mathbf{U},\mathbf{V}) & =g_{m}(\mathbf{V},\mathbf{U}) & \text{symmetry}\\ g_{m}(\mathbf{U},\mathbf{U}) & \geq0 & \text{postive-definiteness}\end{aligned}\end{split}\]

The metric tensor can be expressed as

\[g_{m}=g_{ij}(m)\,dx^{i}\otimes dx^{j}\]

where the components are given by

\[g_{ij}(m)=g_{m}\left(\frac{\partial}{\partial x^{i}},\frac{\partial}{\partial x^{j}}\right)\]

There exists an isomorphism between \(T_{m}M\) and \(T_{m}^{*}M\). Let \(U\) be a vector field (\((1,0)\) tensor field), and let \(\alpha\) be a one-form (\((0,1)\) tensor field). Then the isomorphism is expressed as

\[\alpha_{i}=g_{ij}U^{j}\quad U^{i}=g^{ij}\alpha_{j}\]

where \(g\) with raised indices is the inverse of \(g\) with lowered indices, namely


Thus the flat \((\flat)\) and sharp \((\sharp)\) oprators are easily expressible in terms of the metric

\[\left(U^{\flat}\right)_{i}=g_{ij}U^{j}\qquad U^{i}=g^{ij}\left(U^{\flat}\right)_{j}\]

From now on the convention that we will follow will be that components with raised indices will be the components of vector fields, and the components with lowered indices will be the components of one-forms. Thus \(V^{i}\) will denote the components of a vector field, and \(V_{i}\) will denote the components of a one-form.

Eucledian Space

The flat metric is

\[\begin{split}g=\left[\begin{array}{cc} 1 & 0\\ 0 & 1 \end{array}\right]\end{split}\]

Given two one-forms

\[\begin{split}\mathbf{U}=\left[\begin{array}{c} U_{x}\\ U_{y} \end{array}\right]\qquad\mathbf{V}=\left[\begin{array}{c} V_{x}\\ V_{y} \end{array}\right]\end{split}\]

their inner product is given by


The hodge-star must be such that we must have the following propreties satisfied for one-forms

\[\begin{split}\begin{aligned} \mathbf{U}\cdot\star\mathbf{U} & = & 0\\ \mathbf{U}\cdot\mathbf{U} & = & \star\mathbf{U}\cdot\star\mathbf{U}\end{aligned}\end{split}\]


\[\begin{split}\star\mathbf{U}=\left[\begin{array}{c} -U_{y}\\ U_{x} \end{array}\right]=\left[\begin{array}{cc} 0 & -1\\ 1 & 0 \end{array}\right]\mathbf{U}\end{split}\]

satisfies these properties.

General (Non-Eucledian) Space

In this case the symmetric bilinear form defining the inner product between two vectors

\[\begin{split}g_{ij}=\left[\begin{array}{cc} g_{xx} & g_{xy}\\ g_{xy} & g_{yy} \end{array}\right]\end{split}\]

The correspond inner product between two one-forms on the other hand will be given by

\[\begin{split}g^{ij}=\left[\begin{array}{cc} g^{xx} & g^{xy}\\ g^{xy} & g^{yy} \end{array}\right]\end{split}\]

and the two are inverses of each other

\[\begin{split}\left[\begin{array}{cc} g_{xx} & g_{xy}\\ g_{xy} & g_{yy} \end{array}\right]=\left[\begin{array}{cc} g^{xx} & g^{xy}\\ g^{xy} & g^{yy} \end{array}\right]^{-1}\end{split}\]

The inner product between two one-forms is then given by

\[\begin{split}\mathbf{U}\cdot\mathbf{V}=\mathbf{U}^{T}g\mathbf{V}=\left[\begin{array}{cc} U_{x} & U_{y}\end{array}\right]\left[\begin{array}{cc} g^{xx} & g^{xy}\\ g^{xy} & g^{yy} \end{array}\right]\left[\begin{array}{c} V_{x}\\ V_{y} \end{array}\right]\end{split}\]

Expanding the matrix expressions


The hodge-star must be such that


where \(\mu\) is the normalized volume form, and in 2D it is given by the following two equivalent expressions

\[\mathbf{\mu}=\sqrt{\det g_{ij}}dx\wedge dy\]
\[\mathbf{\mu}=\frac{1}{\sqrt{\det g^{ij}}}dx\wedge dy\]

Since we are working with one-forms the relevant inner product will be the one between one-forms, and we will use the second expression.

The wedge product in 2D is given in components as

\[\mathbf{U}\wedge\mathbf{V}=U_{i}V_{j}\,dx^{i}\wedge dx^{j}\]
\[\mathbf{U}\wedge\mathbf{V}=\left(U_{x}V_{y}-U_{y}V_{x}\right)dx\wedge dy\]

Now we need to look for a \(\mathbf{V}=\star\mathbf{U}\) such that \(\mathbf{U}\wedge\mathbf{V}=(\mathbf{U}\cdot\mathbf{U})\mu\), or

\[U_{x}V_{y}-U_{y}V_{x}=\frac{1}{\sqrt{\det g^{ij}}}\left(g^{xx}U_{x}^{\phantom{x}2}+2g^{xy}U_{x}U_{y}+g^{yy}U_{y}^{\phantom{y}2}\right)\]

Clearly one way to satify that is

\[\begin{split}\begin{aligned} V_{x} & = & -g^{yy}U_{y}-g^{xy}U_{x}\\ V_{y} & = & g^{xx}U_{x}+g^{xy}U_{y}\end{aligned}\end{split}\]

and noting that \(\det g^{ij}=g^{xx}g^{yy}-g^{xy\,2}\), then the expression for the hodge-star is given by

\[\begin{split}\star\mathbf{U}=\mathbf{V}=\frac{1}{\sqrt{g^{xx}g^{yy}-g^{xy\,2}}}\left[\begin{array}{c} -g^{yy}U_{y}-g^{xy}U_{x}\\ g^{xx}U_{x}+g^{xy}U_{y} \end{array}\right]\end{split}\]

which for the flat metric (\(g^{xx}=g^{yy}=1\) and \(g^{xy}=0\)) reverts to our expression in Eucledian space.

To check that this expression of the hodge star is indeed indempotent (up to a sign), we apply it twice

\[\begin{split}\begin{aligned} \star\star\mathbf{U} & =\frac{1}{\left|\det g^{ij}\right|}\left[\begin{array}{c} -g^{yy}\left(g^{xx}U_{x}+\cancel{g^{xy}U_{y}}\right)-g^{xy}\left(-\cancel{g^{yy}U_{y}}-g^{xy}U_{x}\right)\\ g^{xx}\left(-g^{yy}U_{y}-\bcancel{g^{xy}U_{x}}\right)+g^{xy}\left(\bcancel{g^{xx}U_{x}}+g^{xy}U_{y}\right) \end{array}\right]\\ & =\frac{1}{\left|\det g^{ij}\right|}\left[\begin{array}{c} -g^{yy}g^{xx}U_{x}+g^{xy\,2}U_{x}\\ -g^{xx}g^{yy}U_{y}+g^{xy\,2}U_{y} \end{array}\right]\\ & =-\frac{g^{xx}g^{yy}-g^{xy\,2}}{\left|\det g^{ij}\right|}\left[\begin{array}{c} U_{x}\\ U_{y} \end{array}\right]\\ & =-\left[\begin{array}{c} U_{x}\\ U_{y} \end{array}\right]=-\mathbf{U}\end{aligned}\end{split}\]


  • zero-form \(\mathbf{U}=U\,1\)

    \[\star\mathbf{U}=U\,\frac{1}{\sqrt{\det g^{ij}}}\,dx\wedge dy\]
  • one-form \(\mathbf{U}=U_{x}\,dx+U_{y}\,dy\)

    \[\begin{split}\star\mathbf{U}=\frac{1}{\sqrt{\det g^{ij}}}\left[\begin{array}{c} -g^{yy}U_{y}-g^{xy}U_{x}\\ g^{xx}U_{x}+g^{xy}U_{y} \end{array}\right]\end{split}\]
    \[\begin{split}\star\mathbf{U}=\frac{1}{\sqrt{\det g^{ij}}}\left[\begin{array}{cc} -g^{xy} & -g^{yy}\\ g^{xx} & g^{xy} \end{array}\right]\mathbf{U}\end{split}\]
  • two-form \(\mathbf{U}=U\,dx\wedge dy\)

    \[\star\mathbf{U}=U\,\sqrt{\det g^{ij}}\]
Wed, 22 Jul 2015 00:00:00 -0700 <![CDATA[Discrete Contraction and Inner Product]]> Discrete Contraction and Inner Product

Let \(\alpha\) and \(\beta\) be two 1-forms. Their inner product can be expressed as


where \(\langle\cdot,\cdot\rangle\) represents the pairing between a vector field and a one form. In this context one can think of vector fields as row matrices, and forms as row matrices.

The sharp operator is a mapping from forms to vector fields, and the flat operator is its inverse.

\[\begin{split}\begin{aligned} \flat: & \quad & TM\to T^{*}M\\ \sharp: & \quad & T^{*}M\to TM\end{aligned}\end{split}\]

Contraction with a one form

The contraction with a one form is simply given by


Contraction with a two form

Suppose the two form can be expressed as the wedge product of two one forms such as \(\beta\wedge\gamma\). Then

\[\begin{split}\begin{aligned} \langle\alpha^{\sharp},\beta\wedge\gamma\rangle & = & \langle\alpha^{\sharp},\beta\rangle\wedge\gamma-\langle\alpha^{\sharp},\gamma\rangle\wedge\beta\\ & = & \star(\alpha\wedge\star\beta)\wedge\gamma-\star(\alpha\wedge\star\gamma)\wedge\beta\end{aligned}\end{split}\]
Thu, 04 Jun 2015 00:00:00 -0700 <![CDATA[The Hodge Star Operator and the Dual Mesh]]> The Hodge Star Operator and the Dual Mesh

On an \(n\)-dimensional Riemanian manifold the continuous Hodge Star operator establishes an isomorphism between \(k\)-forms and \((n-k)\)-forms.

\[\begin{split}\begin{aligned} \star^{k}: & \quad & \Lambda^{k}\to\Lambda^{n-k}\\ \star^{n-k}: & \quad & \Lambda^{n-k}\to\Lambda^{k}\end{aligned}\end{split}\]

It satisfies the following identity


where \(\mathbf{id}^{k}\) is the identity map between \(k\)-forms (it maps the form to itself)

\[\begin{split}\begin{aligned} \mathbf{id}^{k}: & \quad & \Lambda^{k}\to\Lambda^{k}\\ & & \alpha\mapsto\alpha\end{aligned}\end{split}\]

The hodge star operators is its own inverse up to a sign. Let us consider only the positive case for now (an analogous argument will hold for the negative case).

We want to study under what conditions the discrete Hodge star operator satisfies the discrete counterpart of the above equation


The discrete Hodge star operator acting on primal forms is given by:


The discrete Hodge star acting on dual forms, on the other hand, is given by


Let us assume that eq. (\ref{eq:discrete} ) holds. Then in terms of the basis functions it becomes


The basis functions by their own definition must satisfy


Additionally, we desire the dual basis functions to be linear combinations of the hodge-star of the primal ones, therefore


where the matrix \(\mathbf{A}\) represents the transform between the two basis sets. The hodge star is required in order for the dual basis functions to correspond to the correct simplex dimensions. For example if we take the basis functions for 0-forms in 2d corresponding to vertices, the dual basis functions must correspond to the corresponding elements on the dual mesh wich are faces, and therefore must be 2-forms.

The following sequence of equations are equivalent

\[\begin{split}\begin{aligned} \langle\sigma_{i},\star\tilde{\phi}_{j}\rangle\langle\tilde{\sigma}^{j},\star\phi^{k}\rangle & =\delta_{i}^{k}\\ A_{jn}\langle\sigma_{i},\star\star\phi^{n}\rangle\langle\tilde{\sigma}^{j},\star\phi^{k}\rangle & =\delta_{i}^{k}\\ A_{jn}\delta_{i}^{n}\langle\tilde{\sigma}^{j},\star\phi^{k}\rangle & =\delta_{i}^{k}\\ A_{ji}\langle\tilde{\sigma}^{j},\star\phi^{k}\rangle & =\delta_{i}^{k}\\ A_{ji}H^{jk} & =\delta_{i}^{k}\\ A_{ki}H^{kj} & =\delta_{i}^{j}\quad(k\leftrightarrow j)\end{aligned}\end{split}\]

and for eq. (\ref{eq:linearity} )

\[\begin{split}\begin{aligned} \langle\tilde{\sigma}^{j},\tilde{\phi}_{i}\rangle & =\delta_{i}^{j}\\ A_{ik}\langle\tilde{\sigma}^{j},\star\phi^{k}\rangle & =\delta_{i}^{j}\\ A_{ik}H^{jk} & =\delta_{i}^{j}\end{aligned}\end{split}\]

In matrix form


and since \(\mathbf{H}\) is of full rank, this implies that


We have shown that eq. (\ref{eq:discrete} ) (the Hodge star and its dual are exact inverses) and eq. (\ref{eq:linearity} ) (dual basis functions are linear combinations of primal ones) are completely consistent with each other (can it be shown that one implies the other?).

Therefore, the primal basis functions \(\phi\) and the dual simplices \(\tilde{\sigma}\) completely determine the dual basis functions, and they are implicitly given by



Wed, 22 Apr 2015 00:00:00 -0700 <![CDATA[Invariance of the Discrete Wedge Product]]> Invariance of the Discrete Wedge Product

Consider a manifold \(\mathcal{M}\) and its polyhedrization in the context of discrete exterior calculus, and let \(\sigma_{i}\) denote the simplices and \(\phi^{i}\) the corresponding basis functions, such that the pairing given by integration is natural


Given the projection \(\mathcal{P}\) and reconstruction operators \(\mathcal{R}\),

\[\begin{split}\begin{aligned} \mathcal{P}:\quad & \Lambda^{k}\to\bar{\Lambda}^{k}\\ & \alpha\mapsto\bar{\alpha}_{i}=\int_{\sigma_{i}}\alpha\end{aligned}\end{split}\]
\[\begin{split}\begin{aligned} \mathcal{R}:\quad & \bar{\Lambda}^{k}\to\Lambda^{k}\\ & \bar{\alpha}_{i}\mapsto\alpha=\bar{\alpha}_{i}\phi^{i}\end{aligned}\end{split}\]

the discrete wedge product is defined as a three tensor




Consider now a second manifold \({\cal M}^{\prime}\) and a diffeomorphism

\[\varphi:\quad{\cal M}\to{\cal M}^{'}\]

which induces a new polyderization with the corresponding simplices \(\sigma^{\prime}=\varphi_{*}\sigma\) and basis functions \(\phi^{\prime}=\varphi_{*}\phi\).

Given a simplex \(\sigma\) and a form \(\phi\), the following identities hold for push-forwards and pull-backs:


The discrete wedge product on this new manifold will be elementwise equal to the wedge product on the original manifold. To see why this is the case:

\[\begin{split}\begin{aligned} \mathbf{W}_{\phantom{\prime}i}^{\prime\phantom{i}jk} & =\int_{\sigma_{i}^{\prime}}\phi^{\prime j}\wedge\phi^{\prime k}\\ & =\int_{\varphi_{*}\sigma_{i}}\left(\varphi_{*}\phi^{j}\right)\wedge\left(\varphi_{*}\phi^{k}\right)\\ & =\int_{\varphi_{*}\sigma_{i}}\varphi_{*}\left(\phi^{j}\wedge\phi^{k}\right)\quad(\mbox{eq.}\ref{eq:wedge})\\ & =\int_{\sigma_{i}}\phi^{j}\wedge\phi^{k}\quad(\mbox{eq.}\ref{eq:integration})\\ & =\mathbf{W}_{i}^{\phantom{i}jk}\end{aligned}\end{split}\]

So we have shown that the discrete wedge product is indeed invariant under diffeomorphisms.


This result is consistent with theory in differential geometry, where the wedge product is thought of as a pureley topological operator that is independent of the metric.

Mon, 20 Apr 2015 00:00:00 -0700 <![CDATA[Bitcoin and Silicon Valley]]> Bitcoin and Silicon Valley

This weekend I attended the Bitcoin Job Fair in Sunnyvale, CA. There I had the pleasure of meeting Balaji Srinivasan in person, and he was kind enough to let me take a picture with him.


If you don’t know who he is, I strongly recommend that you watch his talk titled Silicon Valley’s Ultimate Exit. I have been a big fan ever since I saw that talk. One of the observations that Srinivasan makes is about how Silicon Valley technologies are already disrupting major centers of power in the United States, and replacing what until now used to be brick and mortar entities that relied primarily on paper to transmit and store information (appropriately referred to as the Paper Belt), with digital entities based on the internet that have small to no presence in the physical world.

There is hardly an industry that has not been affected by Silicon Valley now. The entertainment industry which until recently derived its profits mostly from the middle men that distributed content from producers to consumers has been transformed by technologies like iTunes and BitTorrent that can serve the same function at close to zero cost. The education industry is being disrupted right now by companies like Khan Academy or Coursera which exist solely on the internet and require no costly physical campus. But the example that was most interesting to me, to which there was only a brief reference in that talk from 2013, was Bitcoin.

Fast forward to 2015, and it is quite exciting to see that Srinivasan has now turned his attention fully onto Bitcoin. His presentation during the fair was entirely focused on the cryptocurrency. One of the more striking facts mentioned was about how Bitcoin is now bigger than Google itself by several significant metrics. In terms of energy consumption, the Bitcoin mining network already draws in more energy than all of Google’s datacenters combined. In terms of computational power, Bitcoin has surpassed Google’s server farms too. While in terms of market capitalization, Bitcoin still trails Google by several orders of magnitude, one has to wonder how long it will be before Google is left behind in that category too by what just a few years ago used to be an obscure open source project, dismissed by most as a toy currency. Given its meteoric rise, one can only wonder where Bitcoin will be in the coming years.

Sun, 19 Apr 2015 00:00:00 -0700 <![CDATA[Wedge Product in Coordinates]]> Wedge Product in Coordinates

Given a manifold \(M\), the wedge product is a map that constructs higher order forms


The wedge product has the following properties:

  • \(\alpha\wedge\beta\) is associative: \(\alpha\wedge(\beta\wedge\gamma)=(\alpha\wedge\beta)\wedge\gamma\)

  • \(\alpha\wedge\beta\) is bilinear in \(\alpha\) and \(\beta\):

  • \(\alpha\wedge\beta\) is anticommutative: \(\alpha\wedge\beta=(-1)^{kl}\beta\wedge\alpha\), where \(\alpha\) is a \(k\)-form and \(\beta\) is an \(l\)-form.

In the discrete setting we will only be able to preserve some of these continuous properties. Namely, the bilinearity and anticommutativity will be preserved exactly, whereas the associativity will be satisfied only in the limit where the mesh size tends to zero (\(h\to0\)) and will not be exact.

The wedge product is a an operator that is independent of the metric, i.e. it is a homomorphism under a pull-back:


Consider the 2D case

\[\begin{split}\begin{aligned} \alpha & =\alpha\\ \beta & =\beta_{x}dx+\beta_{y}dy\\ \alpha\wedge\beta & =\alpha\beta_{x}dx+\alpha\beta_{y}dy\end{aligned}\end{split}\]
\[\begin{split}\begin{aligned} \alpha & =\alpha_{x}dx+\alpha_{y}dy\\ \beta & =\beta_{x}dx+\beta_{y}dy\\ \alpha\wedge\beta & =\left(\alpha_{x}\beta_{y}-\beta_{x}\alpha_{y}\right)dx\wedge dy\end{aligned}\end{split}\]

The 3D case, on the other hand, will be

\[\begin{split}\begin{aligned} \alpha & =\alpha\\ \beta & =\beta_{x}dx+\beta_{y}dy+\beta_{z}dz\\ \alpha\wedge\beta & =\alpha\beta_{x}dx+\alpha\beta_{y}dy+\alpha\beta_{z}dz\end{aligned}\end{split}\]
\[\begin{split}\begin{aligned} \alpha= & \alpha_{x}dx+\alpha_{y}dy+\alpha_{z}dz\\ \beta= & \beta_{x}dx+\beta_{y}dy+\beta_{z}dz\\ \alpha\wedge\beta= & (\alpha_{x}\beta_{y}-\alpha_{y}\beta_{x})dx\wedge dy+\\ + & (\alpha_{y}\beta_{z}-\alpha_{z}\beta_{y})dy\wedge dz+\\ + & (\alpha_{z}\beta_{x}-\alpha_{x}\beta_{z})dz\wedge dx\end{aligned}\end{split}\]
\[\begin{split}\begin{aligned} \alpha & =\alpha_{x}dx+\alpha_{y}dy+\alpha_{z}dz\\ \beta & =\beta_{xy}dx\wedge dy+\beta_{yz}dy\wedge dz+\beta_{zx}dz\wedge dx\\ \alpha\wedge\beta & =(\alpha_{x}\beta_{yz}+\alpha_{y}\beta_{zx}+\alpha_{z}\beta_{xy})dx\wedge dy\wedge dz\end{aligned}\end{split}\]
Sat, 11 Apr 2015 00:00:00 -0700 <![CDATA[Contraction in Coordinates]]> Contraction in Coordinates

Given a manifold \(M\), the interior product is defined as the contraction of a differential form with a vector field.


The contraction of a one-form \(\alpha=\alpha_{i}e^{i}\) with a vector field \(X=X^{i}e_{i}\) is given by

\[\begin{split}\begin{aligned} \mathbf{i}_{X}\alpha & =\langle X,\alpha\rangle\\ & =\alpha_{i}X^{j}\langle e^{i},e_{j}\rangle\\ & =\alpha_{i}X^{j}\delta_{j}^{i}\\ & =\alpha_{i}X^{i}\end{aligned}\end{split}\]

The contraction of a two-form \(\beta=\beta_{ij}e^{i}\wedge e^{j}\) with the same vector field is given by

\[\begin{split}\begin{aligned} \mathbf{i}_{X}\beta & =\beta_{ij}\mathbf{i}_{X}(e^{i}\wedge e^{j})\\ & =\beta_{ij}\left((\mathbf{i}_{X}e^{i})\wedge e^{j}-e^{i}\wedge(\mathbf{i}_{X}e^{j})\right)\\ & =\beta_{ij}\left(X^{i}e^{j}-X^{j}e^{i}\right)\end{aligned}\end{split}\]

The contraction of a three-form \(\gamma=\gamma{}_{ijk}e^{i}\wedge e^{j}\wedge e^{k}\) is going to be

\[\begin{split}\begin{aligned} \mathbf{i}_{X}\gamma & =\gamma_{ijk}\mathbf{i}_{X}(e^{i}\wedge e^{j}\wedge e^{k})\\ & =\gamma_{ijk}\left(\mathbf{i}_{X}e^{i}\wedge e^{j}\wedge e^{k}-e^{i}\wedge\mathbf{i}_{X}e^{j}\wedge e^{k}+e^{i}\wedge e^{j}\wedge\mathbf{i}_{X}e^{k}\right)\\ & =\gamma_{ijk}\left(X^{i}e^{j}\wedge e^{k}-X^{j}e^{i}\wedge e^{k}+X^{k}e^{i}\wedge e^{j}\right)\end{aligned}\end{split}\]

In 2D the explicit formula for the contractions are

\[\begin{split}\begin{aligned} \alpha & =\alpha_{x}dx+\alpha_{y}dy\\ \mathbf{i}_{X}\alpha & =X^{x}\alpha_{x}+X^{y}\alpha_{y}\end{aligned}\end{split}\]
\[\begin{split}\begin{aligned} \beta & =\beta_{xy}\,dx\wedge dy\\ \mathbf{i}_{X}\beta & =-\beta_{xy}X^{y}dx+\beta_{xy}X^{x}dy\end{aligned}\end{split}\]

If 3D, on the other hand, the contractions become

\[\begin{split}\begin{aligned} \alpha & =\alpha_{x}dx+\alpha_{y}dy+\alpha_{z}dz\\ \mathbf{i}_{X}\alpha & =X^{x}\alpha_{x}+X^{y}\alpha_{y}+X^{z}\alpha_{z}\end{aligned}\end{split}\]
\[\begin{split}\begin{aligned} \beta & =\beta_{xy}dx\wedge dy+\beta_{yz}dy\wedge dz+\beta_{zx}dz\wedge dx\\ \mathbf{i}_{X}\beta & =\left(\beta_{zx}X^{z}-\beta_{xy}X^{y}\right)dx+\\ & +\left(\beta_{xy}X^{x}-\beta_{yz}X^{z}\right)dy+\\ & +\left(\beta_{yz}X^{y}-\beta_{zx}X^{x}\right)dz\end{aligned}\end{split}\]
\[\begin{split}\begin{aligned} \gamma & =\gamma_{xyz}dx\wedge dy\wedge dz\\ \mathbf{i}_{X}\gamma & =\gamma_{xyz}\left(X^{z}dx\wedge dy+X^{x}dy\wedge dz+X^{y}dz\wedge dx\right)\end{aligned}\end{split}\]
Thu, 18 Dec 2014 00:00:00 -0800 <![CDATA[Lie Derivative in Coordinates]]> Lie Derivative in Coordinates

Here we will study the coordinate dependent realization of the Lie derivative. Let \(f\) be a scalar (zero-form), \(X\) a vector field, and \(\alpha\) be a one-form, where

\[\begin{split}\begin{aligned} X=X^{i}e_{i}=X^{i}\frac{\partial}{\partial x^{i}}\\ \alpha=\alpha_{i}e^{i}=\alpha_{i}dx^{i}\end{aligned}\end{split}\]

\(\{e_{1},\dots,e_{n}\}=\{\frac{\partial}{\partial x^{1}},\dots,\frac{\partial}{\partial x^{n}}\}\) form the basis that span the space of vector fields, and \(\{e^{1},\dots,e^{n}\}=\{dx^{1},\dots,dx^{n}\}\) are the basis elements for the space of one-forms.

Using Cartan’s magic formula the Lie derivative is given algebraically by


In coordinate form the Lie derivative of a scalar, a one-form and a vector field become

\[\mathfrak{L}_{X}f=\mathbf{i}_{X}\mathbf{d}f=X^{i}\frac{\partial f}{\partial x^{i}}\]
\[\mathfrak{L}_{X}Y=[X,Y]=X^{i}\frac{\partial Y^{j}}{\partial x^{i}}\frac{\partial}{\partial x^{j}}-Y^{i}\frac{\partial X^{j}}{\partial x^{i}}\frac{\partial}{\partial x^{j}}\]

To compute the Lie derivative of a one form we use

\[\begin{aligned} \mathbf{d}\alpha & =\frac{\partial\alpha_{i}}{\partial x^{j}}\,dx^{j}\wedge dx^{i}\end{aligned}\]
\[\begin{split}\begin{aligned} \mathbf{i}_{X}\mathbf{d}\alpha & =X^{k}\mathbf{i}_{\frac{\partial}{\partial x^{k}}}\mathbf{d}\alpha\\ & =X^{k}\frac{\partial\alpha_{i}}{\partial x^{j}}\mathbf{i}_{\frac{\partial}{\partial x^{k}}}\left(dx^{j}\wedge dx^{i}\right)\\ & =X^{k}\frac{\partial\alpha_{i}}{\partial x^{j}}\left(\left(\mathbf{i}_{\frac{\partial}{\partial x^{k}}}dx^{j}\right)\wedge dx^{i}-\left(dx^{j}\wedge\mathbf{i}_{\frac{\partial}{\partial x^{k}}}dx^{i}\right)\right)\\ & =X^{k}\frac{\partial\alpha_{i}}{\partial x^{j}}\left(\delta_{k}^{j}dx^{i}-\delta_{k}^{i}dx^{j}\right)\\ & =X^{j}\frac{\partial\alpha_{i}}{\partial x^{j}}dx^{i}-X^{i}\frac{\partial\alpha_{i}}{\partial x^{j}}dx^{j}\end{aligned}\end{split}\]
\[\begin{split}\begin{aligned} \mathbf{i}_{X}\alpha & =X^{i}\mathbf{i}_{\frac{\partial}{\partial x^{i}}}\alpha\\ & =X^{i}\alpha_{j}\mathbf{i}_{\frac{\partial}{\partial x^{i}}}dx^{j}\\ & =X^{i}\alpha_{j}\delta_{i}^{j}\\ & =X^{i}\alpha_{i}\end{aligned}\end{split}\]
\[\begin{split}\begin{aligned} \mathbf{d}\mathbf{i}_{X}\alpha & =d\mathbf{d}\left(X^{i}\alpha_{i}\right)\\ & =\left(\frac{\partial X^{i}}{\partial x^{j}}\alpha_{i}+X^{i}\frac{\partial\alpha_{i}}{\partial x^{j}}\right)dx^{j}\end{aligned}\end{split}\]

Combining the terms above we obtain the results for the Lie derivative of a one-form

\[\mathfrak{L}_{X}\alpha=X^{j}\frac{\partial\alpha_{i}}{\partial x^{j}}dx^{i}+\frac{\partial X^{i}}{\partial x^{j}}\alpha_{i}dx^{j}\]
\[\mathfrak{L}_{X}\alpha=\left(X^{j}\frac{\partial\alpha_{i}}{\partial x^{j}}+\frac{\partial X^{j}}{\partial x^{i}}\alpha_{j}\right)dx^{i}\]

The Lie derivative is a derivation and so it must satisfy a product rule over the pairing of forms and vector fields


To see that the above expressions for the Lie derivative do indeed satisfy Eq.\ref{eq:derivation}

\[\begin{split}\begin{aligned} \mathfrak{L}_{X}\langle\alpha,Y\rangle & =\mathfrak{L}_{X}\left(\alpha_{i}Y^{i}\right)\\ & =X^{j}\frac{\partial}{\partial x^{j}}\left(\alpha_{i}Y^{i}\right)\\ & =X^{j}\frac{\partial\alpha_{i}}{\partial x^{j}}Y^{i}+X^{j}\alpha_{i}\frac{\partial Y^{i}}{\partial x^{j}}\end{aligned}\end{split}\]
\[\begin{split}\begin{aligned} \langle\mathfrak{L}_{X}\alpha,Y\rangle & =\langle X^{j}\frac{\partial\alpha_{i}}{\partial x^{j}}dx^{i}+\frac{\partial X^{i}}{\partial x^{j}}\alpha_{i}dx^{j},Y^{k}\frac{\partial}{\partial x^{k}}\rangle\\ & =X^{j}\frac{\partial\alpha_{i}}{\partial x^{j}}Y^{k}\langle dx^{i},\frac{\partial}{\partial x^{k}}\rangle+\frac{\partial X^{i}}{\partial x^{j}}\alpha_{i}Y^{k}\langle dx^{j},\frac{\partial}{\partial x^{k}}\rangle\\ & =X^{j}\frac{\partial\alpha_{i}}{\partial x^{j}}Y^{k}\delta_{k}^{i}+\frac{\partial X^{i}}{\partial x^{j}}\alpha_{i}Y^{k}\delta_{k}^{j}\\ & =X^{j}\frac{\partial\alpha_{i}}{\partial x^{j}}Y^{i}+\frac{\partial X^{i}}{\partial x^{j}}\alpha_{i}Y^{j}\end{aligned}\end{split}\]
\[\begin{split}\begin{aligned} \langle\alpha,\mathfrak{L}_{X}Y\rangle & =\langle\alpha_{k}dx^{k},X^{i}\frac{\partial Y^{j}}{\partial x^{i}}\frac{\partial}{\partial x^{j}}-Y^{i}\frac{\partial X^{j}}{\partial x^{i}}\frac{\partial}{\partial x^{j}}\rangle\\ & =\left(\alpha_{k}X^{i}\frac{\partial Y^{j}}{\partial x^{i}}-\alpha_{k}Y^{i}\frac{\partial X^{j}}{\partial x^{i}}\right)\langle dx^{k},\frac{\partial}{\partial x^{j}}\rangle\\ & =\left(\alpha_{k}X^{i}\frac{\partial Y^{j}}{\partial x^{i}}-\alpha_{k}Y^{i}\frac{\partial X^{j}}{\partial x^{i}}\right)\delta_{j}^{k}\\ & =\alpha_{j}X^{i}\frac{\partial Y^{j}}{\partial x^{i}}-\alpha_{j}Y^{i}\frac{\partial X^{j}}{\partial x^{i}}\end{aligned}\end{split}\]

After a simple relabeling of the dummy indices it can be shown that Eq.\ref{eq:derivation} is indeed satisfied.

Thus, in two dimensions it can be shown that

\[\begin{split}\begin{aligned} \mathfrak{L}_{X}f & =X^{x}\frac{\partial f}{\partial x}+X^{y}\frac{\partial f}{\partial y}\\ \mathfrak{L}_{X}Y & =[X^{x}\frac{\partial}{\partial x}+X^{y}\frac{\partial}{\partial y},Y^{x}\frac{\partial}{\partial x}+Y^{y}\frac{\partial}{\partial y}]\\ & =X^{x}\frac{\partial Y^{x}}{\partial x}\frac{\partial}{\partial x}+X^{y}\frac{\partial Y^{x}}{\partial y}\frac{\partial}{\partial x}+X^{x}\frac{\partial Y^{y}}{\partial x}\frac{\partial}{\partial y}+X^{y}\frac{\partial Y^{y}}{\partial y}\frac{\partial}{\partial y}-\\ & -Y^{x}\frac{\partial X^{x}}{\partial x}\frac{\partial}{\partial x}-Y^{y}\frac{\partial X^{x}}{\partial y}\frac{\partial}{\partial x}-Y^{x}\frac{\partial X^{y}}{\partial x}\frac{\partial}{\partial y}-Y^{y}\frac{\partial X^{y}}{\partial y}\frac{\partial}{\partial y}\\ & =\left(X^{x}\frac{\partial Y^{x}}{\partial x}+X^{y}\frac{\partial Y^{x}}{\partial y}-Y^{x}\frac{\partial X^{x}}{\partial x}-Y^{y}\frac{\partial X^{x}}{\partial y}\right)\frac{\partial}{\partial x}+\\ & +\left(X^{x}\frac{\partial Y^{y}}{\partial x}+X^{y}\frac{\partial Y^{y}}{\partial y}-Y^{x}\frac{\partial X^{y}}{\partial x}-Y^{y}\frac{\partial X^{y}}{\partial y}\right)\frac{\partial}{\partial y}\\ \mathfrak{L}_{X}\alpha & =X^{x}\frac{\partial\alpha_{x}}{\partial x}dx+X^{y}\frac{\partial\alpha_{x}}{\partial y}dx+X^{x}\frac{\partial\alpha_{y}}{\partial x}dy+X^{y}\frac{\partial\alpha_{y}}{\partial y}dy+\\ & +\frac{\partial X^{x}}{\partial x}\alpha_{x}dx+\frac{\partial X^{y}}{\partial x}\alpha_{y}dx+\frac{\partial X^{x}}{\partial y}\alpha_{x}dy+\frac{\partial X^{y}}{\partial y}\alpha_{y}dy\\ & =\left(X^{x}\frac{\partial\alpha_{x}}{\partial x}+X^{y}\frac{\partial\alpha_{x}}{\partial y}+\frac{\partial X^{x}}{\partial x}\alpha_{x}+\frac{\partial X^{y}}{\partial x}\alpha_{y}\right)dx+\\ & +\left(X^{x}\frac{\partial\alpha_{y}}{\partial x}+X^{y}\frac{\partial\alpha_{y}}{\partial y}+\frac{\partial X^{x}}{\partial y}\alpha_{x}+\frac{\partial X^{y}}{\partial y}\alpha_{y}\right)dy\end{aligned}\end{split}\]

Given a two form

\[\Omega=\omega\,dx\wedge dy\]

it can be shown that

\[\begin{split}\begin{aligned} \mathfrak{L}_{X}\Omega & =\mathbf{d}\mathbf{i}_{X}\Omega\\ & =\left(\frac{\partial}{\partial x}(\omega X^{x})+\frac{\partial}{\partial y}(\omega X^{y})\right)\,dx\wedge dy\end{aligned}\end{split}\]

Other quantities

\[\begin{split}\begin{aligned} \operatorname{div}(X) & =\frac{\partial X^{x}}{\partial x}+\frac{\partial X^{y}}{\partial y}\\ \operatorname{vort}(X) & =\frac{\partial X^{y}}{\partial x}-\frac{\partial X^{x}}{\partial y}\end{aligned}\end{split}\]
Wed, 10 Dec 2014 00:00:00 -0800 <![CDATA[Sample Probability Problem]]> Sample Probability Problem


Suppose we are given a set \(\Omega\) consisting of \(N\) different elements, and a map


on that set, which maps every item to either \(0\) or \(1\).

The fraction of items that are mapped to \(1\) is given by


Given a subset \(S\subseteq\Omega\) of size \(n\) define the sample proportion

\[\bar{p}=\frac{1}{n}\sum_{x\in S}f(x)\]

For all the subsets of size \(n\) find the expectation value and the variation of the random variable \(\bar{p}\).


Let us first compute

\[\begin{split}\begin{aligned} \operatorname{E}[f(x)] & =p\\ \operatorname{E}[f(x)^{2}] & =\frac{1}{N}\sum_{x\in\Omega}f(x)^{2}\\ & =\frac{1}{N}\sum_{x\in\Omega}f(x)\\ & =p\\\end{aligned}\end{split}\]
\[\begin{split}\begin{aligned} \operatorname{Var}[f(x)] & =\operatorname{E}[f(x)^{2}]-\operatorname{E}[f(x)]^{2}\\ & =p-p^{2}\\ & =p(1-p)\\\end{aligned}\end{split}\]

Let us define the set of all subsets of size \(n\)


where \(2^{\Omega}\) defines the powerset of \(\Omega\), and the size of the set is

\[|\Sigma_{n}|={N \choose n}\]

For example, if \(\Omega=\{1,2,3\}\) then

\[\begin{split}\begin{aligned} \Sigma_{0} & =\{\{\}\} & |\Sigma_{0}|=1\\ \Sigma_{1} & =\{\{1\},\{2\},\{3\}\} & |\Sigma_{1}|=3\\ \Sigma_{2} & =\{\{1,2\},\{1,3\},\{2,3\}\}\quad & |\Sigma_{2}|=3\\ \Sigma_{3} & =\{\{1,2,3\}\} & |\Sigma_{0}|=1\\\end{aligned}\end{split}\]

Now the expectation value and variance that we are required to compute are given in terms of:

\[\begin{split}\begin{aligned} \operatorname{E}[\bar{p}] & =\frac{1}{|\Sigma_{n}|}\sum_{S\in\Sigma_{n}}\frac{1}{n}\sum_{x\in S}f(x)\\ & =\frac{1}{n|\Sigma_{n}|}\sum_{S\in\Sigma_{n}}\sum_{x\in S}f(x)\\ \operatorname{E}[\bar{p}^{2}] & =\frac{1}{|\Sigma_{n}|}\sum_{S\in\Sigma_{n}}\left(\frac{1}{n}\sum_{x\in S}f(x)\right)^{2}\\ & =\frac{1}{n^{2}|\Sigma_{n}|}\sum_{S\in\Sigma_{n}}\left(\sum_{x\in S}f(x)^{2}+2\sum_{x<y\,\in S}f(x)f(y)\right)\\\end{aligned}\end{split}\]

To compute the sums that appear in the expressions above we will use arguments of symmetry to show that the terms on the left and right hand side will be the same up to an integer constant

\[\sum_{S\in\Sigma_{n}}\sum_{x\in S}f(x)=A\sum_{x\in\Omega}f(x)\]

where \(A\) is a constant.

Counting the total number of terms on the left and the right it follows that

\[|\Sigma_{n}|n=A\,N\quad\therefore\quad A=\frac{n}{N}|\Sigma_{n}|\]


\[\sum_{S\in\Sigma_{n}}\sum_{x<y\,\in S}f(x)f(y)=B\sum_{x<y\,\in\Omega}f(x)f(y)\]

and counting the the terms on the left and the right

\[\Sigma_{n}\frac{n(n-1)}{2}=B\,\frac{N(N-1)}{2}\quad\therefore\quad B=|\Sigma_{n}|\frac{n(n-1)}{N(N-1)}\]

Since \(f(x)\) is either \(0\) or \(1\) it can be shown that

\[\begin{split}\begin{aligned} \sum_{x\in\Omega}f(x) & =pN\\ \sum_{x\in\Omega}f(x)^{2} & =pN\\ \sum_{x<y\,\in\Omega}f(x)f(y) & =\frac{pN(pN-1)}{2}\\\end{aligned}\end{split}\]

Finally, we obtain simplified expressions for the expectations

\[\begin{split}\begin{aligned} \operatorname{E}[\bar{p}] & =\frac{1}{n|\Sigma_{n}|}\sum_{S\in\Sigma_{n}}\sum_{x\in S}f(x)\\ & =\frac{1}{n|\Sigma_{n}|}A\sum_{x\in\Omega}f(x)\\ & =\frac{1}{n|\Sigma_{n}|}\frac{n}{N}|\Sigma_{n}|\sum_{x\in\Omega}f(x)\\ & =p\\ \operatorname{E}[\bar{p}^{2}] & =\frac{1}{n^{2}|\Sigma_{n}|}\sum_{S\in\Sigma_{n}}\left(\sum_{x\in S}f(x)^{2}+2\sum_{x<y\,\in S}f(x)f(y)\right)\\ & =\frac{1}{n^{2}|\Sigma_{n}|}\left(A\sum_{x\in\Omega}f(x)^{2}+2B\sum_{x<y\,\in\Omega}f(x)f(y)\right)\\ & =\frac{1}{n^{2}|\Sigma_{n}|}\left(|\Sigma_{n}|\frac{n}{N}\sum_{x\in\Omega}f(x)^{2}+2|\Sigma_{n}|\frac{n(n-1)}{N(N-1)}\sum_{x<y\,\in\Omega}f(x)f(y)\right)\\ & =\frac{p}{n}+\frac{1}{n}\frac{n-1}{N-1}p(pN-1)\end{aligned}\end{split}\]

Thus, we obtain the following final result

\[\begin{split}\begin{aligned} \operatorname{E}[\bar{p}] & =p\\ \operatorname{Var}[\bar{p}] & =\frac{p}{n}+\frac{1}{n}\frac{n-1}{N-1}p(pN-1)-p^{2}\end{aligned}\end{split}\]

It can be seen that this agrees with the Central Limit Theorem in the limit \(N\to\infty\)

Thu, 04 Dec 2014 00:00:00 -0800 <![CDATA[Flow Simulation]]> Flow Simulation

Here we present the simulation results of flows around a number of rigid and variable shape bodies. The different cases of Potential, Stokesian, and Navier-Stokes flow are considered.

The linear Potential and Stokes flows are solved using a simple panel method, whereas the Navier-Stokes flow is solved using the Immersed Boundary Method.

Pulsating Circle

We start with the simple example of a pulsating circle. It’s boundary is described by:

\[r(\theta) = 1+q_0(t)\]

with the periodic shape function

\[q_0(t) = \sin(t)\]

In python code this is expressed as

shapevariables = (
        lambda t: sin(t),

X = lambda s,q: cos(s) (1+q[0])
Y = lambda s,q: sin(s) (1+q[0])
bodies = (

The solution can easily be found from the continuity equaiton alone

\[u_r(\theta) \propto \frac{1}{r}\]

and the numerical simulation is run just for validation purposes. Whereas the panel methods work fine in this case, the immersed boundary method does not converge as it has not been designed to handle cases where the panel size on the boundary changes.

Below is the result of both the stokes and potential solver (the result is the same for both)


Translating Circle

We now consider the periodic translation of a rigid circle. Its boundary is described by:

\[\begin{split}x(\theta) &= \cos(\theta) + q_0(t)\\ y(\theta) &= \sin(\theta)\end{split}\]

with a periodic shape function

\[q_0(t) = \frac{1}{2}\sin(t)\]

The code listing

shapevariables = (
        lambda t: 0.5*sin(t),

X = lambda s,q: cos(s) + q[0]
Y = lambda s,q: sin(s)
bodies = (

Due to Stoke’s Paradox there does not exist a 2D creeping flow around a translating cylinder that vanishes at infinity. The Stokesian solver handles this problem gracefully and outputs a uniform flow for \(Re=0\). The immersed boundary method is used to solve for the \(Re=20\), \(Re=100\) and \(Re=1000\) flows. And the potential solver (assumes irrotational flow) is used to solve for \(Re=\infty\). The results are presented in the table below:

Re=0 Re=1000
im1 im4
Re=20 Potential, irrotatoinal
im2 im5

There is an inverse correlation between boundary thickness and Reynolds number as expected

\[\delta \propto \frac{1}{\sqrt{Re}}\]


You can view the simulations at higher resolution by Right Click->View Image (on Firefox) or using an equivalent method on another browser.

Mason & Burdick Amoeba

Here we take the example given by [MB99]. The boundary is simulated by:

\[r(\theta) = 1 + \epsilon (s_1(t) \cos(2\theta) + s_2(t) \sin(3\theta) )\]

By using an appropriate time evolution of the two shape modes we can achieve forward motion.


In the simulations below, the parameter is set to \(\epsilon=0.1\) and the following evolution in time of the shape modes is used

\begin{eqnarray*} s_{1}(t)&=&1-\cos t\\ s_{2}(t)&=&\sin t \end{eqnarray*}

In order for the immersed boundary method to converge we must also ensure that the area enclosed by the boundary remains constant. This can be achieved by adding a normalization:

\[r(\theta) = \frac{1}{N(t)}\big(1 + \epsilon (s_1(t) \cos(2\theta) + s_2(t) \sin(3\theta) )\big)\]


\[N(t) = \sqrt {1 + \frac{1}{2}\epsilon^2\left ( {s_1(t)}^2 + {s_2(t)}^2 \right)}\]


shapevariables = (
        lambda t: 1-cos(t),
        lambda t: sin(t),

epsilon = 0.1
R = lambda s,q: 1 + epsilon*(q[0]*cos(2*s) + q[1]*cos(3*s))
N = lambda s,q: sqrt(1+0.5*epsilon**2*(q[0]**2+q[1]**2))
X = lambda s,q: R(s,q)*cos(s)/N(s,q)
Y = lambda s,q: R(s,q)*sin(s)/N(s,q)
bodies = (
Re=0 Re=1000
Re=10 Potential, irrotatoinal

Shapere & Wilczek Amoeba

Here we take a different exmample of an amoeba studied by Shapere and Wilczek in [SW89].

This time the boundary is given by:

\begin{eqnarray*} x(\theta)&=&\cos\theta\\ &&\quad+s_{1}(t)(a\cos\theta+b\sin2\theta)\\ &&\quad+s_{2}(t)(-a\cos2\theta+b\sin\theta)\\ y(\theta)&=&\sin\theta\\ &&\quad+s_{1}(t)(-a\sin\theta+b\cos2\theta)\\ &&\quad+s_{2}(t)(a\sin2\theta+b\cos\theta) \end{eqnarray*}

Below we use the shape modes evolving in time as follows

\begin{eqnarray*} s_{1}(t)&=&\cos t\\ s_{2}(t)&=&\sin t \end{eqnarray*}

and the parameters are set to \(a=\frac{3}{10}\quad b=\frac{3}{200}\)

We will not compute the immersed boundary case here, as ensuring the area remains constant in this case is not trivial.

Re=0 Potential, irrotation


[MB99]Richard Mason and Joel Burdick. Propulsion and control of deformable bodies in an ideal fluid. In Robotics and Automation, 1999. Proceedings. 1999 IEEE International Conference on, volume 1, 773–780. IEEE, 1999.
[SW89]Alfred Shapere and Frank Wilczek. Geometry of self-propulsion at low reynolds number. Journal of Fluid Mechanics, 198:557–585, 1989.

Thu, 30 Oct 2014 00:00:00 -0700 <![CDATA[Centrosymmetric Group]]> Centrosymmetric Group

Consider the set of \(N\times N\) matrices that satisfy the property

\[\mathcal{H}=\{H\,|\,H_{ij}=H_{N+1-i,N+1-j},\det H\neq0\}\]

or in matrix forms

\[\begin{split}\begin{pmatrix}a_{1} & a_{2} & \cdots & a_{N-1} & a_{N}\\ b_{1} & b_{2} & \cdots & b_{N-1} & b_{N}\\ \vdots & \vdots & \ddots & \vdots & \vdots\\ b_{N} & b_{N-1} & \cdots & b_{2} & b_{1}\\ a_{N} & a_{N-1} & \cdots & a_{2} & a_{1} \end{pmatrix}\end{split}\]


\[\begin{split}J=\begin{bmatrix}. & . & \dots & . & 1\\ . & . & \dots & 1 & .\\ \\ . & 1 & \dots & . & .\\ 1 & . & \dots & . & . \end{bmatrix}\qquad J_{ij}=\begin{cases} 1 & \text{if }i=N-j+1\\ 0 & \text{otherwise} \end{cases}\end{split}\]

The matrix \(J\) is sometimes referred to as the exchange matrix.

Using the exchange matrix the above set of matrices can also be conveniently defined as

\[\mathcal{H}=\left\{ A\,|\,AJ=JA\right\}\]

Proposition. The set of matrices \(\mathcal{H}\) forms a group.

Proof. It can be easily shown that the identity is in the group, i.e. \(I\in\mathcal{H}\). Also the set is closed under multiplication, i.e. if \(A\in\mathcal{H}\) and \(B\in\mathcal{H}\) then \(AB\in\mathcal{H}\). To see why, consider


If \(A\in\mathcal{H}\), then \(AJ=JA\)


Therefore, \(A^{-1}\in\mathcal{H}\). This completes the proof that \(\mathcal{H}\) forms a group.

Wed, 27 Jul 2011 00:00:00 -0700 <![CDATA[The Orthogonal Complement of the Space of Row-null and Column-null Matrices]]> The Orthogonal Complement of the Space of Row-null and Column-null Matrices

Lemma. Let \(Z\in\text{GL}(n,R)\), and let \(Y\in\mathcal{S}(n,R)\) where \(\mathcal{S}(n,R)\) is the space of row-null column-null \(n\times n\) matrices. Then \(\text{Tr}(ZY)=0\) if and only if \(Z\) has the form


Proof. Consider the space of row-null and column-null matrices

\[\mathcal{S}(n,R)=\{Y_{ij}\in GL(n,R):\sum_{i}Y_{ij}=0,\sum_{j}Y_{ij}=0\}\]

Its dimension is


since the row-nullness and column-nullness are defined by \(2N\) equations, only \(2N-1\) of which are linearly independent.

Consider the following space

\[\mathcal{G}(n,R)=\{Z_{ij}\in GL(n,R):Z_{ij}=\left(p_{j}-p_{i}\right)+\left(q_{j}+q_{i}\right)\}\]

Its dimension is


where \(N-1\) is the contribution from the antisymmetric part and \(N\) is from the symmetric part.

Assume \(Y\in\mathcal{S}\) and \(Z\in\mathcal{G}\), then the Frobenius inner product of two such elements is

\[\begin{split}\begin{aligned} \text{Tr}(ZY) & =\sum_{ij}\left[\left(p_{j}-p_{i}\right)Y_{ji}+\left(q_{j}+q_{i}\right)Y_{ji}\right]\\ & =\sum_{j}(q_{j}+p_{j})\sum_{i}Y_{ji}+\sum_{i}(q_{i}-p_{i})\sum_{j}Y_{ji}\\ & =0\\\end{aligned}\end{split}\]

Since \(\text{dim}(\mathcal{G})+\text{dim}(\mathcal{S})=\text{dim}(GL)\) and \(\mathcal{G}\perp\mathcal{S}\), then \(\mathcal{G}\) and \(\mathcal{S}\) must be complementary in \(GL\). Therefore, if \(Y\) is orthogonal to all the matrices in \(\mathcal{S}\), it must lie in \(\mathcal{G}\).

Fri, 03 Sep 2010 00:00:00 -0700 <![CDATA[Variational Derivation of the Equation of Motion of an Ideal Fluid]]> Variational Derivation of the Equation of Motion of an Ideal Fluid

Continuous Diffeomorphisms

For ideal fluids, the configuration space is the group of volume preserving diffeomorphisms \(\text{Diff}_{\text{Vol}}(\mathcal{F})\) on the fluid container (a region \(\mathcal{F}(t)\) in \(\mathbb{R}^{2}\) or \(\mathbb{R}^{3}\) that changes with time as a result of the boundary motion). A particle located at a point \(x_{0}\in\mathcal{F}(0)\) will travel to a point \(x_{t}=\varphi_{t}(x_{0})\in\mathcal{F}(t)\) at time \(t\).


The kinetic energy is a mapping from the tangent space to the real numbers

\[L:\quad T\text{Diff}_{\text{Vol}}\to\mathbb{R}\]



Notice that the Lagrangian satisfies the particle relabeling symmetry expressed as invariance under right composition:


where \(\phi,\varphi\in\text{TDiff}_{Vol}\).

The action is given by


Due to the particle relabeling symmetry the system can be described in terms of a reduced Lagrangian

\[l:\quad\mathfrak{diff}_{Vol}\to\mathbb{R}\quad\xi\mapsto l(\xi)=L(e,\xi)=L(\phi,\xi\circ\phi)\]

where \(\mathfrak{\xi\in diff}_{Vol}\) is in the Lie algebra of volume preserving diffeomorphisms, and \(e\) is the identity element of the \(\text{TDiff}_{Vol}\) group.

To obtain the reduced action

\[\begin{split}\begin{aligned} S(\varphi(t),\dot{\varphi}(t)) & =\int_{t_{0}}^{t_{1}}L(\varphi(t),\dot{\varphi}(t))\,dt\\ & =\int_{t_{0}}^{t_{1}}L(e,\dot{\varphi}(t)\circ\varphi(t)^{-1})\,dt\\ & =\int_{t_{0}}^{t_{1}}l(\xi(t))\,dt\quad\text{where }\xi=\dot{\varphi}\circ\varphi^{-1}\\ & =:s(\xi(t))\quad\text{Reduced Action Function }\end{aligned}\end{split}\]

Variations must satisfy the Lin Constraint:

\[\begin{split}\begin{aligned} \delta\xi & =\delta(\dot{\varphi}\varphi^{-1})\\ & =\delta\dot{\varphi}\varphi^{-1}-\dot{\varphi}\varphi^{-1}\delta\varphi\varphi^{-1}\\ & =\dot{\eta}+[\xi,\eta],\quad\text{where }\eta=\delta\varphi\circ\varphi^{-1}\\\end{aligned}\end{split}\]

The variation in the reduced action is

\[\begin{split}\begin{aligned} \delta s & =\delta\left(\frac{1}{2}\int_{t_{1}}^{t_{2}}\int_{\mathcal{F}}v^{2}\:dtdV\right)\\ & =\int_{t_{1}}^{t_{2}}\int_{\mathcal{F}}\,\langle v^{\flat},\delta v\rangle dtdV\\ & =\int_{t_{1}}^{t_{2}}\int_{\mathcal{F}}\,\left(\langle v^{\flat},\dot{u}\rangle+\langle v^{\flat},[v,u]\rangle\right)dtdV\\\end{aligned}\end{split}\]
\[\langle v^{\flat},[v,u]\rangle=\langle v^{\flat},L_{v}u\rangle=L_{v}\langle v^{\flat},u\rangle-\langle L_{v}v^{\flat},u\rangle\]

For a divergence free vector field,

\[\int_{\mathcal{F}}L_{v}f\overset{\nabla\cdot v=0}{=}\int_{\partial\mathcal{F}}(v,n)f=0\]


\[\delta s=\int_{t_{1}}^{t_{2}}\int_{\mathcal{F}}\,\left(\langle v^{\flat},\dot{u}\rangle-\langle L_{v}v^{\flat},u\rangle\right)dtdV\]

Setting the variation \(\delta s=0\) and integrating by parts:


This implies that the integrand must be the gradient of a function


which is equivalent to the Euler equation in coordinate form:

Sun, 28 Jun 2009 00:00:00 -0700 <![CDATA[Coupled Massless Body - Fluid System]]> Coupled Massless Body - Fluid System

For a massless body in 2D the force and momentum balance are:

\[F_{x}=0\quad F_{y}=0\quad M_{z}=0\]

where \(F_{x}\) ,\(F_{y}\) and \(M_{z}\) are the force and momentum exerted on the body by the surrounding fluid.

The configuration of the body is described by the following variables:

\[\left\{ x(t),y(t),\theta(t),s_{1}(t),s_{2}(t),\dots\,\right\}\]

where the first three describe its position and orientation in space that are to be determined, and the rest describe its internal shape configuration that is prescribed a priori.

The problem then is to find such \((\dot{x},\dot{y},\dot{\theta})\) that together with the given set \((\dot{s_{1}}.\dot{s_{2}},\dots)\) would be consistent with the momentum balance at every instance in time during the motion. The translation and the rotation of the body must be such that the immersed boundary code would give zero net force.

But since we do not know \((\dot{x},\dot{y},\dot{\theta})\) a priori, we cannot rely on the immersed boundary code to compute the force. This is why we need to rely on geometric mechanics to find the “mechanical connection”, namely how a trajectory in shape space\((s_{1}(t),s_{2}(t),\dots)\) determines the trajectory in position space \((x(t),y(t),\theta(t))\).

Sat, 20 Jun 2009 00:00:00 -0700 <![CDATA[Squirming Sphere in an Ideal Fluid]]> Squirming Sphere in an Ideal Fluid


The purpose of this work is to derive analytically and with an explicit formula the motion due to the variation of the surface of a squirming body in a fluid. Such a derivation has been done before, but only in two dimensions, and there is no literature on any attempts to generalize to three dimensions, which is considerably more involved algebraically. But utilizing software such as Mathematica, the derivation is quite straightforward.

The primary motivation behind this work is the search for innovative propulsion mechanisms. Surface variations may not lead to a very high speed, but they provide higher maneuverability and efficiency and are the primary source of motion of aquatic animals. Two types of fluid are studied - an Eulerian (potential) fluid where all viscous effects are ignored and inertia is the only factor, and a Stokesian (creeping) fluid where all inertial effects are ignored and viscous friction is the only factor. We consider only those cases because they linearize the equations governing the fluid motion, and permit analytic solutions.

Hydromechanical Connections

Consider a body immersed in a fluid. We regard its shape as a shape manifold \(M\). And the configuration manifold

\[Q=M\times G\]

represents the complete position and orientation of the body in space.

\(Q=M\times G\) can be thought of as a trivial bundle with a group structure G, and an action of \(G\) on \(Q\) given by \(\Phi_{h}(r,g)=(r,hg)\) for \((r,g)\in Q\) and \(h\in G.\) An element of the tangent space is given by \(\dot{q}=(\dot{r},\dot{g})\in TQ\).

We are interested in computing a connection on the configuration bundle \(Q\) that will allow us to compute the fiberwise translations resulting from cyclical changes in shape.

A connection is a vector valued one-form \(\Gamma:TQ\to\mathfrak{g}\) on \(Q\).

We follow [Kel98] in defining two types of connections:

  • Mechanical connection:

    where \(\mathbb{I}:\mathfrak{g}\to\mathfrak{g}^{*}\) is the locked inertia tensor such that
    \[ \begin{align}\begin{aligned}\langle\mathbb{I}(q)\xi,\eta\rangle=\ll\xi_{Q}(q),\eta_{Q}(q)\gg_{KE}\\and\end{aligned}\end{align} \]
    and \(J:TQ\to\mathfrak{g}^{*}\) is the momentum map.
  • Stokes connection:


where \(K:TQ\to\mathfrak{g}^{*}\) is a momentum map

\[\langle K(q,\dot{q}),\xi\rangle=\langle F(q,\dot{q}),\xi_{Q}(q)\rangle\]
associated with a dissipative force \(F\), and \(\mathbb{V}:\mathfrak{g}\to\mathfrak{g}^{*}\) is a viscosity tensor such that

\[ \begin{align}\begin{aligned}\langle\mathbb{V}(q)\xi,\eta\rangle=\langle F(\xi_{Q}(q)),\eta_{Q}(q)\rangle\\and\end{aligned}\end{align} \]

Kelly uses the following local expressions

\[\begin{split}\begin{aligned} \Gamma_{\text{mech}} & = & Ad_{g}(g^{-1}\dot{g}+A_{\text{mech}}\dot{r})\\ \Gamma_{\text{mech}} & = & Ad_{g}(g^{-1}\dot{g}+A_{\text{Stokes}}\dot{r})\end{aligned}\end{split}\]

and for zero initial momentum

\[\begin{split}\begin{aligned} g^{-1}\dot{g} & = & -A_{\text{mech}}\dot{r}\\ g^{-1}\dot{g} & = & A_{\text{Stokes}}\dot{r}\end{aligned}\end{split}\]

which completely describe the evolution of a Potential and Stokesian swimmers respectability.

Squirming Sphere

Consider a three dimensional nearly spherical body whose surface is described by the following equation:


Such a body has already been studied previously numerically [Ruf05], and the current method that allows us to obtain explicit formulas for the displacement as will be later seen, are valuable tools for validating numerical simulations.

Coordinate system.

Coordinate system.

As an example we consider the propulsive movement of a roughly spherical device whose boundary shape is modulated by a small amount. We parametrize the boundary in terms of the shape variables \(s_{i}(t)\), and express it as a sum of an unperturbed base shape and a perturbed part:


The configuration manifold is

\[Q=\mathbb{R}^{2}\times SE(3)\]
Boundary perturbations.

Boundary perturbations.

Potential Flow

Irrotational potential flow is best described in terms of the velocity potential, which satisfies


with the Neumann boundary condition


where \(\boldsymbol{u}\big|_{r=F}\) is the prescribed boundary motion.

For an axisymmetric problem, such as this one, the solution of the Laplace equation can be expressed as a series

\[ \begin{align}\begin{aligned}\phi(r,\theta)=\sum_{k=0}^{\infty}\left(A_{k}r^{k}+\frac{B_{k}}{r^{k+1}}\right)P_{k}(\cos\theta)\\where :math:`P_{k}(x)` are the Legendre polynomials satisfying the\end{aligned}\end{align} \]

orthogonality condition:


We require that \(\boldsymbol{u}\to0\) as \(r\to\infty\), and therefore we must take \(A_{k}=0\) for \(\forall k\).


The problem of solving for \(\phi\) is reduced to finding the coefficients \(a_{k}\).

First, in order to enforce the boundary conditions on the perturbed surface, we need to compute the normal \(\boldsymbol{n}\) at each point. The two tangential components on the surface are

\[\begin{split}\begin{aligned} \boldsymbol{t_{1}} & =\frac{\partial}{\partial\theta}(F\boldsymbol{\hat{r}})=\frac{\partial F}{\partial\theta}\hat{\boldsymbol{r}}+F\boldsymbol{\hat{\theta}}\\ \boldsymbol{t_{2}} & =\frac{\partial}{\partial\varphi}(F\boldsymbol{\hat{r}})=\frac{\partial F}{\partial\varphi}\boldsymbol{\hat{r}}+F\sin\theta\,\boldsymbol{\hat{\varphi}}\end{aligned}\end{split}\]

and therefore the normal component can be evaluated as the cross product of the two

\[ \begin{align}\begin{aligned}\begin{split} \begin{aligned} \boldsymbol{n} & =\boldsymbol{t_{1}}\times\boldsymbol{t_{2}}\\ & =F^{2}\sin\theta\:\boldsymbol{\hat{r}}-F\frac{\partial F}{\partial\theta}\sin\theta\,\boldsymbol{\hat{\theta}}-F\frac{\partial F}{\partial\varphi}\,\boldsymbol{\hat{\varphi}}\end{aligned}\end{split}\\Note that we need not normalize the expression for\end{aligned}\end{align} \]

\(\boldsymbol{n}\) since it occurs on both sides of the boundary condition equation Eq.\ref{eq:boundary}. Expanding the spherical coordinates \((r,\theta,\varphi)\) of \(\boldsymbol{n}\) in terms of \(\epsilon\):

\[\begin{split}\begin{aligned} \boldsymbol{n} & =\mathbf{n}^{(0)}+\epsilon\mathbf{n}^{(1)}+\epsilon\mathbf{n}^{(2)}\\ & =\left(\begin{array}{c} 1\\ 0\\ 0 \end{array}\right)+\epsilon\left(\begin{array}{c} 2\left(c[2\theta]s_{1}+c[3\theta]s_{2}\right)\\ 2s[2\theta]s_{1}+3s[3\theta]s_{2}\\ 0 \end{array}\right) \\ & \quad +\epsilon^{2}\left(\begin{array}{c} \left(c[2\theta]s_{1}+c[3\theta]s_{2}\right){}^{2}\\ \left(c[2\theta]s_{1}+c[3\theta]s_{2}\right)\left(2s[2\theta]s_{1}+3s[3\theta]s_{2}\right)\\ 0 \end{array}\right)\end{aligned}\end{split}\]

where we have used the abbreviations \(c[\theta]=\cos(\theta)\) and \(s[\theta]=\sin(\theta)\).

Next, we expand the potential in a power series in terms of the small perturbation parameter \(\epsilon\) .

\[ \begin{align}\begin{aligned}\phi=\phi^{(0)}+\epsilon\phi^{(1)}+\epsilon^{2}\phi^{(2)}+...+\epsilon^{i}\phi^{(i)}+...\\where\end{aligned}\end{align} \]
\[ \begin{align}\begin{aligned}\phi^{(i)}=\sum_{k=0}^{\infty}a_{k}^{\phantom{k}(i)}\Phi_{k}\\\begin{split}Plugging back in Eq.\ `\\ref{eq:boundary} <#eq:boundary>`__, we expand both\end{split}\end{aligned}\end{align} \]

sides in terms of \(\epsilon\)

\[\begin{split}\begin{aligned} \nabla\phi\cdot\boldsymbol{n}\big|_{r=F} & =\sum_{i=0}^{\infty}\epsilon^{i}\sum_{k=0}^{\infty}a_{k}^{(i)}f_{k}^{(i)}\\ \boldsymbol{u}\cdot\boldsymbol{n}\big|_{r=F} & =\sum_{i=0}^{\infty}\epsilon^{i}u^{(i)}\end{aligned}\end{split}\]

and equate terms of equal power


where the first few terms \(f_{k}^{(0)}\) are

\[\begin{split}\begin{aligned} f_{0}^{(0)} & =-1\\ f_{1}^{(0)} & =-2\eta\\ f_{2}^{(0)} & =\frac{3}{2}\left(1-3\eta^{2}\right)\\ f_{3}^{(0)} & =2\left(3\eta-5\eta^{3}\right)\\ f_{4}^{(0)} & =-\frac{5}{8}\left(3-30\eta^{2}+35\eta^{4}\right)\\ f_{5}^{(0)} & =-\frac{3}{4}\left(15\eta-70\eta^{3}+63\eta^{5}\right)\\ f_{6}^{(0)} & =\frac{7}{16}\left(5-105\eta^{2}+315\eta^{4}-231\eta^{6}\right)\\ f_{7}^{(0)} & =\frac{1}{2}\left(35\eta-315\eta^{3}+693\eta^{5}-429\eta^{7}\right)\end{aligned}\end{split}\]

where \(\eta=\cos\theta\).

By the linearity of Laplace’s equation (using superposition) one can write following Kirchoff,


Then, up to second order in \(\epsilon\)

\[\begin{split}\begin{aligned} \phi_{z} & =-\frac{P_{1}}{2r^{2}}+\epsilon\left(\frac{9P_{1}s_{1}}{10r^{2}}+\frac{9P_{2}s_{2}}{7r^{3}}-\frac{6P_{3}s_{1}}{5r^{4}}-\frac{48P_{4}s_{2}}{35r^{5}}\right)\\ & \quad +\epsilon^{2}\bigl(\frac{P_{1}\left(-3339s_{1}^{2}-4990s_{2}^{2}\right)}{2450r^{2}}-\frac{1781P_{2}s_{1}s_{2}}{735r^{3}}+\frac{6P_{3}\left(308s_{1}^{2}-45s_{2}^{2}\right)}{1925r^{4}}\\ & \qquad \qquad +\frac{42264P_{4}s_{1}s_{2}}{13475r^{5}}-\frac{16P_{5}\left(91s_{1}^{2}-90s_{2}^{2}\right)}{637r^{6}}-\frac{64P_{6}s_{1}s_{2}}{11r^{7}}-\frac{512P_{7}s_{2}^{2}}{143r^{8}}\bigg)\\ \\ \phi_{s1} & =\epsilon\left(\frac{P_{0}}{3r}-\frac{4P_{2}}{9r^{3}}\right) + \epsilon^{2}\bigg(-\frac{14P_{0}s_{1}}{15r}-\frac{3P_{1}s_{2}}{35r^{2}} + \frac{8P_{2}s_{1}}{21r^{3}}+\frac{2P_{3}s_{2}}{5r^{4}} \\ & \qquad \qquad \qquad \qquad \qquad \qquad -\frac{512P_{4}s_{1}}{525r^{5}}-\frac{64P_{5}s_{2}}{63r^{6}}\bigg)\\ \\ \phi_{s2} & =\epsilon\left(\frac{3P_{1}}{10r^{2}}-\frac{2P_{3}}{5r^{4}}\right)+\epsilon^{2}\bigg(-\frac{34P_{0}s_{2}}{35r}-\frac{309P_{1}s_{1}}{350r^{2}}-\frac{13P_{2}s_{2}}{315r^{3}}\\ & \qquad \qquad \qquad \qquad \qquad \qquad +\frac{122P_{3}s_{1}}{225r^{4}}+\frac{1128P_{4}s_{2}}{1925r^{5}}-\frac{208P_{5}s_{1}}{189r^{6}}-\frac{1856P_{6}s_{2}}{1617r^{7}}\bigg)\end{aligned}\end{split}\]

The total kinetic energy is given by

\[\begin{split}\begin{aligned} L & =\frac{1}{2}\int_{\mathfrak{F}}dV\,|\nabla\phi|^{2}\\ & \frac{1}{2}\int_{0}^{2\pi}d\phi\,\int_{0}^{\pi}\sin\theta d\theta\int_{r=F}^{\infty}r^{2}dr\:|\nabla\phi|^{2}\\ & =\pi\int_{-1}^{1}d\eta\int_{r=F}^{\infty}r^{2}dr\:|\nabla\phi|^{2}\\ & =\frac{\pi}{5}\dot{z}^{2}-\epsilon\frac{\pi}{5}\text{ }\left(\frac{49}{15}\text{ }s_{1}\dot{z}^{2}+\frac{6}{5}\dot{s}_{2}\dot{z}\right)\\ & \quad +\epsilon^{2}\frac{\pi}{5}\bigg(\frac{230}{189}\dot{s}_{1}^{2}+\frac{241}{225}\dot{s}_{2}^{2}+2\left(\frac{1069}{525}s_{1}\dot{s}_{2}-\frac{628}{245}s_{2}\dot{s}_{1}\right)\dot{z}+\\ & \qquad \qquad + \left(\frac{576}{175}s_{1}^{2}+\frac{144857}{18865}s_{2}^{2}\right)\dot{z}^{2}\bigg)\\ & \quad +O(\epsilon^{3})\end{aligned}\end{split}\]

This Lagrangian is invariant under translation along \(z\), so its momentum map is

\[\begin{split}\begin{aligned} J & = \frac{2\pi}{5}\text{ }\dot{z}-\epsilon\frac{2\pi}{5}\left(\frac{49}{15}\text{ }s_{1}\dot{z}+\frac{3}{5}\dot{s}_{2}\right) \\ & \quad +\epsilon^{2}\frac{2\pi}{5}\left(\frac{1069}{525}s_{1}\dot{s}_{2}-\frac{628}{245}s_{2}\dot{s}_{1}+\left(\frac{576}{175}s_{1}^{2}+\frac{144857}{18865}s_{2}^{2}\right)\dot{z}\right)+O(\epsilon^{3})\end{aligned}\end{split}\]

The locked inertia tensor is given by

\[\mathbb{I}(s_{1},s_{2}):\quad\xi\mapsto\frac{2\pi}{5}\text{ }\xi\left(1-\epsilon\frac{49}{15}\text{ }s_{1}+\epsilon^{2}\left(\frac{576}{175}s_{1}^{2}+\frac{144857}{18865}s_{2}^{2}\right)\right)+O(\epsilon^{3})\]


\[\begin{split}\begin{aligned} \Gamma_{\text{mech}} & =\mathbb{I}^{-1}J \\ & =dz-\epsilon\frac{3}{5}ds_{2}+\epsilon^{2}\left(\frac{8}{105}s_{1}\dot{s}_{2}-\frac{628}{245}s_{2}\dot{s}_{1}\right)+O(\epsilon^{3})\end{aligned}\end{split}\]

Self-propulsion of the sphere requires that its motion remain horizontal with respect to this connection, therefore


Stokes Flow

The solution of creeping-motion can be reduce in two dimensions to the biharmonic equation in terms of the stream function


which satisfies the no-slip and no-penetration boundary condition on the surface


Following [Lea07] p.462 the general solution is


where \(Q_{k}(x)\) are polynomial functions defined in terms of the Legendre polynomials


and satisfying the orthogonality condition


Again we want the velocity to vanish at infinity, so we will neglect the first two terms

\[u_{r}=-\frac{1}{r^{2}}\frac{\partial\psi}{\partial\eta}\qquad u_{\theta}\sqrt{1-\eta^{2}}=-\frac{1}{r}\frac{\partial\psi}{\partial r}\]

Again we will exploit the linearity of the equation by using Kirchoff’s method to express the streamfunction as

\[\begin{split}\begin{aligned} \psi_{z} & =\left(\frac{1}{2r}-\frac{3r}{2}\right)Q_{1}+\\ & +\epsilon\bigg[Q_{3}\left(\frac{12s_{1}}{5r^{3}}-\frac{12s_{1}}{5r}\right)+Q_{1}\left(-\frac{9s_{1}}{10r}+\frac{9rs_{1}}{10}\right)\\ & \qquad \qquad + Q_{4}\left(\frac{24s_{2}}{7r^{4}}-\frac{24s_{2}}{7r^{2}}\right)+Q_{2}\left(\frac{27s_{2}}{14}-\frac{27s_{2}}{14r^{2}}\right)\bigg]\\ & +\epsilon^{2}\bigg[Q_{1}\left(-\frac{3r\left(448s_{1}^{2}+1045s_{2}^{2}\right)}{1225}+\frac{4683s_{1}^{2}+8125s_{2}^{2}}{2450r}\right)\\ & \qquad \qquad +Q_{2}\left(-\frac{1287}{490}s_{1}s_{2}+\frac{2197s_{1}s_{2}}{490r^{2}}\right)\\ & \qquad \qquad +Q_{3}\left(-\frac{6\left(308s_{1}^{2}-135s_{2}^{2}\right)}{1925r}-\frac{4\left(308s_{1}^{2}+45s_{2}^{2}\right)}{1925r^{3}}\right)\\ & \qquad \qquad +Q_{4}\left(-\frac{10764s_{1}s_{2}}{2695r^{4}}-\frac{1836s_{1}s_{2}}{2695r^{2}}\right)+Q_{7}\left(\frac{2048s_{2}^{2}}{143r^{7}}-\frac{1536s_{2}^{2}}{143r^{5}}\right)\\ & \qquad \qquad +Q_{5}\left(\frac{48\left(91s_{1}^{2}-54s_{2}^{2}\right)}{637r^{5}}-\frac{32\left(91s_{1}^{2}-18s_{2}^{2}\right)}{637r^{3}}\right)\\ & \qquad \qquad +Q_{6}\left(\frac{224s_{1}s_{2}}{11r^{6}}-\frac{160s_{1}s_{2}}{11r^{4}}\right)\bigg]\end{aligned}\end{split}\]

Due to the lack of time, I am unable to continue with the calculations. The Mathematica code that I use seems to have trouble computing the \(\psi_{s1}\) and \(\psi_{s2}\) components. This could be related the 3D equivalent of the Stoke’s paradox, or could simply be a bug in the code. In any case, it would be important to carry this calculation to completion at some point in the future in order to compare the displacement resulting from creeping flow to that in potential flow.


A calculation similar to the one done above was carried by Kelly for a two dimensional circle where the algebra is less involved. Figure [Flo:comp] compares the results. In all cases the shape trajectory is a unit circle \((s_{1},s_{2})\) shape space.

Geometric phase for a particular gate. TODO: stokes 3d.

Geometric phase for a particular gate. TODO: stokes 3d.

Explicitly, we quote the results. For the the potential flow case, Kelly obtained:


and for Stoke’s case


If we compare those, to our result


The first thing that we notice is the presence of \(O(\epsilon^{1})\) in the 3d case (our solution). However, this term has no effect on the net displacement and it only results in a more pronounced oscillation of the center of mass during each cycle. Also we notice that the body achieves higher displacement in 3d.

This simplicity with which such explicit equations such as Eq.\ref{eq:magic} can be obtained is a testament to the power of the geometric approach to mechanics.


[Blo07]Anthony Bloch. Nonholonomic Mechanics and Control. Springer, September 2007. ISBN 0387955356.
[Hol08]Darryl D. Holm. GEOMETRIC MECHANICS: Dynamics and Symmetry. Imperial College Press, illustrated edition edition, May 2008. ISBN 1848161956.
[Kel98]Scott D Kelly. The Mechanics and Control of Robotic Locomotion with Appliations to Aquatic Vehicles. PhD thesis, California Institute of Technology, May 1998.
[KP04]M. Kohr and I. Pop. Viscous Incompressible Flow for Low Reynolds Numbers (Advances in Boundary Elements). WIT Press (UK), illustrated edition edition, July 2004. ISBN 1853129917.
[Lea07]L. Gary Leal. Advanced Transport Phenomena: Fluid Mechanics and Convective Transport Processes. Cambridge University Press, 1 edition, June 2007. ISBN 0521849101.
[MR02]Jerrold E. Marsden and Tudor S. Ratiu. Introduction to Mechanics and Symmetry: A Basic Exposition of Classical Mechanical Systems. Springer, 2nd edition, December 2002. ISBN 038798643X.
[MB99]R. Mason and J. Burdick. Propulsion and control of deformable bodies in an ideal fluid. In 1999 IEEE International Conference on Robotics and Automation, 1999. Proceedings, volume 1. 1999.
[MRR06]J. B. Melli, C. W. Rowley, and D. S. Rufat. Motion planning for an articulated body in a perfect planar fluid. SIAM Journal on Applied Dynamical Systems, 5:650, 2006.
[NK07]S. Nair and E. Kanso. Hydrodynamically coupled rigid bodies. Journal of Fluid Mechanics, 592:393–411, 2007.
[Ruf05]Dzhelil S Rufat. Locomotion of a three dimensional variable shape body in inviscid fluid. May 2005.

Mon, 08 Jun 2009 00:00:00 -0700