I want to reproduce the symbolic (algebraic) computations done in §5.A of L.S.'s PhD thesis.
%load_ext watermark
%watermark -v -m -p sympy -g
from sympy import *
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:
var("x y z")
solve(x + 2 + 19)
solve(x + 2*y + 19)
We start by defining $b_1 > 0$ and the other variables.
b1 = symbols("b1", positive=True)
We need a lot of variables.
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.
J = Matrix([[0, -b1], [b1, 0]])
J
h0 = Matrix([h1, h2])
h0
$\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$.
hhat = Lambda(t, (t * J).exp() * h0)
hhat
xhat = Lambda(t, integrate(hhat(s), (s, 0, t)))
xhat
$\hat{z}$ is slightly more complex:
integrand = simplify(b1 * (hhat(s)[0] * xhat(s)[1] - hhat(s)[1] * xhat(s)[0] )/2)
integrand
zhat = Lambda(t, integrate(integrand, (s, 0, t)))
zhat
As we asked $b_1 > 0$, this won't get too complicated:
factor(zhat(t))
Two more expressions:
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))
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:
C10 = 2 * (pi / b1) * Matrix([h1, h2])
C10
A0 = simplify(Matrix([
[diff(exp21, h1), diff(exp21, h2)],
[diff(exp22, h1), diff(exp22, h2)],
]))
A0
rest1
and rest2
already vanished.
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
We need one more variable to solve an equation.
var('dt')
tc = factor(solve(Equality((jacobian + Matrix([[dt,0,0], [0,dt,0], [0,0,0]])).det(), 0), dt)[0])
tc
I can compare by copying the result from the document:
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!
simplify(tc - tc2)
Let's use the one from Ludovic's notebook.
tc = tc2
Next cell.
A12 = simplify(A0 + Matrix([[tc, 0], [0, tc]]))
var("u1 u2 u5")
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)
))
my_psi = factor(simplify(Psi(u1, u2, u5)))
my_psi
Let's compare with the value from Ludovic's notebook:
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))
)
simplify(my_psi - his_psi)
Ok, we have the same result so far! Great!
Next cell.
var("nu0")
Here again, Sympy fails to solve the equation. It can solve on both lines separately, but cannot find a solution that satisfies both lines.
s1 = solve(Eq((Matrix(A12.dot(Matrix([-h2, h1]))) + nu0 * C10)[0]), nu0)[0]
s2 = solve(Eq((Matrix(A12.dot(Matrix([-h2, h1]))) + nu0 * C10)[1]), nu0)[0]
solve(Eq(s1, s2))
It is not possible to impose these constraints. And Sympy fails to solve both equation simultaneously:
solve(Eq(Matrix(A12.dot(Matrix([-h2, h1]))) + nu0 * C10, Matrix([0,0])), nu0)
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
).
v = Matrix([nu, -h2, h1])
v
Let's finish.
var("eta f F");
d = [
lambda F: -eta**2 * diff(F, eta) + eta * t * diff(F, t),
lambda F: diff(F, h1),
lambda F: diff(F, h2)
]
d
Next cell.
var("i1 i2 g1 g2 dt1 dt2");
g = eta * g1(t + eta * dt1, h1, h2) + eta**2 * g2(t, h1, h2)
g
We have a sum to compute:
res = 0
for i1 in range(0, 3):
for i2 in range(0, 3):
res += v[i1] * v[i2] * d[i1](d[i2](g))
simplify(res)
Then a double derivative
d2f = simplify((diff(res, eta, eta) / 2).subs(eta, 0))
Now, we should replace $g_1$ and $g_2$ with the following expression:
g1 = Lambda((t, h1, h2), simplify(Matrix(list(xhat(t)) + [0])))
g1(t, h1, h2)
g2 = Lambda((t, h1, h2), simplify(Matrix([exp21, exp22] + [zhat(t)])))
g2(t, h1, h2)
Let's see if the replacement is possible:
simplify(d2f)
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)]))),
}))
OK. This failed. But we can copy and paste this and the replacement of $g_1$ and $g_2$ will be effective:
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$.
t = 2 * pi / b1
d2f
The replacement does not work apparently, so let's do it manually:
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)]])
d2Exp = simplify(d2f)
It still works well. And we removed the complex exponential, this result is purely real now!
d2Exp
Psi
So now we can call $\Psi$ on the three components of this d2Exp
:
Psi(d2Exp[0], d2Exp[1], d2Exp[2])
simplify(expand(Psi(d2Exp[0], d2Exp[1], d2Exp[2]) / (1 / (1/b1 * 2 * (h1**2 + h2**2) * pi))))
We don't have the value for nu
!
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.