Table of Contents¶

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 [1]:

%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 [2]:

from sympy import *

In [3]:

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 [4]:

var("x y z")

Out[4]:

$$\left ( x, \quad y, \quad z\right )$$

In [5]:

solve(x + 2 + 19)

Out[5]:

$$\left [ -21\right ]$$

In [6]:

solve(x + 2*y + 19)

Out[6]:

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

Computations¶

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

In [7]:

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

We need a lot of variables.

In [8]:

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 [13]:

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

Out[13]:

$$\left[\begin{matrix}0 & - b_{1}\\b_{1} & 0\end{matrix}\right]$$

In [14]:

h0 = Matrix([h1, h2])
h0

Out[14]:

$$\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 [15]:

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

Out[15]:

$$\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 [16]:

xhat = Lambda(t, integrate(hhat(s), (s, 0, t)))
xhat

Out[16]:

$$\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 [20]:

integrand = simplify(b1 * (hhat(s)[0] * xhat(s)[1] - hhat(s)[1] * xhat(s)[0] )/2)
integrand

Out[20]:

$$\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 [21]:

zhat = Lambda(t, integrate(integrand, (s, 0, t)))
zhat

Out[21]:

$$\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 [22]:

factor(zhat(t))

Out[22]:

$$- \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 [23]:

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 [24]:

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 [27]:

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

Out[27]:

$$\left[\begin{matrix}\frac{2 \pi}{b_{1}} h_{1}\\\frac{2 \pi}{b_{1}} h_{2}\end{matrix}\right]$$

In [29]:

A0 = simplify(Matrix([
    [diff(exp21, h1), diff(exp21, h2)],
    [diff(exp22, h1), diff(exp22, h2)],
]))
A0

Out[29]:

$$\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 [32]:

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

Out[32]:

$$\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 [33]:

var('dt')

Out[33]:

$$dt$$

In [34]:

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

In [35]:

tc

Out[35]:

$$- \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 [36]:

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 [37]:

simplify(tc - tc2)

Out[37]:

$$\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 [38]:

tc = tc2

Next cell.

In [39]:

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

In [40]:

var("u1 u2 u5")

Out[40]:

$$\left ( u_{1}, \quad u_{2}, \quad u_{5}\right )$$

In [41]:

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 [43]:

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

Out[43]:

$$- \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 [44]:

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 [45]:

simplify(my_psi - his_psi)

Out[45]:

$$0$$

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

Next cell.

In [46]:

var("nu0")

Out[46]:

$$\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 [47]:

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

In [48]:

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

In [49]:

solve(Eq(s1, s2))

Out[49]:

$$\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 [50]:

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

Out[50]:

$$\left [ \right ]$$

In [51]:

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 [52]:

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

In [53]:

Out[53]:

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

Let's finish.

In [54]:

var("eta f F");

In [55]:

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

In [56]:

Out[56]:

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

Next cell.

In [57]:

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

In [59]:

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

Out[59]:

$$\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 [60]:

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 [61]:

simplify(res)

Out[61]:

$$\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 [62]:

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

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

In [63]:

g1 = Lambda((t, h1, h2), simplify(Matrix(list(xhat(t)) + [0])))

In [64]:

g1(t, h1, h2)

Out[64]:

$$\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 [65]:

g2 = Lambda((t, h1, h2), simplify(Matrix([exp21, exp22] + [zhat(t)])))

In [66]:

g2(t, h1, h2)

Out[66]:

$$\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 [67]:

simplify(d2f)

Out[67]:

$$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 [68]:

simplify(d2f.subs({
    g1: Lambda((t, h1, h2), simplify(Matrix(list(xhat(t)) + [0]))),
    g2: Lambda((t, h1, h2), simplify(Matrix([exp21, exp22] +[zhat(t)]))),
}))

Out[68]:

$$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 [69]:

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 [70]:

t = 2 * pi / b1

In [71]:

d2f

Out[71]:

$$\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 [72]:

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 [73]:

d2Exp = simplify(d2f)

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

In [74]:

d2Exp

Out[74]:

$$\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 [75]:

Psi

Out[75]:

$$\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 [76]:

Psi(d2Exp[0], d2Exp[1], d2Exp[2])

Out[76]:

$$- \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 [77]:

simplify(expand(Psi(d2Exp[0], d2Exp[1], d2Exp[2]) / (1 / (1/b1 * 2 * (h1**2 + h2**2) * pi))))

Out[77]:

$$\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.