# The Lyapunov-Schmidt-Koiter asymptotic method using SymPy

Sébastien Brisard

► To cite this version:

Sébastien Brisard. The Lyapunov-Schmidt-Koiter asymptotic method using SymPy. 2023. [hal-04146928]

HAL Id: hal-04146928

https://univ-eiffel.hal.science/hal-04146928

Preprint submitted on 6 Jul 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# The Lyapunov-Schmidt-Koiter asymptotic method using SymPy

Sébastien Brisard

6/30/23
# Table of contents

<table><tr><td>1</td><td>Welcome!</td><td>4</td></tr><tr><td>2</td><td>Setting-up the mathematical stage</td><td>5</td></tr><tr><td>2.1</td><td>Mathematical setting</td><td>5</td></tr><tr><td>2.2</td><td>Kernel of the hessian of the energy</td><td>6</td></tr><tr><td>2.3</td><td>Two canonical variational problems</td><td>7</td></tr><tr><td>2.4</td><td>Additional symbols</td><td>8</td></tr><tr><td>3</td><td>Stability of the critical point</td><td>9</td></tr><tr><td>4</td><td>Setting-up the computational stage</td><td>11</td></tr><tr><td>4.1</td><td>Importing all necessary modules</td><td>11</td></tr><tr><td>4.2</td><td>Elimination of the derivatives of the jacobian w.r.t. λ</td><td>13</td></tr><tr><td>4.2.1</td><td>The term of order 0</td><td>14</td></tr><tr><td>4.2.2</td><td>The term of order 1</td><td>14</td></tr><tr><td>4.2.3</td><td>The term of order 2</td><td>15</td></tr><tr><td>4.2.4</td><td>The term of order 3</td><td>15</td></tr><tr><td>4.3</td><td>Elimination of the remaining mixed derivatives</td><td>15</td></tr><tr><td>4.4</td><td>Summary: final expression of the energy</td><td>17</td></tr><tr><td>4.5</td><td>Implementation in the lsk.energy module</td><td>19</td></tr><tr><td>5</td><td>Bifurcation equations</td><td>21</td></tr><tr><td>5.1</td><td>Introduction</td><td>21</td></tr><tr><td>5.2</td><td>The variational problem of order 1</td><td>23</td></tr><tr><td>5.3</td><td>The variational problem of order 2</td><td>23</td></tr><tr><td>5.4</td><td>The variational problem of order 3</td><td>25</td></tr><tr><td>6</td><td>Asymptotic expansion of the energy and its hessian</td><td>28</td></tr><tr><td>6.1</td><td>Asymptotic expansion of the energy</td><td>29</td></tr><tr><td>6.2</td><td>Asymptotic expansion of the hessian of the energy</td><td>33</td></tr><tr><td>7</td><td>Eigenmodes of the hessian of the energy</td><td>35</td></tr><tr><td>7.1</td><td>The eigenvalue problem</td><td>36</td></tr><tr><td>7.2</td><td>The variational problem of order 0</td><td>37</td></tr><tr><td>7.3</td><td>The variational problem of order 1</td><td>37</td></tr><tr><td>7.4</td><td>The variational problem of order 2</td><td>39</td></tr><tr><td>8</td><td>Asymmetric bifurcation</td><td>42</td></tr><tr><td>9</td><td>A particular case of symmetric bifurcation</td><td>44</td></tr></table>
<table><tr><td>10 The case of a single mode</td><td>46</td></tr><tr><td>10.1 Asymmetric bifurcation</td><td>46</td></tr><tr><td>10.2 Symmetric bifurcation</td><td>46</td></tr><tr><td>11 The lsk symbols module</td><td>48</td></tr><tr><td>References</td><td>50</td></tr></table>
# 1 Welcome!

These are my notes on the LSK method for the analysis of the stability and bifurcation(s) of a conservative system. These notes are based on several references: the initial PhD thesis of Warner Tjardus Koiter (1945) as well as some graphical illustrations from his lecture notes (W. T. Koiter and Heijden 2009). I enjoyed the concise presentation of Nguyen (2000) as well as the lecture notes of Triantafyllidis (2017). Finally, the chapter by Potier-Ferry (1987) helped me clear some issues. I used the direct approach presented in Appendix A of the paper by Chakrabarti et al. (2018) to derive the bifurcation equations.

When the system under consideration exhibits multiple simultaneous buckling modes, some derivations become a bit tedious. I therefore used the SymPy computer algebra system as much as possible. This required a few "tricks", and I believe that SymPy is now fully part of this work, as reflected by the new title.

I hope the reader will find these notes useful, even though there are still a few points which I do not fully understand (they are clearly indicated in the text). Should you want to report an error or suggest improvements, please do so by creating an issue on the Github repository.

These notes are available in two forms

* as an online website at https://sbrisard.github.io/LSK/

* as a PDF document at ....

The Quarto sources of these notes are available on the Github repo https://github.com/sbrisard/LSK. In particular, all SymPy simulations are available as Jupyter notebooks.

These notes by Sébastien Brisard are licensed under a Creative Commons Attribution 4.0 International License. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.
# 2 Setting-up the mathematical stage

In this chapter, we define the problem mathematically. Symbols and notation are introduced in 2.1, as well as the fundamental assumptions that we make. Then, the *kernel* of the hessian of the energy is introduced in 2.2 as a subspace of the space of admissible states.

## 2.1 Mathematical setting

The space of admissible states of the system under consideration is denoted $U$. It has the structure of a real vector space. The energy of the system is $\mathcal{E}(u, \lambda)$, where $\lambda$ denotes a loading parameter. It is assumed that the fundamental branch of the equilibrium diagram, $u^*(\lambda)$ is known. Then the energy is stationary with respect to the state $u$ along the whole branch. In other words, for all $\hat{u} \in U$

$$ \mathcal{E}_{,u}[u^*(\lambda), \lambda; \hat{u}] = 0, \qquad (2.1) $$

where $\mathcal{E}_{,u}(u, \lambda; \hat{u})$ denotes the (real) value of the differential of the energy $\mathcal{E}$ with respect to the state $u$, evaluated at $(u, \lambda)$, for the test function $\hat{u}$. Similarly, evaluation of the second-, third-, etc., order differential of the energy will be denoted $\mathcal{E}_{,uu}(u, \lambda; \hat{u}, \hat{v})$, $\mathcal{E}_{,uuu}(u, \lambda; \hat{u}, \hat{v}, \hat{w})$, etc. It is assumed that a finite value $\lambda_0 > 0$ of $\lambda$ can be found (critical load), such that

1. $\mathcal{E}_{,uu}[u(\lambda), \lambda] > 0$ for all $0 < \lambda < \lambda_0$,

2. $\mathcal{E}_{,uu}(u_0, \lambda_0) \ge 0$ but $\mathcal{E}_{,uu}(u_0, \lambda_0) \ne 0$.

The load $\lambda_0$ will be referred to as the *critical load*; similarly, the state $u_0$ of the system at the critical load will be referred to as the *critical state*; finally, the pair $(u_0, \lambda_0)$ is the *critical point* of the system. Assumption 1 implies that equilibria along the fundamental branch are *stable below the critical load*. Conversely, it results from assumption 3 that equilibrium points on the fundamental branch are *unstable above the critical load*. Stability at the critical load is yet undetermined.

The goal of these notes is to analyze all equilibrium paths that pass through the critical point $(u_0, \lambda_0)$.

We introduce the following notations

$$ u_0 = u^*(\lambda_0), \quad \dot{u}_0 = \left. \frac{du^*}{d\lambda} \right|_{\lambda=\lambda_0}, \quad \ddot{u}_0 = \left. \frac{d^2 u^*}{d\lambda^2} \right|_{\lambda=\lambda_0}, \quad \ddot{\ddot{u}}_0 = \dots, \quad \dddot{u}_0 = \dots $$

and

$$ \mathcal{E}_2 = \mathcal{E}_{,uu}(u_0, \lambda_0), \quad \mathcal{E}_3 = \mathcal{E}_{,uuu}(u_0, \lambda_0), \quad \mathcal{E}_4 = \mathcal{E}_{,uuuu}(u_0, \lambda_0). $$
Note that $\mathscr{E}_2, \mathscr{E}_3$ and $\mathscr{E}_4$ thus defined are bi-, tri- and quadrilinear forms, respectively. The following derivatives are also introduced

$$ \dot{\mathscr{E}}_2(\hat{u}, \hat{v}) = \left. \frac{d}{d\lambda} \mathscr{E}_{uu}[u^*(\lambda), \lambda; \hat{u}, \hat{v}] \right|_{\lambda=\lambda_0} = \mathscr{E}_{uuu}(u_0, \lambda_0; \dot{u}_0, \hat{u}, \hat{v}) + \mathscr{E}_{uu\lambda}(u_0, \lambda_0; \hat{u}, \hat{v}) \quad (2.2) $$

$$ \ddot{\mathscr{E}}_2(\hat{u}, \hat{v}) = \left. \frac{d^2}{d\lambda^2} \mathscr{E}_{uu}[u^*(\lambda), \lambda; \hat{u}, \hat{v}] \right|_{\lambda=\lambda_0} = \mathscr{E}_{uuuu}(u_0, \lambda_0; \dot{u}_0, \hat{u}, \hat{v}) + 2\mathscr{E}_{uuu\lambda}(u_0, \lambda_0; \dot{u}_0, \hat{u}, \hat{v}) \\ + \mathscr{E}_{uu\lambda\lambda}(u_0, \lambda_0; \hat{u}, \hat{v}) + \mathscr{E}_{uuu}(u_0, \lambda_0; \ddot{u}_0), $$

and, similarly, $\dot{\mathscr{E}}_3, \ddot{\mathscr{E}}_3$, etc.

## 2.2 Kernel of the hessian of the energy

The kernel of the hessian of the energy, $\mathscr{E}_2$, is defined as follows

$$ V = \{u \in U, \mathscr{E}_2(u, u) = 0\}. $$

Since $\mathscr{E}_2$ is a bilinear, symmetric and positive (but not positive definite!) form, $V$ is a vector subspace of $U$.

To prove this result, we must show that, for all $u, v \in V$ and $\alpha \in \mathbb{R}$, $w = u + \alpha v \in V$. From the bilinearity and symmetry of $\mathscr{E}_2$

$$ \mathscr{E}_2(w, w) = \mathscr{E}_2(u + \alpha v, u + \alpha v) = \mathscr{E}_2(u, u) + 2\alpha \mathscr{E}_2(u, v) + \alpha^2 \mathscr{E}_2(v, v), $$

Since $u, v \in \ker \mathscr{E}_2$, the first and the last term vanish, and the above identity reduces to

$$ \mathscr{E}_2(w, w) = 2\alpha \mathscr{E}_2(u, v) $$

The bilinear form $\mathscr{E}_2$ is positive, therefore the left-hand side is positive, for all values of $\alpha \in \mathbb{R}$. The quantity $\mathscr{E}_2(u, v) = 0$ is necessarily null, and $\mathscr{E}_2(w, w) = 0$, which proves that $w \in V$ and that $V$ is a vector subspace of $U$. The following characterization of $V$ holds

$$ v \in V \iff \mathscr{E}_2(v, \hat{u}) = 0 \quad \text{for all } \hat{u} \in U. \quad (2.3) $$

Indeed, if for all $\hat{u} \in U$, $\mathscr{E}_2(v, \hat{u}) = 0$, then in particular $\mathscr{E}_2(v, v) = 0$ and $v \in V$. Conversely, let $v \in V$, $\hat{u} \in U$ and $\alpha \in \mathbb{R}$. Similarly to the previous proof, we write that $\mathscr{E}_2(w, w) \ge 0$, with $w = \hat{u} + \alpha v$

$$ \mathscr{E}_2(w, w) = \mathscr{E}_2(\hat{u}, \hat{u}) + 2\alpha \mathscr{E}_2(\hat{u}, v) + \alpha^2 \mathscr{E}_2(v, v) = 2\alpha \mathscr{E}_2(\hat{u}, v) + \mathscr{E}_2(\hat{u}, \hat{u}) \ge 0, $$

($\mathscr{E}_2(v, v) = 0$ since $v \in V$). The above expression is of degree 1 in $\alpha$, with a constant sign. Therefore the linear term in $\alpha$ must vanish: $\mathscr{E}_2(\hat{u}, v) = 0$, which proves the characterization (2.3) of $V$.
It will be assumed in the remainder of these notes the dimension of V is finite: $m = \dim V < +\infty$; m is the multiplicity of the critical point. A (finite) basis ($v_1, ..., v_m$) of V can therefore be introduced, that is orthonormal in the sense of $\langle \bullet, \bullet \rangle$

$$ \langle v_i, v_j \rangle = \delta_{ij} \quad \text{for all} \quad i, j = 1, ..., m. $$

To close this section, we define the complementary subspace W, orthogonal to V for the scalar product $\langle \bullet, \bullet \rangle$

$$ U = V \otimes_{\perp} W \quad \text{and} \quad \langle \hat{v}, \hat{w} \rangle = 0 \quad \text{for all} \quad \hat{v} \in V \quad \text{and} \quad \hat{w} \in W. $$

### i On the sign of $\dot{E}_2$

It results from definition (2.2) that, for all $\hat{v} \in V$

$$ \begin{aligned} \mathcal{E}_{uu}[u^*(\lambda), \lambda; \hat{v}, \hat{v}] &= \mathcal{E}_2(\hat{v}, \hat{v}) + (\lambda - \lambda_0) \dot{\mathcal{E}}_2(\hat{u}, \hat{u}) + o(\lambda - \lambda_0) \\ &= (\lambda - \lambda_0) \dot{\mathcal{E}}_2(\hat{v}, \hat{v}) + o(\lambda - \lambda_0) \end{aligned} $$

For $\lambda < \lambda_0$, it has been assumed that the hessian is positive definite along the fundamental branch. From the above equation, it therefore results that $\dot{\mathcal{E}}_2(\hat{v}, \hat{v}) < 0$. In other words, $\dot{\mathcal{E}}_2$ is a negative definite form over V.

## 2.3 Two canonical variational problems

The bilinear form $\mathcal{E}_2$ is elliptic over W. Therefore, variational problems of the type: find $w \in W$ such that

$$ \mathcal{E}_2(w, \hat{w}) + \ell(\hat{w}) = 0 \quad \text{for all} \quad \hat{w} \in W $$

are well-posed for any linear form $\ell$ over W. In particular, for $\ell = 0$, the unique solution to the variational problem

$$ \mathcal{E}_2(w, \hat{w}) = 0 \quad \text{for all} \quad \hat{w} \in W $$

is $w = 0$. For $\ell(\bullet) = \dot{\mathcal{E}}_2(v_i, \bullet)$ and $\ell(\bullet) = \mathcal{E}_3(v_i, v_j, \bullet)$, $w_i \in W$ and $w_{ij} \in W$ are defined as the unique solutions in W of the following variational problems

$$ \mathcal{E}_2(w_i, \hat{w}) + \dot{\mathcal{E}}_2(v_i, \hat{w}) = 0 \quad (2.4) $$

and

$$ \mathcal{E}_2(w_{ij}, \hat{w}) + \mathcal{E}_3(v_i, v_j, \hat{w}) = 0, \quad (2.5) $$

for all $\hat{w} \in W$. These variational problems (and their solutions) will pop-up recurrently in what follows.
## 2.4 Additional symbols

We will make use of the following symbols

$$ \dot{E}_{ij} = \dot{\mathcal{E}}_2(v_i, v_j), \qquad (2.6) $$

$$ \ddot{E}_{ij} = \dot{\mathcal{E}}_2(v_i, v_j) + \dot{\mathcal{E}}_2(v_i, w_j) + \dot{\mathcal{E}}_2(v_j, w_i), \qquad (2.7) $$

$$ E_{ijk} = \mathcal{E}_3(v_i, v_j, v_k), \qquad (2.8) $$

$$ \dot{E}_{ijk} = \dot{\mathcal{E}}_3(v_i, v_j, v_k) + \dot{\mathcal{E}}_2(v_i, w_{jk}) + \dot{\mathcal{E}}_2(v_j, w_{ik}) + \dot{\mathcal{E}}_2(v_k, w_{ij}), \qquad (2.9) $$

$$ E_{ijkl} = \mathcal{E}_4(v_i, v_j, v_k, v_l) + \mathcal{E}_3(v_i, v_j, w_{kl}) + \mathcal{E}_3(v_i, v_k, w_{jl}) + \mathcal{E}_3(v_i, v_l, w_{jk}). \qquad (2.10) $$

