Onur Solmaz · Post · /2017/12/22

Nonlinear Finite Elements

This post builds on the formulations I showed in my previous posts by introducing their nonlinear versions.

In a typical nonlinear problem, the variational setting leads to the weak formulation

Find $u\in V$ such that
$\begin{equation} F(u,v) = 0 \htmlId{eq:femnonlinear2}{} \tag{1}\end{equation}$
for all $v\in V$ where the semilinear form $F$ is nonlinear in terms of $u$ and linear in terms of $v$ .

We linearize $F$ :

\begin{equation} \Lin [F(u,v)]_{u=\bar{u}} = F(\bar{u}, v) + \varn{F(u,v)}{u}{\Var u}\evat_{u=\bar{u}} \htmlId{eq:femnonlinear4}{} \tag{2}\end{equation}

Equating (2) to zero yields a linear system in terms of $\Var u$

\begin{equation} \boxed{ a(\Var u, v) = b(v) } \htmlId{eq:femnonlinear6}{} \tag{3}\end{equation}

where

\begin{equation} \boxed{ \begin{aligned} a(\Var u, v) &= \varn{F(u,v)}{u}{\Var u}\evat_{u=\bar{u}} \\ b(v) &= -F(\bar{u}, v) . \end{aligned} } \htmlId{eq:femnonlinear5}{} \tag{4}\end{equation}

We can compute the components of the matrices and vectors according to (3)

\begin{equation} \boxed{ \begin{alignedat}{3} \Aelid{I\!J}{} &= a(N^J,N^I) &&= \varn{F(u,N^I)}{u}{N^J}\evat_{u=\bar{u}} \\ b^{I} &= b(N^I) &&= -F(\bar{u}, N^I). \end{alignedat} } \htmlId{eq:femnonlinear9}{} \tag{5}\end{equation}

Then the update vector $\Var \Bu = [\Var u^1, \Var u^2, \dots, \Var u^{\nnode}]\tra$ is obtained by solving

\begin{equation} \BA \Var \Bu = \Bb \tag{6}\end{equation}

Letting $\Var u$ be the difference between consequent iterates, we obtain the update equation as

\begin{equation} \boxed{ \Bu \leftarrow \bar{\Bu} + \Var\Bu } \tag{7}\end{equation}

Example: Nonlinear Poisson’s Equation

Consider the following nonlinear Poisson’s equation

\begin{equation} \begin{alignedat}{4} - \nabla \dtp (g(u)\nabla u) &= f \quad && \text{in} \quad && \Omega \\ u &= 0 \quad && \text{on} \quad && \del\Omega \end{alignedat} \htmlId{eq:femnonlinear8}{} \tag{8}\end{equation}

The weak formulation reads

Find $u\in V$ such that
$\begin{equation} - \int_\Omega \nabla \dtp (g(u)\nabla u) v \dv= \int_\Omega f v \dv \tag{9}\end{equation}$
for all $v\in V$ where $V=H^1_0(\Omega)$ .

Applying integration by parts and divergence theorem on the left-hand side

\begin{equation} \begin{aligned} \int_\Omega \nabla \dtp (g(u)\nabla u) v \dv &= \int_\Omega \nabla \dtp (g(u)\nabla (u) v) \dv - \int_\Omega g(u)\nabla u\dtp\nabla v \dv \\ &= \underbrace{\int_{\del\Omega} g(u) v (\nabla u\dtp\Bn) \da}_{v = 0 \text{ on } \del\Omega} - \int_\Omega g(u)\nabla u\dtp\nabla v \dv \\ \end{aligned} \tag{10}\end{equation}

Thus we have the semilinear form

\begin{equation} F(u,v) = \int_{\Omega} g(u) \nabla u \dtp \nabla v \dv - \int_{\Omega} f \, v \dv = 0 \tag{11}\end{equation}

The linearized version of this problem is then with (4)

\begin{equation} \begin{aligned} a(\Var u,v) &= \int_{\Omega} \rbr{\deriv{g}{u}\evat_{\bar{u}} \Var u\, \nabla \bar{u} + g(\bar{u})\nabla(\Var u)} \dtp \nabla v \dv \\ b(v) &= \int_{\Omega} [f \, v - g(\bar{u}) \nabla \bar{u} \dtp \nabla v] \dv\\ \end{aligned} \tag{12}\end{equation}

and the matrix and vector components are with (5)

\begin{equation} \begin{aligned} \Aelid{I\!J}{} &= \int_{\Omega} \rbr{\deriv{g}{u}\evat_{\bar{u}} N^J \, \nabla \bar{u} + g(\bar{u})\BB^J} \dtp \BB^I \dv \\ b^{I} &= \int_{\Omega} [f \, N^I - g(\bar{u}) \nabla \bar{u} \dtp \BB^I ] \dv\\ \end{aligned} \tag{13}\end{equation}

where the previous solution and its gradient are computed as

\begin{equation} \bar{u} = \suml{I=1}{\nnode} \bar{u}^I N^I \eqand \nabla \bar{u} = \suml{I=1}{\nnode} \bar{u}^I \BB^I . \tag{14}\end{equation}

Nonlinear Time-Dependent Problems

In the case of a nonlinear time-dependent problem, we have the following weak form:

Find $u \in V$ such that
$\begin{equation} m(\dot{u}, v; t) + F(u,v; t) = 0 \htmlId{eq:nonlineartimedependentweak1}{} \tag{15}\end{equation}$
for all $v \in V$ and $t \in [0,\infty)$ where $F$ is a semilinear form.

Discretization yields the following nonlinear system of equations

\begin{equation} \BM(t)\Bu + \Bf(u; t) = \Bzero \tag{16}\end{equation}

where

\begin{equation} \begin{aligned} M^{I\!J}(t) &= m(N^J, N^I; t) \\ f^{I}(u;t) &= F(u, N^I; t). \end{aligned} \tag{17}\end{equation}

Explicit Euler Scheme

We discretize in time with the finite difference $\dot{u} \approx [u_{n+1}-u_n]/{\Delta t}$ and linearity allows us to write

\begin{equation} \boxed{ m(\dot{u}, v; t) \approx \frac{1}{\Delta t} [m(u_{n+1}, v; t_{n+1}) - m(u_n, v; t_n)] } \label{eq:discretetimedependent1} \end{equation}

We discretize the variational forms in time according to \eqref{eq:discretetimedependent1}, and evaluate the remaining terms at $t_n$ :

\begin{equation} \frac{1}{\Delta t} [m(u_{n+1},v;t_{n+1}) - m(u_{n},v;t_{n})] + F(u_n, v; t_n) = 0 \tag{18}\end{equation}

The corresponding system of equations is

\begin{equation} \frac{1}{\Delta t} [\BM_{n+1}\Bu_{n+1} - \BM_n\Bu_n] + \Bf_n = \Bzero \tag{19}\end{equation}

where $\Bf_n = \Bf(u_n, t_n)$ . This yields the following update equation

\begin{equation} \boxed{ \Bu_{n+1} = \BM_{n+1}\inv [\BM_n\Bu_n - \Delta t \Bf_n] } \tag{20}\end{equation}

For a time-independent $m$ , this becomes

\begin{equation} \Bu_{n+1} = \Bu_n - \Delta t \BM\inv\Bf_n \tag{21}\end{equation}

Implicit Euler Scheme

For the implicit scheme, we evaluate the remaining terms at $t_{n+1}$ and let the result be equal to

\begin{equation} G(u_{n+1}, v) := \frac{1}{\Delta t} [m(u_{n+1},v;t_{n+1}) - m(u_{n},v;t_{n})] + F(u_{n+1}, v; t_{n+1}) = 0 \tag{22}\end{equation}

We will hereon replace $u_{n+1}$ with $u$ for brevity. The update of this nonlinear system requires the linearization of $G(u, v)$ :

\begin{equation} \Lin[G(u,v)]_{u=\bar{u}} = G(\bar{u}, v) + \varn{G}{u}{\Var u}\evat_{u=\bar{u}} = 0 \tag{23}\end{equation}

We thus have the following linear setting for the Newton update $\Var u$ :

\begin{equation} a(\Var u, v) = b(v) \tag{24}\end{equation}

where

\begin{equation} \begin{aligned} a(\Var u, v) &:= \varn{G}{u}{\Var u} \evat_{u=\bar{u}} = \frac{1}{\Delta t} m(\Var u, v; t_{n+1}) + \varn{F(u, v; t_{n+1})}{u}{\Var u} \evat_{u=\bar{u}} \\ b(v) &:= -G(\bar{u}, v) = - F(\bar{u}, v; t_{n+1}) -\frac{1}{\Delta t} [m(\bar{u},v;t_{n+1}) - m(u_{n},v;t_{n})] \end{aligned} \tag{25}\end{equation}

Discretization yields

\begin{equation} \rbr{\frac{1}{\Delta t} \BM_{n+1} + \tilde{\BA}}\Var \Bu = \Bb \tag{26}\end{equation}

where

\begin{equation} \tilde{A}^{I\!J} := \varn{F(u, N^I;t_{n+1})}{u}{N^J} \evat_{u=\bar{u}} \eqand b^I := b(N^I) \tag{27}\end{equation}

The Newton update is rendered

\begin{equation} \boxed{ \Bu \leftarrow \bar{\Bu} + \Var\Bu \eqwith \Var \Bu = [\frac{1}{\Delta t} \BM_{n+1} + \tilde{\BA}]\inv\Bb } \tag{28}\end{equation}

which is repeated until the solution for the next timestep $\Bu$ converges to a satisfactory value.

Nonlinear Coupled Problems

For a nonlinear coupled problem, the weak formulation is as follows

Find $u\in V_1$ , $y\in V_2$ such that
$\begin{equation} \begin{aligned} F(u, y, v) &= 0 \\ G(u, y, w) &= 0 \\ \end{aligned} \htmlId{eq:nonlinearcoupled1}{} \tag{29}\end{equation}$
for all $v\in V_1$ , $w \in V_2$ where $F(\cdot,\cdot, \cdot)$ , $G(\cdot, \cdot, \cdot)$ are nonlinear in terms of $u$ and $y$ and linear in terms of $v$ and $w$ .

We linearize the semilinear forms about the nonlinear terms:

\begin{equation} \begin{alignedat}{4} \Lin[F(u, y, v)]_{\bar{u},\bar{y}} &= F(\bar{u},\bar{y},v) &&+ \varn{F(u, y, v)}{u}{\Var u} \evat_{\bar{u},\bar{y}} &&+ \varn{F(u, y, v)}{y}{\Var y} \evat_{\bar{u},\bar{y}} \\ \Lin[G(u, y, w)]_{\bar{u},\bar{y}} &= G(\bar{u},\bar{y},w) &&+ \varn{G(u, y, w)}{u}{\Var u} \evat_{\bar{u},\bar{y}} &&+ \varn{G(u, y, w)}{y}{\Var y} \evat_{\bar{u},\bar{y}} \end{alignedat} \htmlId{eq:nonlinearcoupled2}{} \tag{30}\end{equation}

