# Presentation¶

I want to reproduce the symbolic (algebraic) computations done in §5.A of L.S.'s PhD thesis.

I want to only use a free and open-source software, so I'm using Python 3 with the Sympy module.

In :
%load_ext watermark
%watermark -v -m -p sympy -g

CPython 3.6.5
IPython 6.4.0

sympy 1.1.1

compiler   : GCC 7.3.0
system     : Linux
release    : 4.15.0-23-generic
machine    : x86_64
processor  : x86_64
CPU cores  : 4
interpreter: 64bit
Git hash   : 3f62e1141f480e65f685cabe81f7a79b17bcd06d

In :
from sympy import *

In :
init_printing(use_latex='mathjax')


Just a small introduction to SymPy: it works by using Python expressions on formal variables. First, we define variables, and then we can solve linear equations for example:

In :
var("x y z")

Out:
$$\left ( x, \quad y, \quad z\right )$$
In :
solve(x + 2 + 19)

Out:
$$\left [ -21\right ]$$
In :
solve(x + 2*y + 19)

Out:
$$\left [ \left \{ x : - 2 y - 19\right \}\right ]$$

# Computations¶

We start by defining $b_1 > 0$ and the other variables.

In :
b1 = symbols("b1", positive=True)


We need a lot of variables.

In :
var("h1 h2 h3 h4 k131 k132 k141 k142 k231 k232 k241 k242 rest1 rest2 t s a b");


Then we can start to follow Ludovic's notebook and define the first expressions.

In :
J = Matrix([[0, -b1], [b1, 0]])
J

Out:
$$\left[\begin{matrix}0 & - b_{1}\\b_{1} & 0\end{matrix}\right]$$
In :
h0 = Matrix([h1, h2])
h0

Out:
$$\left[\begin{matrix}h_{1}\\h_{2}\end{matrix}\right]$$

$\hat{h}$ is defined as $t \mapsto \exp(t J) . h_0$ and $\hat{x}$ as $t \mapsto \int_{s=0}^{s=t} \hat{h}(s) \mathrm{d}s$.

In :
hhat = Lambda(t, (t * J).exp() * h0)
hhat

Out:
$$\left( t \mapsto \left[\begin{matrix}h_{1} \left(\frac{1}{2} e^{i b_{1} t} + \frac{1}{2} e^{- i b_{1} t}\right) + h_{2} \left(\frac{i}{2} e^{i b_{1} t} - \frac{i}{2} e^{- i b_{1} t}\right)\\h_{1} \left(- \frac{i}{2} e^{i b_{1} t} + \frac{i}{2} e^{- i b_{1} t}\right) + h_{2} \left(\frac{1}{2} e^{i b_{1} t} + \frac{1}{2} e^{- i b_{1} t}\right)\end{matrix}\right] \right)$$
In :
xhat = Lambda(t, integrate(hhat(s), (s, 0, t)))
xhat

Out:
$$\left( t \mapsto \left[\begin{matrix}- \frac{h_{2}}{b_{1}} + \frac{1}{4 b_{1}^{2}} \left(\left(- 2 i b_{1} h_{1} + 2 b_{1} h_{2}\right) e^{i b_{1} t} + \left(2 i b_{1} h_{1} + 2 b_{1} h_{2}\right) e^{- i b_{1} t}\right)\\\frac{h_{1}}{b_{1}} + \frac{1}{4 b_{1}^{2}} \left(\left(- 2 b_{1} h_{1} - 2 i b_{1} h_{2}\right) e^{i b_{1} t} + \left(- 2 b_{1} h_{1} + 2 i b_{1} h_{2}\right) e^{- i b_{1} t}\right)\end{matrix}\right] \right)$$

$\hat{z}$ is slightly more complex:

In :
integrand = simplify(b1 * (hhat(s) * xhat(s) - hhat(s) * xhat(s) )/2)
integrand

Out:
$$\frac{h_{1}^{2}}{4} e^{i b_{1} s} - \frac{h_{1}^{2}}{2} + \frac{h_{1}^{2}}{4} e^{- i b_{1} s} + \frac{h_{2}^{2}}{4} e^{i b_{1} s} - \frac{h_{2}^{2}}{2} + \frac{h_{2}^{2}}{4} e^{- i b_{1} s}$$
In :
zhat = Lambda(t, integrate(integrand, (s, 0, t)))
zhat

Out:
$$\left( t \mapsto t \left(- \frac{h_{1}^{2}}{2} - \frac{h_{2}^{2}}{2}\right) + \frac{1}{16 b_{1}^{2}} \left(\left(- 4 i b_{1} h_{1}^{2} - 4 i b_{1} h_{2}^{2}\right) e^{i b_{1} t} + \left(4 i b_{1} h_{1}^{2} + 4 i b_{1} h_{2}^{2}\right) e^{- i b_{1} t}\right) \right)$$

As we asked $b_1 > 0$, this won't get too complicated:

In :
factor(zhat(t))