These symbols define second, third and fourth-order, fully symmetric, tensors over V.

**i** Consistency of the above definitions

Since $\mathcal{E}_2(v, \bullet) = 0$ for all $v \in V$, definitions (2.6) and (2.8) of $\dot{E}_{ij}$ and $E_{ijk}$ can also be written

$$ \dot{E}_{ij} = \dot{\mathcal{E}}_2(v_i, v_j) + \mathcal{E}_2(v_i, w_j) + \mathcal{E}_2(v_j, w_i) $$

and

$$ E_{ijk} = \mathcal{E}_3(v_i, v_j, v_k) + \mathcal{E}_2(v_i, w_{jk}) + \mathcal{E}_2(v_j, w_{ik}) + \mathcal{E}_2(v_k, w_{ij}), $$

which shows the consistency with definitions (2.7) and (2.9) of $\ddot{E}_{ij}$ and $\dot{E}_{ijk}$.
## 3 Stability of the critical point

In this chapter, we discuss the stability of the critical point $(u_0, \lambda_0)$. To this end, we evaluate the potential energy in a neighboring state $u_0 + u$, where $u \in U$ is “small”. We have, to the fourth order

$$ \mathcal{E}(u_0 + u, \lambda_0) - \mathcal{E}(u_0, \lambda_0) = \frac{1}{2}\mathcal{E}_2(u, u) + \frac{1}{6}\mathcal{E}_3(u, u, u) + \frac{1}{24}\mathcal{E}_4(u, u, u, u) + o(\langle u, u \rangle^2), $$

where the linear term has been omitted, $u_0$ being a critical point of the energy. We now expand $u$ as $u = \xi v + \eta w$, with $\xi, \eta \in \mathbb{R}$ and $v \in V$ and $w \in W$ are fixed, orthogonal directions. Since $v \in V$, we have $\mathcal{E}_2(\lambda_0; v, \bullet) = 0$. Owing to the multi-linearity and symmetry of the successive differential of $\mathcal{E}$, the above expression expands as follows

$$ \begin{aligned} \mathcal{E}(u_0 + u, \lambda_0) - \mathcal{E}(u_0, \lambda_0) = & \frac{1}{2}\eta^2 \mathcal{E}_2(w, w) + \frac{1}{6}\xi^3 \mathcal{E}_3(v, v, v) \\ & + \frac{1}{2}\xi^2 \eta \mathcal{E}_3(v, v, w) + \frac{1}{2}\xi \eta^2 \mathcal{E}_3(v, w, w) \\ & + \frac{1}{6}\eta^3 \mathcal{E}_3(w, w, w) + \frac{1}{24}\xi^4 \mathcal{E}_4(v, v, v, v) \\ & + \frac{1}{6}\xi^3 \eta \mathcal{E}_4(v, v, v, w) + \frac{1}{4}\xi^2 \eta^2 \mathcal{E}_4(v, v, w, w) \\ & + \frac{1}{6}\xi \eta^3 \mathcal{E}_4(v, w, w, w) + \frac{1}{24}\eta^4 \mathcal{E}_4(w, w, w, w) \\ & + o\left[\left(\xi^2 + \eta^2\right)^2\right]. \end{aligned} $$

For the equilibrium to be stable, the above expression must be $\ge 0$ for all $\xi$ et $\eta$ small enough. Taking first $\eta = 0$, we get the following necessary conditions

$$ \mathcal{E}_3(v, v, v) = 0 \quad \text{and} \quad \mathcal{E}_4(v, v, v, v) \ge 0 \quad \text{for all} \quad v \in V. \qquad (3.1) $$

**i** Note

Note that the first of these two conditions is equivalent to $E_{ijk} = 0$, for all $i, j, k = 1, \dots, m$.

In other words, if there exists $v \in V$ such that $\mathcal{E}_3(v, v, v) \ne 0$ or $\mathcal{E}_4(v, v, v, v) < 0$, then the equilibrium is *unstable* at the critical point. The above conditions are not sufficient. Indeed, assuming conditions (3.1) to hold, we now take $\eta = \xi^2$

$$ \mathcal{E}(u_0 + u, \lambda_0) - \mathcal{E}(u_0, \lambda_0) = \frac{1}{2}\xi^4\left[\mathcal{E}_2(w, w) + \mathcal{E}_3(v, v, w) + \frac{1}{12}\mathcal{E}_4(v, v, v, v)\right] + o(\xi^4) $$

and we get the further necessary condition

$$ \mathcal{E}_2(w, w) + \mathcal{E}_3(v, v, w) + \frac{1}{12}\mathcal{E}_4(v, v, v, v) \ge 0 \quad \text{for all} \quad v \in V \quad \text{and} \quad w \in W. \qquad (3.2) $$
The direction $v \in V$ being fixed, the above expression is minimal when $w$ satisfies the following variational problem

$$2\mathcal{E}_2(w, \hat{w}) + \mathcal{E}_3(v, v, \hat{w}) = 0 \quad \text{for all} \quad \hat{w} \in W. \qquad (3.3)$$

Expanding $v \in V$ in the $(v_i)$ basis as follows: $v = \xi_i v_i$, it is observed that the solution to the above variational problem is $w = \frac{1}{2}\xi_i \xi_j w_{ij}$, where $w_{ij}$ is the solution to the elementary variational problem (2.5). For this value of $w$, condition (3.2) reads

$$[\mathcal{E}_4(v_i, v_j, v_k, v_l) - 3\mathcal{E}_2(w_{ij}, w_{kl})]\xi_i \xi_j \xi_k \xi_l \geq 0 \quad \text{for all} \quad \xi_1, \dots, \xi_m \in \mathbb{R},$$

which, in view of definition (2.10) of $E_{ijkl}$, is equivalent to

$$E_{ijkl} \xi_i \xi_j \xi_k \xi_l \geq 0 \quad \text{for all} \quad \xi_m, \dots, \xi_m \in \mathbb{R}.$$

Note that Eq. (3.2) implies $\mathcal{E}_4(\lambda_0; v, v, v, v) \geq 0$, which becomes a redundant necessary condition. Indeed, plugging $w = \xi_i \xi_j w_{ij}$ into Eq. (3.2) cancels the first two terms. To sum up, we have the following *necessary* conditions for stability

$$E_{ijk} \xi_i \xi_j \xi_k = 0 \quad \text{and} \quad E_{ijkl} \xi_i \xi_j \xi_k \xi_l \geq 0 \quad \text{for all} \quad \xi_m, \dots, \xi_m \in \mathbb{R}. \qquad (3.4)$$

<figure><img src="image_2.png" /></figure>

Warning

It is claimed in W. T. Koiter and Heijden (2009) that the above condition is also *sufficient* to ensure stability at the critical point. I can't prove this result, though.
# 4 Setting-up the computational stage

This chapter lays the ground for all symbolic calculations that are to follow. The SymPy library is imported and initialized in 4.1. Then, the energy of the system is rewritten with a minimum set of parameters.

## 4.1 Importing all necessary modules

We will use the Python library SymPy for symbolic mathematics and rely also on the IPython library (in particular, IPython display.Math) for pretty LaTeX output. Some useful functions are defined in the lsk.display module, which will be systematically imported. Also, all relevant symbols are defined in the lsk symbols module (see Chapter 11).

Without loss of generality, it will be assumed in all symbolic calculations that $\lambda_0 = 0$ and $u_0 = 0$. The general case $\lambda_0 \neq 0$ and $u_0 \neq 0$ is readily recovered through the substitution $\lambda \leftrightarrow \lambda - \lambda_0$ and $u \leftrightarrow u - u_0$.

The following developments involve the energy $(u, \lambda) \mapsto \mathscr{E}(u, \lambda)$ and its differentials at the critical point $(u_0, \lambda_0)$, as well as the fundamental path $\lambda \mapsto u^*(\lambda)$ and its derivatives at $\lambda = \lambda_0$. It will therefore be convenient to express $\mathscr{E}$ and $u^*$ as Taylor expansions with respect to $u$ and $\lambda$.

```python
from sympy import *
from lsk.display import *
from lsk symbols import *
```

We start with the Taylor expansion of the energy $\mathscr{E}$. We define its differentials at the critical point. These differentials are stored in a dictionary. Values are indexed with the order of the differentials with respect to $u$ and $\lambda$.

```python
d = {
    r"\E_{\{,\lambda\}}(u_0, \lambda_0)": E_\lambda,
    r"\E_{\{,uu\}}(u_0, \lambda_0)": E2,
    r"\E_{\{,u\lambda\}}(u_0, \lambda_0)": E_u\lambda,
    r"\E_{\{,\lambda\lambda\}}(u_0, \lambda_0)": E_\lambda\lambda,
    r"\E_{\{,uuu\}}(u_0, \lambda_0)": E3,
    r"\E_{\{,uu\lambda\}}(u_0, \lambda_0)": E_uuu\lambda,
    r"\E_{\{,u\lambda\lambda\}}(u_0, \lambda_0)": E_u\lambda\lambda,
    r"\E_{\{,\lambda\lambda\lambda\}}(u_0, \lambda_0)": E_\lambda\lambda\lambda,
    r"\E_{\{,uuuu\}}(u_0, \lambda_0)": E4,
```
r"\E_{,uuu\lambda\lambda\lambda}(u_0, \lambda\lambda\lambda_0)": E_{uuu\lambda},
r"\E_{,uu\lambda\lambda\lambda}(u_0, \lambda\lambda\lambda_0)": E_{uu\lambda\lambda},
r"\E_{,u\lambda\lambda\lambda\lambda\lambda}(u_0, \lambda\lambda\lambda_0)": E_{u\lambda\lambda\lambda},
r"\E_{,\lambda\lambda\lambda\lambda\lambda\lambda}(u_0, \lambda\lambda\lambda_0)": E_{\lambda\lambda\lambda\lambda},
}
display_latex_dict(d, num_COLS=3)

$$
\begin{align*}
\mathcal{E}_{,\lambda}(u_0, \lambda_0) &= \mathcal{E}_{\lambda} & \mathcal{E}_{,uu}(u_0, \lambda_0) &= \mathcal{E}_{2} & \mathcal{E}_{,u\lambda}(u_0, \lambda_0) &= \mathcal{E}_{u\lambda} \\
\mathcal{E}_{,\lambda\lambda}(u_0, \lambda_0) &= \mathcal{E}_{\lambda\lambda} & \mathcal{E}_{,uuu}(u_0, \lambda_0) &= \mathcal{E}_{3} & \mathcal{E}_{,uu\lambda}(u_0, \lambda_0) &= \mathcal{E}_{uu\lambda} \\
\mathcal{E}_{,u\lambda\lambda}(u_0, \lambda_0) &= \mathcal{E}_{u\lambda\lambda} & \mathcal{E}_{,\lambda\lambda\lambda}(u_0, \lambda_0) &= \mathcal{E}_{\lambda\lambda\lambda} & \mathcal{E}_{,uuuu}(u_0, \lambda_0) &= \mathcal{E}_{4} \\
\mathcal{E}_{,uuu\lambda}(u_0, \lambda_0) &= \mathcal{E}_{uuu\lambda} & \mathcal{E}_{,uu\lambda\lambda}(u_0, \lambda_0) &= \mathcal{E}_{uu\lambda\lambda} & \mathcal{E}_{,u\lambda\lambda\lambda}(u_0, \lambda_0) &= \mathcal{E}_{u\lambda\lambda\lambda} \\
\mathcal{E}_{,\lambda\lambda\lambda\lambda}(u_0, \lambda_0) &= \mathcal{E}_{\lambda\lambda\lambda\lambda} & & & &
\end{align*}
$$

where

• $\mathcal{E}_{\lambda}$, $\mathcal{E}_{\lambda\lambda}$, $\mathcal{E}_{\lambda\lambda\lambda}$ and $\mathcal{E}_{\lambda\lambda\lambda\lambda}$ are scalar quantities,

• $\mathcal{E}_{u\lambda}$, $\mathcal{E}_{u\lambda\lambda}$ and $\mathcal{E}_{u\lambda\lambda\lambda}$ are linear forms,

• $\mathcal{E}_2$, $\mathcal{E}_{uu\lambda}$ and $\mathcal{E}_{uu\lambda\lambda}$ are bilinear forms,

• $\mathcal{E}_3$ and $\mathcal{E}_{uuu\lambda}$ are trilinear forms,

• $\mathcal{E}_4$ is a quadrilinear form.

Important

Note that all these differentials are defined as SymPy scalars. A definition as a SymPy function (e.g. E2 = Function(r"\E_2")) would be more appropriate. However, SymPy would fail to account for multilinearity or symmetry of these forms. Therefore, we use the following trick: all multi-linear forms are defined as scalars, and the standard multiplication operator * means function application. In other words, E2 * (a * u + b * v) * w (resp. E2 * u * v - E2 * v * u) should be understood as E2(a * u + b * v, w) (resp. E2(u, v) - E2(v, u)). In both cases, the expressions would be correctly simplified.

Whether the symbols in an expression are true scalars or vectors (elements of $U$) should be clear from the context. For example, in the expression: $\lambda * E2 * u * v$, the first * is a true multiplication, while the other * refer to function application.

The energy $\mathcal{E}(u, \lambda)$ is now expressed as a Taylor expansion about the critical point. We use the function `create_E` that is defined in the `lsk` module. We will get back to the optional parameter `simplify_mixed_derivatives` later.

$$
\begin{align*}
u, \lambda &= \text{symbols}(\text{"u \lambda"}) \\
E &= (\lambda * E_{\lambda} + (E2 * u^{**2} + 2 * \lambda * E_{u\lambda} * u + \lambda^{**2} * E_{\lambda\lambda}) / 2 \\
  &\quad + (E3 * u^{**3} + 3 * \lambda * E_{uu\lambda} * u^{**2} + 3 * \lambda^{**2} * E_{u\lambda\lambda} * u + \lambda^{**3} * E_{\lambda\lambda\lambda}) / 6 \\
  &\quad + (E4 * u^{**4} + 4 * \lambda * E_{uuu\lambda} * u^{**3} + 6 * \lambda^{**2} * E_{uu\lambda\lambda} * u^{**2} \\
  &\quad \phantom{+ (E4 * u^{**4} + 4 * \lambda * E_{uuu\lambda} * u^{**3} +} + 4 * \lambda^{**3} * E_{u\lambda\lambda\lambda} * u + \lambda^{**4} * E_{\lambda\lambda\lambda\lambda}) / 24).\text{expand}()
\end{align*}
$$
R_star, terms_per_line=4)

$$
\begin{align*}
\mathcal{R}^*(\lambda; \hat{u}) ={}& \mathcal{E}_2 \dot{u}_0 \hat{u} \lambda + \frac{\mathcal{E}_{u\lambda\lambda\lambda} \hat{u} \lambda^3}{6} + \frac{\mathcal{E}_{u\lambda\lambda} \hat{u} \lambda^2}{2} + \mathcal{E}_{u\lambda} \hat{u} \lambda \\
& + \frac{\mathcal{E}_2 \ddot{u}_0 \hat{u} \lambda^2}{2} + \frac{\mathcal{E}_3 \dot{u}_0^2 \hat{u} \lambda^2}{2} + \frac{\mathcal{E}_{uu\lambda\lambda} \dot{u}_0 \hat{u} \lambda^3}{2} + \mathcal{E}_{uu\lambda} \dot{u}_0 \hat{u} \lambda^2 \\
& + \frac{\mathcal{E}_2 \ddot{u}_0 \hat{u} \lambda^3}{6} + \frac{\mathcal{E}_{uu\lambda\lambda} \ddot{u}_0 \hat{u} \lambda^4}{4} + \frac{\mathcal{E}_{uu\lambda} \ddot{u}_0 \hat{u} \lambda^3}{2} + \frac{\mathcal{E}_{uuu\lambda} \dot{u}_0^2 \hat{u} \lambda^3}{2} \\
& + \frac{\mathcal{E}_2 \ddot{u}_0 \hat{u} \lambda^4}{24} + \frac{\mathcal{E}_3 \dot{u}_0^2 \hat{u} \lambda^4}{8} + \frac{\mathcal{E}_4 \dot{u}_0^3 \hat{u} \lambda^3}{6} + \frac{\mathcal{E}_{uu\lambda} \ddot{u}_0 \hat{u} \lambda^4}{6} \\
& + \frac{\mathcal{E}_3 \ddot{u}_0 \dot{u}_0 \hat{u} \lambda^4}{6} + \frac{\mathcal{E}_3 \ddot{u}_0 \dot{u}_0 \hat{u} \lambda^3}{2} + \frac{\mathcal{E}_4 \ddot{u}_0 \dot{u}_0^2 \hat{u} \lambda^4}{4} + \frac{\mathcal{E}_{uuu\lambda} \ddot{u}_0 \dot{u}_0 \hat{u} \lambda^4}{2}
\end{align*}
$$

