parameters.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from functools import partial\n",
    "\n",
    "import numpy as np\n",
    "from scipy.integrate import solve_ivp\n",
    "import seaborn as sns\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "sns.set_context('paper', font_scale=1.4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The wavelength of a Turing pattern is $2\\pi/k$ where \n",
    "$$ k^2 = \\frac{-d_u d_v(f_u - g_v) + (d_u + d_v)\\sqrt{-d_u d_v f_v g_u}}{d_u d_v(d_v - d_u)}$$\n",
    "\n",
    "The linearized substrate depletion model (Miura & Maini, 2004) is\n",
    "$$ \n",
    "\\begin{eqnarray*}\n",
    "\\dot u & = & 0.5u + 0.6v - u^3 + 2\\tilde\\Delta u \\\\\n",
    "\\dot v & = & -1.9u - 2v + 25\\tilde\\Delta v \n",
    "\\end{eqnarray*}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "23.996077217466947"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fu = 0.5; fv = 0.6; du = 2; gu = -1.9; gv = -2; dv = 25;\n",
    "\n",
    "k = np.sqrt((-du*dv*(fu - gv) + (du + dv)*np.sqrt(-du*dv*fv*gu)) / (du*dv*(dv - du)))\n",
    "2*np.pi/k"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.arange(0, 101, 1)\n",
    "\n",
    "\n",
    "D = 0.1\n",
    "d = 0.01\n",
    "\n",
    "def source(coupling, _, c):\n",
    "    dcdt = D*coupling(c) - d*c\n",
    "    dcdt[0] = 0\n",
    "    return dcdt\n",
    "    \n",
    "    \n",
    "offset = 1\n",
    "\n",
    "def turing(coupling, _, uv):\n",
    "#     assert (uv > 0).all()\n",
    "    u = uv[:101]\n",
    "    v = uv[101:]\n",
    "\n",
    "    dudt = 0.5*(u - offset) + 0.6*(v - offset) - (u - offset)**3 + 2*coupling(u)\n",
    "    dvdt = -1.9*(u - offset) - 2*(v - offset) + 25*coupling(v)\n",
    "    return np.append(dudt, dvdt)\n",
    "\n",
    "\n",
    "def test_diffusion_model(coupling):\n",
    "    \"\"\"Plot gradient from source and Turing pattern\"\"\"\n",
    "    plt.figure(figsize=(12, 3))\n",
    "\n",
    "    analytic_gradient = np.exp(-np.sqrt(d/D)*x)\n",
    "    c_0 = np.append(1, np.zeros(len(x) - 1))\n",
    "    c = solve_ivp(partial(source, coupling), [0, 1e3], c_0, 'BDF')\n",
    "\n",
    "    plt.subplot2grid(shape=(1, 4), loc=(0, 0))\n",
    "    plt.plot(x[:26], analytic_gradient[:26], 'k')\n",
    "    plt.plot(x[:26], c.y[:, -1][:26], 'ko')\n",
    "    plt.legend(['Exponential', coupling.description], frameon=False)\n",
    "    plt.xlabel('Space')\n",
    "    plt.xlim([-1, 26])\n",
    "    plt.ylabel('Concentration')\n",
    "    plt.yticks([])\n",
    "\n",
    "    uv_0 = offset + np.random.rand(len(x)*2)/100\n",
    "    uv = solve_ivp(partial(turing, coupling), [0, 1e3], uv_0, 'BDF')\n",
    "\n",
    "    plt.subplot2grid(shape=(1, 4), loc=(0, 1), colspan=3)\n",
    "    plt.plot(x, uv.y[:, -1][101:], 'ko')\n",
    "    plt.xlabel('Space')\n",
    "    plt.xlim([-1, 101])\n",
    "    plt.yticks([0, offset], ['', '$u_0$'])\n",
    "\n",
    "    sns.despine()\n",
    "    plt.tight_layout()\n",
    "    plt.savefig(coupling.description.replace(' ', '_').replace(\"'\", '').replace('.', '') + '.png', \n",
    "                dpi=300)\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA1cAAADQCAYAAADrn5sNAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3Xt8VPWd//H3BzKg4gVRugNq4hVRwcSCUlywumAl1EotFVcfq9Jt5VfxJ6s2uFZWe8H10gVXUuutbnVVWrVrq1INilYrtt5QEqGtyk9b4kqiCFq13Ab4/P6YSTqZzCQzyZmcubyej8f3kcn3zJnznXOZ8/2c7/d8j7m7AAAAAAC90y/sAgAAAABAKSC4AgAAAIAAEFwBAAAAQAAIrgAAAAAgAARXAAAAABAAgisAAAAACADBFQAAAAAEgOAKAAAAAAJAcAUAAAAAAch7cDVlyhSXRCIVcyo4qcdVfX29DxgwwM3MBwwY4PX19V3mk0jJqbW11ZcuXeqtra3teT3Zp9J9Dqn8EvsTKd+J/YmUj9TNPpI9d89rGjNmjANFLu/HSa4p+bhat26dRyKRDj8CkUjEGxsb0+a3tLS4u3tLS4s3NDS0/4/ytGjRovb9JBKJ+KJFi3q0T6X7HJQf9icELfVcFfT+xLmw/KTb5ln85mRdRyO4AroXejCVmpKPq4aGhrRXYObPn582v6GhgYpLmcq2krJ48eKc9qnFixd3GcijPPTV/kRluHykO1fles7ran/iXFh+crkAlPIbk3UdjXuugCJXU1OjSCTSIS8Siei0005Lmz9s2DDV1dUpFotJkmKxmOrq6tTa2ipJam1t1dKlS9v/R2mor69XZWWlamtrVVlZqfr6ejU1NbXvB23a/s9ln0qeL/lzGhsb2Z9KWOq27Yv9Kd1+jNKRvE+1tLSkPVcNHz48kP3p17/+dZfnQpSeTPvU008/nfE3p0dyicR6kjK1XHHlCUUk9Jaq1JR6XKVeiamvr8+Yn+mqHy1apSvTVbmmpqasr+h2tU+1tLSk/Zz58+ezP5WodL8VmfaDoPan7ro6o7il7gvf/OY3sz5X9WR/ytSi2tDQ4O7UU0tB6jbMVP/JsvdF1nW0UIIrKnAoMqEHU6kp3XGV6USQmk/FpfRle0LpqpKS7nO6yk/9nOTAiv2ptHTVhSZf+1N3F4a6+nwUvnT7VEVFRY+6h2a7P+VyMYB6avEJ6gJQksINrrLs1wgUktCDqdTU23sZe1JxQXHI9YTiHlylNPlzqAiXjlyC9XTvD2q5VIZLV6Z96pvf/GZ3Fd6sZRt0UU8tfkFfAEoo3OCKChyKUOjBVGoKYqCYXCou6d6PwtPTE0o+UBEuDT0J1vuyPFSGi1Mu5598n3tyvXiAwpenC0CFG1yF+aMM9FDowVRqytconNn2Y6cyXJgKrbWIinBxK6RgPRmV4eKW6XwS5j6VjAuNxSfXi8U9VLjBlXv8AOrXr1/oBxCQpdCDqdSUz0ccZDtcNyeW8PXRCSXQMlIRLh6FFqxnQmW4eHR3PimUbcWFxuLRh8F6YQdX7u633367S/JVq1b19ssC+RZ6MJWa+vL5cVSGC1OhX/3NhIpw4SqGYD0TKsPFoZjOJ1xoLHx9HKwXfnD1hz/8wSX5U089FcQXRo6qqqp8xIgRXl1d3Z5WrlwZdrGy8qc//clvueWWDnm1tbX++uuvdztvVVVVT75n6MFUaurL4IrKcOEplqu/mVARLjzFGqwnozJceIo5YE9VTIFhuejjbVL4wdW2bdt8wIABnDxD0sMgoyA8/fTTXl1d3aN5Ca56hspwYSmFkzwV4cJR7MF6JqVwnBSzUgjYk3GhsfD0cbCedR3N3F35NHbsWF+xYkXaadXV1Ro3bpxuv/32vJah0Pz1r3/V66+/nrfPHzlypAYNGtTlew488EA99NBDqqmpac9bs2aNJk2apKefflqHHHKI6uvr9ctf/lJPPvmk7rnnHt1zzz3aZ599tGbNGu2666666667NGLECEnSDTfcoJ/85CcyMx1++OG65ZZbNHToUN11111avHix/u7v/k5NTU0aOHCgHnjgAR188MGSpHvvvVc33XSTtm3bpt13310//OEPVV1d3eV8I0eO1Nq1a3X44Ydr+PDheuyxxzp8nxtuuEH33Xeftm3bpkgkovr6eo0fPz7j986C5bQB+kBXx1W+tLa2qrGxUTU1NYpGo2ppaVFVVVWHp5pHIhE1NzcrGo32adnKTWtrqyorK0tq3S9dulS1tbWd8hsaGjRlypQQSlQ+SnXdl+JxUiy6Oz+knk+KRX19verq6hSLxRSJRLRw4UJddNFFnfIXLFigOXPmhF3ckpNuv8m0TfIg+7pgLpFYT1JXV9jPPvtsHz9+fFARZdFYsWJF2qtpQaUVK1Z0W4Z03QI//PBDf+CBB/yYY47xZ5991isrK9uj/zvvvNMHDBjgq1evdvf4PXPjxo1zd/elS5f6YYcd5hs3bnR397q6Ov/Hf/zH9vn23HNPf/vtt93dfe7cuT5r1ix3d3/uuef8lFNO8S1btri7+7PPPutHHnlkt/Ola7lKbpF6//332/Off/55P/zww9O+Lweht1Slpr5uuUqHq8J9J9vnsxSzYu4uVOxKed334pk26IVSPj/Q6h6OrnrK9NGxnHUdrSK3oC1Yo0aN0pIlS+JNaFZwjQN5M3LkSOWz1WHkyJFZve/+++/v1IJzxhln6Omnn9akSZP0+OOPd7iiNH78eB111FGSpK997WuaPXu2PvnkEz3xxBM666yztPfee0uSZs+erWOPPbbDfAcddJAkady4cbrlllskSQ8//LBWrVqlcePGtb9348aN2rx5c5fzdWflypX693//d23YsEEVFRV64403tHnzZu26665ZzY/s1NTUKBKJdLoymWOrILqR6YronDlzNGPGjKK8+ptONBrVggULOl2BbPtexXqluxClrsvu1n0xy3Sc0NKQX6V8fohGox1adJuamjp8T0mKxWJqbGws6pbfQtLS0tJ+vErx9VtXV6cZM2a0/4YV0roOPbj65JNP9M4776iysjLMovSpQYMGacyYMWEXI63t27dr1apV2nfffbVu3boO07oKgLuatssuu7S/7t+/v7Zv3y4p3mp63nnn6Zprrslpvq5s27ZN06dP1zPPPKMxY8bo448/1l577aWtW7cSXAWMynD+FdsJpbeoCOdfuQTryVKPk+6OK+SunAL2VKUcSBaKYgtg+4W58FGjRkmSVq9eHWYxkOSKK65QZWWlnn/+eX3729/WH//4x/Zpv/vd7/T73/9eknTnnXdqzJgx2mOPPfSFL3xB9913nz766CNJ0q233qovfOEL3S5r2rRpWrx4sZqbmyVJO3fuzKpFb6+99tJf/vKXtNO2bNmirVu3ar/99pMk/fCHP+z289Bzc+bMUXNzsxoaGtTc3Nzez7m+vl6VlZWqra1VZWWl6uvrQy5pcerqhFKq2irCbZWwTBXh1tbWMItZlLpbl6nrvlSV43GVT5l+7zOdH0pNWyAZiUQkKe2FxqVLl/KblYPUddYWwCYr5AA21OCqqqpKgwYNIrgKyZlnnqmampr29PDDD2vJkiW67bbbVFVVpZtuukkzZszQpk2bJEknnniirrvuOh1zzDG68847dffdd0uSTjnlFM2aNUsTJkzQ6NGj9dZbb2nRokXdLn/ChAm6/vrrdfrpp6u6ulpHHXWU7rvvvm7nO/roozV69GiNGjVKU6dO7TBtzz331NVXX61x48bphBNO0MCBA3uwZpALKsP5U2wnlHygIhwc1mUcx1VwCNjjuNAYnHTrrLsAtuDkcoNWT1J3N94fd9xxfs455wRwnxny6c477/Rp06aFXYywhD6ARWoqhAEtMinlG5nzrRwGrshVKQ+20NdYl3/DQBfB4Pc+Mwa6yF2BPxYi6zpaqC1XUrxrIC1XKHVmdo+ZnZv0/zVmdlXi9Y/M7Dkz+7fwShgcrgr3TLl3rcmELjc9l7puiu7qbx7R0hAMfu8zo6U4d92ts6JpCc0lEutJ6u4K+w033OADBw707du3BxJWAnnQ6+NA0h8ljUr6/3FJX5R0rKT6RN5DkqJdfMYsSSskraisrOyzL98TXBXODVc4u5dNqx4Psf6bAhi2uOhwHHaP1vXs0VKcuwJfZ9nX+XJ5c09Sd8HVE0884ZL8jTfeCOSbA3nQ28BqD0mfSuqflPeBpGGSLpR0eiLvUklTs/nMQu4W2IbKcPboWpM7KsKZsW56huOwawTsueNCY9eKLFgvnuBq3bp1LskffPDBQL45kAe9Da4+L+mFpP9HSlqXeD1P0qTE629IOjubzyyG4CoVFb7MCvxqXUGiIpwZ66ZnujsOy7kyzO93z3GhMb0iDNazrveFfs9VNBrVkCFDtGrVqrCLAuTLUEkbJMnM+kv6jqRXEtM+krRX4vWeif9LEv3P/4Z7YXqPez0yY930TFfHYbnfi8Xvd88xom5npT7KZOjBlZkxqAVK3VOS9jGzZZJultRf0quJaS9JOjHx+gT9LegqOVT44hi4IhgMdPE3BOvBSXccUhnm9ztIBKplsA5yaebqScqm+9Ls2bN95MiRvW6vQ/aqqqp8xIgRXl1d7dXV1T5z5kx3d7/yyiv93nvv7XLeroZl//znP+9/+tOfsi7Hd77zHf+Xf/mXrN8fkrweI5J+JOk5SVdmO08xdgt0p/85XWuCV+5dboqwa03RKcdulkV2L0xRoRt40a6D4rnnyt395ptv9v79+/vmzZt7/c1LTb5OjlVVVb5y5coezUtwFX4q1uDKvbwrw+VYSetr5RTAltN3DVORVgR7jIA9/8rtQmOJBOtZ19FC7xYoxZ91tWPHDr3xxhthF6WghNHHe+bMmbrxxhslSdu2bdPcuXM1atQoVVdX66STTur0/tbWVo0fP1633XabJGnIkCHq37+/JOnqq6/WEUccoZqaGtXU1Gjt2rVdLnvVqlU64YQT9NnPflZHHnmkrr766vZp+++/v9atWydJOuusszRmzBhJ0vbt2zVkyBB99FHJ3qpUcsqp/3lqdy261uRfyXc3SVJO3zVMXXWzLLXup6V+L0yhKKfnrJVlV/hcIrGepGyusG/YsMElddsdrZzk+4pkarfA++67z93dzzvvPP/P//xPd3f/7ne/61/60pfaWxTff/99d/9by9Xq1at91KhR/thjj3X6/I0bN/pee+3lmzZtcnf3v/71r2lbJpNbrj7++GPfsmWLu7tv2rTJa2pq/Pnnn3d393PPPdfvuusu37lzpx966KE+evRo/+CDD/y5557zcePGBbJOuhB6S1VqKuaWq1Sl2pqT6epvEV6tKyqlPuJbcvnLrUUlbKn7Tim2uJfq73ExKMWW6BL7TsXVcjVkyBANHz6cQS2S9MUVyfvvv1+NjY1qbGzUmWee2Wn6r371K1188cXaZZddJElDhw5tn/b73/9ep556qu69917V1tZ2mnfPPffUYYcdprPOOku33XabNmzY0P45mWzevFnf+MY3dPTRR+tzn/uc1q5d2/59J0+erGXLlum1117TkUceqVNOOUVPPfWUnnzySZ188sm9WQ0IWXetOcV4Zbirq78lfbWuAJTyiG+p5X/ggQcYuKIPJbfalEqLO63rhaMUW6JL8TtloyCCK0mMGJii0H/g9ttvP0WjUT355JNpp/fv318vvPCCvvWtb+m9997T5z73OS1fvrzLz7ziiiv0mc98Rq+++qqampp04oknasuWLZLiwdVTTz2lZcuW6eSTT9bkyZP15JNP6sknn9TkyZMD/37oO6VYGe7uhELXmvwqxRHfMpV/xowZBOshKIVKY7rfV0aaDE8pXGgkWI8juCpQhfADd9ppp2nRokXtAc769evbpw0ePFjLli3To48+qquuuqrTvJ988onee+89TZw4UVdddZWOP/54rVy5ssvlffjhhxo2bJgqKir0xhtvaNmyZe3Thg0bpn333Vc/+tGPdPLJJ+uEE07QU089pdWrV2v8+PEBfWOEpdgrw5xQCk9qAFvsleGuyk+w3veKvSJM63rhKfYLjQTrSXLpQ9iTlO29IT/5yU9ckn/88ce96hBZavp6tMDke662bt3qdXV1fsQRR3h1dbVPmjTJ3TuOFrh582Y/9dRTfc6cOb5z5872z3nnnXd83LhxPmrUKB89erR/5Stf8Y8++qjT8pLvuXr11Vf9qKOO8rFjx/pXv/pVP+2009rL4u5+ySWX+AEHHND+/0knneS1tbUBrI1uhX6PVWoqpXuuMimWvv/cW1Uciu1erNTycH9V4cl0jBfDvVjF8vtajlKP/WK4b6m7Mhba72sPZV1HK5hK4EsvveSS/IUXXujVNwfyIPRgKjWVQ3BVDJXhMjmhlIxiqQwTsBePYqkIE6wXr+4C4TDOM6nLLJNgvfiCq08//dTNzO+4445efXMgD0IPplJTOQRX7oVXGS7TE0pJybUy3NcVl0IrD3JTiBVhgvXi1lUgHMa5MN0yyyRYL77gyt39kEMO8YsvvrjH3xrIk9CDqdRULsGVe+FUhsv4hFLSuqoM90XFhYC9tIRdES6U30sEK10gHMa27WqZZRCs5ye4kjRa0rmSZienrubJpRI4bdo0nzx5ci+/OxC40IOp1FROwVWqnlSGcz3J5FJBKYMTSknLVBlubGwMtOKS7v0E7KWpLyrC2e5PBOulI5cLMUEF8rle/CnxYD344ErSv0naIulFSU8npV93NV8ulcB58+Z5NBrt7ZcHghZ6MJWayjm4yrUyPH/+/IwnmaAqKCV+Qil56SrDPa24ZLtPEbCXtqAqwr3dn5qamgjWS1TQF4ayeUh2mV/8yUtw9Z6kmlw+3HOsBP7sZz9zSb5+/fpefHcgcKEHU6mpnIMr99wqwxUVFVl30aGCUr6yveG/q4pLLvvU4sWLCdjLSJj7E8F66QrqwlBq/ve//30u/nSWl+DqfyVFcvlwz7ESuGrVKpfkzzzzTC++OxC40IOp1FTuwZV7dpXh/v37pz3JLF68mAoKupVLxSXXfSrT+wmmSleY+xPBeunq7YWhdPmZzp1lfvEn6zpaLg8RvlrSfDOryGGenIwYMUIVFRU8TBhAt1IfnJruYYXf+9730j7oU1LaB7ImT09+f01NDQ/WLEPptnmmh8dKue1T//AP/1CeD9csY2HuTzxounRlcy5cuHChWlpa0u5TS5Ys6ZS/Y8cOVVR0rO5Hkh6Szf7UjWyjMEnrJcUkbZP0fnLqar5cr7CPGDHCp06dWm7RMApb6C1VqYmWq8y66zdeX1+f04hetFAhVZD7VJleAUYS9ifkQ7YtWpm6vafer8y5MD/dAj+fKXU1Xy6VwEWLFnm/fv069QUFQhZ6MJWaCK5yk81N4VRQkAv2KQSJ/Ql9IdtnR7bls591kHUdzTweOOXN2LFjfcWKFd2+r6WlRVVVVR2aJiORiJqbm2l2RNgs7AKkyva4QtdaW1vV2NiompoafmcQCPYpBIn9CUHLtE+xr3Ur67pgTsGVmf0fxZ9ztZ+kdyXdLel27+JDsq0ELl26VLW1tZ3yGxoaNGXKlKzLCOQBwRUAAED5yroumPXgFGY2T9JMSQsk/VnSQZK+JWmo4oNd9ErbTZ2pLVdtN88BAAAAQCHLZbTAr0ua4u63ufvj7n6rpCmJ/F5LHd2kX79+jJwEAAAAoGjkElztoXhXwGQtifxAtA1TOnXqVB144IEMdQwAAACgaOQSXP1G0q1mtrckmdlgSYskPRtkgaLRqM4++2y9/fbbeu+994L8aAAAAADIm1yCqwslVUn6wMz+ImmDpEMlzQ66UBMnTpQk/fa3vw36owEAAAAgL7IOrtz9PXc/SfEA62RJle7+D+7eGnShKisrdcABB2j58uVBfzQAAAAA5EXWowW2cff/lfS/eShLBxMnTtRzzz2X78UAAAAAQCC6bLkys7eTXq83s/fTpXwUbMKECVq5cqU+/fTTfHw8AAAAAASqu5arc5NefzWfBUk1YcIE7dixQy+88IImT57cl4sGAAAAgJx12XLl7sn98vZy99+kJgU4FHuyo446SoMHD6ZrIAAAAICikMtogfdmyL87iIKk6tevn/7+7/+e4AoAAABAUeg2uDKz3cxsUPyl7Zr4vy2NlhTLV+EmTpyo559/XrFY3hYBAAAAAIHIpuXqU0kfS9ot8fqTpPSqpJvzVbgJEyZo06ZNamxszNciAAAAACAQ2QzFfpAkk/RbSccn5buk9e6+OR8Fk6SxY8dq4MCBWr58uY499th8LQYAAAAAeq3blit3X+vuf3b3/RKv21JzPgMrSRo4cKCOO+447rsCAAAAUPByeoiwmdVKOlHSvoq3ZkmS3P2fgy3W30yYMEF33HGH3F1m1v0MAAAAABCCrEcLNLPLJS2WNETSWZJ2SPqy8jighRQf1GL9+vV6880387kYAAAAAOiVXIZi/z+SJrv7+ZK2Jf5Ok7RPXkqWMH78eJkZXQMBAAAAFLRcgqt93P3VxOvtZtbf3ZdLmpyHcrUbPHiwRo8eTXAFAAAAoKDlElw1m9khiddvSpphZl+QtCn4YnU0ceJELV++PN+LAQAAAIAeyyW4ul7SYYnX8yXdIelRSd8LulCpJkyYoLfeekstLS35XhQAAAAA9EhWwZXFh+n7jaRlkuTuDZL2lrS3u9+Wv+LFTZgwQZLoGggAAACgYOXScvVHxR8cLEly923u/mnwReps//3314EHHkhwBQAAAKBgZRVcubtL+oOkyvwWJ7MJEyYQXAEAAAAoWLm0XP1C0qNmdoGZfdHMpralfBUu2cSJE7Vy5Ur94he/UGtra18sEgAAAACyVpHDe2cl/l6Wku+SDg6mOJn9+c9/lrtr+vTpikQiWrBggebMmZPvxQIAAABAVize4y9/xo4d6ytWrOjVZ7S0tKiqqkqxWKw9LxKJqLm5WdFotLdFBLpjYRcgVRDHFQAAALKSdV0w626BZnZ/hvyfZvsZPdXU1NQhsJKkWCymxsbGfC8aAAAAALKSyz1XtRnyTwmiIF2pqalRJBLpkBeJRFRTU5PvRQMAAABAVrq958rMZideRpJetzlE0vrAS5UiGo1qwYIFqqurUywWU//+/bVw4UK6BAIAAAAoGNkMaHFG4m8k6bUUH8jiPUkzAy5TWnPmzNGMGTN0+umna8uWLbrooov6YrEAAAAAkJVugyt3P0mSzOw6d788/0XKLBqNavbs2Tr33HP1zjvv6IADDgizOAAAAADQLut7rtoCK4vbLTnlr3idfelLX1IkEtEvfvGLvlwsAAAAAHQpl9ECx5rZK5K2SvokkT5N/O0zgwcP1qRJk/Tggw/25WIBAAAAoEu5jBZ4p6QnJFUr/tDggyUdpD54gHCq6dOn67nnnlNra2tfLxoAAAAA0soluKqSdIW7/9Hd1yanfBUuk2nTpsnM9NBDD/X1ogEAAAAgrVyCq8clfS5fBcnF0KFDdcIJJ9A1EAAAAEDByGYo9jYbJD1mZkskdeiP5+6XBVqqLEyfPl0XX3yxNmzYoH322aevFw8AAAAAHeTScjVQ0kOSdkgampL63Omnn64dO3ZoyZIlYSweAAAAADrIuuXK3b+Wz4Lkar/99tP48eP14IMPaubMmWEXBwAAAECZy6XlSmZ2uJldaWY3Jf1/dH6K1r3p06friSee0McffxxWEQAAAABAUm7PuTpd0guSDpF0biJ7L0k35KFcWfnKV76ibdu26dFHHw2rCAAAAAAgKbeWq6slneLuMxW/70qSGiWF1nJ10EEH6bOf/SyjBgIAAAAIXS7B1TBJLydee+Lvjhw/I3DTp09XQ0ODNm3aFGYxAAAAAJS5XAKj1ZJOTck7VVJTcMXJ3fTp07Vp0yZde+21am1t7X4GAAAAAMiDXIKrb0m6y8zul7Srmd0t6TZJc/NSsiw9/vjjkqSrr75alZWVqq+vD7M4AAAAAMpU1sGVu78saZSkFZL+S/GWrDHu/mqeytatlpYW1dXVtf8fi8VUV1dHCxYAAACAPpf1c64kyd1bJP1HnsqSs6amJsVisQ55sVhMjY2NmjJlSkilAgAAAFCOchmK/RdmNjEl7wQz+3nwxcpOTU2NIpFIh7xIJKKampqQSgQAAACgXOVyz9UJkp5PyfudpJOCK05uotGoFixY0CHAuuCCCxSNRsMqEgAAAIAylUtwtV3SwJS8XSTtDK44uZszZ46am5v16KOP6pBDDtHatWvDLA4AAACAMpVLcPW0pP8wswpJSvy9TtJv8lGwXESjUU2dOlVz587VI488ojVr1oRdJAAAAABlJpfg6lJJx0lab2avSVovabykf8lHwXrinHPO0ZAhQ7Ro0aKwiwIAAACgzOQyFHuL4sFVraT5kqZIOtbd1+WpbDnbbbfddMEFF+jOO+/Uxo0bwy4OAAAAgDKSS8uV3H2nu7/g7j939xfdPdT7rdKZPXu2YrGYfvzjH4ddFAAAAABlJJeh2A82s/8xs7fN7P3klM8C5mrYsGE6++yzVV9fr23btoVdHAAAAABlIpeWq3slbZV0oaQzUlJBueSSS7Ru3Tr9/OehPYILAAAAQJmpyOG9oyRNdPcd+SpMUKqrqzVp0iTdcMMNOvvss2VmYRcJAAAAQInLpeXqJUkj8lWQoF1yySV69dVXtXz58rCLAgAAAKAM5NJy9aKkBjO7W1Jr8gR3vznQUgWgtrZWhx9+uK655hpt2rRJNTU1ikajYRcLAAAAQInKpeXqeEl/kjRRHe+3+moeytVr/fr1U3V1tR5//HHV1taqsrJS9fX1YRcLAAAAQIkyd8/rAsaOHesrVqzI6zLSaWlpUVVVlWKxWHteJBJRc3MzLVjIVcHdtBfWcQUAAFCGsq4L5vScKzMbYGYnm9lMM5tsZgNyL1vfaGpq6hBYSVIsFlNjY2NIJQIAAABQyrK+58rMjpD0mKQBkpolHSApZmZfdPc/5Kl8PVZTU6NIJNKp5aqmpibEUgEAAAAoVbm0XN0k6ceS9nf38YoHV7dJ+lE+CtZb0WhUCxYsUCQSac8744wz6BIIAAAAIC9yCa6OkfQDT9yklfi7QFJ1PgqNBJEiAAAQBUlEQVQWhDlz5qi5uVkNDQ2aNm2annjiCW3cuDHsYgEAAAAoQbkEV+9JGpOSd4yk94MrTvCi0aimTJmim2++WVu3btUVV1wRdpEAAAAAlKBcgqvrJD1uZovM7BIzu1HS0kR+wRs+fLjmz5+v22+/XS+++GLYxQEAAABQYrIOrtz9vyV9WdJukk6WNEjSV9z9rvwULXgXXnihjj76aF1wwQXasWNH2MUBAAAAUEK6Da7M7Jtm9pAkufsz7n6+u0919/MlXWxm38h7KQNSUVGhm2++WStXrtQtt9wSdnEAAAAAlJBsWq7+WdJVGaZdJen84IqTf8cff7y+/vWva968eXrttde0dOlStba2hl0sAAAAAEUum+DqYHd/Ld0Ed18l6ZBgi5R/1113nbZv365jjjlGtbW1qqysVH19fdjFAgAAAFDEsgmuKsxscLoJifysH0RcKGKxmLZu3aqdO3e2/19XV0cLFgAAAIAeyya4elnSmRmmnSFpRXDF6RtNTU2dBrSIxWJqbGwMqUQAAAAAil02rU7/LulhM9tD0v9IWidpuKSvSpqn+AiCRaWmpkaRSESxWKw9LxKJqKamJsRSAQAAAChm3bZcufszks6VdKGktyRtTvydLek8d/9NPguYD9FoVAsWLFAkEmnP23///bX33nuHWCoAAAAAxSyr51y5+8PufpCkIyRNlHSkux/s7o/ktXR5NGfOHDU3N6uhoUGPPPKI3n33Xc2aNUvuHnbRAAAAABShnAajcPc3Jb2Zp7L0uWg0qilTpkiS7rjjDp177rk68sgjdd5556mxsVE1NTWKRqMhlxIAAABAMciq5aocnHPOOfr2t7+tyy+/XAcccABDtAMAAADICcFVkgsvvFBmpu3bt0tiiHYAAAAA2SO4SrJq1apO91wxRDsAAACAbBBcJWkboj1ZRUUFQ7QDAAAA6BbBVZLUIdrNTDt27NBDDz0UcskAAAAAFDqCqxTJQ7Q3Nzfrggsu0AUXXKB//dd/1c6dO9Xa2qqlS5dyHxYAAACADnIair1cJA/RftNNN+mggw7S3LlztWzZMq1evVqxWEyRSEQLFizQnDlzQi4tAAAAgEJAy1U3zEx1dXW69dZbtXLlSsViMUmMJAgAAACgI4KrLFVVVXXKYyRBAAAAAG0IrrKUbiRBM9PgwYMliXuxAAAAgDJHcJWl1JEEKyoqtNdee+mkk07SaaedpsrKStXW1qqyslL19fUhlxYAAABAXyO4ykHySILvvPOO3nnnHf3TP/2TlixZwr1YAAAAQJkjuMpR20iC0WhUu+++u6ZPn97pPbFYTCtXrmz/ny6DAAAAQOljKPZearsXq63lqs28efO0efNmNTc367LLLmP4dgAAAKDE0XLVS6n3YkUiEc2aNUuDBg3S9OnTdckll9BlEAAAACgDBFcBSL4Xq7m5WbfddpuWL1+uhQsXdnpvLBbTyy+/LInuggAAAEApIbgKSPK9WG3OPvvsTsO3t+Ufd9xxOuCAAzKOMEjgBQAAABQXgqs8Stdl8PLLL9fXv/51vfzyy9q+fbukeGvWpZdeqoaGBm3fvl319fVdDu1O4AUAAAAUHnP3vC5g7NixvmLFirwuo9C1traqsbFRNTU1ikajWrp0qWpra9O+d9CgQdq0aZOSt0skElFzc7Oi0ajq6+tVV1eXdoCM1OUgMBZ2AVJxXAEAAPSZrOuCtFz1gdQug20jDCaLRCL65S9/qenTpys14I3FYjrttNN00UUX6dJLL007QEZPW7toBQMAAACCQXAVgnTdBRcuXKgvf/nLuv766zsFXv369dPuu++un/70p9qxY0eHabFYLOOohC0tLZKUMfDKR0AW5DwAAABAUXH3vKYxY8Y40mtpafGGhgZvaWnpkL9o0SKPRCIuySORiNfX17e/vy2/LfXr188PPfTQDnltqaKiwvfbbz83sw75/fv39yuuuML79+/fIT8SifiaNWt8x44dncqwaNGitGVry+9qWk/m6W4dZcoPep6EvB8nuSaOKwAAgD6TdR2NSmCByiXwShd09e/f37///e/7jBkz0gZeuSYz8zFjxnQK1Pr16+fnn3++X3jhhd6vX79OZbjyyis7BXEVFRV+zz33+OLFi72ioqLTtIaGBn/llVf8sssua58eiUR83rx5vmbNGr/yyis7rIPvfe97/u677/q6det8/vz5HaZdc801vmHDBr/22ms75F9//fX+8ccf+yeffOI/+MEPugzwvACCqdTEcQUAANBnsq6jMaBFEUo3cEXqQBcLFy7URRddpNbWVlVWVrZ3GZTi3RBffPFFjRs3rkN+RUWFbrzxRr3++uu66aabOi131KhRWr16daf84cOHa/v27Xr//ffz8G37XvIAIgkMaAEAAFC+GNCilKV7plbqg4wvuuii9vemu7/rmGOO6ZR/ww036MILL9S8efPSDrixePHitPmvvPKKmpqa0k579dVX0+avWLFCL730UqdpFRUVWrp0qRYtWpT2u8+cOTNt/ne/+1195zvfSTvtzDPPTJs/d+5czZ07t1N+LBZTY2Nj2nkAAACAjHJp5upJovtSYcj1PqRM931lyg96nnRdHSORiDc1NaXNb2lpCXSelPURejfA1MRxBQAA0GeyrqNRCURGYQ8mEXaAlyT0YCo1cVwBAAD0mazraNxzhYKW6cHIXT0wOch5ErjnCgAAoHxlXRckuAK6R3AFAABQvhjQAgAAAAD6EsEVAAAAAAQg790CzWy9pLUZJu8r6YO8FqDwsQ7iCnk9fODuU8IuhJnNkjQr8e/hkt5I87ZCXo/oHbZt6WLblia2a+li25auTNs267pg3oOrLhdutsLdx4ZWgALAOohjPQSD9Vi62Lali21bmtiupYttW7qC2LZ0CwQAAACAABBcAQAAAEAAwg6ubg95+YWAdRDHeggG67F0sW1LF9u2NLFdSxfbtnT1etuGes8VAAAAAJSKsFuuAAAAAKAkEFwBAAAAQAAIrgAAAAAgAARXAAAAABCAUIIrMzvEzJ4zszfNbKWZlcWD2MzsCjN7w8x2mtmXk/L3NbPHzGyNmf3RzLJ6AnQxMrNdzOyhxPdsNLPHzezgxLSyWQ9BK9djqhRxjJQ+M/uamXnbeYDtWvzMbICZ3ZjYhr83s8cS+WzbImdmU83slcS5dZWZnZvIZ9sWkZ7Uwc1sVzO718z+X2L6edksK6yWq1sl3eXuIyRdJmmxmVlIZelLT0maKunZlPxrJb3q7odJ+qqk/zazQX1duD50i7sf4e41kpZIuiORX27rIUjlekyVKo6REmVmVZLOl/RCUjbbtfhdI2mApMPd/ShJX0/ks22LmJn1k/RTSee4+zGSTpX0YzPbXWzbYtOTOvilkuTuh0r6vKTrzKyyuwX1eXBlZkMlfU7Sf0uSuy9LTBrT12Xpa+7+oru/lWbSWZJ+lHjP7yU1Sqrty7L1FXff4u6PJ2W9IOngxOuyWQ9BKudjqhRxjJSuREXtvyRdJGlr0iS2axEzs90kzZJ0ubvvlCR3b0lMZtsWt7aLlEMSfwdL2iBpm9i2RaWHdfDkaeskPSJpRnfLCqPlqlLSOnePJeWtTeSXHTPbR1Ik6YdYKq/1cZGkh1kPvcIxVdo4RkrHpZJ+6+6vtGWwXUvCoZI+lHS5mb1sZr8zs2ls2+Ln7jsknSHpITNbK2m5pHMl7SG2bdHL4hitkvTnDNMyqgiqgL1E96WOymJ9mNm3JY2QNEnSrune0rclKimsuxLAMVI6zOwoxbucTMzm7XkuDoIVUbzC9Za7X2FmIxWvhKe7B4dtW0TMrELSPElfcfdnzexYSQ9LOibd2/u0cMiXrrZjVts4jJarZknDzSySlFeVyC877r5BUszMhiVll/z6MLM6SdMl1br7pnJdDwHhmCpBHCMl5wTFt9caM/uz4l15b1e8iwnbtbitleSS7pEkd39d8a5FB4ltW+xqJA1392clyd1flvSupKPFti16WZxXmyUdmGFaRn0eXLn7ekkvSZopSWZ2suKR4CtdzFbq7pf0f6X2q5s1kpaGWqI8MrNLFe/HerK7f5Q0qazWQ1A4pkoPx0jpcfdb3H2Yux/o7gcqfi/dLHe/RWzXoubuH0h6XImWqkRFbZSkVWLbFrt3FL94eaQkmdmhincDfUNs21LR1Xa8T9KFiWnDJZ0m6efdfaC5e15K2uVCzQ5T/Ob7fSVtUvwE81KfF6SPmdm/SfqmpKGSPpG0RdJYSTsl3S3pEEk7JH3L3R8Nq5z5ZGb7K/5j9bbi60CStrv7WDP7jMpkPQStXI+pUsQxUh7M7BlJN7r7Q2zX4mdmByo+WMlnFD+n/8DdF7Nti5+ZnSXpCsW3az9J17r7T9m2xaUndfDEYDV3SBqneOv0Ne7+k26XFUZwBQAAAAClJqznXAEAAABASSG4AgAAAIAAEFwBAAAAQAAIrgAAAAAgAARXAAAAABAAgisAAAAACADBVcjM7HgzW25mH5nZh2a2wsymhl0uAACScb4CgO5VhF2AcmZme0p6TNKlij/ArL+k4xR/UBkAAAWB8xUAZIeWq3CNkDRQ0p3uvt3dt7r7cnd/zsxONLMPzOxSM3vPzNaZ2WVtM5rZQWb2ZOI9G83sITMbljR9LzO71cz+18z+krjauGti2mFm1pCY9y0zm933Xx0AUEQ4XwFAFgiuwvWmpL9K+qmZTTOzoSnTB0s6XFKVpFpJc83stMQ0k/QfkvaTdLDirZD1SfPeLWm4pM9KGiLpckk7zWw3SU9JekTSMElflHSFmX0h+K8HACgRnK8AIAvmTot+mMxshKTLJJ0saX9Jz0o6P/H6aUn7uPvGxHuvlXSQu/9jms+ZIOmX7j7UzKKSWiQNc/fWlPfNkDTX3Y9NyvuepEp3/1o+viMAoPhxvgKA7tFyFTJ3f9Pdv+HuVYpf0dsm6Z7E5E/bTlQJaxW/8icz+4yZ/SzRjeJjxfvC75t4X6Wkv6aeqBKqJFUnbkj+yMw+UrwP/d8F/+0AAKWC8xUAdI/gqoC4+1pJP5Q0KpG1u5ntnfSWSknvJl5fK2mApBp331NS8ohNzZIGJa4IpmqW9JK7D05Ke7g7Iz4BALLC+QoA0iO4CpGZjTSzuWZWaXGfkfQNSc8n3rJT0jVmtouZjU5M+1li2h6K93//S2K+eW2fm7gCuETSLWY21Mz6J4bQHSjpUUn7m9ksMxtoZhVmNtrM2rtdAACQjPMVAGSH4Cpcn0gaK+l3ideNkj6VdF5i+oeS1ijeveIJSf/p7g8npn1H0hGSPpL0jKRfpXz2eZI2SnpN0gbFrxz2c/dPFe8v/0VJ70haL+l2SXsG/u0AAKWC8xUAZIEBLQqUmZ0o6X/cfd/u3gsAQFg4XwHA39ByBQAAAAABILgCAAAAgADQLRAAAAAAAkDLFQAAAAAEgOAKAAAAAAJAcAUAAAAAASC4AgAAAIAAEFwBAAAAQAD+P8e7I4Hn4Y7eAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1a1b7978d0>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "def ficks_law(u):\n",
    "    \"\"\"Finite difference scheme for Fick's law\"\"\"\n",
    "    lapl = np.empty(len(u))\n",
    "    lapl[0] = u[-1] - 2*u[0] + u[1]\n",
    "    lapl[1:-1] = u[:-2] - 2*u[1:-1] + u[2:]\n",
    "    lapl[-1] = u[-2] - 2*u[-1] + u[0]\n",
    "    return lapl\n",
    "\n",
    "\n",
    "ficks_law.description = \"Fick's law\"\n",
    "test_diffusion_model(ficks_law)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [default]",
   "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.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}