Out:
$$- \frac{e^{- i b_{1} t}}{4 b_{1}} \left(h_{1}^{2} + h_{2}^{2}\right) \left(2 b_{1} t e^{i b_{1} t} + i e^{2 i b_{1} t} - i\right)$$

Two more expressions:

In :
exp21 = (3 * a * h1**2 + a * h2**2 + 2 * b * h1 * h2 + k131 * h1 * h3 + k141 * h1 * h4
+ k231 * h2 * h3 + k241 * h2 * h4 + rest1(h3, h4))

In :
exp22 = (b * h1**2 + 3 * b * h2**2 + 2 * a * h1 * h2 + k132 * h1 * h3 + k142 * h1 * h4
+ k232 * h2 * h3 + k242 * h2 * h4 + rest2(h3, h4))


Both terms rest1 and rest2 are not important, they only depend on $h_3$ and $h_4$ and will vanish as soon as we only differentiate with respect to $h_1$ and $h_2$.

So far so good. Next cell:

In :
C10 = 2 * (pi / b1) * Matrix([h1, h2])
C10

Out:
$$\left[\begin{matrix}\frac{2 \pi}{b_{1}} h_{1}\\\frac{2 \pi}{b_{1}} h_{2}\end{matrix}\right]$$
In :
A0 = simplify(Matrix([
[diff(exp21, h1), diff(exp21, h2)],
[diff(exp22, h1), diff(exp22, h2)],
]))
A0

Out:
$$\left[\begin{matrix}6 a h_{1} + 2 b h_{2} + h_{3} k_{131} + h_{4} k_{141} & 2 a h_{2} + 2 b h_{1} + h_{3} k_{231} + h_{4} k_{241}\\2 a h_{2} + 2 b h_{1} + h_{3} k_{132} + h_{4} k_{142} & 2 a h_{1} + 6 b h_{2} + h_{3} k_{232} + h_{4} k_{242}\end{matrix}\right]$$

rest1 and rest2 already vanished.

In :
jacobian = simplify(Matrix([
[diff(exp21, h1), diff(exp21, h2), C10],
[diff(exp22, h1), diff(exp22, h2), C10],
[diff(zhat(2 * pi / b1), h1), diff(zhat(2 * pi / b1), h2), 0],
]))
jacobian

Out:
$$\left[\begin{matrix}6 a h_{1} + 2 b h_{2} + h_{3} k_{131} + h_{4} k_{141} & 2 a h_{2} + 2 b h_{1} + h_{3} k_{231} + h_{4} k_{241} & \frac{2 \pi}{b_{1}} h_{1}\\2 a h_{2} + 2 b h_{1} + h_{3} k_{132} + h_{4} k_{142} & 2 a h_{1} + 6 b h_{2} + h_{3} k_{232} + h_{4} k_{242} & \frac{2 \pi}{b_{1}} h_{2}\\- \frac{2 \pi}{b_{1}} h_{1} & - \frac{2 \pi}{b_{1}} h_{2} & 0\end{matrix}\right]$$

We need one more variable to solve an equation.

In :
var('dt')

Out:
$$dt$$
In :
tc = factor(solve(Equality((jacobian + Matrix([[dt,0,0], [0,dt,0], [0,0,0]])).det(), 0), dt))

In :
tc

Out:
$$- \frac{1}{h_{1}^{2} + h_{2}^{2}} \left(2 a h_{1}^{3} + 2 a h_{1} h_{2}^{2} + 2 b h_{1}^{2} h_{2} + 2 b h_{2}^{3} + h_{1}^{2} h_{3} k_{232} + h_{1}^{2} h_{4} k_{242} - h_{1} h_{2} h_{3} k_{132} - h_{1} h_{2} h_{3} k_{231} - h_{1} h_{2} h_{4} k_{142} - h_{1} h_{2} h_{4} k_{241} + h_{2}^{2} h_{3} k_{131} + h_{2}^{2} h_{4} k_{141}\right)$$

I can compare by copying the result from the document:

In :
tc2 = (-1 / (h1**2 + h2**2)) * (
2 * a * h1**3
+ 6 * a * h1**2 * h2
- 4 * b * h1**2 * h2
+ 2 * a * h1 * h2**2
+ 2 * b * h2**3
+ h2**2 * h3 * k131
- h1 * h2 * h3 * k132
+ h2**2 * h4 * k141
- h1 * h2 * h4 * k142
- h1 * h2 * h3 * k231
+ h1**2 * h3 * k232
- h1 * h2 * h4 * k241
+ h1**2 * h4 * k242 )


Drat, they are not equal! We might have a mistake somewhere, even though I just CAN'T find it!

In :
simplify(tc - tc2)

Out:
$$\frac{6 h_{1}^{2} h_{2} \left(a - b\right)}{h_{1}^{2} + h_{2}^{2}}$$

Let's use the one from Ludovic's notebook.

In :
tc = tc2


Next cell.

In :
A12 = simplify(A0 + Matrix([[tc, 0], [0, tc]]))

In :
var("u1 u2 u5")