Of course, since $\lambda \mapsto u^*(\lambda)$ is an equilibrium path, we have $\mathcal{R}^*(\lambda; \bullet) = 0$ for all $\lambda$. Therefore, all coefficients of the above polynomial in $\lambda$ are null, which delivers expressions of $\mathcal{E}_{u\lambda}$, $\mathcal{E}_{u\lambda\lambda}$ and $\mathcal{E}_{u\lambda\lambda\lambda}$. Each term is analyzed in term below. Expressions of the mixed derivatives are to be stored in the mixed1 dictionary.

mixed1 = dict()

4.2.1 The term of order 0

This term is uniformly null and therefore delivers no informations.

assert R_starAttach(λ, 0) == 0

4.2.2 The term of order 1

This term delivers the following equation

eq = Eq(R_starAttach(λ, 1), 0)

display(eq)

$$
\mathcal{E}_2 \dot{u}_0 \hat{u} + \mathcal{E}_{u\lambda} \hat{u} = 0
$$

for all $\hat{u} \in \mathcal{U}$. This equation delivers the following expression of $\mathcal{E}_{u\lambda}(u_0, \lambda_0)$

sol = solve(eq, E_u\lambda)

mixed1[E_u\lambda] = sol[0]

display_latex_equation(E_u\lambda, mixed1[E_u\lambda])
$$
\mathcal{E}_{u\lambda} = -\mathcal{E}_2 \dot{u}_0
$$

4.2.3 The term of order 2

eq = Eq(R_star市場.coeff(λ, 2).subs(mixed1), 0)

display(eq)

$$
\frac{\mathcal{E}_2 \ddot{u}_0 \hat{u}}{2} + \frac{\mathcal{E}_3 \dot{u}_0^2 \hat{u}}{2} + \frac{\mathcal{E}_{u\lambda\lambda} \hat{u}}{2} + \mathcal{E}_{uu\lambda} \dot{u}_0 \hat{u} = 0
$$

for all $\hat{u} \in U$. This equation delivers the following expression of $\mathcal{E}_{u\lambda\lambda}(u_0, \lambda_0)$

sol = solve(eq, E_u\lambda\lambda)

mixed1[E_uλλ] = sol[0]

display_latex_equation(E_uλλ, mixed1[E_uλλ])

$$
\mathcal{E}_{u\lambda\lambda} = -\mathcal{E}_2 \ddot{u}_0 - \mathcal{E}_3 \dot{u}_0^2 - 2 \mathcal{E}_{uu\lambda} \dot{u}_0
$$

4.2.4 The term of order 3

eq = Eq(R_star市场市場.coeff(λ, 3).subs(mixed1), 0).expand()

display(eq)

$$
\frac{\mathcal{E}_2 \ddot{u}_0 \hat{u}}{6} + \frac{\mathcal{E}_3 \dot{u}_0 \dot{u}_0 \hat{u}}{2} + \frac{\mathcal{E}_4 \dot{u}_0^3 \hat{u}}{6} + \frac{\mathcal{E}_{u\lambda\lambda\lambda} \hat{u}}{6} + \frac{\mathcal{E}_{uu\lambda\lambda} \dot{u}_0 \hat{u}}{2} + \frac{\mathcal{E}_{uu\lambda} \ddot{u}_0 \hat{u}}{2} + \frac{\mathcal{E}_{uuu\lambda} \dot{u}_0^2 \hat{u}}{2} = 0
$$

for all $\hat{u} \in U$. This equation delivers the following expression of $\mathcal{E}_{u\lambda\lambda\lambda}(u_0, \lambda_0)$

sol = solve(eq, E_uλλλ)

mixed1[E_uλλλ] = sol[0]

display_latex_equation(E_uλλ, mixed1[E_uλλ])

$$
\mathcal{E}_{u\lambda\lambda} = -\mathcal{E}_2 \ddot{u}_0 - \mathcal{E}_3 \dot{u}_0^2 - 2 \mathcal{E}_{uu\lambda} \dot{u}_0
$$

4.3 Elimination of the remaining mixed derivatives

So far, we have found the following expressions
display_latex_dict(mixed1, num_COLS=1)

$$ \mathcal{E}_{u\lambda} = -\mathcal{E}_2 \dot{u}_0 $$

$$ \mathcal{E}_{u\lambda\lambda} = -\mathcal{E}_2 \ddot{u}_0 - \mathcal{E}_3 \dot{u}_0^2 - 2\mathcal{E}_{u\lambda\lambda} \dot{u}_0 $$

$$ \mathcal{E}_{u\lambda\lambda\lambda} = -\mathcal{E}_2 \dddot{u}_0 - 3\mathcal{E}_3 \ddot{u}_0 \dot{u}_0 - \mathcal{E}_4 \dot{u}_0^3 - 3\mathcal{E}_{u\lambda\lambda\lambda} \dot{u}_0 - 3\mathcal{E}_{u\lambda\lambda} \ddot{u}_0 - 3\mathcal{E}_{u\lambda\lambda\lambda} \dot{u}_0^2 $$

We want to get rid of the remaining mixed derivatives, namely: $\mathcal{E}_{u\lambda\lambda}$, $\mathcal{E}_{u\lambda\lambda}$ and $\mathcal{E}_{u\lambda\lambda\lambda}$. To do so, we introduce the derivatives $\dot{\mathcal{E}}_2$, $\ddot{\mathcal{E}}_2$ and $\dot{\mathcal{E}}_3$ defined in Chapter 2.

```python
E_uu_star = E.diff(u, 2).subs(u, u_star).expand()
E_uuu_star = E.diff(u, 3).subs(u, u_star).expand()
mixed2 = dict()
```

The mixed derivative $\mathcal{E}_{u\lambda\lambda}$ can first be expressed as a function of $\dot{\mathcal{E}}_2$.

```python
x = E_uu\lambda
lhs = E2.dot
rhs = E_uu_star市場(λ, 1)
```

```python
display_latex_equation trex, rhs)
```

$$ \dot{\mathcal{E}}_2 = \mathcal{E}_3 \dot{u}_0 + \mathcal{E}_{u\lambda\lambda} $$

```python
sol = solve(Eq(lhs, rhs), x)
mixed2[x] = sol[0]
```

```python
display_latex_equation(x, mixed2[x])
```

$$ \mathcal{E}_{u\lambda\lambda} = -\mathcal{E}_3 \dot{u}_0 + \dot{\mathcal{E}}_2 $$

Then, the expression of $\dot{\mathcal{E}}_3$ delivers an expression of the mixed derivative $\mathcal{E}_{u\lambda\lambda}$.

```python
x = E_uuu\lambda
lhs = E3.dot
rhs = E_uuu_star市場(λ, 1)
```

```python
display_latex_equation trex, rhs)
```

$$ \dot{\mathcal{E}}_3 = \mathcal{E}_4 \dot{u}_0 + \mathcal{E}_{u\lambda\lambda\lambda} $$

```python
sol = solve(Eq(lhs, rhs), x)
mixed2[x] = sol[0]
```

```python
display_latex_equation(x, mixed2[x])
```
$$ \mathcal{E}_{uuu\lambda} = -\mathcal{E}_4\dot{u}_0 + \dot{\mathcal{E}}_3 $$

Finally, $\ddot{\mathcal{E}}_2$ delivers an expression of the mixed derivative $\mathcal{E}_{uu\lambda\lambda}$.

$$ x = E_{uu\lambda\lambda} $$
$$ \text{lhs} = E2\_ddot $$
$$ \text{rhs} = 2 * E_{uu\_star}.coeff(\lambda, 2).subs(mixed2).expand() $$

display_latex_equation(lhs, rhs)

$$ \ddot{\mathcal{E}}_2 = \mathcal{E}_3\ddot{u}_0 - \mathcal{E}_4\dot{u}_0^2 + \mathcal{E}_{uu\lambda\lambda} + 2\dot{\mathcal{E}}_3\dot{u}_0 $$

sol = solve(Eq(lhs, rhs), x)

mixed2[x] = sol[0]

display_latex_equation(x, mixed2[x])

$$ \mathcal{E}_{uu\lambda\lambda} = -\mathcal{E}_3\ddot{u}_0 + \mathcal{E}_4\dot{u}_0^2 + \ddot{\mathcal{E}}_2 - 2\dot{\mathcal{E}}_3\dot{u}_0 $$

## 4.4 Summary: final expression of the energy

The following expressions were derived in 4.2

display_latex_dict(mixed1, num_COLS=1)

$$ \mathcal{E}_{u\lambda} = -\mathcal{E}_2\dot{u}_0 $$

$$ \mathcal{E}_{u\lambda\lambda} = -\mathcal{E}_2\ddot{u}_0 - \mathcal{E}_3\dot{u}_0^2 - 2\mathcal{E}_{u\lambda\lambda}\dot{u}_0 $$

$$ \mathcal{E}_{u\lambda\lambda\lambda} = -\mathcal{E}_2\ddot{u}_0 - 3\mathcal{E}_3\ddot{u}_0\dot{u}_0 - \mathcal{E}_4\dot{u}_0^3 - 3\mathcal{E}_{u\lambda\lambda\lambda}\dot{u}_0 - 3\mathcal{E}_{u\lambda\lambda}\ddot{u}_0 - 3\mathcal{E}_{uuu\lambda}\dot{u}_0^2 $$

and in 4.3

display_latex_dict(mixed2, num_COLS=1)

$$ \mathcal{E}_{uu\lambda} = -\mathcal{E}_3\dot{u}_0 + \dot{\mathcal{E}}_2 $$

$$ \mathcal{E}_{uuu\lambda} = -\mathcal{E}_4\dot{u}_0 + \dot{\mathcal{E}}_3 $$

$$ \mathcal{E}_{uu\lambda\lambda} = -\mathcal{E}_3\ddot{u}_0 + \mathcal{E}_4\dot{u}_0^2 + \ddot{\mathcal{E}}_2 - 2\dot{\mathcal{E}}_3\dot{u}_0 $$

Combining the above results allows to fully eliminate the mixed derivatives

mixed = {k: v.subs(mixed2).expand() for k, v in mixed1.items()}

mixed.update(mixed2)

display_latex_dict(mixed, num_COLS=1)
$$
\mathcal{E}_{u\lambda} = -\mathcal{E}_2 \dot{u}_0
$$

$$
\mathcal{E}_{u\lambda\lambda} = -\mathcal{E}_2 \ddot{u}_0 + \mathcal{E}_3 \dot{u}_0^2 - 2\dot{\mathcal{E}}_2 \dot{u}_0
$$

$$
\mathcal{E}_{u\lambda\lambda\lambda} = -\mathcal{E}_2 \ddot{u}_0 + 3\mathcal{E}_3 \ddot{u}_0 \dot{u}_0 - \mathcal{E}_4 \dot{u}_0^3 - 3\ddot{\mathcal{E}}_2 \dot{u}_0 - 3\ddot{u}_0 \dot{\mathcal{E}}_2 + 3\dot{\mathcal{E}}_3 \dot{u}_0^2
$$

$$
\mathcal{E}_{uu\lambda} = -\mathcal{E}_3 \dot{u}_0 + \dot{\mathcal{E}}_2
$$

$$
\mathcal{E}_{uuu\lambda} = -\mathcal{E}_4 \dot{u}_0 + \dot{\mathcal{E}}_3
$$

$$
\mathcal{E}_{uu\lambda\lambda} = -\mathcal{E}_3 \ddot{u}_0 + \mathcal{E}_4 \dot{u}_0^2 + \ddot{\mathcal{E}}_2 - 2 \dot{\mathcal{E}}_3 \dot{u}_0
$$

These expressions can be plugged into the expansion of the energy.

$$
E = E.subs(mixed).expand()
$$

$$
display_latex_long_equation(r"\E(u, \lambda)", E, terms_per_line=5)
$$

$$
\begin{aligned}
\mathcal{E}(u, \lambda) = & \frac{\mathcal{E}_2 u^2}{2} + \frac{\mathcal{E}_3 u^3}{6} + \frac{\mathcal{E}_{\lambda\lambda\lambda}\lambda^3}{6} + \frac{\mathcal{E}_{\lambda\lambda}\lambda^2}{2} + \mathcal{E}_{\lambda}\lambda \\
& + \frac{\mathcal{E}_4 u^4}{24} + \frac{\mathcal{E}_{\lambda\lambda\lambda\lambda}\lambda^4}{24} + \frac{\ddot{\mathcal{E}}_2 \lambda^2 u^2}{4} + \frac{\dot{\mathcal{E}}_2 \lambda u^2}{2} + \frac{\dot{\mathcal{E}}_3 \lambda u^3}{6} \\
& + -\frac{\mathcal{E}_2 \ddot{u}_0 \lambda^2 u}{2} - \mathcal{E}_2 \dot{u}_0 \lambda u + \frac{\mathcal{E}_3 \dot{u}_0^2 \lambda^2 u}{2} - \dot{\mathcal{E}}_2 \dot{u}_0 \lambda^2 u + \frac{\ddot{\mathcal{E}}_3 \dot{u}_0^2 \lambda^3 u}{2} \\
& + -\frac{\mathcal{E}_3 \ddot{u}_0 \lambda^2 u^2}{4} - \frac{\mathcal{E}_3 \dot{u}_0 \lambda u^2}{2} - \frac{\ddot{\mathcal{E}}_2 \dot{u}_0 \lambda^3 u}{2} - \frac{\ddot{u}_0 \ddot{\mathcal{E}}_2 \lambda^3 u}{2} - \frac{\dot{\mathcal{E}}_3 \dot{u}_0 \lambda^2 u^2}{2} \\
& + -\frac{\mathcal{E}_2 \dddot{u}_0 \lambda^3 u}{6} + \frac{\mathcal{E}_3 \dddot{u}_0 \dot{u}_0 \lambda^3 u}{2} - \frac{\mathcal{E}_4 \dot{u}_0^3 \lambda^3 u}{6} + \frac{\mathcal{E}_4 \dot{u}_0^2 \lambda^2 u^2}{4} - \frac{\mathcal{E}_4 \dot{u}_0 \lambda u^3}{6}
\end{aligned}
$$

From which we deduce the expression of the residual $\mathcal{E}_{,u}$

$$
E_u = E.diff(u)
$$

$$
display_latex_long_equation(r"\E_{\{u\}}(u, \lambda)", E_u, terms_per_line=5)
$$

$$
\begin{aligned}
\mathcal{E}_{,u}(u, \lambda) = & \mathcal{E}_2 u + \frac{\mathcal{E}_3 \dot{u}_0^2 \lambda^2}{2} + \frac{\mathcal{E}_3 u^2}{2} + \frac{\mathcal{E}_4 u^3}{6} + \dot{\mathcal{E}}_2 \lambda u \\
& + -\mathcal{E}_2 \dot{u}_0 \lambda + \frac{\ddot{\mathcal{E}}_2 \lambda^2 u}{2} - \dot{\mathcal{E}}_2 \dot{u}_0 \lambda^2 + \frac{\dot{\mathcal{E}}_3 \dot{u}_0^2 \lambda^3}{2} + \frac{\dot{\mathcal{E}}_3 \lambda u^2}{2} \\
& + -\frac{\mathcal{E}_2 \ddot{u}_0 \lambda^3}{6} - \frac{\mathcal{E}_2 \ddot{u}_0 \lambda^2}{2} - \frac{\mathcal{E}_4 \dot{u}_0^3 \lambda^3}{6} - \frac{\ddot{\mathcal{E}}_2 \dot{u}_0 \lambda^3}{2} - \frac{\ddot{u}_0 \ddot{\mathcal{E}}_2 \lambda^3}{2} \\
& + \frac{\mathcal{E}_3 \ddot{u}_0 \dot{u}_0 \lambda^3}{2} - \frac{\mathcal{E}_3 \ddot{u}_0 \lambda^2 u}{2} - \mathcal{E}_3 \dot{u}_0 \lambda u + \frac{\mathcal{E}_4 \dot{u}_0^2 \lambda^2 u}{2} - \dot{\mathcal{E}}_3 \dot{u}_0 \lambda^2 u \\
& + -\frac{\mathcal{E}_4 \dot{u}_0 \lambda u^2}{2}
\end{aligned}
$$

$$
E_{uu} = E \cdot \text{diff}(u, 2)
$$

