946 lines
91 KiB
Plaintext
946 lines
91 KiB
Plaintext
{
|
||
"cells": [
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 2,
|
||
"metadata": {},
|
||
"outputs": [
|
||
{
|
||
"name": "stdout",
|
||
"output_type": "stream",
|
||
"text": [
|
||
"Populating the interactive namespace from numpy and matplotlib\n"
|
||
]
|
||
}
|
||
],
|
||
"source": [
|
||
"%pylab inline"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 3,
|
||
"metadata": {},
|
||
"outputs": [],
|
||
"source": [
|
||
"import pandas\n",
|
||
"import sympy as Ω\n",
|
||
"import matplotlib.pyplot as plt\n",
|
||
"import numpy as np\n",
|
||
"from scipy.misc import derivative\n",
|
||
"import pickle\n",
|
||
"from tqdm import tqdm_notebook as tqdm\n",
|
||
"\n",
|
||
"from dataclasses import dataclass\n",
|
||
"\n",
|
||
"import scipy.sparse\n",
|
||
"import scipy.sparse.linalg\n",
|
||
"\n",
|
||
"import functools"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 4,
|
||
"metadata": {},
|
||
"outputs": [],
|
||
"source": [
|
||
"Ω.init_printing()"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 5,
|
||
"metadata": {},
|
||
"outputs": [],
|
||
"source": [
|
||
"i, k = Ω.symbols('i k')\n",
|
||
"I, K = Ω.symbols('I K')\n",
|
||
"h_x, h_t = Ω.symbols('h_x h_t')\n",
|
||
"l, T, H, D = Ω.symbols('l T H D')"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 6,
|
||
"metadata": {},
|
||
"outputs": [],
|
||
"source": [
|
||
"u = Ω.IndexedBase('u', shape=(I + 1, K + 1))"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 7,
|
||
"metadata": {},
|
||
"outputs": [],
|
||
"source": [
|
||
"gamma = Ω.symbols('gamma')\n",
|
||
"gamma_val = D * h_t / h_x**2"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"metadata": {
|
||
"heading_collapsed": true
|
||
},
|
||
"source": [
|
||
"## простейшая явная"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 91,
|
||
"metadata": {
|
||
"hidden": true
|
||
},
|
||
"outputs": [
|
||
{
|
||
"data": {
|
||
"image/png": "iVBORw0KGgoAAAANSUhEUgAABCwAAAAZCAYAAAD66r+dAAAABHNCSVQICAgIfAhkiAAADTdJREFUeJztnXmsH1UVxz+FUrshLRAhSCgUpaUh9tFaLYLNA3GLohYhjWhwAgSxVoMEFwhiUQJGjVIMLkDIg4gBpanaGKBEfezWsslWLNuDyCqISxGpyPOPcyedN7/5zXrvnZnfnE/yy++9O3dmzpxzft+ZuffOHVAURVEURVEURVEURWkww8B45PNQrdYoiqIoiqIoaVwOPA/MqNsQxRmLkevykyps44tmG8dZsUhpOqoLSh6qaottXdmdiW0R40mVhs2CUWA1sMrSzhVFGRx2Q4RtHfAI8ArwD+AW4ERgh/pMUxSlxWxArkGGU+r82NQ52YdBLWAJ8DpwWt2GKM5ZBzwDzMyotwfwP+AHsfKfIr+d+fZN845qRTqqC0oRqmiLbV2ZjrRBrAbGyGiwWG1pp4qiDB6nIDrxNHAlcD5wGfB3U34NMKk26xRFaSsvAq+R3iN4N6Izi7xY1Hw2INo7rW5DFOe8A8n9MzPqnWzqHR4rfxD4F4PRqaBakY7qglKEKtriUldG0QYLRek0Adm9E/04AjiKXnHaE3jSbPfjFWxTFJcElM99xR1zkbj8KaXOdOC/wH+AnXwY5YGA8vl4ANKLerFFewaFgMH8nW8GniD95uBa4AVgx0jZDKRn9BZ3pnmjC1oRoLrQJgLarzdltMW1rowSabDw2dK6HBGQjcCcPnXOYeIQrjPM/0cn1N3HLFtn10ylAhrjwed3wHrkhBjlWWQIJvSKtuZFvagv7aG+dMPbzffGlDqLgcnAvYiehHQ1Jicgo9mudriPrvq2qcd9ldn/e/ssfyPSqbAeuZEIGUKu9+8CFiDDuJ8FtgK3Ae90ZK8LVCvSUV1oF03xZRlt8aorPhssngduQIaefDlh+VxTvgm41JQtNt93JNRfYr7vsmijUg2NcbcJLwxei5VrXtSL+tIe6ks35LkJWWq+477vakyORC4c/+BwH131bVOP+1bz3e+m4sPAFHpvcMLHIuYgx7QLcAVwI3AI8GtgZ6uWukO1Ih3VhXbRFF+W0RavujLZ1oZycCvS0/oi2wUnyoWIM1ayvfd2MTL85MmE+uE27rZrplIBjXF3mQwcb/6+LrZM86Je1Jf2UF+6IfTdCnqfvQ8Je2uSbkK6FpMZSO/WZuBlh/vpom+huce9yXwv67N8OZIPG2Ll4U3RUuAwJt4ArUV6d4eAm+2Y6RTViv6oLrSPpviyjLbUpivD+JnD4o/IWwWifMTs+0eRsl1NWfzmJ+QGs3yvWPm3zLImcxGDPVxKY9zsGAfYf97uu2abv0mp4zIv2pAT4D8v1JcTCSif+3X4ssk6YotJbJ+0N8/nbZF1234OCSiXjweY9eI3pjZps28D/P7Owd9xv4IMu44zFZn87hcJy+5FbP5gwrJzzbJ+PatNoitaEaC60CYCysWrab4sqi2udWWUkpNujpFfJMaRZ1mSuMIs39v8Pw14DPgrMDtS70hT79w+2/kb8FxC+XXAd/qsY4NlyDCXpxD7ghLbmE3262N82DSGnZjG0Ri3N8YjJWz7gll3MyLA/XCZF65zAvzkxRh2f5PqS3u5X4cvbeuILcawl6fhRfbtKXVmmzr/ZuJkgnWcQ1YCjyMT+t0JvDvnemPYy8dDTJ2059SL7i8eo7b4doz6f+dQ/riLauFT9D52CTIZ9jhwXKx8KvK45qN9tvczs97cfObWStu04reIf9MYQ3Uhizx+LEKV648x7MWrab4soi0+dGXUbAMo9kjIo8hJJC9P9yl/yHwvAP6CTDiyH3Ai8FKkXjjU5M6EbcxFRCmpVWoIuTFyxUzgfrOPsvt5KbsKIIk+Rvaol7I22YppHI1x/hjnpaxNFwCzYmVDwEeBy5H8inJPQbtWAWuQVxu9BxHYfrjMC9c5AX7ywvZvssu+tJ37dfjSto7YwmaehsNe057TDZ+VvYeJkwn6PoesQPRuJTIz+kpk5vQFJA/pjWIzH18x31NT6lSNUVt824TfebjPMsddVAunsT3+UY4GttE7ynEhcq3fr9d9ETLy8fE8xtZMm7QitKXfzWiI6kI2efxYhCrXHzbj1TRfFtGWWnVlGGnJWO1i4xGONvs5Fdgf+eHchgz1inKVqTcnYRurzLLzYuV7mvIF5v8ZZjt3AftWN72HrRTvGdwbsXF+jrojFI9HGZtsozHOH+MyVI1xgNg3XNGOU8127gPelKO+q7zwnRPQzLxIQn05kYDyue/bl3XkSx18DznOE1PqfMnUuTBW7vscshG4JFb2MHB+xnr9CCiXj3uZ9Vy+prLNvg3w9zsHe7qZpYU7IPM8xXs1d0Seg782YZ3PGttOSli2s9ne72PlxwCvMtEHa8x+90ixzzVt0or9KZ+DQcl1B1EXqvgxDzbumQLK2dgkXxbVFh+6MkpkhIXPt4SERHtZ1yAtNJ+LGmWYjww3eSJW/gbgM+bveCvrENI69GdgHvLM/GvAofS2etXFQmSo2pa6DXGIxnjwY/wV4PtI6/HhyJtAsnCVF23ICagnL9SX9vDtyy7oCOTrNT3YfMcn0fN5DpmC9IjFe5Q2AO9KWc8FzyCP0s1zuI+u+rbJ1yXzkAb+eK/tMmA3kue7SevFPdhsL75sLdIRcZb5/3TgE8AHSB6i7ou2aAWI38czbLXNoOkC1ONHXzTJl0W1xbuu1NFg8QgyTOtY4EPIZHtJM6BuA3YC3hopmwFcBhxk/o+vN4Q442NIz+0lwKdIHuJSF6GNr2dVbDEa48GO8deQSX/uRB4DeSHneq7yog05AfXkhfrSHr59mXSMIxR/5rbJ7IBc3GwDHkipFw7zjt+E+DyH7I70NsUvrJ5Der18Mg7cZGx6i6N9dNW3Tb4uCV/XGe+5XI7oxK8S1lmE9Gren7AsvOmI38CMA2ciOvNV4GzknP1wpM4IfrWoTVoB4ttHgX/mqGuLQdMFqMePvmiSL4tqiytd6UsdDRbbkGdaZiG9smf1qXe9+b4J+CHyvM4WZFK/Z5BnYx6LrTOEBP4yZPj5BX22Hc5emvYZLnRU+VlI/2eazkSGJ4WfTyaU5Z3gq040xv1jXKddNvg08A2k4eFmZMLN1bFP0GddV3mRNyeguXnhCvWlPVz6MomkYwzP2UkTY7WR+cjzxA8g+pDETMSPW9k+SivE5Tmk6aw13+93tP2u+rbJx/0+5NwbvXmYhNy83E5vg88U5MbnPqQnN05aL+kG5FWH5yLzi2yKLfetRW3TikUk+9U1g6YLaX5s+/V0k3xZRFtc6kouhvEzhwXAerJbZqciw8afRobF3oEMkZmFtPaMJqzzEHCpqb88Zdu7I+KX9pme4zjKPPu0BTilz7JdkVbR8LMWeQ4vWjbNgU0u0Bi7s6vOOSxWk32CGE1Z30Ve5M0JaG5euEJ9OZGA8rnv0pdJJB3j3UjvyOze6q3keCQel6bUOZTtPYdxXJ5D4kxBbs6OjZVfBNxYYDtRAsrn4xTkInJjyX1n0WbfBvj7nUP133pImhbugvSi/jJWvgQ51tMS1gmHgP+kzzY3I68rjM8jBXAE8DJyE7MoYblvLWqTVgC8iMynUYYA1YWQND824XoayserKb4sqi0udSXKKL2PjAN+GyxuQYzNczGbl+lmm0uQkQlbyXZGVYom+gzExqVZFQ0jtHPSTdAY541xGZoS4zLYzos6cgKamRdVUV/ao6ovk45xlin7tiUbu0jVuGwELo6VbaH8pJtVOQO5bjs4q6IHBs23ebGpm2la+Hkk1ofFys8z5fuV3GcSC5He3ROQjrPrY8u7oEVV4rofEpMj3JiWyaDogg8/tul62pUvfWlLlq7EGSWjwSL8xIdT2WIS0vLyoOXtLkVa6MMRCN9E3in7Zsv7mYkMyRlCWrjONn/vk2PdQyh2szZCvgaLKja5QGNst6Gmqk1NwUVe+MoJaGZe2ER9aY+qvkw6xqOQN+v4fqZ/kKgalxXIcPSTgAOR3rGtJM/y7oOpyIRt62vaf5RB821eqh53Hi2chvTCXpOw/mbsPh43x+zrbPP/QUhv73CkThe0qEpcj0Huo+Kvv/TFoOiCKz+29XrahS99aUseXQEZNRMftd3Dvkx8Dn2VJSPjHGAMuNLydk9BnBsyCfg58gyNzQvbYZKHwY9E6gSmbN8EG4s0BI2Qr8Eij00+0RjbJ49NTcdFXvjKCciOQUByToR2umoEtoX60h55fBnQ7mNsI1XjArASmY39VbPeMvtmFmIZ8HVkVE6dDKJv81D1uIfJPrcfiFwLJq1vk12RY4kP9b4aeY69S1SJ6/n0vh7SN4OgC678OEw7r6dd+NKHthTRlen0zotXGysQZ55epxGOOQeZGGhy3YbUhMZYSWLQ80Jzwh5d8GUXjrGNaFzc0VXfdvW4Bx2Nqx3Uj/ZQXyqF2ESzZ6ZVqqMxVuJoTtijC77swjG2EY2LO7rq264e96CjcbWD+tEe6ktFURRFURRFURRFURRFURRFURRFURRFURRFURRFURRFURRFURRFURRFURRFURRFURRF6cP/AdoKCcAwJompAAAAAElFTkSuQmCC\n",
|
||
"text/latex": [
|
||
"$$\\left [ - \\gamma {u}_{k - 1,i + 1} - \\gamma {u}_{k - 1,i - 1} + 2 \\gamma {u}_{k - 1,i} - {u}_{k - 1,i} + {u}_{k,i}, \\quad - H {u}_{k,0} + \\frac{- {u}_{k,0} + {u}_{k,1}}{h_{x}}, \\quad H {u}_{k,I} + \\frac{- {u}_{k,I - 1} + {u}_{k,I}}{h_{x}}\\right ]$$"
|
||
],
|
||
"text/plain": [
|
||
"⎡ \n",
|
||
"⎢-γ⋅u[k - 1, i + 1] - γ⋅u[k - 1, i - 1] + 2⋅γ⋅u[k - 1, i] - u[k - 1, i] + u[k,\n",
|
||
"⎣ \n",
|
||
"\n",
|
||
" -u[k, 0] + u[k, 1] -u[k, I - 1] + u[k, I]⎤\n",
|
||
" i], -H⋅u[k, 0] + ──────────────────, H⋅u[k, I] + ──────────────────────⎥\n",
|
||
" hₓ hₓ ⎦"
|
||
]
|
||
},
|
||
"execution_count": 91,
|
||
"metadata": {},
|
||
"output_type": "execute_result"
|
||
}
|
||
],
|
||
"source": [
|
||
"eqs = [\n",
|
||
" (u[k, i] - u[k - 1, i]) / h_t - D * (u[k - 1, i + 1] - 2 * u[k - 1, i] + u[k - 1, i - 1]) / (h_x**2),\n",
|
||
" (u[k, 1] - u[k, 0]) / h_x - H * u[k, 0],\n",
|
||
" (u[k, I] - u[k, I - 1]) / h_x + H * u[k, I],\n",
|
||
"]\n",
|
||
"eqs[0] = (eqs[0] * h_t).expand().subs({gamma_val:gamma})\n",
|
||
"eqs"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"metadata": {
|
||
"hidden": true
|
||
},
|
||
"source": [
|
||
"### $i=1,I-1$"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 92,
|
||
"metadata": {
|
||
"hidden": true
|
||
},
|
||
"outputs": [
|
||
{
|
||
"data": {
|
||
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAY0AAAAZCAYAAADNJ3IwAAAABHNCSVQICAgIfAhkiAAABsVJREFUeJztnWuoVFUUgD/N0usVfCQhJWWWj0Ty9lCUSq5l9aOXWdGPgoaKErKQiB5SdBVJocAselARt0ioSKyMsIIUtah8JGn5SO0qpWWRWaKl1u3H2oNzj3tm9jlnn8fMrA+Gc2c/18za66z9WGcuKIqiKIqiKA3Na8BeoDlrQeqEC4BO4M6sBVEURfHNWOA/4P4E2j4ZuXEuBrYBh4D9wCrgDqB7An0mwY3As8BK4E/EIbxRpc5iYA/QJ1nRFEVR0uVj4A+gKYG2pyE32N3AQmAu8KrprxN4B+iWQL++WY/I+xewCTenMc6Um5msaIrilwIycFuzFUNJkALRdTwcWWW85FGeUi4FruH4FcUgYBci9w0J9V1KgXh2MAkYhji4VtycBoiD2UntrKi8U/rBrweOAF8CZ5QpPwv5cu8CHjF/T7WUO93kLfYmqRKXsPoF1XEtcjtyI3wrofY/BZYgjqmUn4EXzd+tgbwoYy9plgHfmz7D8CYy9i+35DWEvZQ6jb3AJ8gS7EFL2aEmfTXwCnIwBLDGUnasua7zI6bigbD6BdVxLTIZ+Bf4IoO+j5jr0UB6lLGXVz4zV5vTaEh76QkcQGYEQT5ABuOF5v0PwK9l2pmLeNWrfQuoxFqWh9EvqI6zokA0HTcjN+wNnuVxoYfptxO40pIfduxVo4C/bdpW3Len+pqyX1nyGsJegvty/wDfASMD6dcCVyH7pGuAAcAQYG2ZdovKD3rVeciMI888Rx0sIcvgql9QHdcipwEnIBE+aTMPGA18CHxkyQ8z9vLMfuBvZLuplCj2Ugu2Ag728jriEQeb903ADsSD9jdpk02ZOWXa+B34xZK+FHgynLyhmAi8D/yEyFeI0EZ//IbUxZGpw9RxfbU7tOmiX1AdpyVTB/50PMGUqXSeEbY/l9n3fabsJuTmWQ7XsRdX5nYHmUtpxf2zgug5uAUXxV6SthVIwF56WApsNtdRwI/I4c6ZSAz2PpNX3LuzedWhppOllrwWZOAkRR9go+kjaj/7qhcBZGB2AG0JyvQ00C+Q1gJchzy81RHIW+/Qpot+QXUchrzo+JC59qpQZjsyU3Zld5X86cACZBVxGXJzLIfr2AuShB3EoYlj33WRKPaStK1ASvYyFfFIM4CzkAH2OV1jr980ZWyRENNN3hOB9EEmfZR532zaWYcs63xzgPBedTAiY3AJbaOd6g4jSBSZghSIt5frol9QHUclSx2fauqtitm/KzNMfxuAUxzKu449Fwpkc6bRHYkc2x5ID2svadsKeLIXW6xx6WxgAbIauYeuoWkjkUiJnYG6PYG7zd/Bve4WxDtvAUYgB0lHgYs4fqaQFWOAg8DWrAVJEBf9guq4FtmDbPWMSKGvh4D5yKx+EhIhVQ3XsZdnRiBOLriaCWsvtWArYLEXm9PYhkQy3IQcUL0AfB0ocxg4EXk4pkgz8mToaPM+WKcFmZFMQWYXLwO3cvwyL0uKMgZj0OsJF/2C6rgW6QRWAAOBsxPs5zHkEHctsiX1m2M917GXZ8ab67JAelh7qQVbAYu92JzGYSR0rB8ye3jUUqYYHbECeB7ZK9uKHILtQSIMdlg6H4Z8iVORfcpyzKH6YVdrhfpRGUP5/dCZyPKu+LrFknZJAjL5xkW/0Jg6zkomnywyV1vYqw9uA2YjN/+VyCF4W+BVKFPXdeylwRRki7kdeNikTShJe6pMvSuQz/5eID2svdSCrUBle+nCEiqftPdClpe7kaXLGmQJ1g/xSMstdTYjD+4cRJ4QrcRAZLlX6dW7ShtR9u+2Ir+tY2MAMnsrvhYBzwTSqv3WTx7ONKC6fqExdZyVTEEKRNfxSUiUju15CB+0Uf0GtrxCfZex50KBeHbQRuXP0GGp0xdZCbxryQtrL2nbCvi3ly6sQrypiyAu9DbtjUVm6AeA8z21XY6wX1AzIuP4agUN7WRzEO4D3/qF+tRxFPKg4+LPWZyXsRw2khh7aXEv8r1eHLOdLGwFPNmLLeS2G7Ik2YJ4QR+ci3zZG5GfChiJzDjGIfHDvujDsb3c7sgDOC1IGOAuBxkBvvEoT1yZkiAJ/YLqOE86no/MDmcjPy6YF5Iae2nQhDjjRcSPTkvLViAlexmOfKCF0WS0Mg158KdIN+Bt5CDN54yjleoP+xRM2hCLjJtxpx23lYaLTGmShH6hPnXsU6a0mQg8Tr7+CVNSYy8NzkHsfYiHttKyFUjJXm42DTwQQ9A8Mwv4FvsqqxGod/2C6jivNMLYq0XUXqqwmvxHwCjxUB0rijtqL4qiKIqiKIqiKIqiKIqiKIqi1B3/AyKO+5gdgZlpAAAAAElFTkSuQmCC\n",
|
||
"text/latex": [
|
||
"$$\\gamma {u}_{k - 1,i + 1} + \\gamma {u}_{k - 1,i - 1} + \\left(- 2 \\gamma + 1\\right) {u}_{k - 1,i}$$"
|
||
],
|
||
"text/plain": [
|
||
"γ⋅u[k - 1, i + 1] + γ⋅u[k - 1, i - 1] + (-2⋅γ + 1)⋅u[k - 1, i]"
|
||
]
|
||
},
|
||
"execution_count": 92,
|
||
"metadata": {},
|
||
"output_type": "execute_result"
|
||
}
|
||
],
|
||
"source": [
|
||
"Ω.solve(eqs[0], u[k, i])[0].expand().subs({gamma_val:gamma}).collect([\n",
|
||
" u[k-1, i-1],\n",
|
||
" u[k-1, i],\n",
|
||
" u[k-1, i+1],\n",
|
||
"])"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"metadata": {
|
||
"hidden": true
|
||
},
|
||
"source": [
|
||
"### $i=0$"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 93,
|
||
"metadata": {
|
||
"hidden": true
|
||
},
|
||
"outputs": [
|
||
{
|
||
"data": {
|
||
"image/png": "iVBORw0KGgoAAAANSUhEUgAAADYAAAAjCAYAAADfXvn1AAAABHNCSVQICAgIfAhkiAAAAuBJREFUWIXt2E+IV1UUwPHPmATKkFILBWVENFES/c0fkRaJgTuJ2thGCEFoYYtoIZS4KF0YLvy3F1xFyVSkaEEJQlgbJzVCzP5AVCBI6aAWaKSLc19zf2/e2G98v9/MLH5feLx7zzv33nPeufe9ey5dKtmEu3g8ky3GfaycFosSs2q278dl4Vwu+wtXa/Zdi7qONXCxJOvHd/gXL+JQi319jBsYrmkT2hOxSyVZ7uyaiucTcRiv1LTnP+o4NgcrjI/YOmPOFI49gU/w6kP6O4tbNexpYnaNtsvwGH7IZBvEx6NwdhX+xqd4G5/XGG/KWCjW0UupPoTvk6xXRPQGLuCZFvvcaAassWvYhaP4FW/gffyE21iNr9MY/9Qzc2axHW9hACMiigVnsKiizUZtilgnOYwXUnk7jqdyD34RUzXnC1wX/8Df8OwU2NhWVuHAdBvRpUuXevSk+/1ptaJL53lX/HdmLOUt1Rm8V6H3stgDzkv1qjxsqtiAE/hdLKFtVUplxwZwvkJvCD9iNNUbWs+zWuWYyAD+j16RyL4uModKcseWYb6JHRtJ5YVYII4DTuOO2Pg+34JR7eC02HwPi1lUSe7YYFL8pqTTY2wjS0QLXsNBrBVvcEZtk/JEc1A4OlEWmzs2KtbdtSQbxr5OGPio5BEbwClxjpFfu8UiLSLZwEljTsFysQYnwy6RtxXX1grZc5Pss5I/sLNCvl/zUdoV7CjpnNB8GnUO61P5qEhCyzwpXkhxfYgjJVk5tSlz2wRfxWIqLk0DjVTo5OtrLp4W6X5OPz7K6nvxJr4U6/ZgRb9/pqvgVqpPNvKVFFNxMN3LHw7C6MKxNen+bfb8Kc0HOPAZ+rDZ+OjWpVcsh4awvy+V+3Kl3LGfcbPUyRLNkWyIU6k7mU4/7okT4YJ1qd1oetZOhsSMuSCm6jupvKfN44xjkfh5LxUzYHWnB5wK5uArse2BLfhg+szp0qVLp3gArXKVkSmRbzYAAAAASUVORK5CYII=\n",
|
||
"text/latex": [
|
||
"$$\\frac{{u}_{k,1}}{H h_{x} + 1}$$"
|
||
],
|
||
"text/plain": [
|
||
"u[k, 1] \n",
|
||
"────────\n",
|
||
"H⋅hₓ + 1"
|
||
]
|
||
},
|
||
"execution_count": 93,
|
||
"metadata": {},
|
||
"output_type": "execute_result"
|
||
}
|
||
],
|
||
"source": [
|
||
"Ω.solve(eqs[1], u[k, 0])[0].expand().subs({gamma_val:gamma})"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"metadata": {
|
||
"hidden": true
|
||
},
|
||
"source": [
|
||
"### $i=I$"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 94,
|
||
"metadata": {
|
||
"hidden": true
|
||
},
|
||
"outputs": [
|
||
{
|
||
"data": {
|
||
"image/png": "iVBORw0KGgoAAAANSUhEUgAAADYAAAAjCAYAAADfXvn1AAAABHNCSVQICAgIfAhkiAAAAxZJREFUWIXt2FuIVlUUwPHfRAgTQ0U+FBgTQyUKYt9cJAoSg+olIhDyoUIEoYd8kAKhhoguUNFDZj2G0JNUTEVG0kNCFNZLkyYW4Q0EK6HroBVUND2sffj2nDnfbb75ZoS+PxzOPmuvvfdaZ6995X/KHfgLKzLZtZjFmmWxqE0uaZE/im+Ec7nsDxzvlVGLQSvHajhSko3iGP7FvXi5zbYexXMdWVfnXfyKqXYLtNNjX5VkubPrK/IbsU78kIWwB1s7KdDMsUGsNr/HNqg7Uzh2Od7DQ03q68axj3F+gWUrDZkVk0XBxiS7NX0fw1ocwp1N6hoQobSiiU4rNukgFC9tkveTcGICZ9P7tSQ7Knp0FfbhQXzdpK4R/GDuJEREQ5UNd+H71uY3pplj5zCJvXgVn+AN3I8LIiQ/F87906KdRmFY69DeJWE7HscYpjGU5R0UDhdM4sku29ukg1Dshj24J6W3462UHsAZEaoF+7C5i7Y+wo9i/TyLW7qoa8GsxUvL0XCfPn2WjoH0nl1WK/r0nhfE2nLRUt7dHxSLaZkt4vx1RfquOqctFRuxH9+JIbStSqns2Bi+qNCbwEnMpO+a9s9h7fI6nmpDb0jsO3fiz0ZKuWPX40qNHZtO6WtwtdipH8DvOIXb2zBqMTgg9p5TIooqyR0bT4pflnQG1De61HfkO7AbN4k/eFFto/Jjy7hwtNFJNXdsRoy7c0k2hed7YeBCyXtsDB+Ie478eUIM0qIna3hf3Sm4QYzBTpgU57rieaBCdluHdVbyM3ZVyF8096rtWzxc0tlv7m3VIdyc0nvxSEW9V4kfUjxv45WSbLCiXM4FDWbFIhRHUkPTFTr5+LoMN+JwSWcU72Tfz+IxfCrG7e6Ken9JT8H59N1pz1dShOJ4epcnDsLowrH16X00y18pLnzyde1DDONu83u3W4bEcKgJ+4dTejhXyh07jd9KlVxnbk/WcEJM8QWj+FvcGBdsSOVmUt5iMiEi5rAI1adT+plFbmceq8TiPSIiYF2vG1wKBvGZ2PbAfXhz+czp06dPr/gPRuyh+uS9hDQAAAAASUVORK5CYII=\n",
|
||
"text/latex": [
|
||
"$$\\frac{{u}_{k,I - 1}}{H h_{x} + 1}$$"
|
||
],
|
||
"text/plain": [
|
||
"u[k, I - 1]\n",
|
||
"───────────\n",
|
||
" H⋅hₓ + 1 "
|
||
]
|
||
},
|
||
"execution_count": 94,
|
||
"metadata": {},
|
||
"output_type": "execute_result"
|
||
}
|
||
],
|
||
"source": [
|
||
"Ω.solve(eqs[2], u[k, I])[0].expand().subs({gamma_val:gamma})"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"metadata": {
|
||
"heading_collapsed": true
|
||
},
|
||
"source": [
|
||
"## простейшая неявная"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 126,
|
||
"metadata": {
|
||
"hidden": true
|
||
},
|
||
"outputs": [
|
||
{
|
||
"data": {
|
||
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA94AAAAZCAYAAADHaPa6AAAABHNCSVQICAgIfAhkiAAADVdJREFUeJztnXusHUUdxz/lUfpCKDRCkFAoyKMh9tJaLYLNBfFBFLU80vgIngBBrNUgwQeIeNEGjBoFDFGBkAsRA0qD2higRL28rQUKiBR5yIXIUxAfVaSC1z9+s7l79+x7Z2dnz/4+ycm5d3Zmz8x3Zr87uzM7C4qiKIqiKIqiKIqiOGEYmAh9Hm40N4qiKIqidIUrgReA2U1nRKmNJUj/8pQK+/ic2cdHreRI8R31BSUPVb3Ftq/MY+o19URcpGGzYQwYAVZb+nFF6TK7IkZwPfAY8Arwd+B24GRgm+aypihKDaxHzqXDKXF+YOKc6iJDLWAp8D/gjKYzotTO9cCzwJyMeLsBrwPfi4T/CDl2DrSfNeeoV6SjvqAUIY+3uPKVWci19AgwTsaF94ilH1UUBU5DjqtngKuBC4ArgL+Z8OuAaY3lTlEU27wEvEb6CM0m5Phf7CRH/rMe8cSZTWdEqZ23IW3/7Ix4p5p4R0TCHwL+yWDctFavSEd9QSlCHm9pwlfG0AtvRSlEj+y70kkcCRxD/8G8O/CU2e9xFfKmdJMe5dukUh8LkHq5PyXOLOC/wH+A7V1kqmZ6VGuL+yOjWpdays8g0WMwj/PNwJOkd3JvAF4Etg2FzUZGq26vL2vOGHSv6KG+0CZ6DIbXZHlLE74yRujC2+UdwxWIgWwA5ifEOY+pU2rOMv8fGxN3L7PtervZHEhU++b4NbAOOYGEeQ6ZQgb9Rle0vrpQV10ooytUy/p4q/nekBJnCbAd8ABynAd0tV5OQmb9XFvjb3RVW1/LfY35/XcnbH8DctN6HdIhDhhC+q33AguR6aHPAVuAO4G315TfOlCvSEd9oV34omWat3jhKy4vvF8AbkamAnwhZvsCE74RuNyELTHfd8fEX2q+77WYx0FFtfeT4ET6WiS8aH11oa66UEZXqJb1kaczvcx8R/Xvar0chXSCflvjb3RVW1/LfYf5Trrw/gAwnf6OejDdej5Spp2Aq4BbgEOBXwA7Ws1pfahXpKO+0C580TLNW7zwle1s7SgHdyAjeS8xaThhLkYEWcXk6OASZErAUzHxg31sspvNgUS194/tgBPN3zdGthWtry7UVRfK6ArVsj4C/VbS/wxZQHD3PK4z3bV6mY2MNmwG/lXj73RRW/C33BvN9/KE7SuQ9rA+Eh507pcBhzO1I78WGW0bAm6zk81aUa9IRn2hffiiZZq3eOcrw7h5xvt3yKrOYT5ofvv7obBdTFj0oiTgZrN9j0j4N8y2MlzCYE8rqVP7KrqD39r3sP/sy7fNPn+ZEidPfXXhOGmijK4oq2WPcm2yC+2lKaYxuWhins9bQmnb3MZ7lPfH/U3aaEfIJk2cz2zRo7y2vpf7FWQ6Z5QZyCJHP43Z9gCS56Njtq0x25JG0X2iC17RQ32hTfRw6zVQn55x3tKkr4xRcnG1cfKbxAQyRz6Oq8z2Pc3/M4E/AX8B5obiHWXirUnYz1+B52PCbwS+lZAmi7lkv+KiCMuRKQpPI2Xp5Uw3jh2to9SpfRXdob3aj5bI22dN2s2IYSWRp758PU7GsdeGmyhjHsq2sTB1aDmasi9f20sTjGPXZ4PO4l0pceaaOP9m6uIuTdTLKuAJZOGme4B35kgzjl1/PNTES3uOs+hvRuupifOZD9q6LndRP3ya/sesQBYjnaD/fbozkMezHk/Y349NugX5stsobfOKXyH6JjGO+kIesnQsii/9XN/aZJy3NOkrY2YfQLGp5o8jJ5G8PJMQ/rD5Xgj8GXkgfx/kncYvh+IFQ//3xOxjAWJKcXdXhpALljK8nB2FUaTRjuSIOwd40OSnSJ5saR2lTu2r6A75tC9CWe0vBHaOhA0BHwKuROo+zH0F87UauAh5dcG7EENKIk99+Xqc2GzDTZQxD2XbWJg8Wtpsk762lyaw7bPBdLq059iCZ8nuY+riLq7rZSXiQ6uQlVxXIau9LiR+qmCAbX98xXzPSIlTtZ5cn8980dZ1uYv64Uwm6z/MscBW+meDLUL6rEmjoIuRGWJP5Mlsw7TJK4K8JF1UgfpCXrJ0LIov/Vzf2mSct3jpK8PIFflIHTsPcaz5ndOBfZED507632V8jYk3P2Yfq8228yPhu5vwheb/2WY/9wJ7Z+RrT5M26+Xpo5TTaAvlRsRsUpf2VXSH/NqXpar2PSR/wxXzcbrZz++BN+aIn6e+fD1ObOK6jGUo08aqaNmjXJvsQntpiu8gZT05Jc7nTZyLI+Gu62UDcFkk7FHggox0cfQo7497mLR1vh7K9fnMF21dlztMlh9ug6xPEh1l2hZ5TvSGmDSfMvk7JWbbjmZ/vwmFHQ+8ytTyX2R+c7eUvLmgTV6xL+XaYK9kOhhMXyirY16a7Of61CbjvMWmr0BxbxkjNOLtclXzgPAo3kXInYZPhzNlOBAZ/n8yEr4D8Enzd/Ru4RByl+OPwAHIM7KvAYfRfwcnyiJkSs8jOcrQVurSvoru0A3tvwh8F7lzeASycnkWeeqrC8eJ6zK6ogtaduHYDsgzinWI+Y4uluSyXqYjIxTRO/zrgXekpKuDZ5FHZw6o8Tdcns980tb1ebwIByA3kKMjacuBXYlfEyJtVO0Qs7/wtrXITe5zzP9nAh8B3kf8tFeXtMUrQHSfyMirbQbNF6AZHV3hU5uM8xabvgIVvaWJC+/HkGkzJwDvRxaJilvpbiuwPfDmUNhs4ArgYPN/NN0QIsaHkZHBy4CPEz+dKUqQNvq+5UGiLu2r6B5OP6jafwVZROIeZHr5iznT5amvLhwnrsvoii5omVTGUYo9k+Y72yAn6a3AH1LiBdNHo51pl/UyDxkBiHYQnkdGIVwyAdxq8rRfTb/h8nzmk7auz+NFCF6TFR1JWoF4xc9j0ixGRpkejNkWdJ7DHfEJ4GzEY74EnIucRx+NpB3FrRe1yStAtH0c+EeOuLYYNF+AZnR0hU9tMs5bbPoK5PeWWFy+TixgKzJXfj9k1O+chHg3Ie9+uxW5SzEHuWh5ALkbNgtZaCrMEFLxVyCrP99SIF+LiH+O4WzzCdgBEf3MUNjRtOP1FXVpX0V3SNZ+DfDljLRHINM4fOUTwNeQC+jbkIXVoowTv3hFnvry5Tipk7rL2FQ7G0QtoySVMbjpG7e4Uhs5ENFyE3LcxjEH0XELk7NZAlzXi0+sBY4D3ovcbLSN6/OZL/hc7vcg58RwR3ga0gm/i/4bF9ORDvz9yMhalKRRq/XI64XWIIsrbaQf117UNq9YTPxoYN0Mmi8k6TgI/Vyf2mTUW+rwFcjnLZkM4+YZb4B1ZN9hnIFMr30Gmap4NzJlYWfkzsVYTJqHgctN/BUF8/QIcFpM+C7IhU/wWYs8cxMOm5lj/z484w31aF9Fd0jWfh5ykkr7zMqx/yaffRkhe8XIsZT0WfXly3FSJ3WX0UY7K9PGqmjZo1yb9KW9bELuVs+N2dZGTkTq4/KUOIcxOZITxWW9TEcuMk6IhF9CuYutHtWeXZyOdIg2lEyfhcvzmU/auj6Ph0nzw52QUa2fRcKXImU9IyZNMLX0hwn73Iy8Kii6Xs2RyHt7X2dyBDmKay9qk1cAvIQ8b16UHuoLYZJ0HIR+ri9tMs5b6vAVyOctAWP0P9ILuL3wvh3JbJ7GlJdZZp9LgY8hjTBLjIDZJu2yrIi0e3E1sK99Fd2hmPZl8UX7MvhUXy7qygZV22QZiraxLmiZVMadTfg3LeWxi1Rt4xuASyNhj1BuATAbnIX0Pw7JilgzNrzDN23zYNsz0/zwM0hdHx4JP9+E71Phd8MsQlYjPgkZMLkpJk4XvKhK3e6D1MmR9WQtk0HxBRc6tqmfW1ebjPMW274C+bwlzBgZF97BJzq9xRbTkDsID1ne7zLkTnMw+vx15F1ub8qR9lDyX9yMkv/Cew4ynWIIuatzrvl7r5zpbVOH9lV0h2LaF8E37cvgW33VVVe2qdom81KljXVBy6QyHoOs0O/6mddBomobX4lMcz0FOAgZrdhC/Kq0LpiBLMyzrqHfD7DhHb5pmwcb5c7jhzORUbHrYtJvxt6jN/PN75xr/j8YGXkbjsTrghdVqdvjkeuB6KunXDEovlCXjm3t59bRJpO8xaavQH5vmUf/7NY+9kYuKIPPaosZDbO/ycDVlvd7GiJwwDTgJ8jc/KDj1zO/vXdM2rw3GkbJf+E9TPzU4tGc6W1Th/Z5dAc72hdhGL+0L0NT9dXDbV3ZJm+brMow2W2sx+Br2aPdZWwjVeolYBWyvsSrJt1y+9ksxHLgq8hMiaawoSv4p20WNso9TLYfHoT0n5L2YYNdkLJEp49eizzn2TWq1O0F9L/yzTWD4At16ThMO/u5dbRJ37xlFlOvqUdqzFcmK+lfnMwV5yGrSDaxqJwPqPbtoqn60rqyRxe07EIZ24jWSz10VdeulrsLaN1WRzW0i+o5IGykvpfXK+mo9u1B68oeXdCyC2VsI1ov9dBVXbta7i6gdVsd1dAuqqeiKIqiKIqiKIqiKIqiKIqiKIqiKIqiKIqiKIqiKIqiKIqiKIqiKIqiKIqiKIpSmf8Da8CXk/nuMgAAAAAASUVORK5CYII=\n",
|
||
"text/latex": [
|
||
"$$\\left [ - \\gamma {u}_{k,i + 1} - \\gamma {u}_{k,i - 1} + 2 \\gamma {u}_{k,i} - {u}_{k - 1,i} + {u}_{k,i}, \\quad - H {u}_{k,0} + \\frac{- {u}_{k,0} + {u}_{k,1}}{h_{x}}, \\quad H {u}_{k,I} + \\frac{- {u}_{k,I - 1} + {u}_{k,I}}{h_{x}}\\right ]$$"
|
||
],
|
||
"text/plain": [
|
||
"⎡ \n",
|
||
"⎢-γ⋅u[k, i + 1] - γ⋅u[k, i - 1] + 2⋅γ⋅u[k, i] - u[k - 1, i] + u[k, i], -H⋅u[k,\n",
|
||
"⎣ \n",
|
||
"\n",
|
||
" -u[k, 0] + u[k, 1] -u[k, I - 1] + u[k, I]⎤\n",
|
||
" 0] + ──────────────────, H⋅u[k, I] + ──────────────────────⎥\n",
|
||
" hₓ hₓ ⎦"
|
||
]
|
||
},
|
||
"execution_count": 126,
|
||
"metadata": {},
|
||
"output_type": "execute_result"
|
||
}
|
||
],
|
||
"source": [
|
||
"eqs = [\n",
|
||
" (u[k, i] - u[k - 1, i]) / h_t - D * (u[k, i + 1] - 2 * u[k, i] + u[k, i - 1]) / (h_x**2),\n",
|
||
" (u[k, 1] - u[k, 0]) / h_x - H * u[k, 0],\n",
|
||
" (u[k, I] - u[k, I - 1]) / h_x + H * u[k, I],\n",
|
||
"]\n",
|
||
"eqs[0] = (eqs[0] * h_t).expand().subs({gamma_val:gamma})\n",
|
||
"eqs"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"metadata": {
|
||
"hidden": true
|
||
},
|
||
"source": [
|
||
"### $i=1,I-1$"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 127,
|
||
"metadata": {
|
||
"hidden": true
|
||
},
|
||
"outputs": [
|
||
{
|
||
"data": {
|
||
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYEAAAAZCAYAAADXG5K+AAAABHNCSVQICAgIfAhkiAAABx5JREFUeJztnXuIVFUcxz+apesupCYhKWWWj0RwyzaUQsayIiIzK/qjqKGipCwkosheWlJCf5RFFBWxSYKFYiWEFqTpGpWPpIeaaY1RmrZkVmppZX/8zs2ZM2fmnvu+u3M+MNy75zXnfM/v3PO8s+BwOBwOh8PR4LwG7AWas85IDhgHHAVuzTojDofDkQZtwL/APZr7SciDcCmwHTgE7Ac6gFuAninmMQrXAM8Ba4DfkAf86z5xlgK7gZZks+ZwOBzZ8x7wK9CkuU9HHpi7gIXAk8CrKuxRYDHQI71shmYTkt/fgS3YdQLnqXCzks2aw9F9KCKNppBtNhwBGYHMAl4y+F0IXEH1iH8Q8D1S31cnmrtjFAlvX5OA4UiHVcCuEwDpMHbSdWY8OkVcm/SlvHKvAo4AnwCn1Qg/BxH1NuABdT/NEO5U5bc0tpx2X4LqDk77OLkZeTi+YfD7AFiGdBLl/AS8qO4Lml+Y+kyalcA36juDsAixp4sNfo1gg41QxopOYC/wPjINvM8QdphyXwe8gmweAaw3hG1T143xZLNbE1R3cNrHyWTgH+DjgPGOqOvfmnuY+swra9XV1Ak0gg02QhkrOoG1yCjmAHCuIeyzwAnAHcjIaBzQiUyLdbz4n8WW0+5LUN3BaR8XzUArsuxxIEC8XsCN6n655hemPvPKOnWdaPBrBBtshDJWrfX9BWwGRmnuU4DLkXXT9cAAYCiwoUa6nkB6LzkPGSWF4Xm6wdSrBra6g9M+TgYDxyGnYIIwDxgDvAusMPgHqc88sx/4E1n6KCeMDUaxP0jfBrMoY1r4arkAWesaov5uAr4Ffgb6K7fJKszcGmn8AuwxuC8HngqW3//pT7zH1SYC7wA/ImUpBohbUnFsP+0WadroDk77EsG0r7cBOkGFMe0H1OJuFWcL8qCohW19migRv315FLDfGAapI33JK4wNRrE/8LfBEvFqlkUZbYjSdjwqtOxlCLBVXUcDPyCbI6cj56L3KT9vrczUSw5TX6JPk0Gm3gsCZ1nY5x8EkMotAbN9wrUAX6r8BM3TM0A/za0VuBJ56aik+W2ySNNGd8i39rZE0X4HMjq1ZVcdv0Pq2scyrRnAfGSUfxHyIKiFbX2aSMK+wtLEMZ08wthgFPuD9DXLoow2RGk7Hr7teRrSw8wEzkAa3EdUnodepMKYTj/MUH5PaO6DlPto9XezSmcjMu2qxxAVV59em2jHvwPQ+YNwPWo5RaIdR7PRHfKtfRji0D4spyBl67AIO1OF/QI42SK8bX3aUiSe444F7GcCPZF9ix2ae1AbjGJ/EN4Gi4TXLO0yhiFM26nS0nT+t3wEMx+ZLdxJ5fGyUcjpiJ1a3N7A7epeX5NuRUYUXwMjgU+Raeb5VPfSOmOBg8A2n3BdGRvdwWkfJ7uR5ZmRPuHuB55GRo+TkBNAftjWZ54ZiXRa+qg5qA1GsT/IxgbTLmNaVGlp6gS2I0fmrkU2sV6gegf8MHA88gKKRzPyNuUY9bcepxUZRU1FRkQvAzdQPdU04cXN+2mKKNjoDk77ODkKrAYGAmfWCPMwsuG3AVkC6rRM27Y+88x4dV2puQe1wSj2Vx4/TRtMu4xpUaWlaU/gMPAd0ij2Ag8ZwqxAzsmuRnaZW5AG8jkyuuqLbILpXz4cEXEK8GGAjI+l9hreLCpfbe+NNO57y9wuQ343Jc/Y6A750n4u8KBP/EnAqgDflzZLkLd+L0Ue3OXcBDyGPMzXIJvCOiXMm4y29ZkGU9UHZNkCZFO8Xd13UtlePC5Byv625h7UBqPYH9S3waRIsoxZthtrLZdRf+e5DzLF3YVMLdYjU6R+SA+zyhBnK/JyzEHkHHUQtiG/42JiANLQvM8S5Cx2uZv+mzA6edgTAH/dIV/aD0SmzfU+fX3Sz3JPAOTM/h7kDV+d2fifMllVJ22b+rSlSHj7mk39MpQMcU5ERrFvGfyC2mAU+4P6NliPIuE1S7KMcbQbCNd2rLXsQEYANhmxoa9Krw24Hsn8OZZxm1Xc8X4BFe1kszEcB3HrDulqH4Y8aO/9PMDZMaebRH2mxV2IJhdETCeK/UE6NhiVqGUMS9C2Y61lD+TXBjeHypaZ8cgmiTcifxw55zrYIu4EgjWkduw6gRZkCteK9N6PqHv9xZi0SEJ3SFd7W/KmfR9kA3BZjGkmVZ9p0ISMgBfHkFYU+4PkbDBOopYxCFHajrWWI5ARwMLQ2axmOvJyjUcP4E1ks83LUFF971BD3K3Y045dJ1Ag+os3cZKE7pCu9rYUyJf2IC/hPEp8/1QmqfpMg7OQNjQ0hrSi2J8XPwkbjBObMsZFAf+2UyRie76O6o3VNJgDfIV5s7oRyEp3cNonQZb12dVw9hcvXVbPdbjf/c4Kp70jS5z9xYvT0+FwOBwOh8PhcDgcDofD4XA4HA5HFf8B2Sn0Vmz2OcMAAAAASUVORK5CYII=\n",
|
||
"text/latex": [
|
||
"$$\\gamma {u}_{k,i + 1} + \\gamma {u}_{k,i - 1} - \\left(2 \\gamma + 1\\right) {u}_{k,i} + {u}_{k - 1,i}$$"
|
||
],
|
||
"text/plain": [
|
||
"γ⋅u[k, i + 1] + γ⋅u[k, i - 1] - (2⋅γ + 1)⋅u[k, i] + u[k - 1, i]"
|
||
]
|
||
},
|
||
"execution_count": 127,
|
||
"metadata": {},
|
||
"output_type": "execute_result"
|
||
}
|
||
],
|
||
"source": [
|
||
"-eqs[0].collect([\n",
|
||
" u[k, i-1],\n",
|
||
" u[k, i],\n",
|
||
" u[k, i+1],\n",
|
||
"])"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"metadata": {
|
||
"hidden": true
|
||
},
|
||
"source": [
|
||
"### $i=0$"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 128,
|
||
"metadata": {
|
||
"hidden": true
|
||
},
|
||
"outputs": [
|
||
{
|
||
"data": {
|
||
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAL0AAAAZCAYAAABkWi/IAAAABHNCSVQICAgIfAhkiAAABdxJREFUeJztm11sFUUUx39FrP0iKJJoJKGEWCmEhFsqWgSb+m1iVFQM8SPmBoko0QeNJmpM1WjUJ2M1JioG0fAAxgbUB6QvVUSlKl+KFkWwmEgF8QNBEBTw4T+bu507d+/edru9l+4vuZnuzJzZObszZ86c2UJCQgIAbwJ7geoI27wfOAHcGmGbxUYj0nHBUHckoTBmAMeBBxxlHeiltgTIv2Lq3GXlLzP59QPvYmzMBV4CPgb+Qv1flkdmJdAL1Axu1xKipAP4E6h0lP0G/EfwCrAJDY7pVv63wAFgRAR9jIvNSJcDQDfhBv0Fpt6jg9u1hKg4D1n51xxlE9HL3BIgXwX8C/wDnOrLrwaOAeui6WZBpMm/OuXiEqAOKDPyYQY9aILsorQm+LDBfinz0Qte4ah7vkm7AtprBEYCX6HB75Ey99oITEED5xfgIPApcGGhHY+JTmA7GuyFsBwYD1zhKHvEtHejo2y8KVtZ4P1KgaLR2x70lyOLvN5RN8ygbzLpl1a+5+rUmrLRwFvAR8BM4D1gVLgulwSfmNQ16BtNaj8j0H4KZBxONopG75G+v6uRRe4G/nbU9Qb9PLTsu/Astq2Yp3ATMJu+yrWj2Z9CG8aTgS9M2uwoawT2AT85yrxnvGkwOjXEFI3e/kE/DjgFRR5syshYa5f1ssll6dNkz+Zuk1aEaLdU2I/2NeOt/DHABGBNDjnv5dvP6Dk0aMI8+2KkqPT2uzdnmvQPR7065JKsRxPA9Rtj6h4GvvHJVgCTgZ3AakfbE026o/DuZ9GDfEP/7w1T1ukoWxrBPXPxOzDWyvMmv2uJB73gvcBuKz+FIkn9YRHwI5qEG4CLQ8j0kP2sgn75Nvdx692MXOafTf/S/kK/pT9sUpfFzTUT/XiKbUb7Ao9p5j4dAXL70YsZKC8Ap1t5KeB6dODWY5X1dyCFoZLMM/Xw3LwNjvoTgTOADxxlKbQHKpR5QBsa+OtMuhoFE1xuhscONEnCYg9Wm7j1rgG2GtlA+XPQrHCFFZ83ZXcGyD9k6rxo5d9D7lPKUShE2mnlzwWOoI2vRxt6GWcF9MFFmv6HLP20ED5kOQLpZa9ey00btVkScK8pe8bKP9vkTzHX1aadjchlCKILWGzlbQeezSMXNXHr7ecglqX3uze9wK/AJIdgGEvfYNJcm1jXLG9ArpFd1g58DTxmrh8EbgGuBvYE9KFYmIT0sleSehTK3WXlnwYsNH/bzziFVozvTLufowPCWWSvXH7K0bO3V9gO4KJ8CkRMnHrnxT/oTwBrkR96rlWnAThKX1/dJpffNh1Z7a0OGW9C2Ap7J5pp4GGgFbgGWalSwAvd2ivYUXRoV+fLqwaWAFPNtR3BSCEDMAedaSwGbifbdbIZiwITtpHYg6xonMSpd15GWtftwE3AVcAPJq8e+UibTOdd1CCFDgLbfPnlSKkt9D2s8ghaBTpQ6O9p4FoyYcA4mWN+kBkoM8lsgPehVcjmSrSvedfKX4Ni0mvRQUwNcBk6zOtFJ9o7LZkUerZLgOvQ2UapUdR6lyNL4D+AugNZ3tcD5GaRWSn8eF8dvppDrht911LmKLsUnRccI/s7nkJI03+f/gmCoxY9DpnRyBqtcpRVoL3JbuAQWhUXos33ceBDh8w29OwPATcU0Pdy5A7cbOW/TPwTJ069bbJ8ehfecXFDvoqDyDQU0ZmPVp9c8d1i5D70/GZH0FYVmvQzgNvQCyzEAHSR/R3V98S/kS2UgertJ9Sgr0Abjvf7eZOBUossQqu5noqsQcsQ9acQKlHf34movSZkrb0vXp9CsedxIeXnIZd0AToraUODwBVFKSYGqncNco9SaKVoNX/bh4V9aAYeJ9p/IgnDGOTy2O7QCuCzmPvSHyYjl2hCRO3dTebEGuQGvo32QFUmL41Wllz3XITcsCNGzvVpRLExUL1bcLujSwehrwlDwJMoqmYHJU52hqveCSiq1TLUnRgChqveCQkJCQkJCQkJCcOe/wHOT5SuMUTYXgAAAABJRU5ErkJggg==\n",
|
||
"text/latex": [
|
||
"$$\\left(H h_{x} + 1\\right) {u}_{k,0} - {u}_{k,1}$$"
|
||
],
|
||
"text/plain": [
|
||
"(H⋅hₓ + 1)⋅u[k, 0] - u[k, 1]"
|
||
]
|
||
},
|
||
"execution_count": 128,
|
||
"metadata": {},
|
||
"output_type": "execute_result"
|
||
}
|
||
],
|
||
"source": [
|
||
"(-eqs[1] * h_x).expand().collect([\n",
|
||
" u[k, 0],\n",
|
||
" u[k, 1],\n",
|
||
"])"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"metadata": {
|
||
"hidden": true
|
||
},
|
||
"source": [
|
||
"### $i=I$"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 129,
|
||
"metadata": {
|
||
"hidden": true
|
||
},
|
||
"outputs": [
|
||
{
|
||
"data": {
|
||
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAM4AAAAZCAYAAAB5JBFTAAAABHNCSVQICAgIfAhkiAAABb9JREFUeJztm39snVMYxz+dmq7tMmYJsVin1LZmyW5bo7KlaiEh/jDSRIbIzYgR/EFIECmyBX+I+JGIX6FkoRPLhD8wkfltjG0YJTI6yTrG0G02Zlv98T1v+vbcc9/7vnfvvX1v+36Sm9N7fr3Pc+4553nOc95CSkpKbDwP7ATqYuzzZmAIuDzGPpNGG9LxmtEWJKX8zAcOA7c4ytaiidEZ0P4JU+daK3+lyZ995CKWjS7gMeADYDeSf2WBNmuAHUB9aUVLSRprgb+ASY6yXcBBgi3RJjTBWq38b4E9wIQYZCwXm5Eue4A+wi2cM029O0srWkqSOB1Zm6ccZY1oQnwZ0L4W+A/4Bzjal18HHAI+jEfMSGQpbCXzcS7QBFSZ9mEWDmiRbaOyNomUCNg/7FI0SVY56p5h0k8D+msDqoGv0ALyyJhnbQSa0eT7BdgLfAycFVXwMrEO+AEtmCj0AjOA8x1ld5j+LnWUzTBlayI+L2mMeR3thXMesgzrHXXDLJx2k35u5XtuW4MpmwK8ALwHnA28BkwOJ3JF8JFJXQunzaT2GIHOl6ANppIZ8zpW+/6uQ5ahD/jbUddbOJchF8aFZznsAfMGsh1YyMhBW412pgw6hI8FNpi0w1HWBvwO/Owo88Z4UymEKiNjXkf/wpkOHIUiQjZVDFsN1y5qk8/iZMndafpMWhOi30phEJ3zZlj5U4GZwFt52nmTyh6jB9BkDDP2o8140HGEq3a8Sf901GtC7tV6tIhcn6mm7n7gG1/bGmAO8CPwhqPvRpNujS5+Dv3If/Z/njNl6xxlPTE8Mx9/ANOsPG8DcbkwoImzExiw8jMowlcM7wAvBpT3kzsuQZ9CwZEk6hiVDnR82I50ztoV/BZnv0ldO3++XcKPN2Cb0TnJY555ztqAdoPATwF9h+Vh4FgrLwNcjC51+62yYn+oMExieEw9PJf1C0f9RuA44E1HWQadCYuhFVgRUL4VWcew2BPeJok6RqUe2ILkKSjTSWh1uULGD5myqwPa32bqPGrlX0/+2/TJKPy9zsrvAv5FwQSPR9CPfEKADC6yFB+O9tNJ+HD0BKSXbUV7TR8NOS3gRlN2n5V/oslvNt/rTD8bkUsUxKnEo3sUxpqOe3FYHL+rtgP4DZjlaBzG4rSYNF9gwLUDtSA3zy5bDXwN3GW+3wosAS4Afg2QISnMQnrZFm02CtNvs/KPAZaZv+0xziDL9b3p9zN0Cb2AXAtq470CVM4I1njQccTCGQLeR375aVadFuAAI88uNvl821ZkPbY42niLylbau3nPArcD3cBF6E6lEvDC8rYlPYAuhpt8eXXAs8Bc892ONmXQJrIY3Xk9DVxJrhvoog1Zvd1hBY+B8aBjDkvQpL3Bl9dM4RVdj8419is1E9FAbnA1Yvj9NZeVAw3iQeDCQoIHkKV4U74YBRB6kF8+hH4kL+/BPO1eQnKfbOUvN33sAB5H/vN2FDQZQK86VVltelGgYRA4J6L8b5v25SRJOq6gcLCjs0D/TlfNZiJyhfyXnFeZBzwT0G4BwxbLj2dGn8zTrg8tNnsgARah+6RD5L73FoUsxS+cewge9H5Hmylop3zVUVaDzmoDwD5knZehgMZh4F1Hm+/Q2O8DLoko/y509iwnSdJxGnIdgz61BfoPtXBg+HWJlkIVS8g8tPssReedfHcCSeQmNH4LY+irFm0c84Er0I8YdhM5xcixKAY5SknSdQy9cGrQwe71EgoTRAParbrN97lop+ocJXmiMAnJ/kpM/bUjl897U305cnumh2jbhSaVHZ5PGknUsR6duzLICnabv+0L7Rw6gLuJ9x/ZwjAVuW+2a7cK+KTMshTDHOTezYypv+sYfrMC5NK+jKKQnouRRZPHfub9xHOpXGqSqGMnbte8pwTPShkl7kXRzupCFSuY8aBjSpnZQGW4sUfCeNAxJSUlJSUlJSUlJSVl1PgfCDirCadaf2EAAAAASUVORK5CYII=\n",
|
||
"text/latex": [
|
||
"$$\\left(H h_{x} + 1\\right) {u}_{k,I} - {u}_{k,I - 1}$$"
|
||
],
|
||
"text/plain": [
|
||
"(H⋅hₓ + 1)⋅u[k, I] - u[k, I - 1]"
|
||
]
|
||
},
|
||
"execution_count": 129,
|
||
"metadata": {},
|
||
"output_type": "execute_result"
|
||
}
|
||
],
|
||
"source": [
|
||
"(eqs[2] * h_x).expand().collect([\n",
|
||
" u[k, I-1],\n",
|
||
" u[k, I],\n",
|
||
"])"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"metadata": {
|
||
"heading_collapsed": true
|
||
},
|
||
"source": [
|
||
"## КН"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 8,
|
||
"metadata": {
|
||
"hidden": true
|
||
},
|
||
"outputs": [
|
||
{
|
||
"data": {
|
||
"image/png": "iVBORw0KGgoAAAANSUhEUgAABVUAAAAZCAYAAAAxOM4PAAAABHNCSVQICAgIfAhkiAAADtJJREFUeJztnXusHVUVh78WKH0hLTRCkFAoCKUh9tJaLYLN5eGDaNXySCMaPAGCWKtBgg8Q8aIEjBoVDFGBkAsRA0pTtTFAjXqRAtYCBUSKPORC5CmIShGoSP1j7cmdO3fOvM7es/ecWV9ycu7ds2dmzVprfufMPmv2gKIoiqIoiqIoiqIoiqIoilKJQWB77PWgV2sURVEURVEURVGUpnI18Bwww7chijMWI2MHp/Wwjc+ZbZxkxSIldFQXlCL0qi22dWUO48dLt6d1GjQLRoAhYLWlnSuK0j/sjgjbWuAR4BXgX8AG4FRgsj/TFEVRVKMsoX70i/o/PNYj10mDGX1+aPqcXodBDWAJ8AZwlm9DFOesBZ4GZub02wP4H/D9RPuPkXNnvn3Take1IhvVBaUMRbSlLl2ZjoyTDgGj5AyqDlnaqaIo/ccZiE48BVwLXAxcBfzTtN8ATPJmnaIobUc1yg7qR7+o/8PjBeB1siurNiPxWVSLReGzHsnZab4NUZzzDiT3z83pd7rpd2Si/QHgJfrjByPVimxUF5QyFNEWH7oygg6qKkqr6ZD/C2o3jgKWM1Gc9gSeMNs9vgfbFEVRoLpOqUbZQf3YOx30s7ZfmIf4/N6MPtOB/wKvAjvVYZRjOlTPX4ADkWq0yy3Z00906M23obIFeJzsAYwbgeeBHWJtM5Aqsw3uTKuNfteKDqoLTaJDf2hNnrb40JURYoOqdf4atAIRkI3A3C59LmB8Kfw55v/jUvruY5attWum0gMa4/7nt8A65AMxzjPIrSwwUbg1L/xR1vfqdztozvtDNcoOVfwI/e+Xuqjif/W9O95u3jdm9FkM7Ajch+hJRFvjcgpSTX29w3201behHvd1Zv/v6bL8TcgPRuuQwY6IAWRM4m5gAXLL7jPAVuB24J2O7HWBakU2qgvNIhRfZmlLELpS56Dqc8CvkRLeL6Qsn2faNwFXmrbF5v3OlP5LzPvdFm1UekNj3G6iLwavJ9o1L/xR1vfqdztozoeJapQduvkR2u2Xuujmf/W9O4oMlCw170n/tzUuxyAXuH9wuI+2+jbU477NvHcbVP0gMIWJgzDRLfBzkWPaFbgGuAU4DPglsItVS92hWpGN6kKzCMWXWdoShK7saGtDBbgNqQZ5gTHBiXMp4pBVjP0yvxgp5X0ipX+0jc12zVR6QGPcXnYETjZ/35RYpnnhj7K+V7/bQXM+PFSj7JDlR2ivX+oiy//qe3dE/lvJxDnbIqKql7SBkrbFZQZSJbQFeNnhftroWwj3uDeZ92Vdlq9A8mF9oj0auFkKHMH4QZo1SJXcAHCrHTOdolrRHdWF5hGKL7O0JThdGaSeOVX/iDzBNM6HzL5/EGvbzbSlfWkHqSrZDuyVaP+GWRYyl9HfZeca47Bj3MH+/CrfNtv8VUYfl3nRhJwAf3lRxPc+zseQzxMbhKyFofu+g12d8q1R0B/xyvKjavd4OtT3WdvmnHTNJMYeEFbk9bbYuk3+ntuhev4eaNZNXuTapMl606G6b0M/7leQW2yTTEUeGPOzlGX3ITYfm7LsQrOsW/VrSLRBKzqoLjSJDvVqDbjzZ5q2+NSVESo+qGqU4iKxHZm3II1rzPK9zf/TgL8Cfwdmx/odY/pd2GU7/wCeTWm/CfhWl3VssAwpF34Ssa9TYRuzgZkB2DSKnZgm0Rg3N8bDFWz7rFl3CyLA3XCZF65zAurJi1H8nZM+zsdQzhNw4/uQtbDJvh8uaVsIGgVu4zWKG+2Ik+dH1W5/n7Uh5qQvRrF7LkQDAXdk9Jlt+vyH8Q/K8BGXVcBjyENw7gLeXWCdUezm72GmX9a8iWX3mYyTD70Jwbd1H3dZ/XyS9Klhlpv1T0q0T0WmNHm0y/Z+YtabV8xcrzRNK36D+Lcbo6guFCHPj2UJ5Zo/tJxM0xafujJitgGUu/3/UeRDpChPdWl/0LwvAP6GTIC7H3Aq8GKsX1Sye1fKNuYhopQ2cj6AXBi5YiZwv9lH1f28mN8FkGQfJb96uKpNtmKaRGNcPMZFqWrT94BZibYB4MPA1Uh+xbmnpF2rgUuAB4CjEYHthsu8cJ0TUE9e+DwnfZyPoWghuPF9yFoYikaBW50KRaPAbbxcaUdEET+2Wbt9f9aGmJO+sH0uRLc4Zs0bF83ddg/jH5RRd1xWInmyCnni8SrkqcgLSL99M8J2/r5i3qdm9Ok1TnXrTSi+rfu4y+rnNMbiH+c4YBsTq+wXIuMR3aoXFyF3kDxWxFjPNEkrIlu6DZiB6kJR8vxYllCu+UPLyTRtCVJXBpHR1iEXG49xnNnPmcD+yIlzO1IyH+c6029uyjZWm2UXJdr3NO0LzP8zzHbuBvbt3fQJbKV8xcPeiI3zC/Qdpnw8qthkG41x8RhXodcYdxD7Bnu040yznT8Bby7Q31Ve1J0TEGZeZFHE93Wfj6qFY4Ts+yrY8H2H3nUqFI2CsOOVR1E/qnaPp0N9n7Vty8k6+Q5yrKdm9Pm86XNpor3uuGwErki0PQxcnLNeGh2q5+9eZt0NFdYtSt16E4pvfepsnn5ORuYdT1aH7YDMy3hjyjqfMvadlrJsF7O938XaTgBeY/zxX2L2uUeGbXXQJK3Yn2o52Km4HvSnLlT1Y1F8XvOHlJNp2mJTV6C8towQq1SdnNLBNfFKkEuQUeRPx40yzEfKdh9PtO8MfNL8nfwlaAAZwf4LcBAyP9rrwOFMHJ33xUKk5P8h34Y4RGPc/zH+IvBd5FeuI5GnZ+fhKi+akBPgNy+K+L7u87EN50moWtgG34ekUdDceJXxo2q3fYr6v005WTdFqs8ONe/JB8/UGZcpSGVRsjJnPfCujPVc8DQybcpBDvdRp96E5NuQdfYg5EfIZAXcMmB30udgzqqGO9RsL75sDfID03nm/7OBjwLvJ/1W5DppilaA+H17jq226TddAD9+rIuQcjJNW2zqCvSoLT4GVR9Byt1PBD6APPwh7alh24CdgLfG2mYAVwGHmP+T6w0gzvgIUl1yBfBx0m9D8EVk4xt5HRuMxri/Y/wVZBLqu5DbEJ8vuJ6rvGhCToDfvCji+7rPx34/TyBcLex334emUdDMeJX1o2q3Xcr4P5ScHKbcHHChMxm5ANsG/DmjX3RLb3KgpM64zEEqd5IXf88i1UN1sh34vbHpAEf7qFNvQvJtyDq71LwnK8BWIFrxi5R1FiHVYfenLIsGRuKDLNuBcxGN+RJwPvJ5/XBi3WHq1aImaQWIbx8F/l2gry36TRfAjx/rIqScTNMWm7oCxbUllTJzqtpiGzJ/wQHIL+7ndel3M7AEOfnWIvNLHI08yetpYDryAIk4A0jgr0Ke2HtLl21fCHw5x84jkbJe2yyk+xwW55pXxM5IgM+OtR0L3OrALptojLvH2KddNvgE8DVk4OFW5MEZSUZJn/zaVV4UzQkINy9cU8T3Ls/HNFQLxwjJ96pRzYqXK6r4UbXbHmX9H0pORsUaaQ+qaSLzEV9uRvQhjZmIH7cyVu0eUXdcQmINcDzwPuQHK9u41JuQCfm434toVnyQYxIywHIHEwelpyCDM/ciFXFJulWbrQc2IZ8Ly83fSerWoqZpxSLSq/hc02+60M2PTf8uDWHlZFJbXOgKFNOWXAapZ05VgHXk/3o0Fbn97inkFqM7kVLjWcio9EjKOg8CV5r+KzK2PQcRv6zX9ALHUWWei4eAM7os2w25mItea5A5V+Jt0xzY5AKNsTu7fM6vMkT+0wRHMtZ3kRdFcwLCzYs6yPO9y/MxDdXCMULyfQgaBdV1aojwNArcxcsVQ5T3o2r3eDrU91kbSk5uRqpMZpfcXqicjPj6yow+hzNWgZWkzrhMQQaQTky0X0a1gbQOvc0VOAW52N1Ycf08XOpNkpB8W+dxJ8nSz12RarSfJ9qXIMd6Vso60e2+P+qyzS3AS0yc1/wo4GVkkGVRciVD3VrUJK0AeAGZ37UsHVQX4nTzo8/vLHE61Kc14CYn07TFha5AMW2JGGHi9GBAvYOqGxBjiyRTUaabbS4BPoYkYZ4zeqVsos9AbFya19EwTDMfzgIa46IxrkIoMa6C7bzwkRMQZl7kEZLvVQt7o07fVyEU31ehjfHygWq3X1zk5CzT/k1LNraRXs+LjcDlibaHqPYwJRucg1xbHprX0TE29CY03xbBts5m6ednkFgfkWi/yLTv18N+4yxEntp9CvKD+80pfdqgRb3Edj8kJke5MS2XftGFOvzYpO/TrnIyTVts6woU05Y4I+QMqkavZFm6LSYho8MPWN7uUuRXxKh66evAk8BbLO9nJlLaPICMwp9v/t6nwLqHUe5ibZhiAwm92OQCjbHdC/JebQoFF3lRV05AmHlRlNB8r1rYG3X6viih+b4KbYqXb1S7/eIiJ5cDr1L/HJP9RK/nxUrk1uPTgIORKqOtpD+9uQ6mIg85Wedp/xE29CY03xbBxnEX0c9pSDXbDSnrb8He9ClzzX7ON/8fglTMDSb6tUGLeontCchYzyw3puXSL7rgyo9N/T7tIie7aYtNXYHi2jKHiXcNTWBf5KI1eq22aGicA40B11re7hmIgyMmAT9F5kuw+WV4kPTbsIZjfTqmbd8UG8sMVg9TbCChiE11ojG2TxGbQsdFXtSVE5Afgw7pORHZ6eqHqiL48n0H1UJfWthBNaosbYqXb1S7/aI5GSa9xCViFTLf7mtmvWX2zSzFMuCrSIWzL2z4FcLzbR42jnuQ/M/2g5HvaN22YYPdkGNJ3tJ7PTKvYtvoJbYXIw8E8kk/6IIrPw7SzO/TLnIyNG2Zzvjx0iGHduWyEnHm2XkdG8wFyJP/fDwILAQ0xkoa/Z4XIeeEL9+H7JM68Jnzbfd9FTRe7UR9n476JUw0Lm5oq1/betxtQGPbO+pDu6g/lVJsovoEzkoz0BgrSTQnJqI+8Yf6vllovPyhvk9H/RImGhc3tNWvbT3uNqCx7R31oV3Un4qiKIqiKIqiKIqiKIqiKIqiKIqiKIqiKIqiKIqiKIqiKIqiKIqiKIqiKIqiKIqiKIqiKIqiKIqiKIqiKIqiKAHwf+EpYFfXvLtVAAAAAElFTkSuQmCC\n",
|
||
"text/latex": [
|
||
"$$\\left [ - \\gamma {u}_{k - 1,i + 1} - \\gamma {u}_{k - 1,i - 1} + 2 \\gamma {u}_{k - 1,i} - \\gamma {u}_{k,i + 1} - \\gamma {u}_{k,i - 1} + 2 \\gamma {u}_{k,i} - 2 {u}_{k - 1,i} + 2 {u}_{k,i}, \\quad - H {u}_{k,0} + \\frac{- {u}_{k,0} + {u}_{k,1}}{h_{x}}, \\quad H {u}_{k,I} + \\frac{- {u}_{k,I - 1} + {u}_{k,I}}{h_{x}}\\right ]$$"
|
||
],
|
||
"text/plain": [
|
||
"⎡ \n",
|
||
"⎢-γ⋅u[k - 1, i + 1] - γ⋅u[k - 1, i - 1] + 2⋅γ⋅u[k - 1, i] - γ⋅u[k, i + 1] - γ⋅\n",
|
||
"⎣ \n",
|
||
"\n",
|
||
" -u[k, 0] +\n",
|
||
"u[k, i - 1] + 2⋅γ⋅u[k, i] - 2⋅u[k - 1, i] + 2⋅u[k, i], -H⋅u[k, 0] + ──────────\n",
|
||
" hₓ\n",
|
||
"\n",
|
||
" u[k, 1] -u[k, I - 1] + u[k, I]⎤\n",
|
||
"────────, H⋅u[k, I] + ──────────────────────⎥\n",
|
||
" hₓ ⎦"
|
||
]
|
||
},
|
||
"execution_count": 8,
|
||
"metadata": {},
|
||
"output_type": "execute_result"
|
||
}
|
||
],
|
||
"source": [
|
||
"eqs = [\n",
|
||
" (u[k, i] - u[k - 1, i]) / h_t - D / 2 * (\n",
|
||
" (u[k-1, i + 1] - 2 * u[k-1, i] + u[k-1, i - 1]) +\n",
|
||
" (u[k, i + 1] - 2 * u[k, i] + u[k, i - 1])\n",
|
||
" ) / (h_x**2),\n",
|
||
" (u[k, 1] - u[k, 0]) / h_x - H * u[k, 0],\n",
|
||
" (u[k, I] - u[k, I - 1]) / h_x + H * u[k, I],\n",
|
||
"]\n",
|
||
"eqs[0] = (eqs[0] * 2 * h_t).expand().subs({gamma_val:gamma})\n",
|
||
"eqs"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"metadata": {
|
||
"hidden": true
|
||
},
|
||
"source": [
|
||
"### $i=1,I-1$"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 123,
|
||
"metadata": {
|
||
"hidden": true
|
||
},
|
||
"outputs": [
|
||
{
|
||
"data": {
|
||
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAuEAAAAZCAYAAACfHvDPAAAABHNCSVQICAgIfAhkiAAACGxJREFUeJztnXuoZVUdxz8j2sydO+ATERMdzVcieCtHRoq4miZRminRHwltMnIgC4kofKRjSQYFaSCJhlxFQUUZcxJ8gE0zo5SvRFMnsbpGjmaiiTITvqY/fut0z913n7PX3nvtvdc++/uByzln7fX4rd/vu8/97dc6IIQQQgghhBBCCCFECW4EXgWm2zZkQvgEsAv4RtuGRIC0FRZpSwghhJgQ1gAfAN+toe99sWRhA/ACsBN4E9gKnAvsVsOYoSk7hw3Ay8CqBmyMFWlrPNKWEEII0WPuB/4DTNXQ9zrsrN124BbgSuAGN94u4A5gWQ3jhqTsHE5w2y9qxswokbbGI22J3pBgop1t1wxRIwmKcVskyPdtkVDe90diZyqvC2jPMCcDp7P0jN4BwD8wu8+uaexQVJnDc8CLGW27QoK0lUdCeR/1WVuVGJ70l4B3gT8Ch4yofznmzG8CF7r3Z2XUO9ht2xDMUlGVovEFxbhrKMbtId+3y9exM2231dT/g8BGLBkb5hXgWvd+dqi8jB7qpugchrkV0+SpqfI+aDg2bUF8+qpDW9ADfQ0n4a8CD2CXB76fUfcwV/4o8GvspnqAxzLqrnGvT4QxUwSgaHxBMe4ainF7yPftcgrwPvCHFsZ+172+N1RWRg9tkjWHYR5yr+lEqQ8ajk1b0C19ldUW9ENfi1gOvI0dXaX5LSbE493nvwP/HtHPldgRyhdCGygqXTIqEl9QjNsiQTFuiwT5vi0Syvl+GvsH/3Rge3zY3Y27Czgtta2oHtpi3BwG7Om2P5Iq74qGEyZLWxBeXwnhb8Wroi3ojr6C8gj2VOswZ2CT/ZX7vI/7fO+IPh5w2w9Mlf/UbYuZa4j78kZCtR3FJ76gGLdJwuTGWL435PulJJTz/ZGu3f2B7fHh527se0Zs99VDm+TNYcBO7PaCAV3ScMLkaQvC6ishfBJeVltQTl9dyD8gR8c3YRM7yH2eAv6GHY3s7cpOcXWuGNHH68C/MsrvBX5W3F5vPg3cDbyE2ZeU6GNvwi6XU8WmedfG92/Oo0+f+IJi3JRN8/QrxvK9Id+H8/2Jrs64e3aLjnezxxy+4+o+hyUMWfjqoSmby8xhwEssvqUgVg3P0w9tQXl9FbV5zsPmsnOApdqCcvqqO/+AwDnI7hkbt7nXY4B/YjfGH4qt9fiG2za4T+fxjPaHuQGyjl5mMNHUxSrgz26MsuO8kV8FMFHOA+trtOkqYK9U2QzwReyHA+ZT25706NMnvqAYF0Ex9kf7lxGz732Jxfc73euKMXX+Cvy3gH3bc7afD1wNPAt8BksIsvDVQxahbU7jO4cBUyz4GuLVcF+0BeX1Vcd33zBVtQXl9FV3/gEN5CBnYdn9BcBHMHE9zOI1Hm91dQ7JaH++2/aTVPkBrvwY93na9fMEsLqI9Z68TfEjlIMwG4/2qDtHfoKQpoxNaRLMxtmS7X3iC4pxWRTj0Wj/WiBm35ehTd8f6NptrTi+Lxe48Z4G9s+p66uHpikyB7BFHD7AEs4BXdJwwuRpC8LqK6Had9+AENqC4vpqOv+AADlI1rqMw0dWV2Nny7/lGg04Gnvi9cVU2+XAee59+onVGexI5y/AUdi9TO8Bn2TpUVdbHAfsAJ5v25Aa8YkvKMZdJtYYy/cLyPfheBm7/H5UA2P9APgFdmbwJGyVinH46qFJis4BzLfLWHxGtA8ajllbEJ++QmkLiuurC/kHpHSclYS/gD1V+2Xg89jN/X9K1XkH2AM4YqhsGvuFpGPd53SbGezI6EzsSO164ByWXoJok4GN6bUuJwmf+IJi3GVijbF8v4B8H45dwGZgP+DwGsf5Ifbg1+PYJfbXPNr46qEpyswBYK17/d1QWR80HLO2IC59hdQWFNdXF/IPSOk4657wd7BlYQ7HjmIuyahzH7ZG42bsCc9VmNOfwo4cV2IPB6QHPgJz4BnA78cYeQVwcc5ETgI25dQpynGMvvfpIhb/tOpybAf93lDZ54AtgW0KjU98oZ8xbsum0MQS4zTavxaIyfeToPs7sV/kOw1LTELzNeBHWMKzBXvoLM08Sx9g89VDE5SdA8BnXbvfDJXFpOE6iVVbEI++QmsLiuurC/kHeOp4I+Of+lyBXfrYjp1Wfwy7PLAXlt1vymizDVs0fgf2a0/j2A+7FDHub2VOH2Xu1XkeWDdi2z6Y0Ad/dwK/TJVN1WBTmoTq923lxRf6GeO2bEqTMBkxTqP9a4GYfD8Juv8QtlJC1prJIVhP/goSm0a09dFDE6yn3Bz2xM4m3pUqj0nDeSRMprYgnL4SyvtoPWG1BcX11XT+AeFzkP+zFTsy8THCh5WuvzXAVzHDPx6o71EUdc40ZuPavIqOOdp5cCwEoeMLkxnjMijG2Wj/qkaTvi9DDL4f/MT1x1q2I00demiSb2N+/VTFfmLX8Dhi1RZ0W18xaKsKteQgy4C3sKVlQrEWu0F+cCbrx9gaix8OOAbY5YoZ97cDuNS9P9ij7YkUE/IcfklCFZvqoI74wmTGuAmb6iDGGGv/qkaTvvclNt+vwB7i2tjS+FnUpYemmMLORN4RoK8YNexLjNqCbusrFm0VpfYcZPALUbdUMnMx67AF2wcsA27HbuAPuVPNkr/QfOLKVmfYuA1/5vBLEnxsapI64guTGeOQNjVJmzFO0P7Vdd/7Mktcvgf7IY3LsDNOMVCXHprio9h+uDpAXzFquAixaQu6ra+mtRWKWcZ/7yVka3hgZ66Ov8LSB6ImicuBZ8h+KLUPTHp8QTFuM8byvXwvFujD921IpOFiSF/xIQ3n8CjVF6MXcaMYt4d83x7yveg60rDoOtKwEEIIIYQQQgghhBBCCCGEEEIIIYQQQgghhBBC9IT/ATdDtjEGsW6yAAAAAElFTkSuQmCC\n",
|
||
"text/latex": [
|
||
"$$\\gamma {u}_{k - 1,i + 1} + \\gamma {u}_{k - 1,i - 1} + \\gamma {u}_{k,i + 1} + \\gamma {u}_{k,i - 1} + \\left(- 2 \\gamma - 2\\right) {u}_{k,i} + \\left(- 2 \\gamma + 2\\right) {u}_{k - 1,i}$$"
|
||
],
|
||
"text/plain": [
|
||
"γ⋅u[k - 1, i + 1] + γ⋅u[k - 1, i - 1] + γ⋅u[k, i + 1] + γ⋅u[k, i - 1] + (-2⋅γ \n",
|
||
"- 2)⋅u[k, i] + (-2⋅γ + 2)⋅u[k - 1, i]"
|
||
]
|
||
},
|
||
"execution_count": 123,
|
||
"metadata": {},
|
||
"output_type": "execute_result"
|
||
}
|
||
],
|
||
"source": [
|
||
"(-eqs[0]).collect([\n",
|
||
" u[k, i-1],\n",
|
||
" u[k, i],\n",
|
||
" u[k, i+1],\n",
|
||
" u[k-1, i-1],\n",
|
||
" u[k-1, i],\n",
|
||
" u[k-1, i+1],\n",
|
||
"])"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"metadata": {
|
||
"hidden": true
|
||
},
|
||
"source": [
|
||
"### $i=0$"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 124,
|
||
"metadata": {
|
||
"hidden": true
|
||
},
|
||
"outputs": [
|
||
{
|
||
"data": {
|
||
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAL0AAAAZCAYAAABkWi/IAAAABHNCSVQICAgIfAhkiAAABdxJREFUeJztm11sFUUUx39FrP0iKJJoJKGEWCmEhFsqWgSb+m1iVFQM8SPmBoko0QeNJmpM1WjUJ2M1JioG0fAAxgbUB6QvVUSlKl+KFkWwmEgF8QNBEBTw4T+bu507d+/edru9l+4vuZnuzJzZObszZ86c2UJCQgIAbwJ7geoI27wfOAHcGmGbxUYj0nHBUHckoTBmAMeBBxxlHeiltgTIv2Lq3GXlLzP59QPvYmzMBV4CPgb+Qv1flkdmJdAL1Axu1xKipAP4E6h0lP0G/EfwCrAJDY7pVv63wAFgRAR9jIvNSJcDQDfhBv0Fpt6jg9u1hKg4D1n51xxlE9HL3BIgXwX8C/wDnOrLrwaOAeui6WZBpMm/OuXiEqAOKDPyYQY9aILsorQm+LDBfinz0Qte4ah7vkm7AtprBEYCX6HB75Ey99oITEED5xfgIPApcGGhHY+JTmA7GuyFsBwYD1zhKHvEtHejo2y8KVtZ4P1KgaLR2x70lyOLvN5RN8ygbzLpl1a+5+rUmrLRwFvAR8BM4D1gVLgulwSfmNQ16BtNaj8j0H4KZBxONopG75G+v6uRRe4G/nbU9Qb9PLTsu/Astq2Yp3ATMJu+yrWj2Z9CG8aTgS9M2uwoawT2AT85yrxnvGkwOjXEFI3e/kE/DjgFRR5syshYa5f1ssll6dNkz+Zuk1aEaLdU2I/2NeOt/DHABGBNDjnv5dvP6Dk0aMI8+2KkqPT2uzdnmvQPR7065JKsRxPA9Rtj6h4GvvHJVgCTgZ3AakfbE026o/DuZ9GDfEP/7w1T1ukoWxrBPXPxOzDWyvMmv2uJB73gvcBuKz+FIkn9YRHwI5qEG4CLQ8j0kP2sgn75Nvdx692MXOafTf/S/kK/pT9sUpfFzTUT/XiKbUb7Ao9p5j4dAXL70YsZKC8Ap1t5KeB6dODWY5X1dyCFoZLMM/Xw3LwNjvoTgTOADxxlKbQHKpR5QBsa+OtMuhoFE1xuhscONEnCYg9Wm7j1rgG2GtlA+XPQrHCFFZ83ZXcGyD9k6rxo5d9D7lPKUShE2mnlzwWOoI2vRxt6GWcF9MFFmv6HLP20ED5kOQLpZa9ey00btVkScK8pe8bKP9vkTzHX1aadjchlCKILWGzlbQeezSMXNXHr7ecglqX3uze9wK/AJIdgGEvfYNJcm1jXLG9ArpFd1g58DTxmrh8EbgGuBvYE9KFYmIT0sleSehTK3WXlnwYsNH/bzziFVozvTLufowPCWWSvXH7K0bO3V9gO4KJ8CkRMnHrnxT/oTwBrkR96rlWnAThKX1/dJpffNh1Z7a0OGW9C2Ap7J5pp4GGgFbgGWalSwAvd2ivYUXRoV+fLqwaWAFPNtR3BSCEDMAedaSwGbifbdbIZiwITtpHYg6xonMSpd15GWtftwE3AVcAPJq8e+UibTOdd1CCFDgLbfPnlSKkt9D2s8ghaBTpQ6O9p4FoyYcA4mWN+kBkoM8lsgPehVcjmSrSvedfKX4Ni0mvRQUwNcBk6zOtFJ9o7LZkUerZLgOvQ2UapUdR6lyNL4D+AugNZ3tcD5GaRWSn8eF8dvppDrht911LmKLsUnRccI/s7nkJI03+f/gmCoxY9DpnRyBqtcpRVoL3JbuAQWhUXos33ceBDh8w29OwPATcU0Pdy5A7cbOW/TPwTJ069bbJ8ehfecXFDvoqDyDQU0ZmPVp9c8d1i5D70/GZH0FYVmvQzgNvQCyzEAHSR/R3V98S/kS2UgertJ9Sgr0Abjvf7eZOBUossQqu5noqsQcsQ9acQKlHf34movSZkrb0vXp9CsedxIeXnIZd0AToraUODwBVFKSYGqncNco9SaKVoNX/bh4V9aAYeJ9p/IgnDGOTy2O7QCuCzmPvSHyYjl2hCRO3dTebEGuQGvo32QFUmL41Wllz3XITcsCNGzvVpRLExUL1bcLujSwehrwlDwJMoqmYHJU52hqveCSiq1TLUnRgChqveCQkJCQkJCQkJCcOe/wHOT5SuMUTYXgAAAABJRU5ErkJggg==\n",
|
||
"text/latex": [
|
||
"$$\\left(H h_{x} + 1\\right) {u}_{k,0} - {u}_{k,1}$$"
|
||
],
|
||
"text/plain": [
|
||
"(H⋅hₓ + 1)⋅u[k, 0] - u[k, 1]"
|
||
]
|
||
},
|
||
"execution_count": 124,
|
||
"metadata": {},
|
||
"output_type": "execute_result"
|
||
}
|
||
],
|
||
"source": [
|
||
"(-eqs[1] * h_x).expand().collect([\n",
|
||
" u[k, 0],\n",
|
||
" u[k, 1],\n",
|
||
"])"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"metadata": {
|
||
"hidden": true
|
||
},
|
||
"source": [
|
||
"### $i=I$"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 125,
|
||
"metadata": {
|
||
"hidden": true
|
||
},
|
||
"outputs": [
|
||
{
|
||
"data": {
|
||
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAM4AAAAZCAYAAAB5JBFTAAAABHNCSVQICAgIfAhkiAAABb9JREFUeJztm39snVMYxz+dmq7tMmYJsVin1LZmyW5bo7KlaiEh/jDSRIbIzYgR/EFIECmyBX+I+JGIX6FkoRPLhD8wkfltjG0YJTI6yTrG0G02Zlv98T1v+vbcc9/7vnfvvX1v+36Sm9N7fr3Pc+4553nOc95CSkpKbDwP7ATqYuzzZmAIuDzGPpNGG9LxmtEWJKX8zAcOA7c4ytaiidEZ0P4JU+daK3+lyZ995CKWjS7gMeADYDeSf2WBNmuAHUB9aUVLSRprgb+ASY6yXcBBgi3RJjTBWq38b4E9wIQYZCwXm5Eue4A+wi2cM029O0srWkqSOB1Zm6ccZY1oQnwZ0L4W+A/4Bzjal18HHAI+jEfMSGQpbCXzcS7QBFSZ9mEWDmiRbaOyNomUCNg/7FI0SVY56p5h0k8D+msDqoGv0ALyyJhnbQSa0eT7BdgLfAycFVXwMrEO+AEtmCj0AjOA8x1ld5j+LnWUzTBlayI+L2mMeR3thXMesgzrHXXDLJx2k35u5XtuW4MpmwK8ALwHnA28BkwOJ3JF8JFJXQunzaT2GIHOl6ANppIZ8zpW+/6uQ5ahD/jbUddbOJchF8aFZznsAfMGsh1YyMhBW412pgw6hI8FNpi0w1HWBvwO/Owo88Z4UymEKiNjXkf/wpkOHIUiQjZVDFsN1y5qk8/iZMndafpMWhOi30phEJ3zZlj5U4GZwFt52nmTyh6jB9BkDDP2o8140HGEq3a8Sf901GtC7tV6tIhcn6mm7n7gG1/bGmAO8CPwhqPvRpNujS5+Dv3If/Z/njNl6xxlPTE8Mx9/ANOsPG8DcbkwoImzExiw8jMowlcM7wAvBpT3kzsuQZ9CwZEk6hiVDnR82I50ztoV/BZnv0ldO3++XcKPN2Cb0TnJY555ztqAdoPATwF9h+Vh4FgrLwNcjC51+62yYn+oMExieEw9PJf1C0f9RuA44E1HWQadCYuhFVgRUL4VWcew2BPeJok6RqUe2ILkKSjTSWh1uULGD5myqwPa32bqPGrlX0/+2/TJKPy9zsrvAv5FwQSPR9CPfEKADC6yFB+O9tNJ+HD0BKSXbUV7TR8NOS3gRlN2n5V/oslvNt/rTD8bkUsUxKnEo3sUxpqOe3FYHL+rtgP4DZjlaBzG4rSYNF9gwLUDtSA3zy5bDXwN3GW+3wosAS4Afg2QISnMQnrZFm02CtNvs/KPAZaZv+0xziDL9b3p9zN0Cb2AXAtq470CVM4I1njQccTCGQLeR375aVadFuAAI88uNvl821ZkPbY42niLylbau3nPArcD3cBF6E6lEvDC8rYlPYAuhpt8eXXAs8Bc892ONmXQJrIY3Xk9DVxJrhvoog1Zvd1hBY+B8aBjDkvQpL3Bl9dM4RVdj8419is1E9FAbnA1Yvj9NZeVAw3iQeDCQoIHkKV4U74YBRB6kF8+hH4kL+/BPO1eQnKfbOUvN33sAB5H/vN2FDQZQK86VVltelGgYRA4J6L8b5v25SRJOq6gcLCjs0D/TlfNZiJyhfyXnFeZBzwT0G4BwxbLj2dGn8zTrg8tNnsgARah+6RD5L73FoUsxS+cewge9H5Hmylop3zVUVaDzmoDwD5knZehgMZh4F1Hm+/Q2O8DLoko/y509iwnSdJxGnIdgz61BfoPtXBg+HWJlkIVS8g8tPssReedfHcCSeQmNH4LY+irFm0c84Er0I8YdhM5xcixKAY5SknSdQy9cGrQwe71EgoTRAParbrN97lop+ocJXmiMAnJ/kpM/bUjl897U305cnumh2jbhSaVHZ5PGknUsR6duzLICnabv+0L7Rw6gLuJ9x/ZwjAVuW+2a7cK+KTMshTDHOTezYypv+sYfrMC5NK+jKKQnouRRZPHfub9xHOpXGqSqGMnbte8pwTPShkl7kXRzupCFSuY8aBjSpnZQGW4sUfCeNAxJSUlJSUlJSUlJSVl1PgfCDirCadaf2EAAAAASUVORK5CYII=\n",
|
||
"text/latex": [
|
||
"$$\\left(H h_{x} + 1\\right) {u}_{k,I} - {u}_{k,I - 1}$$"
|
||
],
|
||
"text/plain": [
|
||
"(H⋅hₓ + 1)⋅u[k, I] - u[k, I - 1]"
|
||
]
|
||
},
|
||
"execution_count": 125,
|
||
"metadata": {},
|
||
"output_type": "execute_result"
|
||
}
|
||
],
|
||
"source": [
|
||
"(eqs[2] * h_x).expand().collect([\n",
|
||
" u[k, I-1],\n",
|
||
" u[k, I],\n",
|
||
"])"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"metadata": {},
|
||
"source": [
|
||
"## простейшая неявная повышенного порядка аппроксимации"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 9,
|
||
"metadata": {},
|
||
"outputs": [
|
||
{
|
||
"data": {
|
||
"image/png": "iVBORw0KGgoAAAANSUhEUgAABAEAAAApCAYAAABEMzBQAAAABHNCSVQICAgIfAhkiAAAErlJREFUeJztnXvQH1V5xz8BibmASGAqQkq4yM1yScBIwiVdA1bryGg1yCiV/ipIqZeOVpxKRrkURhhtlchgqcX62toKxQxFW7kIilKICCHhUqKB4GuLEbyjFhUC6R/P2Xn33Xfv9/3t9zOz87579pzdc57f2e+ePXvOc0AIIYQQQgghhBCDZTnwGHAv8E23LWspL7sBTwAHtHR9Md5cC7yv7UwIIYQQQgghRMUcx9T7/AbsHX9pXGQP2A4cX1NmbnXn3w5sw17y/wM4KSLuR4HP1JQP0Q3OBe4GfgH8CPgScFhD1z4c+Cmwa0PXE0IIIYQQQoim8Qi84+/QQgaOAs4HXgy8BDgFexH7CvCWQLx5wJnAp5vOYI1MABe0nIeu4QGfBI4FVmIdQ7cACxq49gPAo8AfN3AtIYQQQgghhOgkHvWNBDjAnduLOHY98L3A/iqsc2BWIOwk4GlgdiBsoTvnIVVmtCYmyN4J0PeyFmVn4Fng5EBYnbY4D/ivkucQQgghhBBCiK7i0eJIgKPdxe+NOHYjsA9TX4BPANa7+D5LgIewF8Jg2FPA5qoz2zJDKmuQXbB6+bNAWJ22+BbwcmBuyfMIIYQQQgghROdpoxNgCzb/O4z/gveM+7sI2BqKsxjYGApbAjwIPOf2XwdcljE/KynuGO467EX1CwXTp5FW1nEpZ5g1WLnXBcLqtMVWYCdgryKZFUIIIYQQQog+0XQnwFHY1/0oDsZeyH7p9ucCvwnFWQLcFwoLvyAeEREnjq8CfxtzbOS2ONYAp6ecfzXwq8B2WkTYCTFp08paVTnTyFLOi5ly9hi3eRmu9TFsiMobsSkBPnXa4tfur0YCCCGEEEIIIcaernQC7IQ5CFwbCPsxtkSgz1zgIGZ+EV7K9BdA/4XwBZifgbMS8vNFzEN8EW5jqsMijiuxl1V/+2JE2D0R6bKUtUvlvAw4NGX7Vso5Pg68GftS/2ggvG5b+NNPfpSSvyb4LPBDYH7bGRGdxZ9SdWbB9O916d+SFlGInEi/RFbK6lgRpH3DQnokqmJQeuVRj2PA/dx5V4bCZwFXYE4A9w6En4MN9/Y5zKVfGAhb4cKODYQ9iL103gG8MiVPm7EOiChGJI8EALNVnmHyE2RzDJilrFWVMwse9U4HWAM8jpUnTN22OANbM7NtlmJTG/6y7YyIznMd8APMiWYSL8JG1FweCPsc4+9ctCw3kz566UoXJ6nDcUhIv0ResuhYWMN2xxri1wGPYCP5nsSc+55B8oetcdI+aVQy0iNRNUXbXUU1qy698kh5x0+NUJBVTBVoT+BA7Ov/17Gv/itC8Q/HDLm7298Tu6lf7/ZfBnzHhfk/ylxs/voG4PdS8rMLcH8obDb2xXkj8D9u8/dnMxOPejoB0spatpxgZXowYouaG+9RXyfAFZiPiJVYuf3N/03rtsUE3ViG8mbg52hagkjn5ZiWrk6Jd5aL94pA2EPYyJ42loftCz/BlipN+oK0AbPtUY3kqPtIv0ResuhYWMPOdvtbgX8BLgH+Eat727F2yqyZpwHGS/ukUclIj0TVFG13FdWsuvTKo6VOgEuYmh++DROxdcCFWM9JFOuAdwb2P+DS/S9mzAuBhwPHlwJfxoaGH5ySn2XA5xOOj2hvJAAkl7XKcmbBo75OgDgfAhcE4tRlizlYj9yyMgXA6klWvwdRHIR1anyqZD5ENkaU+726wCZsSdWkB8QNWAfrjm5/PtaxqiUx49kfqxtJPkbmYQ5sf0O5EVZdYkTxe0L6VS8j+q9XcaTpWFjDVmLLB4fj74l9tNmO+RQKM07aNwSNGiE9Et2jSLuriGbVqVceLXUCFOHV2FffHdMiOs4AzmXK90Bw2MatTJ9u8HaSe3RGVN8JUBVVljMLHu2UMwtlbPFOrMe4LCPKNdIudelPrCAvUZzrzv+GiGP7uGPX1XTtLjKi+O/VFVue7671qpjjLwB+C3wmEHacS/MJ4KXYcLPHMQeldwLH1JXZHvEmzEZJDcgTXJywn5Ou1I0ijCh+T9StX9Bv25ZlRLnnS5dtl6RjURqWxGp3rssjjo2T9g1Bo0bUp0d9sUEfGJoti7S7kojTrDr1yiPwjt/1YVE3YsPFF6ZFdByBDWm/F/gkNuQCbKjFSzC/Az6HM93nQF5uAa4FXoPNKV9e4lx5GUo5s1DGFs8A724gj2mchPX6fbOm8x/t/kY5oVzq/t5b07XHja7Y8g73N84HxmuxKUzBB7A/LHQRlv9dgX/CpmQtx5xm7lJ5TvvFy9zfuxLi+COHwnWgK3WjaerWLxiubaugy7ZL0rEoDUvCX156W8SxcdI+aVQyaXo0BBs0xdBsWaTdlUScZrWqVx7dGQlQFYdiy8+NO0MpZxaatMWI4r3W8zEBeKDC/IT5LvGrH/jTdF5b4/W7xojiv1dXbLkr0V96fK7Feo3nBMImXJonmDlPdK07Frdk6VD4KmaHm7He96jtYRfnbaG0XakbRRhR7J5oQr+g37Yty4hyIwG6bLskHYvSsDieh9XBuK90E4yP9g1Bo0bUp0d9sUEfGJoti7S74kjSrAnq0yuPwDv+8wqepG9sYhheQodSziz0xRZ7Y9NdflDT+RcA+wI3xRz3vyqMU29tXXTJlk9i8z33iTg2B5tKdYOL4+M/TEbMzOOmQNqhMospG6WtMgLTv350qW40Sd36BcO1bRV03XZxOhanYXFciq0k9GWiyzou2ieNSiZNj4Zgg6YYoi2LtLviSNKsVvXKI34kwCTxjtyits/VlUkhWmKSfPfARMr5lrt411R4zeB9d5ILuzjm3D/FehvDXAp8JSXvbXMF6cOuJqnu92rDlkll/D7RQ19PxvIZXF92Djb0bEvMuf7Vpdm/WDZbY5LqnkkHuTjrEuLs5uI8xXRfNX26zyap7p6oW7+gX7YtyyTVPl/6YLsoHYvSsDj+wsXdhL2YhBkn7RtHjZqkOT3qqg2aYpLqnpdDtWWedlccSZpVt155Ln2hkQBbyNbD4bM15/mF6DqXAS8MhS0GXgd8FhPZIBtTzvdr9zepV6/MfefP2VofEW9/rMFwY8SxxaTnvQwrgHOw/O0F/CnpDdowH2RqTlUcVf5ebdgyqYxzmao/Qd4APA38ZyDsSEzv4xxhHoX1cn83FP5RzO9GnCOctqnymZTly4XfQ78Rm3fq00bdeAfwfuDFwH8D7wFuz5Cuynuibv2C7mpYEkX1rernS9O2K1LuKB2L0rAo3gWswZbTOpHpPoB8impfF+mbRt2KvQgmvRg1qUddtUEeyrSdqnxejoMtIb8987S7okjTrNb1yiPQSyCESGWE3TNegbR7ubR1LV10tTv/oohj73LHPhxx7HGqFdowr3HXXYV9sRjVeK0wI4r9Xl2y5Q7YMkjh3uIdseVpbgiF/7nL35kR59rFnetrEce+DlxUKqf94WOYjc5IiPN+F+cTofCm68apWOfQ2zH/J5djcxGjhilmYUSxe6Ju/YJu3XdZqVLfRhR/vjRtu7zljtKxOA0L8x4s/w8Av5MQr6j2dZE+aRTAz4D3FUg3oh496rINJsi2hHibbacgXbZlHvLYM2+7K0wWzapbrzxS3vFTIxQkzxAUbdra3vIwcmm8nOnA5vj9kHjnKmXZiPVOhnk+U05JVoWO7enCX+r252OCfy82B6xqfkX+B9lCLI+HFLjeiGK/V9O2TCrjoe7Y2lD4K1z4WaHwq1z4kohzrXDH/iYQNhsra/CeeCglv33nG8TbyMcfind6KLzpunEX8A+hsIcxZ0xFGFHsnqhbv6AfGpZEEX0LMqL486VN22Upd5SOxWlYkL9ycTYAe6RcI6/2rcKW+VoUCFuDNfxflHKtuumTRh1A8Xo7Kpg2TY+6bIMJsnUCBCmrLWXosi2LkmbPvO2uIFk1K69eQT7N8tw5Gl8icJa22rZ9sC92m7Cb700dyFPft6bYjj3Y98CWNKyap4GdgAMDYfOxpRQPc/sbQmkWY8OdvgMcjHlC3YatXTpZQx6LcCTWa7u5wWs2bcukMvpLQH0tFP5HWE/x9aHwo7CHRNRyof6wvuAQ021MLQd6DDbk/LjA8Qms7o7iMt8zdsAeuk9jQ+vj8IfahpdEarJuzMZ+s/BwwZuBYxPS1UHd+gXjq2FN0HXbRelYnIb5fAibV7weG07745Rr5NW+tVg76oNu/xzgzZjTr+Ac5wma1cA+aRSYbbfTrGO4ND0agg2aYoi2zNvu8smjWXn1CrJr1gyGsjrAuLMNG2ayEetpW495nPy/NjMlMrMWeCM27/qRis99E7Ze6zcwB3M7YyJ0P+ZBdx7waCjNYkxQXg98ChsOflnF+SqLn8fnGrxm07ZMKuMfYPM9gw+dWe4665gu/LOxh/J9RPsXiJrb9xz24v9L4G5mjo7xO5CjHOT0kUOw33MD0V83cMcPxL4WfDt0rMm6sQc2/DD8cH8Cc9bUNHXqF4yvhjVB120X1rE4DfP5E+CvXZrbMQdbYSaZmtNbRPu2A6uxub1b3P8nYiNtgjStgX3SKDDbbgF+kTF+VSTp0VBs0ARDtGWedpdPHs0qoleQXbMy4REYKiB6yX3A77adiQExotwwpdmYeNxVUX6CzMGGBW3FvirfA/wZ5ojnOeC2iDRXY85KngR+P+X8F5M+tcJLOUeRIW3/BlyZM43PiGK/V922DBNXxl2x3vR/D4UvxcoVXhrT70H/+5jrbMJe9sMjYD5EvKO5DdjDeLeY433jdMxGVyXEOY6pL01hmqwb/rzXFaHw87AvLEUYUVzD6tQvaP6+C9KWvgUZZbxOFG3aLq3cUToWp2E+F5D+e9wWiF9U+wDuxF7w/zAmbdMa2CeNAvMAf3XOND4j6tGjLtlgNXaP+NszWOdOMCxtLfg2pwN0yZZV6DQk2zNvu8vnggx5u83FLaNXkK5Z0KJPAFGcS4FbMsY9muhhJKLbnEv6XL+m+DbW0HgKG+aUxB7YF4qkbV7KOYo8yDYDZ+dM0wZ5bBkmrozvJlqjP+zC98t5nTjWMtO5FNhD/lngIxVdZ6gUrRuzsQf9KaHwK7BpYW3QJf2CcvddkLb0rU2qsl1auaN0rGoNK8pKbCTls0wNrw8yFA0sUxd+gjkpbIMq9aguGyzApiz4m/+8DYbNTTn/kLQlyZZV6DQk27OpdldR0jTLx0OdAL3jRmY6gohiATZXrOl5oaI8c4DvAV9qOR/zMBFZCpyGiWKSoFRB3gfZfCyPy9IitkwZW8aVcS7W8/6FiDSbqHZJtEeI9lB7MrbM0J4VXmtolL3P7sKGUwbZTHHHgGXpin5BOxqWRJ8a6lXaLqnccTpWtYYV4UjsS+XbsBezmyLiDEEDy9SF/bD3iJX1ZC2VqvSoSRtM0C/HgHnpQ32Ks2eT7a4iZNEsHw91AvSOx4G3psR5PjYELC2e6C4rgPOxF8C2WIZ9ZfR7oC8Cvg/sXfF1dsbmhy3GeoXPc/9nWeJsOfYwydKz2yZlbBlXxkOxhsK+leQwmUms83EvZq7jLMpR9j47FRs6eiZWJ9ZgDZhFSYlqpgv6Bc1pWBJl9K1Nytoua7mb1LE8LMIa++e5/cOwoc1eWxlqkTJ1YRX2HtHmc6MKPWrSBhNk6wQYorbUWZ+y2LOregX5NctDnQC9wl9S41SmnP1twZak8JkFfJ78vYhChDkb69n0mYXNTV9PtS/dHtFzoyYCcUYubN+IPIadHnWRLLYc0d0yngY8hj1Q/q7lvIwbZeqGzzuwjprfunRhHwFDpQrblsUjXd+6SFnbefSz3GAjKTcxcz7uNZjTr6FRpi5cwsy11PtIkzaYIFsb3qOf91hX65NHP+0JxTTLQ50AveLVTDl6eSU2T+h6pi+9cTzWUN8Y2A5vNptCVM6F2PSWcV7FZAhlFMVQ3agP2bY4sp3wUV2QDapEtqwfD3UC9IoPAD9n+tyzt2Jf6IQYZ+5m/IdhDqGMohiqG/Uh2xZHthM+qguyQZXIlvXjoU6AXnE18M+hsAuJXn5DCCGEEEIIIYQI4hF4x9+h1ayILCxm5tyOJbTvjVIIIYQQQgghRM9QJ0C3mQccyPT5/6BOACGEEEIIIYQQBVAnQLc5wv29PxC2O7CQ6Z0AdwDHuP8/Dby3/qwJIYQQQgghhOgb8sDYbRYDD2PLAvosAZ4BHgqEXYQ5ELwdWyXg401lUAghhBBCCCFEv/GQY8A+sh64Fdip7YwIIYQQQgghhOgMHnIMOHYsBRYAT2KjBIQQQgghhBBCiBmoE6D/7A1cBawE9gUOazU3QgghhBBCCCF6xXLgMeA+4B63LW81RyKOucCdwAq3fwpwTXvZEUIIIYQQQgjRAY5n6n3+fuwdf2mrORJCCCGEEEIIIYQQQgghhBBCCCGEEEIIIYQQQgghRFH+H+ez/E7drbKpAAAAAElFTkSuQmCC\n",
|
||
"text/latex": [
|
||
"$$\\left [ - \\frac{D \\left({u}_{k,i + 1} + {u}_{k,i - 1} - 2 {u}_{k,i}\\right)}{h_{x}^{2}} + \\frac{- {u}_{k - 1,i} + {u}_{k,i}}{h_{t}}, \\quad - H {u}_{k,0} + \\frac{- {u}_{k,-1} + {u}_{k,1}}{2 h_{x}}, \\quad H {u}_{k,I} + \\frac{{u}_{k,I + 1} - {u}_{k,I - 1}}{2 h_{x}}\\right ]$$"
|
||
],
|
||
"text/plain": [
|
||
"⎡ D⋅(u[k, i + 1] + u[k, i - 1] - 2⋅u[k, i]) -u[k - 1, i] + u[k, i] \n",
|
||
"⎢- ───────────────────────────────────────── + ──────────────────────, -H⋅u[k,\n",
|
||
"⎢ 2 hₜ \n",
|
||
"⎣ hₓ \n",
|
||
"\n",
|
||
" -u[k, -1] + u[k, 1] u[k, I + 1] - u[k, I - 1]⎤\n",
|
||
" 0] + ───────────────────, H⋅u[k, I] + ─────────────────────────⎥\n",
|
||
" 2⋅hₓ 2⋅hₓ ⎥\n",
|
||
" ⎦"
|
||
]
|
||
},
|
||
"execution_count": 9,
|
||
"metadata": {},
|
||
"output_type": "execute_result"
|
||
}
|
||
],
|
||
"source": [
|
||
"eqs = [\n",
|
||
" (u[k, i] - u[k - 1, i]) / h_t - D * (u[k, i + 1] - 2 * u[k, i] + u[k, i - 1]) / (h_x**2),\n",
|
||
" (u[k, 1] - u[k, -1]) / (2 * h_x) - H * u[k, 0],\n",
|
||
" (u[k, I+1] - u[k, I - 1]) / (2 * h_x) + H * u[k, I],\n",
|
||
"]\n",
|
||
"# eqs[0] = (eqs[0] * h_t).expand().subs({gamma_val:gamma})\n",
|
||
"eqs"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"metadata": {},
|
||
"source": [
|
||
"### $i=1,I-1$"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 10,
|
||
"metadata": {},
|
||
"outputs": [
|
||
{
|
||
"data": {
|
||
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAAqCAYAAACjkKBgAAAABHNCSVQICAgIfAhkiAAACgNJREFUeJztnXuwVVUdxz9X4novD1F0JgQC0QQx0HuhS5B6uxJO1OikpjVmNqexzMb8Q+sPdSaD0TFm7CE12NO6kiVFRGAFPUytyJKnRqAUiRYPs0SKtNSwP75rz9l33/06j733Off8PjN7zj2/vfZav7Xu2Wv/1m/91tpgGIZhtDTHAc8ApxStiDFkWQl8rGglDKOZuR941R2voE77h8CCQLrbgW/kq5qRMzcCG4F/As8C9wEzcix/JvAcMCbHMg1jSHEQuBkYB5wE9ALLUQf/XpdmhEt3dgH6ZUU/sKhgHRqNnwAfQJ34TGA1cAAYm6MOm4BrcizPMIYMp6COuy/k3BrgKff3JchyagukWQC8BLT7ZBNdnqfVU9EM6Cd9h97M9ayFUcD/gAt8sqzb4mbg13XIxzBagqN8f89GN+KWkHTrgUnIOjsH2OzS+ukGdqAb3C97AdhVJ30bgVapZ5DR6Pdy0CfLui0eAeYAnXXIyzCGPMEOfTfymQbxbtiXgcnAvpA0XcC2gKwb2A4ccd/fCdyRQq/5VD8hthp1Ot+r8vokWqWeQZaiej/sk2XdFvuA4cD4SpU1jFbE36HPQpZ3GNPQzfUvZC39JyRNN/BoQBa84c8ISRPGL4DPRJwruSOKpcD7E/K/CTjsOy4PkZ0TcW1e9UwiTT0BbqU80R119CXk8Vk0Z/Iu5HbxyLotXnSfZqEbRgrSdOjDgUuBVe7731HYop9OYCqDrbUeBt7M3s19DPLLXxWh11o0EVcND6IHTxxfQh2Pd6wNkW0Kua7Z6gmyjqcnHI/EXP854DJkQf/ZJ8+jLbwJ2Gdj9MuDu4G/ASML1sNoXDyX9QdryOM6BgagVM0Ul9H8gLwNWIYmQSc42cfRkNrPDHf9RJ+s18ne7JNtRx3IBuC8GH12oQdJGCXiLXSQxVmJK6KfdJOiedYzDX1k63JZiiJbpoecy6MtrgT+WpnKdacHuY+uL1gPo/FZDexHAQRxvBaNdL8QkN9DnQIKLvFlNA44FVnlDyGLvNeXdqZT5nifbBz60V/ovr8ReMLJvMp1Ip/vVuANMbqMBh4LyNqRJbgNeNod3vd2BtNHNh161vUE1Wl7yBHmR+4juw59GZpPmY/q7R1ePfNoi37grmorUCd+CjyPuX2MZOagfvSmhHRXuXTnBuQ70Kj7qEFXVMinGLig6B9o8msxepoEeZjB8cE3uOv+AnzLXftH3/ke4Mdo+D0tRpe5wL0x50sUZ6FDfvVMQx/ZdehR/vZFvjRZtkUHcMidK4qp6AH1lQJ1aCVKpJvTaWR2ohDvuE55HTKUh/lkI5GhXEiY7kJkjQ1LSujjSrT60PPVe1bc/ZTdOQAfIv4JV6L+HXo9qVc909BHcfVMQy1tcQ2yjotkCepg3pphGTe6Mi4OOTfJnVudYfmNRInaOvRGaMtPunLeFnH+GOC/DF5pf5a77vPA6cj9cgAFaPwGeFOawqs17dejIfnEpIQ+zkCugy3AncDXkY/+9chH7zGTwT76Svg52gfkHcj/Oq+GvKqhVeqZhlra4mXg2hx0jGMBspp+m2EZs91n2CR8j/sMWxtiDKYR2nKD+4yaLzofuYmDD5ZZ7nMy0n8MWqX/ELq31yLXZEMzHYXEDXVapZ5paJa2GIncj7/PuJwniY7i8Vyh52esQ6NQojYLvRHacowrJypybCWyujsC8n533TOUO3ePVe5cVCi1YRgJTEU3UZZun7GujPUR53/mzvsnxJc4eaOzjMrdGyWq79CLaMuoOr6I3CVBOtCk58qQc485/d4ecs5bSxIXJWYYQ4o9JC+g8h/3JOQ3z6X7ToZlLnDyWyPyfw5ZbH7Wo91Os6QXDfH3Iv1KVeRxHPHhe3uorO36E8oroi2j6rgXje6CXEB4nHkHcjHujijn2+66k32y29EGegN4Tby+htE07CZ8BXMUYdtX+PFWqQaHxvUs0/P5hi3oOxl1GEGLswv5VrNkFJrTWF5DWQcTzt8BHBuQdaGtIu5GHb6f4AK2IEW0ZVQdOyn/fvxcjLZR+VFAfibqi6NGg7NQxNeTPtkc4JepNTWMFmc8soqyDCNb4cqYHHLuo+7cbT7ZOCc73X0f6fLYgra7zoLDVG6hV7vjZonqXS55t2VUHY9Coa5Ba3sYClVcF5LXR4heZTra5feA+96OHgr+0csOL7FnoQd3TjSMRiW4bXNW7EcTbHHx87VyGhpqPxWQHw182P3tj8roQpbfE06v76NFW2cRbhEWxZnkv/to3m0ZVcdp6DcaHFH0osWYYT73uNFFt8vPO/cKcgduQqGMT6MwSKActthmRybHJBR2tBNFS7y7AXRq9iMvXkVD2hNQmGUWvIS2OzjVJxuJwju9t0Nt9Z3rQr+jC1Fs8leB99FYnTmU9TySlLCO5N2WUXX0FsI9EJBf5NKuCclrFuqUw8KYvc7eexgdAU5Ek6sb0eRrknvLqBMnon88aIi3F9vgqZm4DHXsWb016RaX/34Up78c/UbWIX/78wx8iK1Ak3uHgLck5F2PXTahOpfLd9Fmd5VSonqXS5ZtGUZUHe9FVvTrfLI2ZEmHue88F8rGiHK8/V38I8VPAL+qUF8jJUvQIp80PMrAf7TR2LSjyIjfZZR/B9oAbR8avm9C7oFjkSX2YCD948DXXNqLEvI+Abkh4o4RKXSspkPfBVxd4TVQW4eeZVuGEVbHMcjC/0FA3oPqFbbBm7dL45cjytmJrHH/w2gVWlFqZMB64NMp0s2mtpWhRjF4y8m7C9ZjBFq12oP27z/M4AUoWVBph+7tSVLkHjxJ1NqWUXW8Fv1Wgu9bvs3Jp1SjbAh/orZteo0YDgBXJKQZC/yBgdvKGs1BB5pou69gPeaioby36+MtyKUwIfKK6hlF+b0AL6B3u3ahOaEk5qHOLo31XxS1tmVYHTvR6CBsb6WdJIddVsIeZESOZ3Dop1EDXujTe9Cugv9G4Urn+tIcjSbXkjp9o3HpRZsuFTn/cTXqGDzakB93M/XvPPtIXtxTcrKTQvR8vM761Js0bVkivH7e9cE6Tkc7kYalrzeXo/2bjgBfzKG8lmEh5WiI81A0xBrKs+ltaJJkURHKGUaGLEajzqG6OHGo188I4QY0ez7OJ7uC8lt2zkZP0W2+o9pXzhlGI7GR5t63PImhXj8jhBXANwOyxQyeUTcMw8icml911OJ0obc3+emmvhMghmEYqbAOvXpGoFVpWwNy69ANwzCajLkodMkf+XA8miT1VoduoPzqqLuA63LTzjAMw0hNWOjSArSMt919X4g247ke7RVhGIZhNDGb0UuRhxetiGEYQxvzoWdLD1oleght7WkYhmE0IRPQZlxT0NaXM+KTG4ZhGI1IJ9pjudd9v5T4d1MahmEYhmEYhmEYhmEYhmEYhmEYhmEYhmEYcfwfPEkFnxCCClgAAAAASUVORK5CYII=\n",
|
||
"text/latex": [
|
||
"$$\\frac{D \\left({u}_{k,i + 1} + {u}_{k,i - 1} - 2 {u}_{k,i}\\right)}{h_{x}^{2}} - \\frac{- {u}_{k - 1,i} + {u}_{k,i}}{h_{t}}$$"
|
||
],
|
||
"text/plain": [
|
||
"D⋅(u[k, i + 1] + u[k, i - 1] - 2⋅u[k, i]) -u[k - 1, i] + u[k, i]\n",
|
||
"───────────────────────────────────────── - ──────────────────────\n",
|
||
" 2 hₜ \n",
|
||
" hₓ "
|
||
]
|
||
},
|
||
"execution_count": 10,
|
||
"metadata": {},
|
||
"output_type": "execute_result"
|
||
}
|
||
],
|
||
"source": [
|
||
"(-eqs[0]).collect([\n",
|
||
" u[k, i-1],\n",
|
||
" u[k, i],\n",
|
||
" u[k, i+1],\n",
|
||
" u[k-1, i-1],\n",
|
||
" u[k-1, i],\n",
|
||
" u[k-1, i+1],\n",
|
||
"])"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"metadata": {},
|
||
"source": [
|
||
"### $i=0$"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 15,
|
||
"metadata": {},
|
||
"outputs": [
|
||
{
|
||
"data": {
|
||
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAbcAAAAqCAYAAAAwGBXCAAAABHNCSVQICAgIfAhkiAAADehJREFUeJztnXuwVVUdxz/XqyBwwQRKUhMhMzKVi4qg6PWCNEMP03wxZdr1MYbolDqZj4qkSGlqVHQsy5HMzBRTRysfjKRg6qDhvahBPtBrFGKiZVAq4r398V17zjr77L3PPueesx/nrM/MnXPP2uu5z3r+1m/9FjgcDofD4UiUhcCDaWfC4WgQXHtyNA3b1TCui4Engf8ArwO/A/YN8LcM6Dd/24DXgN8DMwP8tgM9NcxjI7J92hlw1AXXnsrj6r4jER4ATkUNcD/gLmAjMNLn71/APGAMsCfQAdyEGueXfH43AifXLcf551RgatqZcNQF157KswAYmnYmHM1HG/A+cJTl9lHU6DoD/N8NvGJ9H2P8zgbuBf4LrAOmR6Q5E9gKDLLcdjfxTKgo99nnFOD8tDPhSIw02lPW2QtYQm0lUA5HWT6MGtNhltuJQB8wIsD/Wca/NzOdZb6vAD6FKvLdQHdEmhdQKnY5CjXkRmoAB6OZfUvaGXEkRhrtKQ/MAS5LOxOO5mIJajitltsPgRdC/J+OGt9w8/0i4N9oxulxMvB36/t44PPW918DN/rinQestL4fDVxVNvfiLiT2+W1M/0kwCHgWmJh2RhyJkkR7yiPbob3JaWlnxJEt6rWauQLNMI9DohSPA4BVIWE+DmwANpvv7WgTfaPlZy/gRev7p4F9rO+TgNW+eP2b6PsH+AljERL/ZYnzgOeIXwZH/kmiPfkninmhD/g+cB1OwcRRZ64EXiV4j+sNJDr0swPwMnC15fZXYK7P3z0UVl1HAG+imWs3MAxpi/n3ENb74rkddQojkFjmzMjSaD8jKyu3YcAm8r1P4uc8gpUf8ppOrUmqPZ2NVnd5pAVYA5yWdkYcjcsiNDP8RMCzcahzmeFzbwGuRQPVbsZtKJqhHuLzux7osr4/jDTEQFpl/UiBxKPDuB1quT1r8vco2nsoRyfZGdzOpXjl6mcUcAYSp74IvA28BfwJiamiVupLKaiU96MZ8WY0ebgViXPD9vi8sJ0R8V9n/PgnEzdTW4WfXVDduabO6SRBUu0paKKYBsej3+0RdASiH/1ucbgQ5b+1nEeHo1KuRRVyBpLre39t5vnxFDqXMcDHgBOA5Wg10mHFNRU1RruRjTLh2y23l6z/x6AO+Rjz/SAkvuuz8jAE7aF1A5+MWa5OsjO4PQ38OOL5HPSONqD9x8uBxWivpR+VI2yAegO9q/nApebzSuA+4B0TfhmFPRx/2G1Ed4rdJo4DfO5r0CBaKxH5mSYd/+q21unUm6Tb08MUJopp0YPytBlYS2WD2x7G/6z6ZM3RzPSH/F1qnl9uuW1DHeLjqBPdxRfXHCRGsfGr+e+OViQ2F5l416POfT7FG+6TkRr0arQnEYdOajO4dVF+dRPF/gTP1G1mIO1Qfwc+BvibCX9cQLjx5pn/ndvhHzJ+loSEjdoDHAq8hwbJHSz3YajT9f+OA+E+1LnbM/h6pBOHLqr/zZNuTy8xcLoYWB2fjgbpFhNHJYMbqIy3VJm2w5EZpqH9s0o4HVl+8Dbi26xnyyiIcWw6ycbgdjHFq9BKucSk7xfXgVTK+9GEIIwPoiMV/Ugk5g/784iwhxs/T/jcpxn3q5Fi0M1IDLcFeAyYEhFnECOAd4Ff1DmduHQxsN88KYImitXQRe3K20nlg9tP0Go3TLHkYhPnsQHPvJXfXRWklxeastx5EdEE8SxaNTxDscZkFPubcE+hhrDYuLcgzbE3ff4fRAPoZ5DKtH/PIkk60Cp0S5Xh3zOf2wKeHWQ+/xwR/nW0MoBiqyhe2JWE4/n3x++JKMeaZzsh6xrL0bu+h2AxaBifQysRf0OtdTqNxliktJJ3etDv6Bd9exxoPoPq+WTz+VStM5UBmrLceVadfYvCjxaXr1v/32D+QPsWdyAFDJsg+3xpMQVttFfD9hSONNwf8NwboMLUyj3eMJ8fCAg7m3AtTm9l5G9c3u83Fam62w3sDjTTbCd+ub+AVpdL65xOo2FPFGej/ck88oz5PJRSKQGoHmxCIno/Xj3O+6H2IJqy3Hke3GrJWrJtymo0sDPVH7hdiLRJ70WWTWxa0Ey3n/IVfCfz6Q1yXliIp3katnLronTmuNZ87hgjXs/fLAoKMPVKpxGpZqKYRbzOe++AZyORwoy//nt4nbxdPxai9xKnbmeVZi13Cb2Eb2QH/VUiD28meqnsPd5YJr4pxt+3q8jL10zYtZQa3QV1BFHKJDYvGb+edRQv7OOhITQo9wP/o1jJY0ckKl0XEu4WE268z/1HBDfUowg+x1ZtOpXSS21/86zTS33L20nlfUwrUhwKqh8zTXwLQsK+iW5UsLkf1bdqmYvOG76DpCKHxwjTS2374DTKHYcOtB3wD5S/ririiHy//pXbOkpnvVFsqCJDzcBVFIvuQGKvo4FfogpsU+4aktHmc3Okr1LOQWel1gBHUrqnCPH220AD2Tj0m3uakUGzPj/eqqmHYusaE1H984sQ7XBvocprczCyj+jnWKT99wefe7XpVEqtf/Osk8Xyvo/ayOiAZ97KNEj0Ph5Nwvwi+3a0N1sNs1Hbm4uUdeYiqcI+BIsHPWrdBydd7ri0IXH4TVWmVe37ddSBLqrXJPM0Er9aQZhzTZhngA9F+LvC+DuvTHw3GH+XBIQ9PSLcBRQ0FW08w75nBIQZjjRDH7LcBqHBy561entDrWhP4b6AuCpN53ikcTnWcluEOh2/in05usiHtmSt6CJdbUmQFuxzAe63mvjGBjw7xzyzDTB7tyh4ymrDTBxPEe884Ergep/bC+gIR5IkXe5q2ELlK7ey77dW2pKVLKMb5S9JBpvPrTH9X4gOYPcgJY9/RviNo0zyZWTaqJfiowRxVm6TzGeYMklQupPQfp79bBsFbdUpyEq+Zyy3Ax1KDlJnrjSdO9CEwBMBfwP4ItrP84tv6kXadTutel4LtlJoLzYTkHj6FZ/7YAqTRrsetyMFs+fQmdgnUB2cRumq1M8gVO/80oKlFFtLSoIky50Usd5vrQa3lhT/9kAq3WtRp3RiQukmiafFGdRo/XwHbQivQqLITRF+t0MdfB/BA9QoJHu/CVl2OZqCaNQLuxX4S0QanlgySJnkXSSa8OMNSHae+tCAthlZgd9o8gTSkuxDtkKD0q8kHW912oWMAswDPku49f16kKW2kTcGU6r1DKqnO6BD4h7D0HEg74ZzW6GqHb2zY9BZyOvRJC8obj+jkTTBPxl6jeJbGZIgyXInRaz32wjaktuQCK4HFWwVhcsYGwWvLOU0+r4CfA/tPTyClEn89FLY3J+AZN+b0QoFNGjtjMyTHWLSfBRdj2LvS3lhuwlfUbahRrWFYoWVQahhraZw/s4mbLXl3fpgryhaUEN8nNLKXm06S9EAugApqjwZEDYP5K1tHEPBfJ7XSR1Cob5uolBPwxhC8FnQB9CZrhVohd+GJn9PozN+Qym20tKO6u5idFvC8vjFyBT1LPcC4Ftl0p+OTLslTh4Gt4VI/BV25uxVCgdQN6IGMJLsNuBq8K4pCdoot/Esh7SiTi2I5RQ6C0+sOBz4rvn/bWSL8mXgp8g6y2MB8cQRSU5Eg2U3Wll57Idmk2FhD0Qd1PM+93ZKjyscBHyE4Dv6qk1nhsl7C8mJIquh0dpGO5qg2YynoMn6CtGDWyvqrIN+sx8gCzYnoFX5GjQRvA0pWq2geNLUDtyJtG+DtIyj2IQmmP492l0ovnIoCepZ7qsovydaD+WOLL3fAXE/0caCbQ4kWPyUd0agCri4nMcG50VKFUMuQ+9mXKn3qpiItCdPQ/tvYeeDsoBrG8XshupC3MuIw/BuUZgMnIQmQGFWT8JYSalJuudJXqGkEmpR7mqoVqEkb++3hI1IJFaOkWjvJ+kN26TYAPwx7UykTC/qzHeloIa+ltqpmY9F73me+b4vWnF21ij+WuPaRjHeedCzBhjPVCTSHWK+fx+dxwqyPRvGbCSuPwNdWbQIdeJBWotZoRbljksbWiW2ozOw88z/e8QMn8f3W4Snljqbwl7BOkrNPA1Gy+s4DT2v3I7sOzYzJyErLX1IZFpLRqKB8mc+99uIPqSeFq5tlHIaeif7DTCeORQs14DE00vQnuVQ49Zl0tozIp65aEL2rgnbEeE3C8Qpd63oJFg790bLTxfR7zhv77eIWahwK5ApmL2QRpy979IC/IbCVSCNytnoXcSd2TgaG9c2SrkGmYZLQstzPloN50FvIa809Du+CCk32OqzJ1NsY/EwNJPvsf4GOnPLIrsiefipaWfEkQlc2yilm9LrjurFk2RXXN0oNPQ7vhX4lc9tPimplmaApZReFupoTlzbKMZTJjky7Yw4skHW73Nrp3S/YxL5t8tXLVejO8tGpJ0RR+q4tlHMseig/UNpZ8ThKIenluq/IHQ91VmQbhRWoY1fR/Pi2kYpq4m2cepwZIapqAEPs9xGIdFDu+X2KIXLMG+gvAHgvHMEOu/VWs6jo2FxbaOY6cgW4qC0M+JwxGEOpXeMzURnG+xKPAuZlTmfUivRjcpi4My0M+FIDdc2illO+C3wDkeuWQUsQ6aWmoHh6A6jcua4HI5GbxunoKuXHI6GYzKyg3hn2hlJmAnoUsg8Wm53JEOjt429kW1DJ6J3NBy7oY3kccg47r7R3huOA3DiGEcwzdA2vonbZ3M0IEOQtXrP5MoJyFSSw9HsuLbhcDgcDofD4XA4HA6Hw+FwOBwOh8PhcDgcDofD4XA4HDb/BzsvnRFz+4zxAAAAAElFTkSuQmCC\n",
|
||
"text/latex": [
|
||
"$$- \\frac{2 D h_{t} {u}_{k,1}}{h_{x}^{2}} + \\left(\\frac{2 D H h_{t}}{h_{x}} + \\frac{2 D h_{t}}{h_{x}^{2}} + 1\\right) {u}_{k,0} - {u}_{k - 1,0}$$"
|
||
],
|
||
"text/plain": [
|
||
" 2⋅D⋅hₜ⋅u[k, 1] ⎛2⋅D⋅H⋅hₜ 2⋅D⋅hₜ ⎞ \n",
|
||
"- ────────────── + ⎜──────── + ────── + 1⎟⋅u[k, 0] - u[k - 1, 0]\n",
|
||
" 2 ⎜ hₓ 2 ⎟ \n",
|
||
" hₓ ⎝ hₓ ⎠ "
|
||
]
|
||
},
|
||
"execution_count": 15,
|
||
"metadata": {},
|
||
"output_type": "execute_result"
|
||
}
|
||
],
|
||
"source": [
|
||
"solution = Ω.solve(eqs[0].subs({i:0}), u[k, -1])[0]\n",
|
||
"# (-eqs[1]*2*h_x*gamma).subs({u[k, -1]: solution}).expand().collect([\n",
|
||
"# u[k, 0],\n",
|
||
"# u[k, 1],\n",
|
||
"# ])\n",
|
||
"(-eqs[1] * h_x * 2 * gamma_val).subs({u[k, -1]: solution}).expand().collect([\n",
|
||
" u[k, 0],\n",
|
||
" u[k, 1],\n",
|
||
"])"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"metadata": {},
|
||
"source": [
|
||
"### $i=I$"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 17,
|
||
"metadata": {},
|
||
"outputs": [
|
||
{
|
||
"data": {
|
||
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAb0AAAAqCAYAAAAnOoULAAAABHNCSVQICAgIfAhkiAAADgFJREFUeJztnXmwHEUdxz+PQEIOAiRRIkFCImDUAC8JIYFgeAmhKh4ccqUUweUoDEIpUCKHmkokBaG0OAvFg4iIiGCkQORIGUmASHGEvADyoCTwMBqjBBQSUUJ4zz++PbW9szOzs7OzOzu7/ana2p2ePmf7N9396193g8PhcDgcjqZgCfD7rDPhcLQ5Tg4dLcMOKcZ1KfAU8DbwOvBbYGKAvxVAv/lsB/4B3AfMCfDbCXSnmEdHKTtmnQFHXUhbFltRDl3dd9TMQ8DpSLgOAO4GNgEjfP7+BSwARgP7ADOBW5HgfcHndxNwat1y3N6cDkzPOhOOupC2LLaiHC4GhmSdCUdrMQx4HzjacvsIEqiuAP/3AK9Z16ON33nA/cB/gPXArIg05wDbgIGW214mnglV5b61OQ24MOtMOBpGLbKYRA7zwL7AnaSr7XK0OR9CwnK45XYy0AcMD/B/jvHv9UbnmutHgKNQJb0HWBuR5kWUq2GORoLqKrc4BI0EOrLOiKNh1CKLSeQwL8wHrsg6E47W4U4kGAMst6uAP4f4PxMJ1y7m+hLg36in6XEq8FfrejxwjHX9C+AWX7wLgCes62OBayvmXlxIcqG4G6mPfp0wfD0YCDwPHJR1RhwNpRZZjCOHeWUHNPc5I+uMOBpHvUY/V6Ne5QlIreIxGVgTEuajwEZgi7nuRBPwmyw/+wIvW9efAj5uXU8C1vni9U/CHxjgJ4yJqJFIwnVIjdhMXAC8RPzyO/JPrbIYRw79nc+80AdcDtyEM2xx1MA1wN8JnkN7A6kg/ewEvApcb7m9CHzF5+9eiqO0I4A3UW91LTAUWaD55xo2+OK5Cwn8cKSmOTuiLE+iRjIpXTTPSG8osJn8z8XYXECwAVRe00mbNGSxkhwCnItGhHmkA3gBOCPrjDjyyXWoR/ixgHvj0Itjts+9A7gRNWBjjNsQ1Cs91Od3A1CwrlciqzPQqKwfGa54zDRuh1luz5v8rUZzFGF0IPXkwAg/leiieRq98yntnfsZCZyF1LIvA/8F3gIeQ+quKK3Acoqm7/2oB70FdUjuQCrlsDlEL2xXRPw3GT/+DsptpGuktAeqdzfUOZ1GkIYsxpHDoM5nFpyI/rdH0VKNfvS/xeFilP8BlTw6HDY3oso2G+n/vc8wc/9Eii+O0cB+wEnAKjQCmWnFNR0Jmy1AI034TsvtFev3aPSyPc5cH4xUeX1WHgajhmwt8IkK5RmPeoA23ajR9H/2DImji+Zp9J4Fvhdxfz56vhvR3OiVwFI0n9OPyhHWcL2BnvMiYKH5vgZ4APifCb+C4nytP+x2ol+Wa00ck33uL6DGNS01/dkmHf9oOO106k1ashhXDldS7HxmRTfK1xagh+oavb2N/7n1yZqjVekP+Sw096+03Lajl93j6AW5hy+u+UitYuNfjrAXGoXYXGLi3YBe3Isonayfisyu16F5iyiOQQYAtdBFOo1egcqjoSgOJLhnbzMbWbr6X+yjgb+Y8CcEhBtv7vn/Lzv8w8aP/3l6YaPmGIcA76HGcyfLfSh6IfvrQC08gF76do+/HunEoUDy/zwtWYwjh1Da+UxKgdrq+CzUeHeYOKpp9EDlvD1h2g5HQ5iB5ueq4Uy0W4U3iT/MureConoV4DJk+VkLXTRHo3cppSPearnMpO9X+4FM3/tRJyOMD6BlI/1IteYP+6OIsJ80fp70uc8w7tcjY6bbkDpvK/BHYFpEnEEMB94FflrndOJSoLb/vFEEdT6TUCC98nZRfaP3fTQ6DjNoudTEeXzAPW+keHcV6TUj7VDG3KhrgngejRSeo9SCM4oDTbhnUCVfatw7kEXam5bfWiw3QXsV3gV8Gpl3++dFGslMNOLdmjD8e+Z7e8C9g8330xHhX0cjCSjdBcYL+wTheP798XuqzrHm3q5oN5FV6FnfS7A6NYzPotGLX6jTTqfVGIuMZfJON/of/Sp0jynmO6ieTzXfz6SdqQbTDmXMtZnuWxT/pLh8zfp9s/mA5jaWIeMNj1ot9YL2Es2KaWiCPwk7Ulx68WDAfa/hCjN/93jDfO8WEHYe4Val3kjKL4jefz8dmeTbwrgM9VY7iV/uz6HR6PI6p9Nq2J3PeZTPg+eF58z3YZRrFUD1YDNS9fvx6nHeF+y3Qxlz3eilSQ+tuy3XKGB3ki8mXoJGvfejnVxsOlDPuJ/KwrCr+fYaPy8sRFvReoSN9AqU9z57zPfOMeL1/M2laHhTr3RakSSdz2bEe9HvH3BvBDLU8dd/D69BsOvHEvRc4tTtZqAdyhhIL+GT4EGfanTmjnj0Ut1/cEuF+KYZf99KkJevmrA9lG9WDHpBRBmx2Lxi/Hq7wXhhHw8Noca6H3iHUuOSnZHKdX1IuNtNuPE+9+8SLNRHE7wOL2k61dJLuv95s9NLfcvbRfXvpwHIYCmofswx8S0OCfsmOqHC5kFU35KygmjDml7SfVdnUcY4zERTCH9D+SskiKPkWfpHeusp7+lGsTFBBhzRXEupChCkPjsW+Bmq7DaVjnwZZb63RPoq5zy01usF4EhK5zs94szngRq4cai+eJaaQT1HP94oq5vS3UQOQnXXr4q0w72FFlnbHIL2kPRzPLJI/J3PPWk61ZL2f97sNGN530cyMirgnjeSDVLhj0edM7/qvxPN/SZlMuENEKT/rs6ijHEYhlTot9aQVqVn6WhCCiS3bPMsJL9cRZjzTZjngA9G+Lva+LugQnw3G3+XBYQ9MyLcRRQtJ228DZHPCgizC7JUfdhyG4gaNbvn6809DUDzGA8ExFVtOiciC9Cxltt16AXlX5ZTiQL5sN5MiwLZWm+CrHJfCnC/w8Q3NuDeeeaevUevdzKFZ2A31MTxDPHWM0adgFEvGl3GJGyl+pFe2bNMy3qzmmG2+zSWQeZ7W0z/F6OF5d3IuOSfEX7jGLF8EW3x1Evpkoc4I71J5jvMiCUo3UlovtC+t52i9ew0dOqAt8nwTLTgOsgUu9p0lqGOgqdK/jrweTRf6FcN1Yus63ZW9TwNtlGUF5sJSM39ms99EMXOpF2PO5FR3EtoPfCTqA7OoHwUG8QU9PwaaSnZ6DI2irJnmVaj19HEn72ReXkPeiGd3AR5aiSeRWqQMPv5NpqcXoNUmpsj/O6AXvx9BAvnSKTvvxXtgnMsRRWrF3Yb8KeINDz1ZpARy7sELynxGio7T32ooduCdtXfZPIEstrsQ/uwBqVfTTr9aDRbQBslLAA+Q/hpBvWgmeQmbwyi1ILbYxvaGGE/y20oWvLknUhvG3J1omd2HFrL+WPU+QuKO4gpSDvwdtyMp0Cjy9goyp5lO1hvbkfqum40JF9D8TDMdsArZyULwy8B30FzG48iIxY/vRSNCiYgffsWNKIBNWa7oy3eDjVprkZH0bxqxeOFXUv4CHQYEsCtlBrKDERCuI7i+kGbsNGZdwKHPQLpQEL7OOUjsaTpLEcN62JkIPNUQNg8kDe5OY7iFoTeMUiHUqyvmynW0zAGE7yW9SG0Tu0RpBEYhjqFz6I1ikMo3ZWmE9XdpWhnp1XxiwFEn4BRL+pZxsXANyukPwttZ5c2WTzLurMELQSPyzrgw3XKS70okFzHP8mErTSRu5DK6qqVlv/TAu6/gybMV6M5O3ujbwLC/iQiP95OKH7DE09d8cOQcD2oIfaPNJZRPjc41cQVtFwlaTqzUcPwPuELneNQIPl/Hodmk5sCtZV3IdF1t7dC+AGoob8v4N7OaG52I6rjTyOV325IS7DS5/9FVLffQZqEagk7AaOe1LOMo1BHN+ozJEYek8zpZfEs686DRG+kbDOF2nZZySPDkdAvreSxxXmZcoOUK9CzGZdSGgcha84zUCMbtuapGXByU8oYVBfiHjAdhncyxVTgFPSirqbzM87kI2qf3KyptYxJqbbRy8OzTMQmpD6rxAg0fxQ2+mhlNgJ/yDoTGdOLXvJ7UjSX7yE9c/ix6Dl7+7VORD3krpTiTxsnN6V461nPqTGe6WjEONhcX47WmI0JDVGKdwKGf0lHM1FrGathGFKldqJR5QLze+8YYfPwLKvGM5udR3G+YT3lW1oNQmqyOELeityF9r9sZ05Bu9L0AT9IOe4RqAH1q0J/RfTi+6xwclPOGeiZHFBjPPMp7tQDUn/fieaVPBVewaS1T0D4KwnfDKFZiFPGtOgiWF19i+WnQPDzzMOzrJq5FOd9jkKbRt9DqZVRB/BLiseqtCPnoucUp3fkaH2c3JRzA5r/aYTV6SI0em4HQ8JG0FbP8xJ0yOloy+1USveZPBz17rutT629ubyxJ9LBn551RhxNgZObctZSfqxUvXiK5lV755G2ep53AD/3uS2iPqaveWc5tR+K62gNnNyU4hmxHJl1Rhz1J8/n6YEmNP1zJpPI/96E9eB6dGbc8Kwz4sgcJzelHI82EHg464w4HFF4ZrP+w1k3kGwn7nZgDZqEdrQvTm7KWUf0HrAOR1MwHQnvUMttJFJTdFpuqykeRHozlTdHbmWOQOvVBlTy6GhZnNyUMgvtHzkw64w4HJWYT/k5bnPQtlZ2BZ6LttW5EO0P1+4sBc7OOhOOzHByU8oqypdqOBy5Zw06SHCnrDPSBOwCPEbwuWEOh02ry81paLs8h6OlmIo2O/5N1hlpIiagwzrzuBO+ozG0utzsj87bc6p+R0sxBk1Sj0NHwEyM9t5WTMapdRzBtIPcfAM3j+doMQajc55mmuuT0LZQDocjHCc3DofD4XA4HA6Hw+FwOBwOh8PhcDgcDofD4XA4HA6Hw1Fv/g/SlLVoSWO7jQAAAABJRU5ErkJggg==\n",
|
||
"text/latex": [
|
||
"$$- \\frac{2 D h_{t} {u}_{k,I - 1}}{h_{x}^{2}} + \\left(\\frac{2 D H h_{t}}{h_{x}} + \\frac{2 D h_{t}}{h_{x}^{2}} + 1\\right) {u}_{k,I} - {u}_{k - 1,I}$$"
|
||
],
|
||
"text/plain": [
|
||
" 2⋅D⋅hₜ⋅u[k, I - 1] ⎛2⋅D⋅H⋅hₜ 2⋅D⋅hₜ ⎞ \n",
|
||
"- ────────────────── + ⎜──────── + ────── + 1⎟⋅u[k, I] - u[k - 1, I]\n",
|
||
" 2 ⎜ hₓ 2 ⎟ \n",
|
||
" hₓ ⎝ hₓ ⎠ "
|
||
]
|
||
},
|
||
"execution_count": 17,
|
||
"metadata": {},
|
||
"output_type": "execute_result"
|
||
}
|
||
],
|
||
"source": [
|
||
"solution = Ω.solve(eqs[0].subs({i:I}), u[k, I+1])[0]\n",
|
||
"# (eqs[2] * 2 * h_x * gamma).subs({u[k, I+1]: solution}).expand().collect([\n",
|
||
"# u[k, I-1],\n",
|
||
"# u[k, I],\n",
|
||
"# ])\n",
|
||
"(eqs[2] * 2 * h_x * gamma_val).subs({u[k, I+1]: solution}).expand().collect([\n",
|
||
" u[k, I-1],\n",
|
||
" u[k, I],\n",
|
||
"])"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"metadata": {},
|
||
"source": [
|
||
"## КН повышенного порядка аппроксимации"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 134,
|
||
"metadata": {},
|
||
"outputs": [
|
||
{
|
||
"data": {
|
||
"image/png": "iVBORw0KGgoAAAANSUhEUgAABYkAAAAZCAYAAAB0FaHIAAAABHNCSVQICAgIfAhkiAAAD2tJREFUeJztnXuMHVUdxz8tUPpCWmiEIOFREEpj7NJaLYLN8vBBtGp5pBEN3gBBrNUgwQeIuCgBo0YFQ1QgZCFiQGmqNgaoURcpYC1QQKTIQxYi5SGISBGoSP3jd2727uzcmTMz55w5c+/vk2zu3pkzM+d+f7/5zsyZM2dAURRFURRFURRFURRFURRFUYBBYHvH34O11kZRFEVRFEVRFEVpKlcDzwIz6q6IEjWLkPaH0wJu8/NmmycF3KZSL+pHiit6wbPmML79d3taoUEzYwQYAlY52riiKL3D7ogZrgEeAV4BXgTWA6cCk+urmqIoinqUI1THelH942Mdcp00mFHmR6bM6SEq1AAWA28AZ9VdEaURrAGeAmbmlNsD+B/wA/O9rF/+BNlf51WqdTyoR2WjfqS4poxnVTm/c+1Z05F23yFglJxG4iFHG1UUpfc4A/GJLcC1wMXAVcC/zPQbgEm11U5RlH5HPcoNqmO9qP7x8TzwOtk90DYh8VkYpEbxsw7J2Wl1V0RpBO9E9p9zc8qdbsodab6X9csHgJfonZtu6lHZqB8prinjWVXO73x61gjaSKwofU2L/DvN3TgKWMZEc9oTeMKs9/gKdVMURYHyPqUe5QbVsTot9FjbK8xFNL83o8x04L/Aq8BOISrlmRbl8xfgIKTX3uWO6qOMp0W1+MTKZuBxshtBbgSeA3Yw38v45QykZ9/6ivWNhV73qBbqR0qcFPWssud3vj1rhI5G4pB3zpYjxrQB2LdLmQsY/wjEOeb7cSll9zHz1ritplIBjXHv8ztgLXKg7eRp5BEmmHgA17yoj6Laq+5u0JyvD/UoN5TREXpfl1CU0V+198c7zOeGjDKLgB2B+xA/adOvcTkF6Q11vcdt9Ku2LohVu+vM9t/bZf6bkEaWtUiDCZTzywGkHeRuYD7yGPfTwFbgduBdZX9ATahHZZPnR/2gQSj6TcuinlX2/DqoZ4VsJH4W+A3SLfuLKfPnmukbgSvNtEXm886U8ovN590O66hUQ2Pc37RPOF5PTNe8qI+i2qvubtCcjxP1KDd00xH6W5dQdNNftfeHTQPMEvOZ1L9f43IMckH8R4/b6FdtXRCrdreZz24NLh8CpmDfyNTNL9vDLeyLaLArcA1wC3AY8CtgF8ttxIB6VDZ5ftQPGoSi37R06VlZ59dBPWtHVyuy4Dakt87zjBlZJ5ciAq5krGV9EdI1+4mU8u11bHJbTaUCGuP+ZUfgZPP/TYl5mhf1UVR71d0NmvPxoR7lhiwdoX91CUWW/qq9P9r6rWBsHNQk7V48aQ0w/RaXGUivp83Ayx6304/auiJW7Taaz6Vd5i9Hcmqdxbry/BKk4fQIxjdarUZ6QQ4At1psJwbUo7pj40e9rkFI+k1LV55lc34NNXjWIGHGJP4T8ga/Tj5stv3Djmm7mWlpIoH0+tkO7JWY/k0zL2Yuo7e62SfRGMcd4xbuxzH7jlnnrzPK+MyLJuQE1JcXNtrXsT/GvJ+4IGYvjF37Fm59qm6Pgt6IV5aO6t3jaRHuWNvPOembSYy9UMbm7+0dyzb5PLdF+fw9yCxr05BXll72GxtalI9P7Nq9gjxGnWQq8tKmn1uuJ+t4dZ+Zd2zKvAvNvG49A2MjtEfVsR+18OdHTdGgCTT5mFcFF56Vd53i27NGKPniulHszWc7Mk5GGteY+Xub79OAvwH/AGZ3lDvGlLuwy3r+CTyTMv0m4NtdlnHBUqQ795NI/Vol1jEbmBlBnUZxE9MkGuPmxni4RN0+Z5bdjBwcuuEzL3znBITJi1Hq2yfr2B9j2U/Aj/Yxe2GTtR8uWLcYPAr8xmsUP97RSZ6O6t31HWtjzMm6GMXtvtBuYLgjo8xsU+Y/jL1MC+qJy0rgMeTlVHcB77FYZhS3+XuYKZc1HnHRbSbjFKvfZBHLuXxo7Yr+7idJf9x6mVn+JIttZvnlVOSx7ke7LPtTs+xci+3EQGiPqrof/RbRuBujhPWjGDUoSixtQXUc82LQsqpn5Z1fh/CsEbMOoNhwE48iJxy2bOky/UHzOR/4OzK49f7AqcALHeXaXarvSlnHXMTs0u5SDCAXer6YCdxvtlF2Oy/kFwHE9EbJ791dtk6uYppEY2wfY1vK1un7wKzEtAHgI8DVSH51ck/Beq0CLgEeAI5GzL8bPvPCd05AmLyoc5+sY3+MxQvBj/Yxe2EsHgV+fSoWjwK/8fLlHW1sdOxn7677WBtjTtaF632h/Xhs1hiK7bEC72HsZVoQPi4rkDxZibwBfSXyVvX5pD/628Z1/r5iPqdmlKkap1j9JotYzuVDa1f0d09jLIc6OQ7YRvYTQZDvlwuQNpBuPUsXIk8VPWZR1xgI7VFV96OFdG9EhPB+FKMGbYZpVltQHecitlraUkbLKp5lc35dq2cNIq3HQz5W3sFxZjtnAgcgiXk78qhEJ9eZcvumrGOVmXdRYvqeZvp8832GWc/dwH7Vqz6BrRTvkbI3Usd5FmWHKR6PMnVyjcbYPsZlqBrjFlK/wYr1ONOs58/Amy3K+8qL0DkBceZFFjbah94f1QvHiFn7MrjQvkV1n4rFoyDueOVhq6N693hahDvW9ltOhuS7yG89NaPMF0yZSxPTQ8dlA3BFYtrDwMU5y6XRonz+7mWWXV9iWVua4DdZ1HkuX6d2eb97MjL2f7LH3A7IOKc35qzfxi8/bcqcljJvF7P93yemnwC8xnjNLjH13COnTr4J6VFVc+EAyuVtq+RykO9HMWswTLOuf0If88rmky02WlbxLNvz6xCeNUJHT+LJGZXxRWdPnUuQVvHPdFbKMA/pVv14YvrOwKfM/8k7ZgNIK/5fgYOR8QVfBw5n4l2nuliAPOrxUN0V8YjGuPdj/CXge8jd2yOBZy2W8ZUXTcgJqDcvbLQPvT/2w34Sqxf2g/YxeRQ0N15FdFTvdo+t/v2Uk6Gx6aV3qPlMvhAqZFymIL24kj2N1gHvzljOB08hw/Qc7HEbveo3IYhZu4ORG7PJ3qFLgd3JHgfd1i+zejsearafnLcaacg5z3w/G/gY8AHSH58PSUiPqpoLi5Bzqqy6uibPj/pBg1CEPheJQcuynlXk/Dq4Z9XRSPwI8pjDicAHkZfBpL3hcBuwE/DWjmkzgKuAt5nvyeUGEDE+ivT+uQL4BOndv+uiXcc38go2GI1xb8f4q8gA83chj0U8Z7mcr7xoQk5AvXlho33o/bHX9xOI1wt7XfvYPAqaGa+iOqp3u6WI/rHk5DDFxkOMncnIBdg24C8Z5dqPcicbYELGZQ7Scyl58fcM0lMrJNuBP5g6HehpG73qNyGIWbsl5jPZK2454je/7LJcEb9ciPSwuz9lXrsxJtnotB04F/G2LwPnI+cJDyfKDRPWA0N7VNVcWIT0ZPy3ZXkX5PlRP2gQitDnIjFoWcazip5f+/SsVIqMSeyKbch4GQciLebndSl3M7AY2anXIOODHI282e8pYDryQplOBpCkvAp5I/ktXdZ9IfCVnHoeiXS7ds0Cuo+dc675a7MzEuCzO6YdC9zqoV4u0Rh3j3Gd9XLBJ4GvIw0ptyIDrScZJf0lAr7ywjYnIN688I2N9j73xzTUC8eISXv1qGbFyxdldFTvdkdR/WPJyXbnk7QXuDSReYiWmxB/SGMmouNWxp5GaBM6LjGxGjgeeD9yA841Pv0mi6YfI6E+7Wx4H+J7nQ0rk5AGpDtI7wFXxC+nII1V9yI9HpNk9dhbB2xEcmCZ+T9JaA8M7VFVc2Eh6dr6JsuPYtKg6dc/oY95WVqG8uqinlX0/M63Z+UySJgxiQHWkn+XbSryuOcW5JG2O5Fu6rOQVvmRlGUeBK405ZdnrHsOYqpZf9MtfkeZMV8eAs7oMm835OK0/bcaGTuoc9o0D3XygcbYX73qHMdsiPw3n45kLO8jL2xzAuLNixDkae9zf0xDvXCMmLSPwaOgvE8NEZ9Hgb94+WKI4jqqd4+nRbhjbSw5uQnp0TO74Ppi5WRE6yszyhzOWE+1JCHjMgVpmDoxMf0yyjXqtag21uMU5OJ4Q8nl8/DpN1nU6UWdtCgfn7q0g+zfvSvSe/AXiemLkd96VpflhrD3y/bj6T/usq7NwEtMfBcBwFHAy0jjzsKU+RDeA0N7VNVceB4ZH7koLfz5UUwaNP36J/S5SJaWIa5nynjWEMXO73x7VpsRJg6vB4RtJF6PVNbmIGrLdLPOxcDHkaDmiVGVojvhDKSOS/IKGoZp1mDlnWiM/RFLjMvgOi/qyAmIMy/yiEl79cJqhNS+DLFoX4Z+jFcdqHfXi4+cnGWmf8tRHfuRqvvFBuDyxLSHKPfiOhecg1xbHppX0DN1+U0WTTlOutYu63d/FsmXIxLTLzLT96+w3aosAF4ETkEa7m5OKdPrHlg1F/ZH4niU+6pZ4cKPQmswTHOvf2yoomeIfMrTsume1ckIOY3E7b/k4wiumIS0dj/geL1LkDvo7bsr3wCeBN7ieDszkW7xA8gdj/PN//tYLHsYxS4+h7Ezhip18oHG2G0DQ9U6xYKPvAiVExBnXtgSm/bqhdUIqb0tsWlfhn6KV92od9eLj5xcBrxK+PFve4mq+8UK5JHz04BDkB5dW0l/03wIpiIvMFpb0/bbhPSbLJp4nHShnc3vnob0PrwhZfnN1DdcG8j+swWpN8ij328wsVdrr3tg1Vw4AWljmuW+ala48KPQGgzTzOsfW6ro6SufbLXsBc+aw8RezBPYD0nC9t8qp1Ud4yBTgWsdr/cMJCBtJgE/Q8bncHlyP0h6t/DhjjItM22/lDoWaXwfxs4YbOoUEo2xe2zqFDs+8iJUTkB+DFqk50S7nr5uvNlQl/Yt1Avr8sIW6lFF6ad41Y16d71oTsZJlbi0WYmMZ/iaWW6p+2oWYinwNaQHel240NUFgzTvOGnr1S266zdI/u8+BDnPS1u+TnZDfn/yUe/rkfFG+4mquXAx8pKxOqnqR6E1GKaZ1z+2VPFmX/k0iJ2WveBZ0xnf/jvktWY5rECEPjuvYIO5AHnDaB0vBowBjbGSRq/nRcw5UZf2MWsSgjpzvt+1L4PGqz9R7dNRXeJE4+IH1bUaqp/SRnNBNXCJaqn0DBspP7C60gw0xkoSzYmJqCb1odo3C41Xfaj26agucaJx8YPqWg3VT2mjuaAauES1VBRFURRFURRFURRFURRFURRFURRFURRFURRFURRFURRFURRFURRFURRFURRFURRFURRFURRFURRFURRFURRFURRF8cT/AaZ9yPk4uprnAAAAAElFTkSuQmCC\n",
|
||
"text/latex": [
|
||
"$$\\left [ - \\gamma {u}_{k - 1,i + 1} - \\gamma {u}_{k - 1,i - 1} + 2 \\gamma {u}_{k - 1,i} - \\gamma {u}_{k,i + 1} - \\gamma {u}_{k,i - 1} + 2 \\gamma {u}_{k,i} - 2 {u}_{k - 1,i} + 2 {u}_{k,i}, \\quad - H {u}_{k,0} + \\frac{- {u}_{k,-1} + {u}_{k,1}}{2 h_{x}}, \\quad H {u}_{k,I} + \\frac{{u}_{k,I + 1} - {u}_{k,I - 1}}{2 h_{x}}\\right ]$$"
|
||
],
|
||
"text/plain": [
|
||
"⎡ \n",
|
||
"⎢-γ⋅u[k - 1, i + 1] - γ⋅u[k - 1, i - 1] + 2⋅γ⋅u[k - 1, i] - γ⋅u[k, i + 1] - γ⋅\n",
|
||
"⎣ \n",
|
||
"\n",
|
||
" -u[k, -1] \n",
|
||
"u[k, i - 1] + 2⋅γ⋅u[k, i] - 2⋅u[k - 1, i] + 2⋅u[k, i], -H⋅u[k, 0] + ──────────\n",
|
||
" 2⋅\n",
|
||
"\n",
|
||
"+ u[k, 1] u[k, I + 1] - u[k, I - 1]⎤\n",
|
||
"─────────, H⋅u[k, I] + ─────────────────────────⎥\n",
|
||
"hₓ 2⋅hₓ ⎦"
|
||
]
|
||
},
|
||
"execution_count": 134,
|
||
"metadata": {},
|
||
"output_type": "execute_result"
|
||
}
|
||
],
|
||
"source": [
|
||
"eqs = [\n",
|
||
" (u[k, i] - u[k - 1, i]) / h_t - D / 2 * (\n",
|
||
" (u[k-1, i + 1] - 2 * u[k-1, i] + u[k-1, i - 1]) +\n",
|
||
" (u[k, i + 1] - 2 * u[k, i] + u[k, i - 1])\n",
|
||
" ) / (h_x**2),\n",
|
||
" (u[k, 1] - u[k, -1]) / (2 * h_x) - H * u[k, 0],\n",
|
||
" (u[k, I+1] - u[k, I - 1]) / (2 * h_x) + H * u[k, I],\n",
|
||
"]\n",
|
||
"eqs[0] = (eqs[0] * 2 * h_t).expand().subs({gamma_val:gamma})\n",
|
||
"eqs"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"metadata": {},
|
||
"source": [
|
||
"### $i=1,I-1$"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 135,
|
||
"metadata": {},
|
||
"outputs": [
|
||
{
|
||
"data": {
|
||
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAuEAAAAZCAYAAACfHvDPAAAABHNCSVQICAgIfAhkiAAACGxJREFUeJztnXuoZVUdxz8j2sydO+ATERMdzVcieCtHRoq4miZRminRHwltMnIgC4kofKRjSQYFaSCJhlxFQUUZcxJ8gE0zo5SvRFMnsbpGjmaiiTITvqY/fut0z913n7PX3nvtvdc++/uByzln7fX4rd/vu8/97dc6IIQQQgghhBBCCCFECW4EXgWm2zZkQvgEsAv4RtuGRIC0FRZpSwghhJgQ1gAfAN+toe99sWRhA/ACsBN4E9gKnAvsVsOYoSk7hw3Ay8CqBmyMFWlrPNKWEEII0WPuB/4DTNXQ9zrsrN124BbgSuAGN94u4A5gWQ3jhqTsHE5w2y9qxswokbbGI22J3pBgop1t1wxRIwmKcVskyPdtkVDe90diZyqvC2jPMCcDp7P0jN4BwD8wu8+uaexQVJnDc8CLGW27QoK0lUdCeR/1WVuVGJ70l4B3gT8Ch4yofznmzG8CF7r3Z2XUO9ht2xDMUlGVovEFxbhrKMbtId+3y9exM2231dT/g8BGLBkb5hXgWvd+dqi8jB7qpugchrkV0+SpqfI+aDg2bUF8+qpDW9ADfQ0n4a8CD2CXB76fUfcwV/4o8GvspnqAxzLqrnGvT4QxUwSgaHxBMe4ainF7yPftcgrwPvCHFsZ+172+N1RWRg9tkjWHYR5yr+lEqQ8ajk1b0C19ldUW9ENfi1gOvI0dXaX5LSbE493nvwP/HtHPldgRyhdCGygqXTIqEl9QjNsiQTFuiwT5vi0Syvl+GvsH/3Rge3zY3Y27Czgtta2oHtpi3BwG7Om2P5Iq74qGEyZLWxBeXwnhb8Wroi3ojr6C8gj2VOswZ2CT/ZX7vI/7fO+IPh5w2w9Mlf/UbYuZa4j78kZCtR3FJ76gGLdJwuTGWL435PulJJTz/ZGu3f2B7fHh527se0Zs99VDm+TNYcBO7PaCAV3ScMLkaQvC6ishfBJeVltQTl9dyD8gR8c3YRM7yH2eAv6GHY3s7cpOcXWuGNHH68C/MsrvBX5W3F5vPg3cDbyE2ZeU6GNvwi6XU8WmedfG92/Oo0+f+IJi3JRN8/QrxvK9Id+H8/2Jrs64e3aLjnezxxy+4+o+hyUMWfjqoSmby8xhwEssvqUgVg3P0w9tQXl9FbV5zsPmsnOApdqCcvqqO/+AwDnI7hkbt7nXY4B/YjfGH4qt9fiG2za4T+fxjPaHuQGyjl5mMNHUxSrgz26MsuO8kV8FMFHOA+trtOkqYK9U2QzwReyHA+ZT25706NMnvqAYF0Ex9kf7lxGz732Jxfc73euKMXX+Cvy3gH3bc7afD1wNPAt8BksIsvDVQxahbU7jO4cBUyz4GuLVcF+0BeX1Vcd33zBVtQXl9FV3/gEN5CBnYdn9BcBHMHE9zOI1Hm91dQ7JaH++2/aTVPkBrvwY93na9fMEsLqI9Z68TfEjlIMwG4/2qDtHfoKQpoxNaRLMxtmS7X3iC4pxWRTj0Wj/WiBm35ehTd8f6NptrTi+Lxe48Z4G9s+p66uHpikyB7BFHD7AEs4BXdJwwuRpC8LqK6Had9+AENqC4vpqOv+AADlI1rqMw0dWV2Nny7/lGg04Gnvi9cVU2+XAee59+onVGexI5y/AUdi9TO8Bn2TpUVdbHAfsAJ5v25Aa8YkvKMZdJtYYy/cLyPfheBm7/H5UA2P9APgFdmbwJGyVinH46qFJis4BzLfLWHxGtA8ajllbEJ++QmkLiuurC/kHpHSclYS/gD1V+2Xg89jN/X9K1XkH2AM4YqhsGvuFpGPd53SbGezI6EzsSO164ByWXoJok4GN6bUuJwmf+IJi3GVijbF8v4B8H45dwGZgP+DwGsf5Ifbg1+PYJfbXPNr46qEpyswBYK17/d1QWR80HLO2IC59hdQWFNdXF/IPSOk4657wd7BlYQ7HjmIuyahzH7ZG42bsCc9VmNOfwo4cV2IPB6QHPgJz4BnA78cYeQVwcc5ETgI25dQpynGMvvfpIhb/tOpybAf93lDZ54AtgW0KjU98oZ8xbsum0MQS4zTavxaIyfeToPs7sV/kOw1LTELzNeBHWMKzBXvoLM08Sx9g89VDE5SdA8BnXbvfDJXFpOE6iVVbEI++QmsLiuurC/kHeOp4I+Of+lyBXfrYjp1Wfwy7PLAXlt1vymizDVs0fgf2a0/j2A+7FDHub2VOH2Xu1XkeWDdi2z6Y0Ad/dwK/TJVN1WBTmoTq923lxRf6GeO2bEqTMBkxTqP9a4GYfD8Juv8QtlJC1prJIVhP/goSm0a09dFDE6yn3Bz2xM4m3pUqj0nDeSRMprYgnL4SyvtoPWG1BcX11XT+AeFzkP+zFTsy8THCh5WuvzXAVzHDPx6o71EUdc40ZuPavIqOOdp5cCwEoeMLkxnjMijG2Wj/qkaTvi9DDL4f/MT1x1q2I00demiSb2N+/VTFfmLX8Dhi1RZ0W18xaKsKteQgy4C3sKVlQrEWu0F+cCbrx9gaix8OOAbY5YoZ97cDuNS9P9ij7YkUE/IcfklCFZvqoI74wmTGuAmb6iDGGGv/qkaTvvclNt+vwB7i2tjS+FnUpYemmMLORN4RoK8YNexLjNqCbusrFm0VpfYcZPALUbdUMnMx67AF2wcsA27HbuAPuVPNkr/QfOLKVmfYuA1/5vBLEnxsapI64guTGeOQNjVJmzFO0P7Vdd/7Mktcvgf7IY3LsDNOMVCXHprio9h+uDpAXzFquAixaQu6ra+mtRWKWcZ/7yVka3hgZ66Ov8LSB6ImicuBZ8h+KLUPTHp8QTFuM8byvXwvFujD921IpOFiSF/xIQ3n8CjVF6MXcaMYt4d83x7yveg60rDoOtKwEEIIIYQQQgghhBBCCCGEEEIIIYQQQgghhBBC9IT/ATdDtjEGsW6yAAAAAElFTkSuQmCC\n",
|
||
"text/latex": [
|
||
"$$\\gamma {u}_{k - 1,i + 1} + \\gamma {u}_{k - 1,i - 1} + \\gamma {u}_{k,i + 1} + \\gamma {u}_{k,i - 1} + \\left(- 2 \\gamma - 2\\right) {u}_{k,i} + \\left(- 2 \\gamma + 2\\right) {u}_{k - 1,i}$$"
|
||
],
|
||
"text/plain": [
|
||
"γ⋅u[k - 1, i + 1] + γ⋅u[k - 1, i - 1] + γ⋅u[k, i + 1] + γ⋅u[k, i - 1] + (-2⋅γ \n",
|
||
"- 2)⋅u[k, i] + (-2⋅γ + 2)⋅u[k - 1, i]"
|
||
]
|
||
},
|
||
"execution_count": 135,
|
||
"metadata": {},
|
||
"output_type": "execute_result"
|
||
}
|
||
],
|
||
"source": [
|
||
"(-eqs[0]).collect([\n",
|
||
" u[k, i-1],\n",
|
||
" u[k, i],\n",
|
||
" u[k, i+1],\n",
|
||
" u[k-1, i-1],\n",
|
||
" u[k-1, i],\n",
|
||
" u[k-1, i+1],\n",
|
||
"])"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"metadata": {},
|
||
"source": [
|
||
"### $i=0$"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 150,
|
||
"metadata": {},
|
||
"outputs": [
|
||
{
|
||
"data": {
|
||
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAlAAAAAZCAYAAAD+KhDIAAAABHNCSVQICAgIfAhkiAAACXtJREFUeJztnWusHkUZx38teOyNcAmJBoXWIhQaI6fUKgjiq0GDX7wgpkETswLxwiVBgwaNqRCImmiMhXj3QyWgYCAY+ACWD3gp1VIo9UarBDkY2wIKqFSwCtQPz2zOdpmdd3Zm9vKefX7JyZ53Znb2mct/dm67C4qiKIqiKIqiKIH8AHgCWNyxHZ8E9gMf7NiOIbEayfMLujZE6YUOVYPdoDpUlAlkDfAi8KkK/42IsEeOOL5twnw00pbrTTwnRMYzZM4BrgV+CfwLyc/rx5xzK7AHWNKsaYoDlw5Vg5OH6lBRBsBG4B/Awgr/J4HncY+KH0AaiJMjbXkQeAaYHxnPkNmOlMUzwA78Gu43mnCfa9Y0xYFLh6rByUN1qCg9J2P8yNTF8cio97sV/stN/L9xxLEI+B/wH+BlgXaA3BxeADZFxBFCRlwe9o23AccB85A0+TTcII38o+iNM4SM5nQ4BA1C/3SYEWfPnNZhSuPeh1TeLcDSijBXcuD06mfN77MtYY8xfrcmtHFIaHn4cx4i8Jsq/N9gjlsccawGDgZ+i+Q7hJXBNKLLbcBKpLF5DNgLbAbe5E5KJ4Sks2nuBh4y16zDjUhdf4fFb6j6aAuXDkM1CPXrp2owHal12CsNpuxAPQHchUy/fcbiv9y4bwW+b9xWm+N9lvBrzHFbQhuHhJaHP2ciI85fV/j7NN6nmGMx70LKIF96WGriOhS4Dvg5cCpwG3CIw44uCElnX7nHHG0dqKHqoy1cOgzVINSvn6rB7qnSYa80eHDCuO5BesFPMlvZi1wDTAEXItO0IJnxd+AvlvB5HA8ktHFIaHn4sRgZce4A/l0RJk/7WmRK2kY+Ki0KO7QMQG4Gp3NgY3ALMvKaRjZl9oWQdPaVreZ4hsVviPpoi3E6DNUg1K+fqsHuqdLhnNfgvcA/S27vRqbWvlVwO8K43VkRz13G/6iS+5eNn+LHEMojI3yd/nhz7sYK/3nIptb9nn+vt8ThWwYgyw/7gXdZ4rna+NlmR2LJiN97USedbTLCf+8FwHPIkk2REH30QRttktGMDlNoEPzrZ1cahHgdptZgrD1FRsTpsPf3qBn8K2lVRlxn/F5tfi8E/gz8DTi8EO5ME+7qClueAh63uN8JfMUzPSGcgUzR7kLsy1o8f4b4/C/T1/IIzacZ6uXRhjHxnWrCVe1/yhv2XzniONyEeRY4yOLvWwYLkH0MD1dc54cmnuUOW3yYIW0e5vimM9Ye3wY4Z1TzvF3I015FQvTRdFsF8e0VyKzEI8jm6/uBt3icM0N7OkyhQfCrn21pEJrRYagGm7KnyIg4HXbRZ3Bqo7yE97AJ6Mtui9tOc1wJ/BXZ9PUa4Hzg6UK4fJr0fkscy5HCtvU0p5FK0hRLgN+ba4RcJ+b8FPlfpq/lEZpPXwcOs9jwHuQFfDMlv+1j4nvOHBdU+OfTwq519XzPxHZkD0cZ3zI4CdFk1WzYycjo8hGHLT6kzsMc33SWaaLex7CQ2XqRE6KPptsqiG+v1gLrkRvFJnO8AylD2zJJTps6TKFB8KufbWkQmtFhqAabsieGsg7bvkeFaiOKs5Fe4qXAsUjDuBmZhi1yowm31BLHxcbviyX3Vxr3leb3YhPPNmBZvOkvYS9hI7pU56dgEsojNp8ywqeZjzLnVj2y/DXjf74jjk+bMNdU+PuWwSdMONtbeA9B9i7cXXA7B9jHgWW2HumQvMJhr42M8DzM8U1n24zwH/nOR/K5PANRVx9tt1UQpqMtwPdKbg8BXwq4fkYzOkyhQfCrn3U1CP3SYWoNxtpTZEScDtu+R43VRhPvWCj2gNcjvfmLkEQUOQGZKn205P5y4GPm//KIYxrpkf4RWIGs9z4PnMZLe8eKoOXhZg8yvb2iwt9n9LvKHG1PhoB/GbhGWKuQRrDodwvwO+Dz5vdlwLnAWdinspvGN519ZgWSz+XRdV19TII2ppA6V55t2Qi8uWVbXDpMoUHwq591NQj90uFc0CDYddjmPaozbUwhBj6NFNq1FeHuNf7HFdwWAzcwu8Z6bOmcy5Fe4fuRpw0uTWa1nbkwAzUJ5dHlDBTAzeb815bc5yNv0N2H5GMVOzlwlFPGtwy2Uf0SwKpvc70TaVQuRz6VsIYwMuJHmr7pbJsR/iPfj5iwF5fc6+qj7bYK6uson/UpP+m0Drnh1CUjvQ5TaRD86meIBqE/OkytwVh7ioyI02Gb96jU2qhF/uKsx5F3aNi4yoTZA3wTWaPchawx7kaeurAtMz2FrEG/dYwN+dMSrr/RmDjmQgcK+lEeLrruQJ1rzr+o5L7SuLtGvkuQPRfjPvswrgymgP8y+/humfzbXLYR+mak0bQ9NeRLRpqG0qeutcF7kU2uG5B9EfuR5YDc7asV5/0IycujS+519dF2WwWT34Gy6TClBsFdP2M0CP3RYUoNxtqTUodt3qM67UDdbi6eOcIsQKYYdyNPTtyHTMMdhqx9/sxyzk7kJWDPIu+8cHEkMuXn+ls0Jo650oHqQ3m46LoDNYU0NuWX9H3YxOt68dxpJswvxlxjXBnkXyH/ToX/DuQGUW4g3o68N+cF4r7/lZGmA+VT19rgCtwdkhnLOYci0/0/sfjV1UfbbRXU11E+W/GBkvs3kJdG1iUjvQ5TahDc9TNUg9AvHabUYEacPVeQTodt3qNSa6MWm5CK5CN6XxaZONcAH0Iai9gPRo5jrnSg+l4efcin/BMBq8YFDKSJMjgJGVmdh+zD+GnCuENpIp1tcQlSB06PjKeLtgrCN5GXvz33J8I2kadAdRjPJGsQ+qPDTrQxD+mlP5g43lOQHmH+le6rkOm7VyW+zhJk49k00mtdZ/4/pqXzU9PX8uhbPi1ANife3kDcTZTBUmQkts78fh0yChslvEZdmqprbbAQyc+bE8TVVlsF8TpaiyxbXQCciIzw91L9PbWmUR3GMckahH7psBNt5C89uyFxvB9HplBz5gE/Rp6ISNnTHmGfatxQCJMZt2WB57dJl+WRMTn5BLLe/QVkY2JKUpfBEUjel5cabsL9ssGmaaqutcGJyHLDsgRxtdVWwXgdZVRrMOdCZClln7HR9hmbNlEdhjPJGoT2dZjh1kfr2lhrDLqs6Qt1yJXAH0j7LcGm6LI8JimfmmQImoDhpHOSUA3OMoT6OYQ0pkT10QFb6XapZFLQfFKUblENKko1qg9FURRFURRFURRFURRFURRFURRFURRFURRFUQbL/wGv1jgdcREc6gAAAABJRU5ErkJggg==\n",
|
||
"text/latex": [
|
||
"$$- \\gamma {u}_{k - 1,1} - \\gamma {u}_{k,1} + \\left(H \\gamma h_{x} + \\gamma - 1\\right) {u}_{k - 1,0} + \\left(H \\gamma h_{x} + \\gamma + 1\\right) {u}_{k,0}$$"
|
||
],
|
||
"text/plain": [
|
||
"-γ⋅u[k - 1, 1] - γ⋅u[k, 1] + (H⋅γ⋅hₓ + γ - 1)⋅u[k - 1, 0] + (H⋅γ⋅hₓ + γ + 1)⋅u\n",
|
||
"[k, 0]"
|
||
]
|
||
},
|
||
"execution_count": 150,
|
||
"metadata": {},
|
||
"output_type": "execute_result"
|
||
}
|
||
],
|
||
"source": [
|
||
"(eqs[0]/2).subs({i:0}).subs({\n",
|
||
" u[k, -1]: Ω.solve(eqs[1], u[k, -1])[0],\n",
|
||
" u[k-1, -1]: Ω.solve(eqs[1], u[k, -1])[0].subs({k:k-1}),\n",
|
||
"}).expand().collect([\n",
|
||
" u[k - 1, 0],\n",
|
||
" u[k - 1, 1],\n",
|
||
" u[k, 0],\n",
|
||
" u[k, 1],\n",
|
||
"])"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"metadata": {},
|
||
"source": [
|
||
"### $i=I$"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 151,
|
||
"metadata": {},
|
||
"outputs": [
|
||
{
|
||
"data": {
|
||
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAnIAAAAZCAYAAAC1gsMlAAAABHNCSVQICAgIfAhkiAAACLpJREFUeJztnXuoHcUZwH9JNeYmER8UFLUmJtVoKHpjTFUUPYoW/aetaUC0INta6husqGgpPlB8oIgP8IV/pKJWS0Kk+aM2QXzUd2xM6yMRUa9iTH1Ejc+q0fSPb5a7Z++c3Tk7s49z9vvBZe+Zx7ffzLffzO7M7CwoiqIoiqIoiqIolfFn4H1geo06/AHYCpxcow5tYwFS57+rW5EWoj7XTtTnFEUJzkLge+A8S9xKpNHpZOS/3aT5vace9xg5+3rKaTuLgVuAfwKfInV6T0b65cBGYEb5qikG9bnhQn1OUZRaWQl8AoxY4jYBW8geNXgBabgO9NTjFeAzYLKnnLazFrHHZ8A68juVn5o0fyxfNcWgPjdcqM8pilKYiPyn9yz2QUYG7rTEzTay/52RfxrwLfA/YNuCOoB0Wt8BT3jIKEqEXx02jaOAvYFJSJnyOhWQzucttEN3IUJ9zpeIZvlchJ8+Q+VzIRU6AblYnwVm9khzOd3Dyxeb34ssafc0ccsD6qj4oTaun98ijc8DlriDzPHZjPwLgG2A/yC2hGJ2HUXajzXAPKQR/C/wOfAUcHB+USqnSDmr4BHgNXNeV+5H/OdYS5z6XFjU54qjPleBz4W8kXsfWIUMQV5oiZ9twlcDd5mwBeb4vCX9QnNcE1BHxQ+1cf0cgzyVP2OJc+lUDjHHpD2K2DWeIpppZO0A3A08BhwK/A3YPkOPOihSzqbypDnaOhX1ubCozxVHfa4Cn9smoKwnkbvvTYxf3EluBqYAZyLD1CCF/xB425I+lvFCQB0VP9TG9TIdeSpfB3xhiY/r80Rk6sBG/NSebHCK2hWkkzqc7kZqGfJkOoosJm4KRcrZVFab4xGWOPW5cKjP+aE+N6A+9xywORX2c2Ro8bZE2M4m7KEeclaZ+N1S4deYOKU+1MbFiSi+tmMfk3elJW4Sshh7q+Pf/hYZrnYFmSbaChxvkXOlibM9uYYgwm99TD/lrJoObut1AL5CptaSqM9NJEJ9zpeIZvmcrz5JOgyZz43hflH2KvjdJm4P83sEeAP4ANgpke4Yk+7KHrp8BLxnCX8IuM6xPEU4Ahmi3oDoFxWQ8TBwX0N0GsPfpmmabOPQdQ/F63+M/up+SY68Q00621qduMN5OiP/TibNl8APLPGudp2KrHt5vcd57jNyZmfo4soYYesQ3MsZQh8Xf0rS6SPfBuRtySRtb1fHUJ/zZYzB9jkXfZJ0GDyf6/KF9NTq68ibNa68awlbb47zgHeQRYB7AacCHyfSxcPE/7LImI0Y13aHO4pcFGUxA3jJnKPoeQ6kt1GL4KNTCJumabKNQ9c9FK//G4EdU2GjwC+QzUXHUnFrc+R9ZY5TLXHx8H3WOox4jc1aZM1PGle7HoC0HbZRivg8m4E3M3RxJXQdgns5bZThT0UZYfyaiGl7u6o+50/TfK4MfYrSFJ8ro5/rYhFyd3ouMAdp9J5ChqGT3G/SzbTIONvEXZUK39WEzzO/pxs5a4BZ/qpP4HP6f3KcQ7ghXxtFdApNU21cdt2Df/1HFNdxN5PXtv3ADSbu1Iz8F5g0N/eId7XrGSadbafz7ZG1Lo+kwhcDX9N9LdyE3BjtkqGzjQg/O7uWsw46uI0OTEbqOT1Co+3qRCKPvOpzQkSzfM5XnyQdBsvnJvhCGfuhJO+8b0KeIs4yJ06yLzJU/FYqfDvgNPN/+klnFLkbfhWYi8y7bwEOY+JdeV3En/MY5rfCmmrjYa/7jchUxFxLnMvowHxztL1ZBe52zXoCnY80zum4ZcCLwJ/M7/OBk4DjsE87lIlrOZvMXKSe06MP2q6GRX0uDOpz3TS+n5tiFPrYnOyWHumeM/F7J8KmA/cyPtc9J5XnIuQ1718hb8GcG0xrO0WeHK9F9qcpiyaMyDXVxmXXPdQ7Igew1OT/cSJsMrJD+deIbXqxnu6nwDSudl1D781Ns74F+TOkwbsI+SzOQksaFyL86tC1nHXQwW104Dcm3dmpcG1XJxKhPjdsPuerT5IOg+VzVfRzwPhGe+8he93YuMKk2QjciswVbwD+jqwr+QT7VN1HyFqAI3N0iN/iyfrr5Mgo0uCsMno2SacyaIKN05Rd91D/jdxJJv9ZibB55D+hzUDW6OR93ifPrlOAbxh/FT9N/C1I2wgGyHTKFuxv3rkS4d+Iu1y/VfFLZIH2EmQtzVZk+iYOu96S5y9IPf4oFd7GdjWPCPW5YfM5X30GyefS+PhCX6xAChZlpJmKDLG+i7zR8zwyDLkjMg/9qCXPemTTwC+RvWmy+CEy5Jn1Ny1HRpEGZxOyLqJJOpVBE2ycpuy6h/pv5KYgDWFyA9JTjMysDTUPM2kez5GfZ9d4WP+OHvHrkI7LtvblaGQvru/w++ZkhH+n4nL9VsVlZN8YjaXS74BMyzxokdXGdjWPCPW5YfO5CD99LmNwfC6Njy/0xRPIhePSMboyzchcCPwaaQx8P0CcR78Nzl7IRXB0KdoITbmRa5qNq6h7aEb9x5+DmZ+XsABl2BXkrbvNyOeOlgH/CCy/X8oqZxWcg9j/8ACytF11Q33OH/U5YRD6OSYhTwevBJZ7CDKsOWJ+X4EMX+4e+DwzkIWIo8jd8iXm/z0d8i5GKjn9qnSdOpVBE21cVt1D8+p/KrKwdkVguWXZdSbytHqJ+f0T5Em1E/g8rpRVzioYQepyaSB5bW5X+0F9zg/1uXGa2s91EW+SeG9guacjQ8gxk4C/Im/qhLzD72Afal2SSBOZsFmpvFfTe8PGsnWqkjptHFFt3UPz6h9kg9VLkUW1oSjDrjsjNk1PCz1A9kaqZVLW9VsF+yHTQrMCyWtzu9ov6nPFUZ8bp6n9XBcnGgXOr+JkNXE58DJhv1U7SNRp47bXfZm0wXehPeUcNNro2225FttSzlC00RcqZzX1DVG3Ha17RRlO1LcVRVBfUBRFURRFURRFURRFURRFURRFURRFURRFURRFURTg/4+VR9X8IeyTAAAAAElFTkSuQmCC\n",
|
||
"text/latex": [
|
||
"$$- \\gamma {u}_{k - 1,I - 1} - \\gamma {u}_{k,I - 1} + \\left(H \\gamma h_{x} + \\gamma - 1\\right) {u}_{k - 1,I} + \\left(H \\gamma h_{x} + \\gamma + 1\\right) {u}_{k,I}$$"
|
||
],
|
||
"text/plain": [
|
||
"-γ⋅u[k - 1, I - 1] - γ⋅u[k, I - 1] + (H⋅γ⋅hₓ + γ - 1)⋅u[k - 1, I] + (H⋅γ⋅hₓ + \n",
|
||
"γ + 1)⋅u[k, I]"
|
||
]
|
||
},
|
||
"execution_count": 151,
|
||
"metadata": {},
|
||
"output_type": "execute_result"
|
||
}
|
||
],
|
||
"source": [
|
||
"(eqs[0]/2).subs({i:I}).subs({\n",
|
||
" u[k, I+1]: Ω.solve(eqs[2], u[k, I+1])[0],\n",
|
||
" u[k-1, I+1]: Ω.solve(eqs[2], u[k, I+1])[0].subs({k:k-1}),\n",
|
||
"}).expand().collect([\n",
|
||
" u[k - 1, I],\n",
|
||
" u[k - 1, I+1],\n",
|
||
" u[k, I],\n",
|
||
" u[k, I+1],\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.7.2"
|
||
}
|
||
},
|
||
"nbformat": 4,
|
||
"nbformat_minor": 2
|
||
}
|