from sympy import *
from lsk.display import *
from lsk.symbols import *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 Section 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 \E(u, \lambda)\) and its differentials at the critical point \((u_0, \lambda_0)\), as well as the fundamental path \(\lambda \mapsto u^\ast(\lambda)\) and its derivatives at \(\lambda = \lambda_0\). It will therefore be convenient to express \(\E\) and \(u^\star\) as Taylor expansions with respect to \(u\) and \(\lambda\).
We start with the Taylor expansion of the energy \(\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\).
Code
d = {
r"\E_{,\lambda}(u_0, \lambda_0)": E_λ,
r"\E_{,uu}(u_0, \lambda_0)": E2,
r"\E_{,u\lambda}(u_0, \lambda_0)": E_uλ,
r"\E_{,\lambda\lambda}(u_0, \lambda_0)": E_λλ,
r"\E_{,uuu}(u_0, \lambda_0)": E3,
r"\E_{,uu\lambda}(u_0, \lambda_0)": E_uuλ,
r"\E_{,u\lambda\lambda}(u_0, \lambda_0)": E_uλλ,
r"\E_{,\lambda\lambda\lambda}(u_0, \lambda_0)": E_λλλ,
r"\E_{,uuuu}(u_0, \lambda_0)": E4,
r"\E_{,uuu\lambda}(u_0, \lambda_0)": E_uuuλ,
r"\E_{,uu\lambda\lambda}(u_0, \lambda_0)": E_uuλλ,
r"\E_{,u\lambda\lambda\lambda}(u_0, \lambda_0)": E_uλλλ,
r"\E_{,\lambda\lambda\lambda\lambda}(u_0, \lambda_0)": E_λλλλ,
}
display_latex_dict(d, num_cols=3)\(\displaystyle \begin{aligned}\E_{,\lambda}(u_0, \lambda_0)&=\E_{\lambda}&\E_{,uu}(u_0, \lambda_0)&=\E_{2}&\E_{,u\lambda}(u_0, \lambda_0)&=\E_{u\lambda}\\\E_{,\lambda\lambda}(u_0, \lambda_0)&=\E_{\lambda\lambda}&\E_{,uuu}(u_0, \lambda_0)&=\E_{3}&\E_{,uu\lambda}(u_0, \lambda_0)&=\E_{uu\lambda}\\\E_{,u\lambda\lambda}(u_0, \lambda_0)&=\E_{u\lambda\lambda}&\E_{,\lambda\lambda\lambda}(u_0, \lambda_0)&=\E_{\lambda\lambda\lambda}&\E_{,uuuu}(u_0, \lambda_0)&=\E_{4}\\\E_{,uuu\lambda}(u_0, \lambda_0)&=\E_{uuu\lambda}&\E_{,uu\lambda\lambda}(u_0, \lambda_0)&=\E_{uu\lambda\lambda}&\E_{,u\lambda\lambda\lambda}(u_0, \lambda_0)&=\E_{u\lambda\lambda\lambda}\\\E_{,\lambda\lambda\lambda\lambda}(u_0, \lambda_0)&=\E_{\lambda\lambda\lambda\lambda}&\end{aligned}\)
where
- \(\E_\lambda\), \(\E_{\lambda\lambda}\), \(\E_{\lambda\lambda\lambda}\) and \(\E_{\lambda\lambda\lambda\lambda}\) are scalar quantities,
- \(\E_{u\lambda}\), \(\E_{u\lambda\lambda}\) and \(\E_{u\lambda\lambda\lambda}\) are linear forms,
- \(\E_2\), \(\E_{uu\lambda}\) and \(\E_{uu\lambda\lambda}\) are bilinear forms,
- \(\E_3\) and \(\E_{uuu\lambda}\) are trilinear forms,
- \(\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 * (α * u + β * v) * w (resp. E2 * u * v - E2 * v * u) should be understood as E2(α * u + β * v, w) (resp. E2(u, v) - E2(v, u)). In both cases, the expressions would are 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: λ * E2 * u * v, the first * is a true multiplication, while the other * refer to function application.
The energy \(\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.
u, λ = symbols(r"u \lambda")
E = (λ * 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_uuuλ * u**3 + 6 * λ**2 * E_uuλλ * u**2
+ 4 * λ**3 * E_uλλλ * u + λ**4 * E_λλλλ) / 24).expand()Code
display_latex_long_equation(r"\E(u, \lambda)", E, terms_per_line=7)\(\displaystyle \begin{aligned}\E(u, \lambda)={} &\frac{\E_{2} u^{2}}{2} + \frac{\E_{3} u^{3}}{6} + \frac{\E_{4} u^{4}}{24} + \frac{\E_{\lambda\lambda\lambda\lambda} \lambda^{4}}{24} + \frac{\E_{\lambda\lambda\lambda} \lambda^{3}}{6} + \frac{\E_{\lambda\lambda} \lambda^{2}}{2} + \E_{\lambda} \lambda\\ &+\frac{\E_{u\lambda\lambda\lambda} \lambda^{3} u}{6} + \frac{\E_{u\lambda\lambda} \lambda^{2} u}{2} + \E_{u\lambda} \lambda u + \frac{\E_{uu\lambda\lambda} \lambda^{2} u^{2}}{4} + \frac{\E_{uu\lambda} \lambda u^{2}}{2} + \frac{\E_{uuu\lambda} \lambda u^{3}}{6}\end{aligned}\)
The fundamental path \(\lambda \mapsto u^\star(\lambda)\) is also defined through its Taylor expansion.
u_star = (λ * u0_dot
+ λ**2 * u0_ddot / 2
+ λ**3 * u0_dddot / 6
+ λ**4 * u0_ddddot / 24)Code
display_latex_equation(r"u^\star(\lambda)", u_star)\(\displaystyle u^\star(\lambda)=\frac{\ddddot{u}_0 \lambda^{4}}{24} + \frac{\dddot{u}_0 \lambda^{3}}{6} + \frac{\ddot{u}_0 \lambda^{2}}{2} + \dot{u}_0 \lambda\)
Where \(\dot{u}_0\), \(\ddot{u}_0\), etc denote the derivatives of \(u^\star\) with respect to \(\lambda\), at \(\lambda = \lambda_0\).
Code
d = {f"\\frac{{\\D^{k}u^\\star}}{{\\D \\lambda^{k}}}"
"\\biggr \\rvert_{{\\lambda=\\lambda_0}}" : x
for k, x in enumerate([u0_dot, u0_ddot, u0_dddot, u0_ddddot], start=1)}
display_latex_dict(d, num_cols=4)\(\displaystyle \begin{aligned}\frac{\D^1u^\star}{\D \lambda^1}\biggr \rvert_{{\lambda=\lambda_0}}&=\dot{u}_0&\frac{\D^2u^\star}{\D \lambda^2}\biggr \rvert_{{\lambda=\lambda_0}}&=\ddot{u}_0&\frac{\D^3u^\star}{\D \lambda^3}\biggr \rvert_{{\lambda=\lambda_0}}&=\dddot{u}_0&\frac{\D^4u^\star}{\D \lambda^4}\biggr \rvert_{{\lambda=\lambda_0}}&=\ddddot{u}_0\\\end{aligned}\)
4.2 Elimination of the derivatives of the jacobian w.r.t. \(\lambda\)
Since the fundamental path \(\lambda \mapsto u^\star(\lambda)\) is an equilibrium path, the various differentials of the energy at the critical point are not linearly independent. To express the relationships between these forms, we define \(\mathcal R^\star(\lambda; \bullet)\) as the jacobian of the energy along the fundamental path \(u^\star(\lambda)\) \[ \mathcal R^\star(\lambda; \bullet) = \E_{,u}[u^\star(\lambda), \lambda; \bullet]. \]
Combining the expansions of \(\lambda \mapsto u^\star(\lambda)\) and \((u, \lambda) \mapsto \E(u, \lambda)\) delivers and expansion of \(\mathcal{R}^\star\) with respect to the powers of \(\lambda\), up to the fourth order.
R_star = (E.diff(u) * u_hat).subs(u, u_star).series(λ, 0, 5).removeO().expand()Code
display_latex_long_equation(r"\mathcal{R}^\ast(\lambda;"
+ sympy.latex(u_hat) + ")",
R_star, terms_per_line=4)\(\displaystyle \begin{aligned}\mathcal{R}^\ast(\lambda;\hat{u})={} &\E_{2} \dot{u}_0 \hat{u} \lambda + \frac{\E_{u\lambda\lambda\lambda} \hat{u} \lambda^{3}}{6} + \frac{\E_{u\lambda\lambda} \hat{u} \lambda^{2}}{2} + \E_{u\lambda} \hat{u} \lambda\\ &+\frac{\E_{2} \ddot{u}_0 \hat{u} \lambda^{2}}{2} + \frac{\E_{3} \dot{u}_0^{2} \hat{u} \lambda^{2}}{2} + \frac{\E_{uu\lambda\lambda} \dot{u}_0 \hat{u} \lambda^{3}}{2} + \E_{uu\lambda} \dot{u}_0 \hat{u} \lambda^{2}\\ &+\frac{\E_{2} \dddot{u}_0 \hat{u} \lambda^{3}}{6} + \frac{\E_{uu\lambda\lambda} \ddot{u}_0 \hat{u} \lambda^{4}}{4} + \frac{\E_{uu\lambda} \ddot{u}_0 \hat{u} \lambda^{3}}{2} + \frac{\E_{uuu\lambda} \dot{u}_0^{2} \hat{u} \lambda^{3}}{2}\\ &+\frac{\E_{2} \ddddot{u}_0 \hat{u} \lambda^{4}}{24} + \frac{\E_{3} \ddot{u}_0^{2} \hat{u} \lambda^{4}}{8} + \frac{\E_{4} \dot{u}_0^{3} \hat{u} \lambda^{3}}{6} + \frac{\E_{uu\lambda} \dddot{u}_0 \hat{u} \lambda^{4}}{6}\\ &+\frac{\E_{3} \dddot{u}_0 \dot{u}_0 \hat{u} \lambda^{4}}{6} + \frac{\E_{3} \ddot{u}_0 \dot{u}_0 \hat{u} \lambda^{3}}{2} + \frac{\E_{4} \ddot{u}_0 \dot{u}_0^{2} \hat{u} \lambda^{4}}{4} + \frac{\E_{uuu\lambda} \ddot{u}_0 \dot{u}_0 \hat{u} \lambda^{4}}{2}\end{aligned}\)
Of course, since \(\lambda \mapsto u^\ast(\lambda)\) is an equilibrium path, we have \(\mathcal R^\ast(\lambda; \bullet) = 0\) for all \(\lambda\). Therefore, all coefficients of the above polynomial in \(\lambda\) are null, which delivers expressions of \(\E_{u\lambda}\), \(\E_{u\lambda\lambda}\) and \(\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.
Code
assert R_star.coeff(λ, 0) == 04.2.2 The term of order 1
This term delivers the following equation
eq = Eq(R_star.coeff(λ, 1), 0)Code
display(eq)\(\displaystyle \E_{2} \dot{u}_0 \hat{u} + \E_{u\lambda} \hat{u} = 0\)
for all \(\hat{u} \in U\). This equation delivers the following expression of \(\E_{,u\lambda}(u_0, \lambda_0)\)
sol = solve(eq, E_uλ)
mixed1[E_uλ] = sol[0]Code
display_latex_equation(E_uλ, mixed1[E_uλ])\(\displaystyle \E_{u\lambda}=- \E_{2} \dot{u}_0\)
4.2.3 The term of order 2
eq = Eq(R_star.coeff(λ, 2).subs(mixed1), 0)Code
display(eq)\(\displaystyle \frac{\E_{2} \ddot{u}_0 \hat{u}}{2} + \frac{\E_{3} \dot{u}_0^{2} \hat{u}}{2} + \frac{\E_{u\lambda\lambda} \hat{u}}{2} + \E_{uu\lambda} \dot{u}_0 \hat{u} = 0\)
for all \(\hat{u} \in U\). This equation delivers the following expression of \(\E_{,u\lambda\lambda}(u_0, \lambda_0)\)
sol = solve(eq, E_uλλ)
mixed1[E_uλλ] = sol[0]Code
display_latex_equation(E_uλλ, mixed1[E_uλλ])\(\displaystyle \E_{u\lambda\lambda}=- \E_{2} \ddot{u}_0 - \E_{3} \dot{u}_0^{2} - 2 \E_{uu\lambda} \dot{u}_0\)
4.2.4 The term of order 3
eq = Eq(R_star.coeff(λ, 3).subs(mixed1), 0).expand()Code
display(eq)\(\displaystyle \frac{\E_{2} \dddot{u}_0 \hat{u}}{6} + \frac{\E_{3} \ddot{u}_0 \dot{u}_0 \hat{u}}{2} + \frac{\E_{4} \dot{u}_0^{3} \hat{u}}{6} + \frac{\E_{u\lambda\lambda\lambda} \hat{u}}{6} + \frac{\E_{uu\lambda\lambda} \dot{u}_0 \hat{u}}{2} + \frac{\E_{uu\lambda} \ddot{u}_0 \hat{u}}{2} + \frac{\E_{uuu\lambda} \dot{u}_0^{2} \hat{u}}{2} = 0\)
for all \(\hat{u} \in U\). This equation delivers the following expression of \(\E_{,u\lambda\lambda\lambda}(u_0, \lambda_0)\)
sol = solve(eq, E_uλλλ)
mixed1[E_uλλλ] = sol[0]Code
display_latex_equation(E_uλλ, mixed1[E_uλλ])\(\displaystyle \E_{u\lambda\lambda}=- \E_{2} \ddot{u}_0 - \E_{3} \dot{u}_0^{2} - 2 \E_{uu\lambda} \dot{u}_0\)
4.3 Elimination of the remaining mixed derivatives
So far, we have found the following expressions
Code
display_latex_dict(mixed1, num_cols=1)\(\displaystyle \begin{aligned}\E_{u\lambda}&=- \E_{2} \dot{u}_0\\\E_{u\lambda\lambda}&=- \E_{2} \ddot{u}_0 - \E_{3} \dot{u}_0^{2} - 2 \E_{uu\lambda} \dot{u}_0\\\E_{u\lambda\lambda\lambda}&=- \E_{2} \dddot{u}_0 - 3 \E_{3} \ddot{u}_0 \dot{u}_0 - \E_{4} \dot{u}_0^{3} - 3 \E_{uu\lambda\lambda} \dot{u}_0 - 3 \E_{uu\lambda} \ddot{u}_0 - 3 \E_{uuu\lambda} \dot{u}_0^{2}\\\end{aligned}\)
We want to get rid of the remaining mixed derivatives, namely: \(\E_{uu\lambda}\), \(\E_{uuu\lambda}\) and \(\E_{uu\lambda\lambda}\). To do so, we introduce the derivatives \(\dot{\E}_2\), \(\ddot{\E}_2\) and \(\dot{\E}_3\) defined in Chapter 2.
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 \(\E_{uu\lambda}\) can first be expressed as a function of \(\dot{\E}_2\).
x = E_uuλ
lhs = E2_dot
rhs = E_uu_star.coeff(λ, 1)Code
display_latex_equation(lhs, rhs)\(\displaystyle \dot{\E}_2=\E_{3} \dot{u}_0 + \E_{uu\lambda}\)
sol = solve(Eq(lhs, rhs), x)
mixed2[x] = sol[0]Code
display_latex_equation(x, mixed2[x])\(\displaystyle \E_{uu\lambda}=- \E_{3} \dot{u}_0 + \dot{\E}_2\)
Then, the expression of \(\dot{E}_3\) delivers an expression of the mixed derivative \(\E_{uuu\lambda}\).
x = E_uuuλ
lhs = E3_dot
rhs = E_uuu_star.coeff(λ, 1)Code
display_latex_equation(lhs, rhs)\(\displaystyle \dot{\E}_3=\E_{4} \dot{u}_0 + \E_{uuu\lambda}\)
sol = solve(Eq(lhs, rhs), x)
mixed2[x] = sol[0]Code
display_latex_equation(x, mixed2[x])\(\displaystyle \E_{uuu\lambda}=- \E_{4} \dot{u}_0 + \dot{\E}_3\)
Finally, \(\ddot{\E}_2\) delivers an expression of the mixed derivative \(\E_{uu\lambda\lambda}\).
x = E_uuλλ
lhs = E2_ddot
rhs = 2 * E_uu_star.coeff(λ, 2).subs(mixed2).expand()Code
display_latex_equation(lhs, rhs)\(\displaystyle \ddot{\E}_2=\E_{3} \ddot{u}_0 - \E_{4} \dot{u}_0^{2} + \E_{uu\lambda\lambda} + 2 \dot{\E}_3 \dot{u}_0\)
sol = solve(Eq(lhs, rhs), x)
mixed2[x] = sol[0]Code
display_latex_equation(x, mixed2[x])\(\displaystyle \E_{uu\lambda\lambda}=- \E_{3} \ddot{u}_0 + \E_{4} \dot{u}_0^{2} + \ddot{\E}_2 - 2 \dot{\E}_3 \dot{u}_0\)
4.4 Summary: final expression of the energy
The following expressions were derived in Section 4.2
Code
display_latex_dict(mixed1, num_cols=1)\(\displaystyle \begin{aligned}\E_{u\lambda}&=- \E_{2} \dot{u}_0\\\E_{u\lambda\lambda}&=- \E_{2} \ddot{u}_0 - \E_{3} \dot{u}_0^{2} - 2 \E_{uu\lambda} \dot{u}_0\\\E_{u\lambda\lambda\lambda}&=- \E_{2} \dddot{u}_0 - 3 \E_{3} \ddot{u}_0 \dot{u}_0 - \E_{4} \dot{u}_0^{3} - 3 \E_{uu\lambda\lambda} \dot{u}_0 - 3 \E_{uu\lambda} \ddot{u}_0 - 3 \E_{uuu\lambda} \dot{u}_0^{2}\\\end{aligned}\)
and in Section 4.3
Code
display_latex_dict(mixed2, num_cols=1)\(\displaystyle \begin{aligned}\E_{uu\lambda}&=- \E_{3} \dot{u}_0 + \dot{\E}_2\\\E_{uuu\lambda}&=- \E_{4} \dot{u}_0 + \dot{\E}_3\\\E_{uu\lambda\lambda}&=- \E_{3} \ddot{u}_0 + \E_{4} \dot{u}_0^{2} + \ddot{\E}_2 - 2 \dot{\E}_3 \dot{u}_0\\\end{aligned}\)
Combining the above results allows to fully eliminate the mixed derivatives
Code
mixed = {k: v.subs(mixed2).expand() for k, v in mixed1.items()}
mixed.update(mixed2)
display_latex_dict(mixed, num_cols=1)\(\displaystyle \begin{aligned}\E_{u\lambda}&=- \E_{2} \dot{u}_0\\\E_{u\lambda\lambda}&=- \E_{2} \ddot{u}_0 + \E_{3} \dot{u}_0^{2} - 2 \dot{\E}_2 \dot{u}_0\\\E_{u\lambda\lambda\lambda}&=- \E_{2} \dddot{u}_0 + 3 \E_{3} \ddot{u}_0 \dot{u}_0 - \E_{4} \dot{u}_0^{3} - 3 \ddot{\E}_2 \dot{u}_0 - 3 \ddot{u}_0 \dot{\E}_2 + 3 \dot{\E}_3 \dot{u}_0^{2}\\\E_{uu\lambda}&=- \E_{3} \dot{u}_0 + \dot{\E}_2\\\E_{uuu\lambda}&=- \E_{4} \dot{u}_0 + \dot{\E}_3\\\E_{uu\lambda\lambda}&=- \E_{3} \ddot{u}_0 + \E_{4} \dot{u}_0^{2} + \ddot{\E}_2 - 2 \dot{\E}_3 \dot{u}_0\\\end{aligned}\)
These expressions can be plugged into the expansion of the energy.
E = E.subs(mixed).expand()Code
display_latex_long_equation(r"\E(u, \lambda)", E, terms_per_line=5)\(\displaystyle \begin{aligned}\E(u, \lambda)={} &\frac{\E_{2} u^{2}}{2} + \frac{\E_{3} u^{3}}{6} + \frac{\E_{\lambda\lambda\lambda} \lambda^{3}}{6} + \frac{\E_{\lambda\lambda} \lambda^{2}}{2} + \E_{\lambda} \lambda\\ &+\frac{\E_{4} u^{4}}{24} + \frac{\E_{\lambda\lambda\lambda\lambda} \lambda^{4}}{24} + \frac{\ddot{\E}_2 \lambda^{2} u^{2}}{4} + \frac{\dot{\E}_2 \lambda u^{2}}{2} + \frac{\dot{\E}_3 \lambda u^{3}}{6}\\ &+- \frac{\E_{2} \ddot{u}_0 \lambda^{2} u}{2} - \E_{2} \dot{u}_0 \lambda u + \frac{\E_{3} \dot{u}_0^{2} \lambda^{2} u}{2} - \dot{\E}_2 \dot{u}_0 \lambda^{2} u + \frac{\dot{\E}_3 \dot{u}_0^{2} \lambda^{3} u}{2}\\ &+- \frac{\E_{3} \ddot{u}_0 \lambda^{2} u^{2}}{4} - \frac{\E_{3} \dot{u}_0 \lambda u^{2}}{2} - \frac{\ddot{\E}_2 \dot{u}_0 \lambda^{3} u}{2} - \frac{\ddot{u}_0 \dot{\E}_2 \lambda^{3} u}{2} - \frac{\dot{\E}_3 \dot{u}_0 \lambda^{2} u^{2}}{2}\\ &+- \frac{\E_{2} \dddot{u}_0 \lambda^{3} u}{6} + \frac{\E_{3} \ddot{u}_0 \dot{u}_0 \lambda^{3} u}{2} - \frac{\E_{4} \dot{u}_0^{3} \lambda^{3} u}{6} + \frac{\E_{4} \dot{u}_0^{2} \lambda^{2} u^{2}}{4} - \frac{\E_{4} \dot{u}_0 \lambda u^{3}}{6}\end{aligned}\)
From which we deduce the expression of the residual \(\E_{,u}\)
Code
E_u = E.diff(u)
display_latex_long_equation(r"\E_{,u}(u, \lambda)", E_u, terms_per_line=5)\(\displaystyle \begin{aligned}\E_{,u}(u, \lambda)={} &\E_{2} u + \frac{\E_{3} \dot{u}_0^{2} \lambda^{2}}{2} + \frac{\E_{3} u^{2}}{2} + \frac{\E_{4} u^{3}}{6} + \dot{\E}_2 \lambda u\\ &+- \E_{2} \dot{u}_0 \lambda + \frac{\ddot{\E}_2 \lambda^{2} u}{2} - \dot{\E}_2 \dot{u}_0 \lambda^{2} + \frac{\dot{\E}_3 \dot{u}_0^{2} \lambda^{3}}{2} + \frac{\dot{\E}_3 \lambda u^{2}}{2}\\ &+- \frac{\E_{2} \dddot{u}_0 \lambda^{3}}{6} - \frac{\E_{2} \ddot{u}_0 \lambda^{2}}{2} - \frac{\E_{4} \dot{u}_0^{3} \lambda^{3}}{6} - \frac{\ddot{\E}_2 \dot{u}_0 \lambda^{3}}{2} - \frac{\ddot{u}_0 \dot{\E}_2 \lambda^{3}}{2}\\ &+\frac{\E_{3} \ddot{u}_0 \dot{u}_0 \lambda^{3}}{2} - \frac{\E_{3} \ddot{u}_0 \lambda^{2} u}{2} - \E_{3} \dot{u}_0 \lambda u + \frac{\E_{4} \dot{u}_0^{2} \lambda^{2} u}{2} - \dot{\E}_3 \dot{u}_0 \lambda^{2} u\\ &+- \frac{\E_{4} \dot{u}_0 \lambda u^{2}}{2}\end{aligned}\)
Code
E_uu = E.diff(u, 2)
display_latex_long_equation(r"\E_{,uu}(u, \lambda)", E_uu, terms_per_line=7)\(\displaystyle \begin{aligned}\E_{,uu}(u, \lambda)={} &\E_{2} + \E_{3} u + \frac{\E_{4} \dot{u}_0^{2} \lambda^{2}}{2} + \frac{\E_{4} u^{2}}{2} + \frac{\ddot{\E}_2 \lambda^{2}}{2} + \dot{\E}_2 \lambda + \dot{\E}_3 \lambda u\\ &+- \frac{\E_{3} \ddot{u}_0 \lambda^{2}}{2} - \E_{3} \dot{u}_0 \lambda - \E_{4} \dot{u}_0 \lambda u - \dot{\E}_3 \dot{u}_0 \lambda^{2}\end{aligned}\)
4.5 Implementation in the lsk.energy module
This module exposes three functions
create_E(u, λ): asymptotic expansion of the energy \(\E(u, \lambda)\),create_E_u(u, λ): asymptotic expansion of the jacobian \(\E_{,u}(u, \lambda)\),create_E_u(u, λ): asymptotic expansion of the hessian \(\E_{,uu}(u, \lambda)\),create_u_star(λ): asymptotic expansion of the fundamental branch \(u^\star(\lambda)\),
where u and λ must be SymPy expressions.
import lsk.energy
%psource lsk.energyfrom lsk.symbols import * __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_uuuλ : E3_dot - E4 * u0_dot, } 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_uuuλ * 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, λ)