$$
display_latex_long_equation(r"\E_{\{uu\}}(u, \lambda)", E_{uu}, terms_per_line=7)
$$
$$ \mathcal{E}_{uu}(u, \lambda) = \mathcal{E}_2 + \mathcal{E}_3 u + \frac{\mathcal{E}_4 \dot{u}_0^2 \lambda^2}{2} + \frac{\mathcal{E}_4 u^2}{2} + \frac{\ddot{\mathcal{E}}_2 \lambda^2}{2} + \dot{\mathcal{E}}_2 \lambda + \dot{\mathcal{E}}_3 \lambda u \\ + -\frac{\mathcal{E}_3 \ddot{u}_0 \lambda^2}{2} - \mathcal{E}_3 \dot{u}_0 \lambda - \mathcal{E}_4 \dot{u}_0 \lambda u - \dot{\mathcal{E}}_3 \dot{u}_0 \lambda^2 $$

## 4.5 Implementation in the lsk.energy module

This module exposes three functions

* create_E(u, λ) : asymptotic expansion of the energy $\mathcal{E}(u, \lambda)$,

* create_E_u(u, λ) : asymptotic expansion of the jacobian $\mathcal{E}_{u}(u, \lambda)$,

* create_E_u(u, λ) : asymptotic expansion of the hessian $\mathcal{E}_{uu}(u, \lambda)$,

* create_u_star(λ) : asymptotic expansion of the fundamental branch $u^*(\lambda)$,

where u and λ must be SymPy expressions.

```python
import lsk.energy
```

```python
%psource lsk.energy
```

```python
from lsk.symbols import *
```

```python
__mixed_derivatives = {
    E_uλ : -E2 * u0 Dot,
    E_uλλ : -E2 * u0_ddot - 2 * E2 Dot * u0 Dot + E3 * u0 Dot**2,
    E_uλλλ : (
        -E2 * u0_dddot
        - 3 * E2 Dot * u0_ddot
        - 3 * E2_ddot * u0 Dot
        + 3 * E3 * u0 Dot * u0_ddot
        + 3 * E3 Dot * u0 Dot**2
        - E4 * u0 Dot**3
    ),
    E_uuλ : E2 Dot - E3 * u0 Dot,
    E_uuλλ : E4 * u0 Dot**2 - 2 * E3 Dot * u0 Dot - E3 * u0 Ddot + E2 Ddot,
    E_uuλλλ : E3 Dot - E4 * u0 Dot,
}
```

```python
def create_E(u, λ):
    out = (λ * E_λ + (E2 * u**2 + 2 * λ * E_uλ * u + λ**2 * E_λλ) / 2
           + (E3 * u**3 + 3 * λ * E_uuλ * u**2 + 3 * λ**2 * E_uλλ * u
           + λ**3 * E_λλλ) / 6
           + (E4 * u**4 + 4 * λ * E_uuλλ * u**3 + 6 * λ**2 * E_uuλλλ * u**2
           + 4 * λ**3 * E_uλλλ * u + λ**4 * E_λλλλ) / 24)
```
return out.subs(__mixed_derivatives).expand()

def create_E_u(u, λ):
    out = (E2 * u + λ * E_uλ
           + (E3 * u**2 + 2 * λ * E_uuλ * u + λ**2 * E_uλλ) / 2
           + (E4 * u**3 + 3 * λ * E_uuuλ * u**2 + 3 * λ**2 * E_uuλλ * u
              + λ**3 * E_uλλλ) / 6)
    return out.subs(__mixed_derivatives).expand()

def create_E_uu(u, λ):
    out = (E2 + E3 * u + λ * E_uuλ
           + (E4 * u**2 + 2 * λ * E_uuuλ * u + λ**2 * E_uuλλ) / 2)
    return out.subs(__mixed_derivatives).expand()

def create_u_star(λ):
    return (λ * u0.dot + λ**2 * u0_ddot / 2 + λ**3 * u0_dddot / 6
             + λ**4 * u0_ddddot / 24).expand()

And these functions can be tested against the expressions found above.

assert E == lsk.energy.create_E(u, λ)

assert E_u == lsk.energy.create_E_u(u, λ)

assert E_uu == lsk.energy.create_E_uu(u, λ)
# 5 Bifurcation equations

## 5.1 Introduction

In this chapter, the bifurcation analysis of the perfect system is performed symbolically. It is assumed that through the critical point passes a second equilibrium curve $\lambda \mapsto u(\lambda)$, besides the fundamental branch $\lambda \mapsto u^*(\lambda)$. We seek an asymptotic expansion of $u(\lambda)$ for $\lambda \to \lambda_0$ (critical load). It will be convenient to introduce an auxiliary parametrization $\eta$ such that the bifurcated branch is defined as the set of points $(u(\eta), \lambda(\eta))$. The functions $\eta \mapsto \lambda(\eta)$ and $\eta \mapsto u(\eta)$ are expanded as follows

$$ \lambda(\eta) = \lambda^{(0)} + \eta \lambda^{(1)} + \frac{1}{2} \eta^2 \lambda^{(2)} + \dots \qquad (5.1) $$

and

$$ u(\eta) = u^*[\lambda(\eta)] + \eta u^{(1)} + \frac{1}{2} \eta^2 u^{(2)} + \dots \qquad (5.2) $$

It will be shown below that

$$ u^{(1)} = \xi_i^{(1)} v_i \quad \text{and} \quad u^{(2)} = \xi_i^{(2)} v_i + \xi_i^{(1)} \xi_j^{(1)} w_{ij} + 2\lambda^{(1)} \xi_i^{(1)} w_{i\lambda}, \qquad (5.3) $$

where the coefficients $\lambda^{(1)}, \lambda^{(2)}, \xi_i^{(1)}$ and $\xi_i^{(2)}$ solve the first bifurcation equation

$$ \frac{1}{2} E_{ijk} \xi_j^{(1)} \xi_k^{(1)} + \lambda^{(1)} \dot{E}_{ij} \xi_j^{(1)} = 0 \qquad (5.4) $$

and the second bifurcation equation

$$ \frac{1}{3} E_{ijkl} \xi_j^{(1)} \xi_k^{(1)} \xi_l^{(1)} + \lambda^{(1)} (\dot{E}_{ijk} \xi_k^{(1)} + \lambda^{(1)} \ddot{E}_{ij}) \xi_j^{(1)} + (E_{ijk} \xi_k^{(1)} + \lambda^{(1)} \dot{E}_{ij}) \xi_j^{(2)} + \lambda^{(2)} \dot{E}_{ij} \xi_j^{(1)} = 0. \qquad (5.5) $$

The tensors $\dot{E}_{ij}$, $\ddot{E}_{ij}$, $E_{ijk}$, $\dot{E}_{ijk}$ and $E_{ijkl}$ have been defined in Chapter 2 (see 2.4).

The starting point is the symbolic expression of the residual $(u, \lambda) \mapsto \mathcal{E}_{,u}(u, \lambda)$ that was derived in Chapter 4. We plug the postulated expansions (5.1) and (5.2) into the equilibrium equation

$$ \mathcal{E}_{,u}[u(\eta), \lambda(\eta); \hat{u}] = 0 \quad \text{for all} \quad \hat{u} \in U. $$
The coefficients of $\eta^0, \eta^1$, etc deliver a series of variational problems from which $u^{(k)}$ and $\lambda^{(k)}$ are identified.

```python
from sympy import *
from lsk.display import *
from lsk.energy import *
from lsk symbols import *
```

The asymptotic expansion of $\lambda$ according to Eq. (5.1) is first postulated and plugged into the expression of $u^*$. The resulting symbolic expressions are combined to define the asymptotic expansion of $u$ according to Eq. (5.2).

$$ \lambda = \eta * \lambda^1 + \eta^{**2} / 2 * \lambda^2 + \eta^{**3} / 6 * \lambda^3 + \eta^{**4} / 24 * \lambda^4 + O(\eta^{**5}) $$

```python
u_star = create_u_star(\lambda)
```

$$ u = u\_star + \eta * u1 + \eta^{**2} / 2 * u2 + \eta^{**3} / 6 * u3 + \eta^{**4} / 24 * u4 $$

```python
display_latex_equation(r"\lambda(\eta)", \lambda)
```

```python
display_latex_long_equation(r"u\_\lambda(\eta)", u_star, terms_per_line=6)
```

```python
display_latex_long_equation(r"u\_\eta(\eta)", u, terms_per_line=6)
```

$$ \lambda(\eta) = \eta\lambda^{(1)} + \frac{\eta^2\lambda^{(2)}}{2} + \frac{\eta^3\lambda^{(3)}}{6} + \frac{\eta^4\lambda^{(4)}}{24} + O(\eta^5) $$

$$ u^*(\eta) = \frac{\ddot{u}_0\eta^3\lambda^{(1)3}}{6} + \frac{\dot{u}_0\eta^4\lambda^{(2)2}}{8} + \frac{\ddot{u}_0\eta^2\lambda^{(1)2}}{2} + \frac{\dot{u}_0\eta^3\lambda^{(3)}}{6} + \frac{\ddot{u}_0\eta^2\lambda^{(2)}}{2} + \dot{u}_0\eta\lambda^{(1)} \\ + \frac{\dot{u}_0\eta^4\lambda^{(4)}}{24} + \frac{\ddot{u}_0\eta^3\lambda^{(1)\lambda(2)}}{2} + \frac{\ddot{u}_0\eta^4\lambda^{(1)\lambda(3)}}{6} + \frac{\ddot{u}_0\eta^4\lambda^{(1)2}\lambda^{(2)}}{4} + \frac{\ddot{u}_0\eta^4\lambda^{(1)4}}{24} + O(\eta^5) $$

$$ u(\eta) = \frac{\ddot{u}_0\eta^2\lambda^{(1)2}}{2} + \dot{u}_0\eta\lambda^{(1)} + \frac{\eta^4 u^{(4)}}{24} + \frac{\eta^3 u^{(3)}}{6} + \frac{\eta^2 u^{(2)}}{2} + \eta u^{(1)} \\ + \frac{\ddot{u}_0\eta^4\lambda^{(1)4}}{24} + \frac{\ddot{u}_0\eta^3\lambda^{(1)3}}{6} + \frac{\dot{u}_0\eta^4\lambda^{(2)2}}{8} + \frac{\dot{u}_0\eta^4\lambda^{(4)}}{24} + \frac{\dot{u}_0\eta^3\lambda^{(3)}}{6} + \frac{\dot{u}_0\eta^2\lambda^{(2)}}{2} \\ + \frac{\ddot{u}_0\eta^3\lambda^{(1)\lambda(2)}}{2} + \frac{\ddot{u}_0\eta^4\lambda^{(1)\lambda(3)}}{6} + \frac{\ddot{u}_0\eta^4\lambda^{(1)2}\lambda^{(2)}}{4} + O(\eta^5) $$

These expressions are then used to compute an asymptotic expansion of the residual $\mathcal{E}_{,u}$ along the bifurcated branch.

```python
E_u = (create_E_u(u, \lambda) * u\_\hat{u}).expand()
```

The general form of this expansion is

$$ \mathcal{E}_{,u}[u(\eta), \lambda(\eta); \hat{u}] = \mathcal{E}_{,u}^{(0)}(\hat{u}) + \eta \mathcal{E}_{,u}^{(1)}(u^{(1)}; \hat{u}) + \frac{1}{2} \eta^2 \mathcal{E}_{,u}^{(2)}(u^{(1)}, u^{(2)}, \lambda^{(1)}; \hat{u}) \\ + \frac{1}{6} \eta^3 \mathcal{E}_{,u}^{(3)}(u^{(1)}, u^{(2)}, u^{(3)}, \lambda^{(1)}, \lambda^{(2)}; \hat{u}) + \dots, $$
which delivers the following variational problems

$$
\mathcal{E}_{,u}^{(k)}(u^{(1)}, \dots, u^{(k)}, \lambda^{(1)}, \dots, \lambda^{(k-1)}; \hat{u}) \quad \text{for all } \hat{u} \in U, \quad k = 0, 1, 2, \dots
$$

These problems are studied successively in the following sections. Note that the variational problem of order 0 is in fact uninformative, since $\mathcal{E}_{,u}^{(0)} = 0$.

assert E_u市場.coeff(η, 0) == 0

5.2 The variational problem of order 1

This problem reads

$$
\mathcal{E}_2(u^{(1)}, \hat{u}) = 0 \quad \text{for all} \quad \hat{u} \in U.
$$

assert E_u市场市場.coeff(η, 1) == E2 * u1 * u hat

Therefore $u^{(1)} \in V$ and we introduce the following decomposition

$$
u^{(1)} = \xi_i^{(1)} v_i, \tag{5.6}
$$

where $\xi_1^{(1)}, \dots, \xi_m^{(1)}$ are yet unknown scalars.

5.3 The variational problem of order 2

E_u2 = E_u市场.coeff(η, 2)

The term in $\eta^2$ of the residual delivers the following variational problem

display_latex_equation(E_u2, 0)

$$
\frac{\mathcal{E}_2 \hat{u} u^{(2)}}{2} + \frac{\mathcal{E}_3 \hat{u} u^{(1)^2}}{2} + \dot{\mathcal{E}}_2 \hat{u} \lambda^{(1)} u^{(1)} = 0
$$

for all $\hat{u} \in U$. Testing with $\hat{v} \in V$, the term in $\mathcal{E}_2(u^{(2)}, \hat{v})$ vanishes and we get the following variational problem

lhs1 = E_u2.subs(u hat, v hat).subs(E2 * v hat, 0)
display_latex_equation(lhs1, 0)

$$ \frac{\mathcal{E}_3 \hat{v} u^{(1)}^2}{2} + \dot{\mathcal{E}}_2 \hat{v} \lambda^{(1)} u^{(1)} = 0 $$

to be understood as

$$ \frac{1}{2} \mathcal{E}_3(u^{(1)}, u^{(1)}, \hat{v}) + \lambda^{(1)} \dot{\mathcal{E}}_2(u^{(1)}, \hat{v}) = 0, \qquad (5.7) $$

for all $\hat{v} \in V$. The above equation fully defines $u^{(1)}$. Indeed, plugging the decomposition (5.6) delivers the equivalent equations

$$ \frac{1}{2} E_{ijk} \xi_j^{(1)} \xi_k^{(1)} + \lambda^{(1)} \dot{E}_{ij} \xi_j^{(1)} = 0, $$

where we have introduced $\dot{E}_{ij}$ and $E_{ijk}$, defined by Eqs. (2.6) and (2.8), respectively. We finally retrieve the first bifurcation equation (5.4).

We now test the same equation with $\hat{w} \in W$, plugging the expansion (5.6) of $u^{(1)}$

$$ \frac{1}{2} \mathcal{E}_2(u^{(2)}, \hat{w}) + \frac{1}{2} \xi_i^{(1)} \xi_j^{(1)} \mathcal{E}_3(v_i, v_j, \hat{w}) + \lambda^{(1)} \xi_i^{(1)} \dot{E}_2(v_i, \hat{w}) = 0. $$

The second-order term $u^{(2)}$ is projected onto $V$ and $W$

$$ u^{(2)} = u_V^{(2)} + u_W^{(2)}, \quad \text{where} \quad u_V^{(2)} = \xi_i^{(2)} v_i \in V \quad \text{and} \quad u_W^{(2)} \in W. $$

Plugging this decomposition delivers

$$ \frac{1}{2} \mathcal{E}_2(u_W^{(2)}, \hat{w}) + \frac{1}{2} \xi_i^{(1)} \xi_j^{(1)} \mathcal{E}_3(v_i, v_j, \hat{w}) + \lambda^{(1)} \xi_i^{(1)} \dot{E}_2(v_i, \hat{w}) = 0. $$

The unique solution to this variational problem can be expressed as a function of $w_{i\lambda}$ and $w_{ij}$, (see 2.3)

$$ u_W^{(2)} = \xi_i^{(1)} \xi_j^{(1)} w_{ij} + 2\lambda^{(1)} \xi_i^{(1)} w_{i\lambda}, $$

which proves the decomposition (5.3) of $u^{(2)}$.
## 5.4 The variational problem of order 3

We now turn to the third-order term of the residual $\mathcal{E}_{,u}$. It involves $u^{(3)}$, that gets eliminated when testing with $v_i \in V$.

```
lhs2 = E_u市場. coefficient(η, 3).subs(u hat, v[i]).subs(E2 * v[i], 0)
```

We get the following equations for $i = 1, ..., m$

```
display_latex_equation(lhs2, 0)
```

