I want to reproduce the symbolic (algebraic) computations done in §5.A of L.S.'s PhD thesis.
In [1]:
%load_ext watermark
%watermark -v -m -p sympy -g
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]:
In [5]:
solve(x + 2 + 19)
Out[5]:
In [6]:
solve(x + 2*y + 19)
Out[6]:
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]:
In [14]:
h0 = Matrix([h1, h2])
h0
Out[14]:
$\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]:
In [16]:
xhat = Lambda(t, integrate(hhat(s), (s, 0, t)))
xhat
Out[16]:
$\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]:
In [21]:
zhat = Lambda(t, integrate(integrand, (s, 0, t)))
zhat
Out[21]:
As we asked $b_1 > 0$, this won't get too complicated:
In [22]:
factor(zhat(t))
Out[22]:
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]:
In [29]:
A0 = simplify(Matrix([
[diff(exp21, h1), diff(exp21, h2)],
[diff(exp22, h1), diff(exp22, h2)],
]))
A0
Out[29]:
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]:
We need one more variable to solve an equation.
In [33]:
var('dt')
Out[33]:
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]:
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]:
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]:
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]:
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]:
Ok, we have the same result so far! Great!
Next cell.
In [46]:
var("nu0")
Out[46]:
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]:
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]:
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]:
v
Out[53]:
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]:
d
Out[56]:
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]:
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]:
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]:
In [65]:
g2 = Lambda((t, h1, h2), simplify(Matrix([exp21, exp22] + [zhat(t)])))
In [66]:
g2(t, h1, h2)
Out[66]:
Let's see if the replacement is possible:
In [67]:
simplify(d2f)
Out[67]:
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]:
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]:
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]:
In [75]:
Psi
Out[75]:
So now we can call $\Psi$ on the three components of this d2Exp :
In [76]:
Psi(d2Exp[0], d2Exp[1], d2Exp[2])
Out[76]:
In [77]:
simplify(expand(Psi(d2Exp[0], d2Exp[1], d2Exp[2]) / (1 / (1/b1 * 2 * (h1**2 + h2**2) * pi))))
Out[77]:
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.