Algebraic_integration_from_Ludovic_Sacchelli_s_PhD_thesis.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "toc": "true"
   },
   "source": [
    "# Table of Contents\n",
    " <p><div class=\"lev1 toc-item\"><a href=\"#Presentation\" data-toc-modified-id=\"Presentation-1\"><span class=\"toc-item-num\">1&nbsp;&nbsp;</span>Presentation</a></div><div class=\"lev1 toc-item\"><a href=\"#Computations\" data-toc-modified-id=\"Computations-2\"><span class=\"toc-item-num\">2&nbsp;&nbsp;</span>Computations</a></div><div class=\"lev2 toc-item\"><a href=\"#Conclusion\" data-toc-modified-id=\"Conclusion-21\"><span class=\"toc-item-num\">2.1&nbsp;&nbsp;</span>Conclusion</a></div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Presentation\n",
    "\n",
    "I want to reproduce the symbolic (algebraic) computations done in §5.A of [L.S.'s PhD thesis](http://www.cmap.polytechnique.fr/~sacchelli/)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "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."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPython 3.6.5\n",
      "IPython 6.4.0\n",
      "\n",
      "sympy 1.1.1\n",
      "\n",
      "compiler   : GCC 7.3.0\n",
      "system     : Linux\n",
      "release    : 4.15.0-23-generic\n",
      "machine    : x86_64\n",
      "processor  : x86_64\n",
      "CPU cores  : 4\n",
      "interpreter: 64bit\n",
      "Git hash   : 3f62e1141f480e65f685cabe81f7a79b17bcd06d\n"
     ]
    }
   ],
   "source": [
    "%load_ext watermark\n",
    "%watermark -v -m -p sympy -g"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sympy import *"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "init_printing(use_latex='mathjax')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Just a small introduction to SymPy: it works by using Python expressions on formal variables.\n",
    "First, we define variables, and then we can solve linear equations for example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$\\left ( x, \\quad y, \\quad z\\right )$$"
      ],
      "text/plain": [
       "(x, y, z)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "var(\"x y z\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$\\left [ -21\\right ]$$"
      ],
      "text/plain": [
       "[-21]"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "solve(x + 2 + 19)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$\\left [ \\left \\{ x : - 2 y - 19\\right \\}\\right ]$$"
      ],
      "text/plain": [
       "[{x: -2⋅y - 19}]"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "solve(x + 2*y + 19)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Computations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We start by defining $b_1 > 0$ and the other variables."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "b1 = symbols(\"b1\", positive=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We need a lot of variables."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "var(\"h1 h2 h3 h4 k131 k132 k141 k142 k231 k232 k241 k242 rest1 rest2 t s a b\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then we can start to follow Ludovic's notebook and define the first expressions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$\\left[\\begin{matrix}0 & - b_{1}\\\\b_{1} & 0\\end{matrix}\\right]$$"
      ],
      "text/plain": [
       "⎡0   -b₁⎤\n",
       "⎢       ⎥\n",
       "⎣b₁   0 ⎦"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "J = Matrix([[0, -b1], [b1, 0]])\n",
    "J"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$\\left[\\begin{matrix}h_{1}\\\\h_{2}\\end{matrix}\\right]$$"
      ],
      "text/plain": [
       "⎡h₁⎤\n",
       "⎢  ⎥\n",
       "⎣h₂⎦"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "h0 = Matrix([h1, h2])\n",
    "h0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\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$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$\\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)$$"
      ],
      "text/plain": [
       "    ⎡    ⎛ ⅈ⋅b₁⋅t    -ⅈ⋅b₁⋅t⎞      ⎛   ⅈ⋅b₁⋅t      -ⅈ⋅b₁⋅t⎞ ⎤\n",
       "    ⎢    ⎜ℯ         ℯ       ⎟      ⎜ⅈ⋅ℯ         ⅈ⋅ℯ       ⎟ ⎥\n",
       "    ⎢ h₁⋅⎜─────── + ────────⎟ + h₂⋅⎜───────── - ──────────⎟ ⎥\n",
       "    ⎢    ⎝   2         2    ⎠      ⎝    2           2     ⎠ ⎥\n",
       "t ↦ ⎢                                                       ⎥\n",
       "    ⎢   ⎛     ⅈ⋅b₁⋅t      -ⅈ⋅b₁⋅t⎞      ⎛ ⅈ⋅b₁⋅t    -ⅈ⋅b₁⋅t⎞⎥\n",
       "    ⎢   ⎜  ⅈ⋅ℯ         ⅈ⋅ℯ       ⎟      ⎜ℯ         ℯ       ⎟⎥\n",
       "    ⎢h₁⋅⎜- ───────── + ──────────⎟ + h₂⋅⎜─────── + ────────⎟⎥\n",
       "    ⎣   ⎝      2           2     ⎠      ⎝   2         2    ⎠⎦"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "hhat = Lambda(t, (t * J).exp() * h0)\n",
    "hhat"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$\\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)$$"
      ],
      "text/plain": [
       "    ⎡                               ⅈ⋅b₁⋅t                          -ⅈ⋅b₁⋅t⎤\n",
       "    ⎢  h₂   (-2⋅ⅈ⋅b₁⋅h₁ + 2⋅b₁⋅h₂)⋅ℯ       + (2⋅ⅈ⋅b₁⋅h₁ + 2⋅b₁⋅h₂)⋅ℯ       ⎥\n",
       "    ⎢- ── + ───────────────────────────────────────────────────────────────⎥\n",
       "    ⎢  b₁                                    2                             ⎥\n",
       "    ⎢                                    4⋅b₁                              ⎥\n",
       "t ↦ ⎢                                                                      ⎥\n",
       "    ⎢                             ⅈ⋅b₁⋅t                           -ⅈ⋅b₁⋅t ⎥\n",
       "    ⎢h₁   (-2⋅b₁⋅h₁ - 2⋅ⅈ⋅b₁⋅h₂)⋅ℯ       + (-2⋅b₁⋅h₁ + 2⋅ⅈ⋅b₁⋅h₂)⋅ℯ        ⎥\n",
       "    ⎢── + ──────────────────────────────────────────────────────────────── ⎥\n",
       "    ⎢b₁                                    2                               ⎥\n",
       "    ⎣                                  4⋅b₁                                ⎦"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "xhat = Lambda(t, integrate(hhat(s), (s, 0, t)))\n",
    "xhat"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\hat{z}$ is slightly more complex:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$\\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}$$"
      ],
      "text/plain": [
       "  2  ⅈ⋅b₁⋅s     2     2  -ⅈ⋅b₁⋅s     2  ⅈ⋅b₁⋅s     2     2  -ⅈ⋅b₁⋅s\n",
       "h₁ ⋅ℯ         h₁    h₁ ⋅ℯ          h₂ ⋅ℯ         h₂    h₂ ⋅ℯ       \n",
       "─────────── - ─── + ──────────── + ─────────── - ─── + ────────────\n",
       "     4         2         4              4         2         4      "
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "integrand = simplify(b1 * (hhat(s)[0] * xhat(s)[1] - hhat(s)[1] * xhat(s)[0] )/2)\n",
    "integrand"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$\\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)$$"
      ],
      "text/plain": [
       "      ⎛    2     2⎞   ⎛           2            2⎞  ⅈ⋅b₁⋅t   ⎛         2       \n",
       "      ⎜  h₁    h₂ ⎟   ⎝- 4⋅ⅈ⋅b₁⋅h₁  - 4⋅ⅈ⋅b₁⋅h₂ ⎠⋅ℯ       + ⎝4⋅ⅈ⋅b₁⋅h₁  + 4⋅ⅈ⋅\n",
       "t ↦ t⋅⎜- ─── - ───⎟ + ────────────────────────────────────────────────────────\n",
       "      ⎝   2     2 ⎠                                         2                 \n",
       "                                                       16⋅b₁                  \n",
       "\n",
       "     2⎞  -ⅈ⋅b₁⋅t\n",
       "b₁⋅h₂ ⎠⋅ℯ       \n",
       "────────────────\n",
       "                \n",
       "                "
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "zhat = Lambda(t, integrate(integrand, (s, 0, t)))\n",
    "zhat"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we asked $b_1 > 0$, this won't get too complicated:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$- \\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)$$"
      ],
      "text/plain": [
       " ⎛  2     2⎞ ⎛        ⅈ⋅b₁⋅t      2⋅ⅈ⋅b₁⋅t    ⎞  -ⅈ⋅b₁⋅t \n",
       "-⎝h₁  + h₂ ⎠⋅⎝2⋅b₁⋅t⋅ℯ       + ⅈ⋅ℯ         - ⅈ⎠⋅ℯ        \n",
       "─────────────────────────────────────────────────────────\n",
       "                           4⋅b₁                          "
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "factor(zhat(t))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Two more expressions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "exp21 = (3 * a * h1**2 + a * h2**2 + 2 * b * h1 * h2 + k131 * h1 * h3 + k141 * h1 * h4\n",
    "         + k231 * h2 * h3 + k241 * h2 * h4 + rest1(h3, h4))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "exp22 = (b * h1**2 + 3 * b * h2**2 + 2 * a * h1 * h2 + k132 * h1 * h3 + k142 * h1 * h4\n",
    "         + k232 * h2 * h3 + k242 * h2 * h4 + rest2(h3, h4))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "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$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So far so good. Next cell:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$\\left[\\begin{matrix}\\frac{2 \\pi}{b_{1}} h_{1}\\\\\\frac{2 \\pi}{b_{1}} h_{2}\\end{matrix}\\right]$$"
      ],
      "text/plain": [
       "⎡2⋅π⋅h₁⎤\n",
       "⎢──────⎥\n",
       "⎢  b₁  ⎥\n",
       "⎢      ⎥\n",
       "⎢2⋅π⋅h₂⎥\n",
       "⎢──────⎥\n",
       "⎣  b₁  ⎦"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "C10 = 2 * (pi / b1) * Matrix([h1, h2])\n",
    "C10"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$\\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]$$"
      ],
      "text/plain": [
       "⎡6⋅a⋅h₁ + 2⋅b⋅h₂ + h₃⋅k₁₃₁ + h₄⋅k₁₄₁  2⋅a⋅h₂ + 2⋅b⋅h₁ + h₃⋅k₂₃₁ + h₄⋅k₂₄₁⎤\n",
       "⎢                                                                        ⎥\n",
       "⎣2⋅a⋅h₂ + 2⋅b⋅h₁ + h₃⋅k₁₃₂ + h₄⋅k₁₄₂  2⋅a⋅h₁ + 6⋅b⋅h₂ + h₃⋅k₂₃₂ + h₄⋅k₂₄₂⎦"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A0 = simplify(Matrix([\n",
    "    [diff(exp21, h1), diff(exp21, h2)],\n",
    "    [diff(exp22, h1), diff(exp22, h2)],\n",
    "]))\n",
    "A0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`rest1` and `rest2` already vanished."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$\\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]$$"
      ],
      "text/plain": [
       "⎡                                                                          2⋅π\n",
       "⎢6⋅a⋅h₁ + 2⋅b⋅h₂ + h₃⋅k₁₃₁ + h₄⋅k₁₄₁  2⋅a⋅h₂ + 2⋅b⋅h₁ + h₃⋅k₂₃₁ + h₄⋅k₂₄₁  ───\n",
       "⎢                                                                            b\n",
       "⎢                                                                             \n",
       "⎢                                                                          2⋅π\n",
       "⎢2⋅a⋅h₂ + 2⋅b⋅h₁ + h₃⋅k₁₃₂ + h₄⋅k₁₄₂  2⋅a⋅h₁ + 6⋅b⋅h₂ + h₃⋅k₂₃₂ + h₄⋅k₂₄₂  ───\n",
       "⎢                                                                            b\n",
       "⎢                                                                             \n",
       "⎢             -2⋅π⋅h₁                              -2⋅π⋅h₂                    \n",
       "⎢             ────────                             ────────                  0\n",
       "⎣                b₁                                   b₁                      \n",
       "\n",
       "⋅h₁⎤\n",
       "───⎥\n",
       "₁  ⎥\n",
       "   ⎥\n",
       "⋅h₂⎥\n",
       "───⎥\n",
       "₁  ⎥\n",
       "   ⎥\n",
       "   ⎥\n",
       "   ⎥\n",
       "   ⎦"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "jacobian = simplify(Matrix([\n",
    "    [diff(exp21, h1), diff(exp21, h2), C10[0]],\n",
    "    [diff(exp22, h1), diff(exp22, h2), C10[1]],\n",
    "    [diff(zhat(2 * pi / b1), h1), diff(zhat(2 * pi / b1), h2), 0],\n",
    "]))\n",
    "jacobian"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We need one more variable to solve an equation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$dt$$"
      ],
      "text/plain": [
       "dt"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "var('dt')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "tc = factor(solve(Equality((jacobian + Matrix([[dt,0,0], [0,dt,0], [0,0,0]])).det(), 0), dt)[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$- \\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)$$"
      ],
      "text/plain": [
       " ⎛      3            2         2            3     2             2             \n",
       "-⎝2⋅a⋅h₁  + 2⋅a⋅h₁⋅h₂  + 2⋅b⋅h₁ ⋅h₂ + 2⋅b⋅h₂  + h₁ ⋅h₃⋅k₂₃₂ + h₁ ⋅h₄⋅k₂₄₂ - h₁\n",
       "──────────────────────────────────────────────────────────────────────────────\n",
       "                                                                              \n",
       "                                                                              \n",
       "\n",
       "                                                                2             \n",
       "⋅h₂⋅h₃⋅k₁₃₂ - h₁⋅h₂⋅h₃⋅k₂₃₁ - h₁⋅h₂⋅h₄⋅k₁₄₂ - h₁⋅h₂⋅h₄⋅k₂₄₁ + h₂ ⋅h₃⋅k₁₃₁ + h₂\n",
       "──────────────────────────────────────────────────────────────────────────────\n",
       "   2     2                                                                    \n",
       " h₁  + h₂                                                                     \n",
       "\n",
       "2        ⎞ \n",
       " ⋅h₄⋅k₁₄₁⎠ \n",
       "───────────\n",
       "           \n",
       "           "
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tc"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "I can compare by copying the result from the document:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [],
   "source": [
    "tc2 = (-1 / (h1**2 + h2**2)) * (\n",
    "      2 * a * h1**3\n",
    "    + 6 * a * h1**2 * h2\n",
    "    - 4 * b * h1**2 * h2\n",
    "    + 2 * a * h1 * h2**2\n",
    "    + 2 * b * h2**3\n",
    "    + h2**2 * h3 * k131\n",
    "    - h1 * h2 * h3 * k132\n",
    "    + h2**2 * h4 * k141\n",
    "    - h1 * h2 * h4 * k142\n",
    "    - h1 * h2 * h3 * k231\n",
    "    + h1**2 * h3 * k232\n",
    "    - h1 * h2 * h4 * k241\n",
    "    + h1**2 * h4 * k242 )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Drat, they are not equal! We might have a mistake somewhere, even though I just CAN'T find it!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$\\frac{6 h_{1}^{2} h_{2} \\left(a - b\\right)}{h_{1}^{2} + h_{2}^{2}}$$"
      ],
      "text/plain": [
       "    2           \n",
       "6⋅h₁ ⋅h₂⋅(a - b)\n",
       "────────────────\n",
       "     2     2    \n",
       "   h₁  + h₂     "
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "simplify(tc - tc2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's use the one from Ludovic's notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [],
   "source": [
    "tc = tc2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next cell."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [],
   "source": [
    "A12 = simplify(A0 + Matrix([[tc, 0], [0, tc]]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$\\left ( u_{1}, \\quad u_{2}, \\quad u_{5}\\right )$$"
      ],
      "text/plain": [
       "(u₁, u₂, u₅)"
      ]
     },
     "execution_count": 40,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "var(\"u1 u2 u5\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [],
   "source": [
    "Psi = Lambda((u1, u2, u5), simplify(\n",
    "    u5 * Matrix(A12.dot(Matrix([h1, h2]))).dot(Matrix([h2, -h1])) + 2 * pi / b1 * ( h1**2 + h2**2) * (h1 * u2 - h2 * u1)\n",
    "))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$- \\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)$$"
      ],
      "text/plain": [
       " ⎛           2                  3               3                  2          \n",
       "-⎝- 2⋅a⋅b₁⋅h₁ ⋅h₂⋅u₅ - 2⋅a⋅b₁⋅h₂ ⋅u₅ + 2⋅b⋅b₁⋅h₁ ⋅u₅ + 2⋅b⋅b₁⋅h₁⋅h₂ ⋅u₅ + b₁⋅h\n",
       "──────────────────────────────────────────────────────────────────────────────\n",
       "                                                                              \n",
       "\n",
       " 2                   2                                                        \n",
       "₁ ⋅h₃⋅k₁₃₂⋅u₅ + b₁⋅h₁ ⋅h₄⋅k₁₄₂⋅u₅ - b₁⋅h₁⋅h₂⋅h₃⋅k₁₃₁⋅u₅ + b₁⋅h₁⋅h₂⋅h₃⋅k₂₃₂⋅u₅ \n",
       "──────────────────────────────────────────────────────────────────────────────\n",
       "                                                                       b₁     \n",
       "\n",
       "                                                   2                   2      \n",
       "- b₁⋅h₁⋅h₂⋅h₄⋅k₁₄₁⋅u₅ + b₁⋅h₁⋅h₂⋅h₄⋅k₂₄₂⋅u₅ - b₁⋅h₂ ⋅h₃⋅k₂₃₁⋅u₅ - b₁⋅h₂ ⋅h₄⋅k₂\n",
       "──────────────────────────────────────────────────────────────────────────────\n",
       "                                                                              \n",
       "\n",
       "              3            2                  2            3   ⎞ \n",
       "₄₁⋅u₅ - 2⋅π⋅h₁ ⋅u₂ + 2⋅π⋅h₁ ⋅h₂⋅u₁ - 2⋅π⋅h₁⋅h₂ ⋅u₂ + 2⋅π⋅h₂ ⋅u₁⎠ \n",
       "─────────────────────────────────────────────────────────────────\n",
       "                                                                 "
      ]
     },
     "execution_count": 43,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "my_psi = factor(simplify(Psi(u1, u2, u5)))\n",
    "my_psi"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's compare with the value from Ludovic's notebook:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [],
   "source": [
    "his_psi = (1 / b1) * (\n",
    "    2 * h1**3 * (pi * u2 - b * b1 * u5 ) -\n",
    "    h1**2 * (2 * h2 * pi * u1 - 2 * a * b1 * h2 * u5 + b1 * h3 * k132 * u5 +\n",
    "             b1 * h4 * k142 * u5 ) +\n",
    "    h2**2 * ( -2 * h2 * pi * u1 + 2 * a * b1 * h2 * u5 + b1 * h3 * k231 * u5 +\n",
    "            b1 * h4 * k241 * u5 ) +\n",
    "    h1 * h2 * (b1 * (h3 * k131 + h4 * k141 - h3 * k232 - h4 * k242) * u5 +\n",
    "    2 * h2 * (pi * u2 - b * b1 * u5))\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$0$$"
      ],
      "text/plain": [
       "0"
      ]
     },
     "execution_count": 45,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "simplify(my_psi - his_psi)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ok, we have the same result so far! Great!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next cell."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$\\nu_{0}$$"
      ],
      "text/plain": [
       "ν₀"
      ]
     },
     "execution_count": 46,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "var(\"nu0\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here again, Sympy fails to solve the equation. It can solve on both lines separately, but cannot find a solution that satisfies both lines."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [],
   "source": [
    "s1 = solve(Eq((Matrix(A12.dot(Matrix([-h2, h1]))) + nu0 * C10)[0]), nu0)[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [],
   "source": [
    "s2 = solve(Eq((Matrix(A12.dot(Matrix([-h2, h1]))) + nu0 * C10)[1]), nu0)[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$\\left [ \\left \\{ a : b\\right \\}, \\quad \\left \\{ h_{1} : 0\\right \\}\\right ]$$"
      ],
      "text/plain": [
       "[{a: b}, {h₁: 0}]"
      ]
     },
     "execution_count": 49,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "solve(Eq(s1, s2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is not possible to impose these constraints. And Sympy fails to solve both equation simultaneously:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$\\left [ \\right ]$$"
      ],
      "text/plain": [
       "[]"
      ]
     },
     "execution_count": 50,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "solve(Eq(Matrix(A12.dot(Matrix([-h2, h1]))) + nu0 * C10, Matrix([0,0])), nu0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [],
   "source": [
    "nu = simplify(solve(Eq(Matrix(A12.dot(Matrix([-h2, h1]))) + nu0 * C10, Matrix([0,0])), nu0))\n",
    "if nu == []:\n",
    "    nu = var(\"nu\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can continue by using a formal variable for $\\nu$ (`nu`)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [],
   "source": [
    "v = Matrix([nu, -h2, h1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$\\left[\\begin{matrix}\\nu\\\\- h_{2}\\\\h_{1}\\end{matrix}\\right]$$"
      ],
      "text/plain": [
       "⎡ ν ⎤\n",
       "⎢   ⎥\n",
       "⎢-h₂⎥\n",
       "⎢   ⎥\n",
       "⎣h₁ ⎦"
      ]
     },
     "execution_count": 53,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "v"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's finish."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [],
   "source": [
    "var(\"eta f F\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [],
   "source": [
    "d = [\n",
    "    lambda F: -eta**2 * diff(F, eta) + eta * t * diff(F, t),\n",
    "    lambda F: diff(F, h1),\n",
    "    lambda F: diff(F, h2)\n",
    "]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<function __main__.<lambda>(F)>,\n",
       " <function __main__.<lambda>(F)>,\n",
       " <function __main__.<lambda>(F)>]"
      ]
     },
     "execution_count": 56,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "d"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next cell."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [],
   "source": [
    "var(\"i1 i2 g1 g2 dt1 dt2\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$\\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 )}$$"
      ],
      "text/plain": [
       " 2                                        \n",
       "η ⋅g₂(t, h₁, h₂) + η⋅g₁(dt₁⋅η + t, h₁, h₂)"
      ]
     },
     "execution_count": 59,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "g = eta * g1(t + eta * dt1, h1, h2) + eta**2 * g2(t, h1, h2)\n",
    "g"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have a sum to compute:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {},
   "outputs": [],
   "source": [
    "res = 0\n",
    "for i1 in range(0, 3):\n",
    "    for i2 in range(0, 3):\n",
    "        res += v[i1] * v[i2] * d[i1](d[i2](g))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$\\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)$$"
      ],
      "text/plain": [
       "  ⎛      ⎛                                               ⎛       ⎛  2         \n",
       "  ⎜ 2  2 ⎜        ⎛ ∂                 ⎞│                 ⎜   2   ⎜ ∂          \n",
       "η⋅⎜η ⋅ν ⋅⎜2⋅dt₁⋅η⋅⎜───(g₁(ξ₁, h₁, h₂))⎟│             + η⋅⎜dt₁ ⋅η⋅⎜────(g₁(ξ₁, \n",
       "  ⎜      ⎜        ⎝∂ξ₁                ⎠│ξ₁=dt₁⋅η + t     ⎜       ⎜   2        \n",
       "  ⎝      ⎝                                               ⎝       ⎝∂ξ₁         \n",
       "\n",
       "        ⎞│                                                                    \n",
       "        ⎟│                     ⎛ ∂                 ⎞│                         \n",
       "h₁, h₂))⎟│             + 2⋅dt₁⋅⎜───(g₁(ξ₁, h₁, h₂))⎟│             + 2⋅g₂(t, h₁\n",
       "        ⎟│                     ⎝∂ξ₁                ⎠│ξ₁=dt₁⋅η + t             \n",
       "        ⎠│ξ₁=dt₁⋅η + t                                                        \n",
       "\n",
       "     ⎞                                                                        \n",
       "     ⎟                         ⎛  ∂                   ⎛ ∂                 ⎞│  \n",
       ", h₂)⎟ + 4⋅η⋅g₂(t, h₁, h₂) - t⋅⎜η⋅──(g₂(t, h₁, h₂)) + ⎜───(g₁(ξ₁, h₁, h₂))⎟│  \n",
       "     ⎟                         ⎝  ∂t                  ⎝∂ξ₁                ⎠│ξ₁\n",
       "     ⎠                                                                        \n",
       "\n",
       "                ⎛        ⎛  2                 ⎞│                              \n",
       "          ⎞     ⎜        ⎜ ∂                  ⎟│                 ∂            \n",
       "          ⎟ + t⋅⎜- dt₁⋅η⋅⎜────(g₁(ξ₁, h₁, h₂))⎟│             - η⋅──(g₂(t, h₁, \n",
       "=dt₁⋅η + t⎠     ⎜        ⎜   2                ⎟│                 ∂t           \n",
       "                ⎝        ⎝∂ξ₁                 ⎠│ξ₁=dt₁⋅η + t                  \n",
       "\n",
       "         ⎛    2                  ⎛  2                 ⎞│            ⎞⎞     ⎛  \n",
       "         ⎜   ∂                   ⎜ ∂                  ⎟│            ⎟⎟     ⎜  \n",
       "h₂)) + t⋅⎜η⋅───(g₂(t, h₁, h₂)) + ⎜────(g₁(ξ₁, h₁, h₂))⎟│            ⎟⎟ - t⋅⎜dt\n",
       "         ⎜    2                  ⎜   2                ⎟│            ⎟⎟     ⎜  \n",
       "         ⎝  ∂t                   ⎝∂ξ₁                 ⎠│ξ₁=dt₁⋅η + t⎠⎠     ⎝  \n",
       "\n",
       "    ⎛  2                 ⎞│                                                   \n",
       "    ⎜ ∂                  ⎟│                   ∂                   ⎛ ∂         \n",
       "₁⋅η⋅⎜────(g₁(ξ₁, h₁, h₂))⎟│             + 2⋅η⋅──(g₂(t, h₁, h₂)) + ⎜───(g₁(ξ₁, \n",
       "    ⎜   2                ⎟│                   ∂t                  ⎝∂ξ₁        \n",
       "    ⎝∂ξ₁                 ⎠│ξ₁=dt₁⋅η + t                                       \n",
       "\n",
       "                      ⎞                          ⎞            ⎛      ⎛    2   \n",
       "        ⎞│            ⎟                          ⎟            ⎜      ⎜   ∂    \n",
       "h₁, h₂))⎟│            ⎟ + 2⋅g₁(dt₁⋅η + t, h₁, h₂)⎟ - 2⋅η⋅h₁⋅ν⋅⎜dt₁⋅η⋅⎜───────(\n",
       "        ⎠│ξ₁=dt₁⋅η + t⎟                          ⎟            ⎝      ⎝∂h₂ ∂ξ₁ \n",
       "                      ⎠                          ⎠                            \n",
       "\n",
       "               ⎞│                                          ⎛     2            \n",
       "               ⎟│                    ∂                     ⎜    ∂             \n",
       "g₁(ξ₁, h₁, h₂))⎟│             + 2⋅η⋅───(g₂(t, h₁, h₂)) - t⋅⎜η⋅──────(g₂(t, h₁,\n",
       "               ⎠│ξ₁=dt₁⋅η + t       ∂h₂                    ⎝  ∂t ∂h₂          \n",
       "                                                                              \n",
       "\n",
       "        ⎛    2                  ⎞│            ⎞                             ⎞ \n",
       "        ⎜   ∂                   ⎟│            ⎟    ∂                        ⎟ \n",
       " h₂)) + ⎜───────(g₁(ξ₁, h₁, h₂))⎟│            ⎟ + ───(g₁(dt₁⋅η + t, h₁, h₂))⎟ \n",
       "        ⎝∂h₂ ∂ξ₁                ⎠│ξ₁=dt₁⋅η + t⎠   ∂h₂                       ⎠ \n",
       "                                                                              \n",
       "\n",
       "           ⎛      ⎛    2                  ⎞│                                  \n",
       "           ⎜      ⎜   ∂                   ⎟│                    ∂             \n",
       "+ 2⋅η⋅h₂⋅ν⋅⎜dt₁⋅η⋅⎜───────(g₁(ξ₁, h₁, h₂))⎟│             + 2⋅η⋅───(g₂(t, h₁, h\n",
       "           ⎝      ⎝∂h₁ ∂ξ₁                ⎠│ξ₁=dt₁⋅η + t       ∂h₁            \n",
       "                                                                              \n",
       "\n",
       "        ⎛     2                    ⎛    2                  ⎞│            ⎞    \n",
       "        ⎜    ∂                     ⎜   ∂                   ⎟│            ⎟    \n",
       "₂)) - t⋅⎜η⋅──────(g₂(t, h₁, h₂)) + ⎜───────(g₁(ξ₁, h₁, h₂))⎟│            ⎟ + ─\n",
       "        ⎝  ∂t ∂h₁                  ⎝∂h₁ ∂ξ₁                ⎠│ξ₁=dt₁⋅η + t⎠   ∂\n",
       "                                                                              \n",
       "\n",
       "                         ⎞       ⎛    2                     2                 \n",
       "∂                        ⎟     2 ⎜   ∂                     ∂                  \n",
       "──(g₁(dt₁⋅η + t, h₁, h₂))⎟ + h₁ ⋅⎜η⋅────(g₂(t, h₁, h₂)) + ────(g₁(dt₁⋅η + t, h\n",
       "h₁                       ⎠       ⎜     2                     2                \n",
       "                                 ⎝  ∂h₂                   ∂h₂                 \n",
       "\n",
       "       ⎞           ⎛      2                        2                         ⎞\n",
       "       ⎟           ⎜     ∂                        ∂                          ⎟\n",
       "₁, h₂))⎟ - 2⋅h₁⋅h₂⋅⎜η⋅───────(g₂(t, h₁, h₂)) + ───────(g₁(dt₁⋅η + t, h₁, h₂))⎟\n",
       "       ⎟           ⎝  ∂h₂ ∂h₁                  ∂h₂ ∂h₁                       ⎠\n",
       "       ⎠                                                                      \n",
       "\n",
       "       ⎛    2                     2                        ⎞⎞\n",
       "     2 ⎜   ∂                     ∂                         ⎟⎟\n",
       " + h₂ ⋅⎜η⋅────(g₂(t, h₁, h₂)) + ────(g₁(dt₁⋅η + t, h₁, h₂))⎟⎟\n",
       "       ⎜     2                     2                       ⎟⎟\n",
       "       ⎝  ∂h₁                   ∂h₁                        ⎠⎠"
      ]
     },
     "execution_count": 61,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "simplify(res)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then a double derivative"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [],
   "source": [
    "d2f = simplify((diff(res, eta, eta) / 2).subs(eta, 0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we should replace $g_1$ and $g_2$ with the following expression:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {},
   "outputs": [],
   "source": [
    "g1 = Lambda((t, h1, h2), simplify(Matrix(list(xhat(t)) + [0])))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$\\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]$$"
      ],
      "text/plain": [
       "⎡⎛             ⅈ⋅b₁⋅t                      2⋅ⅈ⋅b₁⋅t⎞  -ⅈ⋅b₁⋅t⎤\n",
       "⎢⎝ⅈ⋅h₁ - 2⋅h₂⋅ℯ       + h₂ + (-ⅈ⋅h₁ + h₂)⋅ℯ        ⎠⋅ℯ       ⎥\n",
       "⎢────────────────────────────────────────────────────────────⎥\n",
       "⎢                            2⋅b₁                            ⎥\n",
       "⎢                                                            ⎥\n",
       "⎢⎛      ⅈ⋅b₁⋅t                            2⋅ⅈ⋅b₁⋅t⎞  -ⅈ⋅b₁⋅t ⎥\n",
       "⎢⎝2⋅h₁⋅ℯ       - h₁ + ⅈ⋅h₂ - (h₁ + ⅈ⋅h₂)⋅ℯ        ⎠⋅ℯ        ⎥\n",
       "⎢─────────────────────────────────────────────────────────── ⎥\n",
       "⎢                            2⋅b₁                            ⎥\n",
       "⎢                                                            ⎥\n",
       "⎣                             0                              ⎦"
      ]
     },
     "execution_count": 64,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "g1(t, h1, h2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {},
   "outputs": [],
   "source": [
    "g2 = Lambda((t, h1, h2), simplify(Matrix([exp21, exp22] + [zhat(t)])))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$\\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]$$"
      ],
      "text/plain": [
       "⎡      2       2                                                              \n",
       "⎢3⋅a⋅h₁  + a⋅h₂  + 2⋅b⋅h₁⋅h₂ + h₁⋅h₃⋅k₁₃₁ + h₁⋅h₄⋅k₁₄₁ + h₂⋅h₃⋅k₂₃₁ + h₂⋅h₄⋅k₂\n",
       "⎢                                                                             \n",
       "⎢                2         2                                                  \n",
       "⎢2⋅a⋅h₁⋅h₂ + b⋅h₁  + 3⋅b⋅h₂  + h₁⋅h₃⋅k₁₃₂ + h₁⋅h₄⋅k₁₄₂ + h₂⋅h₃⋅k₂₃₂ + h₂⋅h₄⋅k₂\n",
       "⎢                                                                             \n",
       "⎢       ⎛         ⎛  2     2⎞  ⅈ⋅b₁⋅t       2       2     ⎛  2     2⎞  2⋅ⅈ⋅b₁⋅\n",
       "⎢       ⎝- 2⋅b₁⋅t⋅⎝h₁  + h₂ ⎠⋅ℯ       + ⅈ⋅h₁  + ⅈ⋅h₂  - ⅈ⋅⎝h₁  + h₂ ⎠⋅ℯ       \n",
       "⎢       ──────────────────────────────────────────────────────────────────────\n",
       "⎣                                              4⋅b₁                           \n",
       "\n",
       "                  ⎤\n",
       "₄₁ + rest₁(h₃, h₄)⎥\n",
       "                  ⎥\n",
       "                  ⎥\n",
       "₄₂ + rest₂(h₃, h₄)⎥\n",
       "                  ⎥\n",
       "t⎞  -ⅈ⋅b₁⋅t       ⎥\n",
       " ⎠⋅ℯ              ⎥\n",
       "───────────       ⎥\n",
       "                  ⎦"
      ]
     },
     "execution_count": 66,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "g2(t, h1, h2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's see if the replacement is possible:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$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)$$"
      ],
      "text/plain": [
       "    ⎛    ⎛    3                   ⎞│         2                ⎞           ⎛   \n",
       "  2 ⎜    ⎜   ∂                    ⎟│        ∂                 ⎟           ⎜   \n",
       "h₁ ⋅⎜dt₁⋅⎜────────(g₁(ξ₁, h₁, h₂))⎟│     + ────(g₂(t, h₁, h₂))⎟ - 2⋅h₁⋅h₂⋅⎜dt₁\n",
       "    ⎜    ⎜   2                    ⎟│          2               ⎟           ⎝   \n",
       "    ⎝    ⎝∂h₂  ∂ξ₁                ⎠│ξ₁=t   ∂h₂                ⎠               \n",
       "\n",
       " ⎛      3                    ⎞│           2                 ⎞          ⎛  ⎛   \n",
       " ⎜     ∂                     ⎟│          ∂                  ⎟          ⎜  ⎜   \n",
       "⋅⎜───────────(g₁(ξ₁, h₁, h₂))⎟│     + ───────(g₂(t, h₁, h₂))⎟ + 2⋅h₁⋅ν⋅⎜t⋅⎜───\n",
       " ⎝∂h₂ ∂h₁ ∂ξ₁                ⎠│ξ₁=t   ∂h₂ ∂h₁               ⎠          ⎝  ⎝∂h₂\n",
       "                                                                              \n",
       "\n",
       " 2                  ⎞│                         ⎞       ⎛    ⎛    3            \n",
       "∂                   ⎟│        ∂                ⎟     2 ⎜    ⎜   ∂             \n",
       "────(g₁(ξ₁, h₁, h₂))⎟│     - ───(g₁(t, h₁, h₂))⎟ + h₂ ⋅⎜dt₁⋅⎜────────(g₁(ξ₁, h\n",
       " ∂ξ₁                ⎠│ξ₁=t   ∂h₂               ⎠       ⎜    ⎜   2             \n",
       "                                                       ⎝    ⎝∂h₁  ∂ξ₁         \n",
       "\n",
       "       ⎞│         2                ⎞          ⎛  ⎛    2                  ⎞│   \n",
       "       ⎟│        ∂                 ⎟          ⎜  ⎜   ∂                   ⎟│   \n",
       "₁, h₂))⎟│     + ────(g₂(t, h₁, h₂))⎟ - 2⋅h₂⋅ν⋅⎜t⋅⎜───────(g₁(ξ₁, h₁, h₂))⎟│   \n",
       "       ⎟│          2               ⎟          ⎝  ⎝∂h₁ ∂ξ₁                ⎠│ξ₁=\n",
       "       ⎠│ξ₁=t   ∂h₁                ⎠                                          \n",
       "\n",
       "                      ⎞\n",
       "     ∂                ⎟\n",
       "  - ───(g₁(t, h₁, h₂))⎟\n",
       "t   ∂h₁               ⎠\n",
       "                       "
      ]
     },
     "execution_count": 67,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "simplify(d2f)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$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)$$"
      ],
      "text/plain": [
       "    ⎛    ⎛    3                   ⎞│         2                ⎞           ⎛   \n",
       "  2 ⎜    ⎜   ∂                    ⎟│        ∂                 ⎟           ⎜   \n",
       "h₁ ⋅⎜dt₁⋅⎜────────(g₁(ξ₁, h₁, h₂))⎟│     + ────(g₂(t, h₁, h₂))⎟ - 2⋅h₁⋅h₂⋅⎜dt₁\n",
       "    ⎜    ⎜   2                    ⎟│          2               ⎟           ⎝   \n",
       "    ⎝    ⎝∂h₂  ∂ξ₁                ⎠│ξ₁=t   ∂h₂                ⎠               \n",
       "\n",
       " ⎛      3                    ⎞│           2                 ⎞          ⎛  ⎛   \n",
       " ⎜     ∂                     ⎟│          ∂                  ⎟          ⎜  ⎜   \n",
       "⋅⎜───────────(g₁(ξ₁, h₁, h₂))⎟│     + ───────(g₂(t, h₁, h₂))⎟ + 2⋅h₁⋅ν⋅⎜t⋅⎜───\n",
       " ⎝∂h₂ ∂h₁ ∂ξ₁                ⎠│ξ₁=t   ∂h₂ ∂h₁               ⎠          ⎝  ⎝∂h₂\n",
       "                                                                              \n",
       "\n",
       " 2                  ⎞│                         ⎞       ⎛    ⎛    3            \n",
       "∂                   ⎟│        ∂                ⎟     2 ⎜    ⎜   ∂             \n",
       "────(g₁(ξ₁, h₁, h₂))⎟│     - ───(g₁(t, h₁, h₂))⎟ + h₂ ⋅⎜dt₁⋅⎜────────(g₁(ξ₁, h\n",
       " ∂ξ₁                ⎠│ξ₁=t   ∂h₂               ⎠       ⎜    ⎜   2             \n",
       "                                                       ⎝    ⎝∂h₁  ∂ξ₁         \n",
       "\n",
       "       ⎞│         2                ⎞          ⎛  ⎛    2                  ⎞│   \n",
       "       ⎟│        ∂                 ⎟          ⎜  ⎜   ∂                   ⎟│   \n",
       "₁, h₂))⎟│     + ────(g₂(t, h₁, h₂))⎟ - 2⋅h₂⋅ν⋅⎜t⋅⎜───────(g₁(ξ₁, h₁, h₂))⎟│   \n",
       "       ⎟│          2               ⎟          ⎝  ⎝∂h₁ ∂ξ₁                ⎠│ξ₁=\n",
       "       ⎠│ξ₁=t   ∂h₁                ⎠                                          \n",
       "\n",
       "                      ⎞\n",
       "     ∂                ⎟\n",
       "  - ───(g₁(t, h₁, h₂))⎟\n",
       "t   ∂h₁               ⎠\n",
       "                       "
      ]
     },
     "execution_count": 68,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "simplify(d2f.subs({\n",
    "    g1: Lambda((t, h1, h2), simplify(Matrix(list(xhat(t)) + [0]))),\n",
    "    g2: Lambda((t, h1, h2), simplify(Matrix([exp21, exp22] +[zhat(t)]))),\n",
    "}))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "OK. This failed. But we can copy and paste this and the replacement of $g_1$ and $g_2$ will be effective:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "metadata": {},
   "outputs": [],
   "source": [
    "d2f = (\n",
    "    h1**2*(dt1*diff(g1(t, h1, h2), t, h2, h2)\n",
    "         + diff(g2(t, h1, h2), h2, h2))\n",
    "    - 2*h1*h2*(dt1*diff(g1(t, h1, h2), t, h1, h2)\n",
    "               + diff(g2(t, h1, h2), h1, h2))\n",
    "    + 2*h1*nu*(t*diff(g1(t, h1, h2), t, h2)\n",
    "               - diff(g1(t, h1, h2), h2))\n",
    "    + h2**2*(dt1*diff(g1(t, h1, h2), t, h1, h1)\n",
    "             + diff(g2(t, h1, h2), h1, h1))\n",
    "    - 2*h2*nu*(t*diff(g1(t, h1, h2), t, h1)\n",
    "               - diff(g1(t, h1, h2), h1))\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And the same for $t$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "metadata": {},
   "outputs": [],
   "source": [
    "t = 2 * pi / b1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$\\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]$$"
      ],
      "text/plain": [
       "⎡                                       ⎛    ⎛  ⎛ 2⋅ⅈ⋅b₁⋅t      ⅈ⋅b₁⋅t    ⎞  -\n",
       "⎢      2         2                      ⎜    ⎜  ⎝ℯ         - 2⋅ℯ       + 1⎠⋅ℯ \n",
       "⎢2⋅a⋅h₁  + 6⋅a⋅h₂  - 4⋅b⋅h₁⋅h₂ + 2⋅h₁⋅ν⋅⎜ⅈ⋅t⋅⎜- ──────────────────────────────\n",
       "⎢                                       ⎝    ⎝                   2            \n",
       "⎢                                                                             \n",
       "⎢                                        ⎛  ⎛  ⎛ 2⋅ⅈ⋅b₁⋅t    ⎞  -ⅈ⋅b₁⋅t       \n",
       "⎢                   2         2          ⎜  ⎜  ⎝ℯ         - 1⎠⋅ℯ           ⅈ⋅b\n",
       "⎢-4⋅a⋅h₁⋅h₂ + 6⋅b⋅h₁  + 2⋅b⋅h₂  + 2⋅h₁⋅ν⋅⎜t⋅⎜- ──────────────────────── + ℯ   \n",
       "⎢                                        ⎝  ⎝             2                   \n",
       "⎢                                                                             \n",
       "⎢                                                                 2 ⎛        ⅈ\n",
       "⎢                                                               h₁ ⋅⎝2⋅b₁⋅t⋅ℯ \n",
       "⎢                                                             - ──────────────\n",
       "⎣                                                                             \n",
       "\n",
       "ⅈ⋅b₁⋅t              ⎞   ⎛ 2⋅ⅈ⋅b₁⋅t      ⅈ⋅b₁⋅t    ⎞  -ⅈ⋅b₁⋅t⎞          ⎛  ⎛  ⎛\n",
       "          ⅈ⋅b₁⋅t    ⎟   ⎝ℯ         - 2⋅ℯ       + 1⎠⋅ℯ       ⎟          ⎜  ⎜  ⎝\n",
       "────── + ℯ       - 1⎟ - ────────────────────────────────────⎟ - 2⋅h₂⋅ν⋅⎜t⋅⎜- ─\n",
       "                    ⎠                   2⋅b₁                ⎠          ⎝  ⎝   \n",
       "                                                                              \n",
       "   ⎞   ⎛     2⋅ⅈ⋅b₁⋅t    ⎞  -ⅈ⋅b₁⋅t⎞          ⎛    ⎛⎛ 2⋅ⅈ⋅b₁⋅t      ⅈ⋅b₁⋅t    \n",
       "₁⋅t⎟   ⎝- ⅈ⋅ℯ         + ⅈ⎠⋅ℯ       ⎟          ⎜    ⎜⎝ℯ         - 2⋅ℯ       + 1\n",
       "   ⎟ - ────────────────────────────⎟ - 2⋅h₂⋅ν⋅⎜ⅈ⋅t⋅⎜──────────────────────────\n",
       "   ⎠               2⋅b₁            ⎠          ⎝    ⎝                 2        \n",
       "                                                                              \n",
       "⋅b₁⋅t      2⋅ⅈ⋅b₁⋅t    ⎞  -ⅈ⋅b₁⋅t     2 ⎛        ⅈ⋅b₁⋅t      2⋅ⅈ⋅b₁⋅t    ⎞  -ⅈ\n",
       "      + ⅈ⋅ℯ         - ⅈ⎠⋅ℯ          h₂ ⋅⎝2⋅b₁⋅t⋅ℯ       + ⅈ⋅ℯ         - ⅈ⎠⋅ℯ  \n",
       "───────────────────────────────── - ──────────────────────────────────────────\n",
       "        2⋅b₁                                              2⋅b₁                \n",
       "\n",
       " 2⋅ⅈ⋅b₁⋅t    ⎞  -ⅈ⋅b₁⋅t          ⎞   ⎛     2⋅ⅈ⋅b₁⋅t    ⎞  -ⅈ⋅b₁⋅t⎞ ⎤\n",
       "ℯ         - 1⎠⋅ℯ           ⅈ⋅b₁⋅t⎟   ⎝- ⅈ⋅ℯ         + ⅈ⎠⋅ℯ       ⎟ ⎥\n",
       "─────────────────────── + ℯ      ⎟ - ────────────────────────────⎟ ⎥\n",
       "          2                      ⎠               2⋅b₁            ⎠ ⎥\n",
       "                                                                   ⎥\n",
       "⎞  -ⅈ⋅b₁⋅t              ⎞   ⎛   2⋅ⅈ⋅b₁⋅t      ⅈ⋅b₁⋅t    ⎞  -ⅈ⋅b₁⋅t⎞⎥\n",
       "⎠⋅ℯ           ⅈ⋅b₁⋅t    ⎟   ⎝- ℯ         + 2⋅ℯ       - 1⎠⋅ℯ       ⎟⎥\n",
       "────────── - ℯ       + 1⎟ - ──────────────────────────────────────⎟⎥\n",
       "                        ⎠                    2⋅b₁                 ⎠⎥\n",
       "                                                                   ⎥\n",
       "⋅b₁⋅t                                                              ⎥\n",
       "                                                                   ⎥\n",
       "─────                                                              ⎥\n",
       "                                                                   ⎦"
      ]
     },
     "execution_count": 71,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "d2f"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The replacement does not work apparently, so let's do it manually:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "metadata": {},
   "outputs": [],
   "source": [
    "d2f = Matrix([\n",
    "[ 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))],\n",
    "[-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))],\n",
    "[                                                                                                                                                       -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)]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {},
   "outputs": [],
   "source": [
    "d2Exp = simplify(d2f)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It still works well. And we removed the complex exponential, this result is purely real now!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$\\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]$$"
      ],
      "text/plain": [
       "⎡      2         2               4⋅π⋅h₂⋅ν ⎤\n",
       "⎢2⋅a⋅h₁  + 6⋅a⋅h₂  - 4⋅b⋅h₁⋅h₂ - ──────── ⎥\n",
       "⎢                                   b₁    ⎥\n",
       "⎢                                         ⎥\n",
       "⎢                   2         2   4⋅π⋅h₁⋅ν⎥\n",
       "⎢-4⋅a⋅h₁⋅h₂ + 6⋅b⋅h₁  + 2⋅b⋅h₂  + ────────⎥\n",
       "⎢                                    b₁   ⎥\n",
       "⎢                                         ⎥\n",
       "⎢                 ⎛  2     2⎞             ⎥\n",
       "⎢            -2⋅π⋅⎝h₁  + h₂ ⎠             ⎥\n",
       "⎢            ─────────────────            ⎥\n",
       "⎣                    b₁                   ⎦"
      ]
     },
     "execution_count": 74,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "d2Exp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$\\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)$$"
      ],
      "text/plain": [
       "                                                                              \n",
       "                     2               3            3               2        2  \n",
       "(u₁, u₂, u₅) ↦ 2⋅a⋅h₁ ⋅h₂⋅u₅ + 2⋅a⋅h₂ ⋅u₅ - 2⋅b⋅h₁ ⋅u₅ - 2⋅b⋅h₁⋅h₂ ⋅u₅ - h₁ ⋅h\n",
       "                                                                              \n",
       "\n",
       "                                                                              \n",
       "              2                                                               \n",
       "₃⋅k₁₃₂⋅u₅ - h₁ ⋅h₄⋅k₁₄₂⋅u₅ + h₁⋅h₂⋅h₃⋅k₁₃₁⋅u₅ - h₁⋅h₂⋅h₃⋅k₂₃₂⋅u₅ + h₁⋅h₂⋅h₄⋅k₁\n",
       "                                                                              \n",
       "\n",
       "                                                                   3          \n",
       "                             2                2              2⋅π⋅h₁ ⋅u₂   2⋅π⋅\n",
       "₄₁⋅u₅ - h₁⋅h₂⋅h₄⋅k₂₄₂⋅u₅ + h₂ ⋅h₃⋅k₂₃₁⋅u₅ + h₂ ⋅h₄⋅k₂₄₁⋅u₅ + ────────── - ────\n",
       "                                                                 b₁           \n",
       "\n",
       "  2                  2            3   \n",
       "h₁ ⋅h₂⋅u₁   2⋅π⋅h₁⋅h₂ ⋅u₂   2⋅π⋅h₂ ⋅u₁\n",
       "───────── + ───────────── - ──────────\n",
       "  b₁              b₁            b₁    "
      ]
     },
     "execution_count": 75,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Psi"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So now we can call $\\Psi$ on the three components of this `d2Exp` :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 76,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$- \\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)$$"
      ],
      "text/plain": [
       "                                                                              \n",
       "          2    ⎛  2     2⎞           3 ⎛  2     2⎞           3 ⎛  2     2⎞    \n",
       "  4⋅π⋅a⋅h₁ ⋅h₂⋅⎝h₁  + h₂ ⎠   4⋅π⋅a⋅h₂ ⋅⎝h₁  + h₂ ⎠   4⋅π⋅b⋅h₁ ⋅⎝h₁  + h₂ ⎠   4\n",
       "- ──────────────────────── - ───────────────────── + ───────────────────── + ─\n",
       "             b₁                        b₁                      b₁             \n",
       "\n",
       "                                3 ⎛                   2         2   4⋅π⋅h₁⋅ν⎞ \n",
       "          2 ⎛  2     2⎞   2⋅π⋅h₁ ⋅⎜-4⋅a⋅h₁⋅h₂ + 6⋅b⋅h₁  + 2⋅b⋅h₂  + ────────⎟ \n",
       "⋅π⋅b⋅h₁⋅h₂ ⋅⎝h₁  + h₂ ⎠           ⎝                                    b₁   ⎠ \n",
       "─────────────────────── + ─────────────────────────────────────────────────── \n",
       "          b₁                                       b₁                         \n",
       "\n",
       "        2    ⎛      2         2               4⋅π⋅h₂⋅ν⎞                       \n",
       "  2⋅π⋅h₁ ⋅h₂⋅⎜2⋅a⋅h₁  + 6⋅a⋅h₂  - 4⋅b⋅h₁⋅h₂ - ────────⎟         2         ⎛  2\n",
       "             ⎝                                   b₁   ⎠   2⋅π⋅h₁ ⋅h₃⋅k₁₃₂⋅⎝h₁ \n",
       "- ───────────────────────────────────────────────────── + ────────────────────\n",
       "                            b₁                                         b₁     \n",
       "\n",
       "                                                 2 ⎛                   2      \n",
       "     2⎞         2         ⎛  2     2⎞   2⋅π⋅h₁⋅h₂ ⋅⎜-4⋅a⋅h₁⋅h₂ + 6⋅b⋅h₁  + 2⋅b\n",
       " + h₂ ⎠   2⋅π⋅h₁ ⋅h₄⋅k₁₄₂⋅⎝h₁  + h₂ ⎠              ⎝                          \n",
       "─────── + ─────────────────────────── + ──────────────────────────────────────\n",
       "                       b₁                                         b₁          \n",
       "\n",
       "   2   4⋅π⋅h₁⋅ν⎞                                                              \n",
       "⋅h₂  + ────────⎟                     ⎛  2     2⎞                     ⎛  2     \n",
       "          b₁   ⎠   2⋅π⋅h₁⋅h₂⋅h₃⋅k₁₃₁⋅⎝h₁  + h₂ ⎠   2⋅π⋅h₁⋅h₂⋅h₃⋅k₂₃₂⋅⎝h₁  + h₂\n",
       "──────────────── - ───────────────────────────── + ───────────────────────────\n",
       "                                 b₁                              b₁           \n",
       "\n",
       "                                                                           3 ⎛\n",
       "2⎞                     ⎛  2     2⎞                     ⎛  2     2⎞   2⋅π⋅h₂ ⋅⎜\n",
       " ⎠   2⋅π⋅h₁⋅h₂⋅h₄⋅k₁₄₁⋅⎝h₁  + h₂ ⎠   2⋅π⋅h₁⋅h₂⋅h₄⋅k₂₄₂⋅⎝h₁  + h₂ ⎠           ⎝\n",
       "── - ───────────────────────────── + ───────────────────────────── - ─────────\n",
       "                   b₁                              b₁                         \n",
       "\n",
       "      2         2               4⋅π⋅h₂⋅ν⎞                                     \n",
       "2⋅a⋅h₁  + 6⋅a⋅h₂  - 4⋅b⋅h₁⋅h₂ - ────────⎟         2         ⎛  2     2⎞       \n",
       "                                   b₁   ⎠   2⋅π⋅h₂ ⋅h₃⋅k₂₃₁⋅⎝h₁  + h₂ ⎠   2⋅π⋅\n",
       "───────────────────────────────────────── - ─────────────────────────── - ────\n",
       "               b₁                                        b₁                   \n",
       "\n",
       "                       \n",
       "  2         ⎛  2     2⎞\n",
       "h₂ ⋅h₄⋅k₂₄₁⋅⎝h₁  + h₂ ⎠\n",
       "───────────────────────\n",
       "         b₁            "
      ]
     },
     "execution_count": 76,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Psi(d2Exp[0], d2Exp[1], d2Exp[2])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 77,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$\\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)$$"
      ],
      "text/plain": [
       "   2 ⎛   ⎛        6             4   3          2   5         7         7      \n",
       "4⋅π ⋅⎝b₁⋅⎝- 8⋅a⋅h₁ ⋅h₂ - 24⋅a⋅h₁ ⋅h₂  - 24⋅a⋅h₁ ⋅h₂  - 8⋅a⋅h₂  + 8⋅b⋅h₁  + 24⋅\n",
       "──────────────────────────────────────────────────────────────────────────────\n",
       "                                                                              \n",
       "                                                                              \n",
       "\n",
       "    5   2          3   4            6     6             6             5       \n",
       "b⋅h₁ ⋅h₂  + 24⋅b⋅h₁ ⋅h₂  + 8⋅b⋅h₁⋅h₂  + h₁ ⋅h₃⋅k₁₃₂ + h₁ ⋅h₄⋅k₁₄₂ - h₁ ⋅h₂⋅h₃⋅\n",
       "──────────────────────────────────────────────────────────────────────────────\n",
       "                                                                              \n",
       "                                                                              \n",
       "\n",
       "         5                5                5                  4   2           \n",
       "k₁₃₁ + h₁ ⋅h₂⋅h₃⋅k₂₃₂ - h₁ ⋅h₂⋅h₄⋅k₁₄₁ + h₁ ⋅h₂⋅h₄⋅k₂₄₂ + 2⋅h₁ ⋅h₂ ⋅h₃⋅k₁₃₂ - \n",
       "──────────────────────────────────────────────────────────────────────────────\n",
       "                                                                              \n",
       "                                                                              \n",
       "\n",
       "  4   2               4   2             4   2               3   3             \n",
       "h₁ ⋅h₂ ⋅h₃⋅k₂₃₁ + 2⋅h₁ ⋅h₂ ⋅h₄⋅k₁₄₂ - h₁ ⋅h₂ ⋅h₄⋅k₂₄₁ - 2⋅h₁ ⋅h₂ ⋅h₃⋅k₁₃₁ + 2⋅\n",
       "──────────────────────────────────────────────────────────────────────────────\n",
       "                                                           3                  \n",
       "                                                         b₁                   \n",
       "\n",
       "  3   3               3   3               3   3             2   4             \n",
       "h₁ ⋅h₂ ⋅h₃⋅k₂₃₂ - 2⋅h₁ ⋅h₂ ⋅h₄⋅k₁₄₁ + 2⋅h₁ ⋅h₂ ⋅h₄⋅k₂₄₂ + h₁ ⋅h₂ ⋅h₃⋅k₁₃₂ - 2⋅\n",
       "──────────────────────────────────────────────────────────────────────────────\n",
       "                                                                              \n",
       "                                                                              \n",
       "\n",
       "  2   4             2   4               2   4                5                \n",
       "h₁ ⋅h₂ ⋅h₃⋅k₂₃₁ + h₁ ⋅h₂ ⋅h₄⋅k₁₄₂ - 2⋅h₁ ⋅h₂ ⋅h₄⋅k₂₄₁ - h₁⋅h₂ ⋅h₃⋅k₁₃₁ + h₁⋅h₂\n",
       "──────────────────────────────────────────────────────────────────────────────\n",
       "                                                                              \n",
       "                                                                              \n",
       "\n",
       "5                5                5             6             6        ⎞      \n",
       " ⋅h₃⋅k₂₃₂ - h₁⋅h₂ ⋅h₄⋅k₁₄₁ + h₁⋅h₂ ⋅h₄⋅k₂₄₂ - h₂ ⋅h₃⋅k₂₃₁ - h₂ ⋅h₄⋅k₂₄₁⎠ + 4⋅π\n",
       "──────────────────────────────────────────────────────────────────────────────\n",
       "                                                                              \n",
       "                                                                              \n",
       "\n",
       "   ⎛  6       4   2       2   4     6⎞⎞\n",
       "⋅ν⋅⎝h₁  + 3⋅h₁ ⋅h₂  + 3⋅h₁ ⋅h₂  + h₂ ⎠⎠\n",
       "───────────────────────────────────────\n",
       "                                       \n",
       "                                       "
      ]
     },
     "execution_count": 77,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "simplify(expand(Psi(d2Exp[0], d2Exp[1], d2Exp[2]) / (1 / (1/b1 * 2 * (h1**2 + h2**2) * pi))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We don't have the value for `nu` !"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Conclusion"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We do NOT obtain the same result as the document. Everything failed at the end.\n",
    "\n",
    "Too bad, but still, it was interesting. I guess?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> See [here](https://github.com/Naereen/notebooks) for other notebooks I wrote."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.5"
  },
  "toc": {
   "colors": {
    "hover_highlight": "#DAA520",
    "running_highlight": "#FF0000",
    "selected_highlight": "#FFD700"
   },
   "moveMenuLeft": true,
   "nav_menu": {
    "height": "67px",
    "width": "252px"
   },
   "navigate_menu": true,
   "number_sections": true,
   "sideBar": false,
   "threshold": 4,
   "toc_cell": true,
   "toc_position": {
    "height": "42px",
    "left": "967.188px",
    "right": "20px",
    "top": "120px",
    "width": "159px"
   },
   "toc_section_display": "block",
   "toc_window_display": true
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "position": {
    "height": "382px",
    "left": "782px",
    "right": "20px",
    "top": "120px",
    "width": "350px"
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}