$$ \frac{\mathcal{E}_3 u^{(1)} u^{(2)} v_i}{2} + \frac{\mathcal{E}_4 u^{(1)} v_i^3}{6} + \frac{\dot{\mathcal{E}}_2 \lambda^{(1)} u^{(1)} v_i}{2} + \frac{\dot{\mathcal{E}}_2 \lambda^{(1)} u^{(2)} v_i}{2} + \frac{\dot{\mathcal{E}}_2 u^{(1)} \lambda^{(2)} v_i}{2} + \frac{\dot{\mathcal{E}}_3 \lambda^{(1)} u^{(1)} v_i^2}{2} = 0 $$

For the terms that involve only $v_i$ and $u^{(1)}$, we use Eqs. (2.6), (2.7) and (2.9) together with the decomposition $u^{(1)} = \xi_i^{(1)} v_i$.

```
d = {
    E2.dot * u1 * v[i]: ξ1[j] * E dot[i, j],
    E2_ddot * u1 * v[i]: (ξ1[j] * (E_ddot[i, j]
        - E2.dot * v[i] * w[j]))
        - E2.dot * u1 * w[i]),
    E3.dot * u1 * u1 * v[i]: (ξ1[j] * ξ1[k] * (E dot[i, j, k]
        - E2.dot * v[i] * w[j, k]))
        - 2 * ξ1[j] * E2.dot * u1 * w[i, j]),
}
lhs2 = rhs2.subs(d).expand()
display_latex_equation(lhs2, 0)
```

$$ \frac{\mathcal{E}_3 u^{(1)} u^{(2)} v_i}{2} + \frac{\mathcal{E}_4 u^{(1)} v_i^3}{6} - \frac{\dot{\mathcal{E}}_2 \lambda^{(1)} u^{(1)} w_i}{2} - \frac{\dot{\mathcal{E}}_2 \lambda^{(1)} v_i w_j \xi_j^{(1)}}{2} - \dot{\mathcal{E}}_2 \lambda^{(1)} u^{(1)} w_{i,j} \xi_j^{(1)} + \frac{\dot{\mathcal{E}}_2 \lambda^{(1)} u^{(2)} v_i}{2} - \frac{\dot{\mathcal{E}}_2 \lambda^{(1)} v_i w_{j,k} \xi_j^{(1)} \xi_k^{(1)}}{2} + \frac{\lambda^{(1)} \dot{\mathcal{E}}_{i,j} \xi_j^{(1)}}{2} + \frac{\lambda^{(1)} \dot{\mathcal{E}}_{i,j,k} \xi_j^{(1)} \xi_k^{(1)}}{2} + \frac{\lambda^{(2)} \dot{\mathcal{E}}_{i,j} \xi_j^{(1)}}{2} = 0 $$

In the remainder of this section, we simplify this equation. To do so, we rely heavily on Eqs. (2.5) and (2.4) as well as the definitions of 2.4. We first apply the following simplification

$$ \dot{\mathcal{E}}_2(u^{(2)}, v_i) = \dot{E}_{ij} \xi_j^{(2)} + \xi_j^{(1)} \xi_k^{(1)} \dot{\mathcal{E}}_2(w_{jk}, v_i) + 2\lambda^{(1)} \dot{\mathcal{E}}_2(w_i, u^{(1)}). $$
Proof

$$
\begin{align*}
\dot{\mathscr{E}}_2(u^{(2)}, v_i) &= \dot{E}_{ij} \xi_j^{(2)} + \xi_j^{(1)} \xi_k^{(1)} \dot{\mathscr{E}}_2(w_{jk}, v_i) + 2\lambda^{(1)} \xi_j^{(1)} \dot{\mathscr{E}}_2(w_j, v_i) \\
&= \dot{E}_{ij} \xi_j^{(2)} + \xi_j^{(1)} \xi_k^{(1)} \dot{\mathscr{E}}_2(w_{jk}, v_i) - 2\lambda^{(1)} \xi_j^{(1)} \mathscr{E}_2(w_i, w_j) \\
&= \dot{E}_{ij} \xi_j^{(2)} + \xi_j^{(1)} \xi_k^{(1)} \dot{\mathscr{E}}_2(w_{jk}, v_i) + 2\lambda^{(1)} \xi_j^{(1)} \dot{\mathscr{E}}_2(w_i, v_j) \\
&= \dot{E}_{ij} \xi_j^{(2)} + \xi_j^{(1)} \xi_k^{(1)} \dot{\mathscr{E}}_2(w_{jk}, v_i) + 2\lambda^{(1)} \dot{\mathscr{E}}_2(w_i, u^{(1)}).
\end{align*}
$$

$$
\text{expr} = (\text{E\_dot}[i, j] * \xi_2[j] + \xi_1[j] * \xi_1[k] * \text{E2\_dot} * \text{v}[i] * \text{w}[j, k] \\
+ 2 * \lambda_1 * \text{E2\_dot} * \text{u1} * \text{w}[i])
$$

$$
lhs2 = \text{lhs2.subs}(\text{E2\_dot} * \text{u2} * \text{v}[i], \text{expr}).\text{expand()}
$$

$$
\text{display\_latex\_equation}(\text{lhs2}, \emptyset)
$$

$$
\frac{\mathscr{E}_3 u^{(1)} u^{(2)} v_i}{2} + \frac{\mathscr{E}_4 u^{(1)} v_i}{6} + \frac{\dot{\mathscr{E}}_2 \lambda^{(1)}^2 u^{(1)} w_i}{2} - \frac{\dot{\mathscr{E}}_2 \lambda^{(1)}^2 v_i w_j \xi_j^{(1)}}{2} - \dot{\mathscr{E}}_2 \lambda^{(1)} u^{(1)} w_{i,j} \xi_j^{(1)} + \frac{\lambda^{(1)}^2 \ddot{E}_{i,j} \xi_j^{(1)}}{2} + \\
\frac{\lambda^{(1)} \dot{E}_{i,j,k} \xi_j^{(1)} \xi_k^{(1)}}{2} + \frac{\lambda^{(1)} \dot{E}_{i,j} \xi_j^{(2)}}{2} + \frac{\lambda^{(2)} \dot{E}_{i,j} \xi_j^{(1)}}{2} = 0
$$

Then,

$$
\mathscr{E}_3(u^{(1)}, u^{(2)}, v_i) = E_{ijk} \xi_j^{(2)} \xi_k^{(1)} + \frac{1}{3} \xi_j^{(1)} \xi_k^{(1)} \xi_l^{(1)} E_{ijkl} \\
- \frac{1}{3} \mathscr{E}_4(u^{(1)}, u^{(1)}, u^{(1)}, v_i) + 2\lambda^{(1)} \xi_j^{(1)} \dot{\mathscr{E}}_2(u^{(1)}, w_{ij}).
$$

Proof

$$
\begin{align*}
\mathscr{E}_3(u^{(1)}, u^{(2)}, v_i) &= \xi_j^{(2)} \xi_k^{(1)} \mathscr{E}_3(v_i, v_j, v_k) + \xi_j^{(1)} \xi_k^{(1)} \xi_l^{(1)} \mathscr{E}_3(v_i, v_j, w_{kl}) \\
&\quad + 2\lambda^{(1)} \xi_j^{(2)} \xi_k^{(1)} \mathscr{E}_3(v_i, v_j, w_k) \\
&= E_{ijk} \xi_j^{(2)} \xi_k^{(1)} \\
&\quad + \frac{1}{3} \xi_j^{(1)} \xi_k^{(1)} \xi_l^{(1)} \left[ \mathscr{E}_3(v_i, v_j, w_{kl}) + \mathscr{E}_3(v_i, v_k, w_{jl}) + \mathscr{E}_3(v_i, v_l, w_{jk}) \right] \\
&\quad - 2\lambda^{(1)} \xi_j^{(1)} \xi_k^{(1)} \mathscr{E}_2(w_{ij}, w_k) \\
&= E_{ijk} \xi_j^{(2)} \xi_k^{(1)} + \frac{1}{3} \xi_j^{(1)} \xi_k^{(1)} \xi_l^{(1)} \left[ E_{ijkl} - \mathscr{E}_4(v_i, v_j, v_k, v_l) \right] \\
&\quad + 2\lambda^{(1)} \xi_j^{(1)} \xi_k^{(1)} \dot{\mathscr{E}}_2(w_{ij}, v_k)
\end{align*}
$$

and the above identity is retrieved.
expr = (E[i, j, k] * ξ2[j] * ξ1[k] + E[i, j, k, l] * ξ1[j] * ξ1[k] * ξ1[l] / 3 - E4 * v[i] * u1**3 / 3 + 2 * λ1 * ξ1[j] * E2 dot * u1 * w[i, j])

lhs2 =lhs2.subs(E3 * u1 * u2 * v[i], expr).expand()

display_latex_equation(lhs2, 0)

$$
\frac{\dot{E}_2 \lambda^{(1)} u^{(1)} w_i}{\frac{2}{6}} - \frac{\dot{E}_2 \lambda^{(1)} v_i w_j \xi_j^{(1)}}{\frac{2}{2}} + \frac{\lambda^{(1)} \ddot{E}_{i,j} \xi_j^{(1)}}{2} + \frac{\lambda^{(1)} \dot{E}_{i,j,k} \xi_j^{(1)} \xi_k^{(1)}}{2} + \frac{\lambda^{(1)} \dot{E}_{i,j} \xi_j^{(2)}}{2} + \frac{\lambda^{(2)} \dot{E}_{i,j} \xi_j^{(1)}}{2} + \frac{\lambda^{(2)} \dot{E}_{i,j} \xi_j^{(1)}}{2} + \frac{\lambda^{(2)} \dot{E}_{i,j,k} \xi_k^{(1)} \xi_l^{(1)}}{2} = 0
$$

Finally,

$$
\xi_j^{(1)} \dot{E}_2(w_j, v_i) = -\xi_j^{(1)} E_2(w_i, w_j) = \xi_j^{(1)} \dot{E}_2(v_j, w_i) = \dot{E}_2(u^{(1)}, w_i).
$$

lhs2 =lhs2.subs(ξ1[j] * E2 dot * v[i] * w[j], E2 dot * u1 * w[i])

display_latex_equation(lhs2, 0)

$$
\frac{\lambda^{(1)} \dot{E}_{i,j} \xi_j^{(1)}}{2} + \frac{\lambda^{(1)} \dot{E}_{i,j,k} \xi_j^{(1)} \xi_k^{(1)}}{2} + \frac{\lambda^{(1)} \dot{E}_{i,j} \xi_j^{(2)}}{2} + \frac{\lambda^{(2)} \dot{E}_{i,j} \xi_j^{(1)}}{2} + \frac{E_{i,j,k,l} \xi_j^{(1)} \xi_k^{(1)} \xi_l^{(1)}}{6} + \frac{E_{i,j,k} \xi_k^{(1)} \xi_j^{(2)}}{2} = 0
$$

And the second bifurcation equation (5.5) is retrieved.

expected = (E[i, j, k, l] * ξ1[j] * ξ1[k] * ξ1[l] / 3
+ λ1 * (E Dot[i, j, k] * ξ1[k] + λ1 * E Ddot[i, j]) * ξ1[j]
+ (E[i, j, k] * ξ1[k] + λ1 * E Dot[i, j]) * ξ2[j]
+ λ2 * E Dot[i, j] * ξ1[j])

assert expand(2 *lhs2 - expected) == 0
# 6 Asymptotic expansion of the energy and its hessian

The present chapter is organized as follows. In 6.1, we derive the following asymptotic expansion of the energy along the bifurcated branch

$$ \mathcal{E}[u(\eta), \lambda(\eta)] = \mathcal{E}[u^* \circ \lambda(\eta), \lambda(\eta)] + \frac{1}{6}\eta^3 \lambda^{(1)} \xi_i^{(1)} \xi_j^{(1)} \dot{E}_{ij} \\ + \frac{1}{24}\eta^4 \{E_{ijkl} \xi_i^{(1)} \xi_j^{(1)} \xi_k^{(1)} \xi_l^{(1)} + 4\lambda^{(1)} \dot{E}_{ijk} \xi_i^{(1)} \xi_j^{(1)} \xi_k^{(1)} \\ + 6[(\lambda^{(1)})^2 \ddot{E}_{ij} + \lambda^{(2)} \dot{E}_{ij}] \xi_i^{(1)} \xi_j^{(1)}\} + o(\eta^4). \quad (6.1) $$

Then, in 6.2, we derive the asymptotic expansion of the hessian of the energy along the bifurcated branch

$$ \mathcal{E}_{uu}[u(\eta), \lambda(\eta); \bullet, \bullet] = \mathcal{E}_2(\bullet, \bullet) + \eta [\mathcal{E}_3(u^{(1)}, \bullet, \bullet) + \lambda^{(1)} \dot{\mathcal{E}}_2(\bullet, \bullet)] \\ + \frac{1}{2}\eta^2 [\mathcal{E}_4(u^{(1)}, u^{(1)}, \bullet, \bullet) + \mathcal{E}_3(u^{(2)}, \bullet, \bullet) \\ + 2\lambda^{(1)} \dot{\mathcal{E}}_3(u^{(1)}, \bullet, \bullet) + (\lambda^{(1)})^2 \ddot{\mathcal{E}}_2(\bullet, \bullet) + \lambda^{(2)} \dot{\mathcal{E}}_2(\bullet, \bullet)]. \quad (6.2) $$

from sympy import *
from lsk.display import *
from lsk.energy import *
from lsk symbols import *

As in Chapter 5, $\lambda$ and $u^*$ are defined as asymptotic expansions of the powers of $\eta$.

$\lambda = \eta * \lambda_1 + \eta^{**2} / 2 * \lambda_2 + \eta^{**3} / 6 * \lambda_3 + \eta^{**4} / 24 * \lambda_4 + O(\eta^{**5})$
$u\_star = create\_u\_star(\lambda)$

display_latex_equation(r"\lambda\eta(\lambda)\eta", \lambda)
display_latex_long_equation(r"u\_\lambda\eta(\lambda)\eta", u\_star, terms\_per\_line=5)

$$ \lambda(\eta) = \eta\lambda^{(1)} + \frac{\eta^2\lambda^{(2)}}{2} + \frac{\eta^3\lambda^{(3)}}{6} + \frac{\eta^4\lambda^{(4)}}{24} + O(\eta^5) $$
$$
\begin{aligned}
u^*(\eta) = & \frac{\ddot{u}_0 \eta^3 \lambda^{(1)3}}{6} + \frac{\dot{u}_0 \eta^2 \lambda^{(1)2}}{2} + \frac{\dot{u}_0 \eta^3 \lambda^{(3)}}{6} + \frac{\dot{u}_0 \eta^2 \lambda^{(2)}}{2} + \dot{u}_0 \eta \lambda^{(1)} \\
& + \frac{\dddot{u}_0 \eta^4 \lambda^{(1)4}}{24} + \frac{\dddot{u}_0 \eta^4 \lambda^{(1)2} \lambda^{(2)}}{4} + \frac{\dot{u}_0 \eta^4 \lambda^{(2)2}}{8} + \frac{\ddot{u}_0 \eta^3 \lambda^{(1)} \lambda^{(2)}}{2} + \frac{\dot{u}_0 \eta^4 \lambda^{(4)}}{24} \\
& + \frac{\dot{u}_0 \eta^4 \lambda^{(1)} \lambda^{(3)}}{6} + O(\eta^5)
\end{aligned}
$$

The bifurcated branch $u(\eta)$ is also expanded. Moreover, the second-order term, $u^{(2)}$ is expressed as the orthogonal decomposition

$$u^{(2)} = u_{V}^{(2)} + u_{W}^{(2)}, \quad \text{with} \quad u_{V}^{(2)} \in V \quad \text{and} \quad u_{W}^{(2)} \in W.$$

It was in fact shown in Chapter 5 that

$$u_{W}^{(2)} = \xi_{i}^{(1)} \xi_{j}^{(1)} w_{ij} + 2\lambda^{(1)} \xi_{i}^{(1)} w_{i},$$

while $u_V^{(2)}$ is expanded as follows:

$$u_V^{(2)} = \xi_i^{(2)} v_i.$$

u = expand(u_star + \eta * u1 + \eta^{**2} / 2 * (u2_V + u2_W) + \eta^{**3} / 6 * u3 + \eta^{**4} / 24 * u4)

display_latex_long_equation(r"u(\eta)", u, terms_per_line=5)

