# # Presentation # # I want to reproduce the symbolic (algebraic) computations done in §5.A of [L.S.'s PhD thesis](http://www.cmap.polytechnique.fr/~sacchelli/). # I want to only use a free and open-source software, so I'm using [Python 3](https://docs.python.org/3) with the [Sympy](https://sympy.org) module. # In: get_ipython().run_line_magic('load_ext', 'watermark') get_ipython().run_line_magic('watermark', '-v -m -p sympy -g') # 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") # In: solve(x + 2 + 19) # In: solve(x + 2*y + 19) # # 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 # In: 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$. # In: hhat = Lambda(t, (t * J).exp() * h0) hhat # In: xhat = Lambda(t, integrate(hhat(s), (s, 0, t))) xhat # $\hat{z}$ is slightly more complex: # In: integrand = simplify(b1 * (hhat(s) * xhat(s) - hhat(s) * xhat(s) )/2) integrand # In: zhat = Lambda(t, integrate(integrand, (s, 0, t))) zhat # As we asked $b_1 > 0$, this won't get too complicated: # In: factor(zhat(t)) # 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 # In: A0 = simplify(Matrix([ [diff(exp21, h1), diff(exp21, h2)], [diff(exp22, h1), diff(exp22, h2)], ])) A0 # 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 # We need one more variable to solve an equation. # In: var('dt') # In: tc = factor(solve(Equality((jacobian + Matrix([[dt,0,0], [0,dt,0], [0,0,0]])).det(), 0), dt)) # In: tc # 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) # 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") # 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 # 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) # Ok, we have the same result so far! Great! # Next cell. # In: 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. # 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)) # 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) # 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 # 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 # 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 # 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) # 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) # In: g2 = Lambda((t, h1, h2), simplify(Matrix([exp21, exp22] + [zhat(t)]))) # In: g2(t, h1, h2) # Let's see if the replacement is possible: # In: simplify(d2f) # 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)]))), })) # 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 # 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 # In: Psi # So now we can call $\Psi$ on the three components of this d2Exp : # In: Psi(d2Exp, d2Exp, d2Exp) # In: simplify(expand(Psi(d2Exp, d2Exp, d2Exp) / (1 / (1/b1 * 2 * (h1**2 + h2**2) * pi)))) # 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](https://github.com/Naereen/notebooks) for other notebooks I wrote.