notebook_diffusion_mixed_formulation.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 70,
   "metadata": {},
   "outputs": [],
   "source": [
    "from nutils import *\n",
    "import numpy, unittest\n",
    "\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib import collections"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "metadata": {},
   "outputs": [],
   "source": [
    "nelems = 4\n",
    "degree = 2\n",
    "tol  = 1e-5\n",
    "ndims = 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "metadata": {},
   "outputs": [],
   "source": [
    "domain, geom = mesh.rectilinear([numpy.linspace(0,1,nelems+1)]*ndims)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {},
   "outputs": [],
   "source": [
    "# prepare residual\n",
    "ns = function.Namespace()\n",
    "ns.pi = numpy.pi\n",
    "ns.x = geom\n",
    "\n",
    "# construct bases\n",
    "ns.alpha, ns.beta = function.chain([\n",
    "    domain.basis( 'std', degree=degree+1),\n",
    "    domain.basis( 'discont', degree=degree),\n",
    "])\n",
    "\n",
    "# solution variables\n",
    "ns.phi = 'alpha_n ?lhs_n'\n",
    "ns.u   = 'beta_n ?lhs_n'\n",
    "\n",
    "# analytical solution\n",
    "ns.g  = 'cos( 2 pi x_0 )'\n",
    "ns.dg = '-(2 pi) sin( 2 pi x_0 )'\n",
    "ns.f  = '(2 pi)^2 cos( 2 pi x_0 )'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "metadata": {},
   "outputs": [],
   "source": [
    "# construct residual and inertia vector\n",
    "res = domain.integral('- alpha_m,0 u + alpha_m phi' @ ns, geometry=ns.x, degree=degree*2)\n",
    "\n",
    "# construct residual and inertia vector\n",
    "res += domain.integral('- beta_n phi_,0 + beta_n f' @ ns, geometry=ns.x, degree=degree*2)\n",
    "\n",
    "# boundary integral\n",
    "res += domain.boundary.integral('alpha_m g n_0' @ ns, geometry=ns.x, degree=degree*2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "solve > solver returned with residual 1e-15\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x1517795dd8>]"
      ]
     },
     "execution_count": 75,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAD8CAYAAABzTgP2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3Xd4U3X///HnO+miBdpCy7BlCggICFiWCOJgIwVlO3DiAge4bvV2gHuPr4qIA0UFRLmZgoCAKEPKKDIEyi4ttGwKncnn90fD/SvchRaS9iTp+3FduZqcnJPzOoy8erYYY1BKKaVOs1kdQCmllHfRYlBKKXUGLQallFJn0GJQSil1Bi0GpZRSZ9BiUEopdQYtBqWUUmfQYlBKKXUGLQallFJnCLA6wMWIiooytWvXtjqGUkr5lNWrVx80xkQXNZ5PFkPt2rVJSEiwOoZSSvkUEdldnPF0U5JSSqkzaDEopZQ6gxaDUkqpM2gxKKWUOoMWg1JKqTNoMSillDqDFoNSSqkzeKQYRORLEUkTkQ3neF9E5EMRSRKR9SLSssB7Q0Vkm+sx1BN5zuU/a/cxJWEvmTmOkpyNUkp53IZ9xxg9cxO5DmeJz8tTJ7h9Dfwf8M053u8O1Hc92gCfAm1EpBLwAhAHGGC1iMwwxhzxUK4zzEhM4bd/0nhl9mb6XxnLLW1rUScqrCRmpZRSbsvKdTB7fSrfrtjNur1HKRdo56aWMTSJCS/R+XqkGIwxv4tI7fOMEg98Y4wxwAoRiRCR6kAnYL4x5jCAiMwHugE/eCLX2b5os5/tsWn8sLcS3yzbwfg/dtKhfhS3tq3F9Q2rEGDXLWtKKevtPnSS71fuYUrCXo6cyqVFlIPPrsrgmoqphFTpVOLzL61LYsQAewu8TnYNO9fw/yEiw4BhADVr1ryoEPLXOOrt/J1/A8+FliM1tD5/7othwfc1mBh6GW1bt6d/2zpUqRByUZ+vlFIXy+E0LPonjRl/riFj12qayi6+Ck+hUbmdBGekwBrXiJd1hurNSjRLaRWDFDLMnGf4/w40ZhwwDiAuLq7QcYp06zQ4uBVSE5HURC5JTaRf1lL6OzMgF7L/CGDL0ppsiWxC3V6jiKnf/KJmo5RSxZWXnUnCT+/gTPqNZo7t3CBHIRAMgoTUg+pXQfUrXI9mUC6yxDOVVjEkAzUKvI4FUlzDO501fHGJpbAHQNXG+Y/mgwEQpxMO74DUdWTuSKBC0irqHZ1H0MRZbKo1mIYDX8EWVvJ/EUqpMsYYUldNg7nP0NaZyt6AmuTWuAbHZW2xx7RAqjWB4AqWRCutYpgBDBeRSeTvfD5mjEkVkXnAqyJy+pu3C/CvUsqUz2aDqHoQVY+Ipv2IAA6k7mXDd0/Taff3ZLwzi7xrnqFSh3vBZi/VaEop/+Q4sJnUyY8Re3g524lhWbtxtOsyAJHCNqKUPsnfH+zmh4j8QP5v/lHAAfKPNAoEMMaMlfyl/T/ydyyfAu40xiS4pr0LeMb1Ua8YY74qan5xcXGmpC+7bYxh7sL5VF76Aq1lE0cqNCDipneQOh1LdL5KKT+WeYTjc8cQlvgVGSaEOZXv5Prb/kWVyNJZMxCR1caYuCLH80QxlLbSKIbTkg+fZOrET+h3aCyxcpDM+r0o1+NViKxVKvNXSvkBRx7O1V+TM38MgTnHmCo3UK7L89zYrmmpriVoMXiQ02n4YdkWDv36Dvfa/kOQDWztH0Y6jIQgPQ9CKXUeO38nZ9aTBB3azApnI2Ze8gjDB/eheni5Uo+ixVACdh86yeuTFtJ1/6f0sS8jL/JSAoZOh4gaRU+slCpbjMH89jKy9G32mWje5jba9LiDga1rWrYvobjFoGd0XYBalcP4+IEbOdjlY4Y6niPzSCq547vAwSSroymlvInTiXPOE8jSt5mc14lnYr5g1KNPMKhNLa/ZwXw+WgwXyGYT7ulQl389dB8PBLzEiYwMcsd3hf2FXiZKKVXWOPJwTHsA26rPGZfXk+QOb/DVPR2JjQy1OlmxaTFcpIbVKvLqA7cyIuQ1DmYa8r7sDnv/sjqWUspKednkTb4N+9+TeDu3P7YuYxjVtSE2m/evJRSkxeCGmpVDee/BfjxZ8U2Ss0PJ+zoeti+yOpZSygo5J8n9tj8BW+fwUt7t1OzzAvd0vNTqVBdFi8FNVSqG8NEDvXkx6m2Scivj+K4//DPb6lhKqdKUeZTcr3pj272UJ/Pup/XAZxjQyncPStFi8ICI0CA+Htadd2LeZX1eLZyTb4PESVbHUkqVhow0cr7oAanreMzxKL1uf5zuTatbncotWgweEhYcwEd3Xc8Xdd9neV5DmHYf5q/PrY6llCpJR/eSM74rzoPbeIinGHrPCDo2iLY6ldu0GDwoJNDO+7e1Z/rl7zHfcSUy53HM7+9YHUspVRIOJpHzeReyj+7nIfvzPDrsfq6sVcnqVB6hxeBhAXYbrw9ozfK495jmaI/8Nhrn8k+tjqWU8qQTB8j5sicnMjIYETSG5x64i8aXVLQ6lcdoMZQAm034d+9m7OzwLr86rsTMexZ2/Wl1LKWUJzhyOfX9rThOHuapsNG8+uAtfneLYC2GEiIijOzSkMQrX2eXswqZ398Gx1OsjqWUclPm7GcITf2LVwMe5OX7BnNJROlf86ikaTGUsJG9W/FlzBgc2Sc58c0QyMu2OpJS6iLlrv2BcmvG8bWzO/3vfIxq4f55G2AthhJmtwlPD+3Du2GPUuHgWo5PG2V1JKXURTCpiZgZj7DC2YgqN71Js9gIqyOVGI8Ug4h0E5EtIpIkIk8X8v57IrLO9dgqIkcLvOco8N4MT+TxNhVCArnznkf4SvpQceO3nFrxtdWRlFIX4tRhTkwYzCFnGH+3e58ezWtanahEuX1rTxGxAx8Dncm/h/MqEZlhjNl0ehxjzGMFxh8BtCjwEZnGmObu5vB2NSqF0uS2t/jz6yRazX2cvJgmBNQo8uq3SimrOR0cnHAbFTIPMK7Wh4zq1sbqRCXOE2sMrYEkY8wOY0wOMAmIP8/4g4EfPDBfn9OqbhUOdvuUNBNOxjdD4ORBqyMppYqQPvMFog78wecVHmT4bYN84rLZ7vJEMcQAewu8TnYN+x8iUguoA/xWYHCIiCSIyAoR6eOBPF4t/qpm/Nr0bUJyDpP6xWBw5FkdSSl1DkfX/Ez02o+YbruBAcOeJSTQbnWkUuGJYiisPs91W7hBwFRjjKPAsJquOwoNAd4XkUIvRygiw1wFkpCenu5eYosNvakP30c/RvXDf7F36lNWx1FKFSI7dTNBMx9kvbmUekM/oUpF/zwCqTCeKIZkoOBlBGOBcx2wP4izNiMZY1JcP3cAizlz/0PB8cYZY+KMMXHR0b59LRK7TRhw71PMCOpJjc3j2f/nd1ZHUkoVYLKOceSrAZxyBnK45xdcXquq1ZFKlSeKYRVQX0TqiEgQ+V/+/3N0kYhcBkQCywsMixSRYNfzKKA9sOnsaf1R+eAAWg77hEQuI3z+YxzblWh1JKUUgDHs+HwoUdnJ/HHFm3RqXejvqn7N7WIwxuQBw4F5wGZgijFmo4iMFpHeBUYdDEwyxhTczNQISBCRRGAR8HrBo5n8XWxUBDLwG06YcpycOARnTqbVkZQq87bMeJtLDy1iVtX7ie87yOo4lpAzv6d9Q1xcnElISLA6hscsmj2Ja1fdx7rad9H8jvesjqNUmXV43zbKfd6eDQFNaPrkr4QEuX1Ev1cRkdWufbrnpWc+e4FOPQbyR/muNNn5NXs2LLM6jlJlknE62T9xGA4jRA782O9K4UJoMXgBEaHh0I84IhXJnfYQuTl6PSWlStvq6R/TOHMNay97lHr1G1kdx1JaDF4iKroqe9uO5lLHDv767kWr4yhVpuzft4sG615jc+DlXDXwSavjWE6LwYu07DaUxArXELfrc7Zs8J99KEp5M6fTsGfigwSTQ/jAsdjtZeMktvPRYvAydW7/hCwJJmfacLJycq2Oo5TfW/Kf8bTO/JPNlz3IJfWaWR3HK2gxeJmK0bHsb/c8TR2bWTTxNavjKOXXdu7dS5PEMewKqs8VA/9tdRyvocXghS7rMoykCq3psPtjVifqiW9KlYQ8h5MdEx8lUk5QccBYxB5odSSvocXgjUSIue0z7AI50x/hRGaO1YmU8juzpn3H9dkL2HnZvVSqp5fAL0iLwUuVq1KXg22fpp1zLbMmvm91HKX8yqZdKcT9/SIHgmpQv99oq+N4HS0GL1ajy8Psq9CMbskfsGTNRqvjKOUXsnIdbPn+SS6RQ4T2+xQCy85VU4tLi8Gb2exE3zKO8pJF9szHOZShJ74p5a7JP08lPnsWKfVvpUKDDlbH8UpaDF4uqFojjrZ6jC5mGZMnfoovXttKKW/x17YU2m98kWNBVYjtp0f9nYsWgw+I7vYUh8Lqc3Pqe8xN2GJ1HKV8Ulaug82Tn6eeLYVyfT+C4ApWR/JaWgy+wB5IxODPiZLjHJ87mlM5ejtQpS7Uj78uZnDuz6TV6UNI465Wx/FqWgw+wh7bgkMNBnJT3lwm//Jb0RMopf4r9Vgm1f96FactiCo3vWF1HK/nkWIQkW4iskVEkkTk6ULev0NE0kVknetxT4H3horINtdjqCfy+KsqvUfjsAdTa/Vr7D18yuo4SvmMn6Z+zw2SQFbbR6BCNavjeD23i0FE7MDHQHegMTBYRBoXMupkY0xz12O8a9pKwAtAG6A18IKIRLqbyW+Vr0LOVSO5zraGn6d+a3UapXxCwo50rt/9PseCqxNx3aNWx/EJnlhjaA0kGWN2GGNygElAfDGn7QrMN8YcNsYcAeYD3TyQyW9V7PQwx0Ji6Jr8Ecu27bc6jlJezeE0LPvpAxrZ9hDSfQwElrM6kk/wRDHEAHsLvE52DTvbzSKyXkSmikiNC5xWnRYQTLmer9DQtpeEnz4gz+G0OpFSXmva8s0MzpjAoUotCL6in9VxfIYnikEKGXb2wfYzgdrGmGbAAmDCBUybP6LIMBFJEJGE9PT0iw7rD4Ka9OFwVBy3ZH7Lj3/qGdFKFeZYZi4nFrxJtByn0k3vgBT2daMK44liSAZqFHgdC6QUHMEYc8gYc/q03c+BK4s7bYHPGGeMiTPGxEVHR3sgtg8TIbLv20RKBtm/vcGRk3qRPaXO9s2cJQxxzuRIvZuQ2CuLnkD9lyeKYRVQX0TqiEgQMAiYUXAEEale4GVvYLPr+Tygi4hEunY6d3ENU0WQmBacuGwAQ8wcvpq10Oo4SnmVpLQTXLruTcQWQOSNL1sdx+e4XQzGmDxgOPlf6JuBKcaYjSIyWkR6u0Z7WEQ2ikgi8DBwh2vaw8AY8stlFTDaNUwVQ3iv0Rh7EJdveIdNKcetjqOUVzDG8MPUyfSwrySv3QgI192WF0p88do7cXFxJiFB74kMkLnwdcotfY2XKr3B8yPuQ3Q7qirj5m9Mperk7tQtd4ryo9ZBUKjVkbyGiKw2xhR58wk989nHlev4CBkh1el38BPmrN9ndRylLJWd52DV9E9oZttJSLeXtBQukhaDrwssR2iPl7nctpvEmZ+QmeOwOpFSlpmweCN3ZX/L8cpXEHDFQKvj+CwtBj9ga3ozJ6Jbcm/uRL78bb3VcZSyxIHjWeQtfZ9qcoSK8W+BTb/eLpb+yfkDESrEv0W0HMO+7H32Hc20OpFSpe6zGUu4ixmcrB8PNdtYHcenaTH4i9g4TjXsx50ym8+mL7I6jVKlau2eIzTb8j52m42wnnp4qru0GPxIaPfR2Ow2WiV9wMaUY1bHUarUTJs5nT72ZZh2D0FETavj+DwtBn8SHoOj7XButK9g2uzZVqdRqlSs3H6Q7vs/IzOwEkHXjLI6jl/QYvAzIR0fITMgnKv3jGXDPl1rUP5vwezJtLNvIuDaJyG4vNVx/IIWg78JCUc6PEoneyKzZ061Oo1SJWpZUjq9Do4nI6Q6ga3vsjqO39Bi8EMh7e7nZFAU16Z8RuKeI1bHUapEGGP4c9YErrDtIPiGZyAg2OpIfkOLwR8FhWK/9ila27awcNZ3VqdRqkQs25ZG/OEvORZam8AWQ6yO41e0GPxUSKs7OBYSQ5f9n7Nm9yGr4yjlUcYYVs8aRwPbPsp1ex7sAVZH8itaDP4qIIjgzs/RxLaLP2d+aXUapTzqj39S6HN0AocqNiKoSV+r4/gdLQY/FtJiIIfDLqVH2hes3plmdRylPMIYw8bZH1PTlk6Fni/ppS9KgP6J+jObnbBuL3CpLZXVMz61Oo1SHvH7pj30PfEdaZEtCWrQxeo4fskjxSAi3URki4gkicjThbw/UkQ2ich6EVkoIrUKvOcQkXWux4yzp1XuCW7Sm7QKl9Pz8ARWJaVaHUcptxhj2Dn7ParKUSJvHKP3cS4hbheDiNiBj4HuQGNgsIg0Pmu0tUCcMaYZMBV4s8B7mcaY5q5Hb5RniRDeawwxcoiNMz6wOo1Sblmyfjt9Tk4hJfpqAutebXUcv+WJNYbWQJIxZocxJgeYBMQXHMEYs8gYc8r1cgUQ64H5qmIKvux6UiJb0/PYd6z8Z4/VcZS6KMYYUue+RYScJDpeL5RXkjxRDDHA3gKvk13DzuVu4JcCr0NEJEFEVohIHw/kUYWoHD+GaDnOjtnv4Iu3c1VqyZpN9D41jb3VuxIY28LqOH7NE8VQ2Ea+Qr95RORWIA54q8Dgmq57kA4B3heRS88x7TBXgSSkp6e7m7nMCa7dlj3Rneh5fAp/bdphdRylLogxhiO/vkGI5FK9z2ir4/g9TxRDMlCjwOtYIOXskUTkBuBZoLcxJvv0cGNMiuvnDmAxUOivAsaYccaYOGNMXHR0tAdilz1V48dQXjLZN+cNXWtQPmXJqnX0yJrD3hrxBFRtaHUcv+eJYlgF1BeROiISBAwCzji6SERaAJ+RXwppBYZHikiw63kU0B7Y5IFMqhDBsc3YXb073TL+w8q/N1sdR6licToNmQtexSaG2D4vWh2nTHC7GIwxecBwYB6wGZhijNkoIqNF5PRRRm8B5YEfzzostRGQICKJwCLgdWOMFkMJuqTvaIIkl8O/vKprDconLF2xnM7ZC9hTZyABlWtbHadM8MgFRowxc4A5Zw17vsDzG84x3TKgqScyqOIJrlKfbTVu4oY901i5dh1tW+pOPOW9jDE4F71GrgRSu+8LVscpM/TM5zKoVt+XMGIja8HrVkdR6rxW/fUn1+QsZU+927BXrGp1nDJDi6EMCqoUy5aYm7j65K8kbfnb6jhKnZNZ/DqZEkKd3v+yOkqZosVQRtXq/Rx52Dn8y6tWR1GqUDs2rKRN5lI21xpCUMUoq+OUKVoMZVR41ZqsrdKHlkfmcmi3HqGkvE/GvJc5YcpRP17XFkqbFkMZFtvrGfKwkzprjNVRlDrD4aQEmp34ndXVBxFeSc9bKm1aDGVYjVp1WRrRm0bpc8jav9XqOEr916HZozluQql745NWRymTtBjKuKhuT5FjAkmZ8ZLVUZQCIGvPWuofWcKSSv2pGXOJ1XHKJC2GMq55wwbMDe1JrZQ5ONN0rUFZL33WixwzocR0H2l1lDJLi6GMExHKXTuKbBPIAd3XoCzmTF5LjbTFzAq7iRb1axU9gSoRWgyKG668nGkB3am6Zxak61qDss6h2S9x1IRR6bpHEL07m2W0GBQBdhvOdiPINIEcnas3QFEW2bea6NRFTA6Mp3OLelanKdO0GBQA8VdfwQ90o+L2GZD2j9VxVBl0fO7LHDHlCb7qAQLs+tVkJf3TVwBUDAnkaPP7OWWCyVygZ0OrUpacQMW9vzGBG7npqkZWpynztBjUfw28pjnfOLoQsnUGpOnZ0Kr0ZC14mcOmPNkt76ZiSKDVcco8LQb1XzUqhbKrwV2cJIS8316zOo4qK/b+RciuRXzu6MWQDpdbnUbhoWIQkW4iskVEkkTk6ULeDxaRya73V4pI7QLv/cs1fIuIdPVEHnXxBndqzld5XQn4Zzoc2Gh1HFUG5P32KoepwIGGt1OjUqjVcRQeKAYRsQMfA92BxsBgEWl81mh3A0eMMfWA94A3XNM2Jv9WoJcD3YBPXJ+nLNKiZiQJ1YeQQSgHZr7E+wv08FVVMjbsO8ZnE78nYOcixub24raOZ39tKKt4Yo2hNZBkjNlhjMkBJgHxZ40TD0xwPZ8KXC/5BynHA5OMMdnGmJ1AkuvzlIUGdWzGF3ldqZo8j3kLF7A59bjVkZQfemPuPzTe8jGHCGdDzABa1Iy0OpJy8UQxxAB7C7xOdg0rdBzXPaKPAZWLOa0qZV0ur8bcCjdx3ITySMDPjF+60+pIys9sTj1OZtIfdLBv4NPcXtzeUY9E8iaeKIbCTk88+y7z5xqnONPmf4DIMBFJEJGE9PT0C4yoLoTdJgy4uilfOrrRzb6KbYl/knY8y+pYyo988cdOHgv4iXQTzuKKN9K5cTWrI6kCPFEMyUCNAq9jgZRzjSMiAUA4cLiY0wJgjBlnjIkzxsRFR+v12Uta/7gaTAm4keMmlOG2qXyzfLfVkZSfSDueRUriAtrbNzI270Zuubohdpte/sKbeKIYVgH1RaSOiASRvzN5xlnjzACGup73A34zxhjX8EGuo5bqAPWBvzyQSbmpfHAA93ZuwYqqg+hiX03CikWcysmzOpbyA98s380I21Qyg6PYUqM/A+JqFD2RKlVuF4Nrn8FwYB6wGZhijNkoIqNFpLdrtC+AyiKSBIwEnnZNuxGYAmwC5gIPGWMc7mZSnnFn+zp0uetF8oIqcnfeZH5as8/qSMrHZeY4+GfFHNrZNlHu2seZeH8nwoIDrI6lziL5v7j7lri4OJOQkGB1jDLDLHkTWfQK95V7m0+fuAebrvarizRx+S7q/zKQFmGHCRq5HgLLWR2pTBGR1caYuKLG0zOfVZGkzf3kBIYzIOM7Fv6TZnUc5aOcTsPaJTNoY/uHwGtGaSl4MS0GVbSQitivHsH19rUsWviL1WmUj/pt8wEGnppIZkhV5Mo7rI6jzkOLQRWLve39ZAVG0DntS/5OPmZ1HOWDVvz2E61tWwjqNAoCQ6yOo85Di0EVT3AFuGoE19oTmT9/ptVplI/ZkHyUbulfkxFcFXvcHVbHUUXQYlDFFnLV/ZwMiCBu52ekHM20Oo7yIUvnTSHOthX7NaMgINjqOKoIWgyq+ILLk9tmOB1t61nwq641qOJJPXqKtrvHcSyoKuVa32F1HFUMWgzqgkRc8yAn7BHU2/QRGdl6wpsq2u+/TKKFbRuO9iN1bcFHaDGoCxMUxvGWD3IV61my4OwT3JU608msXBpt+ZjDAVWp1P4uq+OoYtJiUBcspvNwjtoiqLr6PRxO3ztBUpWeZXN/oBlJZLR+BAKCrI6jikmLQV24oDD2NxlGnHM9q5bMsjqN8lIOh5PY9R+QZqtCzevutTqOugBaDOqi1O/5KIeIIGzZW1ZHUV5q3cJJNHImkdZihK4t+BgtBnVR7MFhJDW4m6a5iWxdqWdDq7MYQ8Rf75IiVWnU7X6r06gLpMWgLlqT3o+RTgQsfs3qKMrLbP/zRy7N28aORg9gD9S1BV+jxaAuWlj5CqyrdScNMhNJS/zV6jjKWzidBC59k71UpUUvXVvwRVoMyi1Nej/MflOJzHkvgQ9ewl153v6VU6iZvY3EuvcRFqpXUPVFWgzKLdUrV2JFjbupdWoDB9foeQ1lniMPFr3CdhNDuz4PWJ1GXSS3ikFEKonIfBHZ5voZWcg4zUVkuYhsFJH1IjKwwHtfi8hOEVnnejR3J4+yRuu+D7PbVCV3/kvgdFodR1nowJ8TqJazh8QGI6hcMdTqOOoiubvG8DSw0BhTH1joen22U8DtxpjLgW7A+yISUeD9J4wxzV2PdW7mURa4pHJF/qp9P9WztnPwr0lWx1FWycsm4Pc32GDq0ilez3L2Ze4WQzwwwfV8AtDn7BGMMVuNMdtcz1OANCDazfkqL9Oh7/38Y2rCb6+AI9fqOMoCBxZ/RuW8A2xu/AiVyus1kXyZu8VQ1RiTCuD6WeV8I4tIayAI2F5g8CuuTUzvicg5/zWJyDARSRCRhPT0dDdjK0+rFhHKmkuHE5WTzKE/vrI6jiptOScJWf4uq0xjOvcabHUa5aYii0FEFojIhkIe8RcyIxGpDnwL3GmMOb0h+l9AQ6AVUAl46lzTG2PGGWPijDFx0dG6wuGNro+/nbWmPvalb0JultVxVClKW/Ah4Y4jJDUdSUSYri34uiKLwRhzgzGmSSGP6cAB1xf+6S/+Qu8ULyIVgdnAc8aYFQU+O9Xkywa+Alp7YqGUNaqGl+Pvyx4hIi+dQ0s+tTqOKi2ZRwhL+D8W05IePf5na7LyQe5uSpoBDHU9HwpMP3sEEQkCpgHfGGN+POu906Ui5O+f2OBmHmWxbjf250/TlODl70H2CavjqFKQNu9twpwZJDcfRXhooNVxlAe4WwyvA51FZBvQ2fUaEYkTkfGucQYAHYE7Cjks9TsR+Rv4G4gCXnYzj7JYlQohbG38COUdxzi88H2r46iSlpFGeOJ4fuEqenfranUa5SEB7kxsjDkEXF/I8ATgHtfzicDEc0x/nTvzV97pxp69mb+pFVcnfAqdHoLQSlZHUiXk4C+vEuHMIT1uFBVDdG3BX+iZz8rjosoHs6vZYwQ7TnHk1zesjqNKytE9hG+cyHS5lr6dr7E6jfIgLQZVIm7u1pmZXE1Y4pdwPNXqOKoEHJozBqcxZLQZSQVdW/ArWgyqRFQKCyK1xWOI08HRea9aHUd52sFtRG6dyo/SlZuva2t1GuVhWgyqxAzs3IGfuI4KG7+DwzutjqM86PCsFzhlgsi9aiTlg93aVam8kBaDKjGRYUEciXuUHGPn6C9jrI6jPCU1kUq7ZvO9rRcDrtHrXvojLQZVogZf35of6EbFbT9D2mar4ygPODr7BY6aMAI7jCBM1xb8khaDKlERoUFktRlBhgnh+JwXrI6j3LV7ORHJi/jG3peBVzexOo0qIVoMqsTdcm0LJkhvKu6aB7v+tDqOulhOJxkzn2K/iaRCxwcJDdK1BX+lxaBKXHi5QAKuGsE+U5mMGU/ozXx8lHP9FMofTGRswK0Muqqh1XFUCdISz1duAAATT0lEQVRiUKXizmsb83nQ7ZQ/vBHHuu+tjqMuVM5Jsub+m0RnXZr3vI9yQXarE6kSpMWgSkVIoJ3WNw5jjbMe2XNfgOwMqyOpC5D1+/uEZqUxpfIDxLeItTqOKmFaDKrUdG9anWlVhhOac5Csxe9YHUcV17F92Jd9wCxHWwbdPJD8iyErf6bFoEqNiHBLv5uZ7rgK+4r/g6N7rY6kiuH47H/jdDjZ2HgkTWPDrY6jSoEWgypVDatVJKnZKBxOw7FZz1odRxXBJCdQcetPfCM9ufvGTlbHUaVEi0GVurt7dmSC9CY8aTpmz0qr46hzMYZj054g3YQTdM3jRJXXW3aWFW4Vg4hUEpH5IrLN9TPyHOM5CtykZ0aB4XVEZKVr+smuu70pPxcRGkTF6x/ngIng2H/08FVvlbt+KhGH1jAh5DaGdLzc6jiqFLm7xvA0sNAYUx9Y6HpdmExjTHPXo3eB4W8A77mmPwLc7WYe5SP6X9WQb8PuIOJwIjmJU6yOo86Wm0nWnOfY6KxFq74jCLTrxoWyxN2/7Xhgguv5BPLv21wsrvs8XwdMvZjplW8LsNtof9Nw1jvrkP3L85BzyupIqoCMxR9QIXs/cy4ZwTUNq1kdR5Uyd4uhqjEmFcD1s8o5xgsRkQQRWSEip7/8KwNHjTF5rtfJQMy5ZiQiw1yfkZCenu5mbOUN2tWLZn7NR6iQc4Djv71rdRx12on9BC57n1+drejfb4jVaZQFiiwGEVkgIhsKecRfwHxqGmPigCHA+yJyKVDYwdDmXB9gjBlnjIkzxsRFR0dfwKyVNxt480DmOtsQvPIjOJ5idRwFHJrxb8SZw86WT1E7KszqOMoCRRaDMeYGY0yTQh7TgQMiUh3A9TPtHJ+R4vq5A1gMtAAOAhEicvpKXLGAfjOUMbGRoexr9S9w5nFwuh6+ajXnvrVEbvuRyfae3NL9WqvjKIu4uylpBjDU9XwoMP3sEUQkUkSCXc+jgPbAJmOMARYB/c43vfJ/Q7p2ZErAjURt/xnH3tVWxym7jOHQz09wxJQnvOszeme2MszdYngd6Cwi24DOrteISJyIjHeN0whIEJFE8ovgdWPMJtd7TwEjRSSJ/H0OX7iZR/mgckF2ons8Q7qpyKGfR4E55xZFVYIy//4P0YdW8WPF2+nVSq+eWpaJ8cH/hHFxcSYhIcHqGMqDjDGMe/8F7jv2Aad6fUZo3CCrI5UtOac48s6VpGXayLp7CVfUirI6kSoBIrLatb/3vPTgZOUVRIQO/R8j0VkX5y9PwanDVkcqUw7OepHI7BQW1X1cS0FpMSjv0Tg2kjXNxxCcd4LUyY9aHafMyN6dQOT6z/nZ1plBA26xOo7yAloMyqvc0rsHU8oNoPru6RxfP8fqOP4vL4djk4aRbsKJ7vsGEaF6VRqlxaC8TFCAjbjbxrDNxJI34xFM1nGrI/m15FmvUSVzO7/WeYoOTS+1Oo7yEloMyutcFhPN31e+TERuOjsnP2V1HL91ct8Gqqz7kAX2Dtw8+F6r4ygvosWgvFJ8rz7MCYun7s7vObhxkdVx/I/TwaHvhpFhQojq9y5hes6CKkCLQXklu01odvvbJJtocqc9hDNbL7LnSdtmvkPNUxv5o97jNG/UwOo4ystoMSivVbNaNFtav0L1vH1snKSXy/CUY/u2Ebv2bVYGxNF10Air4ygvpMWgvNp1PQawJKwbjXZ8TfKm5VbH8X3GkPLdfTiMEDHg/wgO1E1I6n9pMSivJiI0GvoBR6QiOT8/SF5OttWRfNramR/T6NRqEho8ymUNGlkdR3kpLQbl9apUqcauNqOpm7eDVT+8ZHUcn5W2bzeXrnmVjYGXc/XAJ6yOo7yYFoPyCa26D2VN+WtouWMcWzfqFVgvlDGG3RMfJMjkUHHAWAICdBOSOjctBuUz6t7+CVkSTM7PD5Gdm2t1HJ/y+/QvaJX5Bxsve4ga9ZtZHUd5OS0G5TMiqsSyr+3zNHFsZvHE162O4zP2JCfTeO1odgbWo8WA56yOo3yAFoPyKY27DmNrhdZ02PURS/9YbHUcr3ciM5s9X99DpJyg/ICx2AICrY6kfIBbxSAilURkvohsc/2MLGSca0VkXYFHloj0cb33tYjsLPBec3fyqDJAhFp3fU22PZSa84excfseqxN5LYfTMO+zp7k6bzl7Wv6L6PqtrI6kfIS7awxPAwuNMfWBha7XZzDGLDLGNDfGNAeuA04BvxYY5YnT7xtj1rmZR5UBwZExyMBviJGDHJl4BweO6VnRhZk86WtuOvIVO6t1o+6NehSSKj53iyEemOB6PgHoU8T4/YBfjDH6P1m5JeKyjhy8+iWuNqtZOHYUmTkOqyN5lVlLltFjy7OklatLnbu+BBGrIykf4m4xVDXGpAK4flYpYvxBwA9nDXtFRNaLyHsiEuxmHlWGVLt+OCm1+zIk83u++upTnE7fu01tSfhrazJ1f7ufQJsQdfePEBRmdSTlY4osBhFZICIbCnnEX8iMRKQ60BSYV2Dwv4CGQCugEnDOayyLyDARSRCRhPT09AuZtfJXIlxyy6ekV2jErSmv8PWshVYnstzugxmkfX8/DWUP3DyegGi9x4K6cEUWgzHmBmNMk0Ie04EDri/801/8aef5qAHANGPMfw9AN8akmnzZwFdA6/PkGGeMiTPGxEVHRxd3+ZS/CyxH1F2TsQcE0j7hEWYlbLU6kWWOZ+Uyc/xL9GIpx9o8TliT7lZHUj7K3U1JM4ChrudDgennGXcwZ21GKlAqQv7+iQ1u5lFlkETWImjg19S3pWCbMYK1uw9bHanU5TmcfPjlN9yX+QWHY28gsuszVkdSPszdYngd6Cwi24DOrteISJyIjD89kojUBmoAS86a/jsR+Rv4G4gCXnYzjyqjAhtcR1bH5+hhW8HvE14g5Wim1ZFK1Yf/+Z37DrxEZlgslW79Emx6ipK6eGKM7+2wi4uLMwkJCVbHUN7GGE58ewuh2+fwXIWX+ffD9xMa5P/XBPph2TYumzuIJgH7CLp/EVTRq6aqwonIamNMXFHj6a8Vyn+IUGHgZ2SFX8rjJ15n9MR5fn+k0rKkgzh/eZqWtiTsN43VUlAeocWg/EtwBcJun0SFAMOQXc/yws+rcfhpOazefYS5E9/iFvsCstuMwN6kqNOIlCoeLQblf6LqE9h/HM1sO7kucSRP/LCCXIfT6lQetXRbOhPGf8i/GU9WjQ4Ed3nR6kjKj2gxKL8kDXvCjR/Qyb6egVse5eGvl/jN2dFzN6Qye8LbvGd/Hy5pQciQiWD3/30pqvRoMSj/deUdyM3jaWVP4oHdjzF8/HyOZ/n2fRx+TNjLX5Ne4fWAsThrdyTwjulQLsLqWMrPaDEo/9a0H7bB33N5QApP7x/JiLGzOJThm/eN/mLpDvZOe5HnA74lr0FPAm+dope7UCVCi0H5vwZdsd/2E3WCjvLykSd55NNppB7znfMcjDG8++sWnPOeZWTgVBxNBxEw8BsI0EuLqZKhxaDKhjodCLhzJtVDcng342me+HgyOw+etDpVkZxOw+gZf3PJ709yb8AcnK2GYe/7qe5TUCVKi0GVHTFXEnDXL0SGBfJ/Oc/ywqffsinluNWpzinP4eSpHxO4MuEJBgUsxnR4HFuPN/WsZlXi9F+YKluqNibwnnmEVYjkU8eLvDXuC1Z74bWVsnIdPPrtMnpsGEUv+wpM5zHI9f/W+yqoUqHFoMqeSnUJvPdXgirFMpZX+ezzTxi7ZLvXnAi3KeU4t3w8n9t2jKKTfT3c+AHS/mGrY6kyRItBlU0VLyHw7nnYqzZiXMBbVFnwMMM+nsH29AzLIuU6nHy4YAsTPnmZz47eRyt7EnLzeLjyDssyqbJJL6Knyrack5jf38G57EOynTbGOvpQ6YaR3N7hMmy20ttss/XACcZ9N4nbjn7CFbYd5FW/koBeb0HMlaWWQfm/4l5ET4tBKYDDO8ma8wwhSXPY64xmcqX76H/rA9SKKl+is3U4Dd/9upzwZa8Qb/uDrJAqhHR/GZr2153MyuO0GJS6CGb7Yo7/53HCT2xjuWnCwfYv0vOGG0pk7WFHajrLJ46mb8ZkAm1OctsMJ/TaxyG4ZMtIlV2lctltEekvIhtFxCki55yZiHQTkS0ikiQiTxcYXkdEVorINhGZLCJB7uRRyl1yaSfCH13B0Wtfo4l9Nz3+7M/Cd24jJSXZY/NwOpws+GkcQWPbccvJbzh6SUcCRqwitNuLWgrKK7i1xiAijQAn8BnwuDHmf36NFxE7sJX8O7wlA6uAwcaYTSIyBfjZGDNJRMYCicaYT4uar64xqNJgTh1m2+RnqLtrMicpR3LVTsQ2bkd43VZQtUnxv8SdTjiyE8e+tezesJy8Hb/TIG8rewPrUD7+LSKbdC7ZBVHKpbhrDG6dPmmM2eya2flGaw0kGWN2uMadBMSLyGbgOmCIa7wJwItAkcWgVGmQ0Eo0uHMs+7cNI3na89Q+sJTwtNmwGAwCUfWR6lfA6Ue1ZhBcAQ5ug9R1kJoIqYk4UxOx5WRgB2JMADttNVnT9Dla9HkUsQdavZhK/Y/SOK8+Bthb4HUy0AaoDBw1xuQVGB5TCnmUuiDV6rek2pOz2HvoJP/3xxq2rvuD2jlJtD6yl+YnllL+7x///8j2YHDkX6TPYQ9hd2BdlmW2I9FZh+AaLejUvgPXXh6LvRSPeFLqQhVZDCKyAKhWyFvPGmOmF2Mehf0PMOcZfq4cw4BhADVr1izGbJXyrBqVwxge34GsHlfxy4ZU3l2+mzV7jlI9MIO7656gV5U0IuUEK09ewoRdESw+FE6F0BAGtKvBQ61rUjtKr4SqfEORxWCMucHNeSQDNQq8jgVSgINAhIgEuNYaTg8/V45xwDjI38fgZialLlpIoJ2+LWLp2yKWDfuO8d3K3byzNoWXt1T/7zhX1Ijgjetq0atZdUIC7RamVerClcampFVAfRGpA+wDBgFDjDFGRBYB/YBJwFCgOGsgSnmNJjHhvHZTM57u3oif1ySz+9ApbmoZQ7NYvXmO8l1uFYOI9AU+AqKB2SKyzhjTVUQuAcYbY3oYY/JEZDgwD7ADXxpjNro+4ilgkoi8DKwFvnAnj1JWCS8XyJ3t61gdQymP0BPclFKqjCiVE9yUUkr5Hy0GpZRSZ9BiUEopdQYtBqWUUmfQYlBKKXUGLQallFJn0GJQSil1Bp88j0FE0oHdFzl5FPmX4yhLdJnLBl1m/+fu8tYyxkQXNZJPFoM7RCShOCd4+BNd5rJBl9n/ldby6qYkpZRSZ9BiUEopdYayWAzjrA5gAV3mskGX2f+VyvKWuX0MSimlzq8srjEopZQ6D78tBhHpJiJbRCRJRJ4u5P1gEZnsen+liNQu/ZSeVYxlHikim0RkvYgsFJFaVuT0pKKWucB4/UTEiIhPH8FSnOUVkQGuv+eNIvJ9aWf0tGL8u64pIotEZK3r33YPK3J6koh8KSJpIrLhHO+LiHzo+jNZLyItPRrAGON3D/JvCLQdqAsEAYlA47PGeRAY63o+CJhsde5SWOZrgVDX8wfKwjK7xqsA/A6sAOKszl3Cf8f1yb/pVaTrdRWrc5fCMo8DHnA9bwzssjq3B5a7I9AS2HCO93sAvwACtAVWenL+/rrG0BpIMsbsMMbkkH/r0PizxokHJrieTwWuFxEpxYyeVuQyG2MWGWNOuV6uIP8+276sOH/PAGOAN4Gs0gxXAoqzvPcCHxtjjgAYY9JKOaOnFWeZDVDR9Tyc89w73lcYY34HDp9nlHjgG5NvBRAhItXPM/4F8ddiiAH2Fnid7BpW6DjGmDzgGFC5VNKVjOIsc0F3k/8bhy8rcplFpAVQwxgzqzSDlZDi/B03ABqIyJ8iskJEupVaupJRnGV+EbhVRJKBOcCI0olmqQv9/35B3Lrnsxcr7Df/sw+/Ks44vqTYyyMitwJxwDUlmqjknXeZRcQGvAfcUVqBSlhx/o4DyN+c1In8NcKlItLEGHO0hLOVlOIs82Dga2PMOyLSDvjWtczOko9nmRL9/vLXNYZkoEaB17H87+rlf8cRkQDyV0HPt+rm7YqzzIjIDcCzQG9jTHYpZSspRS1zBaAJsFhEdpG/LXaGD++ALu6/6+nGmFxjzE5gC/lF4auKs8x3A1MAjDHLgRDyrynkz4r1//1i+WsxrALqi0gdEQkif+fyjLPGmQEMdT3vB/xmXHt1fFSRy+zarPIZ+aXg69ueoYhlNsYcM8ZEGWNqG2Nqk79fpbcxJsGauG4rzr/r/5B/kAEiEkX+pqUdpZrSs4qzzHuA6wFEpBH5xZBeqilL3wzgdtfRSW2BY8aYVE99uF9uSjLG5InIcGAe+Uc1fGmM2Sgio4EEY8wM4AvyVzmTyF9TGGRdYvcVc5nfAsoDP7r2s+8xxvS2LLSbirnMfqOYyzsP6CIimwAH8IQx5pB1qd1TzGUeBXwuIo+RvznlDh//JQ8R+YH8zYFRrn0nLwCBAMaYseTvS+kBJAGngDs9On8f//NTSinlYf66KUkppdRF0mJQSil1Bi0GpZRSZ9BiUEopdQYtBqWUUmfQYlBKKXUGLQallFJn0GJQSil1hv8Hi5BWzfP61iAAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x15177442b0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# find lhs such that res == 0 and substitute this lhs in the namespace\n",
    "lhs = solver.solve_linear('lhs', res)\n",
    "ns = ns(lhs=lhs)\n",
    "\n",
    "# plot solution\n",
    "bezier = domain.sample('bezier', 9)\n",
    "x, u = bezier.eval([ns.x, ns.u])\n",
    "x, g = bezier.eval([ns.x, ns.g])\n",
    "plt.plot(x, u)\n",
    "plt.plot(x, g)\n",
    "# plt.ylim(-2,2)"
   ]
  }
 ],
 "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"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}