$$
\begin{aligned}
u(\eta) = & \frac{\eta^4 u^{(4)}}{24} + \frac{\eta^3 u^{(3)}}{6} + \frac{\eta^2 u_V^{(2)}}{2} + \frac{\eta^2 u_W^{(2)}}{2} + \eta u^{(1)} \\
& + \frac{\dddot{u}_0 \eta^3 \lambda^{(1)3}}{6} + \frac{\dddot{u}_0 \eta^2 \lambda^{(1)2}}{2} + \frac{\dot{u}_0 \eta^3 \lambda^{(3)}}{6} + \frac{\dot{u}_0 \eta^2 \lambda^{(2)}}{2} + \dot{u}_0 \eta \lambda^{(1)} \\
& + \frac{\dddot{u}_0 \eta^4 \lambda^{(1)4}}{24} + \frac{\dddot{u}_0 \eta^4 \lambda^{(1)2} \lambda^{(2)}}{4} + \frac{\dot{u}_0 \eta^4 \lambda^{(2)2}}{8} + \frac{\ddot{u}_0 \eta^3 \lambda^{(1)} \lambda^{(2)}}{2} + \frac{\dot{u}_0 \eta^4 \lambda^{(4)}}{24} \\
& + \frac{\dot{u}_0 \eta^4 \lambda^{(1)} \lambda^{(3)}}{6} + O(\eta^5)
\end{aligned}
$$

## 6.1 Asymptotic expansion of the energy

We expand the following quantity at the critical point

$$\Delta\mathcal{E}(\eta) = \mathcal{E}[u(\eta), \lambda(\eta)] - \mathcal{E}[u^* \circ \lambda(\eta), \lambda(\eta)],$$
where $\eta$ is the parametrization of the bifurcated branch introduced in Chapter 5 ($\eta = 0$ at the critical point).

The asymptotic expansion of $\Delta\mathcal{E}$ results from plugging the expansions of $\lambda(\eta)$ and $u(\eta)$ defined above into the general expression of the energy that was derived in Chapter 4.

We then use these expansions to evaluate the energy along the **fundamental** and **bifurcated** branches, as well as the difference $\Delta\mathcal{E}$ of these two quantities. The resulting expressions are too long to be displayed. We first apply some elementary simplifications.

$$ \begin{aligned} \Delta E = & (\text{create}_E(u, \lambda) - \text{create}_E(u_{\text{star}}, \lambda)).\text{expand()}.\text{subs}(\{ \\ & \text{E2} * u1 : 0, \\ & \text{E2} * u2\_V : 0 \}) \end{aligned} $$

```python
assert \Delta E \cdot \text{coeff}(\eta, 0) == 0
assert \Delta E \cdot \text{coeff}(\eta, 1) == 0
assert \Delta E \cdot \text{coeff}(\eta, 2) == 0
```

Simplification of the third- and fourth- order terms is performed below, first observing that, for all $v \in V$

$$ \mathcal{E}_3(u^{(1)}, u^{(1)}, v) = -2\lambda^{(1)} \dot{\mathcal{E}}_2(u^{(1)}, v). \qquad (6.3) $$

### Proof

Let $v = \xi_i v_i \in V$

$$ \begin{aligned} \mathcal{E}_3(u^{(1)}, u^{(1)}, v) &= \xi_i \xi_j^{(1)} \xi_k^{(1)} \mathcal{E}_3(v_i, v_j, v_k) = E_{ijk} \xi_i \xi_j^{(1)} \xi_k^{(1)} \\ &= -2\lambda^{(1)} \dot{E}_{ij} \xi_i \xi_j^{(1)} = -2\lambda^{(1)} \dot{\mathcal{E}}_2(u^{(1)}, v), \end{aligned} $$

where we have used the definitions (2.6) and (2.8) of $\dot{E}_{ij}$ and $E_{ijk}$ and the first bifurcation equation (5.4).

### Note

Note that the first bifurcation equation was used to prove Eq. (6.3).

$$ \begin{aligned} d = & \{ \\ & \text{E3} * u1 * u1 * u1 : -2 * \lambda1 * \text{E2\_dot} * u1 * u1, \\ & \text{E3} * u1 * u1 * u2\_V : -2 * \lambda1 * \text{E2\_dot} * u1 * u2\_V, \end{aligned} $$

$$ \Delta E = \Delta E \cdot \text{subs}(d) \cdot \text{expand()} $$
The subsequent simplifications rely heavily on definitions (2.5) and (2.4) of $w_{ij}$ and $w_i$, definitions (2.6), (2.7), (2.8), (2.9) and (2.10) of $\dot{E}_{ij}$, $\ddot{E}_{ij}$, $E_{ijk}$, $\dot{E}_{ijk}$ and $E_{ijkl}$. Finally, expressions (5.3) of $u^{(1)}$ and $u^{(2)}$ are used.

$$ 
\begin{aligned}
\mathcal{E}_2(u_W^{(2)}, u_W^{(2)}) &= \xi_i^{(1)} \xi_j^{(1)} \xi_k^{(1)} \xi_l^{(1)} \mathcal{E}_2(w_{ij}, w_{kl}) + 4 \lambda^{(1)} \xi_i^{(1)} \xi_j^{(1)} \xi_k^{(1)} \mathcal{E}_2(w_{ij}, w_k) \\
&\quad + 4 (\lambda^{(1)})^2 \xi_i^{(1)} \xi_j^{(1)} \mathcal{E}_2(w_i, w_j) \\
&= \frac{1}{3} [\mathcal{E}_4(u^{(1)}, u^{(1)}, u^{(1)}, u^{(1)}) - E_{ijkl} \xi_i^{(1)} \xi_j^{(1)} \xi_k^{(1)} \xi_l^{(1)}] \\
&\quad + \frac{4}{3} \lambda^{(1)} [\dot{\mathcal{E}}_3(u^{(1)}, u^{(1)}, u^{(1)}) - \dot{E}_{ijk} \xi_i^{(1)} \xi_j^{(1)} \xi_k^{(1)}] \\
&\quad + 2 (\lambda^{(1)})^2 [\ddot{\mathcal{E}}_2(u^{(1)}, u^{(1)}) - \ddot{E}_{ij} \xi_i^{(1)} \xi_j^{(1)}]
\end{aligned}
\quad (6.4) $$

$$ 
\begin{aligned}
\text{expr} = & ((E4 * u1 * u1 * u1 * u1 \\
& - E[i, j, k, l] * \xi1[i] * \xi1[j] * \xi1[k] * \xi1[l]) / 3 \\
& + 4 * \lambda1 / 3 * (E3\_dot * u1 * u1 * u1 \\
& - E\_dot[i, j, k] * \xi1[i] * \xi1[j] * \xi1[k]) \\
& + 2 * \lambda1**2 * (E2\_ddot * u1 * u1 \\
& - E\_ddot[i, j] * \xi1[i] * \xi1[j]))
\end{aligned}
$$

$$ \Delta E = \Delta E.\text{subs}(E2 * u2_W * u2_W, \text{expr}).\text{expand()} $$

$$ 
\begin{aligned}
\dot{\mathcal{E}}_2(u^{(1)}, u_W^{(2)}) = & \frac{1}{3} [\dot{E}_{ijk} \xi_i^{(1)} \xi_j^{(1)} \xi_k^{(1)} - \dot{\mathcal{E}}_3(u^{(1)}, u^{(1)}, u^{(1)})] \\
& + \lambda^{(1)} [\ddot{E}_{ij} \xi_i^{(1)} \xi_j^{(1)} - \ddot{\mathcal{E}}_2(u^{(1)}, u^{(1)})]
\end{aligned}
$$

<figure><img src="image_3.png" /></figure>

Proof

$$ 
\begin{aligned}
\dot{\mathcal{E}}_2(u^{(1)}, u_W^{(2)}) &= \xi_i^{(1)} \xi_j^{(1)} \xi_k^{(1)} \dot{\mathcal{E}}_2(v_i, w_{jk}) + 2\lambda^{(1)} \xi_i^{(1)} \xi_j^{(1)} \dot{\mathcal{E}}_2(v_i, w_{j\lambda}) \\
&= \frac{1}{3} \xi_i^{(1)} \xi_j^{(1)} \xi_k^{(1)} [\dot{\mathcal{E}}_2(v_i, w_{jk}) + \dot{\mathcal{E}}_2(v_j, w_{ik}) + \dot{\mathcal{E}}_2(v_k, w_{ij})] \\
&\quad + \lambda^{(1)} \xi_i^{(1)} \xi_j^{(1)} [\dot{\mathcal{E}}_2(v_i, w_{j\lambda}) + \dot{\mathcal{E}}_2(v_j, w_{i\lambda})] \\
&= \frac{1}{3} \xi_i^{(1)} \xi_j^{(1)} \xi_k^{(1)} [\dot{E}_{ijk} - \dot{\mathcal{E}}_3(v_i, v_j, v_k)] \\
&\quad + \lambda^{(1)} \xi_i^{(1)} \xi_j^{(1)} [\ddot{E}_{ij} - \ddot{\mathcal{E}}_2(v_i, v_j)]
\end{aligned}
$$

$$ 
\begin{aligned}
\text{expr} = & ((E\_dot[i, j, k] * \xi1[i] * \xi1[j] * \xi1[k] \\
& - E3\_dot * u1 * u1 * u1) / 3 \\
& + \lambda1 * (E\_ddot[i, j] * \xi1[i] * \xi1[j] \\
& - E2\_ddot * u1 * u1))
\end{aligned}
$$

$$ \Delta E = \Delta E.\text{subs}(E2\_dot * u1 * u2\_W, \text{expr}).\text{expand()} $$
$$ \Delta \mathcal{E}^{(3)} = \dot{\mathcal{E}}_2 \lambda^{(1)} u^{(1)} $$

$$ \Delta \mathcal{E}^{(4)} = 6 \dot{\mathcal{E}}_2 u^{(1)^2} \lambda^{(2)} + 6 \lambda^{(1)^2} \ddot{E}_{i,j} \xi^{(1)}_i \xi^{(1)}_j + 4 \lambda^{(1)} \dot{E}_{i,j,k} \xi^{(1)}_i \xi^{(1)}_j \xi^{(1)}_k + E_{i,j,k,l} \xi^{(1)}_i \xi^{(1)}_j \xi^{(1)}_k \xi^{(1)}_l $$

## 6.2 Asymptotic expansion of the hessian of the energy

We now turn to the hessian of the energy, which is expanded to second order in $\eta$

$$ \mathcal{E}_{uu}[u(\eta), \lambda(\eta)] = \mathcal{H}^{(0)} + \eta \mathcal{H}^{(1)} + \frac{1}{2} \eta^2 \mathcal{H}^{(2)} $$

E_uu = create_E_uu(u, \lambda).expand() + O(\eta^{3})

H0 = E_uuolscoff(η, 0)

H1 = E_uuolscoff(η, 1)

H2 = expand(2 * E_uuolscoff(η, 2))

d = {
  r"\order[0]{\mathcal{H}}": H0,
  r"\order[1]{\mathcal{H}}": H1,
  r"\order[2]{\mathcal{H}}": H2,
}

display_latex_dict(d, num_COLS=1)

$$ \mathcal{H}^{(0)} = \mathcal{E}_2 $$

$$ \mathcal{H}^{(1)} = \mathcal{E}_3 u^{(1)} + \dot{\mathcal{E}}_2 \lambda^{(1)} $$

$$ \mathcal{H}^{(2)} = \mathcal{E}_3 u_V^{(2)} + \mathcal{E}_3 u_W^{(2)} + \mathcal{E}_4 u^{(1)^2} + \ddot{\mathcal{E}}_2 \lambda^{(1)^2} + \dot{\mathcal{E}}_2 \lambda^{(2)} + 2 \dot{\mathcal{E}}_3 \lambda^{(1)} u^{(1)} $$

assert H0 == E2

assert H1 == E3 * u1 + \lambda1 * E2.dot

expected = (E4 * u1 * u1

+ E3 * u2_V

+ E3 * u2_W

+ 2 * \lambda1 * E3.dot * u1

+ \lambda1**2 * E2.ddot

+ \lambda2 * E2.dot)

assert H2 == expected

Using Eq. (5.3), the above expressions can be expanded as follows

$$ \mathcal{H}^{(0)} = \mathcal{E}_2, \quad \mathcal{H}^{(1)} = \xi_i^{(1)} \mathcal{E}_3(v_i, \bullet, \bullet) + \lambda^{(1)} \dot{\mathcal{E}}_2, $$
and

$$ \mathcal{H}^{(2)} = \xi_i^{(1)} \xi_j^{(1)} \mathcal{E}_4(v_i, v_j, \bullet, \bullet) + \xi_i^{(2)} \mathcal{E}_3(v_i, \bullet, \bullet) \\ + \xi_i^{(1)} \xi_j^{(1)} \mathcal{E}_3(w_{ij}, \bullet, \bullet) + 2\lambda^{(1)} \xi_i^{(1)} \mathcal{E}_3(w_i, \bullet, \bullet) \\ + 2\lambda^{(1)} \xi_i^{(1)} \dot{\mathcal{E}}_3(v_i, \bullet, \bullet) + (\lambda^{(1)})^2 \ddot{\mathcal{E}}_2(\bullet, \bullet) + \lambda^{(2)} \dot{\mathcal{E}}_2(\bullet, \bullet), $$

and Eq. (6.2) is retrieved.
# 7 Eigenmodes of the hessian of the energy

In view of stability analysis, the eigenvalues $\alpha$ and eigenvectors $x$ of the hessian of the energy are expanded in this chapter to second order in $\eta$

```python
from sympy import *
from lsk.display import *
from lsk symbols import *
from lsk.energy import *

x0_V, x0_W = symbols(r"{\order[0]{x}_V} {\order[0]{x}_W}")
x1_V, x1_W = symbols(r"{\order[1]{x}_V} {\order[1]{x}_W}")
x2_V, x2_W = symbols(r"{\order[2]{x}_V} {\order[2]{x}_W}")
a0, a1, a2 = symbols(r"{\order[0]{\alpha}} {\order[1]{\alpha}} {\order[2]{\alpha}}")

\chi_0 = IndexedBase(r"\order[0]{\chi}")
\chi_1 = IndexedBase(r"\order[1]{\chi}")

x = (x0_V + x0_W) + \eta * (x1_V + x1_W) + \eta^{2/2} * (x2_V + x2_W) + O(\eta^{3})
\alpha = a0 + \eta * a1 + \eta^{2/2} * a2 + O(\eta^{3})

display_latex_dict({"x": x, "\alpha": \alpha}, num_COLS=1)
```

$$x = x_{W}^{(0)} + x_{V}^{(0)} + \eta (x_{V}^{(1)} + x_{W}^{(1)}) + \frac{\eta^{2} (x_{V}^{(2)} + x_{W}^{(2)})}{2} + O(\eta^{3})$$

$$\alpha = \alpha^{(0)} + \eta\alpha^{(1)} + \frac{\eta^2\alpha^{(2)}}{2} + O(\eta^3)$$

Note that the eigenvector $x$ is projected onto $V$ and $W$: $x_V^{(k)} \in V$ and $x_W^{(k)} \in W$. It will be convenient to expand the $V$-component in the $(v_1, \dots, v_m)$ basis

$$x_V^{(k)} = \chi_i^{(k)} v_i.$$

We immediately have the following simplification rules

```python
rules = {
    E2 * v[i]: 0,
    E2 * x0_V: 0,
    E2 * x1_V: 0,
    E2 * x2_V: 0
}
```
We focus in this chapter on potentially unstable eigenmodes, for which the eigenvalue $\alpha$ might be negative in the vicinity of $\eta = 0$. It will be shown that these eigenmodes are necessarily such that $\alpha^{(0)} = 0$. In that case, $x_W^{(0)} = 0$ and

$$x_W^{(1)} = \lambda^{(1)} \chi_j^{(0)} w_j + \chi_j^{(0)} \xi_k^{(1)} w_{jk}. \qquad (7.1)$$

The coefficients $\alpha^{(1)}, \alpha^{(2)}, \chi_i^{(0)}$ and $\chi_i^{(1)}$ solve the following equations

$$(E_{ijk} \xi_k^{(1)} + \lambda^{(1)} \dot{E}_{ij}) \chi_j^{(0)} = \alpha^{(1)} \chi_i^{(0)}, \qquad (7.2)$$

and

$$[E_{ijkl} \xi_k^{(1)} \xi_l^{(1)} + \lambda^{(1)} (2\dot{E}_{ijk} \xi_k^{(1)} + \lambda^{(1)} \ddot{E}_{ij}) + E_{ijk} \xi_k^{(2)} + \lambda^{(2)} \dot{E}_{ij}] \chi_j^{(0)} + 2(E_{ijk} \xi_k^{(1)} + \lambda^{(1)} \dot{E}_{ij}) \chi_j^{(1)} = 2\alpha^{(1)} \chi_i^{(1)} + \alpha^{(2)} \chi_i^{(0)}. \qquad (7.3)$$

## 7.1 The eigenvalue problem

The eigenvalues $\alpha \in \mathbb{R}$ and eigenvectors $x \in U$ of the hessian are such that

$$\mathcal{E}_{uu}[u(\eta), \lambda(\eta); x, \hat{u}] = \alpha \langle x, \hat{u} \rangle \quad \text{for all } \hat{u} \in U, \qquad (7.4)$$

where $\eta \mapsto \lambda(\eta)$ and $\eta \mapsto u(\eta)$ define the bifurcated branch.

In Chapter 5 and in Chapter 6, the following asymptotic expansions of $\lambda(\eta)$, $u(\eta)$ and $\mathcal{E}_{uu}[u(\eta), \lambda(\eta); \bullet, \bullet]$ were derived

H0 = E2
H1 = E3 * u1 + \lambda1 * E2 Dot
H2 = (E4 * u1 * u1 + E3 * (u2_V + u2_W) + 2 * \lambda1 * E3 Dot * u1 + \lambda1**2 * E2 Ddot + \lambda2 * E2 Dot)

H = expand(H0 + \eta * H1 + \eta**2 / 2 * H2 + O(\eta**3))

display_latex_long_equation(r"\E_{uu}[u(\eta), \lambda(\eta)]", H, terms_per_line=6)

$$\mathcal{E}_{uu}[u(\eta), \lambda(\eta)] = \mathcal{E}_2 + \frac{\mathcal{E}_3 \eta^2 u_V^{(2)}}{2} + \frac{\mathcal{E}_3 \eta^2 u_W^{(2)}}{2} + \mathcal{E}_3 \eta u^{(1)} + \frac{\mathcal{E}_4 \eta^2 u^{(1)}^2}{2} + \dot{\mathcal{E}}_2 \eta \lambda^{(1)} \\ + \dot{\mathcal{E}}_3 \eta^2 \lambda^{(1)} u^{(1)} + \frac{\dot{\mathcal{E}}_2 \eta^2 \lambda^{(2)}}{2} + \frac{\ddot{\mathcal{E}}_2 \eta^2 \lambda^{(1)}^2}{2} + O(\eta^3)$$
We define the left-hand side and right-hand side of the variational problem that defines the eigenpairs $(\alpha, x)$ at all orders.

lhs = (H * x * uHat).subs rules. expand()

rhs = (α * x * uHat).subs rules. expand()

The terms of order 0, 1 and 2 in $\eta$ are identified below.

## 7.2 The variational problem of order 0

lhs0 =lhs .coeff (η, 0) .subs (rules)

rhs0 = rhs .coeff (η, 0)

The lowest-order problem reads: find $x_V^{(0)} \in V$, $x_W^{(0)} \in W$ and $\alpha^{(0)} \in \mathbb{R}$ such that, for all $\hat{u} \in U$

display_latex_equation (lhs0, rhs0)

$$ \mathcal{E}_2 \hat{u} x_W^{(0)} = \hat{u} \alpha^{(0)} x_V^{(0)} + \hat{u} \alpha^{(0)} x_W^{(0)} $$

to be understood as

$$ \mathcal{E}_2(x_W^{(0)}, \hat{u}) = \alpha^{(0)} \langle x_V^{(0)}, \hat{u} \rangle + \alpha^{(0)} \langle x_W^{(0)}, \hat{u} \rangle. $$

This variational equation is tested with an element $\hat{w} \in W$. Then, $\langle x_W^{(0)}, \hat{w} \rangle$ vanished and we are left with finding $x_W^{(0)} \in W$ such that, for all $\hat{w} \in W$,

$$ \mathcal{E}_2(x_W^{(0)}, \hat{w}) = \alpha^{(0)} \langle x_W^{(0)}, \hat{w} \rangle. $$

If $x_W^{(0)} \neq 0$, then $\alpha^{(0)}$ is an eigenvalue of $\mathcal{E}_2$ over $W$. Therefore, $\alpha^{(0)} > 0$ and $\alpha > 0$ in the vicinity of the critical point $\lambda = \lambda_0$: the resulting eigenmode is stable. So, in order to find potentially unstable modes, we need to consider that $x_W^{(0)} = 0$ and $\alpha^{(0)} = 0$.

rules [α0] = 0

rules [x0_W] = 0

## 7.3 The variational problem of order 1

lhs1 =lhs .coeff (η, 1) .subs (rules)

rhs1 =rhs .coeff (η, 1) .subs (rules)
The problem reads: find $x_V^{(0)} \in V$, $x_W^{(1)} \in W$ and $\alpha^{(1)} \in \mathbb{R}$ such that, for all $\hat{u} \in U$

display_latex_equation(lhs1, rhs1)

$$ \mathcal{E}_2 \hat{u} x_W^{(1)} + \mathcal{E}_3 \hat{u} x_V^{(0)} u^{(1)} + \dot{\mathcal{E}}_2 \hat{u} x_V^{(0)} \lambda^{(1)} = \hat{u} x_V^{(0)} \alpha^{(1)} $$

Testing first with $\hat{u} = v_i$ delivers the following variational problem: find $x_V^{(0)} \in V$ such that, for all $i = 1, \dots, m$

d = {
  uHat: v[i],
  x0V : x0[j] * v[j],
  u1 : x1[k] * v[k]
}
lhs1a = rhs1.subs(d).subs rules
rhs1a = rhs1.subs(d).subs rules.subs(x0[j] * v[j] * v[i], x0[i])
display_latex_equation(lhs1a, rhs1a)

$$ \mathcal{E}_3 \chi^{(0)}_j v_i v_j v_k \xi^{(1)}_k + \dot{\mathcal{E}}_2 \lambda^{(1)} \chi^{(0)}_j v_i v_j = \alpha^{(1)} \chi^{(0)}_i $$

and Eq. (7.2) is retrieved.

The test function is now picked in $W = V^\perp$ and we get the following variational problem: find $x_W^{(1)} \in W$ such that, for all $\hat{w} \in W$,

d = {
  uHat: wHat,
  x0V : x0[i] * v[i],
  u1 : x1[j] * v[j]
}
lhs1b = rhs1.subs(d).subs rules
rhs1b = rhs1.subs(d).subs rules.subs(v[i] * wHat, 0)
display_latex_equation(lhs1b, rhs1b)

$$ \mathcal{E}_2 \hat{w} x_W^{(1)} + \mathcal{E}_3 \hat{w} \chi^{(0)}_i v_i v_j \xi^{(1)}_j + \dot{\mathcal{E}}_2 \hat{w} \lambda^{(1)} \chi^{(0)}_i v_i = 0 $$

(observe that, in the RHS, $\langle v_i, \hat{w} \rangle = 0$ since $V$ and $W$ are orthogonal subspaces). The solution to the above problem is expressed as a linear combination of the $w_{ij}$ and $w_i$ -defined by the variational problems (2.5) and (2.4), respectively–, delivering Eq. (7.1).

$$ \text{rules}[x1_W] = x0[k] * x1[l] * w[k, l] + \lambda1 * x0[k] * w[k] $$
## 7.4 The variational problem of order 2

The terms of second order are tested against elements of V only.

```R
lhs2 = expand(2 * lhs.coeff(η, 2).subs(u hat, v[i]).subs rules))
rhs2 = expand(2 * rhs.coeff(η, 2).subs(u hat, v[i]).subs rules))
```

<figure><img src="image_4.png" /></figure>

### Simplification of these expressions

Use some orthogonality conditions in the right-hand side.

```R
rhs2 = rhs2.subs({
  x0_V * v[i]: χ0[i],
  x1_V * v[i]: χ1[i],
  v[i] * w[k]: 0,
  v[i] * w[k, 1]: 0
})
```

Plug expansions

$$u^{(1)} = \xi_i^{(1)} v_i, \quad u_V^{(2)} = \xi_i^{(2)} v_i, \quad x_V^{(0)} = \chi_i^{(0)} v_i, \quad \text{and} \quad x_V^{(1)} = \chi_i^{(1)} v_i.$$

```R
d = dict()
d[E2.dot * v[i] * x0_V] = χ0[j] * E2.dot * v[i] * v[j]
d[E2.dot * v[i] * x1_V] = χ1[j] * E2.dot * v[i] * v[j]
d[E2.ddot * x0_V * v[i]] = χ0[j] * E2.ddot * v[i] * v[j]
d[E3 * v[i] * x1_V * u1] = χ1[j] * ξ1[k] * E3 * v[i] * v[j] * v[k]
d[E3 * v[i] * x0_V * u2_V] = χ0[j] * ξ2[k] * E3 * v[i] * v[j] * v[k]
d[E3.dot * v[i] * x0_V * u1] = χ0[j] * ξ1[k] * E3.dot * v[i] * v[j] * v[k]
d[E4 * v[i] * x0_V * u1 * u1] = (χ0[j] * ξ1[k] * ξ1[l] * E4 * v[i] * v[j] * v[k] * v[l])
```

```R
lhs2 = lhs2.subs(d)
```

Rename some indices.

```R
lhs2 = lhs2.subs({
  χ0[k] * E2.dot * v[i] * w[k]: χ0[j] * E2.dot * v[i] * w[j],
  χ0[k] * ξ1[l] * w[k, 1]: χ0[j] * ξ1[k] * w[j, k]
})
```

$$ \begin{aligned}
\mathcal{E}_3(x_V^{(0)}, u_W^{(2)}, v_i) &= \chi_j^{(0)} \xi_k^{(1)} [\xi_l^{(1)} \mathcal{E}_3(v_i, v_j, w_{kl}) + 2\lambda^{(1)} \mathcal{E}_3(v_i, v_j, w_k)] \\
&= \chi_j^{(0)} \xi_k^{(1)} \xi_l^{(1)} \mathcal{E}_3(v_i, v_j, w_{kl}) + 2\lambda^{(1)} \chi_j^{(0)} \xi_k^{(1)} \dot{\mathcal{E}}_2(v_k, w_{ij})
\end{aligned} $$
lhs2 =lhs2.subs(E3 * v[i] * x0_V * u2_W,
(χ0[j] * ξ1[k] * ξ1[l]
* E3 * v[i] * v[j] * w[k,l]
+ 2 * λ1 * χ0[j] * ξ1[k]
* E2 Dot * v[k] * w[i, j]).expand())

$$
\chi_k^{(0)} \mathcal{E}_3(v_i, u^{(1)}, w_k) = \chi_j^{(0)} \xi_k^{(1)} \mathcal{E}_3(v_i, v_k, w_j) = \chi_j^{(0)} \xi_k^{(1)} \dot{\mathcal{E}}_2(v_j, w_{ik})
$$

lhs2 =lhs2.subs(χ0[k] * E3 * v[i] * u1 * w[k],
χ0[j] * ξ1[k] * E2 dot * v[j] * w[i, k]).expand())

