piston.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# piston example: some symbolic computations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sympy\n",
    "sympy.init_printing()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## GENERIC formulation\n",
    "\n",
    "$$\n",
    "\t\t\\begin{bmatrix}\n",
    "\t\t\t\\dot{q} \\\\\n",
    "\t\t\t\\dot{p} \\\\\n",
    "\t\t\t\\dot{s}_1 \\\\\n",
    "\t\t\t\\dot{s}_2\n",
    "\t\t\\end{bmatrix}\n",
    "\t\t\\: = \\:\n",
    "\t\t\\begin{bmatrix}\n",
    "\t\t\t0 && 1 && 0 && 0 \\\\\n",
    "\t\t\t-1 && 0 && 0 && 0 \\\\\n",
    "\t\t\t0 && 0 && 0 && 0 \\\\\n",
    "\t\t\t0 && 0 && 0 && 0 \\\\\n",
    "\t\t\\end{bmatrix}\n",
    "\t\t\\,\n",
    "\t\t\\begin{bmatrix}\n",
    "\t\t\t\\left(\\pi_2 - \\pi_1\\right) A \\\\\n",
    "\t\t\tv \\\\\n",
    "\t\t\t\\theta_1 \\\\\n",
    "\t\t\t\\theta_2\n",
    "\t\t\\end{bmatrix}\n",
    "\t\t\\: + \\:\n",
    "\t\t\\begin{bmatrix}\n",
    "\t\t\t0 && 0 && 0 && 0 \\\\\n",
    "\t\t\t0 && 0 && 0 && 0 \\\\\n",
    "\t\t\t0 && 0 && \\alpha \\, \\frac{\\theta_2}{\\theta_1} && -\\alpha \\\\\n",
    "\t\t\t0 && 0 && -\\alpha && \\alpha \\, \\frac{\\theta_1}{\\theta_2}\n",
    "\t\t\\end{bmatrix}\n",
    "\t\t\\,\n",
    "\t\t\\begin{bmatrix}\n",
    "\t\t\t0 \\\\ 0 \\\\ 1 \\\\ 1\n",
    "\t\t\\end{bmatrix}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "θ_1 = sympy.Symbol('θ_1', real=True, positive=True)\n",
    "θ_2 = sympy.Symbol('θ_2', real=True, positive=True)\n",
    "π_1 = sympy.Symbol('π_1', real=True, positive=True)\n",
    "π_2 = sympy.Symbol('π_2', real=True, positive=True)\n",
    "υ = sympy.Symbol('υ', real=True)\n",
    "\n",
    "A = sympy.Symbol('A', real=True, positive=True)\n",
    "α = sympy.Symbol('α', real=True, positive=True)\n",
    "d = sympy.Symbol('d', real=True, positive=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "dE = sympy.Matrix([\n",
    "    [(π_2 - π_1) * A],\n",
    "    [υ],\n",
    "    [θ_1],\n",
    "    [θ_2]\n",
    "])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "dS = sympy.Matrix([\n",
    "    [0],\n",
    "    [0],\n",
    "    [1],\n",
    "    [1]\n",
    "])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "J = sympy.Matrix([\n",
    "    [0, 1, 0, 0],\n",
    "    [-1, 0, 0, 0],\n",
    "    [0, 0, 0, 0],\n",
    "    [0, 0, 0, 0]\n",
    "])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "assert J.is_anti_symmetric()\n",
    "assert (J @ dS).applyfunc(sympy.simplify).is_zero"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$\\left[\\begin{matrix}0 & 0 & 0 & 0\\\\0 & 0 & 0 & 0\\\\0 & 0 & \\frac{α θ_{2}}{θ_{1}} & - α\\\\0 & 0 & - α & \\frac{α θ_{1}}{θ_{2}}\\end{matrix}\\right]$$"
      ],
      "text/plain": [
       "⎡0  0   0     0  ⎤\n",
       "⎢                ⎥\n",
       "⎢0  0   0     0  ⎥\n",
       "⎢                ⎥\n",
       "⎢      α⋅θ₂      ⎥\n",
       "⎢0  0  ────   -α ⎥\n",
       "⎢       θ₁       ⎥\n",
       "⎢                ⎥\n",
       "⎢            α⋅θ₁⎥\n",
       "⎢0  0   -α   ────⎥\n",
       "⎣             θ₂ ⎦"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# heat conduction\n",
    "R_α = sympy.Matrix([\n",
    "    [0, 0, 0, 0],\n",
    "    [0, 0, 0, 0],\n",
    "    [0, 0, α*θ_2/θ_1, -α],\n",
    "    [0, 0, -α, α*θ_1/θ_2]\n",
    "])\n",
    "R_α"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "assert R_α.is_symmetric()\n",
    "assert all([val >= 0 for val in R_α.eigenvals()])\n",
    "assert (R_α @ dE).applyfunc(sympy.simplify).is_zero"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$\\left[\\begin{matrix}0 & 0 & 0 & 0\\\\0 & \\frac{d \\left(θ_{1} + θ_{2}\\right)}{2} & - \\frac{d υ}{2} & - \\frac{d υ}{2}\\\\0 & - \\frac{d υ}{2} & \\frac{d υ^{2}}{2 θ_{1}} & 0\\\\0 & - \\frac{d υ}{2} & 0 & \\frac{d υ^{2}}{2 θ_{2}}\\end{matrix}\\right]$$"
      ],
      "text/plain": [
       "⎡0       0         0      0  ⎤\n",
       "⎢                            ⎥\n",
       "⎢   d⋅(θ₁ + θ₂)  -d⋅υ   -d⋅υ ⎥\n",
       "⎢0  ───────────  ─────  ─────⎥\n",
       "⎢        2         2      2  ⎥\n",
       "⎢                            ⎥\n",
       "⎢                   2        ⎥\n",
       "⎢      -d⋅υ      d⋅υ         ⎥\n",
       "⎢0     ─────     ────     0  ⎥\n",
       "⎢        2       2⋅θ₁        ⎥\n",
       "⎢                            ⎥\n",
       "⎢                          2 ⎥\n",
       "⎢      -d⋅υ             d⋅υ  ⎥\n",
       "⎢0     ─────       0    ──── ⎥\n",
       "⎣        2              2⋅θ₂ ⎦"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# friction between piston and cylinder\n",
    "# dissipated kinetic energy goes to internal energy of both compartments with equal shares\n",
    "R_d = sympy.Matrix([\n",
    "    [0, 0, 0, 0],\n",
    "    [0, d*(θ_1+θ_2)/2, -d*υ/2, -d*υ/2],\n",
    "    [0, -d*υ/2, d*υ**2/(2*θ_1), 0],\n",
    "    [0, -d*υ/2, 0, d*υ**2/(2*θ_2)]\n",
    "])\n",
    "R_d"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "assert R_d.is_symmetric()\n",
    "assert (R_d @ dE).applyfunc(sympy.simplify).is_zero"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABmoAAAAtCAYAAACuyTQAAAAABHNCSVQICAgIfAhkiAAAGWVJREFUeJzt3XuQNUV9xvHvIhdfRe4B0aBvySsY44tWhLyAEgdMFBJRy1JSlliuKbXKBGJI8BKMBjVqItGAUglgxFFD1GgCYiLiBbfUGIwhojGiBvQgV/GGF0RekTd/9Iw7O2dmTk93z/TMmedTdWp355ye7Z7u6enf6bmAiIiIiIiIiIiIiIiIRLFSsWwbcBFwC3BHtuxk4Dt9ZUpERERERERERERERGSJvBw4Jvv9wcAa8EKAnSs+vAk4GDgduKSHzIkMyS7A9tiZEBERkaVSdXKUiCw3xRUiIiISmuKK8Xtt4fcZcGD+R9VEjciUJZgryC6KnA8RERERERmvBMUVIiIiImJpp9gZEBmYE4APxc6EiIiIiIiMmuIKEREREbGmiRqRjfYCvu+Qbifgc8D7W6Y7CHMvwi8DXwCe5vC/+/QZTD6/BLwyQvq+xMrnENpDqDy47hN98i2rT/oh1HWofIyhrn349AdDqeehi9mGtA+ISFdc4wqYRmwxlbgC4uR1KG1hKrFFzLgiRPoQNKay49ofDKGOx2DscQVMYz+QFhJgB/DUyPkQ6dthwGmOaf8AeA/tO9IDgUdmv+8P3ADcxzEPfdgj+3kv4D+AR/Wcvi+x8jmE9hAqD677RJ98y+qTfgh1HSofY6hrHz79wVDqeehitiHtAyLSBZ+4AqYRW0wlroA4eR1KW5hKbBEzrgiRPgSNqey49gdDqOMxGHtcAdPYD2TdDLgk/0NX1IitKczenwhcWvPeKcC1wJ3AZcB+hff2x2yPCxrWXZf+Fsz2BLgNc9bdfnOp+9NUToAfZj93zV6h0/epKa82+VxUVhd9tgffNum7T/SlKZ82ZfVNX2cIdR0iH2Opa5+0Pv3W2Pr4rtI26aMNTWUfEBkLxRXN/emifqXLPi20LsfbY4krwD+GcjGmsabvPtGnrmIo3/rqq76nMqbqcszs2neNqX/vMm2TvtpQl9/vDWk/kIFI0BU1tn4FuH/sTPSkbmZ4F+AxsTJl4RTgGsvPXliz/NXA14FjMdvg2tJn3wEcjdl3qma8F6XPHZ7ldcUyv7b+FHPZ5A+BbwMfBB7hkc/PAj8C/ipw+hBClrUpn7Zl9dFVewD/Num7T4R0BuaYdW7Fe23qqaqsvultDaGuXfPRdV2H7r+q+PYHtsbSx1fpss/rug2NfR+QaZlKbNF0xumQY4s+4gpo7lf6Gps0GdJ4eyxxRVNepxJX1OXDd58ILURs4RpDLUpvq6v6HvuYakhjZt++awz9e50xxxXQ/X6s2GJ6ZhSuqKmS0G6ipjxg3Rv4FnBwxee6mDFtwyUP7wf+uGL5fYFTA/2PkOn78kXgQdnvTwceEjEvTd4KvNvic/cHXlux/NHAz4EjCstOw8yOA/wGkGa/J8x3pIvS5/bFnFV4tEVecylwpsXnLgeeiznAbgUuBm4F9nHIZ26PbL35Qds3fSihy1qVz7ZlddG2PaTYtQXwb5O++0STFPtyABwJfANz1ko5mGpTT1Vl9U1vq8t9v+sy+NS1rS76ryLf/sCWTxtZZEjbqC3XNpQynX1gamzjivyzMcfMIeMK6Ca2iL2NbBXjChhubNF1XAHN/UqXfVqK/RhsaOPtscQVVXkdYlwB/R1nffeJJint4goIE1u4xlCL0tvqat9fhjHV0MbMrn2X4opqffQVXe/Hii2maUbgiZrygPUszAxgUR9niSzimofDgO8Be5aWn8HGgMLnf4RK35fyzPAuwJvjZafRfwIvtfjc84GjKpa/B/hEadkzgXuy318G3ITZsW4F7gDe1iI9wG7AJ4FnW+SzKKX94BNgd8zB5cTCMpt8lr0UOD1Q+q6EKGs5ny5lzT0dc3ZKvu88ELiLjX2JS3tIsW8Lvm3Sd59okmJfjj2B64DjMLdPKQdTtvVUV1bf9LC4vrve933LsCj/PnXtKlT/5Zq23B90tU/7iL2NirpqQyn97APQXIYY+8Cys4krIP6YOXRcAeFji9jbyFbVGadDjS26jiuguV8J0afVSXGLK2AY4+2xxBXgH0PluhyDpPRznPXdJ5qktGvTIWILnxhqUfrYceQyjqmGMGZuG1sorogXV0D3+7Fii2maEXiipjhgvQ9wO3BM4f2uZkxTws961rkK8zCn3ArmEsCQ/2MI28lG3cxwSvyz9A4FrgB+CvwP8OvAT4AnWqR9K/PPbNotS/+HpeW/h7nXZFnCxhlvm/QrmC8kzrTIY1nqmO5AzD6e31bCtpx7Ar+U/X5v4OPAkwKk75JLWZvy2bZNlL0G+HRp2eeA389+d20PqWUa3zbpu08sktb83yrvZf1S8TU2BlO2+awrq2/6XFN9d73vhyjDovZalNDPGT+u/VcV3/4AutunffS9jZp01YZS+tkHwL4MCTrrLYRFcQXEHzOHjisgfGwRexvZajrjNCVubBE7roCN/UqoPq1O6pAmF2O8PZa4AvxjqDpdjkFSi3R9jbeLErq7osY3tvCNoXzHI13u+8s6pooxZvaJLRRXxIsroJ/9uChBscVUzPCYqFk0YH0G5iyx4hlRPjOmTVLCzno2+XPg3wt/Hw6cH/h/xN5OvmdonAac7JrJAA4BfgCcA2wBfge4HtOWD8g+szPw18BJpbSbgLdUrPOoLP2dwI8Lr7uAKys+n7CxI7VJ/1hMHV9deG1dWFojxe0g/V7g88C9WuQTYDPmy4UvAl8CXhkofZdcytqUz7ZtouxS5s/O+m/MWbTg3h5S7NqCb5v03SdCleP5mDrKH764xsbtapvPurL6ps811XfX+36IMixqr0UJ/QwkXfuvKr79AXS3T/voexs16aoNpfSzD4B9GRIUTLXlEldA/DFz6LgCwscWsbdRiDO/Y8YWQ4grYGO/EqpPq5Pi/uVfjPF2U/quxIqh6nQ5Bknp/jjru0+EKEMuRGzhG0P5jke63PeXdUwVY8y8GffYQnFFvLgC+tmPixIUW0zFDMeJGpsB6znARwtpfGdMm6SEm/Vc5HhgO2bgDfAc4M8C/o8hbCffMzSeAby+de7C+Sjz94y+ELil8PdvAh8A3lf63JOBp1Ss8zmYLw+2lF6fxwRmi/imLzuDjQeEn2HaZXFZ+azTsrMw22RLwHyGLueZmH6l6ZVYrCdGWQ/CDOy/jLm38dNK6a/HBAK5nYDvMB/kL+LaFoZW1y7lOBTzYMOHFZatsXHAM5RyxqzvEGUIlf8zibtPX4o5lpYHutpG/ttoyH3e0PZj2cglroD4Y+Yu4goIG1vE3kYQ5szvmLGF4gq7uAKmE1vEKmfTcTbk8SnGcXYobXossYXGVOvOJO4+XTVu1jYKs4366PNifZcSsgyyXGY4TtTYDFg/wMb7SNvOmNZ9QdAkJcys56IvGcDcT3oH6w8yfTEbgx/fcobIY50Uu+3ke4bGE4h378SDMNvviNLyCzAPTi3aAvwIE8Tm3oK5vUbZqcCNpWUHAHcD27K/m9ruovRt63UfNh4Q/hlz/+7isk21qeGNmAfyPrxlPnN1ZbVJ36as+2EGyk2vqvoq8imrT50eiLkPPMD+wA2FvO7FfDs9FrOfV92rvolrW1iGNr2K2Y53F147MP3T3Zh9ewhtOnZ9+26DkPmPvU8fi7mvctt67rOOx7qNhtznDW0/lo1c4gqwG3O7xBUQ7qziRe2qHFeAW2zhGlfY5LFOSj9xBcSLLWLGFdDfeNs3roB44+227df3GBszhqo7zoY+PsUYaw6lTa8SJrbosk2PYUwFzdtgSONm3326atw8tHHnWLdRH31eV9+lKK4QVzMcJmpsB6yXs/GyfdsZx7ovCIq6mvVs+pIh91BM+R+R/X068KKA5QyRx5zrdvKd2T0ecz/mGJ6M6Rx3LS2/iuoz8f4XOCH7fYX6fD8Tc9uG4j2mz8F0vrmmtrsofZt6rZJifzn3m6k+0NrkM1dXVpv0vmVtw7esPnVa9kXWb/Pxa5h+ZI/s7xXMPWn/piZtGyl2bWEZ2vRemL64+Poc8I/Z7ysW+cx12aZj17fvNugy/22F6L8S2tdzzDpuK9Y2KhtSnzf0/XjKXOMKsBtz28QV0N1ZxYvaVTmuALfYwjWusMljLlZcAfFii5hxBcQbb6e0u/VZzPH2mOIK8KvTsvw42/XxKaX7seZQ2nSo2KLLNj2GMRU0b4OhjKlC7NMwP25epnFnzG1U1keflxLmu5Sx1K8MzwyHiRrbAetFwD8V/radbc3/b1NA1eUZ5EXFLxly29h4K4ZnA69w/B8JbmeTLMpjzmU7hZjZPQn4C8vPhnYiJv/3Kyx7XLbsdys+/3rgvOz3bcALatZ7AOZsxFdg7iX6EsxDWLeUPpdQ3XZt0+ea6rVKit0B5W8xB5TjgPsXXrs75DNhvqxtywnty2orVFkT/Ov0cOAa1u/P/quYNvkUzC1f3pa9v3tF2rZS7NrCsrTpsjU2nrk7hDYdu759t0GX+W+jy33atz1rG200tD5v6PvxlLnGFWA/Zk5YPFHT5VnFReV2VY4rwD22SPCPK6rymIsVV0C82CJ2XAFxxtsp9mOwIY23YfhxBYSp0+JxtuvjU0r3Y80htemyNdxii4Ru2vSYxlQJ1dtgCGOqLvfpZRl3Dmkb9dXnpYT9LiU3xPqVYZrhMFFjO2A9HfNArJzPbOsiKeHOIM+Vv2TIPQ+4qfD3YZjbM7j8jwS/s2nq8lgnZfF2CjGz+xLMg0NjeADmzMHzMbeRKN7n/NCKzx8J3Iwp52swM991TsrW9RPgw2y8d20uob7t2qSH9vUK9vtA3T1Ji2lt85lQXVbb9OBWVluhyprgV6f7Yi55Pbq0/O8x9+T/LvB2zJkWIaTYByLL0KbL1pi/xcoQ2nTs+vbdBl3lv42u92nf9qxtZAy1zxvyfty1X8aMVW/GfEE+A84G9o6Yp5xrXAF+V4gtkhLurOJcVbsqxxXgHlsk+J+l3+Z4ndJPXAHxYovYcQXEGW+n2I/BhjLehnHEFeBfp1XH2S6PTyn9jDWH0qbL1nCLLRK6a9NjGVMl1G+D2GOqrvfpZRh3DmUb9dnnpYT9LgWGW7+uhhxXQH27rXutRsllvRkOEzW2A9atwM8x90OEMDPvdVLCznrWfckA8E7m75H8rw7/A/zOJmnKY52UxdspxMzuu1gPyGJ4FvBNzD1R14BXA3ewMUjNrWDugb6NMPe+Tmj/ZUCRS73GkjCNsia4l3M34JOYs2PHKmEa9QzTKmudBL9tMAYJ7mVchjq2keC2jZahz1u2Oj4Yc8uKHZhB/18CV2R/fwVT3phc4wrwP6u5SUrYs4rr2lVVXAFusUWC31nqbdt+Sj9xBcSNLWLGFTCdsUnCNMoJii0SplHXCdMoZ5OE5Y8rQLGFjYRpxhbLVr9DjyvAjE3Lr9sxeTy74r1H9Zq7xWY4TNSA/YD1M8Aphb9DzLxXSQk369nUEWzCnJV2ZGn5H2GCD9v/kUtwO5vEtbNKsdtOPjO7m4A3tMxXbBdg7jv78gDrSpjOwDthGmVNcCvnCubhyGeGzEwECdOoZ5hWWeskLH9AlTDNQKGNhPbbaBn6vNh1/DbMQ+QfGnCdl2PG8qeWlr8pW37eXIr+ucYV4H9Wc52UcGcV17WrurgC3GKLBPez1F3afkr3cQWEiS1WMe098VyPjZBxBUxnbJIwjXKCYouEadR1wjTK2SRh+eMKUGxhI2F6scUQ6jd0bDGGuKLKDJO/zXGzYWWG40SNrScCXwPu1TJdQpwOfVFHcArwkYrlu2ICqrYSlq+zehbNl/kP0ZMw7fywAOtKmM7AO2EaZU1wK+djgXuAqwuvreGy1ZuEadQzTKusdRKWO6D6GPBtzBeWNwJHWaZbpjpexHUbjb3PG0Idf431y/CvBd4C/DbVz1208ZBsXd9gftLjfpiHv98B3Ndx/X0bW1wBze2qLq4At9giYTmPYSFii1X6m6gJGVfA8tZrWcI0ygmKLRKmUdcJ0yhnk4TljitAsYWNKcYWQ6nfkLHFmOOKGc0TNZuz91PMiVLvxVyBfg+mH0uYv+1fef2zmve2YfrBW4HtwA2YOwk8oGFdnU7UgJm9e3CLz7vuxCEs6gheQPX9iMFUbJsgYhk7q10w96Iem3szf99zFz5td8j1WmUqZY3ZHw3BVOoZplXWOlNv702WpY6l3hDqeBMmeDoXuI71wOpO4DLanxH3vCz9+TXv52fFPd4xvzGMKa6A5nbVFFdAu9hiWY9hoWKLVfqbqAkVV8Dy1mvZVMoJ8fuk2KZS11MpZ5Opt/VFlqWepdpQ6jdkbDHmuGKG3UTNpzBX+H8W8zzF8zDj0AS3iZrnAndjJrDejblC/GLM7ZxvBh5Us67OJ2pExmAsD+4SERGRaTgUOA34KOZhneUz4rYtSH9W9vk/qXn/3Oz9F4bIrMiArdLfRA0orhAREZHh8YktxhxXzLCbqNkBvK7i/YT2EzWHYK6guRZ4YOm94zCTNRfXrOsXEzVVD0QUmYrbYmdAREREpOCrmLO5fgvYB3gR8D3MgzxPYf7+0GV7Zj9/UPN+vnwvv2yKSIniChERERkan9hiCnHFt4BXBVrXCzFXhr8IuKn03hXApcCJmNvG1do5UGZE2prR7jYWFwEnd5MVERERkUHYFTgGOCF7PTxb/jPgMyy4LN7CSvZzh+d6RIZkRn1c8YmKZe/AXHEjIiIissy6jC2WIa74AuZKoxDy2z0+Djii4v39Mc/dPAS4qm4lmqiRWK4Dftri8zd3lRERERGRiDYDx2OCp8ez/kDOm4ELgQ9hblfwQ4t15We27Vnz/h6lz4ksg7OZP5vzUcBTMJMys9J7V/eQJxEREZEYNhMmtphCXHFrwHXtm/188YLP7d70pu1EzZhnxySulZrlIR82pfYpIiIiQ1Y3HlrDnHUF5sGTV2KCp8tw+zL5q9nPQ2rezx8e+jWHdYeicZu4qtuPzq5YtoqZqEkx+5kttU8REREZsrrxEISNLcYQV/iqG/fdk/2smzfZk/kJquLEls0JdpVsJ2qaGoFIbGqfIiIiMka7Yb5Ivgz4CHC75/ry2zw9AfMsynsK790PeAxwJyZoi0XjNhkytU8REREZq5CxxRjiiq58P/t5UMV7WzBXcpcnaq4EHo251dy/hcxMgplRemrIlY7YKcC1mMZ3GbBfT2lFRKQb6tdFZJldjhnLlx8O+qZs+XkVadLsvdUuMyY6/vRoFdOmk7jZEJElp35dRJbZWOOKWZaHzTXvb87eT2ve3wUzEXM75tkyuU2YK5R2MH9r3YcB2zFXGFVdhZQ/L6gqr43PCkqY1kTNGZjynlvx3quBrwPHAo/EHEQvtFyvT1oREXGnfl1ExuLjwFdavN5gsc6DgW9h+sFLgNcDV2R/f5X1+ycXvTN7/2T3ogg6/gzJKpqoERF/6tdFZExCxxZjjStm+E3UgOmndwA3YY4B52Get/6pbNmsIs3JmMmanwEfBN4IvBmz7b6L2eZVedVETeZI4BvAF5g/8D4a+DlwRGHZacBtFuv1SSsiIu7Ur4vImMww427b1z9Yrvcg4O3ALZhg4XrgHGCfms9/HnMv5b0dyiCGjj/DsoomakTEj/p1ERmbGeFjizHGFTP8J2pWgJdhJme2A9/ETGzdJ1v/rCbd1my91wN3Ad8DvgScDxxXk1dN1GAe7nMdZiOtMX/gfQ/r9+PLPZON9+Sr45NWRETcqF8XEWlvL8wXRjZX60g1HX9ERJaL+nURkfYUV7Q3ozBRs1O8fER3AfB+zCVbZbsBTwYuLi3fxPzDgkKmFRERd+rXRUTaOwZzef6bYmdkxHT8ERFZLurXRUTaU1zRgYTlv6Lm+cBVmAf5wPwZEkdhtsGdwI8Lr7uAKxes2yetiIi4Ub8uIiIx6PgjIrJc1K+LiEhfZhSuqNk5Xj6iORR4HWaWb3vNZw7BHCi3lpa/D/h09vul2To+Djy9RdqDgHcB+2NmGV8F/ItDOURExFC/LiIiMYQ6/kD1McgmrY5BIiLhqF8XEZFBSVjuK2pWMeW7u/Dagbkf6N2YS1FPBW4spTsge39b9vexwImYy2GLFqU9EHhktnx/4AbMg4hERMTNKurXRUSkf6uEOf5A9THIJq2OQSIi4ayifl1ERPozY+JX1FwC/Fdp2duB/8OcObEd+A6wB+YZPvkD3c7AnOHw2ezvT2Amtcps0t6S/bwN+D6wH/BNx/KIiEyd+nUREYkh1PEHqo9BtscfHYNERMJQvy4iIoOSsNxX1FRZY+M9Rw/A3DP0FcBm4CWYA+SWUrqE+TOvbdMCHA5cA6y4Z11ERCqsoX5dRET6t4bb8Qfmj0Ft0oKOQSIiXVhD/bqIiHRjRuGKmioJmqgBOAm4HvgJ8GHgYRXpEua/0LNNuy/wZeBolwyLiEijNdSvi4hI/9ZwO/5A9THINq2OQSIi3VhD/bqIiHRjhiZqgkmo/kJvkd2ATwLPDpobERHxlaB+XURE4kjQMUhEZJkkqF8XEZF6MwoTNTvFy8dkrQApcAXwrrhZERGRANSvi4hILDoGiYgsF/XrIiLyC9uAa4FPYS7D/DCwd9QcDdPHgG9jLle9ETjKMt1jMQ+Nu7rw2tpFBkVEpBX16yIiEouOQSIiy0X9uoiIVHk563Mu1wB/Fzc7IiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIjIx/w9+Jj52PEpTvgAAAABJRU5ErkJggg==\n",
      "text/latex": [
       "$$\\left [ \\frac{d \\left(θ_{1} + θ_{2}\\right) \\left(θ_{1} θ_{2} + υ^{2}\\right)}{4 θ_{1} θ_{2}} - \\frac{d \\sqrt{θ_{1}^{4} θ_{2}^{2} + 2 θ_{1}^{3} θ_{2}^{3} - 2 θ_{1}^{3} θ_{2} υ^{2} + θ_{1}^{2} θ_{2}^{4} + 4 θ_{1}^{2} θ_{2}^{2} υ^{2} + θ_{1}^{2} υ^{4} - 2 θ_{1} θ_{2}^{3} υ^{2} - 2 θ_{1} θ_{2} υ^{4} + θ_{2}^{2} υ^{4}}}{4 θ_{1} θ_{2}} \\geq 0, \\quad \\frac{d \\left(θ_{1} + θ_{2}\\right) \\left(θ_{1} θ_{2} + υ^{2}\\right)}{4 θ_{1} θ_{2}} + \\frac{d \\sqrt{θ_{1}^{4} θ_{2}^{2} + 2 θ_{1}^{3} θ_{2}^{3} - 2 θ_{1}^{3} θ_{2} υ^{2} + θ_{1}^{2} θ_{2}^{4} + 4 θ_{1}^{2} θ_{2}^{2} υ^{2} + θ_{1}^{2} υ^{4} - 2 θ_{1} θ_{2}^{3} υ^{2} - 2 θ_{1} θ_{2} υ^{4} + θ_{2}^{2} υ^{4}}}{4 θ_{1} θ_{2}} \\geq 0, \\quad \\mathrm{True}\\right ]$$"
      ],
      "text/plain": [
       "⎡                                _____________________________________________\n",
       "⎢            ⎛         2⎞       ╱   4   2       3   3       3     2     2   4 \n",
       "⎢d⋅(θ₁ + θ₂)⋅⎝θ₁⋅θ₂ + υ ⎠   d⋅╲╱  θ₁ ⋅θ₂  + 2⋅θ₁ ⋅θ₂  - 2⋅θ₁ ⋅θ₂⋅υ  + θ₁ ⋅θ₂  \n",
       "⎢──────────────────────── - ──────────────────────────────────────────────────\n",
       "⎣        4⋅θ₁⋅θ₂                                                              \n",
       "\n",
       "____________________________________________________________                  \n",
       "      2   2  2     2  4          3  2            4     2  4                   \n",
       "+ 4⋅θ₁ ⋅θ₂ ⋅υ  + θ₁ ⋅υ  - 2⋅θ₁⋅θ₂ ⋅υ  - 2⋅θ₁⋅θ₂⋅υ  + θ₂ ⋅υ        d⋅(θ₁ + θ₂)⋅\n",
       "──────────────────────────────────────────────────────────── ≥ 0, ────────────\n",
       " 4⋅θ₁⋅θ₂                                                                  4⋅θ₁\n",
       "\n",
       "                    __________________________________________________________\n",
       "⎛         2⎞       ╱   4   2       3   3       3     2     2   4       2   2  \n",
       "⎝θ₁⋅θ₂ + υ ⎠   d⋅╲╱  θ₁ ⋅θ₂  + 2⋅θ₁ ⋅θ₂  - 2⋅θ₁ ⋅θ₂⋅υ  + θ₁ ⋅θ₂  + 4⋅θ₁ ⋅θ₂ ⋅υ\n",
       "──────────── + ───────────────────────────────────────────────────────────────\n",
       "⋅θ₂                                                               4⋅θ₁⋅θ₂     \n",
       "\n",
       "_______________________________________________          ⎤\n",
       "2     2  4          3  2            4     2  4           ⎥\n",
       "  + θ₁ ⋅υ  - 2⋅θ₁⋅θ₂ ⋅υ  - 2⋅θ₁⋅θ₂⋅υ  + θ₂ ⋅υ            ⎥\n",
       "─────────────────────────────────────────────── ≥ 0, True⎥\n",
       "                                                         ⎦"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "[val >= 0 for val in R_d.eigenvals()]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$\\left[\\begin{matrix}υ\\\\A \\left(π_{1} - π_{2}\\right) - d υ\\\\\\frac{\\frac{d υ^{2}}{2} - α θ_{1} + α θ_{2}}{θ_{1}}\\\\\\frac{\\frac{d υ^{2}}{2} + α θ_{1} - α θ_{2}}{θ_{2}}\\end{matrix}\\right]$$"
      ],
      "text/plain": [
       "⎡        υ         ⎤\n",
       "⎢                  ⎥\n",
       "⎢A⋅(π₁ - π₂) - d⋅υ ⎥\n",
       "⎢                  ⎥\n",
       "⎢   2              ⎥\n",
       "⎢d⋅υ               ⎥\n",
       "⎢──── - α⋅θ₁ + α⋅θ₂⎥\n",
       "⎢ 2                ⎥\n",
       "⎢──────────────────⎥\n",
       "⎢        θ₁        ⎥\n",
       "⎢                  ⎥\n",
       "⎢   2              ⎥\n",
       "⎢d⋅υ               ⎥\n",
       "⎢──── + α⋅θ₁ - α⋅θ₂⎥\n",
       "⎢ 2                ⎥\n",
       "⎢──────────────────⎥\n",
       "⎣        θ₂        ⎦"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "R = R_α + R_d\n",
    "(J @ dE + R @ dS).applyfunc(sympy.simplify)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Port-Hamiltonian reformulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "θ_0 = sympy.Symbol('θ_0', real=True, positive=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$\\left[\\begin{matrix}A \\left(- π_{1} + π_{2}\\right)\\\\υ\\\\- θ_{0} + θ_{1}\\\\- θ_{0} + θ_{2}\\end{matrix}\\right]$$"
      ],
      "text/plain": [
       "⎡A⋅(-π₁ + π₂)⎤\n",
       "⎢            ⎥\n",
       "⎢     υ      ⎥\n",
       "⎢            ⎥\n",
       "⎢  -θ₀ + θ₁  ⎥\n",
       "⎢            ⎥\n",
       "⎣  -θ₀ + θ₂  ⎦"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dH = dE - θ_0 * dS\n",
    "dH"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$\\left[\\begin{matrix}υ\\\\- A \\left(- π_{1} + π_{2}\\right)\\\\0\\\\0\\end{matrix}\\right]$$"
      ],
      "text/plain": [
       "⎡      υ      ⎤\n",
       "⎢             ⎥\n",
       "⎢-A⋅(-π₁ + π₂)⎥\n",
       "⎢             ⎥\n",
       "⎢      0      ⎥\n",
       "⎢             ⎥\n",
       "⎣      0      ⎦"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "J @ dH"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$\\left[\\begin{matrix}0\\\\0\\\\- α + \\frac{α θ_{2}}{θ_{1}}\\\\\\frac{α \\left(θ_{1} - θ_{2}\\right)}{θ_{2}}\\end{matrix}\\right]$$"
      ],
      "text/plain": [
       "⎡     0     ⎤\n",
       "⎢           ⎥\n",
       "⎢     0     ⎥\n",
       "⎢           ⎥\n",
       "⎢      α⋅θ₂ ⎥\n",
       "⎢ -α + ──── ⎥\n",
       "⎢       θ₁  ⎥\n",
       "⎢           ⎥\n",
       "⎢α⋅(θ₁ - θ₂)⎥\n",
       "⎢───────────⎥\n",
       "⎣     θ₂    ⎦"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "C = sympy.eye(4)\n",
    "D_α = R_α / θ_0\n",
    "y_α = C.T @ dH\n",
    "u_α = -D_α @ y_α\n",
    "u_α = u_α.applyfunc(sympy.simplify)\n",
    "u_α"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$\\left[\\begin{matrix}0\\\\- d υ\\\\\\frac{d υ^{2}}{2 θ_{1}}\\\\\\frac{d υ^{2}}{2 θ_{2}}\\end{matrix}\\right]$$"
      ],
      "text/plain": [
       "⎡ 0  ⎤\n",
       "⎢    ⎥\n",
       "⎢-d⋅υ⎥\n",
       "⎢    ⎥\n",
       "⎢   2⎥\n",
       "⎢d⋅υ ⎥\n",
       "⎢────⎥\n",
       "⎢2⋅θ₁⎥\n",
       "⎢    ⎥\n",
       "⎢   2⎥\n",
       "⎢d⋅υ ⎥\n",
       "⎢────⎥\n",
       "⎣2⋅θ₂⎦"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "C = sympy.eye(4)\n",
    "D_d = R_d / θ_0\n",
    "y_d = C.T @ dH\n",
    "u_d = -D_d @ y_d\n",
    "u_d = u_d.applyfunc(sympy.simplify)\n",
    "u_d"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "D = D_α + D_d\n",
    "assert (θ_0 * C @ D @ C.T - R).applyfunc(sympy.simplify).is_zero"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$\\left[\\begin{matrix}υ\\\\- A \\left(- π_{1} + π_{2}\\right) - d υ\\\\\\frac{d υ^{2}}{2 θ_{1}} - α + \\frac{α θ_{2}}{θ_{1}}\\\\\\frac{d υ^{2}}{2 θ_{2}} + \\frac{α \\left(θ_{1} - θ_{2}\\right)}{θ_{2}}\\end{matrix}\\right]$$"
      ],
      "text/plain": [
       "⎡         υ         ⎤\n",
       "⎢                   ⎥\n",
       "⎢-A⋅(-π₁ + π₂) - d⋅υ⎥\n",
       "⎢                   ⎥\n",
       "⎢     2             ⎥\n",
       "⎢  d⋅υ        α⋅θ₂  ⎥\n",
       "⎢  ──── - α + ────  ⎥\n",
       "⎢  2⋅θ₁        θ₁   ⎥\n",
       "⎢                   ⎥\n",
       "⎢   2               ⎥\n",
       "⎢d⋅υ    α⋅(θ₁ - θ₂) ⎥\n",
       "⎢──── + ─────────── ⎥\n",
       "⎣2⋅θ₂        θ₂     ⎦"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(J @ dH + C @ u_α + C @ u_d)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$\\left ( \\left[\\begin{matrix}0\\\\0\\\\- α + \\frac{α θ_{2}}{θ_{1}}\\\\\\frac{α \\left(θ_{1} - θ_{2}\\right)}{θ_{2}}\\end{matrix}\\right], \\quad \\left[\\begin{matrix}A \\left(- π_{1} + π_{2}\\right)\\\\υ\\\\- θ_{0} + θ_{1}\\\\- θ_{0} + θ_{2}\\end{matrix}\\right], \\quad \\frac{α \\left(θ_{1} - θ_{2}\\right) \\left(- θ_{1} \\left(θ_{0} - θ_{2}\\right) + θ_{2} \\left(θ_{0} - θ_{1}\\right)\\right)}{θ_{1} θ_{2}}\\right )$$"
      ],
      "text/plain": [
       "⎛⎡     0     ⎤                                                            ⎞\n",
       "⎜⎢           ⎥                                                            ⎟\n",
       "⎜⎢     0     ⎥  ⎡A⋅(-π₁ + π₂)⎤                                            ⎟\n",
       "⎜⎢           ⎥  ⎢            ⎥                                            ⎟\n",
       "⎜⎢      α⋅θ₂ ⎥  ⎢     υ      ⎥  α⋅(θ₁ - θ₂)⋅(-θ₁⋅(θ₀ - θ₂) + θ₂⋅(θ₀ - θ₁))⎟\n",
       "⎜⎢ -α + ──── ⎥, ⎢            ⎥, ──────────────────────────────────────────⎟\n",
       "⎜⎢       θ₁  ⎥  ⎢  -θ₀ + θ₁  ⎥                    θ₁⋅θ₂                   ⎟\n",
       "⎜⎢           ⎥  ⎢            ⎥                                            ⎟\n",
       "⎜⎢α⋅(θ₁ - θ₂)⎥  ⎣  -θ₀ + θ₂  ⎦                                            ⎟\n",
       "⎜⎢───────────⎥                                                            ⎟\n",
       "⎝⎣     θ₂    ⎦                                                            ⎠"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "u_α, y_α, (y_α.T @ u_α).applyfunc(sympy.simplify)[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "y_{\\alpha}^\\mathrm{T} \\, u_{\\alpha}\n",
    "\\: = \\:\n",
    "- \\theta_0 \\:\n",
    "\\alpha \\, \\frac{{\\left( \\theta_1 - \\theta_2 \\right)}^2}{\\theta_1 \\, \\theta_2}\n",
    "\\: = \\:\n",
    "- \\theta_0 \\:\n",
    "\\alpha \\, \\left( \\frac{\\theta_1 - \\theta_2}{\\theta_2} + \\frac{\\theta_2 - \\theta_1}{\\theta_1} \\right)\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$\\left ( \\left[\\begin{matrix}0\\\\- d υ\\\\\\frac{d υ^{2}}{2 θ_{1}}\\\\\\frac{d υ^{2}}{2 θ_{2}}\\end{matrix}\\right], \\quad \\left[\\begin{matrix}A \\left(- π_{1} + π_{2}\\right)\\\\υ\\\\- θ_{0} + θ_{1}\\\\- θ_{0} + θ_{2}\\end{matrix}\\right], \\quad - \\frac{d θ_{0} υ^{2} \\left(θ_{1} + θ_{2}\\right)}{2 θ_{1} θ_{2}}\\right )$$"
      ],
      "text/plain": [
       "⎛⎡ 0  ⎤                                     ⎞\n",
       "⎜⎢    ⎥                                     ⎟\n",
       "⎜⎢-d⋅υ⎥                                     ⎟\n",
       "⎜⎢    ⎥  ⎡A⋅(-π₁ + π₂)⎤                     ⎟\n",
       "⎜⎢   2⎥  ⎢            ⎥         2           ⎟\n",
       "⎜⎢d⋅υ ⎥  ⎢     υ      ⎥  -d⋅θ₀⋅υ ⋅(θ₁ + θ₂) ⎟\n",
       "⎜⎢────⎥, ⎢            ⎥, ───────────────────⎟\n",
       "⎜⎢2⋅θ₁⎥  ⎢  -θ₀ + θ₁  ⎥        2⋅θ₁⋅θ₂      ⎟\n",
       "⎜⎢    ⎥  ⎢            ⎥                     ⎟\n",
       "⎜⎢   2⎥  ⎣  -θ₀ + θ₂  ⎦                     ⎟\n",
       "⎜⎢d⋅υ ⎥                                     ⎟\n",
       "⎜⎢────⎥                                     ⎟\n",
       "⎝⎣2⋅θ₂⎦                                     ⎠"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "u_d, y_d, (y_d.T @ u_d).applyfunc(sympy.simplify)[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "y_{d}^\\mathrm{T} \\, u_{d}\n",
    "\\: = \\:\n",
    "- \\theta_0 \\:\n",
    "\\frac{d \\, \\upsilon^2}{2} \\, \\frac{\\theta_1 + \\theta_2}{\\theta_1 \\, \\theta_2}\n",
    "\\: = \\:\n",
    "- \\theta_0 \\:\n",
    "\\frac{d \\, \\upsilon^2}{2} \\, \\left( \\frac{1}{\\theta_1} + \\frac{1}{\\theta_2} \\right)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Port-Hamiltonian reformulation 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "f_6 = θ_1 - θ_0\n",
    "f_9 = θ_2 - θ_0\n",
    "e_6 = 1/(θ_0 + f_6) * α * (f_9 - f_6)\n",
    "e_9 = 1/(θ_0 + f_9) * α * (f_6 - f_9)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAe0AAAAXCAYAAAAmyhcwAAAABHNCSVQICAgIfAhkiAAAB8RJREFUeJztnVmMFUUUhj9AccYZo9EHXAIuMRoIMRiFuEQcCcaoUYPREHmxXWIE1xh9kAc1PmgMElEDUXxwNG5EXOISDBqNWyIuEeMuEgdx33AHUcGHU51p+1Zv1VXdl+F8yc2dqequ+uecOufW7erqAUVRFEVRtjvuBb4D+gK13wtcA3wEbALWAzcAOwbqzxc+dR8GbAXOK3FskT/UntXs2YSebQGNczvdGuc+2N594iNP+NTjhcOBLcAVgdrfC3gP2AwsA24C3kcMeWegPn0QQvdjwNdAf84xRf5Qew5Txp4A44B/gdsC6+lmNM7tdGuc+0B9ItTJEyH01GYl8DMyk/DNWOBN4Hfg6ER5P7AOMdCeAfqtSyjd0xBHz885Js8fas//U8aeABeY444LrKeb0TjvpFvj3Afqk2Fc80QoPbU4CJntLQ3U/nzECBdZ6haZujMD9R2Z9gcczg2p+0PE2WMsdUX+aNOedWjLnjErgB8Sx2yr49IVjXM73RrnPmjLJxH1xnco3S55wrue0ZayHYHLgdXARuAL4BZktrAz8C1wf+L4c4FRyFd+3/QCVyGXJWyD80fz3m2zvdC6HwImADMtdXn+UHvaybMnwK7ADOBJZFa8rdoxicZ5fbo1zn2gPumkap4Ioif9ob078CoSvOuR6/K/IMF9KXChOebaxDkzjcDXynZagVnAbsADwN+W+h7zvjlA33UIrftV8368pS7PH2pPO3n2BDgZ+TB7tCE9odE490O3xrkP1CedVM0TQfTskPr9IWAqcBnDC+kLkFn4icAkYBD41NT1AVOQywZ/lO20Aieb932A6yz18YxnfYC+6xBa9xvmfXqqvMgfak87WfaMmYXY89mG9IRG49wP3RrnPlCfdFI1T4TWw0zk2vpLyGWXJB+Zuk3A+ET5QaZ8pUuHJVhn2i967Reo/wi3tZUmdG8EvkmVFfmjbXu60pY9QWbCvwEPN6wnjwj3NT+N804iRlac+6BNn0S4j+/QuqvkieB67jEnn26pW23qbk2VH2nKs9ZVhkoKjl/3Jc7tM2XvZrTdj1xSSM9Q5gGfIYnnLeCYjPPrah3MaMdV99XITO5X4HtkXWRyjt4vgX9SZXn+cNE13ej4ypx7Ro6eJEO4+92H7qq2BLs9AU41/c+pocdVE/gblzEa5+3GedWYqhrnMUM065OmcsVgTltN5IqyecJVT6Edk5fHj0WuuT+T0cGfyEbwJBvNew921iJBVZavEj/vYylLcgJyM83TibLZSMKZB7wCzEXu5psEfF7Q9yJk7SHJFOA0JNENpepWZ7TjohtkVrkEGUCjgOuB5xDtP1na6WXY/jF5/nDR1Qe8A9wNPJJxno06fk/jonuAarYEuz1BLnltTrTflH9jfI3LGI3zduO8akxVjfOYpn3imit8ju8mckXZPOGqp7Qde5HtA2ssdQcgzn/eUrc3Mht4Ja9xRyabtp/IqH/K1E9LlK0C7kodtwa40VFDZPoYqHCOi24b/ciNJqdY6kYj/lqbKs/zR11dVWbPPvFhzzxbQrY9xyDbN1Z41lNGUx4R1cclaJxnEdFOnBfFlEuc+6DtXBHhNr5D54oqecKHHqsd47vHe5FZxhbLibcAO2G/JPA1cknh4IxO6xCvG9huhT8COAkx0uumbCzyuLn0Os9K4KgA+rKoqjuLXRD/2GZ7ByP+Ss868/zhS1fT+NCdZ0vItud0YA/kaUg+9ZTRFAKNc380EU8uce4DzRXV8q4tT/jSk8koZBH9X+CQRPlchtcSshpebuoPdOm4gA8smvZFLpls4P+L9/HsM31n3zXAx479R7jN+KrozmIZ8Db2jfznGF0XW+ry/FFHV1vftKG+PfNsCdn2vN30O86znjKa8ohwG5ca53Yi2onzophyjXMftJkrItz8AWFzRdU8UVdPoR1vNwd9ByxGrqdvAR4HXjB1dyBbRZKcZepsT3upyxzT9vfIN4GlyGP7NtA5q46DOX1DyrXIXbEuRLgNniq6bSxAZtNZAfkg8o1ovKUuzx91dLX5oV1Hd5EtIdue64GXPespqymPCPekpnHeSUQ7cV4UU65x7oM2c0WE+/gOmSuq5om6egrt2APcbAT8jQT1QmShfBrD20HSG8vHIpcBVhUIcOVs5MHqm5B9pEsZXuBP6/iHzsfBLQZedOw7wn3wlNWdZiHyNKpJGfW7IjdCPJ5RX+QPV11tfmiDm+4iW0K2Pacif3PWP2MI5d8yRLiPS43zTiKaj3PIj6m6ce6DtnJFhLs/IEyucM0TrnogcM692nRwaKgOSrKKzkfFfYL7DSpNcxvFCf0S7N80koTwR9sf2lUpY0vItucNpnz/FjR1KxrnfsmLqbbi3AcjMVc0mSdigtqxB9lE/mSoDkoyG7nt/nxgIrKd4HdkDaHbWYLsFZyB3LgQv5L/Cq4X2T6wvKAtX/7oR7ZfTEEG0JXm5wk12w1NGVtCvj0/pHj7VAhN3YzGeX3KxFTTce6DkZwrmswTjdpxOrKuFPKfsZdhHrLI/xfy0IWsR851G1szXtcljploft+vRHs+/DGQoWmwRptNUMaWUM2eTWnqdjTO6zFAcUw1Hec+GGDk5oom88RAhp7BBvpWFEVRFEVRFEVRFEVRFEVRFEVRFEVRlFr8B7HPebaSj56DAAAAAElFTkSuQmCC\n",
      "text/latex": [
       "$$\\frac{α \\left(- θ_{0} + θ_{2}\\right) \\left(θ_{1} - θ_{2}\\right)}{θ_{2}} + \\frac{α \\left(- θ_{0} + θ_{1}\\right) \\left(- θ_{1} + θ_{2}\\right)}{θ_{1}}$$"
      ],
      "text/plain": [
       "α⋅(-θ₀ + θ₂)⋅(θ₁ - θ₂)   α⋅(-θ₀ + θ₁)⋅(-θ₁ + θ₂)\n",
       "────────────────────── + ───────────────────────\n",
       "          θ₂                        θ₁          "
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(e_6*f_6 + e_9*f_9)\n",
    "# this is always negative\n",
    "# net outflow is negative == net inflow is positive == exergy destruction rate is postive"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "f_4 = υ\n",
    "f_7 = θ_1 - θ_0\n",
    "f_10 = θ_2 - θ_0\n",
    "e_4 = -d * f_4\n",
    "e_7 = 1/(θ_0 + f_7) * sympy.Rational(1,2) * d * f_4**2\n",
    "e_10 = 1/(θ_0 + f_10) * sympy.Rational(1,2) * d * f_4**2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAcUAAAAaCAYAAADG3kfKAAAABHNCSVQICAgIfAhkiAAACLlJREFUeJztnXmoVUUcxz+pyROFotWSVsKyDVtptbJFQlpsQYqoEy1EO9GCUvYKaY/MUsigXtmCtFjZLkGLgqZlVGo7mu22WBmaLa8/fnPwNO+sc2fO8u584HKfM3Nmfr+Z850zd2bOCB6Px+PxeBrJOGAB8BuwEpgF7F6pRW5pN389niJ4fXjanleAs5Ebfw9gJvAdsEmVRjmk3fz1eIrg9eFpPNsC3cDTlvIbBPwDHGcpP5vY9hWy/X0Y+AEYaLFMnQHABOAjYC2wArgJ2NBhmTawafc+SNueY826fDS17pOooj9wrZGmtpHXR0WciFTWtZby20rld7Cl/Gxi21dI93df4F/gCovlxZX/IbAOmAHcCixWNt3nsNxWcWH3TOBbpCNOYkukk56s/r0pcK669jNgDfArMAfpQPqU7EPVlN0fuNZIU9uoKn2APY00te65ETHyWEv5zQAWAX0t5WcT275Cur+vAquQ0ZIL+gMLgdX8v9MZBCxHbuzBjspuBVd274+07/iUNOerNEeof1+g/v0N8ChwM/AA0m7dwJPABiX6UDVl9wcuNdLUNqpSH2BHI02tewBeQBzb0kJetyMjkZ1ayCNQ9hxuwR4dm75Cur9DkRHwNEtlxTEe8eeimLhJKu5UB+UGtNZGLu1eioguqRN+CfgxEj8SmdrTR7uDgS+VLSfH5NPUus+izP7AtUaqaiNorZ2q1AfY0UiVdZ+LfsBlwPvIT9/lwNXI0/1b4Cst/RTE6JFa+HTkJt4tpow7ge+BXVu0NaA10ZfhK2T7e4vK98hi5udmAPALMnqLm5+/TpV/iYOyA8zbyLXd16vrR8XEbQT8CTyYM69Q2Pdo4U2t+5C69AcuNVJlG4F5O1WpD7CjEes+pK1hmNAfeBF5Ov8D3Au8BtyAjNAGA+9q1+ylvt/TwqchwhmrhU8GzkB+bi+xZbgBZfgK+fw9Stkwr5AH+RkDbAw8BvwVE9+hvtc5Kt8U13bPVd9Hx8SNRu6RvJtIQvv+1sKbWvdQr/7ApUaa2kZV6gPsaMS6D/3yJszJFKQCJgATkSc0QBfwhvo7KoI+wJ7I6PFnLa/F6js6MpyKCOBElT6cJ16tPmXi2lfI5+9AYDgyVfGHkSfZjFbfQ4DOmPij1PcKR+Wb4truBep7REzcGKQ9ZufIpx9wpvr7ZS2uqXUP9ekPXGukqW1UpT7AjkZqXffhwuobCfFLVHx0u/QuKmxmTPrNVdyzkbDuhE+noc0BZtMOZfgK+fwdqsJeLeJAQZan2BL9bO+g7ADzKbwy7F6DvBsXpQP4HXgiZx53KDteiIlrat3XqT9wrZEq2wjM26kqfYA9jTj3YVnOAsLPI5Frp6uwwxLyfkvFD4mEnabCrotJf4iKu9fIk54so5hvXSl51cnXA9W1M1LSLMO8XQeqsA8S8h6ETE3YGIkVtbMrJa+y7P6anlOex6uyT89x/aUq7VJ6vnTe1LqHZmlkGc3Qh4mtXQn5VKkPsKMRJz7o06efIy895uWbyN/HAD8Bbyak3RFZDP86EhauHyyKSR+OIOfGxJkwCZl7jjIcOAF4CLnZouhrGlHq5Osa9d2RkqaVdh0SExZlFLLArY/gLgSuQt4fWgxcjnSEadhsIxO7xwEnATsjGwDmqbAPU8oZwPo2CBmDiDHul1+Ui4C7kV9NR9JzyrAsH8Bu3UOzNNIUfYC9djK1ewRi9z6I7acir0kkEacPsKORsnwwogN5YuuL5iEHqfgXtfDZ9Bwtgiy+fom8m+LydJaA4tMOdfN1a5XvHINr87C7yv+5hPjnVfz+kbCxyKL3ecAwZDPEauQEk6IEFG8jMLO76LFhfZAdkZ9HwvoiW8xfyrDvctaPcreo0Ic0Aszqvp00UrU+wKydTOwGead0IjLw6gZOSSkjTh9gTyNl+GDMhshPZH17dUh4s0/Uwr9C5pV1rlDpx9kyMIGA4jdT3XzdADm2aqXh9Vlshtj3dkzcAchNr3du84H7tbBPkRdxixJg1jGb2K2TdWzYMFXGU5GwI1TY+Sn5XqPSLFJ2JlGGD2kEmNV9O2mkan2AWTvZuLeyHihx+gB7GinDh5ZYqgrQxRc61438ZI6yEjnhPjqNOwqZypiP+zPrAsxEXzdfn1RltnKQQRpLkI51z0jYdshUzS/8fxG7P9Ih6i/LTiF500UaAWZtBMXsjiPr2LCzVfzFkbB7VJlJL6SH700tJN+vN9c+pBFgXvftpJEq9QHm7dTqvZX1QInTB9jViGsfWuJMVcA65NDd25An+BrWn0SwnXbNLBX+JrK7aBbi4PskTynZJMDsZqqbr+EGhbgTHWxwusp/JXAX8s7YKuSmO0hLG05V6duwJwAfG5QdYN4xF7E7jqxjwx5HOrhtImErSF4bOkvZ87eypzPmE5TsQxoB5nXfThqpUh9g3k6t3ltZD5Q4fYBdjbj2oWUuA75A5su/R342741sm/0xJv0OyDbpP9TnPWQkmbZpxCYB5qKvk6/9kXWj+RbySuIsZEPAWmSaaxo9135gvegP1cKvR06vL0qAeRtBfrt1so4N2wjp4J+JhO2nbE06cLqT7N2Cr5foQxYBrdV9O2mkKn1Aa+1kem9B+gMlTh/gRiOufPA0nHFIA++VldAxLqaHyibPsWGX0LNzu0mF7eDOtNzYOgqxN1EHjfQGfYSkPVDi9AH10gj4h2KvpgMZgc+q2hBkNK4fvPwJ5hsJymQy2Q+TAcg2cH0b91KyX1Uogzw+tCN10UiT9REl6YGSpA+oj0ZC/EOxlzMCmYZx+RpLHsYia0nnIjvQJiFbzvU1pLoxFdnoMRI5Miz86P8n3DBkmmf7Em3LS14f2pU6aKSp+gC5j4arTzdwpfo7+jpJnfUB+XzweKxzIbID7E/gHZLPP6wTSWsYnRXaVJTe4EM70ER9gKxbxt1fXdWZVJjDab4PHo/H4/F4PB6Px+PxeDwej8fj8Xg8Ho+nIv4D/PI4FjG458kAAAAASUVORK5CYII=\n",
      "text/latex": [
       "$$- d υ^{2} + \\frac{d υ^{2} \\left(- θ_{0} + θ_{2}\\right)}{2 θ_{2}} + \\frac{d υ^{2} \\left(- θ_{0} + θ_{1}\\right)}{2 θ_{1}}$$"
      ],
      "text/plain": [
       "            2                 2           \n",
       "     2   d⋅υ ⋅(-θ₀ + θ₂)   d⋅υ ⋅(-θ₀ + θ₁)\n",
       "- d⋅υ  + ─────────────── + ───────────────\n",
       "               2⋅θ₂              2⋅θ₁     "
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(e_4*f_4 + e_7*f_7 + e_10*f_10)"
   ]
  }
 ],
 "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}