From 1338b6c06d004cb912ba8a0c5aab98b9af27d52f Mon Sep 17 00:00:00 2001 From: Gerardo Marx Date: Tue, 12 Oct 2021 20:33:02 -0500 Subject: [PATCH] Channel Flow --- Channel-Flow.ipynb | 1369 +++++++++++++++++++++++++++++++++ ft07_navier_stokes_channel.py | 127 +++ 2 files changed, 1496 insertions(+) create mode 100644 Channel-Flow.ipynb create mode 100644 ft07_navier_stokes_channel.py diff --git a/Channel-Flow.ipynb b/Channel-Flow.ipynb new file mode 100644 index 0000000..fe7f4b1 --- /dev/null +++ b/Channel-Flow.ipynb @@ -0,0 +1,1369 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 29, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "[,\n", + " ]" + ] + }, + "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/h5sWrgkvTVgbVr1zqPf67VAddkrg5UVVU5z5G5OuCjwEnm6oBrMkuT33BD5nflLp3MiUEfBU4aGxtZv349ZWVlXs5ZiJgIdInQjuFYIvRRfMR3afJsrA5kSmDUqFHOc0CERKASsCMs9wmoBC6MSIhAJWCHSsCOsEkAIiAClYAdKgE7wigBCLkIVAJ2qATsCKsEIMQiUAnYoRKwI8wSgJCKQCVgh0rAjrBLAEIoApWAHSoBO6IgAQiZCFQCdqgE7IiKBCBEIlAJ2KESsCNKEoCQiEAlYIdKwI6oSQBCIAKVgB0qATuiKAHIcxGoBOxQCdgRVQlAHotAJWCHSsCOKEsA8lQEKgE7VAJ2RF0CkIciUAnYoRKwQyWQIK9EoBKwQyVgh0rgDHkjApWAHSoBO1QCg8npPQsDfEsAyIoENm7cqBIYgng8rhIYBnJeBM3NzRw6dMi7BHp7e71JAOCZZ57hxIkT3iTQ1dWVleIjfX193iQAsGHDBk6dOuVNArFYLCsSaGxszBsJQI6LIChwEo/Hueqqqzhw4IDT+EePHgUS79aFhYVUVFTwwQcfOM1RV5eoIN/T00NlZSW9vb3s2rXLaY69e/cCiaIX1dXVHD58mMOHDzvN0d/fn3o8Z84cPvzwQ6fxg63AT506RVlZGcXFxc7H6Z133gESxVSuu+46jhw5wpEjR5zmSC9wsmDBAvbt2+cl/u7du1m4cKGzuFrgRFHylMgVOHnkkUecx9+6dWvqXW3BggXcdtttTuP39fXxu9/9LlWias2aNZSWljrN0dnZybPPPptq+xin+vr61Lv1pEmTnG+nHY/H2bBhQ+rqb9myZcybN89pjmBOIHjT8zFOjY2NvPTSS6m2jxybN29OFeVxXQTGSgQishz4b0Ah8L+NMf90jmNWAb8CDPAXY8y9l9q5oMCJy8IgwerAhx9+SFVVFfv27WPGjBnOi49s3LiREydOUFlZSVNTE6WlpU5zdHR0sGnTJkpKShg7diz9/f1O48MZCVRWVtLc3MzcuXOd5ggKkra3tzN37lwOHDhAWVmZ0xzd3d1s2rSJoqIipk6dSmtrq/NxamxsZPPmzVRUVNDT08P06dOdn7OBBObNm8f+/fsZN26cs/hgsXwoIoXA40AtMB9YLSLzM46pBH4KfN0Y8zXgh0576YjMJUIfZamGY4lw4sSJznNko+5AelXibKwOVFRUOM/R2NjIc889R3l5Offdd5/zicFAAjt27OCGG25g2bJlTuMH2NxHsAg4aIxpMcb0Ac8Cd2Yc833gcWPMMQBjjNsZGAfofQL2ZOM+gTCUJs+UwOjRo53GP5cEfKxogZ0IpgEfp7Vbk8+lMweYIyJvi8j25EeJsxCRh0SkXkTq02dXfaMSsEclYEeYJADu7iwcAVQC3wRWA/8uImd9SDLGPGGMqTbGVPuoEXcuVAL2qATsCJsEwE4EbcCMtPb05HPptAJ1xph+Y8whoJGEGIYVlYA9KgE7wigBsBPB+0CliMwWkSLgHqAu45jnSVwNICKlJD4qtDjs5wWjErBHJWBHWCUAFiIwxsSAh4HNwD5gvTFmj4g8KiIrk4dtBo6KyF7gdeDHxpijvjo9FCoBe1QCdoRZAmB5H4Ex5kXgxYznfpH22AA/Sv4bVlQC9qgE7Ai7BCCPvoZsg0rAHpWAHVGQAIRIBCoBe1QCdkRFAhASEagE7FEJ2BElCUAIRKASsEclYEfUJAB5LgKVgD0qATuiKAHIYxGoBOxRCdgRVQlAnopAJWCPSsCOKEsA8lAEKgF7VAJ2RF0CkGciUAnYoxKwQyWQIG9EoBKwRyVgh0rgDHkhApWAPSoBO1QCg8nZzUsDsiGBvr6+rEigrq5OJTAE3d3d/OlPf8prCQB5JQHIcREYYzh48CAHDx70WnzktddeQ0S8SQDgueeeo6SkxJsEurq6vFcgeu+99wC8SQBg06ZNFBcXe5NALBbzLoGglkG+SAByuK6BMYZHH3001Z43b57zAW1qaiIWi6XaVVVVTuMDgwpclJeXe9loND3HnDlznG80GovFaGpqSrV9j1NJSYkXCaTnmD17tpcKROk5fIzToUOHOH36NBCRugbpgiotLU1VJXJFX1/fIAlMmTKFrq4upzky92UcGBjwmmP8+PEcO3bMafx4PD5o7H2MU3d396B2cXGx13EqKirixIkTnDhxwlsOH+N0+vTplARck7MiKCgooLy8nMsuu4w1a9Y4jR1MDAYnYG1tLYsWLXKaI5gYPHXqFJAoeOGj7sALL7wAJN5Ff/hDt7vIB3MCgQhqampYunSp0xzBxGBfXx8A69atY/bs2U5zBAVJBwYGAPjpT3/qNH4wMRiIoKqqilWrVjnNcfr0aZ5++ulUwZyVK1cO8T8ujJxeNSgoKPC6OrBixQqnsQPSVwduvvlmLznSJwbnzp1LUVGR0/jpE4PLli1j5MiRTuPD4NWBJUuWOI8Pg6sSX3vttc7nTjJXB3xsyhtIoKOjw0stDshxEbgmc4nQx2e44Vgi9Fl8REuTn59sFB9Jl8Ddd9/NnDlznOeACIlA7xOwQyVgRzbuE8iWBCAiIlAJ2KESsCNsEoAIiEAlYIdKwI4wSgBCLgKVgB0qATvCKgEIsQhUAnaoBOwIswQgpCJQCdihErAj7BKAEIpAJWCHSsCOKEgAQiYClYAdKgE7oiIBCJEIVAJ2qATsiJIEICQiUAnYoRKwI2oSgBCIQCVgh0rAjihKAPJcBCoBO1QCdkRVAmApAhFZLiIHROSgiPzkS477jogYERlyI4RLRSVgh0rAjihLACxEICKFwONALTAfWC0iZ/3Vicg44BFgh+tOZqISsEMlYEfUJQB2VwSLgIPGmBZjTB/wLHDnOY57DPhnwM8WKklUAnaoBOxQCSSwEcE04OO0dmvyuRQici0wwxjzgsO+nYVKwA6VgB0qgTNc8mShiBQA/wr8g8WxD4lIvYjUZ+7nNxQqATtUAnaoBAZjI4I2YEZae3ryuYBxwALgDRH5CFgM1J1rwtAY84QxptoYU30hWzqpBOxQCdihEjgbmzPxfaBSRGaTEMA9wL3Bi8aYbqA0aIvIG8B/NMacf6/yCyAbEujs7OSNN97wLoG3337ba92BsEjg/fffVwlkmSHPRmNMTEQeBjYDhcBvjDF7RORRoN4YU+erc729vbS3t3svPrJz506vxUcArxKIxWJ89tlnfPbZZ94k0N/fz969e70WHwHYvn07FRUVXiQAiasmnxLo7u7m008/paCgIG8kADlc4CQej/PYY4+l2q7/eIBBdQ2ykUNEnG80Go/HicfjqXYYxikbOcLwM0SiwElBwZnpixtvvNF5/M7OTpqbm73mSBfdpEmTmDdvntP4wbtbgI+fobu7mz179njN0djYOKgYiOscxhjefffdVLumpsbLx4GdO3em2j7Gqb29nY8++giAMWPGOI2dsyIAqKioYOzYsdx+++1O43Z0dLBr1y5GjBhBLBbzUuCkvr4+VbQDYO3atU4LnAQTgwGlpaXOxymYEyguLqa3t5eamhrnORobG9m+fXuq7brASTAnEFBYWOi8SEswJ1BQUEA8Hqeqqsr5OLW1tdHQ0ICIYIxxXgcir79rcDFkrg74IH11wEcRlczVAR9zJ5kTgz4KnKRXJb777rudx8+cGLzxxhu9TgyuWrXKS4GTtrY2nn76aUaPHs0DDzzgPD5ETASZEvBRkFSXCO3ILE2e78VHVq1axdy5c53nyJTA+PHjneeACIlA7xOwYzgk4Lo0ebaXCPNdAhAREagE7FAJ2BE2CUAERKASsEMlYEcYJQAhF4FKwA6VgB1hlQCEWAQqATtUAnaEWQIQUhGoBOxQCdgRdglACEWgErBDJWBHFCQAIROBSsAOlYAdUZEAhEgEKgE7VAJ2REkCEBIRqATsUAnYETUJQAhEoBKwQyVgRxQlAHkuApWAHSoBO6IqAchjEagE7FAJ2BFlCUCeikAlYIdKwI6oSwDyUAQqATtUAnaoBBLklQhUAnaoBOxQCZwhb0SgErBDJWCHSmAweSEClYAdYZAAoBIYBnJ681LIjgQaGho4cuSINwkAbNiwgU8++STvJbBjxw7vEujs7PQmgVgs5l0C7e3teSUByHERtLe3px4vWbKE48ePc/z4cWfxP//8cwCOHDnCyJEjWbRoEa2trc7iA7z11lsAfPLJJ8ydO5fy8vLUltSu2Lt3LwDHjh1jyZIl9Pf3O8/R39/P0aNHAbjpppvo7Ox0Gv+DDz4AEtvMl5WVMXfuXA4fPuw0xzvvvAMkzqvFixdTXFzsfJzSa3rW1tZy7Ngxjh075ix+8Dtoa2tj4cKFzuLmTYETRVEGE7kCJz62HW9oaGD37t0AjBo1yvl22sGcwKlTpwBYtmwZ5eXlTnN0d3fz/PPPp9o+xqmxsXFQcRDXOYwxbNmyhY6ODgCuvvpqp+90kJgTeO6551JtH+PU1tbGq6++6jXHjh072L9/PwDXXXed09g5KwI4U+Bk1qxZTuPW19eze/dupk2bRltbG7feeqvTHOkSKCsro7Ozk3nz5jktcNLd3U1dXR3FxcWMGjWKkSNHOh+nYE6goqKC9vZ2ampqnOYIVgc6OjpS43TNNdc4zZFefGTSpEl8/vnnzsepra2Nbdu2MWHCBE6ePMkVV1zhPMf27dvZv39/6nfheg4oL1YNXJK5OuCa9NWBpUuX5vXEYPrqgOsCJ+lLhIsWLWL58uVO48PZqwM+CpKmrw7cf//9XiYGt2/fzubNm6mqquI73/mO8/gQMRFkSsB1QdJMCdTU1DiND+FYIhwOCfheIrz//vudXvEFZErA9TkbEBkRZPM+AZXA+TmXBPL9PoF8lwBERAQqATtUAnaETQJgKQIRWS4iB0TkoIj85Byv/0hE9orIX0XkNRH5ivuuXhwqATtUAnaEUQJgIQIRKQQeB2qB+cBqEcksv7sLqDbGXA38Efgvrjt6MagE7FAJ2BFWCYDdFcEi4KAxpsUY0wc8C9yZfoAx5nVjzMlkczsw3W03LxyVgB0qATvCLAGwE8E04OO0dmvyufPxIPDSuV4QkYdEpF5E6tNvxXSNSsAOlYAdYZcAOJ4sFJH7gGrgX871ujHmCWNMtTGmesqUKS5Tp1AJ2KESsCMKEgC7OwvbgBlp7enJ5wYhIkuAnwO3GGN63XTvwlAJ2KESsCMqEgC7K4L3gUoRmS0iRcA9QF36ASKyEPg1sNIYc8R9N4dGJWCHSsCOKEkALERgjIkBDwObgX3AemPMHhF5VERWJg/7F+Ay4A8i8oGI1J0nnBdUAnaoBOyImgTA8ktHxpgXgRcznvtF2uMljvtljUrADpWAHVGUAOT5nYUqATtUAnZEVQKQxyJQCdihErAjyhKAPBWBSsAOlYAdUZcA5KEIVAJ2qATsUAkkyCsRqATsUAnYoRI4Q96IQCVgh0rADpXAYHJ6z8IA3xKIx+NZkcCmTZtUAkOgEhgecl4ETU1NNDU1eS0+snnzZgBvEgB48sknKSoq8iaBrq4u7xWItm/fDuBNAgDr16+noKDAmwSCAic+JbBv3z727duXNxKAHBdBeoGTiRMnsnXrVm/xRYSenh62bNniNEewFbgxhvLycvbs2cOePXuc5ggKnMTjccrKyti2bZvT+JAocBJQWFjIK6+84jR++pbpEydO5PDhw84LnAQ5ent7ueKKK3jvvfecxofBBU5KSkp47bXXnMYPztlXX33V6bbvOSuCeDyeeiwi7Nq1y2n8gYGBQTlGjBjBlxVcuRjS/3gg8UtMl4+PHEGtBlcYY4jFYql2Nsbp+PHj3nM0NTU5jZ+Zo6CggJ07dzqNH4vFCAoSnTx5coijL4ycFUFBQQGTJ0+mpKSEdevWOY09MDDAxo0bU++ktbW1LFq0yGmOYGIwKHf1yCOPOL8MbWxsZP369QwMDDBy5Eh+9rOfOY2fPicAsHjxYpYtW+Y0RzAnEAhy3bp1zJ4922mO9vZ2nnrqKXp7E1+KtakQdKEEcwIAV155JWvWrHEaPxaL8Yc//IHGxkYAVq5cOcT/uDByetWguLjYy+pAIIGbbrrJaeyA9NWB6uohq01dFIEEysrK+OpXv+p8P/3M0uQjR470OjHoa24mkMDo0aNZsGCBlzmm9InByZMnO68BkS4BX+dsTovANekSyPeqxIEE1q5dy6hRo5zGz5SA66sAyE7xkXQJPPDAA5SUlDjPkbk6kF6qzwXpElixYoW3N5bIiEAlYMe5JJCPS4SZEshGBSLXqwPZkgBERAQqATtUAvaESQIQARGoBOxQCdgTNglAyEWgErBDJWBPGCUAIRaBSsAOlYA9YZUAhFQEKgE7VAL2hFkCEEIRqATsUAnYE3YJQMhEoBKwQyVgTxQkACESgUrADpWAPVGRAIREBCoBO1QC9kRJAhACEagE7FAJ2BM1CUCei0AlYIdKwJ4oSgDyWAQqATtUAvZEVQKQpyJQCdihErAnyhKAPBSBSsAOlYA9UZcA5JkIVAJ2qATsUQkkyBsRqATsUAnYoxI4g5UIRGS5iBwQkYMi8pNzvF4sIs8lX98hIrNcdlIlYIdKwB6VwGCG3MBNRAqBx4HbgVbgfRGpM8bsTTvsQeCYMeZKEbkH+GfgbhcdDJMEtmzZ4k0CQGgksG3btryWwMDAQF5JAOx2MV4EHDTGtACIyLPAnUC6CO4EfpV8/Efgv4uImGDv5YtkYGCAlpYWAG8SAHjnnXcoLi72JgGAl156iYqKCm8S6Orqoqury5sE+vv72bVrl9fiI5DYr3/ChAneJBCLxbxLINhpOF8kACBD/a2KyHeB5caY/5BsrwVuMMY8nHbM7uQxrcl2c/KYrvPFra6uNl+2d308Huexxx5LtadMmWL1A10I6cUospFj3LhxXiSQnqO0tNS5BHp7ezl+/Hiq7XucioqKvEggPcekSZO8SOCzzz5LtX2MU3d3N319fVRWVnLvvfcOebyINBhjhrRRVusaiMhDwEMAM2fOHOrY1OP58+d761NwcvjKEcS/7LLLmDFjhtccVVVVXsqQDQwMpETge5wgURfAZ465c+d6KUNmjEmJwNc4jRkzhsOHD1NeXu40ro0I2oD0M3h68rlzHdMqIiOA8cDRzEDGmCeAJyBxRfBlSUXESyEKRVHOxmbV4H2gUkRmi0gRcA9Ql3FMHXB/8vF3ga2XOj+gKEr2GPKKwBgTE5GHgc1AIfAbY8weEXkUqDfG1AH/B3haRA4Cn5GQhaIoeYLVHIEx5kXgxYznfpH2+DRwl9uuKYqSLfLmzkJFUfyhIlAURUWgKIqKQFEUVASKomBxi7G3xCKfAoctDi0Fznur8jCTy30D7d+lkMt9A/v+fcUYM+S9zsMmAltEpN7mXunhIJf7Btq/SyGX+wbu+6cfDRRFURGItSETAAADGklEQVQoipIfInhiuDvwJeRy30D7dynkct/Acf9yfo5AURT/5MMVgaIonskZEQz3BqmX2LcficheEfmriLwmIl/JVt9s+pd23HdExIhI1mbDbfomIquS47dHRP5vtvpm0z8RmSkir4vIruTv944s9u03InIkuQPYuV4XEfm3ZN//KiLXXnQyY8yw/yPx9eZm4KtAEfAXYH7GMX8L/K/k43uA53Kob7cCY5KP/yZbfbPtX/K4ccCbwHagOlf6BlQCu4CJyfbluTR2JD6L/03y8Xzgoyz27xvAtcDu87x+B/ASIMBiYMfF5sqVK4LUBqnGmD4g2CA1nTuBJ5OP/wjcJj725bqIvhljXjfGnEw2t5PYxSlb2IwdwGMkdpc+nWN9+z7wuDHmGIAx5kiO9c8AJcnH44H2bHXOGPMmif09zsedwFMmwXZggohMvZhcuSKCacDHae3W5HPnPMYYEwO6gck50rd0HiRh6WwxZP+Sl4wzjDEvZLFfYDd2c4A5IvK2iGwXkeVZ651d/34F3CcirST25Pj77HTNigs9N89LVjcvDTsich9QDdwy3H0JEJEC4F+BB4a5K+djBImPB98kcSX1pohcZYz5fFh7dYbVwG+NMf9VRGpI7MS1wBgTH+6OuSRXrgguZINUvmyD1GHqGyKyBPg5sNIY05uFfgUM1b9xwALgDRH5iMRnybosTRjajF0rUGeM6TfGHAIaSYghG9j070FgPYAx5l1gFIn7/HMBq3PTimxNfAwxKTICaAFmc2bS5msZx/wdgycL1+dQ3xaSmHSqzMWxyzj+DbI3WWgzdsuBJ5OPS0lc6k7Oof69BDyQfFxFYo5Asvj7ncX5JwtXMHiy8L2LzpOtH8jiB76DxLtBM/Dz5HOPkniHhYSJ/wAcBN4DvppDfXsV6AQ+SP6ry6Wxyzg2ayKwHDsh8dFlL/AhcE8ujR2JlYK3k5L4AFiaxb79HvgE6Cdx5fQg8APgB2lj93iy7x9eyu9V7yxUFCVn5ggURRlGVASKoqgIFEVRESiKgopAURRUBIqioCJQFAUVgaIowP8HSys+AHlS1+kAAAAASUVORK5CYII=\n", + "text/plain": [ + "
" + ] + }, + "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\u001b[0m in \u001b[0;36m\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\u001b[0m in \u001b[0;36m\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/NMcVhctIztorEF4Ty6uX1XI9SUrmJ4TWtcfzhMk8LMbl3N5YehYg/AxPNy1ahbXly8lNy0jgt7fRqNw89mFBayZsZzzikLnqQj1ZABWfnbFUi6bPS9krEF4Tw7+ZcFMrl2wnLzR2SH6QfWgp+vLR3PNrAuZPX2KSU+Kn65YwBWLFpMUVKI7GCPQk51/Pq+caxddxLhxuRH0g60suLl6UgarZ13E/JlTTfZZg/9YOo8rF58fMtYgvCcHd0yfyOrFKyksHHtC+4mpLFQBzxQa8DAhM48J6aEz/oS2wdfGQmFaHilR6/SV7+HdueNS8xidFL5Of2hLr6dpE/IpTMuPWac/6CmFwtR80hJCJ1wJ9eRlfFoe+Snh6/SHogGDay4bTVF6rklPBpDIhIxcMpKi1+kP8ZSRy9hMs548XDU9nQk5eTHGDvj/r4FUxmXmkhlm9qXw7jyMS89lzCgznsA71MdC4eh8EiOMHRh8f38fTGRcVh5ZmaGTwAS38ffZgoxcxuWOjthnh7byHYtRY0hOCd9nQz3BuKx8crLDj3eJiZnJD4bjEc8EJ+U/uF97DI+2u81NXlH+g/u1y+3WAy5zE0uU/+B+Xf6D+7XT7dIuj7lJO8p/cL82DEPb3OZjOFyuuD25Dbd2uF2m22it44rRZ7Npm8tpatKOIZ5c5j0ZhqFtJiemKf/B/bqtu1vbHfF7cjrMe/J4PNpmj+9Y2J3OiBOJhGvjNjza4TDfZ90n4cnhcmq3O7InTucEJ8OBDDoShOHnE/PVQBCE4UcSgSAIkggEQZBEIAgCkggEQUASgSAISCIQBAFJBIIgYPJ25kqpy4AHACvwqNb6J0HrJwLrgByf5m6t9cunaq7suR+GeVUDivevWEt+emgtefg2HsBK9XXfxaISTOi9Zae7r/k6qQmhpapD2wwWZL116XWMzwwd/xDN0/7V38FqSTKh9273rqvvIC0xN2RteE+avy+/nKn580168m531ep/JcGSbELv9bT96tvITBwTw5Nf7x1vUHPd91BB8z9E81S5+lskWlJNe9p61U1kJ4XO/xDek3dcRvV1/xeLssbQD3rad+03SbKmm9B7PW35zJcYnVJswhOcbJ/dc81dpCRkh1kfnZhXBEopK/AgcDkwA1irlJoRJPs34Fmt9RzgeuB/4nZiGm+HXvLyH9jf+XqY0WPBywpvbtJMff7HuAybiRje3XLeX++jy3HIhB/vY/mG59ja/EJcnqb9+ac4PX0mPHm3e9bffkGbvdqkJwufeesVyp//D7QOHtEWzpN3u6f/+SfYPWZGzXk9zf3br2ka2G1S793u8uf/E0MHj3iM7GnGn3+Mzd1u2lPFC4/Q0PeRCb0F7/lQM/X5H+HRwSMeI3s65y//SZ+rybSnhX//PQc6N5nUD/ZZt+Ew0cbr6dy//pRux1ET+qCIsUqMlVKLgXu01pf6lr8LoLX+cYDmYaBOa/1Tn/5nWusl0d73ZEqM/Zkw0pk6kv6ty65jfEboDDOR9EDYM3W0NpHO1JH0f19+GVPzK+LyFO5MHa3NjqtvIyPMmTpajHBn6mj6cGfqaG0inamjxQh3po6mD3emjtYm0pk6uqfQM3W0NmbP1H79u1dcz5i0KfH12evuxqqGDlIyW2JsJhFcB1ymtf6Kb/mLwEKt9Z0BmvHAa8AoIB24RGu9Ldr7ylgDQRh+zvRYg7XA41rrQuAK4EkV5rSilLpVKbVVKbW1tbX1NIUWBOFUMZMIGoGigOVC32uB3Aw8C6C1/gBIAUImk9NaP6K1rtBaV+Tn55+cY0EQTjtmEsFHQJlSqlQplYT3x8D1QZojwMUASqnpeBOBnPIF4WNCzESgtXYDdwIbgCq8/zqwTyl1r1JqlU/2L8AtSqldwNPAjfps3ehAEIS4MVVH4KsJeDnote8HPK8Elp5ea4IgnCmkslAQBEkEgiBIIhAEAZO/EYwUdnbU87MD955YfnLhY1H1DX3dtNr7uP/g903ptdZsbDzI7xu9RZMPn/cAaakZUdvs7jzKfTX34K8njxWjzd5Dw0AHPzvwQ1OeAN5uOshv6390YjlWm70dx/npgX8z7anLZqeqpZX/aboHgCcWPBqzom3L3np+3W/+WBw83Mo9zXeb9mQbcFK5+yi/UN7tXjf/t1gs0c9bO7Yd5n73YKVdrBhH6tv4XtN3THtyudzs3FLHLxN/6vP0CBZL9IrHPdvr+S+X+f10vKGDf238tmlPHo/Btjf38ausXwDw+4qHSLBGuwV6eEb0XYy/uOXmCGsG5x+AoTv3SxHbeNFoX6eymowRekAitxlkqP4mhs4iEIzB4/N/i9WSMKI8PTr31yQnpo8AT/4Y3oE1/jYn0z/MeApMhF/c8mXCXziH3+7h8hSYCOPxZLay8GN1RTBIaKcIXBOMd/d4d84TCx4zMcmEQeBgIjNn7cF2oZ7Cd259Yt0TC35nauKLQF+n3VM/kD7oad38x6KfgTvxFpQHbEdMT71AZhye2oC8wY69bv6j0c/AbcBowBL+AxeWLrxjZiN6CtoHQZ68Z+AoY1LafO+fELnPhtDh246IVwXRPcXXP7yM6CsCQRBODZnXQBAE00giEARBEoEgCJIIBEFAEoEgCEgiEAQBSQSCICCJQBAEPmaJwOnx3v563YFFrDuwKKbe4XGjtUZrhym91hqnx43WhukY8XpyxukJwGl426w7sIR1B6LeHBrwbjeAw9FjSu/0eDC0RmunKT2A0+XfbnOenE6v3jDspvRutwePx4jPk/PkPHmPRWy9x23gcRto7R45njwGbpcHrT2m+2A4RnRlYfSN8pejev3fULY5fBtf1eVQNSfaRNKfTIxT1cfvycINZe+HxtCB0jg8AQwpdR6Mse7Aosj125E8VS8aWjV8ogI2iicdTR/saQkEztfgb2v5X+gJwBL8JnBj+RapLDw1zNT+xy89qTc45fc/SwT7Hu7t8A8PGUlE9BQ84Uxsol4hnOp2ey9Tz/xj3rx5Ol4GXA6ttdaP1yzUj9csNKF3asMwtOEZMKU3DMPXxhNHjPg82dxO7TE82jDMefLHMAxDP16zWD9eszimvt/p1Fprbbe1m9LbXS7t9ni0YdhM6bXWesDmjWHW08CAdz953H2m9E6nW7tcbm0YdvOeBvzHwrwnb//oN6V3udza6XRrw3CajmGL05PthKcBU3q326Mddqc2DHfYPghs1SY+jyP6q4EgCKeGDDoSBME0kggEQZBEIAiCJAJBEJBEIAgCkggEQUASgSAISCIQBIERfjvzyqMFQ5b3G+D0wHlh5m+YUXQsbJvX2uHiXLAytNJ7UD9hiP4fPTA9CwpPKP0oZhQ1hsRow3tH7Clh/EeKUQ1MPbE0ODBg8P2H6j9wwMRkmBDV02CbY2440g2Lcs3HeHUfLJsKaQlDY0Tahu17ID8HiorMeerrgppDMHeOeU9vvQ7z50BarjlP+3dCWipMnGrOk+6AXdtg9jIg2Zyn916EWXMgY4I5T/Ufef8WzzfnCTts3QQVy4C0SJ6G9vEPn4Op0yF7JifNiK4sDN7gWMwoOkbV0YKoZdeBU2XMKDpmIkbgxBTegxGPr5OJEX50z2n2dKQgRn36WfBUXxDjGvUseDpc4D2LmPaECV+Dbc15GowRr6dzJh43VVk4ohNBMFtaDzMuNYui5GPsb/rUidf92TiYo33dNPZ1My8vG+14i9qOr0bVa635x7FqFuQXk2l8SHXrWt+awWwczL7uOjIT05iQ0El188UxPXU5Ozlub2RKeik4NnGg/ZaoeoC3mis5N2ci2eykuuWzMT3tam4iyWqlPLuH6qZl+IezRdJ3DdjY3dDEgpLxJHg2UtN2c0xP7+2qY2rxGEal7qW65dqYnmoOt2AYmrIiJ9VNi2J6GrA52b77CPPOKyKRjdS03RQzxofbDlFclEteTg3Vzati6uvrWhnod1A+LZH+vtdp6Pt2VL3L5WHLBweYW1FKsuVtqltviBlj54d1jC3IYeyYo+xvviLmdjce7aCjvZfpMzJx2TdxqPtOIPKx8HgMNr+yk1kXTCM1eTPVrV84sW5G0THTJcamEoFS6jLgAbx56FGt9U/CaD4L3IN3S3dprT8f7T1lrIEgDD+nbcozpZQVeBBYCTQAHyml1mutKwM0ZcB3gaVa606l1JiTty4IwpnGzL8aLABqtdZ1Wmsn8AxwVZDmFuBBrXUngNa65fTaFARhODGTCCYARwOWG/D/gD1IOVCulHpPKbXZ91UiBKXUrUqprUqpra2trSfnWBCE087pqiNIAMqAFcBa4LdKqZxgkdb6Ea11hda6Ij8//zSFFgThVDGTCBqBooDlQt9rgTQA67XWLq31IaAGb2IQBOFjgJlE8BFQppQqVUolAdcD64M0f8N7NYBSKg/vV4W60+hTEIRhJGYi0Fq7gTuBDUAV8KzWep9S6l6llP8fazcA7UqpSmAj8G2tdftwmRYE4fTysSooEgQhPj6R9yzc29pMt92OdtfjcVTSf7yY/uPFEfXH+3o52NmB1k60c/sJfaQ2Wms+aj2Cy/CgXXuwtT8XM0Z9/1F6XL1oTyOGa39MfZezh6MDjRiGG+38KKYngO0d9d6JUVyV2Hs3xtTXtLXR2t+P9jSb8tRts1N5vAXDMEx72nPgGHanC+3aj8v+Yez91NBOS1sv2ujAcMY+dgM2J1XVxzEMjXZuM+WpqrIR24AT7a7Fbd8VU3+soYOmY11ooxfDuSem3uXysG/XUe/EK84dpjzV7Gukv9eOdh/CY2K7m4910VjfjtZ2tHNnzBgej8HeLQfxuD1o5y5TnsIxoq8I+o8XRy7bjlQnb4BSEVYqUCi07001GgOwxBEjfXw9/cci7ORotfsRJwcJxQNYh9uTEWVdmNddQOIwe9IeUJFOTWHaOIGkePtHnPp4PZ1MjHiPheEGS6SxBkH6jIIjn4wrAocKeNh9jwg7LX18PU4L2JX2PgyN3a6xuzRaaYL/Sx9fjyVCDE+Ezg3gVIMPhwKHJ7onTZA+xjZYw3lygCuKp7D7KUKHTB9fj8cSvk2kD3VipBhRkp8znN4e2ZOyBug9sT0lRfLkiRzDdSqejNieiOQpSgy3JYI+QgxLgnlPZhnRVwSCIJwan8jfCARBGB4kEQiCIIlAEARJBIIgIIlAEAQkEQiCgCQCQRCQRCAIAh+zRGD4ip+0NtDawGgqj0OvMZrKo7bRWmNo40Qbw+00EWNQb86TX69P6OPZDsPjOqntjqo3/Po4PBkBngxPXHrAlF4Hbnd/Z1wxTB2LuD0ZQz05B4bVk5k+G+LJ7YgZIxwjurLQ3VSO7xh5CbRqDX9rdzeg/Dojtt4GJAW+4C9PtXorNoMzpWVcDa7AHR2gD/gzBFfw+wS0CafvB1LCbbfy1r2H9dQY4eBHiOHSYDHCrMBbxx5crdoPpEQo3Q0XI6onC1jDjR1wgTXCqSmcJweQEKcnZ0M5kYaiWMPsKIcbEiLow3lyAyoOTwBuN3F5sgOJEWIEe7KOP3B67mJ8NjEI6vQmaqmVJ2CnBuo1aBX6FlZABX4g/IKQ4L6Xm8rj0uN7OWwb/He5H0pC8Gsxtttoity5I223xYjc+cKhiKIPsxEn5UnF58kNJEYbyBN0PIymcixRroHDHQsr8XkyiJw4IvWRuD25QEUYdKR1fH79jOhEkDSuZshyh2OA0clpaKMbTRK0zAK8md5P4D7tddlJsSaSoDToAXRLRYg+JUCvtabbZSMnKQ1tdKA9Lmi/YIg+NEYvGQkZoPvQJMT0ZPPYSVBWEpQVdDe6ZWGIPjVoPwxudyfaY0D74qgxOm02clJSQPejSYSWc6PqBxwuLEqRnGgF3YluWRSizwjy1NlnY1RGqteTYYW20H0bGKOn10ZGegoKGygLuvm8EH3glZnT6cbjMUhNTfIei/5D0L92iD4zyFN31wDZOb79RCK0zI3qqbfHRlp6MhaLE7SBbpkd1ZPL5cbpcJOekeL15OqEzssj9qehnrp9nkJjBHrq67WTmpqExeoG7USH2YbkAL3HYzDQ5yAzO9XXZ23QfmGA3lxWGNFfDQRBODVk0JEgCKaRRCAIgiQCQRAkEQiCgCQCQRCQRCAIApIIBEFAEoEgCIzwysJg3mjaCUqxOMtJsmcn9t7/OLHOfyvpQA71trOp6QAXji2h2PIXBnp/OmR9cButNevq3mFmTiGzU3fj6v4Ag2eixtjeuR2H4eC89DRS3R9h6/1hVH2ro4tNLdtYMGoaRfwdW++Po3oCeLTyI6bm5LEoZz9u204M50NR9e8cPkxzXz8XFSeSozZh6/3PqPqufhvPvb+HZTNKKMt+AVvAfo3U5vnXdzI+P4v5kw+hPHtxDtwfVb+jqoEjxzo4f04WoxJfxdb346h6m83J83/ZyoL5kygr2ICt9/sxPb380k6ys9OYN6uZBL0XR/9Pouqr9zVSs6+RRRcUkJuxHlvfT6PqXU43zz/2NrMXT6F8yrvYe78b09Ob63eQmGSlYvEASZY92Pui99nDB1vY8WEdiy4oZlzOX2P2WY/H4Lnfv8PMOcVMn7oNe+83Q97TDCO6svCkJjiJdwILAgYpmWhzWicTiaA/I5OuxDmphltHOWucJU8OIPl09Y8IbeKd4CTe/jTcnsxOcDKirwiCs9/mtmpyktIoT+lDuw9h7749ohbgSF8HNT2tLMkvJNW9gYGu/xNVr7XmpcadzM0tYbxlF/bu7RieB0L06QWDzyt7KkmyJDEpGXAfwNb11agx2hxdVPUcZm5OGWkmPAG8UFfJ7PzxTEyqxG6rwbDdE9Im0NPWxkYMQzN3nMbi3hvTU1e/jfer6zl/ajGZltcY6Lorpqc3PqxhSlEeE3NrcTvrcfZ+K6qnyoNN9NsczC5Pwer+CFv316LGsNmcvPNeDQvnTyIr6U0GukL16UFt3n2nmgkTRlE8oRHDdRRH79ejxqiraaK9tZdZ83JJMN7B1h19u11ON5te3sW8peWMynyfgc5bo+oBPtxYxeixWUya0ol212PvviNqm6OH22g80s7s+QUk8zoDXd+Iqvd4DN58cSezF07vU4iDAAAbGUlEQVQmL2crA503BellrIEg/K9HxhoIgmAaU4lAKXWZUqpaKVWrlLo7im61UkorpWJmIEEQRg4xE4FSygo8CFwOzADWKqVmhNFlAv8MbDndJgVBGF7MXBEsAGq11nVaayfwDHBVGN0PgZ/ivZOSIAgfI8wkggnA0YDlBt9rJ1BKzQWKtNYvnUZvgiCcIU75x0KllAW4H/gXE9pblVJblVJbW1tbTzW0IAinCTOJoBEoClgu9L3mJxOYCWxSSh0GFgHrw/1gqLV+RGtdobWuyM/PP3nXgiCcVswkgo+AMqVUqVIqCbgeWO9fqbXu1lrnaa1LtNYlwGZgldZaigQE4WNCzESgtXYDdwIbgCrgWa31PqXUvUqpVcNtMJB2xxHchhPDXU/l0QIqjxZE1Xc4e+hwdKONPlyOvb42EyLqtdbU9Bz3Tizh2h9TD3CkvwW7x4nhbqDy6ISYbbqc/TTbu9DagdtRZSpGZXsLhtYYrhpT+sOdnfQ5nWjPcVOeeux2Gjq70dqNx1lpKsaBI614DAPDVWvqWBw/3kVfnx3taTXlyW5z0lDfjtYaj3Ofue2ubcHt9mC4D56IEY2W4130dA2gjQ56e9+JGcPlcnP4YIuvf5jbT/V1rTidbgz3YVP7qa2lh872PrTRg31gR8w2Ho9BXXWTz5O5/hSOEV1Z6N0g8/5mFB2LuaNPRQ+KGUWNwx4jnm0+aU/1BXH8QmTGU6AmgqdwN+kP9HS4YOgMIAETu5jyFKIP4yl4u4fbk5kY8Xo6VDB0cEAUT+dMPP7Jqyw0GDp5USAzio6Fb+OK1CbcLCg+vTu8fkZRY8irhieS/uQ8eWOE8eSKw5NfH3Y2HN97Bx55w0yMCJ5O9HsdpPfh8T5O7NcIMQDvB87te/j00Y+FGtS7A/Qn2gx+KE94sgToTXmKECOiJ4bqo8Q40T/CeNJhPflIMO/JLCP6ikAQhFNDxhoIgmAaSQSCIEgiEARBEoEgCEgiEAQBSQSCICCJQBAEJBEIgsDHLBH4i5/WHVgSl15rzboDi+LUxxfj5DyZb3MynhwOh6k2Z3K7DcMYcZ7Mxjg1T7H7X2iMePusuRjhGNGVhd6dHa4YVzOYw7zPbyh739cm0s5QAXq4oWxzeL2HgNpyFaI3HyOKpxMxIumDttvtk55o4z9mEWL4S4ujxgjQa18bhbek9sSLgfspjCeABP+xiODJ38T/NyRGBL0OeMuEods9xJMR4Mn/9Kx6WjS0y/qfW3XAZARR2kT0FEEfw9MnpLLQP+NF8CPQtiI0WYRrE7gu+MMZoLMGt/Hrl4TXR4wRxZM1vH4wRtB2JwS3GdpBQmJYI8cIq1eBMQLXEd1TQuCxUAG6AHsW5X0kRIoR9EH2662BMSIdC2NQbxkpngjvSQX32QBMeQros4H6qJ7MXyGM6CuCYGzu46RYx/JE7fmAMeQsHY4+dy8JKpEk+nmy7tMAUdtorWmydzEuJZsnapcQeGaIRLujg1FJOTxp0lOPy4ZCkW618eTBy03FaOjuZkJWFk/ULvVtQ3R9U28vuWlpuJ31PNvw+Zie+h1OXB6DrBQnTx78lClPx1u7GZfn96S5oeyDqPq2tl6ys9NQqpmnDl0bM4bd7sI+4CQ7R/HEwYtMeWo+1sWY8dk+T7GPRUdbLxmZKVgTOvhD3VUxY7hcbno6Bxidn8QTtStMeWo53kXe2CyePOhNFLE8dXX2k5KSSFJSD0/WXRmzjcdj0NHaQ97YDN/nYqje7BXBxyoRCIIQH5+QrwaCIJwJJBEIgiCJQBAESQSCICCJQBAEJBEIgoAkAkEQkEQgCAJD744+4rHbXuGW3c/jv/H7uvm/xWKJnMu6nUdptu3mh/v/jL/s8smFj0XUa635R/Mmnqz/A97aUitPLnw0qqdtHfv4xYH7Azw9gsVijahvGuhhU1MtL7X+2pQngGd27OYl5wM+T5aY+i11R/l16z0nPP2+4iESrIkR9V39Nl7ZUc3Lib/EX2sfK8bLb+7l6fSfn4gRS79vTwM/absH0gZLlaPtW9uAkw3rd/Dn0od8r8Q+Fm9u2MPvc35xwtMTCx5FKRVRf6DqGPce+j7kmvPkcnl48S9b+UvJw5jdT++8WcUj6f9t2tORQ618b9/dkD84uUG0GB6PwYtPfcDzU3+H2f4RjhFdWfjFLTdHWBM48GbojordRvsSiNWn/zLhL4wG9d6dO9hBzPga6ukmwt9j3q83eHz+b7FaEkxvQ/CHbzg8PTr31yQnpsfhafj3U2BCMNc/VNyeAj+ssT0NPRbD02cZctIz12e9nj6RJcZ226vcsvs5zF8RNNBs28MP9z+H2SuC15vf4on6JzF7RbC9s5Kf1/yMeK4I3m4+yN9bfmXKE8Cfdu7hRccvMJvxPzx0lF+13EM8VwQbdtbwYsIDvldin+FfeXMvf4zjiqBybwM/brsHUr1n39/M/DkZ6VkR9bYBJ6+9uJPni3/jeyX2sdj42l5+lz3oKdbZt7bqGD84/H0Ybf6K4OW/beP5iQ9h9org3Y1VPJwW5xVB1d2Qa/6K4OWnN/Ns2WOE6x+fyEQgCEJ8yFgDQRBMI4lAEARJBIIgSCIQBAGTiUApdZlSqlopVauUujvM+m8qpSqVUruVUm8opYpPv1VBEIaLmIlAKWUFHgQuB2YAa5VSM4JkO4AKrfV5wPPAf51uo4IgDB9mrggWALVa6zqttRN4BrgqUKC13qi1HvAtbgYKT69NQRCGEzOJYAJwNGC5wfdaJG4GXgm3Qil1q1Jqq1Jqa2trq3mXgiAMK6d1rIFS6gtABbA83Hqt9SPAI+AtKIr3/cue+yEANdd9D6Vi5zC/vnL1t0i0pJrWb7vqK2QljY/LU/V1/xeLilxRGKzfd+03SbKmm9Z/tOrL5CSbu9Aa9PRdLCr2Ifbr91xzFykJ2ab1H3z6S+Slmfs56GQ97b7mTlITRpnWv33pWsZlTo5avRfcZv/q72C1JJnW77r6dtIS82LqA9vUXPdvcXmqWv2vJFiSTeu3X30bmYljTHkKR8zKQqXUYuAerfWlvuXvAmitfxykuwT4FbBca90SK7CZykL/RobHw8srL2FK9tIhO7jsuXsJXxvuLdmsWv1tEiwpJmJ4a9u9H8ChnT1yGw/PLV/ErPxL4/IUnBRibffmz3yB3JQpQZ4ixfCWnQYnz1ie9lx7FynW7AB9NE8u3rvyesakDv3pKJan4KQQ/Vgodl9zB6kJuSY9uXnrimsoSJ8dl6f9192NVSUG6KN72nH1LWQkjjPpycOGlZcyKWexSU/+Pjs0KUT3pNm66iayk4uA01tZ+BFQppQqVUolAdcD6wMFSqk5wMPAKjNJ4NTxjgMoyZoXJssGLuuAv4qtV908JAlExj/1j5XspKK4PJ2Td0Fcnj5cdYOpKwMvbsDKqOSSMOsid+59134jzBVUZE/vX7F2SBKIjgdIJC+lPC5Pu6+5M8aVgR7y961LV5NiHR2HpwTGpgb/ph3Jk3fcwK6r7xiSBKJ707xw/krSE8aa9OTtHxOz5sThSbH96ttMXRkMDsqykJUU7Zt7pOa+efWiPYArgBrgIPA932v34v3gA7wONAM7fY/1sd5z3rx5Ol6mPHuvNgxjWPUej2fEeXK5XMPuKZ7t9nsaScdiyrP3apvNNqI8nWyM0+kJ2KpNfMZl0JEgfIKRQUeCIJhGEoEgCJIIBEGQRCAIApIIBEFAEoEgCEgiEAQBSQSCIDDCJziZeu/PgcA71LuBdl6963yKs1ZiCVN6GdhGoYEefnVzIUvGrCIjzKCM0BhOvnBlEqsnX8b0rGlhB4oMjeEB4L6bF3Dp2EUkJ4QOXhmqB7Dxtc9NY3XRCiakhZbNhnpyAYn87uYrWVIwxYQnb6npDz6zgM/MqCA9OdZ+ArBz65Jyrpu1jOL8fBOe3EACv7v2UhbPmI7FEsuTtwT2B8vncWXFIjLSQweBhcawc+vMYq6dfyGlRaGDwCJ5emjVxSw7dyZWa+h5LnS7nfzb0jmsmr+U7KzQUu/QGA5umD6e1RUXMbU0tPw8kqdfX3YBK2bPJjEp9CMXrn/cPX8WVy28gNGjQ2/5Hq7Pri3OZs3ClcyYOsnU4KZgRvwVgT4xYYNGYwXy+FPty9R3vUD4qkjta+NvlcnXH6thw7HHcXp6TcRI5A8vOniq/i8c6D0QUa9P6C1o4DfbX+cvjZswtBElhvevQTIP/mk3D9e+QIejz4SnBMDFY/tfZPOxOhOeFBqDf//7Rp7e9i4eI7InFeDpkfcreWTLS7R09ZjwZAXcPLH7Fd7fUxVBT4An0Cj+/a23+eNbb+B2e6J68j5L5pG9dTz27noaj7eb9vSHfa+ycfvuGJ78+yyB/3jvPZ7cuAGnw20iRhLrqhr5/QfrOXykyZQnjYcn977Gax9sjdBngz0l85OPtvD4Gy8x0G83ESORp+u7+P37f6Om5khYfSw+FiXG/gz46jcWU5wxG0uMgUN+/SvfuIz81BwyE6MPKfbr//qVq9HpA8zImhEzq/rbPPG1lZyXXUZqQvSBISdi3HU1mYmZFKVHH8bq1wOs+8qVLBg/BYtJT+u+fDkz80rISDW3n/5+x3VY3UlMHhd9AE2gp6fWXsGcSWVhz7rh2vzxhiuYMrqQ7MzoA6xOHLs71uDpV5QVRx9AE+jpibWXUTFpqmlPf7rx00zMKmB0ToYp/atf/xy2difTJxeb7x+f+xRzSstJSoo+mMmvf/7LVzIuJZ/8/BxT+pduX4OtzcbMaeGvFGWCE0EQZKyBIAjmkUQgCIIkAkEQJBEIgoAkAkEQkEQgCAKSCARBQBKBIAiM8ESwbNV9LFt1HxuqaxlwujDsb2E4d6EjlPEGtjnW00u/6whN/a/jNgZi6nccO47HMNjWsYkOZ/Q7si9bdR8vV9XQ53Cyt7uK6t7aiKXFgTGOdnXTYu/hteN76HWFLx0N1H90tAG3YfC36ioOd3XG9PTC3kp67HZ21Daytfoobk9sT4c6OunutfHa25X09Npi6j+oP4LL42HTxkrq69silsz62zy9bQedNhv79x9n69ZDuFzhS4sDY9S0tdE/4OCVN/bS1R372L1ddxiH280771Zz8GBLTE/LVt1HW38/B2ub2fJBbcTS4kB9ZXMLTqebDet30N4WvlQ9UP/GgYPYXW62vFvDgapjpjw19fZx9HArH2yswj7gjKnfc7wJl9vDa0+9R2tDR0S9GUZ0ZeGyVfedGFjhBK66rY6vnvsm6Sn5JGTejSXl0yFllResug/w3uE95VP9fP6C15mQPUBx5lqmjLqdREtGxBh9wHlfa2DcuEYmZczgMwU3UZBaEuLLH8MFZKx1MaroCHmpGVxfdC0X5C/CEjSHgF+vgQHAc0MfyWmaqwrn8bWpK8lJSovoqQfwfCqZfu1i1rhx3LPsImaPDS2Z9sdwA56FaTgbnIxOTeGrq5awetm5WC1DPZ2/6j7fOAOvp6TRFpL7YeWy6dzxxeWMzhlaChy8n7ImJGLdP8DkSWO48+srmTU7dMYjvycPoGcmk7RjgKzEZL7whSWsvraChARrWL32xUgbbSWx08OKJVO546YVjM0fOgAn0JMNSCtJxlo5wMTC0XztjktYsGBSVE/GuUmkbO0ngwQ+t3Yxn1u7mKTkhLB6v6dRCW6SjjpZeH45X/2XSykoHDpoLNCTA0ibmEDi1h4KxuZw6zc+xZIVoQPZAvusZ6oi7YNeMpwWrv78YtbeupzUtOSwMbxD6iCvt53EXe3MWT6d2398PcXTCk5ozVYWjujRhzA49cPzT9xKQfK7wGewJF+AsmRG1f/yfz7PlHwnHfbJjElbTnrixKgxvvWjK1k+uYSdvRsoSiujOH0a1ghTmPljPPP4V2ijHofhYFbOTDISwtfR+/X/9avrmZibweste1k2Zhol6XkRa9YVcOcPLmXltDL+Wl1JeW4eCwoKSbJG9/TE726itbmX9p5+zp9ZyqjMtLB6f1r40QNrmDYqj1ff3MeSeZOYNDG6py9+aynXzpvNe69XMe7GHGbPLg758AR7euzRG3Eet9GwtJNFiyaTmxu+tt+vv/f+1cwaO44XX9nFwnmTKJ80JqongJeeup3t7x4k+/o05s4tISUlfG2/X//ww1/E2umhdn4TixaXMWZs6Ci/QP2/3XcNi4oKeemZD5m7YBLlMyaEHXEZ2OavT95G9Qf1JF2bwNyFk0hLDz8exa//1W/+iewBxZ7Zh1m0fBrjCiNP9aaA7/zkM1xQUsyrj2zinB9PYfr8yTHHWUR8v5F8RSAIwqkhYw0EQTCNJAJBECQRCIIgiUAQBCQRCIKAJAJBEDCZCJRSlymlqpVStUqpu8OsT1ZK/cm3fotSquR0GxUEYfiImQiUUlbgQeByYAawVik1I0h2M9CptZ4C/Bz46ek2KgjC8GHmimABUKu1rtNaO4FngKuCNFcB63zPnwcuVidzc/UYaKMjas12MB6jB69lk++vNQPu7rg89bj64vLU57Lj9ESubQ9Hhy1yvX04OvtscXmy2V3YHa64YnT1xOepu9eGYZj3ZHe6GLCbP3Yn46mn14YnyniMYFwuD319kceIhKM7yliJcPT12iPe6j0cHo+Hno7IYx/MYqbEeAJwNGC5AVgYSaO1diuluoFcoO1UzK20rAHgtZ5LANCOd8BoRidfiEq+GJIWDCk99ev/1LUMAI/RRk//M6SlrCA95VNkpF6KJWiswUrLGq56uZSKCm/x1faO9WQnjmVy5kLKM5cyOnlCiB7g9iNfAaCy5wCNtmbmjZrJ/NHncW721JCxBista/js/rsA6HM7eKpuC0vHTGbF2KlcNH4a2UmpYWNc+9G3AHhm3x5GpaZyUckkLp9cTtno3LD6ta98G4B9h5vZW9/EBTNLWX7uJJacUxIy1mClZQ1f+rP3/Z0uD0+v38rcmUUsnTeZZQumMCo7LWyML/7V2+bVd6pITLRy/txJLJ9fRnnJmPD6v3n1tUfa2Fp5hKWzJ3H+3EksnT0pZKzBSssavvCCdxs8hsGTL33EeWUFnD9nMhdWTCEvJ/TYBXp6c0sNLrfB+XMnsaxiCudMGR9Wf+Mz/wJAY2Mn726uZdH8SSxZOJnFC6eQHDQByUrLGr781De9Cxqe+dNmysrGsXjxFJZdMJUxY7LCxvC32bz5IN3dAyxePIWlS8qYObMwbJ+96cn/A0Bray9vvL6P+QsmsXhJGUuWlpGamhQS4+Yfff7E8l8eeInC8gIWXTmPC1YvYvyk6LelD8cZHWuglLoVuBVg4sTItf8hJJR5/7r3gyULlVAOCaEzuvzDeI6VljUkJU4FwOVOIjFhMkmJU0lKLEep8GMBLlt6LeAdxJGRMJq85GLykovJTMwNqwcoSvN2suP2FizKwsS0AgpTx4UkAT9TsrwflBZbD8UZuUzJHMPkrHwyEyPPh1A22jv3wZj0dEpzRjF1dB4FGeHHWABMHu/129zZx+TxuUwpyGNyQW5IEvBTWuR9/95+OwVjsyktyqO0KJesjMjzIZQWetvkjkonb1QGpYV5FIzJjqnv6rVROiGXSUV5lBbmhSQBP5MKvdtgd7gZn5fFpMI8Jk3IJScjdGakEE/Vx0hPSWRSYR6F4yLPC1A60at3OtyUTMyltDiP0uJ8khLDeyop9urdHoMx+VmUluRRWpLPqFGR52jwt6mpbmL0qHRKivMoLBwd0mfXPLSS5776D4pLfbNLKUVR0WhKS/MpKc2POGai+BzvLEva0OQWjKbknCJKZk5k9Pjo8yFEIuZYA6XUYuAerfWlvuXvAmitfxyg2eDTfKCUSgCagHwd5c1PZqyB1gYqwgft9Oi9duP5VmNoI+KH/3TovW10zMlNhugNHXFAzOnQn4kYI9WTUnH2j2H3ZKCUijwo6zSONfgIKFNKlSqlkoDrgfVBmvXADb7n1wFvRksCJ0s8H+qT00feoZGI90Mdr97bJk5PcX6A4tWfiRgj1VPc/WPYPVlOaq7DYGJ+NfB9578T2ABYgd9prfcppe4Ftmqt1wOPAU8qpWqBDrzJQhCEjwmmfiPQWr8MvBz02vcDntuBNafXmiAIZwqpLBQEQRKBIAiSCARBQBKBIAhIIhAEgbN481KlVCtQb0KaxymWKg8jI9kbiL9TYSR7A/P+irXW+bFEZy0RmEUptdVMZdTZYCR7A/F3Koxkb3D6/clXA0EQJBEIgvDxSASPnG0DURjJ3kD8nQoj2RucZn8j/jcCQRCGn4/DFYEgCMPMiEkEI/kGqSa8fVMpVamU2q2UekMpFTo18Fn0F6BbrZTSSqkz9mu4GW9Kqc/69t8+pdQfz5Q3M/6UUhOVUhuVUjt8x/eKM+jtd0qpFqXU3gjrlVLqlz7vu5VSc086mNb6rD/wDm8+CEwCkoBdwIwgzR3AQ77n1wN/GkHeLgTSfM9vP1PezPrz6TKBt4HNQMVI8QaUATuAUb7lMSNp3+H9Ln677/kM4PAZ9LcMmAvsjbD+CuAVvJMjLwK2nGyskXJFMGJukHoy3rTWG7XW/rtUbgYKz4Av0/58/BDv3aXju/vm8Hu7BXhQa90JoLVuGWH+NOC/MWE2cOxMmdNav433/h6RuAp4QnvZDOQopcZH0UdkpCSCcDdInRBJo7V2A/4bpI4Eb4HcjDdLnyli+vNdMhZprV86g77A3L4rB8qVUu8ppTYrpS47Y+7M+bsH+IJSqgHvPTm+fmasmSLevhmRM3rz0k86SqkvABXA8rPtxY/y3q/tfuDGs2wlEgl4vx6swHsl9bZS6lytdddZdTXIWuBxrfXPfPfvfFIpNVNrbf4+6B8DRsoVQSNQFLBc6HstrMZ3g9RsoH2EeEMpdQnwPWCV1tpxBnz5ieUvE5gJbFJKHcb7XXL9GfrB0My+awDWa61dWutDQA3exHAmMOPvZuBZAK31B0AK3jr/kYCpvmmKM/XDR4wfRRKAOqCUwR9tzgnSfI2hPxY+O4K8zcH7o1PZSNx3QfpNnLkfC83su8uAdb7neXgvdXNHkL9XgBt9z6fj/Y1AncHjW0LkHws/zdAfCz886ThnaoNMbPAVeM8GB4Hv+V67F+8ZFryZ+DmgFvgQmDSCvL0ONAM7fY/1I2nfBWnPWCIwue8U3q8ulcAe4PqRtO/w/kvBe74ksRP41Bn09jRwHHDhvXK6Gfgq8NWAffegz/ueUzmuUlkoCMKI+Y1AEISziCQCQRAkEQiCIIlAEAQkEQiCgCQCQRCQRCAIApIIBEEA/j9MCH2JrH9c9QAAAABJRU5ErkJggg==\n", + "text/plain": [ + "
" + ] + }, + "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 +} diff --git a/ft07_navier_stokes_channel.py b/ft07_navier_stokes_channel.py new file mode 100644 index 0000000..a05cdea --- /dev/null +++ b/ft07_navier_stokes_channel.py @@ -0,0 +1,127 @@ +""" +FEniCS tutorial demo program: Incompressible Navier-Stokes equations +for channel flow (Poisseuille) on the unit square using the +Incremental Pressure Correction Scheme (IPCS). + + u' + u . nabla(u)) - div(sigma(u, p)) = f + div(u) = 0 +""" + +from __future__ import print_function +from fenics import * +import numpy as np + +T = 10.0 # final time +num_steps = 500 # number of time steps +dt = T / num_steps # time step size +mu = 1 # kinematic viscosity +rho = 1 # density + +# Create mesh and define function spaces +mesh = UnitSquareMesh(16, 16) +V = VectorFunctionSpace(mesh, 'P', 2) +Q = FunctionSpace(mesh, 'P', 1) + +# Define boundaries +inflow = 'near(x[0], 0)' +outflow = 'near(x[0], 1)' +walls = 'near(x[1], 0) || near(x[1], 1)' + +# Define boundary conditions +bcu_noslip = DirichletBC(V, Constant((0, 0)), walls) +bcp_inflow = DirichletBC(Q, Constant(8), inflow) +bcp_outflow = DirichletBC(Q, Constant(0), outflow) +bcu = [bcu_noslip] +bcp = [bcp_inflow, bcp_outflow] + +# Define trial and test functions +u = TrialFunction(V) +v = TestFunction(V) +p = TrialFunction(Q) +q = TestFunction(Q) + +# Define functions for solutions at previous and current time steps +u_n = Function(V) +u_ = Function(V) +p_n = Function(Q) +p_ = Function(Q) + +# Define expressions used in variational forms +U = 0.5*(u_n + u) +n = FacetNormal(mesh) +f = Constant((0, 0)) +k = Constant(dt) +mu = Constant(mu) +rho = Constant(rho) + +# Define strain-rate tensor +def epsilon(u): + return sym(nabla_grad(u)) + +# Define stress tensor +def sigma(u, p): + return 2*mu*epsilon(u) - p*Identity(len(u)) + +# Define variational problem for step 1 +F1 = rho*dot((u - u_n) / k, v)*dx + \ + rho*dot(dot(u_n, nabla_grad(u_n)), v)*dx \ + + inner(sigma(U, p_n), epsilon(v))*dx \ + + dot(p_n*n, v)*ds - dot(mu*nabla_grad(U)*n, v)*ds \ + - dot(f, v)*dx +a1 = lhs(F1) +L1 = rhs(F1) + +# Define variational problem for step 2 +a2 = dot(nabla_grad(p), nabla_grad(q))*dx +L2 = dot(nabla_grad(p_n), nabla_grad(q))*dx - (1/k)*div(u_)*q*dx + +# Define variational problem for step 3 +a3 = dot(u, v)*dx +L3 = dot(u_, v)*dx - k*dot(nabla_grad(p_ - p_n), v)*dx + +# Assemble matrices +A1 = assemble(a1) +A2 = assemble(a2) +A3 = assemble(a3) + +# Apply boundary conditions to matrices +[bc.apply(A1) for bc in bcu] +[bc.apply(A2) for bc in bcp] + +# Time-stepping +t = 0 +for n in range(num_steps): + + # Update current time + t += dt + + # Step 1: Tentative velocity step + b1 = assemble(L1) + [bc.apply(b1) for bc in bcu] + solve(A1, u_.vector(), b1) + + # Step 2: Pressure correction step + b2 = assemble(L2) + [bc.apply(b2) for bc in bcp] + solve(A2, p_.vector(), b2) + + # Step 3: Velocity correction step + b3 = assemble(L3) + solve(A3, u_.vector(), b3) + + # Plot solution + plot(u_) + + # Compute error + u_e = Expression(('4*x[1]*(1.0 - x[1])', '0'), degree=2) + u_e = interpolate(u_e, V) + error = np.abs(u_e.vector().array() - u_.vector().array()).max() + print('t = %.2f: error = %.3g' % (t, error)) + print('max u:', u_.vector().array().max()) + + # Update previous solution + u_n.assign(u_) + p_n.assign(p_) + +# Hold plot +interactive()