where the evaluations take place at $u=\bar{u}$ and $y=\bar{y}$ .

Equating the linearized residuals to zero, we obtain a linear system of the form

\begin{equation} \begin{alignedat}{3} a(u, v) &+ b(y, v) &&= c(v) \\ d(u, w) &+ e(y, w) &&= f(w) \\ \end{alignedat} \htmlId{eq:coupledweakform1}{} \tag{31}\end{equation}

with the bilinear forms $a$ , $b$ , $d$ , $e$ and the linear forms $c$ , $f$ which are defined as

\begin{equation} \begin{gathered} \begin{alignedat}{4} a(\Var u, v) &:= \varn{F(u, y, v)}{u}{\Var u} \evat_{\bar{u},\bar{y}} \quad & b(\Var y, v) &:= \varn{F(u, y, v)}{y}{\Var y} \evat_{\bar{u},\bar{y}} \\ d(\Var u, w) &:= \varn{G(u, y, w)}{u}{\Var u} \evat_{\bar{u},\bar{y}} \quad & e(\Var y, w) &:= \varn{G(u, y, w)}{y}{\Var y} \evat_{\bar{u},\bar{y}} \end{alignedat} \\ \text{and} \\ \begin{aligned} c(v) &:= -F(\bar{u},\bar{y},v) \\ f(w) &:= -G(\bar{u}, \bar{y}, w) \end{aligned} \end{gathered} \tag{32}\end{equation}

Discretizing as done in the previous section, we obtain the following linear system of equations

\begin{equation} \begin{bmatrix} \BA & \BB \\ \BD & \BE \end{bmatrix} \begin{bmatrix} \Var \Bu \\ \Var \By \end{bmatrix} = \begin{bmatrix} \Bc \\ \Bf \end{bmatrix} \tag{33}\end{equation}

whose solution yields the update values $\Var \Bu$ and $\Var \By$ . Thus the Newton update equations are

\begin{equation} \begin{alignedat}{3} \Bu &\leftarrow \bar{\Bu} &&+ \Var\Bu \\ \By &\leftarrow \bar{\By} &&+ \Var\By . \end{alignedat} \tag{34}\end{equation}

Example: Cahn-Hilliard Equation

The Cahn-Hilliard equation describes the process of phase separation, by which the two components of a binary fluid spontaneously separate and form domains pure in each component. The problem is nonlinear, coupled and time-dependent. The IBVP reads

\begin{equation} \begin{alignedat}{4} \partd{c}{t} &= \nabla\dtp(\BM\nabla \mu) \qquad&& \text{in} \qquad&& \Omega\times I \\ \nabla c\dtp\Bn &= 0 && \text{on} && \del\Omega\times I\\ \nabla \mu\dtp\Bn &= 0 && \text{on} && \del\Omega\times I\\ c &= c_0 && \text{in} && \Omega, t = 0 \\ \mu &= 0 && \text{in} && \Omega, t = 0 \\ \end{alignedat} \htmlId{eq:cahnhilliard1}{} \tag{35}\end{equation}

where

\begin{equation} \mu = \deriv{f}{c} - \nabla\dtp(\BLambda\nabla c) \htmlId{eq:cahnhilliard2}{} \tag{36}\end{equation}

and $t\in I = [0,\infty)$ . Here,

$c$ is the scalar variable for concentration,
$\mu$ is the scalar variable for the chemical potential,
$f: c \mapsto f(c)$ is the function representing chemical free energy,
$\BM$ is a second-order tensor describing the mobility of the chemical,
$\BLambda$ is a second-order tensor describing both the interface thickness and direction of phase transition.

The fourth-order PDE governing the problem can be formulated as a coupled system of two second-order PDEs with the variables $c$ and $\mu$ , as demonstrated in (35) and (36).

The weak formulation then reads

Find $c \in V_1$ , $\mu\in V_2$ such that
$\begin{equation} \begin{aligned} \int_\Omega \partd{c}{t} v \dx - \int_\Omega \nabla\dtp(\BM\nabla \mu) v \dx &=0 \\ \int_\Omega \sbr{\mu - \deriv{f}{c}} w \dx + \int_\Omega \nabla\dtp(\BLambda\nabla c) w\dx &= 0 \end{aligned} \tag{37}\end{equation}$
for all $v \in V_1$ , $w \in V_2$ and $t \in I$ .

We discretize in time implicitly with $\del c/\del t \approx (c_{n+1}-c_n)/\Var t$ . We also denote the values for the next timestep $c_{n+1}$ and $\mu_{n+1}$ as $c$ and $\mu$ for brevity. Using integration-by-parts, the divergence theorem, and the given boundary conditions, we arrive at the following nonlinear forms

\begin{equation} \begin{alignedat}{3} F(c,\mu,v) &= \int_\Omega \frac{1}{\Var t} (c-c_n) v \dx + \int_\Omega (\BM\nabla \mu)\dtp \nabla v \dx &&= 0 \\ G(c,\mu,w) &= \int_\Omega \sbr{\mu - \deriv{f}{c}} w \dx - \int_\Omega (\BLambda\nabla c)\dtp \nabla w\dx &&= 0 \end{alignedat} \tag{38}\end{equation}

which is a nonlinear coupled system of the form (29).

We linearize the forms according to (30) and obtain the following variations

\begin{align*} \varn{F}{c}{\Var c} &= \int_\Omega \frac{1}{\Var t} \Var c\, v \dx \\ \varn{F}{\mu}{\Var \mu} &= \int_\Omega (\BM\nabla (\Var\mu))\dtp \nabla v \dx \\ \varn{G}{c}{\Var c} &= - \int_\Omega \dderiv{f}{c}\Var c \, w \dx - \int_\Omega (\BLambda\nabla (\Var c))\dtp \nabla w\dx \\ \varn{G}{\mu}{\Var \mu} &= \int_\Omega \Var\mu \, w \dx \end{align*}

We substitute basis functions and obtain our system matrix and vectors

\begin{align*} P^{I\!J} &= \int_\Omega \frac{1}{\Var t} N^JN^I \dx \\ Q^{IL} &= \int_\Omega (\BM\BB^L)\dtp\BB^I \dx \\ r^{I} &= \int_\Omega \frac{1}{\Var t}(\bar{c}-c_n)N^I \dx + \int_\Omega (\BM\nabla\bar{\mu})\dtp\BB^I\dx \\ S^{K\!J} &= - \int_\Omega \dderiv{f}{c}\evat_{c=\bar{c}} N^J N^K \dx - \int_\Omega (\BLambda \BB^J)\dtp \BB^K\dx \\ T^{K\!L} &= \int_\Omega N^L N^K \dx \\ u^{K} &= \int_\Omega \sbr{\bar{\mu} - \deriv{f}{c}\evat_{c=\bar{c}}} N^K \dx - \int_\Omega (\BLambda\nabla \bar{c})\dtp \BB^K\dx \end{align*}

which constitute the system

\begin{equation} \begin{bmatrix} \BP & \BQ \\ \BS & \BT \end{bmatrix} \begin{bmatrix} \Var \Bc \\ \Var \Bmu \end{bmatrix} = \begin{bmatrix} \Br \\ \Bu \end{bmatrix} \tag{39}\end{equation}

Solution yields the update values $\Var \Bc$ and $\Var \Bmu$ . The Newton update equations are then

\begin{equation} \begin{alignedat}{3} \Bc &\leftarrow \bar{\Bc} &&+ \Var\Bc \\ \Bmu &\leftarrow \bar{\Bmu} &&+ \Var\Bmu . \end{alignedat} \tag{40}\end{equation}

The system is solved for $c_{n+1}$ and $\mu_{n+1}$ at each $t=t_n$ to obtain the evolutions of the concentration and chemical potential.

Onur Solmaz · Post · /2017/12/12

Coupled Finite Elements

In this post, I’ll introduce the FE formulation of a generalized linear and coupled weak form. Said weak formulation has the form

Find $u\in V_1$ , $y\in V_2$ such that
$\begin{equation} \begin{alignedat}{3} a(u, v) &+ b(y, v) &&= c(v) \\ d(u, w) &+ e(y, w) &&= f(w) \\ \end{alignedat} \htmlId{eq:coupledweakform1}{} \tag{1}\end{equation}$
for all $v\in V_1$ , $w \in V_2$ where $a(\cdot, \cdot): V_1\times V_1 \to \IR$ , $b(\cdot, \cdot): V_2\times V_1 \to \IR$ , $d(\cdot, \cdot): V_1\times V_2 \to \IR$ , $e(\cdot, \cdot): V_2\times V_2 \to \IR$ are bilinear forms and $c(\cdot): V_1\to \IR$ , $f(\cdot): V_2\to \IR$ are linear forms.

Here, the objective is to solve for the two unknown functions $u$ and $y$ . One can also imagine an arbitrary degree of coupling between $n$ variables with $n$ equations.

We introduce the following discretizations

\begin{equation} \begin{alignedat}{3} u_h &= \suml{J=1}{n_n^1} u^J N^J \qquad\qquad & v_h &= \suml{I=1}{n_n^1} u^I N^I \qquad\qquad & u_h, v_h\in V_{h1} \\ y_h &= \suml{L=1}{n_n^2} u^L N^L & w_h &= \suml{K=1}{n_n^2} u^K N^K & y_h, w_h\in V_{h2} \\ \end{alignedat} \tag{2}\end{equation}

where the corresponding number of shape functions are $n_n^1$ and $n_n^2$ , respectively.

Substituting the discretizations in (1), we obtain two linear systems of equations

\begin{equation} \begin{alignedat}{3} \suml{J=1}{n_n^1} a(N^J, N^I) \,u^J &+ \suml{L=1}{n_n^2} b(N^L, N^I) \, y^L &&= c(N^I) \\ \suml{J=1}{n_n^1} d(N^J, N^K) \,u^J &+ \suml{L=1}{n_n^2} e(N^L, N^K) \, y^L &&= f(N^K) \\ \end{alignedat} \tag{3}\end{equation}

for $I=1,\dots,n_n^1$ and $K=1,\dots,n_n^2$ .

We write this system as