Out:
$$\left ( u_{1}, \quad u_{2}, \quad u_{5}\right )$$
In :
Psi = Lambda((u1, u2, u5), simplify(
u5 * Matrix(A12.dot(Matrix([h1, h2]))).dot(Matrix([h2, -h1])) + 2 * pi / b1 * ( h1**2 + h2**2) * (h1 * u2 - h2 * u1)
))

In :
my_psi = factor(simplify(Psi(u1, u2, u5)))
my_psi

Out:
$$- \frac{1}{b_{1}} \left(- 2 a b_{1} h_{1}^{2} h_{2} u_{5} - 2 a b_{1} h_{2}^{3} u_{5} + 2 b b_{1} h_{1}^{3} u_{5} + 2 b b_{1} h_{1} h_{2}^{2} u_{5} + b_{1} h_{1}^{2} h_{3} k_{132} u_{5} + b_{1} h_{1}^{2} h_{4} k_{142} u_{5} - b_{1} h_{1} h_{2} h_{3} k_{131} u_{5} + b_{1} h_{1} h_{2} h_{3} k_{232} u_{5} - b_{1} h_{1} h_{2} h_{4} k_{141} u_{5} + b_{1} h_{1} h_{2} h_{4} k_{242} u_{5} - b_{1} h_{2}^{2} h_{3} k_{231} u_{5} - b_{1} h_{2}^{2} h_{4} k_{241} u_{5} - 2 \pi h_{1}^{3} u_{2} + 2 \pi h_{1}^{2} h_{2} u_{1} - 2 \pi h_{1} h_{2}^{2} u_{2} + 2 \pi h_{2}^{3} u_{1}\right)$$

Let's compare with the value from Ludovic's notebook:

In :
his_psi = (1 / b1) * (
2 * h1**3 * (pi * u2 - b * b1 * u5 ) -
h1**2 * (2 * h2 * pi * u1 - 2 * a * b1 * h2 * u5 + b1 * h3 * k132 * u5 +
b1 * h4 * k142 * u5 ) +
h2**2 * ( -2 * h2 * pi * u1 + 2 * a * b1 * h2 * u5 + b1 * h3 * k231 * u5 +
b1 * h4 * k241 * u5 ) +
h1 * h2 * (b1 * (h3 * k131 + h4 * k141 - h3 * k232 - h4 * k242) * u5 +
2 * h2 * (pi * u2 - b * b1 * u5))
)

In :
simplify(my_psi - his_psi)

Out:
$$0$$

Ok, we have the same result so far! Great!

Next cell.

In :
var("nu0")

Out:
$$\nu_{0}$$

Here again, Sympy fails to solve the equation. It can solve on both lines separately, but cannot find a solution that satisfies both lines.

In :
s1 = solve(Eq((Matrix(A12.dot(Matrix([-h2, h1]))) + nu0 * C10)), nu0)

In :
s2 = solve(Eq((Matrix(A12.dot(Matrix([-h2, h1]))) + nu0 * C10)), nu0)

In :
solve(Eq(s1, s2))

Out:
$$\left [ \left \{ a : b\right \}, \quad \left \{ h_{1} : 0\right \}\right ]$$

It is not possible to impose these constraints. And Sympy fails to solve both equation simultaneously:

In :
solve(Eq(Matrix(A12.dot(Matrix([-h2, h1]))) + nu0 * C10, Matrix([0,0])), nu0)

Out:
$$\left [ \right ]$$
In :
nu = simplify(solve(Eq(Matrix(A12.dot(Matrix([-h2, h1]))) + nu0 * C10, Matrix([0,0])), nu0))
if nu == []:
nu = var("nu")


We can continue by using a formal variable for $\nu$ (nu).

In :
v = Matrix([nu, -h2, h1])

In :
v

Out:
$$\left[\begin{matrix}\nu\\- h_{2}\\h_{1}\end{matrix}\right]$$

Let's finish.

In :
var("eta f F");

In :
d = [
lambda F: -eta**2 * diff(F, eta) + eta * t * diff(F, t),
lambda F: diff(F, h1),
lambda F: diff(F, h2)
]

In :
d

Out:
[<function __main__.<lambda>(F)>,
<function __main__.<lambda>(F)>,
<function __main__.<lambda>(F)>]

Next cell.

In :
var("i1 i2 g1 g2 dt1 dt2");

In :
g = eta * g1(t + eta * dt1, h1, h2) + eta**2 * g2(t, h1, h2)
g

Out:
$$\eta^{2} \operatorname{g_{2}}{\left (t,h_{1},h_{2} \right )} + \eta \operatorname{g_{1}}{\left (dt_{1} \eta + t,h_{1},h_{2} \right )}$$

We have a sum to compute:

In :
res = 0
for i1 in range(0, 3):
for i2 in range(0, 3):
res += v[i1] * v[i2] * d[i1](d[i2](g))

In :
simplify(res)