$$
\begin{align*}
\chi_j^{(0)} \xi_k^{(1)} \mathcal{E}_3(v_i, w_{j,k}, u^{(1)}) &= \chi_j^{(0)} \xi_k^{(1)} \xi_l^{(1)} \mathcal{E}_3(v_i, v_l, w_{j,k}) \\
&= \frac{1}{2} \chi_j^{(0)} \xi_k^{(1)} \xi_l^{(1)} [\mathcal{E}_3(v_i, v_k, w_{j,l}) + \mathcal{E}_3(v_i, v_l, w_{j,k})]
\end{align*}
$$

lhs2 =lhs2.subs(χ0[j] * ξ1[k] * E3 * v[i] * w[j, k] * u1,
χ0[j] * ξ1[k] * ξ1[l] / 2
* (E3 * v[i] * v[k] * w[j, l]
+ E3 * v[i] * v[l] * w[j, k]).expand())

lhs2 =lhs2.subs({
E2 dot * v[i] * v[j]: E dot[i, j],
E2 ddot * v[i] * v[j]: (E ddot[i, j] - 2 * E2 dot * v[i] * w[j]),
E3 * v[i] * v[j] * v[k]: E[i, j, k],
E3 dot * v[i] * v[j] * v[k]: (E dot[i, j, k]
- E2 dot * v[i] * w[j, k]
- E2 dot * v[j] * w[i, k]
- E2 dot * v[k] * w[i, j]),
E4 *v[i] * v[j] * v[k] * v[l]: (E[i, j, k, l]
- E3 * v[i] * (v[j] * w[k, l]
+ v[k] * w[j, l]
+ v[l] * w[j, k]))
}).expand())

And Eq. (7.3) is finally retrieved.

expected =((E[i, j, k, l] * ξ1[k] * ξ1[l]
+ λ1 * (2 * E Dot[i, j, k] * ξ1[k] + λ1 * E Ddot[i, j])
+ E[i, j, k] * ξ2[k]
+ λ2 * E Dot[i, j]) * χ0[j]
+ 2 * (E[i, j, k] * ξ1[k] + λ1 * E Dot[i, j]) * χ1[j])
assert expand(lhs2 - expected) == 0

expected = 2 * α1 * χ1[i] + α2 * χ0[i]

assert expandTre RHS2 - expected) == 0
# 8 Asymmetric bifurcation

In this chapter, we consider a bifurcated branch for which $\lambda^{(1)} \neq 0$. The bifurcation equation (5.4) shows that necessarily, $E_{ijk}$ is not identically nul. This equation has at most $(2^m - 1)$ pairs of real solutions $(\lambda^{(1)}, u^{(1)})$ et $(-\lambda^{(1)}, -u^{(1)})$; furthermore, multiplication by $\xi_i^{(1)}$ shows that

$$ \lambda^{(1)} = - \frac{E_{ijk} \xi_i^{(1)} \xi_j^{(1)} \xi_k^{(1)}}{2 \dot{E}_{ij} \xi_i^{(1)} \xi_j^{(1)}}. \qquad (8.1) $$

<figure><img src="image_5.png" /></figure>

## Warning

I can't prove that the bifurcation equation (5.4) has at most $(2^m - 1)$ pairs of real solutions.

Along the bifurcated branch, we have $\lambda = \lambda_0 + \eta \lambda^{(1)} + o(\eta)$, and $\eta$ can be eliminated. In other words, $\eta = \lambda (\lambda^{(1)} = 1$ and $\lambda^{(2)} = \lambda^{(3)} = \dots = 0)$ can be selected as a parameter. It is therefore possible to express the bifurcated branch as a function of $\lambda$: $u(\lambda)$. For example, combining Eqs. (5.7) and (6.2), we find that

$$ \begin{aligned} \mathcal{E}_{,uu}[u(\eta), \lambda(\eta); u^{(1)}, u^{(1)}] &= \eta [\mathcal{E}_3(u^{(1)}, u^{(1)}, u^{(1)}) + \lambda^{(1)} \dot{\mathcal{E}}_2(u^{(1)}, u^{(1)})] + o(\eta) \\ &= -\eta \lambda^{(1)} \dot{\mathcal{E}}_2(u^{(1)}, u^{(1)}) + o(\eta), \end{aligned} $$

or

$$ \mathcal{E}_{,uu}[u(\lambda), \lambda; u^{(1)}, u^{(1)}] = -(\lambda - \lambda_0) \dot{\mathcal{E}}_2(u^{(1)}, u^{(1)}) + o(\lambda - \lambda_0). \qquad (8.2) $$

For $\lambda < \lambda_0$, the above quantity is negative (since $\dot{\mathcal{E}}_2$ is negative definite). In other words: **for asymmetric bifurcations, below the critical load, the bifurcated branch is unstable.**

To investigate the stability above the critical load, we need to analyse the sign of the eigenvalues $\alpha$ of the Hessian. At first order, $\alpha = \eta \alpha^{(1)} + o(\eta)$, where $\alpha^{(1)}$ is an eigenvalue of $(E_{ijk} \xi_k^{(1)} + \lambda^{(1)} \dot{E}_{ij})$ (see Chapter 7). Let $\alpha_{min}$ and $\alpha_{max}$ be the minimum and maximum eigenvalues of this second-order tensor. Three cases must be discussed

1. If $\alpha_{min} \alpha_{max} > 0$, then $(E_{ijk} \xi_k^{(1)} + \lambda^{(1)} \dot{E}_{ij})$ is positive or negative definite: all eigenvalues have the same sign, $\epsilon \in \{-1, +1\}$. Then the sign of the eigenvalues $\alpha$ of the Hessian is $\epsilon \eta$ and there is a stability switch at the critical load. Since the bifurcated branch is unstable below the critical load, this means that it is stable above the critical load.
2. If $\alpha_{min} \alpha_{max} < 0$, then the extremal eigenvalues of the Hessian are $\eta \alpha_{min}$ and $\eta \alpha_{max}$, the product of which is $\eta^2 \alpha_{min} \alpha_{max} < 0$. The bifurcated branch is unstable for all values of $\lambda$.