\begin{equation} \boxed{ \begin{alignedat}{3} \BA \Bu &+ \BB\By &&= \Bc \\ \BD \Bu &+ \BE\By &&= \Bf \\ \end{alignedat} \eqor \begin{bmatrix} \BA & \BB \\ \BD & \BE \end{bmatrix} \begin{bmatrix} \Bu \\ \By \end{bmatrix} = \begin{bmatrix} \Bc \\ \Bf \end{bmatrix} } \htmlId{eq:coupledsystem1}{} \tag{4}\end{equation}

where the components of given matrices and vectors are defined as

\begin{equation} \begin{alignedat}{6} A^{I\!J} &:= a(N^J, N^I) \qquad & B^{I\!L} &:= b(N^L, N^I) \qquad & c^{I} &:= c(N^I)\\ D^{K\!J} &:= d(N^J, N^K) & E^{K\!L} &:= e(N^L, N^K) & f^{K} &:= f(N^K)\\ \end{alignedat} \tag{5}\end{equation}

Solution of (4) yields the unknown vectors $\Bu$ and $\By$ .

Onur Solmaz · Post · /2017/12/07

Time-Dependent Finite Elements

Time dependent problems are commonplace in physics, chemistry and many other disciplines. In this post, I’ll introduce the FE formulation of linear time-dependent problems and derive formulas for explicit and implicit Euler integration.

The weak formulation of a first order time-dependent problem reads:

Find $u \in V$ such that
$\begin{equation} m(\dot{u}, v; t) + a(u,v; t) = b(v; t) \htmlId{eq:timedependentweak1}{} \tag{1}\end{equation}$
for all $v \in V$ and $t \in [0,\infty)$ .

We can convert (1) into a system of equations

\begin{equation} \BM(t)\dot{\Bu} + \BA(t)\Bu = \Bb(t) \tag{2}\end{equation}

where the components of the matrices and vectors involved are calculated as

\begin{equation} \begin{aligned} M^{I\!J}(t) &= m(N^J, N^I; t) \\ A^{I\!J}(t) &= a(N^J, N^I; t) \\ b^{I}(t) &= b(N^I; t). \end{aligned} \tag{3}\end{equation}

If we further discretize in time with the finite difference $\dot{u} \approx [u_{n+1}-u_n]/{\Delta t}$ , linearity allows us to write

\begin{equation} \boxed{ m(\dot{u}, v; t) \approx \frac{1}{\Delta t} [m(u_{n+1}, v; t_{n+1}) - m(u_n, v; t_n)] } \htmlId{eq:discretetimedependent1}{} \tag{4}\end{equation}

This reflects on the system as

\begin{equation} \BM(t)\dot{\Bu} \approx \frac{1}{\Delta t} [\BM_{n+1}\Bu_{n+1} - \BM_n\Bu_n] \htmlId{eq:discretetimedependent2}{} \tag{5}\end{equation}

Here, $u_{n+1}:= u(x, t_{n+1})$ , $\BM_{n+1} = \BM(t_{n+1})$ and vice versa for $u_n$ and $\BM_n$ .

Explicit Euler Scheme

For the explicit Euler scheme, we substitute evaluate the remaining terms at $t_n$

\begin{equation} \frac{1}{\Delta t} [m(u_{n+1}, v; t_{n+1}) - m(u_n, v; t_n)] + a(u_n,v; t_n) = b(v; t_n) \quad \forall v \in V\,. \tag{6}\end{equation}

The corresponding system is

\begin{equation} \frac{1}{\Delta t} [\BM_{n+1}\Bu_{n+1} - \BM_n\Bu_n] + \BA_n\Bu_n = \Bb_n \tag{7}\end{equation}

The update equation becomes

\begin{equation} \boxed{ \Bu_{n+1} = \BM_{n+1}\inv [\BM_n\Bu_n + \Delta t(\Bb_n - \BA_n\Bu_n)] } \tag{8}\end{equation}

If $m$ is time-independent, that is $m(\dot{u}, v;t) = m(\dot{u}, v)$ , we have

\begin{equation} \Bu_{n+1} = \Bu_n + \Delta t\, \BM\inv(\Bb_n - \BA_n\Bu_n) \tag{9}\end{equation}

Implicit Euler Scheme

For the implicit Euler scheme, we substitute evaluate the remaining terms at $t_{n+1}$

\begin{equation} \frac{1}{\Delta t} [m(u_{n+1}, v; t_{n+1}) - m(u_n, v; t_n)] + a(u_n,v; t_{n+1}) = b(v; t_{n+1}) \quad \forall v \in V\,. \tag{10}\end{equation}

The corresponding system is

\begin{equation} \frac{1}{\Delta t} [\BM_{n+1}\Bu_{n+1} - \BM_n\Bu_n] + \BA_{n+1}\Bu_{n+1} = \Bb_{n+1} \tag{11}\end{equation}

The update equation becomes

\begin{equation} \boxed{ \Bu_{n+1} = [\BM_{n+1}+\Delta t \BA_{n+1}]\inv [\BM_n\Bu_n + \Delta t \,\Bb_{n+1}] } \tag{12}\end{equation}

If $m$ is time-independent, one can just substitute $\BM=\BM_{n+1}=\BM_n$ .

Example: Reaction-Advection-Diffusion Equation

The IBVP of a linear reaction-advection-diffusion problem reads

\begin{equation} \begin{alignedat}{4} \partd{u}{t} &= \nabla\dtp(\BD\nabla u) - \nabla\dtp(\Bc u) + ru + f \qquad&& \text{in} \qquad&& \Omega\times I\\ u &= \bar{u} && \text{on} && \del\Omega\times I\\ u &= u_0 && \text{in} && \Omega, t = 0 \\ \end{alignedat} \tag{13}\end{equation}

where $t\in I = [0,\infty)$ ,

$\BD$ is a second-order tensor describing the diffusivity of $u$ ,
$\Bc$ is a vector describing the velocity of advection,
$r$ is a scalar describing the rate of reaction,
and $f$ is a source term for $u$ .

The weak formulation is then

Find $u \in V$ such that
$\begin{equation} \int_\Omega \dot{u} v \dv = \int_\Omega [\nabla\dtp(\BD\nabla u) - \nabla\dtp(\Bc u) + ru + f] v \dv \tag{14}\end{equation}$
for all $v \in V$ and $t \in I$ .

We have the following integration by parts relationships:

\begin{equation} \int_\Omega \nabla \dtp(\BD\nabla u) v \dv = \cancel{\int_\Omega \nabla\dtp(v\BD\nabla u) \dv} - \int_\Omega (\BD\nabla u)\dtp\nabla v \dv \end{equation}

for the diffusive part and

\begin{equation} \int_\Omega \nabla\dtp(\Bc u) v \dv = \cancel{\int_\Omega \nabla \dtp (\Bc u v) \dv} - \int_\Omega u \Bc \dtp \nabla v \dv \tag{15}\end{equation}

for the advective part. The canceled terms are due to divergence theorem and the fact that $v=0$ on the boundary. Then our variational formulation is of the form (1) where

\begin{align*} m(\dot{u}, v) &= \int_\Omega \dot{u} v \dv \\ a(u, v) &= \int_\Omega (\BD\nabla u) \dtp \nabla v \dv - \int_\Omega u\Bc \dtp \nabla v \dv - \int_\Omega ruv \dv \\ b(v) &= \int_\Omega fv \dv \end{align*}

From these forms, we obtain the following system matrices and vector

\begin{align*} M^{I\!J} &= \int_\Omega N^J N^I \dv \\ A^{I\!J} &= \int_\Omega (\BD\BB^J) \dtp \BB^I \dv - \int_\Omega N^J\Bc \dtp \BB^I \dv - \int_\Omega r N^JN^I \dv \\ b^I &= \int_\Omega f N^I \dv \end{align*}

where $\BM$ is constant through time.

Onur Solmaz · Post · /2017/11/21

Vectorial Finite Elements

Many initial boundary value problems require solving for unknown vector fields, such as solving for displacements in a mechanical problem. Discretization of weak forms of such problems leads to higher-order linear systems which need to be reshaped to be solved by regular linear solvers. There are also more indices involved than a scalar problem, which can be confusing. In this post, I’ll try to elucidate the procedure by deriving for a basic higher-order system and giving an example.

The weak formulation of a linear vectorial problem reads

Find $\Bu\in V$ such that
$\begin{equation} a(\Bu, \Bv) = b(\Bv) \tag{1}\end{equation}$
for all $\Bv\in V$ .

Discretizations of vectorial problems requires the expansion of vectorial quantities as linear combinations of the basis vectors $\Be_i$ :

\begin{equation} \Bu = \suml{i=1}{\ndim} u_i \,\Be_i \htmlId{eq:discrete6}{} \tag{2}\end{equation}

where $\cbr{u_i}_{i=1}^{\ndim}$ are the components corresponding to the basis vectors and $\ndim=\dim V$ . Here, we chose Cartesian basis vectors for simplicity.

We can therefore express its discretization as

\begin{equation} \Bu_h = \suml{I=1}{\nnode} \Bu^I N^I = \suml{I=1}{\nnode}\suml{i=1}{\ndim} u^I_i \Be_i N^I. \htmlId{eq:discrete7}{} \tag{3}\end{equation}

Substituting discretized functions in the weak formulation, we obtain

\begin{equation} \begin{aligned} a(\Bu_h, \Bv_h) &= \suml{i=1}{\ndim}\suml{j=1}{\ndim} \, a(u_{h,j}\,\Be_j, v_{h,i}\,\Be_i)\\ &= \suml{I=1}{\nnode}\suml{J=1}{\nnode} \suml{i=1}{\ndim}\suml{j=1}{\ndim} u^J_j v^I_i \,a(\Be_j N^J, \Be_i N^I) \end{aligned} \tag{4}\end{equation}

and

\begin{equation} b(\Bv_h) = \suml{i=1}{\ndim} b(v_{h,i}\,\Be_i) = \suml{I=1}{\nnode} \suml{i=1}{\ndim} v^I_i b(\Be_i N^I). \tag{5}\end{equation}

We define the following arrays

\begin{equation} \boxed{ \begin{aligned} A^{I\!J}_{ij} &= a(\Be_j N^J, \Be_i N^I) \\ b^{I}_{i} &= b(\Be_i N^I). \end{aligned} } \tag{6}\end{equation}

Hence we can express the linear system

\begin{equation} a(\Bu_h,\Bv_h) = b(\Bv_h) \tag{7}\end{equation}

as

\begin{equation} \suml{I=1}{\nnode}\suml{J=1}{\nnode} \suml{i=1}{\ndim}\suml{j=1}{\ndim} u^J_j v^I_i \,A^{I\!J}_{ij} = \suml{I=1}{\nnode} \suml{i=1}{\ndim} v^I_i b^{I}_{i}. \tag{8}\end{equation}

For arbitrary $\Bv_h$ , this yields the following system of equations