Out:
$$\eta \left(\eta^{2} \nu^{2} \left(2 dt_{1} \eta \left. \frac{\partial}{\partial \xi_{1}} \operatorname{g_{1}}{\left (\xi_{1},h_{1},h_{2} \right )} \right|_{\substack{ \xi_{1}=dt_{1} \eta + t }} + \eta \left(dt_{1}^{2} \eta \left. \frac{\partial^{2}}{\partial \xi_{1}^{2}} \operatorname{g_{1}}{\left (\xi_{1},h_{1},h_{2} \right )} \right|_{\substack{ \xi_{1}=dt_{1} \eta + t }} + 2 dt_{1} \left. \frac{\partial}{\partial \xi_{1}} \operatorname{g_{1}}{\left (\xi_{1},h_{1},h_{2} \right )} \right|_{\substack{ \xi_{1}=dt_{1} \eta + t }} + 2 \operatorname{g_{2}}{\left (t,h_{1},h_{2} \right )}\right) + 4 \eta \operatorname{g_{2}}{\left (t,h_{1},h_{2} \right )} - t \left(\eta \frac{\partial}{\partial t} \operatorname{g_{2}}{\left (t,h_{1},h_{2} \right )} + \left. \frac{\partial}{\partial \xi_{1}} \operatorname{g_{1}}{\left (\xi_{1},h_{1},h_{2} \right )} \right|_{\substack{ \xi_{1}=dt_{1} \eta + t }}\right) + t \left(- dt_{1} \eta \left. \frac{\partial^{2}}{\partial \xi_{1}^{2}} \operatorname{g_{1}}{\left (\xi_{1},h_{1},h_{2} \right )} \right|_{\substack{ \xi_{1}=dt_{1} \eta + t }} - \eta \frac{\partial}{\partial t} \operatorname{g_{2}}{\left (t,h_{1},h_{2} \right )} + t \left(\eta \frac{\partial^{2}}{\partial t^{2}} \operatorname{g_{2}}{\left (t,h_{1},h_{2} \right )} + \left. \frac{\partial^{2}}{\partial \xi_{1}^{2}} \operatorname{g_{1}}{\left (\xi_{1},h_{1},h_{2} \right )} \right|_{\substack{ \xi_{1}=dt_{1} \eta + t }}\right)\right) - t \left(dt_{1} \eta \left. \frac{\partial^{2}}{\partial \xi_{1}^{2}} \operatorname{g_{1}}{\left (\xi_{1},h_{1},h_{2} \right )} \right|_{\substack{ \xi_{1}=dt_{1} \eta + t }} + 2 \eta \frac{\partial}{\partial t} \operatorname{g_{2}}{\left (t,h_{1},h_{2} \right )} + \left. \frac{\partial}{\partial \xi_{1}} \operatorname{g_{1}}{\left (\xi_{1},h_{1},h_{2} \right )} \right|_{\substack{ \xi_{1}=dt_{1} \eta + t }}\right) + 2 \operatorname{g_{1}}{\left (dt_{1} \eta + t,h_{1},h_{2} \right )}\right) - 2 \eta h_{1} \nu \left(dt_{1} \eta \left. \frac{\partial^{2}}{\partial \xi_{1}\partial h_{2}} \operatorname{g_{1}}{\left (\xi_{1},h_{1},h_{2} \right )} \right|_{\substack{ \xi_{1}=dt_{1} \eta + t }} + 2 \eta \frac{\partial}{\partial h_{2}} \operatorname{g_{2}}{\left (t,h_{1},h_{2} \right )} - t \left(\eta \frac{\partial^{2}}{\partial h_{2}\partial t} \operatorname{g_{2}}{\left (t,h_{1},h_{2} \right )} + \left. \frac{\partial^{2}}{\partial \xi_{1}\partial h_{2}} \operatorname{g_{1}}{\left (\xi_{1},h_{1},h_{2} \right )} \right|_{\substack{ \xi_{1}=dt_{1} \eta + t }}\right) + \frac{\partial}{\partial h_{2}} \operatorname{g_{1}}{\left (dt_{1} \eta + t,h_{1},h_{2} \right )}\right) + 2 \eta h_{2} \nu \left(dt_{1} \eta \left. \frac{\partial^{2}}{\partial \xi_{1}\partial h_{1}} \operatorname{g_{1}}{\left (\xi_{1},h_{1},h_{2} \right )} \right|_{\substack{ \xi_{1}=dt_{1} \eta + t }} + 2 \eta \frac{\partial}{\partial h_{1}} \operatorname{g_{2}}{\left (t,h_{1},h_{2} \right )} - t \left(\eta \frac{\partial^{2}}{\partial h_{1}\partial t} \operatorname{g_{2}}{\left (t,h_{1},h_{2} \right )} + \left. \frac{\partial^{2}}{\partial \xi_{1}\partial h_{1}} \operatorname{g_{1}}{\left (\xi_{1},h_{1},h_{2} \right )} \right|_{\substack{ \xi_{1}=dt_{1} \eta + t }}\right) + \frac{\partial}{\partial h_{1}} \operatorname{g_{1}}{\left (dt_{1} \eta + t,h_{1},h_{2} \right )}\right) + h_{1}^{2} \left(\eta \frac{\partial^{2}}{\partial h_{2}^{2}} \operatorname{g_{2}}{\left (t,h_{1},h_{2} \right )} + \frac{\partial^{2}}{\partial h_{2}^{2}} \operatorname{g_{1}}{\left (dt_{1} \eta + t,h_{1},h_{2} \right )}\right) - 2 h_{1} h_{2} \left(\eta \frac{\partial^{2}}{\partial h_{1}\partial h_{2}} \operatorname{g_{2}}{\left (t,h_{1},h_{2} \right )} + \frac{\partial^{2}}{\partial h_{1}\partial h_{2}} \operatorname{g_{1}}{\left (dt_{1} \eta + t,h_{1},h_{2} \right )}\right) + h_{2}^{2} \left(\eta \frac{\partial^{2}}{\partial h_{1}^{2}} \operatorname{g_{2}}{\left (t,h_{1},h_{2} \right )} + \frac{\partial^{2}}{\partial h_{1}^{2}} \operatorname{g_{1}}{\left (dt_{1} \eta + t,h_{1},h_{2} \right )}\right)\right)$$

