FEniCS
/
Navier-Stokes-New
generated from gmarx/git-template

{ "cells": [  {   "cell_type": "code",   "execution_count": 29,   "metadata": {},   "outputs": [    {     "data": {      "text/plain": [       "[<matplotlib.lines.Line2D at 0x7f62a08eee10>,\n",       " <matplotlib.lines.Line2D at 0x7f62a0912128>]"      ]     },     "execution_count": 29,     "metadata": {},     "output_type": "execute_result"    },    {     "data": {      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQIAAAD8CAYAAACcoKqNAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJztnX1wlFW64H9PAgkfEr6CCeFjYDRAGLREIxJHx3FEIFKFVTMjigg65Y41915vOTV3p3Y+tmbm6j/37q29VXtr3d3x7k6NOu4oM1CSGT9ARUv8AE3AGflMSJA1iQkGMQSBJJ0++0f323QaMAc4p9P9vs+viqo+3S/Pc3Ly5tdvn/P2ecQYg6Io0aZguDugKMrwoyJQFEVFoCiKikBRFFQEiqKgIlAUBRWBoiioCBRFQUWgKAowYrgSl5aWmlmzZg1XekWJBA0NDV3GmClDHTdsIpg1axb19fXDlV5RIoGIHLY5Tj8aKIqiIlAURUWgKAoqAkVRUBEoioKFCETkNyJyRER2n+d1EZF/E5GDIvJXEbnWfTcVRfGJzRXBb4HlX/J6LVCZ/PcQ8D8vvVuKomSTIe8jMMa8KSKzvuSQO4GnTGLPs+0iMkFEphpjPrnUzv3617/m9OnTXH311Zca6pz09PSwa9cuAL7xjW94ybFt2zaMMXzta19j8uTJXnK8+eabgL+fwRjDtm3bvOZobGyko6ODyy+/nHnz5nnJEYzTzTffjIh4zeFrnE6ePEl9fT3f/OY3ueWWW5zFFZs9C5Mi+LMxZsE5Xvsz8E/GmLeS7deA/2SMOetuIRF5iMRVAzNnzrzu8OEvv9fhH//xH4f+CRQlovzyl78c8hgRaTDGVA91XFbvLDTGPAE8AVBdXT2kgSoqKhg7diz33nuv037E43Gef/55PvzwQwBqa2tZtGiR0xw9PT08+eSTHD16FIBHHnmECRMmOM3R3NzMs88+SywWY8yYMfz4xz92Gt8Yw9atW3nrrbcAWLx4McuWLXOao7e3l2eeeYaPP/4YgHXr1jF79mynOTo6Onjqqac4deoUYPcHdKHU19fzwgsvADBnzhxWr17tNP7AwAAbN25k7969AKxcudJpfBerBm3AjLT29ORzOUm6BBYvXuwlRyCBnp4eFi5c6CVHIIHJkycza9YsxowZ4zR+ugSuu+46RowY4fxyOpBAa2sr119/vdPYAYEERo4cSVVVFSNGuH/vCyRQWVnJpEmTKCwsdBo/XQK+zlkXIqgD1iVXDxYD3S7mB3yQLoFvfetb3HTTTc5zpEtgzZo1zJw503mOdAmsW7fOuwRWrFjhVQLf/e53mT9/vtP4MFgC999/PxMnTnSeI10Cq1at8iqBZcuWDZ8IROT3wLvAXBFpFZEHReQHIvKD5CEvAi3AQeDfgb/10tNLJFMCN998s/McKgE7hkMCkyZNcp4jUwKurzayJQGwWzX40g87ydWCv3PWIw+oBOxQCdgTJglABO4sVAnYoRKwJ2wSgJCLQCVgh0rAnjBKAEIsApWAHSoBe8IqAQipCFQCdqgE7AmzBCCEIlAJ2KESsCfsEoCQiUAlYIdKwJ4oSABCJAKVgB0qAXuiIgEIiQhUAnaoBOyJkgQgBCJQCdihErAnahKAPBeBSsAOlYA9UZQA5LEIVAJ2qATsiaoEIE9FoBKwQyVgT5QlAHkoApWAHSoBe6IuAcgzEagE7FAJ2KMSSJA3IlAJ2KESsEclcIa8EIFKwA6VgD0qgcEMW1l0W7IhgRMnTmRFAi+//LI3CQChkEBnZydvvvlmXksgHo/nlQQgx0UQj8dpbm4G8CYBSNQeKCoq8iYBgD//+c+UlZV5k0BXV5dXCfT397Nz505ExJsEADZv3kxJSYk3CcRiMa8SiMViHDhwACBvJACWdQ18UF1dberrzyp9kCIej/PYY4+l2lOnTnXeh08+GbzHqu8c48eP9yKB9Bzl5eXOJXDixAl6enpSbd/jNGrUKC8bjabnuPzyy51vNBqLxfj0009Tbd/jlLd1DS6E9JO5srLSefwTJ04MavvIEVzNBFx++eXOczQ1NaUeX3nllc4l0NfXN0gCPsaps7NzUHvGjBnnOfLiSR+nmTNnUlxc7DR++tUr+Bmn7u5u5zEDcloEU6dO5bLLLnNe4CSYGAzwUeCkubmZQ4cOEY/HAfcFToKJweAEnzx5MmvWrHEWH87MCYgIxhhqampYunSp0xzBxGCAjwIn9fX1g0Twve99z2n8gYEBNmzYkGpXVVWxatUqpzm6u7sHnbO5WODEGz7q06WvDrj+ZQWkrw7cfvvtzuNnrg5UVVV5nxgcOXKk0/gweHXA9YkdkD4xeMMNN3hZHdiwYQP79u1j6dKlTJkyxWl8OCOBkydPctdddzmPDzkuAtfoEqEdWnzEjkwJ1NTUOI0PgyWwdu1apk2b5jwHREgEKgE79D4BO8IkAYiICFQCdqgE7AibBCACIlAJ2KESsCOMEoCQi0AlYIdKwI6wSgBCLAKVgB0qATvCLAEIqQhUAnaoBOwIuwQghCJQCdihErAjChIASxGIyHIROSAiB0XkJ+d4faaIvC4iu0TkryJyh/uuDo1KwA6VgB1RkQBYiEBECoHHgVpgPrBaRDLPnP8MrDfGLATuAf6H644OhUrADpWAHVGSANhdESwCDhpjWowxfcCzwJ0ZxxigJPl4PNDurotDoxKwQyVgR9QkAHZfOpoGfJzWbgVuyDjmV8AWEfl7YCywxEnvLFAJ2KESsCOKEgB3k4Wrgd8aY6YDdwBPi8hZsUXkIRGpF5H69O9tXywqATtUAnZEVQJgJ4I2IP0L4tOTz6XzILAewBjzLjAKKM0MZIx5whhTbYypvtRvaakE7FAJ2BFlCYCdCN4HKkVktogUkZgMrMs45v8BtwGISBUJEVz6W/55UAnYoRKwI+oSAAsRGGNiwMPAZmAfidWBPSLyqIgEXyL/B+D7IvIX4PfAA8bTHmgqATtUAnaoBBJYjaox5kXgxYznfpH2eC/wdbddOxuVgB0qATtUAmfImzsLVQJ2qATsUAkMJi9EoBKwQyVgh0rgbHJ289KAbEigpaWFV155xXvxkd27d3uTQF9fX1aKj/iWQENDA3v27FEJZJmcFkF7+5kbFO+66y4mTZp01jbkl0IQ68CBA4wdO5Zvf/vbxONxpzmC+yV2797NvHnzuOWWW/jiiy+cxYdE/wcGBjh+/Di1tbXMnDnT6c8AiQInLS0tFBQUsHr1aoqKipzmOHbsGAB79uyhoqKCFStWcPr0aWfxAVpbWzHGsG/fPr7+9a9z1VVXOR+n9Ptj7rnnHsaPH+80R7C1/KlTp5zFhDwqcKIoymAiUeCkoODM9MUdd7j/MmNLSwv79+/3mmPr1q2pd7UFCxY4/1jT19fHq6++mmr7+Bk6OztpaGjwmqOhoSFV5GT06NHceuutTuPH43FefvnlVNvHz9Dd3c3bb7/tNUdTU1OqPoPrcylnRQBQUVHB2LFjuf76653GbW5uZsuWLYwbN46enh5qa2ud5ggmBk+fPk1hYSEDAwPcdtttTgucpBcfgUSBE9fj1NHRweuvv05JSQnHjx+npqbGeY76+no6OzspKiqir6+Pu+66y2mBk6AqMUBhYSEi4vxnCOYEiouL6e3tpaqqynmOxsZGWlpaGDt2LF988QXXXHON0/h5sWrgkvTVgbVr1zqPf67VAddkrg5UVVU5z5G5OuCjwEnm6oBrMkuT33BD5nflLp3MiUEfBU4aGxtZv349ZWVlXs5ZiJgIdInQjuFYIvRRfMR3afJsrA5kSmDUqFHOc0CERKASsCMs9wmoBC6MSIhAJWCHSsCOsEkAIiAClYAdKgE7wigBCLkIVAJ2qATsCKsEIMQiUAnYoRKwI8wSgJCKQCVgh0rAjrBLAEIoApWAHSoBO6IgAQiZCFQCdqgE7IiKBCBEIlAJ2KESsCNKEoCQiEAlYIdKwI6oSQBCIAKVgB0qATuiKAHIcxGoBOxQCdgRVQlAHotAJWCHSsCOKEsA8lQEKgE7VAJ2RF0CkIciUAnYoRK      "text/plain": [       "<Figure size 432x288 with 1 Axes>"      ]     },     "metadata": {      "needs_background": "light"     },     "output_type": "display_data"    }   ],   "source": [    "from fenics import *\n",    "from mshr import *\n",    "import numpy as np\n",    "\n",    "t = 0\n",    "T = 10\n",    "num_steps = 500\n",    "dt = T/num_steps\n",    "\n",    "mesh = UnitSquareMesh(8,8)\n",    "plot(mesh)"   ]  },  {   "cell_type": "code",   "execution_count": 23,   "metadata": {},   "outputs": [],   "source": [    "V = VectorFunctionSpace(mesh, 'P', 2)\n",    "Q = FunctionSpace(mesh, 'P', 1)\n",    "#V = VectorFunctionSpace(mesh, 'P', 2, dim=10)\n",    "\n",    "u = TrialFunction(V)\n",    "v = TestFunction(V)\n",    "p = TrialFunction(Q)\n",    "q = TestFunction(Q)"   ]  },  {   "cell_type": "code",   "execution_count": 24,   "metadata": {},   "outputs": [],   "source": [    "def boundary(x, on_boundary):\n",    "    return near(x[0], 0)\n",    "\n",    "boundary = 'near(x[0], 0)'"   ]  },  {   "cell_type": "code",   "execution_count": 25,   "metadata": {},   "outputs": [],   "source": [    "# Define boundaries\n",    "inflow   = 'near(x[0], 0)'\n",    "outflow  = 'near(x[0], 1)'\n",    "walls   = 'near(x[1], 0) || near(x[1], 1)'\n",    "\n",    "# Define boundary conditions\n",    "bcu_noslip  = DirichletBC(V, Constant((0, 0)), walls)\n",    "bcp_inflow  = DirichletBC(Q, Constant(8), inflow)\n",    "bcp_outflow = DirichletBC(Q, Constant(0), outflow)\n",    "bcu = [bcu_noslip]\n",    "bcp = [bcp_inflow, bcp_outflow]"   ]  },  {   "cell_type": "code",   "execution_count": 42,   "metadata": {},   "outputs": [],   "source": [    "u_n = Function(V)\n",    "#u_n.name = \"u_n\"\n",    "U = 0.5*(u_n + u)\n",    "\n",    "p_n = Function(Q)\n",    "#p_n.name = \"p_n\""   ]  },  {   "cell_type": "code",   "execution_count": 43,   "metadata": {},   "outputs": [],   "source": [    "U   = 0.5*(u_n + u)\n",    "mu = Constant(1)\n",    "rho = Constant(1)\n",    "n   = FacetNormal(mesh)\n",    "f   = Constant((0, 0))\n",    "k   = Constant(dt)\n",    "#mu  = Constant(mu)\n",    "#rho = Constant(rho)"   ]  },  {   "cell_type": "code",   "execution_count": 44,   "metadata": {},   "outputs": [    {     "ename": "TypeError",     "evalue": "'<' not supported between instances of 'Mesh' and 'Mesh'",     "output_type": "error",     "traceback": [      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",      "\u001b[0;31mTypeError\u001b[0m                                 Traceback (most recent call last)",      "\u001b[0;32m<ipython-input-44-0cd64d6e94f6>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m     12\u001b[0m     \u001b[0;34m+\u001b[0m \u001b[0minner\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0msigma\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mU\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mp_n\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mepsilon\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mv\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mdx\u001b[0m\u001b[0;31m \u001b[0m\u001b[0;31m\\\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     13\u001b[0m     \u001b[0;34m+\u001b[0m \u001b[0mdot\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mp_n\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mn\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mv\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mds\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0mdot\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmu\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mnabla_grad\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mU\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mn\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mv\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mds\u001b[0m\u001b[0;31m \u001b[0m\u001b[0;31m\\\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 14\u001b[0;31m     \u001b[0;34m-\u001b[0m \u001b[0mdot\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mf\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mv\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mdx\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     15\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     16\u001b[0m \u001b[0ma1\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mlhs\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mF1\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",      "\u001b[0;32m/usr/local/lib/python3.6/dist-packages/ufl/measure.py\u001b[0m in \u001b[0;36m__rmul__\u001b[0;34m(self, integrand)\u001b[0m\n\u001b[1;32m    435\u001b[0m         \u001b[0mdomain\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mufl_domain\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    436\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0mdomain\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 437\u001b[0;31m             \u001b[0mdomains\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mextract_domains\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mintegrand\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    438\u001b[0m             \u001b[0;32mif\u001b[0m \u001b[0mlen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mdomains\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    439\u001b[0m                 \u001b[0mdomain\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mdomains\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",      "\u001b[0;32m/usr/local/lib/python3.6/dist-packages/ufl/domain.py\u001b[0m in \u001b[0;36mextract_domains\u001b[0;34m(expr)\u001b[0m\n\u001b[1;32m    353\u001b[0m     \u001b[0;32mfor\u001b[0m \u001b[0mt\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mtraverse_unique_terminals\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mexpr\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    354\u001b[0m         \u001b[0mdomainlist\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mextend\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mufl_domains\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 355\u001b[0;31m     \u001b[0;32mreturn\u001b[0m \u001b[0msorted\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mjoin_domains\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mdomainlist\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    356\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    357\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",      "\u001b[0;31mTypeError\u001b[0m: '<' not supported between instances of 'Mesh' and 'Mesh'"     ]    }   ],   "source": [    "# Define strain-rate tensor\n",    "def epsilon(u):\n",    "    return sym(nabla_grad(u))\n",    "\n",    "# Define stress tensor\n",    "def sigma(u, p):\n",    "    return 2*mu*epsilon(u) - p*Identity(len(u))\n",    "\n",    "# Define variational problem for step 1\n",    "F1 = rho*dot((u - u_n) / k, v)*dx + \\\n",    "    rho*dot(dot(u_n, nabla_grad(u_n)), v)*dx \\\n",    "    + inner(sigma(U, p_n), epsilon(v))*dx \\\n",    "    + dot(p_n*n, v)*ds - dot(mu*nabla_grad(U)*n, v)*ds \\\n",    "    - dot(f,v)*dx\n",    "\n",    "a1 = lhs(F1)\n",    "L1 = rhs(F1)"   ]  },  {   "cell_type": "code",   "execution_count": 50,   "metadata": {},   "outputs": [    {     "name": "stdout",     "output_type": "stream",     "text": [      "t = 0.02: error = 0.84\n",      "max u: 0.16\n",      "t = 0.04: error = 0.684\n",      "max u: 0.315688161445\n",      "t = 0.06: error = 0.548\n",      "max u: 0.451967018333\n",      "t = 0.08: error = 0.444\n",      "max u: 0.555932937013\n",      "t = 0.10: error = 0.365\n",      "max u: 0.634821697889\n",      "t = 0.12: error = 0.3\n",      "max u: 0.700458664738\n",      "t = 0.14: error = 0.246\n",      "max u: 0.754470057841\n",      "t = 0.16: error = 0.202\n",      "max u: 0.798401077724\n",      "t = 0.18: error = 0.165\n",      "max u: 0.834736754803\n",      "t = 0.20: error = 0.136\n",      "max u: 0.864370379848\n",      "t = 0.22: error = 0.111\n",      "max u: 0.888757133065\n",      "t = 0.24: error = 0.0913\n",      "max u: 0.908746179289\n",      "t = 0.26: error = 0.0749\n",      "max u: 0.925128943939\n",      "t = 0.28: error = 0.0614\n",      "max u: 0.938595908082\n",      "t = 0.30: error = 0.0504\n",      "max u: 0.949614547115\n",      "t = 0.32: error = 0.0413\n",      "max u: 0.958678413639\n",      "t = 0.34: error = 0.0339\n",      "max u: 0.966095207289\n",      "t = 0.36: error = 0.0278\n",      "max u: 0.97219191582\n",      "t = 0.38: error = 0.0228\n",      "max u: 0.977186437734\n",      "t = 0.40: error = 0.0187\n",      "max u: 0.981286118546\n",      "t = 0.42: error = 0.0154\n",      "max u: 0.984650083709\n",      "t = 0.44: error = 0.0126\n",      "max u: 0.987406612232\n",      "t = 0.46: error = 0.0103\n",      "max u: 0.989672435868\n",      "t = 0.48: error = 0.00848\n",      "max u: 0.991525862918\n",      "t = 0.50: error = 0.00696\n",      "max u: 0.993052001726\n",      "t = 0.52: error = 0.00571\n",      "max u: 0.99429840471\n",      "t = 0.54: error = 0.00468\n",      "max u: 0.9953261608\n",      "t = 0.56: error = 0.00384\n",      "max u: 0.996164570633\n",      "t = 0.58: error = 0.00315\n",      "max u: 0.996856517384\n",      "t = 0.60: error = 0.00259\n",      "max u: 0.997420565004\n",      "t = 0.62: error = 0.00212\n",      "max u: 0.997886381956\n",      "t = 0.64: error = 0.00175\n",      "max u: 0.99826587716\n",      "t = 0.66: error = 0.00143\n",      "max u: 0.998579465891\n",      "t = 0.68: error = 0.00118\n",      "max u: 0.998834779083\n",      "t = 0.70: error = 0.000965\n",      "max u: 0.999045922679\n",      "t = 0.72: error = 0.000799\n",      "max u: 0.999217645778\n",      "t = 0.74: error = 0.000651\n",      "max u: 0.99935987071\n",      "t = 0.76: error = 0.000644\n",      "max u: 0.999475305263\n",      "t = 0.78: error = 0.000445\n",      "max u: 0.999571183133\n",      "t = 0.80: error = 0.000587\n",      "max u: 0.999648697626\n",      "t = 0.82: error = 0.000447\n",      "max u: 0.999713420101\n",      "t = 0.84: error = 0.00054\n",      "max u: 0.999765376725\n",      "t = 0.86: error = 0.000443\n",      "max u: 0.999809165564\n",      "t = 0.88: error = 0.000503\n",      "max u: 0.999843886616\n",      "t = 0.90: error = 0.000434\n",      "max u: 0.999873617869\n",      "t = 0.92: error = 0.000471\n",      "max u: 0.999896707363\n",      "t = 0.94: error = 0.000421\n",      "max u: 0.99991700559\n",      "t = 0.96: error = 0.000443\n",      "max u: 0.99993223784\n",      "t = 0.98: error = 0.000407\n",      "max u: 0.999946212975\n",      "t = 1.00: error = 0.000419\n",      "max u: 0.999956130296\n",      "t = 1.02: error = 0.000391\n",      "max u: 0.999965873668\n",      "t = 1.04: error = 0.000397\n",      "max u: 0.999972188576\n",      "t = 1.06: error = 0.000376\n",      "max u: 0.999979106649\n",      "t = 1.08: error = 0.000377\n",      "max u: 0.999982972641\n",      "t = 1.10: error = 0.00036\n",      "max u: 0.999988011482\n",      "t = 1.12: error = 0.000359\n",      "max u: 0.999990205399\n",      "t = 1.14: error = 0.000345\n",      "max u: 0.999994001581\n",      "t = 1.16: error = 0.000342\n",      "max u: 0.999995046491\n",      "t = 1.18: error = 0.00033\n",      "max u: 0.999998028507\n",      "t = 1.20: error = 0.000326\n",      "max u: 0.999998276527\n",      "t = 1.22: error = 0.000316\n",      "max u: 1.00000073292\n",      "t = 1.24: error = 0.000311\n",      "max u: 1.00000042106\n",      "t = 1.26: error = 0.000302\n",      "max u: 1.0000025462\n",      "t = 1.28: error = 0.000297\n",      "max u: 1.000001834\n",      "t = 1.30: error = 0.000289\n",      "max u: 1.0000037588\n",      "t = 1.32: error = 0.000284\n",      "max u: 1.00000275375\n",      "t = 1.34: error = 0.000277\n",      "max u: 1.00000456632\n",      "t = 1.36: error = 0.000272\n",      "max u: 1.00000334098\n",      "t = 1.38: error = 0.000265\n",      "max u: 1.00000510051\n",      "t = 1.40: error = 0.00026\n",      "max u: 1.00000370399\n",      "t = 1.42: error = 0.000254\n",      "max u: 1.00000545008\n",      "t = 1.44: error = 0.000249\n",      "max u: 1.00000391588\n",      "t = 1.46: error = 0.000243\n",      "max u: 1.0000056748\n",      "t = 1.48: error = 0.000238\n",      "max u: 1.00000402601\n",      "t = 1.50: error = 0.000233\n",      "max u: 1.00000581495\n",      "t = 1.52: error = 0.000229\n",      "max u: 1.00000406771\n",      "t = 1.54: error = 0.000223\n",      "max u: 1.00000589766\n",      "t = 1.56: error = 0.000219\n",      "max u: 1.00000406359\n",      "t = 1.58: error = 0.000214\n",      "max u: 1.00000594124\n",      "t = 1.60: error = 0.00021\n",      "max u: 1.00000402893\n",      "t = 1.62: error = 0.000206\n",      "max u: 1.00000595803\n",      "t = 1.64: error = 0.000202\n",      "max u: 1.00000397412\n",      "t = 1.66: error = 0.000197\n",      "max u: 1.00000595638\n",      "t = 1.68: error = 0.000194\n",      "max u: 1.00000390623\n",      "t = 1.70: error = 0.00019\n",      "max u: 1.00000594193\n",      "t = 1.72: error = 0.000186\n",      "max u: 1.00000383007\n",      "t = 1.74: error = 0.000182\n",      "max u: 1.00000591851\n",      "t = 1.76: error = 0.000179\n",      "max u: 1.00000374894\n",      "t = 1.78: error = 0.000175\n",      "max u: 1.00000588871\n",      "t = 1.80: error = 0.000172\n",      "max u: 1.00000366508\n",      "t = 1.82: error = 0.000168\n",      "max u: 1.00000585431\n",      "t = 1.84: error = 0.000165\n",      "max u: 1.00000358003\n",      "t = 1.86: error = 0.000162\n",      "max u: 1.00000581654\n",      "t = 1.88: error = 0.000159\n",      "max u: 1.00000349485\n",      "t = 1.90: error = 0.000156\n",      "max u: 1.00000577623\n",      "t = 1.92: error = 0.000153\n",      "max u: 1.00000341026\n",      "t = 1.94: error = 0.00015\n",      "max u: 1.00000573399\n",      "t = 1.96: error = 0.000148\n",      "max u: 1.00000332675\n",      "t = 1.98: error = 0.000145\n",      "max u: 1.00000569024\n",      "t = 2.00: error = 0.000142\n",      "max u: 1.00000324465\n",      "t = 2.02: error = 0.000139\n",      "max u: 1.00000564528\n",      "t = 2.04: error = 0.000137\n",      "max u: 1.00000316419\n",      "t = 2.06: error = 0.000134\n",      "max u: 1.00000559933\n",      "t = 2.08: error = 0.000132\n",      "max u: 1.00000308552\n",      "t = 2.10: error = 0.00013\n",      "max u: 1.00000555257\n",      "t = 2.12: error = 0.000127\n",      "max u: 1.0000030087\n",      "t = 2.14: error = 0.000125\n",      "max u: 1.00000550513\n",      "t = 2.16: error = 0.000123\n",      "max u: 1.00000293381\n",      "t = 2.18: error = 0.000121\n",      "max u: 1.00000545711\n",      "t = 2.20: error = 0.000119\n",      "max u: 1.00000286085\n",      "t = 2.22: error = 0.000116\n",      "max u: 1.00000540861\n",      "t = 2.24: error = 0.000115\n",      "max u: 1.00000278982\n",      "t = 2.26: error = 0.000112\n",      "max u: 1.00000535969\n",      "t = 2.28: error = 0.000111\n",      "max u: 1.00000272072\n",      "t = 2.30: error = 0.000109\n",      "max u: 1.00000531042\n",      "t = 2.32: error = 0.000107\n",      "max u: 1.00000265352\n",      "t = 2.34: error = 0.000105\n",      "max u: 1.00000526086\n",      "t = 2.36: error = 0.000103\n",      "max u: 1.00000258819\n",      "t = 2.38: error = 0.000101\n",      "max u: 1.00000521104\n",      "t = 2.40: error = 9.99e-05\n",      "max u: 1.00000252468\n",      "t = 2.42: error = 9.81e-05\n",      "max u: 1.00000516102\n",      "t = 2.44: error = 9.66e-05\n",      "max u: 1.00000246297\n",      "t = 2.46: error = 9.49e-05\n",      "max u: 1.00000511083\n",      "t = 2.48: error = 9.35e-05\n",      "max u: 1.000002403\n",      "t = 2.50: error = 9.18e-05\n",      "max u: 1.00000506051\n",      "t = 2.52: error = 9.05e-05\n",      "max u: 1.00000234473\n",      "t = 2.54: error = 8.89e-05\n",      "max u: 1.00000501009\n",      "t = 2.56: error = 8.76e-05\n",      "max u: 1.00000228812\n",      "t = 2.58: error = 8.6e-05\n",      "max u: 1.0000049596\n",      "t = 2.60: error = 8.48e-05\n",      "max u: 1.00000223312\n",      "t = 2.62: error = 8.33e-05\n",      "max u: 1.00000490907\n",      "t = 2.64: error = 8.21e-05\n",      "max u: 1.00000217968\n",      "t = 2.66: error = 8.07e-05\n",      "max u: 1.00000485851\n",      "t = 2.68: error = 7.96e-05\n",      "max u: 1.00000212776\n",      "t = 2.70: error = 7.82e-05\n",      "max u: 1.00000480797\n",      "t = 2.72: error = 7.71e-05\n",      "max u: 1.00000207732\n",      "t = 2.74: error = 7.58e-05\n",      "max u: 1.00000475744\n",      "t = 2.76: error = 7.48e-05\n",      "max u: 1.0000020283\n",      "t = 2.78: error = 7.35e-05\n",      "max u: 1.00000470696\n",      "t = 2.80: error = 7.25e-05\n",      "max u: 1.00000198067\n",      "t = 2.82: error = 7.13e-05\n",      "max u: 1.00000465655\n",      "t = 2.84: error = 7.04e-05\n",      "max u: 1.00000193438\n",      "t = 2.86: error = 6.92e-05\n",      "max u: 1.00000460621\n",      "t = 2.88: error = 6.83e-05\n",      "max u: 1.00000188939\n",      "t = 2.90: error = 6.72e-05\n",      "max u: 1.00000455597\n",      "t = 2.92: error = 6.63e-05\n",      "max u: 1.00000184566\n",      "t = 2.94: error = 6.52e-05\n",      "max u: 1.00000450584\n",      "t = 2.96: error = 6.43e-05\n",      "max u: 1.00000180315\n",      "t = 2.98: error = 6.33e-05\n",      "max u: 1.00000445583\n",      "t = 3.00: error = 6.25e-05\n",      "max u: 1.00000176237\n",      "t = 3.02: error = 6.15e-05\n",      "max u: 1.00000440597\n",      "t = 3.04: error = 6.07e-05\n",      "max u: 1.0000017347\n",      "t = 3.06: error = 5.97e-05\n",      "max u: 1.00000435625\n",      "t = 3.08: error = 5.9e-05\n",      "max u: 1.00000170731\n",      "t = 3.10: error = 5.8e-05\n",      "max u: 1.0000043067\n",      "t = 3.12: error = 5.73e-05\n",      "max u: 1.00000168021\n",      "t = 3.14: error = 5.64e-05\n",      "max u: 1.00000425731\n",      "t = 3.16: error = 5.57e-05\n",      "max u: 1.0000016534\n",      "t = 3.18: error = 5.49e-05\n",      "max u: 1.00000420812\n",      "t = 3.20: error = 5.42e-05\n",      "max u: 1.00000162688\n",      "t = 3.22: error = 5.33e-05\n",      "max u: 1.00000415911\n",      "t = 3.24: error = 5.27e-05\n",      "max u: 1.00000160066\n",      "t = 3.26: error = 5.19e-05\n",      "max u: 1.00000411031\n",      "t = 3.28: error = 5.13e-05\n",      "max u: 1.00000157472\n",      "t = 3.30: error = 5.05e-05\n",      "max u: 1.00000406172\n",      "t = 3.32: error = 4.99e-05\n",      "max u: 1.00000154907\n",      "t = 3.34: error = 4.91e-05\n",      "max u: 1.00000401334\n",      "t = 3.36: error = 4.86e-05\n",      "max u: 1.00000152372\n",      "t = 3.38: error = 4.78e-05\n",      "max u: 1.0000039652\n",      "t = 3.40: error = 4.73e-05\n",      "max u: 1.00000149865\n",      "t = 3.42: error = 4.66e-05\n",      "max u: 1.00000391729\n",      "t = 3.44: error = 4.61e-05\n",      "max u: 1.00000147387\n",      "t = 3.46: error = 4.54e-05\n",      "max u: 1.00000386963\n"     ]    },    {     "name": "stdout",     "output_type": "stream",     "text": [      "t = 3.48: error = 4.49e-05\n",      "max u: 1.00000144937\n",      "t = 3.50: error = 4.42e-05\n",      "max u: 1.00000382221\n",      "t = 3.52: error = 4.37e-05\n",      "max u: 1.00000142516\n",      "t = 3.54: error = 4.31e-05\n",      "max u: 1.00000377505\n",      "t = 3.56: error = 4.26e-05\n",      "max u: 1.00000140124\n",      "t = 3.58: error = 4.2e-05\n",      "max u: 1.00000372815\n",      "t = 3.60: error = 4.15e-05\n",      "max u: 1.0000013776\n",      "t = 3.62: error = 4.09e-05\n",      "max u: 1.00000368153\n",      "t = 3.64: error = 4.05e-05\n",      "max u: 1.00000135424\n",      "t = 3.66: error = 3.99e-05\n",      "max u: 1.00000363517\n",      "t = 3.68: error = 3.95e-05\n",      "max u: 1.00000133116\n",      "t = 3.70: error = 3.89e-05\n",      "max u: 1.0000035891\n",      "t = 3.72: error = 3.85e-05\n",      "max u: 1.00000130836\n",      "t = 3.74: error = 3.8e-05\n",      "max u: 1.00000354332\n",      "t = 3.76: error = 3.76e-05\n",      "max u: 1.00000128584\n",      "t = 3.78: error = 3.7e-05\n",      "max u: 1.00000349782\n",      "t = 3.80: error = 3.67e-05\n",      "max u: 1.0000012636\n",      "t = 3.82: error = 3.61e-05\n",      "max u: 1.00000345262\n",      "t = 3.84: error = 3.58e-05\n",      "max u: 1.00000124164\n",      "t = 3.86: error = 3.53e-05\n",      "max u: 1.00000340772\n",      "t = 3.88: error = 3.49e-05\n",      "max u: 1.00000121994\n",      "t = 3.90: error = 3.44e-05\n",      "max u: 1.00000336312\n",      "t = 3.92: error = 3.41e-05\n",      "max u: 1.00000119852\n",      "t = 3.94: error = 3.36e-05\n",      "max u: 1.00000331884\n",      "t = 3.96: error = 3.33e-05\n",      "max u: 1.00000117737\n",      "t = 3.98: error = 3.29e-05\n",      "max u: 1.00000327486\n",      "t = 4.00: error = 3.25e-05\n",      "max u: 1.00000115649\n",      "t = 4.02: error = 3.21e-05\n",      "max u: 1.0000032312\n",      "t = 4.04: error = 3.18e-05\n",      "max u: 1.00000113587\n",      "t = 4.06: error = 3.14e-05\n",      "max u: 1.00000318786\n",      "t = 4.08: error = 3.11e-05\n",      "max u: 1.00000111552\n",      "t = 4.10: error = 3.06e-05\n",      "max u: 1.00000314484\n",      "t = 4.12: error = 3.04e-05\n",      "max u: 1.00000109544\n",      "t = 4.14: error = 3e-05\n",      "max u: 1.00000310215\n",      "t = 4.16: error = 2.97e-05\n",      "max u: 1.00000107561\n",      "t = 4.18: error = 2.93e-05\n",      "max u: 1.00000305978\n",      "t = 4.20: error = 2.9e-05\n",      "max u: 1.00000105605\n",      "t = 4.22: error = 2.86e-05\n",      "max u: 1.00000301775\n",      "t = 4.24: error = 2.84e-05\n",      "max u: 1.00000103674\n",      "t = 4.26: error = 2.8e-05\n",      "max u: 1.00000297605\n",      "t = 4.28: error = 2.78e-05\n",      "max u: 1.00000101769\n",      "t = 4.30: error = 2.74e-05\n",      "max u: 1.00000293469\n",      "t = 4.32: error = 2.72e-05\n",      "max u: 1.00000099889\n",      "t = 4.34: error = 2.68e-05\n",      "max u: 1.00000289366\n",      "t = 4.36: error = 2.66e-05\n",      "max u: 1.00000098034\n",      "t = 4.38: error = 2.62e-05\n",      "max u: 1.00000285298\n",      "t = 4.40: error = 2.6e-05\n",      "max u: 1.00000096204\n",      "t = 4.42: error = 2.57e-05\n",      "max u: 1.00000281263\n",      "t = 4.44: error = 2.55e-05\n",      "max u: 1.00000094399\n",      "t = 4.46: error = 2.51e-05\n",      "max u: 1.00000277263\n",      "t = 4.48: error = 2.49e-05\n",      "max u: 1.00000092619\n",      "t = 4.50: error = 2.46e-05\n",      "max u: 1.00000273298\n",      "t = 4.52: error = 2.44e-05\n",      "max u: 1.00000090863\n",      "t = 4.54: error = 2.41e-05\n",      "max u: 1.00000269367\n",      "t = 4.56: error = 2.39e-05\n",      "max u: 1.00000089131\n",      "t = 4.58: error = 2.36e-05\n",      "max u: 1.0000026547\n",      "t = 4.60: error = 2.34e-05\n",      "max u: 1.00000087423\n",      "t = 4.62: error = 2.31e-05\n",      "max u: 1.00000261609\n",      "t = 4.64: error = 2.29e-05\n",      "max u: 1.00000085739\n",      "t = 4.66: error = 2.27e-05\n",      "max u: 1.00000257783\n",      "t = 4.68: error = 2.25e-05\n",      "max u: 1.00000084078\n",      "t = 4.70: error = 2.22e-05\n",      "max u: 1.00000253991\n",      "t = 4.72: error = 2.2e-05\n",      "max u: 1.00000082441\n",      "t = 4.74: error = 2.18e-05\n",      "max u: 1.00000250235\n",      "t = 4.76: error = 2.16e-05\n",      "max u: 1.00000080827\n",      "t = 4.78: error = 2.13e-05\n",      "max u: 1.00000246514\n",      "t = 4.80: error = 2.12e-05\n",      "max u: 1.00000079236\n",      "t = 4.82: error = 2.09e-05\n",      "max u: 1.00000242828\n",      "t = 4.84: error = 2.08e-05\n",      "max u: 1.00000077667\n",      "t = 4.86: error = 2.05e-05\n",      "max u: 1.00000239178\n",      "t = 4.88: error = 2.04e-05\n",      "max u: 1.00000076121\n",      "t = 4.90: error = 2.01e-05\n",      "max u: 1.00000235562\n",      "t = 4.92: error = 2e-05\n",      "max u: 1.00000074597\n",      "t = 4.94: error = 1.97e-05\n",      "max u: 1.00000231982\n",      "t = 4.96: error = 1.96e-05\n",      "max u: 1.00000073096\n",      "t = 4.98: error = 1.94e-05\n",      "max u: 1.00000228437\n",      "t = 5.00: error = 1.92e-05\n",      "max u: 1.00000071616\n",      "t = 5.02: error = 1.9e-05\n",      "max u: 1.00000224927\n",      "t = 5.04: error = 1.89e-05\n",      "max u: 1.00000070157\n",      "t = 5.06: error = 1.86e-05\n",      "max u: 1.00000221453\n",      "t = 5.08: error = 1.85e-05\n",      "max u: 1.00000068721\n",      "t = 5.10: error = 1.83e-05\n",      "max u: 1.00000218014\n",      "t = 5.12: error = 1.82e-05\n",      "max u: 1.00000067305\n",      "t = 5.14: error = 1.8e-05\n",      "max u: 1.00000214609\n",      "t = 5.16: error = 1.78e-05\n",      "max u: 1.0000006591\n",      "t = 5.18: error = 1.76e-05\n",      "max u: 1.0000021124\n",      "t = 5.20: error = 1.75e-05\n",      "max u: 1.00000064537\n",      "t = 5.22: error = 1.73e-05\n",      "max u: 1.00000207906\n",      "t = 5.24: error = 1.72e-05\n",      "max u: 1.00000063183\n",      "t = 5.26: error = 1.7e-05\n",      "max u: 1.00000204607\n",      "t = 5.28: error = 1.69e-05\n",      "max u: 1.0000006185\n",      "t = 5.30: error = 1.67e-05\n",      "max u: 1.00000201343\n",      "t = 5.32: error = 1.66e-05\n",      "max u: 1.00000060538\n",      "t = 5.34: error = 1.64e-05\n",      "max u: 1.00000198113\n",      "t = 5.36: error = 1.63e-05\n",      "max u: 1.00000059245\n",      "t = 5.38: error = 1.61e-05\n",      "max u: 1.00000194918\n",      "t = 5.40: error = 1.6e-05\n",      "max u: 1.00000058035\n",      "t = 5.42: error = 1.58e-05\n",      "max u: 1.00000191758\n",      "t = 5.44: error = 1.57e-05\n",      "max u: 1.00000057923\n",      "t = 5.46: error = 1.56e-05\n",      "max u: 1.00000188632\n",      "t = 5.48: error = 1.55e-05\n",      "max u: 1.00000057807\n",      "t = 5.50: error = 1.53e-05\n",      "max u: 1.00000185541\n",      "t = 5.52: error = 1.52e-05\n",      "max u: 1.00000057687\n",      "t = 5.54: error = 1.5e-05\n",      "max u: 1.00000182483\n",      "t = 5.56: error = 1.49e-05\n",      "max u: 1.00000057564\n",      "t = 5.58: error = 1.48e-05\n",      "max u: 1.0000017946\n",      "t = 5.60: error = 1.47e-05\n",      "max u: 1.00000057438\n",      "t = 5.62: error = 1.45e-05\n",      "max u: 1.00000176471\n",      "t = 5.64: error = 1.44e-05\n",      "max u: 1.00000057308\n",      "t = 5.66: error = 1.43e-05\n",      "max u: 1.00000173515\n",      "t = 5.68: error = 1.42e-05\n",      "max u: 1.00000057174\n",      "t = 5.70: error = 1.4e-05\n",      "max u: 1.00000170594\n",      "t = 5.72: error = 1.4e-05\n",      "max u: 1.00000057038\n",      "t = 5.74: error = 1.38e-05\n",      "max u: 1.00000167705\n",      "t = 5.76: error = 1.37e-05\n",      "max u: 1.00000056898\n",      "t = 5.78: error = 1.36e-05\n",      "max u: 1.00000164851\n",      "t = 5.80: error = 1.35e-05\n",      "max u: 1.00000056755\n",      "t = 5.82: error = 1.34e-05\n",      "max u: 1.00000162029\n",      "t = 5.84: error = 1.33e-05\n",      "max u: 1.0000005661\n",      "t = 5.86: error = 1.32e-05\n",      "max u: 1.0000015924\n",      "t = 5.88: error = 1.31e-05\n",      "max u: 1.00000056461\n",      "t = 5.90: error = 1.29e-05\n",      "max u: 1.00000156484\n",      "t = 5.92: error = 1.29e-05\n",      "max u: 1.00000056309\n",      "t = 5.94: error = 1.27e-05\n",      "max u: 1.00000153761\n",      "t = 5.96: error = 1.27e-05\n",      "max u: 1.00000056154\n",      "t = 5.98: error = 1.25e-05\n",      "max u: 1.0000015107\n",      "t = 6.00: error = 1.25e-05\n",      "max u: 1.00000055997\n",      "t = 6.02: error = 1.23e-05\n",      "max u: 1.00000148411\n",      "t = 6.04: error = 1.23e-05\n",      "max u: 1.00000055837\n",      "t = 6.06: error = 1.21e-05\n",      "max u: 1.00000145785\n",      "t = 6.08: error = 1.21e-05\n",      "max u: 1.00000055674\n",      "t = 6.10: error = 1.2e-05\n",      "max u: 1.0000014319\n",      "t = 6.12: error = 1.19e-05\n",      "max u: 1.00000055509\n",      "t = 6.14: error = 1.18e-05\n",      "max u: 1.00000140628\n",      "t = 6.16: error = 1.17e-05\n",      "max u: 1.00000055341\n",      "t = 6.18: error = 1.16e-05\n",      "max u: 1.00000138096\n",      "t = 6.20: error = 1.15e-05\n",      "max u: 1.00000055171\n",      "t = 6.22: error = 1.14e-05\n",      "max u: 1.00000135596\n",      "t = 6.24: error = 1.14e-05\n",      "max u: 1.00000054998\n",      "t = 6.26: error = 1.12e-05\n",      "max u: 1.00000133127\n",      "t = 6.28: error = 1.12e-05\n",      "max u: 1.00000054823\n",      "t = 6.30: error = 1.11e-05\n",      "max u: 1.00000130689\n",      "t = 6.32: error = 1.1e-05\n",      "max u: 1.00000054645\n",      "t = 6.34: error = 1.09e-05\n",      "max u: 1.00000128281\n",      "t = 6.36: error = 1.09e-05\n",      "max u: 1.00000054465\n",      "t = 6.38: error = 1.07e-05\n",      "max u: 1.00000125904\n",      "t = 6.40: error = 1.07e-05\n",      "max u: 1.00000054283\n",      "t = 6.42: error = 1.06e-05\n",      "max u: 1.00000123557\n",      "t = 6.44: error = 1.05e-05\n",      "max u: 1.00000054099\n",      "t = 6.46: error = 1.04e-05\n",      "max u: 1.0000012124\n",      "t = 6.48: error = 1.04e-05\n",      "max u: 1.00000053913\n",      "t = 6.50: error = 1.03e-05\n",      "max u: 1.00000118953\n",      "t = 6.52: error = 1.02e-05\n",      "max u: 1.00000053725\n",      "t = 6.54: error = 1.01e-05\n",      "max u: 1.00000116695\n",      "t = 6.56: error = 1.01e-05\n",      "max u: 1.00000053535\n",      "t = 6.58: error = 9.98e-06\n",      "max u: 1.00000114467\n",      "t = 6.60: error = 9.93e-06\n",      "max u: 1.00000053343\n",      "t = 6.62: error = 9.84e-06\n",      "max u: 1.00000112267\n",      "t = 6.64: error = 9.79e-06\n",      "max u: 1.00000053149\n",      "t = 6.66: error = 9.7e-06\n",      "max u: 1.00000110097\n",      "t = 6.68: error = 9.65e-06\n",      "max u: 1.00000052953\n",      "t = 6.70: error = 9.56e-06\n",      "max u: 1.00000107955\n",      "t = 6.72: error = 9.52e-06\n",      "max u: 1.00000052756\n",      "t = 6.74: error = 9.43e-06\n",      "max u: 1.00000105841\n",      "t = 6.76: error = 9.38e-06\n",      "max u: 1.00000052556\n",      "t = 6.78: error = 9.29e-06\n",      "max u: 1.00000103755\n",      "t = 6.80: error = 9.25e-06\n",      "max u: 1.00000052356\n",      "t = 6.82: error = 9.17e-06\n",      "max u: 1.00000101698\n",      "t = 6.84: error = 9.12e-06\n",      "max u: 1.00000052153\n",      "t = 6.86: error = 9.04e-06\n",      "max u: 1.00000099668\n",      "t = 6.88: error = 9e-06\n",      "max u: 1.00000051949\n",      "t = 6.90: error = 8.91e-06\n",      "max u: 1.00000097665\n",      "t = 6.92: error = 8.87e-06\n",      "max u: 1.00000051744\n",      "t = 6.94: error = 8.79e-06\n",      "max u: 1.0000009569\n",      "t = 6.96: error = 8.75e-06\n",      "max u: 1.00000051537\n"     ]    },    {     "name": "stdout",     "output_type": "stream",     "text": [      "t = 6.98: error = 8.67e-06\n",      "max u: 1.00000093741\n",      "t = 7.00: error = 8.63e-06\n",      "max u: 1.00000051328\n",      "t = 7.02: error = 8.55e-06\n",      "max u: 1.00000091819\n",      "t = 7.04: error = 8.51e-06\n",      "max u: 1.00000051118\n",      "t = 7.06: error = 8.44e-06\n",      "max u: 1.00000089924\n",      "t = 7.08: error = 8.4e-06\n",      "max u: 1.00000050907\n",      "t = 7.10: error = 8.32e-06\n",      "max u: 1.00000088054\n",      "t = 7.12: error = 8.29e-06\n",      "max u: 1.00000050695\n",      "t = 7.14: error = 8.21e-06\n",      "max u: 1.00000086211\n",      "t = 7.16: error = 8.17e-06\n",      "max u: 1.00000050481\n",      "t = 7.18: error = 8.1e-06\n",      "max u: 1.00000084393\n",      "t = 7.20: error = 8.07e-06\n",      "max u: 1.00000050267\n",      "t = 7.22: error = 7.99e-06\n",      "max u: 1.00000082601\n",      "t = 7.24: error = 7.96e-06\n",      "max u: 1.00000050051\n",      "t = 7.26: error = 7.89e-06\n",      "max u: 1.00000080834\n",      "t = 7.28: error = 7.85e-06\n",      "max u: 1.00000049834\n",      "t = 7.30: error = 7.78e-06\n",      "max u: 1.00000079092\n",      "t = 7.32: error = 7.75e-06\n",      "max u: 1.00000049616\n",      "t = 7.34: error = 7.68e-06\n",      "max u: 1.00000077375\n",      "t = 7.36: error = 7.65e-06\n",      "max u: 1.00000049397\n",      "t = 7.38: error = 7.58e-06\n",      "max u: 1.00000075682\n",      "t = 7.40: error = 7.55e-06\n",      "max u: 1.00000049176\n",      "t = 7.42: error = 7.48e-06\n",      "max u: 1.00000074013\n",      "t = 7.44: error = 7.45e-06\n",      "max u: 1.00000048956\n",      "t = 7.46: error = 7.39e-06\n",      "max u: 1.00000072368\n",      "t = 7.48: error = 7.36e-06\n",      "max u: 1.00000048734\n",      "t = 7.50: error = 7.29e-06\n",      "max u: 1.00000070747\n",      "t = 7.52: error = 7.26e-06\n",      "max u: 1.00000048511\n",      "t = 7.54: error = 7.2e-06\n",      "max u: 1.0000006915\n",      "t = 7.56: error = 7.17e-06\n",      "max u: 1.00000048288\n",      "t = 7.58: error = 7.11e-06\n",      "max u: 1.00000067576\n",      "t = 7.60: error = 7.08e-06\n",      "max u: 1.00000048063\n",      "t = 7.62: error = 7.02e-06\n",      "max u: 1.00000066024\n",      "t = 7.64: error = 6.99e-06\n",      "max u: 1.00000047838\n",      "t = 7.66: error = 6.93e-06\n",      "max u: 1.00000064496\n",      "t = 7.68: error = 6.9e-06\n",      "max u: 1.00000047613\n",      "t = 7.70: error = 6.84e-06\n",      "max u: 1.0000006299\n",      "t = 7.72: error = 6.81e-06\n",      "max u: 1.00000047386\n",      "t = 7.74: error = 6.76e-06\n",      "max u: 1.00000061506\n",      "t = 7.76: error = 6.73e-06\n",      "max u: 1.00000047159\n",      "t = 7.78: error = 6.67e-06\n",      "max u: 1.00000060045\n",      "t = 7.80: error = 6.65e-06\n",      "max u: 1.00000046932\n",      "t = 7.82: error = 6.59e-06\n",      "max u: 1.00000058605\n",      "t = 7.84: error = 6.56e-06\n",      "max u: 1.00000046704\n",      "t = 7.86: error = 6.51e-06\n",      "max u: 1.00000057187\n",      "t = 7.88: error = 6.48e-06\n",      "max u: 1.00000046475\n",      "t = 7.90: error = 6.43e-06\n",      "max u: 1.0000005579\n",      "t = 7.92: error = 6.4e-06\n",      "max u: 1.00000046246\n",      "t = 7.94: error = 6.35e-06\n",      "max u: 1.00000054414\n",      "t = 7.96: error = 6.32e-06\n",      "max u: 1.00000046017\n",      "t = 7.98: error = 6.27e-06\n",      "max u: 1.00000053059\n",      "t = 8.00: error = 6.25e-06\n",      "max u: 1.00000045787\n",      "t = 8.02: error = 6.19e-06\n",      "max u: 1.00000051725\n",      "t = 8.04: error = 6.17e-06\n",      "max u: 1.00000045556\n",      "t = 8.06: error = 6.12e-06\n",      "max u: 1.00000050411\n",      "t = 8.08: error = 6.1e-06\n",      "max u: 1.00000045326\n",      "t = 8.10: error = 6.05e-06\n",      "max u: 1.00000049117\n",      "t = 8.12: error = 6.02e-06\n",      "max u: 1.00000045095\n",      "t = 8.14: error = 5.97e-06\n",      "max u: 1.00000047843\n",      "t = 8.16: error = 5.95e-06\n",      "max u: 1.00000044864\n",      "t = 8.18: error = 5.9e-06\n",      "max u: 1.00000046589\n",      "t = 8.20: error = 5.88e-06\n",      "max u: 1.00000044632\n",      "t = 8.22: error = 5.83e-06\n",      "max u: 1.00000045354\n",      "t = 8.24: error = 5.81e-06\n",      "max u: 1.00000044401\n",      "t = 8.26: error = 5.76e-06\n",      "max u: 1.00000044139\n",      "t = 8.28: error = 5.74e-06\n",      "max u: 1.00000044169\n",      "t = 8.30: error = 5.69e-06\n",      "max u: 1.00000042942\n",      "t = 8.32: error = 5.67e-06\n",      "max u: 1.00000043937\n",      "t = 8.34: error = 5.63e-06\n",      "max u: 1.00000041765\n",      "t = 8.36: error = 5.6e-06\n",      "max u: 1.00000043705\n",      "t = 8.38: error = 5.56e-06\n",      "max u: 1.00000040882\n",      "t = 8.40: error = 5.54e-06\n",      "max u: 1.00000043472\n",      "t = 8.42: error = 5.49e-06\n",      "max u: 1.00000040465\n",      "t = 8.44: error = 5.47e-06\n",      "max u: 1.0000004324\n",      "t = 8.46: error = 5.43e-06\n",      "max u: 1.00000040053\n",      "t = 8.48: error = 5.41e-06\n",      "max u: 1.00000043008\n",      "t = 8.50: error = 5.37e-06\n",      "max u: 1.00000039643\n",      "t = 8.52: error = 5.35e-06\n",      "max u: 1.00000042775\n",      "t = 8.54: error = 5.3e-06\n",      "max u: 1.00000039238\n",      "t = 8.56: error = 5.29e-06\n",      "max u: 1.00000042543\n",      "t = 8.58: error = 5.24e-06\n",      "max u: 1.00000038837\n",      "t = 8.60: error = 5.22e-06\n",      "max u: 1.00000042311\n",      "t = 8.62: error = 5.18e-06\n",      "max u: 1.0000003844\n",      "t = 8.64: error = 5.16e-06\n",      "max u: 1.00000042078\n",      "t = 8.66: error = 5.12e-06\n",      "max u: 1.00000038046\n",      "t = 8.68: error = 5.1e-06\n",      "max u: 1.00000041846\n",      "t = 8.70: error = 5.06e-06\n",      "max u: 1.00000037656\n",      "t = 8.72: error = 5.05e-06\n",      "max u: 1.00000041614\n",      "t = 8.74: error = 5.01e-06\n",      "max u: 1.0000003727\n",      "t = 8.76: error = 4.99e-06\n",      "max u: 1.00000041382\n",      "t = 8.78: error = 4.95e-06\n",      "max u: 1.00000036887\n",      "t = 8.80: error = 4.93e-06\n",      "max u: 1.00000041151\n",      "t = 8.82: error = 4.89e-06\n",      "max u: 1.00000036509\n",      "t = 8.84: error = 4.88e-06\n",      "max u: 1.00000040919\n",      "t = 8.86: error = 4.84e-06\n",      "max u: 1.00000036134\n",      "t = 8.88: error = 4.82e-06\n",      "max u: 1.00000040688\n",      "t = 8.90: error = 4.78e-06\n",      "max u: 1.00000035762\n",      "t = 8.92: error = 4.77e-06\n",      "max u: 1.00000040456\n",      "t = 8.94: error = 4.73e-06\n",      "max u: 1.00000035394\n",      "t = 8.96: error = 4.71e-06\n",      "max u: 1.00000040226\n",      "t = 8.98: error = 4.68e-06\n",      "max u: 1.0000003503\n",      "t = 9.00: error = 4.66e-06\n",      "max u: 1.00000039995\n",      "t = 9.02: error = 4.62e-06\n",      "max u: 1.00000034669\n",      "t = 9.04: error = 4.61e-06\n",      "max u: 1.00000039765\n",      "t = 9.06: error = 4.57e-06\n",      "max u: 1.00000034312\n",      "t = 9.08: error = 4.56e-06\n",      "max u: 1.00000039535\n",      "t = 9.10: error = 4.52e-06\n",      "max u: 1.00000033959\n",      "t = 9.12: error = 4.5e-06\n",      "max u: 1.00000039305\n",      "t = 9.14: error = 4.47e-06\n",      "max u: 1.00000033608\n",      "t = 9.16: error = 4.45e-06\n",      "max u: 1.00000039075\n",      "t = 9.18: error = 4.42e-06\n",      "max u: 1.00000033262\n",      "t = 9.20: error = 4.41e-06\n",      "max u: 1.00000038846\n",      "t = 9.22: error = 4.37e-06\n",      "max u: 1.00000032918\n",      "t = 9.24: error = 4.36e-06\n",      "max u: 1.00000038618\n",      "t = 9.26: error = 4.32e-06\n",      "max u: 1.00000032578\n",      "t = 9.28: error = 4.31e-06\n",      "max u: 1.0000003839\n",      "t = 9.30: error = 4.28e-06\n",      "max u: 1.00000032242\n",      "t = 9.32: error = 4.26e-06\n",      "max u: 1.00000038162\n",      "t = 9.34: error = 4.23e-06\n",      "max u: 1.00000031908\n",      "t = 9.36: error = 4.21e-06\n",      "max u: 1.00000037935\n",      "t = 9.38: error = 4.18e-06\n",      "max u: 1.00000031578\n",      "t = 9.40: error = 4.17e-06\n",      "max u: 1.00000037708\n",      "t = 9.42: error = 4.14e-06\n",      "max u: 1.00000031252\n",      "t = 9.44: error = 4.12e-06\n",      "max u: 1.00000037481\n",      "t = 9.46: error = 4.09e-06\n",      "max u: 1.00000030928\n",      "t = 9.48: error = 4.08e-06\n",      "max u: 1.00000037256\n",      "t = 9.50: error = 4.05e-06\n",      "max u: 1.00000030608\n",      "t = 9.52: error = 4.03e-06\n",      "max u: 1.0000003703\n",      "t = 9.54: error = 4e-06\n",      "max u: 1.00000030291\n",      "t = 9.56: error = 3.99e-06\n",      "max u: 1.00000036805\n",      "t = 9.58: error = 3.96e-06\n",      "max u: 1.00000029977\n",      "t = 9.60: error = 3.95e-06\n",      "max u: 1.00000036581\n",      "t = 9.62: error = 3.92e-06\n",      "max u: 1.00000029666\n",      "t = 9.64: error = 3.91e-06\n",      "max u: 1.00000036357\n",      "t = 9.66: error = 3.88e-06\n",      "max u: 1.00000029359\n",      "t = 9.68: error = 3.86e-06\n",      "max u: 1.00000036134\n",      "t = 9.70: error = 3.83e-06\n",      "max u: 1.00000029054\n",      "t = 9.72: error = 3.82e-06\n",      "max u: 1.00000035911\n",      "t = 9.74: error = 3.79e-06\n",      "max u: 1.00000028753\n",      "t = 9.76: error = 3.78e-06\n",      "max u: 1.00000035689\n",      "t = 9.78: error = 3.75e-06\n",      "max u: 1.00000028454\n",      "t = 9.80: error = 3.74e-06\n",      "max u: 1.00000035468\n",      "t = 9.82: error = 3.71e-06\n",      "max u: 1.00000028159\n",      "t = 9.84: error = 3.7e-06\n",      "max u: 1.00000035247\n",      "t = 9.86: error = 3.67e-06\n",      "max u: 1.00000027867\n",      "t = 9.88: error = 3.66e-06\n",      "max u: 1.00000035027\n",      "t = 9.90: error = 3.63e-06\n",      "max u: 1.00000027577\n",      "t = 9.92: error = 3.62e-06\n",      "max u: 1.00000034808\n",      "t = 9.94: error = 3.6e-06\n",      "max u: 1.00000027291\n",      "t = 9.96: error = 3.58e-06\n",      "max u: 1.00000034589\n",      "t = 9.98: error = 3.56e-06\n",      "max u: 1.00000027007\n",      "t = 10.00: error = 3.55e-06\n",      "max u: 1.00000034371\n"     ]    },    {     "ename": "NameError",     "evalue": "name 'interactive' is not defined",     "output_type": "error",     "traceback": [      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",      "\u001b[0;32m<ipython-input-50-00a50a4c75fe>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m    126\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    127\u001b[0m \u001b[0;31m# Hold plot\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 128\u001b[0;31m \u001b[0minteractive\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",      "\u001b[0;31mNameError\u001b[0m: name 'interactive' is not defined"     ]    },    {     "data": {      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQIAAAD8CAYAAACcoKqNAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJztnXl8VOW9/9/PTPYdkrCEhCRAwiLKFnYFXKhLLS5IK71ttVq1Wuvtr7e9tb/eX6+193a53tra1lu12orWatUuUjesCu6g7EtCQggEEsi+J7Of5/fHzJDJrGeAQPR+377GzJnzeeb7Oec88z1nhu9zHqW1RhCE/91YzrYBQRDOPpIIBEGQRCAIgiQCQRCQRCAIApIIBEFAEoEgCEgiEAQBSQSCIAAJZytwXl6eLikpOVvhBeF/Bdu2bWvTWufH0p21RFBSUsLWrVvPVnhB+F+BUqrejE6+GgiCIIlAEARJBIIgIIlAEAQkEQiCgIlEoJT6nVKqRSm1N8J6pZT6pVKqVim1Wyk19/TbFARhODFzRfA4cFmU9ZcDZb7HrcBvTt2WIAhnkph1BFrrt5VSJVEkVwFPaO89zzYrpXKUUuO11sdP1dxKyxoAftv+CAAvVFbh9Hi4cPIk5owfj9USmsdWWtbw206v/lB3J38/WMVFEyezrKiUrKTk8Pour16jeWDne8zMHcclRVMoyRoV2VP3wwC8erSaDvsAF0+Ywrz8IhIiefLpm2y9PH1wOxeNL2P5+MnkJKdG1QM8VPU+pZm5XDyhjMmZuSilwnvybcfGo3XU93ZxSdFkFowrIslqDR/Dt1+77HYe27aN5aWlXFhaSm5aWlj9Y82/PbH8xNvbGZOVzvJzJlE+Pi+ip8ePeLdj1856DtQ0sXhJGbPnlJCUHNr1VlrWsO7ggwDY+h0888AG5q2YzoJLZjJ6bHZYPcC6ww8B8OKftpCansyi5dOYMqMgxNOJGIe856rqPQ3s3lLHwgunMXvxFFJSk8Lr6/4HAJfLw1O/eJXzFpWx8OIZ5BdE7h/rDnrb/OO5DwFYcMk5TJ09EUuE/rGu1rvdh6saef/lHSy6dBZzLzyHtMyUqPvJ8Gieuv9lps0tYeGnzmXcxLwQvRlOR0HRBOBowHKD77WQRKCUuhXvVQMTJ040HeDXH2wBYH9rKwc7OqhqbeW6mTO5vLxsyMH2H4Rf7vgAgH6Xkw2HD1DZ0UptVzs3n1tBZphkcEvOrZz71o0AvN14iLcbD1HV0cLN58xnxugx4T3tfQ+Agz3tVHY2UdXVzOrS81hVcg6WMB3w15XvAuD0eHjpaCX7u1o40NPKV6YuCpsMbsm+jVkffBGA95sP8VpjNfu7mrmxfD6zcieE9fTLne8DcKS3m+0tjVR1tHB1XzefLTs3bNL89RbvfjUMg1dqatjX0kJNWxu3VlSQl54eor957C0s+MPNAGyra6R7wE71sVbWnj+bismFYT39YZ13u9tae9m54zB1tS0ca+zkyqvmkpAQmqD++PNXvU+0ZvOG3dTsrOdw1TFW334J+RNCP3gAf3x4IwBVu4/S3tzDoerjXPHZhVQsLRui8/ePP/7mTQC6O/rZ+k4NddXHaTjcxqp/WkxScmLo+//ytRPPP3yzispthzm0/xirb72QguLwH7w//nIDAAf3NHDscBt1VY1cdv1iFn1qZtg++8f/fhGAvu4BNr+yi0N7Gzh6oImrb7uY1IzQZHBiPwFbN1ay691q6iobufa2iymeOj6sp6horWM+gBJgb4R1LwLnByy/AVTEes958+bpeKlqadF9DodpfVNfrz7a02VabxiG3tl6TLs9HtNtqjtbdI/TblrfauvTh3s7TOu11npXW6N2xeHpQGeb7rQPmNZ32Wy6tr1dG4Zhus3eI03a6XKb1jccbdednX2m9bZ+u66rbIjL04HKRu2wO03rjx9t1+3N3ab1TodLH9hzJC5PB/c1aNuA+T7b0tCuW462m9a73R69f8dh7YnQP4Ct2sRnXGkTdzH2fTV4UWs9M8y6h4FNWuunfcvVwAod46tBRUWFlhJjQRhelFLbtNYVsXSn458P1wNf8v3rwSKgO1YSEARhZBHzNwKl1NPACiBPKdUA/DuQCKC1fgh4GbgCqAUGgC8Pl1lBEIYHM/9qsDbGeg187bQ5EgThjCOVhYIgSCIQBEESgSAISCIQBAFJBIIgMMITwbJV97Fs1X28e+gwTo+H95rrqO5uJloR1LJV9/GHDz6kvX+AI71dbGqow+52x4xR1dLqLbVt3E2bvTemr421dTjcbrZ31FHd0xjT07JV99HS10ezrYdNzfuxeZwx9XuON2FozauN+2iy9cT0tKG6lgGni22tjexoPYZhwtOxnh46B2xsqD5AvzO2p+2Nx/AYBht31nKsvTumpxcr99PrcFBT08SePUfxeIyYMY50dWEbcPLWP/bR32ePqf/oaANuw2Dz+wdoONoR09MLeyvptts5XNPEri0Hcbs8MWPUtXfgcrrZ9NJOejr7Y+rfP1yPy+Nh6/u11B9sidk/nt22k44BGw2H29j2fi1OZ+w+W9PahtvtYeOft9DVGr1/xMJUZeFwYKaycNmq+/C7cwKuqxPozexjYuYovnnOhVxROCNkYMkFq+4DwAD6lkB3nkFaSiJry2fxjdnnh4w1CIzRDSSu9eDK6GdBXinfPudypmaF1m37Y7gArrOg89sYn57FLZNXcnnBXCzKElav8RZaOL/YT1q6YlXRHO6cegk5SUMH+QR66gHUGrCnDzA3t4i7z72U80aHjjXwx3ADtiuS6UqxkZ+ezl3nLuHz5bNDxhqcv+o+lM9TLzCwOIFEq5UrZ0zlXy88P2SsQaCnPiBhTiquVhczJo7hm9cup6I8dKyB35MH0DOTSdoxwKiUVNauXcTqaytCxhoE76cs7SKt1cP5F03nK3etZMy4oQOPAj3ZgPSJiSTs6ae0JJ9b77iYBQsnR/XEVAtpGzsYZU3i6i+dz5qvLAsZaxDoqR8Y7bGRctjOghXTuOU7n6ZgYm5ETw4gtdBK0ofdTCzM5aa7VrLkwmlR+6ynXJG2uZccI5FPr5nP529bQWpa+D7rP3a5fe0k7epg7orp3Pafn6N4WsEJrdnKwrN2F2Oz+HfZU7+/icq+Rjxolo8rY3Ry6Ai5QP1PHljDmOx0Xm+o5ZKiKZTnhI6QC2zz5e8s4zOzZvKX4x9xTvYEKnJLSLaGDkAJjPH4ozdy3Gii121naf40cpMzo+rvvX81k8dm8/Lx3SwfO5VpWeOjegJ47anb+fvxXUzOzGNR/iRSE6J7+u0jX6LF3cux/l4uLpzM+PSssHp/Wvh/913DzHFjeG7XXlZMmcS548eGHTQVGOOlP9zGpm0HGTcqk0XTJpKWEjpqL1D/4ENfQDc7ODy/jSWLpzBuXE5U/d0/WcXCggm89PSHVCyZwrSZhVit4S9e/W3+8sQt7H6rjoxrU6hYMImMMAN1AvUPPLiW9HYP+885wsIV05lQEn7wkF//rR99mvMnFvPS4+8ze9FkZswpJiExdNBUYJs/Pf4Vat+tx3qlhYqlU8jKjt5n//tXnyPfnsj2GQdYuHwaEyflR+0ft3//YlZOncqG/3mDGf9vCucuLScx6eQ+0iP6ikAQhFPjTI41EAThY44kAkEQJBEIgiCJQBAEJBEIgoAkAkEQkEQgCAKSCARB4GOSCKbe+3NcHjfNtui17YH6qff+nBZbFx4jcm17sH7AbafHFbmOPLiN0+Wi2dYZl6d2ezduI3Jte7De7nHQ6TDvyeZ00jTQGbW2PcSTrQdHlPEYwXqHx0lHv3lP3f39tPR0xeWpa6APh91lWu/0uOjsjj5GJLhNa3d8+6nH3sfAQOSxD8FtXB4X7V3x9dm2ni6MOPpsn3OAvt4BUzGiMaIrC6fe+3MG3Wk0YMHDhcvT+VzZClaMnx5SgjnYxt/SQKG4+/q5fLZkCelBYw1C9RqFk5uvm8iaiRdRmhE61iC4jcLDnIoEbp59IRePnxVSohvek4vv/NMcVhctIztorEF4Ty6uX1X      "text/plain": [       "<Figure size 432x288 with 1 Axes>"      ]     },     "metadata": {      "needs_background": "light"     },     "output_type": "display_data"    }   ],   "source": [    "# %load ft07_navier_stokes_channel.py\n",    "\"\"\"\n",    "FEniCS tutorial demo program: Incompressible Navier-Stokes equations\n",    "for channel flow (Poisseuille) on the unit square using the\n",    "Incremental Pressure Correction Scheme (IPCS).\n",    "\n",    "  u' + u . nabla(u)) - div(sigma(u, p)) = f\n",    "                                 div(u) = 0\n",    "\"\"\"\n",    "\n",    "from __future__ import print_function\n",    "from fenics import *\n",    "import numpy as np\n",    "\n",    "T = 10.0           # final time\n",    "num_steps = 500    # number of time steps\n",    "dt = T / num_steps # time step size\n",    "mu = 1             # kinematic viscosity\n",    "rho = 1            # density\n",    "\n",    "# Create mesh and define function spaces\n",    "mesh = UnitSquareMesh(16, 16)\n",    "V = VectorFunctionSpace(mesh, 'P', 2)\n",    "Q = FunctionSpace(mesh, 'P', 1)\n",    "\n",    "# Define boundaries\n",    "inflow  = 'near(x[0], 0)'\n",    "outflow = 'near(x[0], 1)'\n",    "walls   = 'near(x[1], 0) || near(x[1], 1)'\n",    "\n",    "# Define boundary conditions\n",    "bcu_noslip  = DirichletBC(V, Constant((0, 0)), walls)\n",    "bcp_inflow  = DirichletBC(Q, Constant(8), inflow)\n",    "bcp_outflow = DirichletBC(Q, Constant(0), outflow)\n",    "bcu = [bcu_noslip]\n",    "bcp = [bcp_inflow, bcp_outflow]\n",    "\n",    "# Define trial and test functions\n",    "u = TrialFunction(V)\n",    "v = TestFunction(V)\n",    "p = TrialFunction(Q)\n",    "q = TestFunction(Q)\n",    "\n",    "# Define functions for solutions at previous and current time steps\n",    "u_n = Function(V)\n",    "u_  = Function(V)\n",    "p_n = Function(Q)\n",    "p_  = Function(Q)\n",    "\n",    "# Define expressions used in variational forms\n",    "U   = 0.5*(u_n + u)\n",    "n   = FacetNormal(mesh)\n",    "f   = Constant((0, 0))\n",    "k   = Constant(dt)\n",    "mu  = Constant(mu)\n",    "rho = Constant(rho)\n",    "\n",    "# Define strain-rate tensor\n",    "def epsilon(u):\n",    "    return sym(nabla_grad(u))\n",    "\n",    "# Define stress tensor\n",    "def sigma(u, p):\n",    "    return 2*mu*epsilon(u) - p*Identity(len(u))\n",    "\n",    "# Define variational problem for step 1\n",    "F1 = rho*dot((u - u_n) / k, v)*dx + \\\n",    "     rho*dot(dot(u_n, nabla_grad(u_n)), v)*dx \\\n",    "   + inner(sigma(U, p_n), epsilon(v))*dx \\\n",    "   + dot(p_n*n, v)*ds - dot(mu*nabla_grad(U)*n, v)*ds \\\n",    "   - dot(f, v)*dx\n",    "a1 = lhs(F1)\n",    "L1 = rhs(F1)\n",    "\n",    "# Define variational problem for step 2\n",    "a2 = dot(nabla_grad(p), nabla_grad(q))*dx\n",    "L2 = dot(nabla_grad(p_n), nabla_grad(q))*dx - (1/k)*div(u_)*q*dx\n",    "\n",    "# Define variational problem for step 3\n",    "a3 = dot(u, v)*dx\n",    "L3 = dot(u_, v)*dx - k*dot(nabla_grad(p_ - p_n), v)*dx\n",    "\n",    "# Assemble matrices\n",    "A1 = assemble(a1)\n",    "A2 = assemble(a2)\n",    "A3 = assemble(a3)\n",    "\n",    "# Apply boundary conditions to matrices\n",    "[bc.apply(A1) for bc in bcu]\n",    "[bc.apply(A2) for bc in bcp]\n",    "\n",    "# Time-stepping\n",    "t = 0\n",    "for n in range(num_steps):\n",    "\n",    "    # Update current time\n",    "    t += dt\n",    "\n",    "    # Step 1: Tentative velocity step\n",    "    b1 = assemble(L1)\n",    "    [bc.apply(b1) for bc in bcu]\n",    "    solve(A1, u_.vector(), b1)\n",    "\n",    "    # Step 2: Pressure correction step\n",    "    b2 = assemble(L2)\n",    "    [bc.apply(b2) for bc in bcp]\n",    "    solve(A2, p_.vector(), b2)\n",    "\n",    "    # Step 3: Velocity correction step\n",    "    b3 = assemble(L3)\n",    "    solve(A3, u_.vector(), b3)\n",    "\n",    "    # Plot solution\n",    "    plot(u_)\n",    "\n",    "    # Compute error\n",    "    u_e = Expression(('4*x[1]*(1.0 - x[1])', '0'), degree=2)\n",    "    u_e = interpolate(u_e, V)\n",    "    error = np.abs(u_e.vector().get_local() - u_.vector().get_local()).max()\n",    "    print('t = %.2f: error = %.3g' % (t, error))\n",    "    print('max u:', u_.vector().get_local().max())\n",    "\n",    "    # Update previous solution\n",    "    u_n.assign(u_)\n",    "    p_n.assign(p_)\n",    "\n",    "# Hold plot\n",    "#interactive()\n"   ]  },  {   "cell_type": "code",   "execution_count": null,   "metadata": {},   "outputs": [],   "source": []  } ], "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.7"  } }, "nbformat": 4, "nbformat_minor": 2}