\begin{equation} \boxed{ \suml{J=1}{\nnode} \suml{j=1}{\ndim} A^{I\!J}_{ij} \,u^J_j = b^{I}_{i} } \htmlId{eq:discrete8}{} \tag{9}\end{equation}

for $i=1,\dots,\ndim$ and $I=1,\dots,\nnode$ .

We reshape this higher-order system as shown in the previous post Reshaping Higher Order Linear Systems:

\begin{equation} \BAhat \Buhat = \Bbhat \tag{10}\end{equation}

by defining a map $i_d$ that maps original indices to the reshaped indices

\begin{equation} i_d := \left\{ \begin{array}{rl} [1,\nnode]\times[1,\ndim] & \to [1,\nnode\ndim]\\[1ex] (I,i) & \mapsto \ndim(I-1)+i\\ \end{array} \right. \tag{11}\end{equation}

where we used 1-based indexing of the arrays. We set

\begin{equation} \boxed{ \begin{alignedat}{3} \alpha &:= i_d(I,i) &&= \ndim(I-1) + i \\ \beta &:= i_d(J,j) &&= \ndim(J-1) + j \\ \end{alignedat} } \tag{12}\end{equation}

and write

\begin{equation} \hat{A}_{\alpha\beta} = A^{I\!J}_{ij} \quad,\quad \hat{u}_{\beta} = u^{J}_{j} \eqand \hat{b}_{\alpha} = b^{I}_{i} \tag{13}\end{equation}

The inverse index mapping can be obtained as shown in the previous post.

Example: Linear Elasticity

Our initial-boundary value problem is

\begin{equation} \begin{alignedat}{4} -\div\Bsigma &= \rho\Bgamma \qquad&& \text{in} \qquad&& \Omega\\ \Bu &= \bar{\Bu} && \text{on} && \del\Omega_u \\ \Bt &= \bar{\Bt} && \text{on} && \del\Omega_t \\ \end{alignedat} \tag{14}\end{equation}

The weak formulation reads

Find $\Bu\in V$ such that
$\begin{equation} -\int_\Omega \div\Bsigma \dtp \Bv \dv = \int_\Omega \rho \Bgamma\dtp\Bv \dv \tag{15}\end{equation}$
for all $\Bv\in V$ where $V=H^1(\Omega)$ .

We apply integration by parts on the left-hand side

\begin{equation} \int_\Omega \div\Bsigma\dtp\Bv \dv = \int_\Omega \div(\Bsigma\Bv) \dv - \int_\Omega \Bsigma : \nabla\Bv \dv \tag{16}\end{equation}

and apply the divergence theorem to the first resulting term:

\begin{equation} \int_\Omega \div(\Bsigma\Bv) \dv = \int_{\del\Omega_t} \bar{\Bt}\dtp\Bv \da \tag{17}\end{equation}

Substituting the linear stress $\Bsigma=\IC:\Bvareps=\IC:\nabla\Bu$ , we obtain the following variational forms:

\begin{align} \htmlId{eq:linelastdiscretebilinear}{} a(\Bu,\Bv) &= \int_\Omega \nabla\Bv:\IC:\nabla\Bu \dv \tag{18}\\ \htmlId{eq:linelastdiscretelinear}{} b(\Bv) &= \int_\Omega \rho \Bgamma\dtp\Bv \dv + \int_{\del\Omega_t} \bar{\Bt}\dtp\Bv \da \tag{19}\end{align}

We have the following discretizations of the unknown function and test function

\begin{equation} \Bu_h = \suml{J=1}{\nnode} \Bu^J N^J \eqand \Bv_h = \suml{I=1}{\nnode} \Bv^I N^I. \tag{20}\end{equation}

With the given discretizations, the matrix corresponding to (18) can be calculated as

\begin{equation} \begin{aligned} A^{I\!J}_{ij} = a(\Be_j N^J, \Be_i N^I) &= \int_\Omega \nabla(\Be_iN^I) : \IC : \nabla(\Be_jN^J) \dv \\ &= \int_\Omega (\Be_i\dyd \nabla N^I) : \IC : (\Be_j \dyd \nabla N^J) \dv \\ &= \int_\Omega \partd{N^I}{x_k} \, C_{ikjl} \, \partd{N^J}{x_l} \dv, \end{aligned} \tag{21}\end{equation}

and finally obtain

\begin{equation} \boxed{ A^{I\!J}_{ij} = \int_\Omega B^I_k \, C_{ikjl} \, B^J_l \dv \,. } \tag{22}\end{equation}

The vector corresponding to (19) is calculated as

\begin{equation} b^{I}_{i} = b(\Be_i N^I) = \int_{\del\Omega_t} \bar{\Bt}\dtp(\Be_iN^I) \da + \int_\Omega\rho\Bgamma\dtp(\Be_iN^I) \dv \tag{23}\end{equation}

which yields

\begin{equation} \boxed{ b^{I}_{i} = \int_{\del\Omega_t} \bar{t}_i N^I \da + \int_\Omega \rho \gamma_i N^I \dv } \tag{24}\end{equation}

Onur Solmaz · Post · /2017/11/20

Reshaping Higher Order Linear Systems

In vectorial problems, we end up with linear systems of higher order, such as

\begin{equation} \suml{k=1}{N} \suml{l=1}{M} A_{ijkl} \, u_{kl} = b_{ij} \htmlId{eq:1}{} \tag{1}\end{equation}

for $i=1,\dots,N$ and $j=1,\dots,M$ .

These systems cannot be solved readily with existing software. In order to be able to solve them with existing software, we need to reshape them by defining a matrix of matrices $\BAhat$ and vector of vectors $\Buhat$ and $\Bbhat$ :

\begin{equation} \underbrace{ \left[ \begin{array}{ccc|c|ccc} A_{1111} & \cdots & A_{111M} & & A_{11N1} & \cdots & A_{11NM} \\ \vdots & \ddots & \vdots & \cdots & \vdots & \ddots & \vdots \\ A_{1M11} & \cdots & A_{1M1M} & & A_{1MN1} & \cdots & A_{1MNM} \\[1ex] \hline & \vdots & &\ddots & & \vdots & \\ \hline &&&&&& \\[-1.5ex] A_{N111} & \cdots & A_{N11M} & & A_{N1N1} & \cdots & A_{N1NM} \\ \vdots & \ddots & \vdots & \cdots & \vdots & \ddots & \vdots \\ A_{NM11} & \cdots & A_{NM1M} & & A_{NMN1} & \cdots & A_{NMNM} \end{array} \right] }_{\BAhat} \underbrace{ \left[ \begin{array}{c} u_{11}\\ \vdots \\ u_{1M} \\[1ex] \hline \vdots\\ \hline \\[-1.5ex] u_{N1} \\ \vdots \\ u_{NM} \end{array} \right] }_{\Buhat} = \underbrace{ \left[ \begin{array}{c} b_{11}\\ \vdots \\ b_{1M} \\[1ex] \hline \vdots\\ \hline \\[-1.5ex] b_{N1} \\ \vdots \\ b_{NM} \end{array} \right]}_{\Bbhat} \tag{2}\end{equation}

This allows us to express the linear system as

\begin{equation} \BAhat \Buhat = \Bbhat \htmlId{eq:3}{} \tag{3}\end{equation}

Here, we reshape the system by defining a map $i_d$ that maps original indices to the reshaped indices

\begin{equation} i_d := \left\{ \begin{array}{rl} [1,N]\times[1,M] & \to [1,NM]\\[1ex] (i,j) & \mapsto M(i-1)+j\\ \end{array} \right. \tag{4}\end{equation}

where we used 1-based indexing of the arrays. We set

\begin{equation} \boxed{ \begin{alignedat}{3} \alpha &:= i_d(i,j) &&= M(i-1) + j \\ \beta &:= i_d(k,l) &&= M(k-1) + l \\ \end{alignedat} } \tag{5}\end{equation}

and write

\begin{equation} \hat{A}_{\alpha\beta} = A_{ijkl} \quad,\quad \hat{u}_{\beta} = u_{kl} \eqand \hat{b}_{\alpha} = b_{ij} \tag{6}\end{equation}

For reference, the inverse of the index mapping reads

\begin{equation} i_d^{-1} := \left\{ \begin{array}{rl} [1,NM] & \to [1,N]\times[1,M] \\[1ex] \alpha & \mapsto (1+(\alpha-\modop(\alpha,M))/M\;,\; \modop(\alpha,M)) \end{array} \right. \tag{7}\end{equation}

Thus, we have for our reshaped indices,

\begin{equation} \begin{aligned} j &= \modop(\alpha,M) \\ i &= 1+(\alpha-j)/M \end{aligned} \quad\eqand\quad \begin{aligned} l &= \modop(\beta,M) \\ k &= 1+(\beta-l)/M \end{aligned} \tag{8}\end{equation}

Expressed as a regular linear system (3), the higher-order system (1) can be solved with a linear solver such as LAPACK.

Onur Solmaz · Post · /2017/11/14

Linear Finite Elements

Beginning with this post, I’ll be publishing about the basics of finite element formulations, from personal notes that accumulated over the years. This one is about linear and scalar problems which came to be the “Hello World” for FE. Details regarding spaces and discretization are omitted for the sake of brevity. For those who want to delve into theory, I recommend “The Finite Element Method: Theory, Implementation, and Applications” by Larson and Bengzon.

The weak formulation of a canonical linear problem reads

Find $u\in V$ such that
$\begin{equation} a(u, v) = b(v) \htmlId{eq:femlinear1}{} \tag{1}\end{equation}$
for all $v \in V$ where $a(\cdot, \cdot)$ is a bilinear form and $b(\cdot)$ is a linear form.

We define the discretization of $u$ as

\begin{equation} u_h := \suml{J=1}{\nnode} u^J N^J ,\quad u_h \in V_h \quad\text{where}\quad V_h\subset V \tag{2}\end{equation}

The discretization $u_h$ is a linear combination of basis functions $N^J$ and corresponding scalars $u^J$ , $J=1,\dots,\nnode$ so that $V_h$ is a subset of $V$ . The discretization of (1) then reads

\begin{equation} a(u_h, v_h) = b(v_h) \quad \forall v_h \in V_h . \tag{3}\end{equation}

We then have

\begin{equation} a\rbr{\suml{J=1}{\nnode} u^J N^J, \suml{I=1}{\nnode} v^I N^I} = b\rbr{\suml{I=1}{\nnode} v^I N^I} \tag{4}\end{equation}

Using the linearity properties,

\begin{equation} a(\alpha u, \beta v) = \alpha\beta\, a(u,v) \eqand b(\alpha v) = \alpha b(v) \tag{5}\end{equation}

we obtain

\begin{equation} \suml{I=1}{\nnode} \suml{J=1}{\nnode} u^J v^I a(N^J, N^I) = \suml{I=1}{\nnode} v^I b(N^I) . \htmlId{eq:femlinear2}{} \tag{6}\end{equation}

For arbitrary test function values $v^I$ , we can express (6) as a system of $\nnode$ equations

\begin{equation} \suml{J=1}{\nnode} u^J a(N^J, N^I) = b(N^I) \htmlId{eq:femlinear3}{} \tag{7}\end{equation}

for $I = 1,2,\dots,\nnode$ . If we expand the summations as

\begin{alignat*}{6} & a(N^1, N^1) u^1 &&+ a(N^2, N^1) u^2 &&+ \cdots &&+ a(N^{\nnode}, N^1) u^{\nnode} &&\quad=\quad b(N^1) \\ & a(N^1, N^2) u^1 &&+ a(N^2, N^2) u^2 &&+ \cdots &&+ a(N^{\nnode}, N^2) u^{\nnode} &&\quad=\quad b(N^2) \\ & \qquad\vdots && \qquad\quad\;\vdots && \quad\;\;\vdots && \qquad\qquad\vdots && \qquad\qquad\vdots \\ & a(N^1, N^{\nnode}) u^1 &&+ a(N^2, N^{\nnode}) u^2 &&+ \cdots &&+ a(N^{\nnode}, N^{\nnode}) u^{\nnode} &&\quad=\quad b(N^{\nnode}) \end{alignat*}

we can see that the terms with $a$ constitute a matrix $\BA$ and the terms with $b$ constitute a vector $\Bb$ , allowing us to write

\begin{equation} \BA\Bu = \Bb \htmlId{eq:discrete9}{} \tag{8}\end{equation}

where we chose to express the unknown coefficients $u^I$ as a vector $\Bu = [u^1,u^2,\dots,u^{\nnode}]\tra$ .

\It can be seen that the components of the $\BA$ and $\Bb$ are defined as

\begin{equation} \boxed{ \Aelid{I\!J}{} = a(N^J,N^I) \eqand b^I = b(N^I), } \tag{9}\end{equation}

we can express the linear system as

\begin{alignat*}{6} & \Aelid{11}{} u^1 &&+ \Aelid{12}{} u^2 &&+ \cdots &&+ \Aelid{1\nnode}{} u^{\nnode} &&\quad=\quad b^1 \\ & \Aelid{21}{} u^1 &&+ \Aelid{22}{} u^2 &&+ \cdots &&+ \Aelid{2\nnode}{} u^{\nnode} &&\quad=\quad b^2 \\ & \quad\vdots && \qquad\;\vdots && \quad\;\;\vdots && \qquad\;\vdots && \qquad\quad\;\;\vdots \\ & \Aelid{\nnode 1}{} u^1 &&+ \Aelid{\nnode 2}{} u^2 &&+ \cdots &&+ \Aelid{\nnode\nnode}{} u^{\nnode} &&\quad=\quad b^{\nnode} \end{alignat*}

Note that with the given definitions, (7) becomes

\begin{equation} \boxed{ \suml{J=1}{\nnode} \Aelid{I\!J}{} \,u^J = b^I \quad\text{for}\quad I=1,2,\dots\nnode. } \htmlId{eq:discrete10}{} \tag{10}\end{equation}

Example: Poisson’s Equation

In the weak form of Poisson’s equation

\begin{equation} \begin{alignedat}{4} - \Var u &= f \quad && \text{in} \quad && \Omega \\ u &= 0 \quad && \text{on} \quad && \del\Omega \end{alignedat} \tag{11}\end{equation}

The weak formulation reads

Find $u\in V$ such that
$\begin{equation} - \int_\Omega \Delta(u) v \dv= \int_\Omega f v \dv \tag{12}\end{equation}$
for all $v\in V$ where $V=H^1_0(\Omega)$ .

Applying integration by parts and divergence theorem on the left-hand side

\begin{equation} \begin{aligned} \int_\Omega \Delta(u) v \dv &= \int_\Omega \nabla \dtp (\nabla (u) v) \dv - \int_\Omega \nabla u\dtp\nabla v \dv \\ &= \underbrace{\int_{\del\Omega} v (\nabla u\dtp\Bn) \da}_{v = 0 \text{ on } \del\Omega} - \int_\Omega \nabla u\dtp\nabla v \dv \\ \end{aligned} \tag{13}\end{equation}

We have the following variational forms:

\begin{equation} \begin{aligned} a(u,v) &= \int_{\Omega} \nabla u \dtp \nabla v \dv\\ b(v) &= \int_{\Omega} f \, v \dv\\ \end{aligned} \tag{14}\end{equation}

Following (7), we can calculate the stiffness matrix $\BA$ as

\begin{equation} \begin{aligned} \Aelid{I\!J}{} = a(N^J, N^I) &= \int_{\Omega} \nabla N^J \dtp \nabla N^I \dv \\ &= \int_{\Omega} \BB^J \dtp \BB^I \dv \end{aligned} \tag{15}\end{equation}

where we have defined the gradient of the basis functions as

\begin{equation} \BB^I := \nabla N^I\,. \tag{16}\end{equation}

Similarly, we integrate the force term into a vector $\Bb$ as

\begin{equation} \begin{aligned} b^I &= \int_{\Omega} f N^I \dv \end{aligned} \tag{17}\end{equation}

Onur Solmaz · Post · /2017/10/12

Balance Laws

Calculus is all about relating the change in one quantity to another quantity.

\Var A = B

Imagine it this way: You have a box full of marbles, and you decide to put some more in. $A$ is the variable representing the amount of marbles, while $B$ is the variable representing the amount of marbles that you put in. If you had $A_1$ marbles at the beginning, you have

A_2 = A_1+\Var A = A_1 + B

marbles following your action. This is the most fundamental algebraic pattern that characterizes balance laws.

Take the first law of thermodynamics for example—a.k.a. balance of energy. We have

\Var U = Q + W

where $U$ is the internal energy of a closed system, $Q$ is the amount of heat supplied to the system, and $W$ is the amount of work done on the system on its surroundings. Here, $A\equiv U$ and $B\equiv Q+W$ . Despite having three quantities, it is the combined effect of two which is related to the remaining quantity. Balance laws derived by physicists and chemists can get quite complex and hard to understand.

It’s always that the change in one quantity is related to the combined effect of remaining quantities. Keeping separate track of your main variable $A$ and affecting variables that compose $B$ gives you a mental model which helps you remember and even build your own balance laws.

Introducing Time

Let $A: t \to A(t)$ be a function of time. We can rewrite the equation in terms of the change in $A$ in a time period $\Var t$ :

\frac{\Var A}{\Var t} = C

where the new variable $C$ represents the change in quantity in $\Var t$ amount of time. In our previous analogy, $C$ is the amount of marbles put in, say, a minute. As $\Var t \to 0$ , we have

\deriv{A}{t} = C(t).

This prototypical balance law allows us to relate the rates of change of quantities.

Let’s introduce time into the balance of energy. The equation becomes

\deriv{U}{t} = P_T(t) + P_M(t)

where the new quantities $P_T$ and $P_M$ are called thermal power and mechanical power, representing the thermal and mechanical work done on the system per unit time, respectively. Given power functions and initial conditions, integrating them would give us the evolution of the internal energy through time.

Introducing Space

Let’s say we are not satisfied with an abstract box where the amount of stuff that goes in is measured automatically. We want to write a balance law over different shapes of bodies and we need to specify exactly where the stuff goes in and out.

To do that, we need to rephrase our laws to work over a continuous domain. The branch of physics that focuses on such problems is called continuum mechanics.

We introduce our spatial domain $\Omega$ and its boundary $\del\Omega$ . Our quantities now vary over both space and time, so we need to integrate them over the whole domain in order to relate them:

\ddt\int_\Omega a \dx = \int_\Omega b \dx + \int_{\del\Omega} \Bc \dtp \Bn \ds

where

$a(x,t)$ is the variable representing the main continuous quantity,
$b(x,t)$ is the variable representing the rate of change of the quantity inside the domain,
and $\Bc(x,t)$ is the variable representing the negative rate of change of the quantity on the boundary of the domain—negative due to surface normals $\Bn$ having outward direction by definition.

Notice that when we introduce space, our prototypical balance law needs an additional vectorial quantity, $\Bc$ . In physical laws, one needs to differentiate actions inside a body from actions on the surface of the body. That’s because one is over a volume and the other over an area, and they have to be integrated separately.

The area integral is actually a flux where the vectorial quantity $\Bc$ is penetrating the surface with a given direction. Given that it’s positive when stuff exits the domain, it’s called the efflux of the underlying quantity. Similarly, we name the rate of change field $b$ as the supply of the underlying quantity, because it being positive results in an increase.

The idea is to get rid of the integrals by a process called “localization”. In order to localize, we have to convert the surface integral into a volume integral using the divergence theorem:

\int_{\del\Omega} \Bc \dtp \Bn \ds = \int_\Omega \nabla\dtp\Bc \dx

Assuming $\Omega$ doesn’t move, we can also write

\ddt\int_\Omega a \dx = \int_\Omega \deriv{a}{t} \dx

Collecting the results, we have

\int_\Omega \deriv{a}{t} \dx = \int_\Omega b \dx + \int_\Omega \nabla\dtp\Bc \dx

Notice that all integrals are over $\Omega$ now. This allows us to make the balance law more strict by enforcing it point-wise:

\deriv{a}{t} = b + \nabla\dtp \Bc \quad \forall x \in \Omega

This is the localized version of the prototypical balance law that is used everywhere in continuum mechanics. Unfortunately, I can’t give the energy balance example, because it would require too many additional definitions. For that, I recommend the excellent Mathematical Foundations of Elasticity by Marsden and Hughes.

Conclusion

In physics and chemistry, one shouldn’t blindly memorize formulas, but try to see the underlying logic. In this case, I tried to elucidate balance laws, which all build upon the same algebraic and geometrical concepts. I went from discrete to continuous by introducing time and space to the equations, which became more complex but retained the same idea: putting things in a box and trying to calculate how that changes the contents.

Onur Solmaz · Post · /2017/09/20

Taylor and Volterra Series

In the theory of computational mechanics, there are a few operations used that are not taught in Calculus 101, which can be confusing without taking a lecture in calculus of variations. One of them is taking variations (a.k.a. Gateaux derivatives), akin to taking directional derivatives, but with functions of functions called functionals.

You need to take variations when you are linearizing a nonlinear problem for the purpose of solving with a numerical scheme. Linearization is the process of expanding a function or functional into a series, and discarding terms that are of order higher than linear—i.e. quadratic, cubic, quartic, etc. These expansions are called Taylor for functions, and Volterra for functionals.

Taylor Series

A function $f:\IR\to\IR$ can be expanded about a point $\bar{x}$ as a power series:

\begin{equation} \begin{aligned} f(x) &= f(\bar{x}) + \frac{\dif f}{\dif x}\evat_{\bar{x}} \frac{(x-\bar{x})}{1!} + \frac{\dif^2 f}{\dif x^2}\evat_{\bar{x}}\frac{(x-\bar{x})^2}{2!} + \frac{\dif^3 f}{\dif x^3}\evat_{\bar{x}}\frac{(x-\bar{x})^3}{3!} + \cdots \\ &= \suml{n=0}{\infty} \frac{\dif^n f}{\dif x^n} \evat_{\bar{x}}\frac{(x-\bar{x})^n}{n!} \end{aligned} \tag{1}\end{equation}

Letting $x$ be a perturbation $\var x$ from the expansion point $\bar{x}$ , that is $x\to\bar{x}+\var x$ , the series can also be phrased as follows

\begin{equation} \begin{aligned} f(\bar{x}+\var x) &= f(\bar{x}) + \frac{\dif f}{\dif x}\evat_{\bar{x}} \frac{\var x}{1!} + \frac{\dif^2 f}{\dif x^2}\evat_{\bar{x}}\frac{\var x^2}{2!} + \frac{\dif^3 f}{\dif x^3}\evat_{\bar{x}}\frac{\var x^3}{3!} + \cdots \\ &= \suml{n=0}{\infty} \frac{\dif^n f}{\dif x^n} \evat_{\bar{x}}\frac{\var x^n}{n!} \end{aligned} \htmlId{eq:2}{} \tag{2}\end{equation}

This is what is taught in Calculus 101 and everyone knows. Now for the part that you may have missed:

Variation

Let $X$ be the space of functions $\IR\to\IR$ . The variation of a functional $F\in X$ is defined as

\begin{equation} \boxed{ \varn{F(u)}{u}{v} := \lim_{\eps\to 0} \frac{F(u+\eps v) - F(u)}{\eps} \equiv \deriv{}{\eps} F(u + \eps v) \evat_{\eps = 0} } \tag{3}\end{equation}

where $v \in X$ is called the perturbation of the variation. This operation is analogous to taking the directional derivative of a function.

Shorthand notation

When working with variational formulations, writing out variations can be a bit of a hassle if there are many symbols involved. Therefore we use the following shorthand for variations:

\begin{equation} \Var F := \varn{F(u)}{u}{v} \tag{4}\end{equation}

Here, we assume that there is no chance of confusing the varied function or perturbation. We use this shorthand in contexts where the perturbation does not play an important role.

The shorthand for evaluation is

\begin{equation} \bar{F} := F(\bar{u}) \eqand \bar{\Var} F := \varn{F(u)}{u}{v}\evat_{\bar{u}} \tag{5}\end{equation}

where there is no risk of confusion for $\bar{u}\in X$ .

Volterra Series

Let $X$ be the space of functions $\IR\to\IR$ . Analogous to the Taylor series, a functional $F\in X$ can be expanded about a point $\bar{u}$ as a power series:

\begin{equation} \boxed{ \begin{aligned} F(\bar{u}+v) &= F(\bar{u}) + \frac{1}{1!} \varn{F(u)}{u}{v}\evat_{\bar{u}} + \frac{1}{2!} D^2_u F(u) \dtp v^2 \evat_{\bar{u}} + \frac{1}{3!} D^3_u F(u) \dtp v^3 \evat_{\bar{u}} + \cdots \\ &= \suml{n=0}{\infty} \frac{1}{n!} D^n_u F(u) \dtp v^n \evat_{\bar{u}} \end{aligned} } \htmlId{eq:6}{} \tag{6}\end{equation}

where $v\in X$ is the perturbation of the expansion. This is called the Volterra series expansion of $F$ . Verbally, the Volterra series expansion of a functional about a function is the infinite sum of the variations of the functional with increasing degree, evaluated at that function, each divided by the factorial of the degree.

In shorthand notation, the expansion is rendered

\begin{equation} \boxed{ F = \bar{F} + \frac{\bar{\Var} F}{1!} + \frac{\bar{\Var}^2 F}{2!} + \frac{\bar{\Var}^3 F}{3!} + \cdots } \htmlId{eq:7}{} \tag{7}\end{equation}

To me, there is an elegance in (7) that is not reflected in (2) or (6).

Onur Solmaz · Post · /2017/07/22

Isomorphisms in Linear Mappings between Vector Spaces

Equipping a vector space with an inner product results in a natural isomorphism $\CV\to\CV^\ast$ , where the metric tensor can be interpreted as the linear mapping $\Bg:\CV\to\CV^\ast$ and its inverse $\Bg\inv:\CV^\ast\to\CV$ .

Notation: Given two real vector spaces $\CV$ and $\CW$ , we denote their inner products as $\dabrn{\cdot,\cdot}_{\CV}$ and $\dabrn{\cdot,\cdot}_{\CW}$ respectively. Given vectors $\Bv\in\CV$ and $\Bw\in\CW$ , we define their lengths as

\begin{equation} \Norm{\Bv}_{\CV} = \sqrt{\dabrn{\Bv,\Bv}_{\CV}} \eqand \Norm{\Bw}_{\CW} = \sqrt{\dabrn{\Bw,\Bw}_{\CW}}. \tag{1}\end{equation}

Regarding $\CV$ and $\CW$ ,

their bases are denoted $\cbrn{\BE_A}$ and $\cbrn{\Be_a}$ ,
their dual bases are denoted $\cbrn{\BE^A}$ and $\cbrn{\Be^a}$ ,
their metrics are denoted $\BG$ and $\Bg$ with the components $G_{AB}=\dabrn{\BE_A,\BE_B}_{\CV}$ and $g_{ab}=\dabrn{\Be_a,\Be_b}_{\CW}$ ,

respectively. Here, the indices pertaining to $\CV$ are uppercase $(ABC\dots)$ and the indices pertaining to $\CW$ are lowercase $(abc\dots)$ .

Definition: Let $\BP:\CV\to\CW$ be a linear mapping. Then the transpose, or adjoint of $\BP$ , written $\BP\tra$ , is the linear mapping

\begin{equation} \boxed{ \BP\tra: \CW\to\CV \quad\text{such that}\quad \dabrn{\Bv,\BP\tra\Bw}_{\CV} = \dabrn{\BP\Bv,\Bw}_{\CW} } \tag{2}\end{equation}

for all $\Bv\in\CV$ and $\Bw\in\CW$ . Carrying out the products,

\begin{equation} G_{BA} v^B (P\tra){}^{A}{}_{d} w^d = g_{ab} P{}^{b}{}_{C}v^Cw^a. \tag{3}\end{equation}

For arbitrary $\Bv$ and $\Bw$ ,

\begin{equation} G_{BA} (P\tra){}^{A}{}_{a} = g_{ab} P{}^{b}{}_{A} \tag{4}\end{equation}

from which we can obtain the components of the transpose as

\begin{equation} \boxed{ (P\tra){}^{A}{}_{a} = g_{ab} P{}^{b}{}_{B} G^{AB} \eqwith \BP\tra = (P\tra){}^{A}{}_{a} \BE_A\dyd\Be^a . } \tag{5}\end{equation}

If $\BB:\CV\to\CV$ is a linear mapping, it is called symmetric if $\BB=\BB\tra$ .

Definition: Let $\BP:\CV\to\CW$ be a linear mapping. Then the dual of $\BP$ is a metric independent mapping

\begin{equation} \boxed{ \BP^\ast: \CW^\ast\to\CV^\ast \quad\text{such that}\quad \abrn{\Bv,\BP^\ast\Bbeta}_{\CV} = \abrn{\BP\Bv,\Bbeta}_{\CW} } \tag{6}\end{equation}

defined through natural pairings for all $\Bv\in\CV$ and $\Bbeta\in\CW^\ast$ . Carrying out the products,

\begin{equation} v^A (P^\ast){}_{A}{}^{a} \beta_a = P{}^{b}{}_{B} v^B \beta_b. \tag{7}\end{equation}

For arbitrary $\Bv$ and $\Bbeta$ , we obtain the components of the dual mapping as

\begin{equation} \boxed{ (P^\ast){}_{A}{}^{a} = P{}^{a}{}_{A} \eqwith \BP^\ast = (P^\ast){}_{A}{}^{a} \BE^A\dyd\Be_a = P{}^{a}{}_{A} \BE^A\dyd\Be_a . } \tag{8}\end{equation}

To fully appreciate the symmetry that originates from the duality, we can think of not just the mappings between $\CV$ and $\CW$ , but also between their dual spaces. To this end we can enumerate four mappings corresponding to $\cbr{\CV,\CV^\ast}\to\cbr{\CW,\CW^\ast}$ and their duals, corresponding to $\cbr{\CW,\CW^\ast}\to\cbr{\CV,\CV^\ast}$ . Their definitions can be found in the table below.

Tensors $\BP$ , $\BQ$ , $\BR$ and $\BS$ as linear mappings (top), and their duals $\BP^\ast$ , $\BQ^\ast$ , $\BR^\ast$ and $\BS^\ast$ (bottom). In the respective tables, the first row displays the tensor spaces, basis vectors and components of the subsequent mappings, and the second and third row display the representations of the tensor as linear and bilinear mappings respectively. The results of the mappings are given in the mapping, matrix and index representations respectively. The mappings are over vectors $\Bv\in\CV$ , $\Bw\in\CW$ and one-forms $\Balpha\in\CV^\ast$ , $\Bbeta\in\CW^\ast$ .
$\BP\in\CW \dyd\CV^\ast$ $P^a_{\idxsep A}=\BP(\Be^a,\BE_A)$ $\BP = P^{a}_{\idxsep A}\, \Be_a \dyd \BE^A$	$\BQ\in\CW^\ast \dyd\CV^\ast$ $Q_{aA}=\BQ(\Be_a,\BE_A)$ $\BQ = Q_{aA}\, \Be^a \dyd \BE^A$	$\BR\in\CW \dyd\CV$ $R^{aA}=\BR(\Be^a,\BE^A)$ $\BR = R^{aA}\, \Be_a \dyd \BE_A$	$\BS\in\CW^\ast \dyd\CV$ $S_a^{\idxsep A}=\BS(\Be_a,\BE^A)$ $\BS = S_{a}^{\idxsep A}\, \Be^a \dyd \BE_A$
Mappings
$\BP: \CV \to \CW$ $\begin{aligned} \Bv &\mapsto \BP(\Be^a,\Bv) \Be_a \\ &= \BP\Bv \\ v^A\BE_A &\mapsto P^a_{\idxsep A} v^A \Be_a \end{aligned}$	$\BQ: \CV \to \CW^\ast$ $\begin{aligned} \Bv &\mapsto \BQ(\Be_a,\Bv) \Be^a \\ &= \BQ\Bv \\ v^A\BE_A &\mapsto Q_{aA} v^A \Be^a \end{aligned}$	$\BR: \CV^\ast \to \CW$ $\begin{aligned} \Balpha &\mapsto \BR(\Be^a,\Balpha) \Be_a \\ &= \BR\Balpha\tra \\ \alpha_A \BE^A &\mapsto R^{aA} \alpha_A \Be_a \end{aligned}$	$\BS: \CV^\ast \to \CW^\ast$ $\begin{aligned} \Balpha &\mapsto \BS(\Be_a,\Balpha) \Be^a \\ &= \BS\Balpha\tra \\ \alpha_A \BE^A &\mapsto S_a^{\idxsep A} \alpha_A \Be^a \end{aligned}$
$\BP: \CW^\ast \times \CV \to \IR$ $\begin{aligned} (\Bbeta,\Bv) &\mapsto \BP(\Bbeta,\Bv) \\ &=\Bbeta\BP\Bv \\ &= \beta_a P^{a}_{\idxsep A} v^A \end{aligned}$	$\BQ: \CW \times \CV \to \IR$ $\begin{aligned} (\Bw,\Bv) &\mapsto \BQ(\Bw,\Bv) \\ &= \Bw\tra\BQ\Bv \\ &= w^a Q_{aA} v^A \end{aligned}$	$\BR: \CW^\ast \times \CV^\ast \to \IR$ $\begin{aligned} (\Bbeta,\Balpha) &\mapsto \BR(\Bbeta,\Balpha) \\ &= \Bbeta\BR\Balpha\tra \\ &= \beta_a R^{aA} \alpha_A \end{aligned}$	$\BS: \CW \times \CV^\ast \to \IR$ $\begin{aligned} (\Bw,\Balpha) &\mapsto \BS(\Bw,\Balpha) \\ &= \Bw\tra\BS\Balpha\tra \\ &= w^a S_a^{\idxsep A} \alpha_A \end{aligned}$
Duals
$\BP^\ast\in \CV^\ast \dyd\CW$ $P^{\ast \, a}_A=\BP^\ast(\BE_A,\Be^a)$ $\BP^\ast = P^{\ast \, a}_A \, \BE^A \dyd \Be_a$	$\BQ^\ast\in \CV^\ast \dyd\CW^\ast$ $Q^\ast_{Aa}=\BQ^\ast(\BE_A,\Be_a)$ $\BQ^\ast = Q^\ast_{Aa}\, \BE^A \dyd \Be^a$	$\BR^\ast\in\CV \dyd \CW$ $R^{\ast Aa}=\BR^\ast(\BE^A,\Be^a)$ $\BR^\ast = R^{\ast Aa}\, \BE_A \dyd \Be_a$	$\BS^\ast\in\CV \dyd\CW^\ast$ $S^{\ast A}{}_{a}=\BS^\ast(\BE^A,\Be_a)$ $\BS^\ast = S^{\ast A}{}_{a}\, \BE_A \dyd \Be^a$
$\BP^\ast: \CW^\ast \to \CV^\ast$ $\begin{aligned} \Bbeta &\mapsto \BP^\ast(\BE_A,\Bbeta) \BE^A \\ &= \BP^\ast\Bbeta\tra \\ \beta_a\Be^a &\mapsto P^{\ast \, a}_A \beta_a \BE^A \end{aligned}$	$\BQ^\ast: \CW \to\CV^\ast$ $\begin{aligned} \Bw &\mapsto \BQ^\ast(\BE_A,\Bw) \BE^A \\ &= \BQ^\ast\Bw \\ w^a\Be_a &\mapsto Q^\ast_{Aa} w^a \BE^A \end{aligned}$	$\BR^\ast: \CW^\ast \to \CV$ $\begin{aligned} \Bbeta &\mapsto \BR^\ast(\BE^A,\Bbeta) \BE_A \\ &= \BR^\ast\Bbeta\tra \\ \beta_a\Be^a &\mapsto R^{\ast Aa} w^a \BE_A \end{aligned}$	$\BS^\ast: \CW \to\CV$ $\begin{aligned} \Bw &\mapsto \BS^\ast(\BE^A,\Bw) \BE_A \\ &= \BS^\ast\Bw \\ w^a\Be_a &\mapsto S^{\ast A}{}_{a} w^a \BE_A \end{aligned}$
$\BP^\ast: \CV \times \CW^\ast \to \IR$ $\begin{aligned} (\Bv,\Bbeta) &\mapsto \BP^\ast(\Bv,\Bbeta) \\ &= \Bv\tra\BP^\ast\Bbeta\tra \\ &= v^A P^{\ast \, a}_A \beta_a \end{aligned}$	$\BQ^\ast: \CV \times \CW \to \IR$ $\begin{aligned} (\Bv,\Bw) &\mapsto \BR^\ast(\Bv,\Bw) \\ &= \Bv\tra\BQ^\ast\Bw \\ &= v^A Q^\ast_{Aa} w^a \end{aligned}$	$\BR^\ast: \CV^\ast \times \CW^\ast \to \IR$ $\begin{aligned} (\Balpha,\Bbeta) &\mapsto \BS^\ast(\Balpha,\Bbeta) \\ &= \Balpha\BR^\ast\Bbeta\tra \\ &= \alpha_A R^{\ast Aa} \beta_a \end{aligned}$	$\BS^\ast: \CV^\ast \times \CW \to \IR$ $\begin{aligned} (\Balpha,\Bw) &\mapsto \BS^\ast(\Balpha,\Bw) \\ &= \Balpha\BS^\ast\Bw \\ &= \alpha_A S^{\ast A}{}_{a} w^a \end{aligned}$

The commutative diagrams pertaining to these mappings can be found in the figure below

Commutative diagrams involving the linear mappings

\BP,\BQ,\BR,\BS

and their dual

\BP^\ast,\BQ^\ast,\BR^\ast,\BS^\ast

based on the metrics

\BG

and

\Bg

of

\CV

and

\CW

.

Onur Solmaz · Post · /2017/07/13

Metrics and Natural Isomorphisms

The assignment of an inner product to a non-degenerate and finite-dimensional vector space $\CV$ , results in emergence of the natural isomorphism to its dual $\CV\to\CV^\ast$ , which means that the morphisms $\CV\to\CV^\ast$ and $\CV^\ast\to\CV$ are of the same structure and one is the inverse of the other. The notion of naturality (of an isomorphism) becomes most clear in the context of category theory; however it should be sufficient for now to say that a natural isomorphism between a vector space an its dual is one that is basis-independent. As the origin of the isomorphism, the inner product is encapsulated in an object called the metric, defined below, in order to make the resulting symmetry of the mappings more obvious.

In the context of differential geometry, the metric object is used synonymously with the inner product of a vector space. More specifically, the metric tensor

\begin{equation} \Bg:= \left\{ \begin{aligned} \CV \times \CV &\to \IR \\ (\Bv,\Bw) &\mapsto \dabrn{\Bv, \Bw} \end{aligned} \right. \tag{1}\end{equation}

of a real vector space $\CV$ is an object whose components contain the information necessary to linearly transform a vector to its covector. This operation is denoted by the symbol $\flat$ and reads

\begin{equation} \flat := \left\{ \begin{aligned} (\CV^\ast \to \IR) &\to (\CV \to \IR) \\ \text{or}\quad \CV &\to\CV^\ast \\ \Bv(\cdot) &\mapsto \Bg(\Bv, \cdot). \end{aligned} \right. \tag{2}\end{equation}

We simply define the one-form $\Bv^\flat$ as

\begin{equation} \Bv^\flat(\Bw) \equiv \Bg(\Bv,\Bw) = \dabrn{\Bv,\Bw}. \tag{3}\end{equation}

We input the basis vectors $\Be_a$

\begin{equation} \Bv^\flat(\Be_a) = \dabrn{v^b\Be_b, \Be_a} = \dabrn{\Be_b, \Be_a} v^b = \dabrn{\Be_a, \Be_b} v^b \tag{4}\end{equation}

and define the components of the metric tensor as

\begin{equation} \boxed{ g_{ab} = \dabrn{\Be_a, \Be_b} \eqwith \Bg = g_{ab}\, \Be^a\dyd \Be^b. } \tag{5}\end{equation}

We then simply say that the operator $\flat$ denotes an index lowering¹ through

\begin{equation} \Bv^\flat = \Bg\Bv \quad\text{and component-wise }\quad v_a = g_{ab} v^b. \tag{6}\end{equation}

Moreover, we can define the inverse of the metric tensor as

\begin{equation} \Bg\inv:= \left\{ \begin{aligned} \CV^\ast\times\CV^\ast &\to \IR \\ (\Balpha,\Bbeta) &\mapsto \dabrn{\Balpha, \Bbeta} \end{aligned} \right. \tag{7}\end{equation}

The operation of transforming a covector to its corresponding vector is denoted by the symbol $\sharp$ and reads

\begin{equation} \sharp := \left\{ \begin{aligned} (\CV \to \IR) &\to (\CV^\ast \to \IR) \\ \text{or}\quad \CV^\ast &\to\CV \\ \Balpha(\cdot) &\mapsto \Bg\inv(\cdot,\Balpha). \end{aligned} \right. \tag{8}\end{equation}

Here, the vector corresponding to the covector $\Balpha$ is denoted $\Balpha^\sharp$ and reads

\begin{equation} \Balpha^\sharp(\Bbeta) = \Bg\inv(\Bbeta,\Balpha) = \dabrn{\Bbeta, \Balpha} \tag{9}\end{equation}

We input the dual basis vectors $\Be^a$

\begin{equation} \Balpha^\sharp(\Be^a) = \dabrn{\Be^a, \alpha_b\Be^b} = \dabrn{\Be^a, \Be^b} \alpha_b \tag{10}\end{equation}

and define the components of the inverse metric $\Bg\inv$ as

\begin{equation} \boxed{ g^{ab} = \dabrn{\Be^a,\Be^b} \eqwith \Bg\inv = g^{ab}\,\Be_a\dyd\Be_b. } \tag{11}\end{equation}

Then the operator $\sharp$ denotes an index raising through

\begin{equation} \Balpha^\sharp = \Bg\inv\Balpha \quad\text{and component-wise }\quad \alpha^a = g^{ab}\alpha_b. \tag{12}\end{equation}

In some literature, the natural isomorphism $\CV\to\CV^\ast$ is called the musical isomorphism—which is also the origin of the notation introduced above—because the process of transforming a vector to its dual space and a covector to the original space is analogous to lowering and raising notes.

With the given definition of the metric, we can elaborate on the advantage of denoting inner products of different objects with different symbols. Whereas $\abrn{\cdot,\cdot}$ always denotes a natural pairing between a vector space and its dual, one can write $\dabrn{\cdot,\cdot}_{\CV}:\CV\times\CV\to\IR$ to denote an inner product of vectors and $\dabrn{\cdot,\cdot}_{\CV^\ast}:\CV^\ast\times\CV^\ast\to\IR$ to denote an inner product of covectors. Using the metric, we can link these notations as

\begin{equation} \begin{alignedat}{5} &\dabrn{\Bv,\Bw}_{\CV} &&= \abrn{\Bv, \Bw^\flat} &&= \abrn{\Bv^\flat, \Bw} &&= g_{ab} v^a w^b\\ &\dabrn{\Balpha,\Bbeta}_{\CV^\ast} &&= \abrn{\Balpha, \Bbeta^\sharp} &&= \abrn{\Balpha^\sharp, \Bbeta} &&= g^{ab}\alpha_a\beta_b\\ \end{alignedat} \tag{13}\end{equation}

for all $\Bv,\Bw\in\CV$ and $\Balpha,\Bbeta\in\CV^\ast$ . Similarly,

\begin{equation} \begin{alignedat}{4} &\abrn{\Bv, \Balpha} &&= \dabrn{\Bv^\flat, \Balpha}_{\CV^\ast} &&= \dabrn{\Bv, \Balpha^\sharp}_{\CV}. \end{alignedat} \tag{14}\end{equation}

Despite the symmetricity of the inner product, we choose to think of the first operand as a vector and the second as a covector in a natural pairing, as a convention.

The metric tensor has the following properties:

For orthonormal bases, the metric tensor equals the identity tensor, that is, $g_{ij}=\delta_{ij}$ .
The diagonal terms equal to the square of the lengths of the basis vectors, that is, $g_{ii}=\Norm{\Be_i}^2$ (no summation).
The off-diagonal terms are zero if the basis vectors are orthogonal. Specifically, $g_{ij}=0$ iff $\Be_i$ and $\Be_j$ are orthogonal.

In musical notation, the flat symbol $\flat$ is used to lower a note by one semitone, whereas the sharp symbol $\sharp$ is used to raise a note by one semitone. It is recommended to pronounce $\Bv^\flat$ as v-flat and $\Balpha^\sharp$ alpha-sharp. ↩

Onur Solmaz · Post · /2017/07/06 · HN

Duality of Vector Spaces

When I was learning about Continuum Mechanics for the first time, the covariance and contravariance of vectors confused the hell out of me. The concepts gain meaning in the context of Riemannian Geometry, but it was surprising to find that one doesn’t need to learn an entire subject to grasp the logic behind co-/contravariance. An intermediate knowledge of linear algebra is enough—that is, one has to be acquainted with the concept of vector spaces and one-forms.

The duality of co-/contravariance arises when one has to define vectors in terms of a non-orthonormal basis. The reason such terminology doesn’t show up in engineering education is that Cartesian coordinates are enough for most engineering problems. But every now and then, a complex problem with funky geometrical requirements show up, like one that requires measuring distances and areas on non-flat surfaces. Then you end up with dual vector spaces. I’ll try to give the basics of duality below.

Definition: Let $\CV$ be a finite-dimensional real vector space. The space $\CV^\ast = \CL(\CV,\IR)$ , defined as the the space of all one-forms $\Balpha:\CV\to\IR$ , is called the dual space to $\CV$ .

Let $B=\cbr{\Be_1,\dots,\Be_n}$ be a basis of $\CV$ . Any vector $\Bv\in\CV$ can be written in terms of $B$ as

\begin{equation} \Bv = a_1 \Be_1 + \cdots + a_n\Be_n \htmlId{eq:vectorrep1}{} \tag{1}\end{equation}

with the components $a_1,\dots,a_n\in\IR$ . For any $i=1,\dots,n$ , we can define the $i$ -th component $a_i$ by a one-form as

\begin{equation} \Be^i := \left\{ \begin{aligned} \CV &\to \IR \\ \Bv &\mapsto \Be^i(\Bv) = a_i \end{aligned}\right. \tag{2}\end{equation}

These elements are linear and thus are in the space $\CL(\CV,\IR)$ ¹. Given any basis $B=\setveci{\Be}$ , we call $B^\ast = \setveciup{\Be}$ the basis of $\CV^\ast$ dual to $B$ . The fact that $B^\ast$ really is a basis of $\CV^\ast$ can be proved by showing that $\Be^i$ are linearly independent. Then $\Bv$ has the following representation

\begin{equation} \Bv = \Be^1(\Bv)\, \Be_1 + \cdots + \Be^n(\Bv)\, \Be_n. \htmlId{eq:vectorrep2}{} \tag{3}\end{equation}

Instead of $a_i$ , it is practical to denote the components of $\Bv$ as $v^i$ , lightface of the same symbol with a raised index corresponding to the raised index of the dual basis:

\begin{equation} \Bv = v^1 \Be_1 + \cdots + v^n \Be_n \eqwith v^i = \Be^i(\Bv). \tag{4}\end{equation}

In fact, this convention is more compatible with the symmetry caused by the duality. This point will be more clear after the introduction of dual basis representation of one-forms.

Proposition: Each $\Be^i \in \CL(\CV,\IR)$ can be identified by its action on the basis $B$ :

\begin{equation} \Be^i(\Be_j) = \begin{cases} 1 & \text{if } i=j \\ 0 & \text{otherwise}. \end{cases} \htmlId{eq:dualbasis2}{} \tag{5}\end{equation}

Proof: For any $\Bv\in\CV$ , $\Be^i(\Bv)$ must give $v^i$ , the $i$ -th component of $\Bv$ . Setting $\Bv = \Be_j$ , one sees that $\Be^i(\Bv)=v^i = 1$ when $i=j$ , and is zero otherwise.

Geometrically, (5) implies that a basis vector is perpendicular to all the dual basis vectors, except its own dual.

Dual Basis Representation of One-Forms

Let $\Balpha$ be a one form in $\CV^\ast$ with the corresponding dual basis $\setveciup{\Be}$ . Then similar to a vector, $\Balpha$ has the following representation

\begin{equation} \begin{aligned} \Balpha(\cdot) &= \Balpha(\Be_1)\,\Be^1(\cdot) + \dots + \Balpha(\Be_n)\,\Be^n(\cdot) \\ &= \alpha_1 \Be^1(\cdot) + \dots + \alpha_n \Be^n(\cdot) \end{aligned} \tag{6}\end{equation}

where the components of the one-form $\Balpha$ are defined as

\begin{equation} \alpha_i = \Balpha(\Be_i). \tag{7}\end{equation}

Proof: We substitute (3) and obtain

\begin{equation} \begin{aligned} \Balpha(\Bv) &= \Balpha\rbr{\suml{i=1}{n} \Be^i(\Bv)\, \Be_i} = \suml{i=1}{n} \Balpha(\Be_i)\, \Be^i(\Bv) \\ \end{aligned} \tag{8}\end{equation}

using $\Balpha$ ‘s linearity.

Notation: Let $\CV$ be a finite-dimensional real vector space. For $\Bv\in\CV$ and $\Balpha\in\CV^\ast$

\begin{equation} \abrn{\cdot,\cdot} := \left\{\begin{aligned} \CV\times\CV^\ast &\to \IR \\ (\Bv, \Balpha) &\mapsto \Balpha(\Bv) \end{aligned}\right. \tag{9}\end{equation}

denotes the action of $\Balpha$ on $\Bv$ , and is called a natural pairing or dual pairing between a vector space and its dual. It is of the essence to understand that $\abrn{\cdot,\cdot}$ does not denote an inner product in $\CV$ ; that is, $\abr{\Bv,\Balpha}$ means $\Balpha(\Bv)$ .

With this notation, (3) can be written as

\begin{equation} \Bv = \abrn{\Bv,\Be^1}\, \Be_1 + \cdots + \abrn{\Bv,\Be^n}\, \Be_n. \tag{10}\end{equation}

and (5) as

\begin{equation} \abrn{\Be_i, \Be^j} = \delta_{ij}. \htmlId{eq:dualbasis1}{} \tag{11}\end{equation}

Using the convention that $\Be_i$ are column vectors and $\Be^i$ are row vectors, (11) can be rearranged in the following manner

\begin{equation} \left[ \begin{array}{ c|c|c|c } \Be_1&\Be_2&\cdots&\Be_n \end{array} \right]\inv = \left[ \begin{array}{ c } \Be^1 \\ \hline \Be^2 \\ \hline \vdots \\ \hline \Be^n \end{array} \right] \htmlId{eq:computedualbasis1}{} \tag{12}\end{equation}

which can be used to compute a dual basis.

Example: Given a two-dimensional vector space $\CV$ with a basis $\Be_1=[2,-0.5]\tra$ , $\Be_2=[1,1]\tra$ , we use (12) to compute

\begin{equation} \begin{bmatrix} 2 & 1 \\ -0.5 & 1 \end{bmatrix}\inv = \begin{bmatrix} 0.4 & -0.4 \\ 0.2 & 0.8 \end{bmatrix} \tag{13}\end{equation}

and obtain the dual basis vectors as $\Be^1=[0.4,-0.4]$ and $\Be^2=[0.2,0.8]$ . The result is given in the following figure,

where one can see that $\Be_1\perp\Be^2$ , $\Be^1\perp\Be_2$ .

A body

\CB

embedded in

\IR^2

with curvilinear coordinates. Every point

\CP

at

\BX

has an associated two-dimensional vector space, called

\CB

's tangent space at

\BX

, denoted

\tang_{\BX}\CB

. The basis

\Be_i

corresponding to coordinates

\theta_i

are not necessarily orthogonal and can admit corresponding duals

\Be^i

, due to curvilinearity. The coordinates appear to be affine at the point's immediate vicinity, and thus in the tangent space.

The introduction of the dual space allows us to reinterpret a one-form $\Balpha$ as an object residing in the dual space. In fact, the canonical duality $\CV^{\ast\ast}=\CV$ states that every vector $\Bv$ can be interpreted as a functional on the space $\CV^\ast$ via

\begin{equation} \Bv:= \left\{ \begin{aligned} \CV^\ast &\to \IR \\ \Balpha &\mapsto \Bv(\Balpha) \text{ or } \abrn{\Bv, \Balpha} \end{aligned}\right. \tag{14}\end{equation}

Despite being denoted with bold letters, one-forms should not be confused with vectors. ↩