231 lines
20 KiB
Plaintext
231 lines
20 KiB
Plaintext
{
|
|
"cells": [
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"# imports\n",
|
|
"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": 113,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"# environment\n",
|
|
"dim = 2\n",
|
|
"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",
|
|
"alpha_x = np.ones((grid_size['x'],grid_size['y']))\n",
|
|
"alpha_y = np.ones((grid_size['x'],grid_size['y']))\n",
|
|
"\n",
|
|
"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": 114,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"data": {
|
|
"text/plain": [
|
|
"<matplotlib.image.AxesImage at 0x122e2fa60>"
|
|
]
|
|
},
|
|
"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": [
|
|
"<Figure size 432x288 with 1 Axes>"
|
|
]
|
|
},
|
|
"metadata": {
|
|
"needs_background": "light"
|
|
},
|
|
"output_type": "display_data"
|
|
}
|
|
],
|
|
"source": [
|
|
"# initialization\n",
|
|
"C_t[10,10] = 2000\n",
|
|
"plt.imshow(C_t, vmin=0, vmax=20)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 5,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"def alpha_interblock(alpha1, alpha2, harmonic=False):\n",
|
|
" if not harmonic:\n",
|
|
" return 0.5 * (alpha1 + alpha2)\n",
|
|
" else:\n",
|
|
" return 2 / ((1/alpha1) + (1/alpha2))\n"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 51,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"# simulating one time step\n",
|
|
"def simulate(C_t): \n",
|
|
" C_t1 = np.zeros((grid_size['x'],grid_size['y']))\n",
|
|
"\n",
|
|
" # inner cells\n",
|
|
" for i in range(1, grid_size['x']-1):\n",
|
|
" for j in range(1, grid_size['y']-1):\n",
|
|
" C_t1[i,j] = C_t[i,j] \\\n",
|
|
" + max_stable_time_step/delta_x**2 * (alpha_interblock(alpha_x[i+1,j], alpha_x[i,j]) * C_t[i+1,j]\n",
|
|
" - (alpha_interblock(alpha_x[i+1,j], alpha_x[i,j]) + alpha_interblock(alpha_x[i-1,j], alpha_x[i,j])) * C_t[i,j] \n",
|
|
" + alpha_interblock(alpha_x[i-1,j], alpha_x[i,j]) * C_t[i-1,j]) \\\n",
|
|
" + 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": 115,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"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": [
|
|
"<Figure size 432x288 with 2 Axes>"
|
|
]
|
|
},
|
|
"metadata": {
|
|
"needs_background": "light"
|
|
},
|
|
"output_type": "display_data"
|
|
}
|
|
],
|
|
"source": [
|
|
"records = []\n",
|
|
"\n",
|
|
"for i in range(1000):\n",
|
|
" C_t = simulate(C_t)\n",
|
|
" 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": 116,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"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))"
|
|
]
|
|
},
|
|
{
|
|
"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.10.1"
|
|
},
|
|
"orig_nbformat": 4
|
|
},
|
|
"nbformat": 4,
|
|
"nbformat_minor": 2
|
|
}
|