Then a double derivative

In :
d2f = simplify((diff(res, eta, eta) / 2).subs(eta, 0))


Now, we should replace $g_1$ and $g_2$ with the following expression:

In :
g1 = Lambda((t, h1, h2), simplify(Matrix(list(xhat(t)) + )))

In :
g1(t, h1, h2)

Out:
$$\left[\begin{matrix}\frac{e^{- i b_{1} t}}{2 b_{1}} \left(i h_{1} - 2 h_{2} e^{i b_{1} t} + h_{2} + \left(- i h_{1} + h_{2}\right) e^{2 i b_{1} t}\right)\\\frac{e^{- i b_{1} t}}{2 b_{1}} \left(2 h_{1} e^{i b_{1} t} - h_{1} + i h_{2} - \left(h_{1} + i h_{2}\right) e^{2 i b_{1} t}\right)\\0\end{matrix}\right]$$
In :
g2 = Lambda((t, h1, h2), simplify(Matrix([exp21, exp22] + [zhat(t)])))

In :
g2(t, h1, h2)

Out:
$$\left[\begin{matrix}3 a h_{1}^{2} + a h_{2}^{2} + 2 b h_{1} h_{2} + h_{1} h_{3} k_{131} + h_{1} h_{4} k_{141} + h_{2} h_{3} k_{231} + h_{2} h_{4} k_{241} + \operatorname{rest_{1}}{\left (h_{3},h_{4} \right )}\\2 a h_{1} h_{2} + b h_{1}^{2} + 3 b h_{2}^{2} + h_{1} h_{3} k_{132} + h_{1} h_{4} k_{142} + h_{2} h_{3} k_{232} + h_{2} h_{4} k_{242} + \operatorname{rest_{2}}{\left (h_{3},h_{4} \right )}\\\frac{e^{- i b_{1} t}}{4 b_{1}} \left(- 2 b_{1} t \left(h_{1}^{2} + h_{2}^{2}\right) e^{i b_{1} t} + i h_{1}^{2} + i h_{2}^{2} - i \left(h_{1}^{2} + h_{2}^{2}\right) e^{2 i b_{1} t}\right)\end{matrix}\right]$$

Let's see if the replacement is possible:

In :
simplify(d2f)

Out:
$$h_{1}^{2} \left(dt_{1} \left. \frac{\partial^{3}}{\partial \xi_{1}\partial h_{2}^{2}} \operatorname{g_{1}}{\left (\xi_{1},h_{1},h_{2} \right )} \right|_{\substack{ \xi_{1}=t }} + \frac{\partial^{2}}{\partial h_{2}^{2}} \operatorname{g_{2}}{\left (t,h_{1},h_{2} \right )}\right) - 2 h_{1} h_{2} \left(dt_{1} \left. \frac{\partial^{3}}{\partial \xi_{1}\partial h_{1}\partial h_{2}} \operatorname{g_{1}}{\left (\xi_{1},h_{1},h_{2} \right )} \right|_{\substack{ \xi_{1}=t }} + \frac{\partial^{2}}{\partial h_{1}\partial h_{2}} \operatorname{g_{2}}{\left (t,h_{1},h_{2} \right )}\right) + 2 h_{1} \nu \left(t \left. \frac{\partial^{2}}{\partial \xi_{1}\partial h_{2}} \operatorname{g_{1}}{\left (\xi_{1},h_{1},h_{2} \right )} \right|_{\substack{ \xi_{1}=t }} - \frac{\partial}{\partial h_{2}} \operatorname{g_{1}}{\left (t,h_{1},h_{2} \right )}\right) + h_{2}^{2} \left(dt_{1} \left. \frac{\partial^{3}}{\partial \xi_{1}\partial h_{1}^{2}} \operatorname{g_{1}}{\left (\xi_{1},h_{1},h_{2} \right )} \right|_{\substack{ \xi_{1}=t }} + \frac{\partial^{2}}{\partial h_{1}^{2}} \operatorname{g_{2}}{\left (t,h_{1},h_{2} \right )}\right) - 2 h_{2} \nu \left(t \left. \frac{\partial^{2}}{\partial \xi_{1}\partial h_{1}} \operatorname{g_{1}}{\left (\xi_{1},h_{1},h_{2} \right )} \right|_{\substack{ \xi_{1}=t }} - \frac{\partial}{\partial h_{1}} \operatorname{g_{1}}{\left (t,h_{1},h_{2} \right )}\right)$$
In :
simplify(d2f.subs({
g1: Lambda((t, h1, h2), simplify(Matrix(list(xhat(t)) + ))),
g2: Lambda((t, h1, h2), simplify(Matrix([exp21, exp22] +[zhat(t)]))),
}))

