diff --git a/proto/FTCS.ipynb b/proto/FTCS.ipynb index 857b8e8..f47d159 100644 --- a/proto/FTCS.ipynb +++ b/proto/FTCS.ipynb @@ -2,17 +2,22 @@ "cells": [ { "cell_type": "code", - "execution_count": 1, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# imports\n", - "import numpy as np" + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "\n", + "%matplotlib notebook\n", + "%matplotlib inline\n", + "from ipywidgets import interact, IntSlider" ] }, { "cell_type": "code", - "execution_count": 25, + "execution_count": 113, "metadata": {}, "outputs": [], "source": [ @@ -21,6 +26,13 @@ "grid_size = {'x': 20, 'y': 20} # to start out\n", "domain_size = {'x': 20, 'y': 20} # to start out\n", "\n", + "# boundary_condition\n", + "# constant boundary condition\n", + "bc_left = 0\n", + "bc_right = 0\n", + "bc_top = 0\n", + "bc_bottom = 0\n", + "\n", "delta_x = domain_size['x'] / grid_size['x']\n", "delta_y = domain_size['y'] / grid_size['y'] \n", "\n", @@ -30,23 +42,48 @@ "max_stable_time_step_x = delta_x**2 / (2 * np.max(alpha_x))\n", "max_stable_time_step_y = delta_y**2 / (2 * np.max(alpha_y))\n", "max_stable_time_step = max_stable_time_step_x if max_stable_time_step_x <= max_stable_time_step_y else max_stable_time_step_y\n", + "max_stable_time_step = 0.25\n", "\n", "C_t = np.zeros((grid_size['x'],grid_size['y']))" ] }, { "cell_type": "code", - "execution_count": 26, + "execution_count": 114, "metadata": {}, - "outputs": [], + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 114, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQUAAAD4CAYAAADl7fPiAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy88F64QAAAACXBIWXMAAAsTAAALEwEAmpwYAAANiklEQVR4nO3df6xk5V3H8ffH5UcTBMuKbPllS+qGhDa6NpvFRjQgLQVC3NY0dYlRVAzYlMQmJgY1KU39p8Yg0dC02dYN1LRAo67dpMuPzWpCm7SUhSwFWpCV0LC3lLXdyhZbi9t+/eOea+5zd+7udc7MnbnT9yu5mXOe55lznpNJPnvOmdnzTVUhSQt+YtITkDRdDAVJDUNBUsNQkNQwFCQ1Tpr0BAY5JafWazht0tOQZtZ/81+8Wj/IoL6pDIXXcBqX5IpJT0OaWQ/X3mX7vHyQ1OgVCkmuSvJMkgNJbhnQf2qSe7v+h5O8oc/+JI3f0KGQZB3wEeBq4GLguiQXLxl2A/Cdqvo54HbgL4fdn6TV0edMYQtwoKqeq6pXgXuArUvGbAXu6pb/AbgiycCbG5KmQ59QOA94YdH6wa5t4JiqOgq8DPz0oI0luTHJviT7/ocf9JiWpD6m5kZjVW2vqs1VtflkTp30dKQfW31CYQ64YNH6+V3bwDFJTgJ+Cvh2j31KGrM+ofAIsDHJhUlOAbYBu5aM2QVc3y2/G/iX8v9qS1Nt6B8vVdXRJDcDDwDrgB1V9VSSDwH7qmoX8HfA3yc5ABxmPjgkTbFM4z/cZ2R9+YtGaXwerr0cqcMDvwmcmhuNkqaDoSCpYShIahgKkhqGgqSGoSCpYShIahgKkhqGgqSGoSCpYShIahgKkhqGgqSGoSCpYShIahgKkhqGgqSGoSCp0adC1AVJ/jXJV5M8leSPBoy5LMnLSfZ3fx/oN11J49an6vRR4I+r6rEkpwOPJtlTVV9dMu7zVXVtj/1IWkVDnylU1YtV9Vi3/F3gaxxbIUrSGjOSewpdNelfBB4e0P3WJI8nuS/Jm46zDcvGSVOgz+UDAEl+EvhH4P1VdWRJ92PA66vqlSTXAP8MbBy0naraDmyH+Ue8952XpOH0OlNIcjLzgfCpqvqnpf1VdaSqXumWdwMnJzmrzz4ljVefbx/CfAWor1XVXy8z5nULpeeTbOn2Zy1JaYr1uXz4ZeC3gSeS7O/a/gz4WYCq+hjz9SPfm+Qo8H1gm7UkpenWp5bkF4CBZacWjbkDuGPYfUhaff6iUVLDUJDUMBQkNQwFSQ1DQVLDUJDUMBQkNQwFSQ1DQVLDUJDUMBQkNQwFSQ1DQVLDUJDUMBQkNQwFSQ1DQVLDUJDU6B0KSZ5P8kRXFm7fgP4k+dskB5J8Jclb+u5T0vj0rvvQubyqvrVM39XM13rYCFwCfLR7lTSFVuPyYSvwyZr3JeC1Sc5Zhf1KGsIoQqGAB5M8muTGAf3nAS8sWj/IgJqTlo2TpsMoLh8uraq5JGcDe5I8XVUP/X83Ytk4aTr0PlOoqrnu9RCwE9iyZMgccMGi9fO7NklTqG8tydOSnL6wDFwJPLlk2C7gd7pvIX4JeLmqXuyzX0nj0/fyYQOwsysXeRLw6aq6P8kfwv+VjtsNXAMcAL4H/F7PfUoao16hUFXPAb8woP1ji5YLeF+f/UhaPf6iUVLDUJDUMBQkNQwFSQ1DQVLDUJDUMBQkNQwFSQ1DQVLDUJDUMBQkNQwFSQ1DQVLDUJDUMBQkNQwFSQ1DQVLDUJDUGDoUklzUlYpb+DuS5P1LxlyW5OVFYz7Qe8aSxmroZzRW1TPAJoAk65h/bPvOAUM/X1XXDrsfSatrVJcPVwD/XlVfH9H2JE3IqEJhG3D3Mn1vTfJ4kvuSvGm5DVg2TpoOmX8Ce48NJKcA3wDeVFUvLek7A/hRVb2S5Brgb6pq44m2eUbW1yW5ote8JC3v4drLkTqcQX2jOFO4GnhsaSAAVNWRqnqlW94NnJzkrBHsU9KYjCIUrmOZS4ckr0tXPirJlm5/3x7BPiWNSa8KUV39yLcDNy1qW1wy7t3Ae5McBb4PbKu+1yuSxqr3PYVx8J6CNF7jvqcgaYYYCpIahoKkhqEgqWEoSGoYCpIahoKkhqEgqWEoSGoYCpIahoKkhqEgqWEoSGoYCpIahoKkhqEgqdHryUvSgge+sX/FY99x7qaxzUP9eaYgqbGiUEiyI8mhJE8ualufZE+SZ7vXM5d57/XdmGeTXD+qiUsaj5WeKdwJXLWk7RZgb1fHYW+33kiyHrgVuATYAty6XHhImg4rCoWqegg4vKR5K3BXt3wX8M4Bb30HsKeqDlfVd4A9HBsukqZIn3sKG6rqxW75m8CGAWPOA15YtH6wa5M0pUZyo7Gr5dDrWfHWkpSmQ59QeCnJOQDd66EBY+aACxatn9+1HaOqtlfV5qrafDKn9piWpD76hMIuYOHbhOuBzw4Y8wBwZZIzuxuMV3ZtkqbUSr+SvBv4InBRkoNJbgA+DLw9ybPA27p1kmxO8gmAqjoM/AXwSPf3oa5N0pRa0S8aq+q6ZbqOqe1WVfuAP1i0vgPYMdTsJK06f+askfCny7PDnzlLahgKkhqGgqSGoSCpYShIahgKkhqGgqSGoSCpYShIahgKkhqGgqSGoSCpYShIahgKkhqGgqSGoSCpYShIahgKkhonDIVl6kj+VZKnk3wlyc4kr13mvc8neSLJ/iT7RjhvSWOykjOFOzm21Nse4M1V9fPAvwF/epz3X15Vm6pq83BTlLSaThgKg+pIVtWDVXW0W/0S80VeJM2AUdxT+H3gvmX6CngwyaNJbjzeRiwbJ02HXo94T/LnwFHgU8sMubSq5pKcDexJ8nR35nGMqtoObAc4I+t71aWUNLyhzxSS/C5wLfBbXYHZY1TVXPd6CNgJbBl2f5JWx1ChkOQq4E+AX6+q7y0z5rQkpy8sM19H8slBYyVNj5V8JTmojuQdwOnMXxLsT/Kxbuy5SXZ3b90AfCHJ48CXgc9V1f1jOQpJI5Nlzvwn6oysr0tyTJlKSSPycO3lSB3OoD5/0SipYShIahgKkhqGgqSGoSCpYShIahgKkhqGgqSGoSCpYShIahgKkhqGgqSGoSCpYShIahgKkhqGgqSGoSCpYShIagxbNu6DSea65zPuT3LNMu+9KskzSQ4kuWWUE5c0HsOWjQO4vSsHt6mqdi/tTLIO+AhwNXAxcF2Si/tMVtL4DVU2boW2AAeq6rmqehW4B9g6xHYkraI+9xRu7qpO70hy5oD+84AXFq0f7NoGsmycNB2GDYWPAm8ENgEvArf1nUhVba+qzVW1+WRO7bs5SUMaKhSq6qWq+mFV/Qj4OIPLwc0BFyxaP79rkzTFhi0bd86i1XcxuBzcI8DGJBcmOQXYBuwaZn+SVs8Jq053ZeMuA85KchC4FbgsySbmS80/D9zUjT0X+ERVXVNVR5PcDDwArAN2VNVT4zgISaNj2Tjpx5Bl4yStmKEgqWEoSGoYCpIahoKkhqEgqWEoSGoYCpIahoKkhqEgqWEoSGoYCpIahoKkhqEgqWEoSGoYCpIahoKkhqEgqbGSZzTuAK4FDlXVm7u2e4GLuiGvBf6zqjYNeO/zwHeBHwJHq2rzSGYtaWxOGArMl427A/jkQkNV/ebCcpLbgJeP8/7Lq+pbw05Q0uo6YShU1UNJ3jCoL0mA9wC/NuJ5SZqQvvcUfgV4qaqeXaa/gAeTPJrkxuNtyLJx0nRYyeXD8VwH3H2c/kurai7J2cCeJE93BWuPUVXbge0w/4j3nvOSNKShzxSSnAT8BnDvcmOqaq57PQTsZHB5OUlTpM/lw9uAp6vq4KDOJKclOX1hGbiSweXlJE2RE4ZCVzbui8BFSQ4muaHr2saSS4ck5ybZ3a1uAL6Q5HHgy8Dnqur+0U1d0jhYNk76MWTZOEkrZihIahgKkhqGgqSGoSCpYShIahgKkhqGgqSGoSCpYShIahgKkhqGgqSGoSCpYShIahgKkhqGgqTGVD5kJcl/AF9f0nwWMIv1I2b1uGB2j20Wjuv1VfUzgzqmMhQGSbJvFitMzepxwewe26we1wIvHyQ1DAVJjbUUCtsnPYExmdXjgtk9tlk9LmAN3VOQtDrW0pmCpFVgKEhqrIlQSHJVkmeSHEhyy6TnMypJnk/yRJL9SfZNej59JNmR5FCSJxe1rU+yJ8mz3euZk5zjMJY5rg8mmes+t/1JrpnkHEdt6kMhyTrgI8DVwMXAdUkunuysRuryqto0A9973wlctaTtFmBvVW0E9nbra82dHHtcALd3n9umqto9oH/NmvpQYL5S9YGqeq6qXgXuAbZOeE5aoqoeAg4vad4K3NUt3wW8czXnNArLHNdMWwuhcB7wwqL1g13bLCjgwSSPJrlx0pMZgw1V9WK3/E3miw7PipuTfKW7vFhzl0XHsxZCYZZdWlVvYf7S6H1JfnXSExqXmv/ue1a+//4o8EZgE/AicNtEZzNiayEU5oALFq2f37WteVU1170eAnYyf6k0S15Kcg5A93powvMZiap6qap+WFU/Aj7OjH1uayEUHgE2JrkwySnANmDXhOfUW5LTkpy+sAxcCTx5/HetObuA67vl64HPTnAuI7MQdJ13MWOf20mTnsCJVNXRJDcDDwDrgB1V9dSEpzUKG4CdSWD+c/h0Vd0/2SkNL8ndwGXAWUkOArcCHwY+k+QG5v8r/HsmN8PhLHNclyXZxPzl0PPATZOa3zj4M2dJjbVw+SBpFRkKkhqGgqSGoSCpYShIahgKkhqGgqTG/wLOTMif1Yv+LQAAAABJRU5ErkJggg==", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], "source": [ "# initialization\n", - "C_t[1,1] = 1" + "C_t[10,10] = 2000\n", + "plt.imshow(C_t, vmin=0, vmax=20)" ] }, { "cell_type": "code", - "execution_count": 8, + "execution_count": 5, "metadata": {}, "outputs": [], "source": [ @@ -59,7 +96,7 @@ }, { "cell_type": "code", - "execution_count": 23, + "execution_count": 51, "metadata": {}, "outputs": [], "source": [ @@ -77,143 +114,87 @@ " + max_stable_time_step/delta_y**2 * (alpha_interblock(alpha_y[i,j+1], alpha_y[i,j]) * C_t[i,j+1]\n", " - (alpha_interblock(alpha_y[i,j+1], alpha_y[i,j]) + alpha_interblock(alpha_y[i,j-1], alpha_y[i,j])) * C_t[i,j] \n", " + alpha_interblock(alpha_y[i,j-1], alpha_y[i,j]) * C_t[i,j-1])\n", + " \n", + " # boundary conditions\n", + " # left\n", + " for i in range(0, grid_size['y']):\n", + " C_t1[i,0] = C_t[i,0] \\\n", + " + max_stable_time_step/delta_x**2 * (alpha_interblock(alpha_x[i,1], alpha_x[i,0]) * C_t[i,1] \n", + " - (alpha_interblock(alpha_x[i,1], alpha_x[i,0]) + 2 * alpha_x[i,0]) * C_t[i,0]\n", + " + 2 * alpha_x[i,0] * bc_left)\n", + "\n", + " # right\n", + " n = grid_size['x'] - 1 # maximum index in x-direction (columns)\n", + " for i in range(0, grid_size['y']):\n", + " C_t1[i,n] = C_t[i,n] \\\n", + " + max_stable_time_step/delta_x**2 * (2 * alpha_x[i,n] * bc_right \n", + " - (alpha_interblock(alpha_x[i,n-1], alpha_x[i,n]) + 2 * alpha_x[i,n]) * C_t[i,n]\n", + " + alpha_interblock(alpha_x[i,n-1], alpha_x[i,n]) * C_t[i,n-1])\n", + "\n", + " # top\n", + " for j in range(0, grid_size['x']):\n", + " C_t1[0,j] = C_t[0,j] \\\n", + " + max_stable_time_step/delta_y**2 * (alpha_interblock(alpha_y[1,j], alpha_x[0,j]) * C_t[1,j] \n", + " - (alpha_interblock(alpha_y[1,j], alpha_y[0,j]) + 2 * alpha_y[0,j]) * C_t[0,j]\n", + " + 2 * alpha_y[0,j] * bc_top)\n", + "\n", + " # bottom\n", + " m = grid_size['y'] - 1 # maximum index in y-direction (rows)\n", + " for j in range(0, grid_size['x']):\n", + " C_t1[m,j] = C_t[m,j] \\\n", + " + max_stable_time_step/delta_y**2 * (2 * alpha_y[m,j] * bc_bottom \n", + " - (alpha_interblock(alpha_y[m,j], alpha_y[m-1,j]) + 2 * alpha_y[m,j]) * C_t[m,j]\n", + " + alpha_interblock(alpha_y[m,j], alpha_y[m-1,j]) * C_t[m-1,j])\n", + "\n", "\n", " return C_t1" ] }, { "cell_type": "code", - "execution_count": 27, + "execution_count": 115, "metadata": {}, "outputs": [ { - "name": "stdout", - "output_type": "stream", - "text": [ - "[[ 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", - " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", - " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", - " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", - " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00]\n", - " [ 0.00000000e+00 4.93210938e+02 -6.85390625e+02 5.61357422e+02\n", - " -3.21421875e+02 1.35512695e+02 -4.28554688e+01 1.01855469e+01\n", - " -1.79687500e+00 2.28515625e-01 -1.95312500e-02 9.76562500e-04\n", - " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", - " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00]\n", - " [ 0.00000000e+00 -6.85390625e+02 9.41660156e+02 -7.56210938e+02\n", - " 4.20468750e+02 -1.70039062e+02 5.07128906e+01 -1.10742188e+01\n", - " 1.71875000e+00 -1.75781250e-01 9.76562500e-03 0.00000000e+00\n", - " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", - " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00]\n", - " [ 0.00000000e+00 5.61357422e+02 -7.56210938e+02 5.86933594e+02\n", - " -3.10078125e+02 1.16542969e+02 -3.12890625e+01 5.84472656e+00\n", - " -7.03125000e-01 4.39453125e-02 0.00000000e+00 0.00000000e+00\n", - " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", - " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00]\n", - " [ 0.00000000e+00 -3.21421875e+02 4.20468750e+02 -3.10078125e+02\n", - " 1.51593750e+02 -5.08593750e+01 1.15312500e+01 -1.64062500e+00\n", - " 1.17187500e-01 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", - " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", - " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00]\n", - " [ 0.00000000e+00 1.35512695e+02 -1.70039062e+02 1.16542969e+02\n", - " -5.08593750e+01 1.43554688e+01 -2.46093750e+00 2.05078125e-01\n", - " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", - " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", - " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00]\n", - " [ 0.00000000e+00 -4.28554688e+01 5.07128906e+01 -3.12890625e+01\n", - " 1.15312500e+01 -2.46093750e+00 2.46093750e-01 0.00000000e+00\n", - " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", - " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", - " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00]\n", - " [ 0.00000000e+00 1.01855469e+01 -1.10742188e+01 5.84472656e+00\n", - " -1.64062500e+00 2.05078125e-01 0.00000000e+00 0.00000000e+00\n", - " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", - " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", - " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00]\n", - " [ 0.00000000e+00 -1.79687500e+00 1.71875000e+00 -7.03125000e-01\n", - " 1.17187500e-01 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", - " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", - " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", - " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00]\n", - " [ 0.00000000e+00 2.28515625e-01 -1.75781250e-01 4.39453125e-02\n", - " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", - " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", - " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", - " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00]\n", - " [ 0.00000000e+00 -1.95312500e-02 9.76562500e-03 0.00000000e+00\n", - " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", - " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", - " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", - " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00]\n", - " [ 0.00000000e+00 9.76562500e-04 0.00000000e+00 0.00000000e+00\n", - " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", - " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", - " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", - " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00]\n", - " [ 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", - " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", - " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", - " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", - " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00]\n", - " [ 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", - " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", - " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", - " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", - " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00]\n", - " [ 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", - " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", - " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", - " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", - " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00]\n", - " [ 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", - " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", - " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", - " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", - " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00]\n", - " [ 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", - " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", - " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", - " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", - " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00]\n", - " [ 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", - " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", - " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", - " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", - " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00]\n", - " [ 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", - " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", - " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", - " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", - " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00]\n", - " [ 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", - " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", - " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", - " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00\n", - " 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00]]\n" - ] + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAT4AAAD8CAYAAADub8g7AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy88F64QAAAACXBIWXMAAAsTAAALEwEAmpwYAAAYa0lEQVR4nO3dfawd1X3u8e8TY0BxoOC4cXgriVoLiaLitJZpFNqaEIixuCGpuKndKoWEyiEKUpBa9dJWgoj8k6ol6PaSG18HLENFHO5N4gQ1DuDSSA5SQjCWSWxeioOIsHFsGRNe8gL4nOf+MXPSne29zxnPnnPO3meeTzQ6M2vWzFobKz+tl5lZsk1ERJu8abYrEBEx0xL4IqJ1EvgionUS+CKidRL4IqJ1EvgionUS+CJiVkk6S9K3JT0uabekT5XpCyVtlfR0+ffUPtdfVeZ5WtJVlcrMc3wRMZsknQacZnuHpJOAR4EPAlcDh21/VtINwKm2/0fXtQuB7cAywOW1f2D7xcnKTIsvImaV7f22d5T7rwBPAGcAVwB3ltnupAiG3d4PbLV9uAx2W4GVU5V5XAP1btzxOsEnsmC2qxExZ/2Sn/G6X9Mg93j/RQv8wuGxSnkf/cFru4FfdiStt72+O5+kdwDvAh4GFtveX576CbC4x63PAJ7rON5bpk1qKAPfiSzgAl0829WImLMe9oMD3+PQ4TEevv/MSnnnn/ajX9peNlkeSW8Bvgpcb/tl6b/ism1LamxcLl3diKjJjHm80jYVSfMpgt7dtr9WJh8ox/8mxgEP9rh0H3BWx/GZZdqkBgp8klZKekrSnnLwsfv8CZLuKc8/XDZjI2IOMDCOK22TUdG0uwN4wvbnOk7dC0zM0l4FfKPH5fcDl0o6tZz1vbRMm1TtwCdpHvB54DLgXGCNpHO7sl0DvGj7d4BbgX+sW15EDJ/xiv+bwnuAjwDvlbSz3FYBnwUukfQ08L7yGEnLJN0OYPsw8BngkXK7uUyb1CBjfMuBPbafKSvzZYpZmMc78lwBfLrc/wpwmyQ5z9BEjDxj3qjQjZ3yPvZDQL+JlqMG+21vB/6q43gDsOFYyhykq1tlNuVXeWwfAV4C3trrZpLWStouafsbvDZAtSJiJhgYw5W2YTM0s7rl1PZ6gJO1cPj+S0XEUaYavxtWgwS+KrMpE3n2SjoO+A3ghQHKjIghYWBsREetBunqPgIskfROSccDqylmYTp1zspcCfxHxvci5o7xituwqd3is31E0nUUU8fzgA22d0u6Gdhu+16KKep/lbQHOEwRHCNiDvCQjt9VMdAYn+0twJautBs79n8J/PdByoiI4WTDG6MZ94ZnciMiRo0Y6/sUynBL4IuIWgyMp8UXEW2TFl9EtErxAHMCX0S0iIE3PJofeErgi4hajBgb0S/bJfBFRG3jTlc3IlokY3wR0UJiLGN8EdEmxReYE/giokVs8brnzXY1akngi4jaxjPGFxFtUkxuNNPVlbQBuBw4aPu8Mu0e4JwyyynAT20v7XHts8ArwBhwZKplLCGBLyJqa3RyYyNwG3DXRILtP/tVSdItFEtX9HOR7UNVC0vgi4hampzcsL2t3/Kz5fKTHwbe20hhZEHxiBjAmFVpG9AfAQdsP93nvIEHJD0qaW2VG6bFFxG1GPGGK4eQRZK2dxyvLxcYq2INsGmS8xfa3ifpbcBWSU/a3jbZDRP4IqKWY5zcOFRl0qFbuUjZnwJ/0Lce9r7y70FJmynW/J408NXu6ko6S9K3JT0uabekT/XIs0LSSx2ro9/Y614RMXpMtW7ugF3d9wFP2t7b66SkBZJOmtgHLgV2TXXTQVp8R4C/tr2jLPhRSVttP96V7zu2Lx+gnIgYUk1NbkjaBKyg6BLvBW6yfQfFAmWbuvKeDtxuexWwGNhczH9wHPAl2/dNVd4gq6ztB/aX+69IegI4A+gOfBExB9k09jiL7TV90q/ukfY8sKrcfwY4/1jLa2SMr5yGfhfwcI/T75b0GPA88De2d/e5x1pgLcCJvLmJakXENComN1r6ypqktwBfBa63/XLX6R3A2bZflbQK+DqwpNd9yhme9QAna+GILmES0S6j+iHSgWotaT5F0Lvb9te6z9t+2far5f4WYL6kRYOUGRHDwYhxV9uGTe0WX/k09R3AE7Y/1yfP2ykePLSk5RSB9oW6ZUbEcBnVFt8gXd33AB8BfihpZ5n298BvAdheB1wJfELSEeAXwGrb6cZGzAHFurotC3y2H4LJv0lj+zaKF48jYs5RPj0fEe1SLC/Z0lndiGgnW+3r6kZEZLGhiGiV4nt8GeOLiFbJ8pIR0TLF4yxp8UVEi7T6Xd2IaK8sKB4RrVJ8lipd3YhomYzxRUSrFF9nSVc3IlqkeGUtgS8iWmV0W3yjWeuIGArjqNI2FUkbJB2UtKsj7dOS9nWs0riqz7UrJT0laY+kG6rUO4EvImqZmNVtaHnJjcDKHum32l5ablu6T0qaB3weuAw4F1gj6dypCkvgi4jaxv2mSttUbG8DDteownJgj+1nbL8OfBm4YqqLEvgiopZjXHNjkaTtHdvaisVcJ+kHZVf41B7nzwCe6zjeW6ZNKpMbEVGLgSPVJzcO2V52jEV8AfhMWdRngFuAjx3jPXpqYnnJZ4FXgDHgSPePKxcl+p8UCwD/HLja9o5By42I2Teds7q2D0zsS/oi8G89su0Dzuo4PrNMm1RTLb6LbB/qc+4yirV0lwAXUETxCxoqNyJmyzQvHSnpNNv7y8MPAbt6ZHsEWCLpnRQBbzXw51Pdeya6ulcAd5Wrq31P0ildPygiRlCTHyKVtAlYQTEWuBe4CVghaWlZ1LPAx8u8pwO3215l+4ik64D7gXnABtu7pyqvicBn4AFJBv6P7fVd5/sNPv5a4CsHO9cCnMibG6hWREy3plp8ttf0SL6jT97nKYbOJo63AEc96jKZJgLfhbb3SXobsFXSk+XU9DEpA+Z6gJO1MGvvRgy5Vn+I1Pa+8u9BSZspnqvpDHy1Bh8jYrgZcWR8NJ+IG6jWkhZIOmliH7iUowcg7wX+UoU/BF7K+F7E3NDUK2szbdAW32Jgc/HECscBX7J9n6RrAWyvo+h7rwL2UDzO8tEBy4yIYeCWdnVtPwOc3yN9Xce+gU8OUk5EDJ9Wj/FFRHsl8EVEqxgxNqKTGwl8EVHbME5cVJHAFxG1uK2TGxHRbk7gi4h2md6PFEynBL6IqC0tvohoFRvGxhP4IqJlMqsbEa1i0tWNiNbJ5EZEtJBH9MuZCXwRUVu6uhHRKsWsbjPv6kraAFwOHLR9Xpn2T8B/A14HfgR81PZPe1z7LJOs9NjLaL5hHBFDwa62VbARWNmVthU4z/bvAf8J/N0k119ke2nVtXsT+CKiNluVtqnv423A4a60B2wfKQ+/R7FsRSMS+CKiFlMt6JWBb5Gk7R3b2mMs7mPAt/pWpVjp8dGq980YX0TUdgyTuoeqdkO7SfoH4Ahwd58sx7zSY+0Wn6RzJO3s2F6WdH1XnhWSXurIc2Pd8iJiyBg8rkpbXZKuppj0+ItyGYujq9Gx0iMwsdLjpGq3+Gw/BSwtKzePYsnIzT2yfsf25XXLiYjhNZ2Ps0haCfwt8Ce2f94nzwLgTbZf6Vjp8eap7t3UGN/FwI9s/7ih+0XECGhqVlfSJuC7wDmS9kq6BrgNOImi+7pT0roy7+mStpSXLgYekvQY8H3gm7bvm6q8psb4VgOb+px7d1mp54G/sb27V6ZyUHItwIm8uaFqRcR0afJdXdtreiTf0Sfv8xRL1vZd6XEqA7f4JB0PfAD4fz1O7wDOtn0+8L+Ar/e7j+31tpfZXjafEwatVkRMNwNWtW3INNHVvQzYYftA9wnbL9t+tdzfAsyXtKiBMiNiCDT4APOMaqKru4Y+3VxJbwcO2Lak5RSB9oUGyoyIWTfYjO1sGijwlbMolwAf70i7FsD2OuBK4BOSjgC/AFb3m5KOiBE0ov9vHijw2f4Z8NautHUd+7dRzMxExFzjfJ0lItqojS2+iGi7tPgiom3GZ7sC9STwRUQ9E8/xjaAEvoiobVSf0Ujgi4j6EvgionXS1Y2ItlFafBHRKha08ZW1iGi5tPgionUS+CKidRL4IqJVRvgB5qyrGxG1ydW2Ke8jbZB0UNKujrSFkrZKerr8e2qfa68q8zwt6aoq9U7gi4j6XHGb2kZgZVfaDcCDtpcAD5bHv0bSQuAm4AKKZSVv6hcgOyXwRURtTbX4ygXAD3clXwHcWe7fCXywx6XvB7baPmz7RWArRwfQo2SMLyLqqz7Gt0jS9o7j9bbXT3HNYtv7y/2fUCwl2e0M4LmO471l2qQS+CKinurdWIBDtpfVLqpYt6exOeRKXd2ZHniMiBHR3BhfLwcknQZQ/j3YI88+4KyO4zPLtElVHePbyAwOPEbEaNB4ta2me4GJxtJVwDd65LkfuFTSqWVsubRMm1SlwDfTA48RMSIaavFJ2gR8FzhH0l5J1wCfBS6R9DTwvvIYScsk3Q5g+zDwGeCRcru5TJvUIGN80zbwGBHDr+qMbRW21/Q5dXGPvNuBv+o43gBsOJbyGpncaGLgUdJaYC3Aiby5iWpFxHRr4ZsbjQ482l5ve5ntZfM5YYBqRcSMmd7JjWkzSOCbtoHHiBgNTT3APNOqPs4yowOPETECPO2zutOm0hjfTA88RsSIGMLWXBV5cyMi6kvgi4i2GcbxuyrydZaIaJ20+CKivhFt8SXwRUQ9Hs4Z2yoS+CKivrT4IqJNxOhObiTwRUR9CXwR0SpD+jpaFQl8EVFfJjciom3S4ouI9kngi4hWGdJv7VWRV9YiorYmvscn6RxJOzu2lyVd35VnhaSXOvLcOEi90+KLiPoaaPHZfgpYCiBpHsVX2jf3yPod25cPXmICX0QMYBpeWbsY+JHtHzd+5w7p6kZEPVXX2yhahYskbe/Y1va562pgU59z75b0mKRvSfrdQaqeFl9E1KJyq+iQ7WWT3k86HvgA8Hc9Tu8Azrb9qqRVwNeBJdWL/3Vp8UVEfc2usnYZsMP2gaOKsV+2/Wq5vwWYL2lR3WpPGfgkbZB0UNKujrR/kvSkpB9I2izplD7XPivph+UszPa6lYyI4dTwKmtr6NPNlfR2SSr3l1PErhfq1rtKi28jsLIrbStwnu3fA/6T3k3TCRfZXjpVMzciRlBDLT5JC4BLgK91pF0r6dry8Epgl6THgH8BVtuuPac85Rif7W2S3tGV9kDH4ffKSkVEmzT4IVLbPwPe2pW2rmP/NuC2ZkprZozvY8C3+pwz8ICkRyeZxQFA0tqJGZ83eK2BakXEtGt2jG/GDDSrK+kfgCPA3X2yXGh7n6S3AVslPWl7W6+MttcD6wFO1sIh/E8VEd1G9SMFtVt8kq4GLgf+ol9f2/a+8u9Biiexl9ctLyKG0Ii2+GoFPkkrgb8FPmD7533yLJB00sQ+cCmwq1feiBhNDc/qzpgqj7NsAr4LnCNpr6RrKAYZT6Lovu6UtK7Me7qkLeWli4GHylmY7wPftH3ftPyKiJh5pvgQaZVtyFSZ1V3TI/mOPnmfB1aV+88A5w9Uu4gYWllsKCLaKYEvItpG9Z8hnlUJfBFRz5DO2FaRwBcRtWWMLyJaZxo+RDojEvgior60+CKiVYb04eQqEvgior4EvohokzzAHBGtpPHRjHwJfBFRT57ji4g2GtXHWbLKWkTU19yaG5MuTKbCv0jaUy5y9vuDVDstvoioreHJjYtsH+pz7jKKdXSXABcAXyj/1pIWX0TUY8Cutg3uCuAuF74HnCLptLo3S+CLiNo0Xm0DFk0sJlZu3YuPTbUw2RnAcx3He8u0WtLVjYhajvE5vkNTrK1deWGyJqTFFxH1VO3mVujqVliYbB9wVsfxmWVaLVXW3Ngg6aCkXR1pn5a0r5yB2SlpVZ9rV0p6qpyJuaFuJSNiODWx2FDFhcnuBf6ynN39Q+Al2/vr1rtKV3cjxeJCd3Wl32r7n/tdJGke8HngEor++COS7rX9eM26RsSwaWZWdzGwWRIUMelLtu+TdC2A7XXAFor1fPYAPwc+OkiBVRYb2ibpHTXuvRzYUy46hKQvU8zMJPBFzBFNPM7Sb2GyMuBN7Bv45OClFQYZ47uufJBwg6RTe5w/plkYSWsnZnze4LUBqhURM8LAmKttQ6Zu4PsC8NvAUmA/cMugFbG93vYy28vmc8Kgt4uIGTCqC4rXepzF9oGJfUlfBP6tR7ZGZ2EiYgiN6CprtVp8XU9Mf4ijZ2AAHgGWSHqnpOOB1RQzMxExR8zZFp+kTcAKiiev9wI3ASskLaXo5T8LfLzMezpwu+1Vto9Iug64H5gHbLC9ezp+RETMgrn8WSrba3ok39En7/MUU84Tx1sopqEjYo4RoCGcuKgir6xFRG0a0TG+BL6IqGcud3UjInpr7JNTMy6BLyJqG8YZ2yoS+CKivrT4IqJVnFndiGij0Yx7CXwRUV8eZ4mI9kngi4hWMTCiC4on8EVELcLp6kZEC42PZpMvq6xFRD0TXd0q2yQknSXp25Iel7Rb0qd65Fkh6aWOBc5uHKTqafFFRG0NdXWPAH9te0e52tqjkrb2WJjsO7Yvb6LABL6IqK+BwFcuE7m/3H9F0hMU6/NM28Jk6epGRE3NLSg+oVzR8V3Awz1Ov1vSY5K+Jel3B6l5WnwRUc/EKmvVLJK0veN4ve31nRkkvQX4KnC97Ze7rt8BnG37VUmrgK8DS2rVmwS+iBjAMYzxHbK9rO99pPkUQe9u21/rPt8ZCG1vkfS/JS2yfehY6wzV1tzYAFwOHLR9Xpl2D3BOmeUU4Ke2l/a49lngFWAMODLZD4+IEdTAGJ8kUSxn8YTtz/XJ83bggG1LWk4xTPdC3TKrtPg2ArcBd00k2P6zjgrdArw0yfUX1Y3KETHEDIw3Mqv7HuAjwA8l7SzT/h74LQDb64ArgU9IOgL8Alht14+6VRYb2lYOOB6ljNQfBt5btwIRMaqa+QKz7Yco1i6aLM9tFA2wRgw6q/tHFM3Pp/ucN/CApEclrZ3sRpLWStouafsbvDZgtSJiRjQ8qztTBp3cWANsmuT8hbb3SXobsFXSk7a39cpYzvCsBzhZC4fvv1RE/DoDYy17ZU3SccCfAvf0y2N7X/n3ILAZWF63vIgYNgaPV9uGzCBd3fcBT9re2+ukpAXl6ydIWgBcCuwaoLyIGDYj2tWdMvBJ2gR8FzhH0l5J15SnVtPVzZV0uqQt5eFi4CFJjwHfB75p+77mqh4Rs2piVrfKNmSqzOqu6ZN+dY+054FV5f4zwPkD1i8ihtkQtuaqyJsbEVFfAl9EtIoNY2OzXYtaEvgior60+CKidRL4IqJdhnPGtooEvoiox+AhfDi5igS+iKhvRF9ZS+CLiHrskV1eMoEvIurL5EZEtI3T4ouIdhnODxBUkcAXEfU09+n5GZfAFxG1GPCIvrKWBcUjoh439yFSSSslPSVpj6Qbepw/QdI95fmH+60DVFUCX0TU5nFX2iYjaR7weeAy4FxgjaRzu7JdA7xo+3eAW4F/HKTeCXwRUV8zLb7lwB7bz9h+HfgycEVXniuAO8v9rwAXl6s81jKUY3yv8OKhf/dXftyVvAiYi+vzztXfBXP3t82F33X2oDd4hRfv/3d/ZVHF7CdK2t5xvL5cYAzgDOC5jnN7gQu6rv9VHttHJL0EvJWa/w5DGfhs/2Z3mqTttpfNRn2m01z9XTB3f9tc/V3HyvbK2a5DXenqRsRs2wec1XF8ZpnWM0+5wuNvAC/ULTCBLyJm2yPAEknvlHQ8xUJm93bluRe4qty/EvgPu/7T00PZ1e1j/dRZRtJc/V0wd3/bXP1ds6Ics7sOuB+YB2ywvVvSzcB22/cCdwD/KmkPcJgiONamAYJmRMRISlc3IlongS8iWmckAt9Ur7OMKknPSvqhpJ1dzziNHEkbJB2UtKsjbaGkrZKeLv+eOpt1rKPP7/q0pH3lv9tOSatms45x7IY+8FV8nWWUXWR76Rx4Lmwj0P1c1w3Ag7aXAA+Wx6NmI0f/LoBby3+3pba3zHCdYkBDH/io9jpLzDLb2yhm2zp1vmZ0J/DBmaxTE/r8rhhxoxD4er3OcsYs1aVpBh6Q9KiktbNdmWmw2Pb+cv8nwOLZrEzDrpP0g7IrPHJd+LYbhcA3l11o+/cpuvGflPTHs12h6VI+bDpXnp36AvDbwFJgP3DLrNYmjtkoBL4qr7OMJNv7yr8Hgc0U3fq55ICk0wDKvwdnuT6NsH3A9piLRWW/yNz7d5vzRiHwVXmdZeRIWiDppIl94FJg1+RXjZzO14yuAr4xi3VpzEQwL32IuffvNucN/Str/V5nmeVqNWExsLn8pNhxwJds3ze7VapP0iZgBbBI0l7gJuCzwP+VdA3wY+DDs1fDevr8rhWSllJ03Z8FPj5b9Yt68spaRLTOKHR1IyIalcAXEa2TwBcRrZPAFxGtk8AXEa2TwBcRrZPAFxGt8/8BNnHxTGa4MTYAAAAASUVORK5CYII=", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" } ], "source": [ - "for i in range(10):\n", + "records = []\n", + "\n", + "for i in range(1000):\n", " C_t = simulate(C_t)\n", - "print(C_t)" + " records.append(C_t)\n", + "\n", + "# with np.printoptions(precision=2,floatmode='fixed',linewidth=200):\n", + "# print(C_t)\n", + "\n", + "plt.imshow(C_t, vmin=0, vmax=20)\n", + "plt.colorbar()\n", + "plt.show()\n" ] }, { "cell_type": "code", - "execution_count": 28, + "execution_count": 116, "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "0.5\n" - ] - } - ], + "outputs": [], "source": [ - "print(max_stable_time_step)" + "def update(w = 1):\n", + " fig = plt.imshow(figsize = (15,10))\n", + " y = records[w]\n", + " plt.imshow(y)\n", + " \n", + "interact(update, w = IntSlider(min=0, max = 99, step = 1, value = 0))" ] }, {