3. If $\alpha_{min} \alpha_{max} = 0$, the analysis is inconclusive.

To close this section, it is observed that the dominant term of the expansion (6.1) of the potential energy along the bifurcated branch is of the third order in $\eta$

$$ \mathcal{E}[u(\eta), \lambda(\eta)] = \mathcal{E}[u^* \circ \lambda(\eta), \lambda(\eta)] + \frac{1}{6} \lambda^{(1)} \eta^3 \dot{E}_{ij} \xi_i^{(1)} \xi_j^{(1)} + o(\eta^3). $$

Eliminating $\lambda$ and plugging expression (8.1) of $\lambda^{(1)}$ delivers the expression of the potential energy, where $\lambda$ is the parameter

$$ \begin{aligned} \mathcal{E}[u(\lambda), \lambda] &= \mathcal{E}[u^*(\lambda), \lambda] + \frac{(\lambda - \lambda_0)^3}{6(\lambda^{(1)})^2} \dot{E}_{ij} \xi_i^{(1)} \xi_j^{(1)} + o(\lambda^3) \\ &= \mathcal{E}[u^*(\lambda), \lambda] + \frac{2(\dot{E}_{ij} \xi_i^{(1)} \xi_j^{(1)})^3}{3(E_{ijk} \xi_i^{(1)} \xi_j^{(1)} \xi_k^{(1)})^2} (\lambda - \lambda_0)^3 + o(\lambda^3). \end{aligned} $$

Recalling that $\dot{E}_{ij} \xi_i^{(1)} \xi_j^{(1)} < 0$, it is found that, above the critical load, the potential energy is smaller along the bifurcated branch than along the fundamental branch.

**i** Note

As expected, the above expression does not depend on the scaling of $u^{(1)}$ (of the $\xi_i^{(1)}$).

**i** Note

It has been shown in Chapter 3 that, when $E_{ijk}$ is not identically null, the critical point is unstable.
9 A particular case of symmetric bifurcation

In this chapter, we consider the case $E_{ijk} = 0$ for all $i, j, k = 1, \dots, m$. Then, from Eq. (5.4), $\lambda^{(1)} = 0$ on all bifurcated branches. It is assumed that, on the bifurcated branch under consideration, the next term of the expansion of $\lambda$ is non-zero: $\lambda^{(2)} \neq 0$. The bifurcation is symmetric, and the bifurcation equation (5.5) reduces to

$$ \frac{1}{3} E_{ijkl} \xi_j^{(1)} \xi_k^{(1)} \xi_l^{(1)} + \lambda^{(2)} \dot{E}_{ij} \xi_j^{(1)} = 0, \qquad (9.1) $$

which has at most $(3^m - 1)/2$ pairs of real solutions $(\lambda^{(2)}, u^{(1)})$ and $(-\lambda^{(2)}, -u^{(1)})$.

<figure><img src="image_6.png" /></figure>

Warning

I can't prove that the bifurcation equation (9.1) has at most $(3^m - 1)/2$ pairs of real solutions.

Upon multiplication by $\xi_i^{(1)}$, the above equation delivers the following expression of $\lambda^{(2)}$

$$ \lambda^{(2)} = - \frac{E_{ijkl} \xi_i^{(1)} \xi_j^{(1)} \xi_k^{(1)} \xi_l^{(1)}}{3 \dot{E}_{ij} \xi_i^{(1)} \xi_j^{(1)}}. \qquad (9.2) $$

Since $\dot{E}_{ij} \xi_i^{(1)} \xi_j^{(1)} < 0$, $\lambda^{(2)}$ has the same sign as $E_{ijkl} \xi_i^{(1)} \xi_j^{(1)} \xi_k^{(1)} \xi_l^{(1)}$. In other words, if $E_{ijkl} \xi_i^{(1)} \xi_j^{(1)} \xi_k^{(1)} \xi_l^{(1)} > 0$, (resp. $< 0$) then the bifurcated branch exists above (resp. below) the critical load $\lambda_0$ only.

Turning now to the eigenpairs of the Hessian of the energy along the bifurcated branch, Eq. (7.2) shows that $\alpha^{(1)} = 0$. Then $\alpha = \alpha^{(2)} \eta^2/2 + o(\eta^2)$ and, from Eq. (7.3)

$$ (E_{ijkl} \xi_k^{(1)} \xi_l^{(1)} + \lambda^{(2)} \dot{E}_{ij}) \chi_j^{(0)} = \alpha^{(2)} \chi_i^{(0)}. $$

If $(E_{ijkl} \xi_k^{(1)} \xi_l^{(1)} + \lambda^{(2)} \dot{E}_{ij})$ is positive definite, then the bifurcated branch is stable (note that, in that case, the bifurcated branch exists above the critical load only). If one of the eigenvalues of this tensor is $< 0$, then the bifurcated branch is unstable. The stability is undecided when all eigenvalues are $\ge 0$.

i Note

Note that, from Eq. (9.1),

$$ E_{ijkl} \xi_i^{(1)} \xi_j^{(1)} \xi_k^{(1)} \xi_l^{(1)} + \lambda^{(2)} \dot{E}_{ij} \xi_i^{(1)} \xi_j^{(1)} \xi_k^{(1)} \xi_l^{(1)} = \frac{2}{3} E_{ijkl} \xi_i^{(1)} \xi_j^{(1)} \xi_k^{(1)} \xi_l^{(1)}. $$
To conclude this section, it is observed that, when $\lambda^{(1)} = 0$, the dominant term of the potential energy along the bifurcated branch is of the fourth order, see Eq. (6.1). Combining with Eq. (9.2),

$$ \begin{aligned} \mathcal{E}[u(\eta), \lambda(\eta)] &= \mathcal{E}[u^* \circ \lambda(\eta), \lambda(\eta)] + \frac{1}{24}\eta^4 (E_{ijkl} \xi_i^{(1)} \xi_j^{(1)} \xi_k^{(1)} \xi_l^{(1)} \\ & \qquad + 6\lambda^{(2)} \dot{E}_{ij} \xi_i^{(1)} \xi_j^{(1)}) + o(\eta^4) \\ &= \mathcal{E}[u^* \circ \lambda(\eta), \lambda(\eta)] - \frac{1}{24}\eta^4 E_{ijkl} \xi_i^{(1)} \xi_j^{(1)} \xi_k^{(1)} \xi_l^{(1)} + o(\eta^4). \end{aligned} \quad (9.3) $$

The expansion $\lambda = \lambda_0 + \lambda^{(2)} \eta^2/2 + o(\eta^2)$ can be inverted as follows

$$ \eta^4 = \frac{4(\lambda - \lambda_0)^2}{(\lambda^{(2)})^2} + o(\lambda^2) = \frac{36(\dot{E}_{ij} \xi_i^{(1)} \xi_j^{(1)})^2}{(E_{ijkl} \xi_i^{(1)} \xi_j^{(1)} \xi_k^{(1)} \xi_l^{(1)})^2} (\lambda - \lambda_0)^2 $$

and expression (9.3) reads

$$ \mathcal{E}[u(\eta), \lambda(\eta)] = \mathcal{E}[u^* \circ \lambda(\eta), \lambda(\eta)] - \frac{3(\dot{E}_{ij} \xi_i^{(1)} \xi_j^{(1)})^2}{2E_{ijkl} \xi_i^{(1)} \xi_j^{(1)} \xi_k^{(1)} \xi_l^{(1)}} (\lambda - \lambda_0)^2 + o(\lambda^2). $$

Again, the above expression does not depend on the scaling of $u^{(1)}$ (of the $\xi_i^{(1)}$). Note that, if $E_{ijkl} \xi_i^{(1)} \xi_j^{(1)} \xi_k^{(1)} \xi_l^{(1)} > 0$, then only loads that are greater than the critical load can be reached on the bifurcated branch, where the energy is lower than the fundamental branch.
# 10 The case of a single mode

In this chapter, we discuss the case $m = 1$; all tensors considered previously ($\dot{E}_{ij}$, $E_{ijk}$, $E_{ijkl}$, etc) then reduce to simple scalars. To avoid ambiguity, indices are kept: $\dot{E}_{11}$, $E_{111}$, $E_{1111}$. Since $\dot{E}_2$ is negative definite over $V$ (see note in 2.2), we have $\dot{E}_{11} < 0$.

It is first observed that the following conditions are necessary to ensure stability of the critical point

$$E_{111} = 0 \quad \text{and} \quad E_{1111} \ge 0,$$

which shows that asymmetric critical points are always unstable.

## 10.1 Asymmetric bifurcation

We first consider the case $E_{111} \neq 0$. Owing to the discussion above, the bifurcation point is unstable. It follows from Eq. (8.1) that

$$\lambda^{(1)} = -\frac{E_{111} \xi_1^{(1)}}{\dot{E}_{11}}.$$

Then, Eq. (8.2) delivers, with $u^{(1)} = \xi_1^{(1)} v_1$

$$(\xi_1^{(1)})^2 \mathcal{E}_{,uu}[u(\lambda), \lambda; v_1, v_1] = -(\lambda - \lambda_0)(\xi_1^{(1)})^2 \dot{\mathcal{E}}_2(v_1, v_1) + o(\lambda - \lambda_0)$$

or

$$\mathcal{E}_{,uu}[u(\lambda), \lambda; v_1, v_1] = -(\lambda - \lambda_0) \dot{E}_{11} + o(\lambda - \lambda_0).$$

Since $\dot{E}_{11} < 0$ we conclude that asymmetric bifurcations branches are unstable for $\lambda \le \lambda_0$ and stable for $\lambda > \lambda_0$ (stability switch).

## 10.2 Symmetric bifurcation

We now consider the case $E_{111} = 0$. From the general discussion of Chapter 3, the critical point is stable if $E_{1111} > 0$ and unstable if $E_{1111} < 0$.
<figure><img src="image_7.png" /></figure>

Warning

Is this really true? Even for $m = 1$ the proof that conditions (3.4) are sufficient to ensure stability does not seem so obvious to me.

The bifurcation equation (9.1) reduces to

$$ \frac{1}{3} E_{1111} (\xi_1^{(1)})^2 + \lambda^{(2)} \dot{E}_{11} = 0, $$

which in particular shows that $\lambda^{(2)}$ has the same sign as $E_{1111}$. Since the expansion of $\lambda$ reads: $\lambda = \lambda_0 + \lambda^{(2)}\eta^2/2 + o(\eta^2)$, the bifurcation branch exists only for loads above the critical load ($\lambda \ge \lambda_0$) if $E_{1111} > 0$ and only for loads below the critical load ($\lambda \le \lambda_0$) if $E_{1111} < 0$.

From Eq. (6.2), the hessian of the energy along the bifurcated branch reads

$$ \begin{aligned} \mathcal{E}_{uu}[u(\eta), \lambda(\eta); v_1, v_1] &= \frac{1}{2}\eta^2 \left[ (\xi_1^{(1)})^2 E_{1111} + \lambda^{(2)} \dot{E}_{11} \right] + o(\eta^2) \\ &= -\eta^2 \lambda^{(2)} \dot{E}_{11} + o(\eta^2), \end{aligned} $$

which has the sign of $\lambda^{(2)}$. Therefore the Hessian is positive (resp. negative) definite if $E_{1111} > 0$ (resp. $< 0$).

To sum up, if $E_{1111} > 0$, then the bifurcation branch (including the critical point) is stable and exists only for loads greater than the critical load. Conversely, if $E_{1111} < 0$, then the bifurcation branch (including the critical point) is unstable and exists only for loads lower than the critical load.
# 11 The lsk.symbols module

This module defines all symbols required for the SymPy derivations (see below).

```python
import lsk.symbols
```

```python
%psource lsk.symbols
```

```python
from sympy import Idx, IndexedBase, Symbol
```

```python
E2 = Symbol(r"\E_2")
```

```python
E3 = Symbol(r"\E_3")
```

```python
E4 = Symbol(r"\E_4")
```

```python
E2Dot = Symbol(r"\dot{\E}_2")
```

```python
E3Dot = Symbol(r"\dot{\E}_3")
```

```python
E2DDot = Symbol(r"\ddot{\E}_2")
```

```python
E_λ = Symbol(r"\E_{\lambda}")
```

```python
E_λλ = Symbol(r"\E_{\lambda\lambda}")
```

```python
E_λλλ = Symbol(r"\E_{\lambda\lambda\lambda}")
```

```python
E_λλλλ = Symbol(r"\E_{\lambda\lambda\lambda\lambda}")
```

```python
E_uλ = Symbol(r"\E_{u\lambda}")
```

```python
E_uλλ = Symbol(r"\E_{u\lambda\lambda}")
```

```python
E_uλλλ = Symbol(r"\E_{u\lambda\lambda\lambda}")
```

```python
E_uuλ = Symbol(r"\E_{uu\lambda}")
```

```python
E_uuλλ = Symbol(r"\E_{uu\lambda\lambda}")
```

```python
E_uuuλ = Symbol(r"\E_{uuu\lambda}")
```

```python
u0Dot = Symbol(r"\dot{u}_0")
```

```python
u0DDot = Symbol(r"\ddot{u}_0")
```

```python
u0DDdot = Symbol(r"\dddot{u}_0")
```

```python
u0DDddot = Symbol(r"\ddddot{u}_0")
```

```python
i = Idx("i")
```

```python
j = Idx("j")
```

```python
k = Idx("k")
```

```python
l = Idx("l")
```

```python
λ = IndexedBase(r"\lambda")
```

```python
ξ = IndexedBase(r"\xi")
```
ξ₁ = IndexedBase(r"{\order[1]{\xi}}")
ξ₂ = IndexedBase(r"{\order[2]{\xi}}")
v = IndexedBase("v")
w = IndexedBase("w")
E = IndexedBase("E")
E Dot = IndexedBase("\dot{E}")
E Ddot = IndexedBase("\ddot{E}")

u₁ = Symbol(r"{\order[1]{u}}")
u₂ = Symbol(r"{\order[2]{u}}")
u₂_V = Symbol(r"{\order[2]{u}_V}")
u₂_W = Symbol(r"{\order[2]{u}_W}")
u₃ = Symbol(r"{\order[3]{u}}")
u₄ = Symbol(r"{\order[4]{u}}")

λ₁ = Symbol(r"{\order[1]{\lambda}}")
λ₂ = Symbol(r"{\order[2]{\lambda}}")
λ₃ = Symbol(r"{\order[3]{\lambda}}")
λ₄ = Symbol(r"{\order[4]{\lambda}}")

η = Symbol(r"\eta")

u_hat = Symbol(r"\hat{u}")

v_hat = Symbol(r"\hat{v}")

w_hat = Symbol(r"\hat{w}")

Note that the symbols: λ, ξ, v, w, E, E and E_Ddot are all instances of the class IndexedBase.
In other words, they can be indexed with symbolic indices

$$
\bullet v[i] \rightarrow v_i
$$

$$
\bullet E[i, j, k] \rightarrow E_{ijk}
$$

etc.
References

Chakrabarti, Aditi, Serge Mora, Franck Richard, Ty Phou, Jean-Marc Fromental, Yves Pomeau, and Basile Audoly. 2018. "Selection of Hexagonal Buckling Patterns by the Elastic Rayleigh-Taylor Instability." *Journal of the Mechanics and Physics of Solids* 121 (December): 234–57. https://doi.org/10.1016/j.jmps.2018.07.024.

Koiter, W. T., and A. M. A. van der Heijden. 2009. *W.T. Koiter's Elastic Stability of Solids and Structures*. Cambridge ; New York: Cambridge University Press.

Koiter, Warner Tjardus. 1945. "The Stability of Elastic Equilibrium." PhD thesis, Delft: Technische Hooge School Delft.

Nguyen, Quoc-Son. 2000. *Stabilité et mécanique non linéaire*. Paris: Hermès Science Publications.

Potier-Ferry, Michel. 1987. "Foundations of Elastic Postbuckling Theory." In *Buckling and Post-Buckling*, 1–82. Lecture Notes in Physics. Springer, Berlin, Heidelberg. https://doi.org/10.1007/BFb0009197.

Triantafyllidis, Nicolas. 2017. "Stability of Solids: From Structures to Materials." Lecture Notes. Palaiseau, France: École polytechnique.