Out:
$$h_{1}^{2} \left(dt_{1} \left. \frac{\partial^{3}}{\partial \xi_{1}\partial h_{2}^{2}} \operatorname{g_{1}}{\left (\xi_{1},h_{1},h_{2} \right )} \right|_{\substack{ \xi_{1}=t }} + \frac{\partial^{2}}{\partial h_{2}^{2}} \operatorname{g_{2}}{\left (t,h_{1},h_{2} \right )}\right) - 2 h_{1} h_{2} \left(dt_{1} \left. \frac{\partial^{3}}{\partial \xi_{1}\partial h_{1}\partial h_{2}} \operatorname{g_{1}}{\left (\xi_{1},h_{1},h_{2} \right )} \right|_{\substack{ \xi_{1}=t }} + \frac{\partial^{2}}{\partial h_{1}\partial h_{2}} \operatorname{g_{2}}{\left (t,h_{1},h_{2} \right )}\right) + 2 h_{1} \nu \left(t \left. \frac{\partial^{2}}{\partial \xi_{1}\partial h_{2}} \operatorname{g_{1}}{\left (\xi_{1},h_{1},h_{2} \right )} \right|_{\substack{ \xi_{1}=t }} - \frac{\partial}{\partial h_{2}} \operatorname{g_{1}}{\left (t,h_{1},h_{2} \right )}\right) + h_{2}^{2} \left(dt_{1} \left. \frac{\partial^{3}}{\partial \xi_{1}\partial h_{1}^{2}} \operatorname{g_{1}}{\left (\xi_{1},h_{1},h_{2} \right )} \right|_{\substack{ \xi_{1}=t }} + \frac{\partial^{2}}{\partial h_{1}^{2}} \operatorname{g_{2}}{\left (t,h_{1},h_{2} \right )}\right) - 2 h_{2} \nu \left(t \left. \frac{\partial^{2}}{\partial \xi_{1}\partial h_{1}} \operatorname{g_{1}}{\left (\xi_{1},h_{1},h_{2} \right )} \right|_{\substack{ \xi_{1}=t }} - \frac{\partial}{\partial h_{1}} \operatorname{g_{1}}{\left (t,h_{1},h_{2} \right )}\right)$$

OK. This failed. But we can copy and paste this and the replacement of $g_1$ and $g_2$ will be effective:

In :
d2f = (
h1**2*(dt1*diff(g1(t, h1, h2), t, h2, h2)
+ diff(g2(t, h1, h2), h2, h2))
- 2*h1*h2*(dt1*diff(g1(t, h1, h2), t, h1, h2)
+ diff(g2(t, h1, h2), h1, h2))
+ 2*h1*nu*(t*diff(g1(t, h1, h2), t, h2)
- diff(g1(t, h1, h2), h2))
+ h2**2*(dt1*diff(g1(t, h1, h2), t, h1, h1)
+ diff(g2(t, h1, h2), h1, h1))
- 2*h2*nu*(t*diff(g1(t, h1, h2), t, h1)
- diff(g1(t, h1, h2), h1))
)


And the same for $t$.

In :
t = 2 * pi / b1

In :
d2f

Out:
$$\left[\begin{matrix}2 a h_{1}^{2} + 6 a h_{2}^{2} - 4 b h_{1} h_{2} + 2 h_{1} \nu \left(i t \left(- \frac{1}{2} \left(e^{2 i b_{1} t} - 2 e^{i b_{1} t} + 1\right) e^{- i b_{1} t} + e^{i b_{1} t} - 1\right) - \frac{e^{- i b_{1} t}}{2 b_{1}} \left(e^{2 i b_{1} t} - 2 e^{i b_{1} t} + 1\right)\right) - 2 h_{2} \nu \left(t \left(- \frac{1}{2} \left(e^{2 i b_{1} t} - 1\right) e^{- i b_{1} t} + e^{i b_{1} t}\right) - \frac{e^{- i b_{1} t}}{2 b_{1}} \left(- i e^{2 i b_{1} t} + i\right)\right)\\- 4 a h_{1} h_{2} + 6 b h_{1}^{2} + 2 b h_{2}^{2} + 2 h_{1} \nu \left(t \left(- \frac{1}{2} \left(e^{2 i b_{1} t} - 1\right) e^{- i b_{1} t} + e^{i b_{1} t}\right) - \frac{e^{- i b_{1} t}}{2 b_{1}} \left(- i e^{2 i b_{1} t} + i\right)\right) - 2 h_{2} \nu \left(i t \left(\frac{1}{2} \left(e^{2 i b_{1} t} - 2 e^{i b_{1} t} + 1\right) e^{- i b_{1} t} - e^{i b_{1} t} + 1\right) - \frac{e^{- i b_{1} t}}{2 b_{1}} \left(- e^{2 i b_{1} t} + 2 e^{i b_{1} t} - 1\right)\right)\\- \frac{h_{1}^{2}}{2 b_{1}} \left(2 b_{1} t e^{i b_{1} t} + i e^{2 i b_{1} t} - i\right) e^{- i b_{1} t} - \frac{h_{2}^{2}}{2 b_{1}} \left(2 b_{1} t e^{i b_{1} t} + i e^{2 i b_{1} t} - i\right) e^{- i b_{1} t}\end{matrix}\right]$$

