{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 30724 characters]",
      "text/plain": [
       "<Figure size 640x480 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 36572 characters]",
      "text/plain": [
       "<Figure size 640x480 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#8.5.1\n",
    "\n",
    "dt = .01\n",
    "p = 10\n",
    "D = 1\n",
    "tau = 5\n",
    "xmax = 100\n",
    "tmaxs = 50\n",
    "tmax = int(tmaxs / dt)\n",
    "\n",
    "deltaT = np.zeros(xmax)\n",
    "deltaR = np.zeros(xmax)\n",
    "T = np.zeros((xmax, tmax))\n",
    "R = np.zeros((xmax, tmax))\n",
    "\n",
    "\n",
    "S = np.zeros(xmax)\n",
    "S[int(xmax / 2) - 1] = 1\n",
    "\n",
    "for t in range(0, tmax):\n",
    "\n",
    "    if t > 0:\n",
    "        for x in range(0,xmax):\n",
    "            T[x,t] = T[x, t-1]\n",
    "            R[x,t] = R[x, t-1]\n",
    "\n",
    "    for x in range(0,xmax):\n",
    "        deltaT = S[x] + p*T[x,t]**2 / (1 + T[x,t]**2) - R[x,t] * T[x,t] / (0.1 + T[x,t]) - T[x,t] \n",
    "        if x > 0 and x < xmax - 1:\n",
    "            deltaT += D*(T[x-1,t] + T[x+1,t] - 2*T[x,t])\n",
    "\n",
    "        deltaR= T[x,t] - R[x,t] / tau\n",
    "\n",
    "        T[x,t] += deltaT* dt\n",
    "        R[x,t] += deltaR* dt\n",
    "\n",
    "plt.title(\"T\")\n",
    "plt.imshow(T, cmap='viridis', interpolation='nearest', extent=[0, tmaxs, xmax, 0])\n",
    "plt.xlabel(\"Time\")\n",
    "plt.ylabel(\"Position\")\n",
    "plt.colorbar()\n",
    "plt.show()\n",
    "\n",
    "plt.title(\"R\")\n",
    "plt.imshow(R, cmap='viridis', interpolation='nearest', extent=[0, tmaxs, xmax, 0])\n",
    "plt.xlabel(\"Time\")\n",
    "plt.ylabel(\"Position\")\n",
    "plt.colorbar()\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 21592 characters]",
      "text/plain": [
       "<Figure size 640x480 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 18600 characters]",
      "text/plain": [
       "<Figure size 640x480 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#8.5.6\n",
    "\n",
    "dt = .001\n",
    "D = 1\n",
    "xmax = 20\n",
    "tmaxs = 100\n",
    "tmax = int(tmaxs / dt)\n",
    "\n",
    "a = np.zeros((xmax, tmax))\n",
    "b = np.zeros((xmax, tmax))\n",
    "\n",
    "a[:,0] = .1 # influences number of stripes\n",
    "b[:,0] = .1\n",
    "a[xmax - 1,0] = 1\n",
    "\n",
    "for t in range(0, tmax):\n",
    "\n",
    "    if t > 0:\n",
    "        for x in range(0,xmax):\n",
    "            a[x,t] = a[x, t-1]\n",
    "            b[x,t] = b[x, t-1]\n",
    "\n",
    "    for x in range(0,xmax):\n",
    "\n",
    "        delta_a = a[x,t] ** 2 / b[x,t] - a[x,t]\n",
    "        delta_b = a[x,t] ** 2  - b[x,t]\n",
    "        \n",
    "        if x > 0 and x < xmax - 1:\n",
    "            delta_b += D*(b[x-1,t] + b[x+1,t] - 2*b[x,t])\n",
    "\n",
    "        if x == 0:\n",
    "            delta_b += D*(b[x+1,t] - b[x,t])\n",
    "\n",
    "        if x == xmax - 1:\n",
    "            delta_b += D*(b[x-1,t] - b[x,t])\n",
    "            \n",
    "        a[x,t] += delta_a * dt\n",
    "        b[x,t] += delta_b * dt\n",
    "\n",
    "plt.title(\"a\")\n",
    "plt.imshow(a, cmap='viridis', interpolation='nearest', extent=[0, tmaxs, xmax, 0])\n",
    "plt.xlabel(\"Time\")\n",
    "plt.ylabel(\"Position\")\n",
    "plt.colorbar()\n",
    "plt.show()\n",
    "\n",
    "plt.title(\"b\")\n",
    "plt.imshow(b, cmap='viridis', interpolation='nearest', extent=[0, tmaxs, xmax, 0])\n",
    "plt.xlabel(\"Time\")\n",
    "plt.ylabel(\"Position\")\n",
    "plt.colorbar()\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "venv",
   "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.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