The replacement does not work apparently, so let's do it manually:

In :
d2f = Matrix([
[ 2*a*h1**2 + 6*a*h2**2 - 4*b*h1*h2 + 2*h1*nu*(I*t*(-(exp(2*I*b1*t) - 2*exp(I*b1*t) + 1)*exp(-I*b1*t)/2 + exp(I*b1*t) - 1) - (exp(2*I*b1*t) - 2*exp(I*b1*t) + 1)*exp(-I*b1*t)/(2*b1)) - 2*h2*nu*(t*(-(exp(2*I*b1*t) - 1)*exp(-I*b1*t)/2 + exp(I*b1*t)) - (-I*exp(2*I*b1*t) + I)*exp(-I*b1*t)/(2*b1))],
[-4*a*h1*h2 + 6*b*h1**2 + 2*b*h2**2 + 2*h1*nu*(t*(-(exp(2*I*b1*t) - 1)*exp(-I*b1*t)/2 + exp(I*b1*t)) - (-I*exp(2*I*b1*t) + I)*exp(-I*b1*t)/(2*b1)) - 2*h2*nu*(I*t*((exp(2*I*b1*t) - 2*exp(I*b1*t) + 1)*exp(-I*b1*t)/2 - exp(I*b1*t) + 1) - (-exp(2*I*b1*t) + 2*exp(I*b1*t) - 1)*exp(-I*b1*t)/(2*b1))],
[                                                                                                                                                       -h1**2*(2*b1*t*exp(I*b1*t) + I*exp(2*I*b1*t) - I)*exp(-I*b1*t)/(2*b1) - h2**2*(2*b1*t*exp(I*b1*t) + I*exp(2*I*b1*t) - I)*exp(-I*b1*t)/(2*b1)]])

In :
d2Exp = simplify(d2f)


It still works well. And we removed the complex exponential, this result is purely real now!

In :
d2Exp

Out:
$$\left[\begin{matrix}2 a h_{1}^{2} + 6 a h_{2}^{2} - 4 b h_{1} h_{2} - \frac{4 \pi}{b_{1}} h_{2} \nu\\- 4 a h_{1} h_{2} + 6 b h_{1}^{2} + 2 b h_{2}^{2} + \frac{4 \pi}{b_{1}} h_{1} \nu\\- \frac{2 \pi}{b_{1}} \left(h_{1}^{2} + h_{2}^{2}\right)\end{matrix}\right]$$
In :
Psi

Out:
$$\left( \left ( u_{1}, \quad u_{2}, \quad u_{5}\right ) \mapsto 2 a h_{1}^{2} h_{2} u_{5} + 2 a h_{2}^{3} u_{5} - 2 b h_{1}^{3} u_{5} - 2 b h_{1} h_{2}^{2} u_{5} - h_{1}^{2} h_{3} k_{132} u_{5} - h_{1}^{2} h_{4} k_{142} u_{5} + h_{1} h_{2} h_{3} k_{131} u_{5} - h_{1} h_{2} h_{3} k_{232} u_{5} + h_{1} h_{2} h_{4} k_{141} u_{5} - h_{1} h_{2} h_{4} k_{242} u_{5} + h_{2}^{2} h_{3} k_{231} u_{5} + h_{2}^{2} h_{4} k_{241} u_{5} + \frac{2 \pi}{b_{1}} h_{1}^{3} u_{2} - \frac{2 \pi}{b_{1}} h_{1}^{2} h_{2} u_{1} + \frac{2 \pi}{b_{1}} h_{1} h_{2}^{2} u_{2} - \frac{2 \pi}{b_{1}} h_{2}^{3} u_{1} \right)$$

So now we can call $\Psi$ on the three components of this d2Exp :

In :
Psi(d2Exp, d2Exp, d2Exp)

Out:
$$- \frac{4 \pi}{b_{1}} a h_{1}^{2} h_{2} \left(h_{1}^{2} + h_{2}^{2}\right) - \frac{4 \pi}{b_{1}} a h_{2}^{3} \left(h_{1}^{2} + h_{2}^{2}\right) + \frac{4 \pi}{b_{1}} b h_{1}^{3} \left(h_{1}^{2} + h_{2}^{2}\right) + \frac{4 \pi}{b_{1}} b h_{1} h_{2}^{2} \left(h_{1}^{2} + h_{2}^{2}\right) + \frac{2 \pi}{b_{1}} h_{1}^{3} \left(- 4 a h_{1} h_{2} + 6 b h_{1}^{2} + 2 b h_{2}^{2} + \frac{4 \pi}{b_{1}} h_{1} \nu\right) - \frac{2 \pi}{b_{1}} h_{1}^{2} h_{2} \left(2 a h_{1}^{2} + 6 a h_{2}^{2} - 4 b h_{1} h_{2} - \frac{4 \pi}{b_{1}} h_{2} \nu\right) + \frac{2 \pi}{b_{1}} h_{1}^{2} h_{3} k_{132} \left(h_{1}^{2} + h_{2}^{2}\right) + \frac{2 \pi}{b_{1}} h_{1}^{2} h_{4} k_{142} \left(h_{1}^{2} + h_{2}^{2}\right) + \frac{2 \pi}{b_{1}} h_{1} h_{2}^{2} \left(- 4 a h_{1} h_{2} + 6 b h_{1}^{2} + 2 b h_{2}^{2} + \frac{4 \pi}{b_{1}} h_{1} \nu\right) - \frac{2 \pi}{b_{1}} h_{1} h_{2} h_{3} k_{131} \left(h_{1}^{2} + h_{2}^{2}\right) + \frac{2 \pi}{b_{1}} h_{1} h_{2} h_{3} k_{232} \left(h_{1}^{2} + h_{2}^{2}\right) - \frac{2 \pi}{b_{1}} h_{1} h_{2} h_{4} k_{141} \left(h_{1}^{2} + h_{2}^{2}\right) + \frac{2 \pi}{b_{1}} h_{1} h_{2} h_{4} k_{242} \left(h_{1}^{2} + h_{2}^{2}\right) - \frac{2 \pi}{b_{1}} h_{2}^{3} \left(2 a h_{1}^{2} + 6 a h_{2}^{2} - 4 b h_{1} h_{2} - \frac{4 \pi}{b_{1}} h_{2} \nu\right) - \frac{2 \pi}{b_{1}} h_{2}^{2} h_{3} k_{231} \left(h_{1}^{2} + h_{2}^{2}\right) - \frac{2 \pi}{b_{1}} h_{2}^{2} h_{4} k_{241} \left(h_{1}^{2} + h_{2}^{2}\right)$$
In :
simplify(expand(Psi(d2Exp, d2Exp, d2Exp) / (1 / (1/b1 * 2 * (h1**2 + h2**2) * pi))))

Out:
$$\frac{4 \pi^{2}}{b_{1}^{3}} \left(b_{1} \left(- 8 a h_{1}^{6} h_{2} - 24 a h_{1}^{4} h_{2}^{3} - 24 a h_{1}^{2} h_{2}^{5} - 8 a h_{2}^{7} + 8 b h_{1}^{7} + 24 b h_{1}^{5} h_{2}^{2} + 24 b h_{1}^{3} h_{2}^{4} + 8 b h_{1} h_{2}^{6} + h_{1}^{6} h_{3} k_{132} + h_{1}^{6} h_{4} k_{142} - h_{1}^{5} h_{2} h_{3} k_{131} + h_{1}^{5} h_{2} h_{3} k_{232} - h_{1}^{5} h_{2} h_{4} k_{141} + h_{1}^{5} h_{2} h_{4} k_{242} + 2 h_{1}^{4} h_{2}^{2} h_{3} k_{132} - h_{1}^{4} h_{2}^{2} h_{3} k_{231} + 2 h_{1}^{4} h_{2}^{2} h_{4} k_{142} - h_{1}^{4} h_{2}^{2} h_{4} k_{241} - 2 h_{1}^{3} h_{2}^{3} h_{3} k_{131} + 2 h_{1}^{3} h_{2}^{3} h_{3} k_{232} - 2 h_{1}^{3} h_{2}^{3} h_{4} k_{141} + 2 h_{1}^{3} h_{2}^{3} h_{4} k_{242} + h_{1}^{2} h_{2}^{4} h_{3} k_{132} - 2 h_{1}^{2} h_{2}^{4} h_{3} k_{231} + h_{1}^{2} h_{2}^{4} h_{4} k_{142} - 2 h_{1}^{2} h_{2}^{4} h_{4} k_{241} - h_{1} h_{2}^{5} h_{3} k_{131} + h_{1} h_{2}^{5} h_{3} k_{232} - h_{1} h_{2}^{5} h_{4} k_{141} + h_{1} h_{2}^{5} h_{4} k_{242} - h_{2}^{6} h_{3} k_{231} - h_{2}^{6} h_{4} k_{241}\right) + 4 \pi \nu \left(h_{1}^{6} + 3 h_{1}^{4} h_{2}^{2} + 3 h_{1}^{2} h_{2}^{4} + h_{2}^{6}\right)\right)$$

We don't have the value for nu !

## Conclusion¶

We do NOT obtain the same result as the document. Everything failed at the end.

Too bad, but still, it was interesting. I guess?

See here for other notebooks I wrote.