{ "cells": [ { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "# A script to illustrate the backward/implicit Euler method for systems or scalar ODEs.\n", "# Key points:\n", "# The spatial variables are y[0],...,y[d-1] for a system with d dependent variables.\n", "# t0 is the initial time, T is the final time\n", "# The vector \"f\" is an expression involving t,y[0],...,y[d-1] that defines y'=f(t,y)\n", "# Define f as below (scalar and vector examples shown).\n", "# Also Jf is the Jacobian matrix of f with respect to y, supplied by user.\n", "# Solution times are returned in row vector \"tv\"\n", "# Solution values are return in matrix \"yv\". Rows correspond to times, columns\n", "# to variables y[0],...,y[d-1].\n", "\n", "var('t y');" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "# SYSTEM EXAMPLE: Define f in y' = f(t,y). Initial condition is y(t0) = ic.\n", "# You must also define \"Jf\", the Jacobian matrix for f with respect to y (used in Newton's\n", "# method for solving the implicit equation).\n", "def f(t,y):\n", " return vector([-sin(t)*y[1],-y[0]+cos(y[1])])\n", "def Jf(t,y):\n", " return(matrix([[0,-sin(t)],[-1,-sin(y[1])]]))\n", "ic = vector(SR,[1.0,2.0]);\n", "\n", "#SCALAR EXAMPLE. Uncomment to test.\n", "#def f(t,y):\n", "# return vector([-sin(y[0])])\n", "#def Jf(t,y):\n", "# return(matrix([[-cos(y[0])]]))\n", "#ic = vector(SR,[1.0]);" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "# Set solution interval t0 <= t <= T, and stepsize.\n", "t0 = 0.0;\n", "T = 1.0;\n", "stepsize = 0.1;" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "#Call routine for implicit Euler's method.\n", "load('implicit_euler_sys.sage');\n", "[tv4,yv4] = implicit_euler_sys(f,Jf,ic,t0,T,stepsize)" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Implicit Euler Method\n", "time 1.00000000000000\n", "y[ 0 ] = 0.418858816516853\n", "y[ 1 ] = 1.29371030048058\n" ] } ], "source": [ "#Print solution values at final time.\n", "N = len(tv4)-1;\n", "d = len(ic);\n", "print('Implicit Euler Method');\n", "print('time',tv4[N]);\n", "for j in range(d):\n", " print('y[',j,'] = ',yv4[N,j]);" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlAAAAGSCAYAAAAsKv9nAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA3bUlEQVR4nO3dd5xcZb3H8c8vCVlCIAsEAgGCFKUkIEQIHQy9iEgVEUWKKFUQkWssiCjEQrvoInpBkGsBREVAgihIMUiXelFQipQA0nYhkE177h/PbNhsdjc7s+XM7Hzer9e8ZnPmnJnfckj2u895zu+JlBKSJEnquSFFFyBJklRrDFCSJEllMkBJkiSVyQAlSZJUJgOUJElSmQxQkiRJZTJASZIklckAJUmSVCYDlCRJUpkMUJIkSWUyQEmSJJWppgNURIyJiLE93PeLEbFnf9ckSZIGv5oNUBGxBtAEvNXDQ84Bjo6I3fqtKEmSVBcipVTch0cMBaYAHwVeB5YGpgHfSinN6ua4ZUr7fSyl9Fwnr28BNKSUbu2wfXngdmCHlNJLffaNSJKkulJ0gPoh8GFg85TS8xExCrgVmAF8KHVRXET8D3B3Sul/unj9ReDWlNKBnbz2NeB9KaVD+ur7kCRJ9aWwS3gRsTVwFPCdlNLzACmlFuDrwO7AIuGndNxEYHvgki5eXx9YCbiji4/+AfDhiHhvr74BSZJUt4qcA3VY6fm3HbZPA1qBT3dx3CnAJSmluV28vk3p+bbOXkwpvQ7cDBzR81IlSZLeVWSAmgy81XEOU0ppDvAEsFVpjtQCpUt8ewPXd/O+2wLNwIPd7HMbsE/5JUuSJPVDgIqI5SPiooj4bUT8KiKW6PD6dyLi18AaQEsXb/MmMIJ8Ka69nYB5wEMd3vNjEXFvRNwLfAKYA9xd2rZ5J+9/G7BuRHR8f0mSpMXqjxGo04Gvlh77AwvaBkREAEcCywFDyUGpM23bl++wfSvgnpTSvPYbU0qXp5Q2JY8qBXBWSmnT0uOuTt7/6dLzRj39piRJktr0aYAqTeCekVJ6kRx2AF5ut8uG5PB0c+nPXbUqmFN6buywfX3g+W5K2L70fGs3+wC8QR7JWmMx+0mSJC2ir0egVgT+t/T1IcA/gbvbvb5d6bnTCd7tLF167hiwViWHn65MBmYC93b35qX2CM0sGtAkSZIWa1hfvllK6TaAiFidfDfcVzv0ctqOPO9pOvA2XQe4ZUrPr3fYvjQ5+HRlMjC9mzv02msFGnqwnyRJ0kL66y68fUvPV3XYvh1wW2kO01Pky3mdWYEccJ7psH0eMLyzAyJiHLAmi79812Z5er4MjCRJ0gL9FaAmkedC/aNtQ0SsS76rrm3+0zRgtYhYKERFxLLA6sBfO04WJ1++6zixvM0i858i4vOliesLiYiR5NGn7uZTSZIkdaq/AtRoFh092qn0/OfS89Wl5/067HdA6fmyTt73aboOUJPII1T3AETEWsAaXSwH09a+4B+dvCZJktSt/gpQ9wKrtzXCjIj3A98EXqPU4DKlNB24BjgtIt5T2m8M8CXgFuDSTt73IfJlus68BryRUppdarg5FTiti323JM/FeqTM70uSJKl/FhOOiKWAC4Cx5Dvx3gKOAaallD7abr8G4FRgL/Lk8BWBPwFTSuvidXzfbcnhanRK6Y0Ory0PXAG8Sm6D8M2U0uNd1PcjYFRK6aBefaOSJKku9XmAKs05WjKl9E67bXuT17z7cErpul689xDgBeDTlb5P6T3+CRyXUupuSRhJkqRO9cclvBuAl0oTtdsC1SnAH3oTngBSSvOBi4CP9eJt9if3iprWm1okSVL96o8ANYncPPOd0hyoc0ufc2Afvf9/A7tFRFdzoRbnS8AJXUwulyRJWqz+uIS3M7AzsBQwhhym/julNKfbA8v7jKOAzVJKh5d53BTy/KmT+6oWSZJUf/plEvlAiIifAjenlH7aw/13BQ4Cjuikv5QkSVKP9Vcbg4FwBLBFaQHjnpgBHGZ4kiRJvVWzI1CSJElFqeURKEmSpEIYoCRJkspkgJIkSSqTAUqSJKlMBihJkqQyGaAkSZLKZICSJEkqkwFKkiSpTAYoSZKkMhmgJEmSymSAkiRJKlNNBKjIRkVEFF2LJEnSsKILABa7mnFzczONjY00NzcPRD2SJKl2FDK4UhMjUJIkSdWk7AAVEdtFxLUR8UJEpIjYuwfHfDAi7ouIWRHxZEQcVVG1kiRJVaCSEaiRwIPAcT3ZOSLWBK4HbgcmAmcC50fEfhV8tiRJUuHKngOVUpoGTAPo4Zzuo4B/p5ROLP35sYjYFDgZ+HVnB7S2ttLa2rrgz2/d9Vb+4l5g6XIrrkKrAqtQ0FVbSZLUWwMxiXxL4MYO2/4AHBERS6S06BzyqVOn8o1vfGPBn+/l3vzFjv1W48BbDtgQ2KD0vCEwAVi2wJokSVKPRGcBpscHRyRgn5TS1d3s8zhwaUrpzHbbtgKmA6uklF7oeMwiI1B/e4tVJ69K8/RmRi09quJ6q8J84BngYeCR0vM/gHml18excKjaAFgfaBjwSiWpZs2ePZu5c+cWXYb6yLBhwxg+fHhXLxdyPWeg2hh0TGnRxXYAGhoaaGholxgmlp43AGo8PwGwMfCRdn9uJYeo9qHqcuC7pdeHAuvwbrBqe14L76OUpA5mz57No48+yvz584suRX1kyJAhTJgwobsQNeAGIkC9CKzcYdsYYC7w6gB8fvVrAN5ferTXDDzKwsHqPOC10utLAeNZeLRqQ2AlnF8lqW7NnTuX+fPns8YaazBixIiiy1EvvfPOOzz99NO8+uqrCwZYRo4cWXRZAxKg/gp8uMO2XYB7U0pzBuDza1cjsFXp0SaRI2n7UPUQecTqndI+K7DoZcANgGUGpGpJqgojRoxgqaWWKroM9ZHrrruOmTNnMmrUKA488MDCQ1TZASoilgbe227TmhGxMfBaSunfETEVWDWldEjp9QuB4yLiHOB/yJPKjwAO6lXl9SqAsaXHLu22zwOeZOFg9UegiTzvCuA9LDpatS5QPSOikiR1auTIkcyZM4eWlhZaW1trL0ABmwJ/bvfnc0rPPwUOJf9oX73txZTSUxGxB3AucCzwAvC5lFKnLQxUoaHA+0qPfdttnwU8xruh6mHgf4HnSq8PI4eojncEvgfnV0mSqkbb5bvZs2cXXQpQWR+oW+hmhk1K6dBOtt0KfKDcz1IfWJI8CX9ih+2v8+78qrZRqxuAN0qvL01uq9AxWK3Y7xVLklT1qmExYRVhOWCb0qNNAp7n3dGqR4D7yCNWbV0lxrBwqNoBWHNgSpYkqVoYoPSuAFYrPXZrt30u8C8WHq26HjifHLo+AOwH7E9utyBJ0iDnLBctXts8qf2Bb5AX4Hmc3GbhSvItBWfy7lyqb5AvD1beo1WS6sall17Ksssuu+DPp512GhtvvHGPj48Irr766j6vqycmT57MiSeeWMhnF62qA1RTUxPjx49n0qRJRZeiziwDHABcAfwH+C2wEXA273ZQ/yrwAIYpSTXn0EMPZe+99x7wzz355JO56aaberz/jBkz2H333QF4+umniQgeeOCBbo9p26+zx5133tmb8utGVV/CO/bYYzn22GNpaWmhsbGx6HLUnRHA3qVHK/An4CrgAuAMctf0/UuPTbHRpyR1Yemll2bppZfu8f4rr9yxV3XP/elPf2LChAkLbRs9enTF71euefPmEREMGVLV4zmdqr2KVf0agA8BlwAvkZeO3hH4CbAZsAZwEnAH7/aokqQqN3nyZI4//nhOPPFElltuOVZaaSV+/OMfM3PmTA477DCWWWYZ1l57baZNm7bgmFtuuYWI4Pe//z0bbbQRSy65JJtvvjkPP/xwl5/T2SW8n/zkJ0yYMIGGhgbGjh3Lcccdt+C19pfw1lwz39UzceJEIoLJkyd3+z2NHj2alVdeeaHHEkssAXQ+AnfiiSd2+56zZ8/mlFNOYdVVV2XkyJFsvvnm3HLLLQteb7tced111zF+/HgaGhp45plnuq2xWhmg1L+WIDf8/DEwA7gZ2BP4JbA1efHk44FbeXdBZUmqUj/96U9ZYYUVuPvuuzn++OM5+uijOeCAA9hqq624//772XXXXfnkJz/J22+/vdBxX/ziFznrrLO45557GDNmDHvttRdz5vRsMY4f/vCHHHvssXzmM5/h4Ycf5pprruG9731vp/vefffdQB5ZmjFjBr/5zW969w2X6bDDDmP69OlcfvnlPPTQQxxwwAHstttuPPHEEwv2efvtt5k6dSoXXXQRjz76KGPGjBnQGvtKVV/C0yAzDNi+9DifvMjPVeRJ6T8gt0jYh3xH32Ry+JKkvvA28Pce7LceeZ3RLmy00UZ89atfBWDKlCl8+9vfZoUVVuDII48E4NRTT+WHP/whDz30EFtsscWC477+9a+z8847AzmErbbaavz2t7/lox/96GJL+ta3vsUXvvAFTjjhhAXbupobvOKKuVlf28jS4my11VaLXD5rbm5m6NChiz22o3/961/88pe/5LnnnmOVVVYB8nyuG264gUsuuYQzzzwTgDlz5nDBBRew0UYblf0Z1cQApWIM5d0+VOcA9/BumPoRsDx5PtV+wE643Iyk3vk7sEkP9ruPbts+v//97676PnToUEaPHs2GG264YNtKK60EwMsvv7zQcVtuueWCr5dffnnWXXddHnvsscWW8/LLL/PCCy+w44479qD48l1xxRWsv/76C22rJDwB3H///aSUWGedhfvZtLa2LjSvavjw4Qv9d6xVBigVbwiweenxXeBv5DB1FXneVCOwFzlM7UKesC5J5ViPHI56sl832uYHtYmIhbZF5Dtk5s9f/ATPtn27M2JE//6DN27cuC4vBw4ZMoSUFr6FurvLjvPnz2fo0KHcd999i4Sw9pPiR4wY0aPvvdoZoFRdgvzb3wfId+89Qh6VuorcEX1p8gT1/YHdgWLXkpRUK5ai0AXF7rzzTlZfPS8T+/rrr/P444+z3nqLSWvAMssswxprrMFNN93E9ttvv9j9hw/Pw/Xz5vV+UumKK67II488stC2Bx54YJEQ2WbixInMmzePl19+mW233bbXn1/tDFCqXsG7a/CdRh6CbwtTB5BHonYnh6kPAaMKqVKSFuv0009n9OjRrLTSSnzlK19hhRVW6HGPqdNOO42jjjqKMWPGsPvuu/Pmm28yffp0jj/++EX2HTNmDCNGjOCGG25gtdVWY8kll+y2DdCrr77Kiy++uNC2ZZddliWXXJIddtiB733ve1x22WVsueWW/OxnP+ORRx5h4sSOi6tm66yzDgcffDCHHHIIZ599NhMnTuSVV17h5ptvZsMNN2SPPfbo0fdbK7wLT7VjPeAr5Et8/ySHqmeBj5MnoO8FXEZeKFmSqsi3v/1tTjjhBDbZZBNmzJjBNddcs2C0aHE+9alPcd5553HBBRcwYcIE9txzz4Xuamtv2LBhnH/++fzoRz9ilVVW4SMf+Ui3773TTjsxduzYhR5tLRF23XVXvva1r3HKKacwadIk3nzzTQ455JBu3++SSy7hkEMO4Qtf+ALrrrsue+21F3fddRfjxo3r0fdaS6Lj9c0CdFlAU1MTTU1NzJs3j8cff5zm5mZGjXKYQR08A/yGPDJ1B/nuvR3JI1MfAVYorjRJA+/tt9/mscceY/3112eppbq5pW4A3HLLLWy//fa8/vrrCy3Xop5rO59PP/00b7zxxoIgt/zyy7ftUsiEqqoegTr22GP5v//7P+65556iS1E1ew/weWA68Bz5rr53gM8AK5Pv4ruQ3NRTkqQ+UNUBSirbqsBxwC3AC0AT+XeT44CxwAeB7wPPF1SfJGlQMEBp8FoJ+CzwR/Lo00Xku/a+AKwGbEUerarNVQQk1YDJkyeTUvLy3SBkgFJ9GA0cDlwPvEyebL4i8GXy2nyTyA083ymoPklSTTFAqf4sC3wS+B3wH/K6fKsCR5PnU30LeK2o4iRJtcAApfq2DPAx4GrgH+Ru598CVgdOAJ4uqjBJUjUzQElt3gf8EPg3cBK58/l7gYOBB4orS5JUfQxQUkdjgNPJQeoccnuEieR1+P5EN53LJEn1wgAldWVp4HPkrue/AF4Bdiavp/VLYG5xpUmSilXVa+G170QuFWYYcBB5rtRNwHfJy8d8mdzA8whc1FiqQu+84221g0G1nseqXsqlTUtLC42NjS7lourxAPA94AqgETiW3KxzTIE1SQJg9uzZPProo8yfP7/oUtRHUko8+eSTvPnmm1WzlEtVj0BJVWtj4OfAmcC5wNnkQHUouVHne4sqTNLw4cOZMGECc+fOpbm5meuuu46RI0fS0NBQdGmq0Lx585g7t7rmTRigpN54D3AecCr5Dr7zyQ059wO+CGxWWGVSXRs+fDjDhw9n1qxZzJw5kzlz5higBoHW1taiS1jAACX1heWBr5DbH1xGHpHanLz23inA7hQ0yCzVt4aGBkaNGkVLSwuzZ88uuhz1gVGjRlVFGHYOlNQf5pE7nX8XuAuYQB6ROggYXmBdUh2aOXNmVY1cqHcaGhoYOXKhO3cK+fXUACX1pwT8hRykriMvGXMi8BnA/5UlqS8UEqDsAyX1pwC2Ba4FHiU34/wyMA74EvBCcaVJkipngJIGynjgJ8BTwGfJk87XIPeReqy4siRJ5TNASQNtVfIlvX8DZwA3kMPVXuTLfYVfVZckLU5VB6impibGjx/PpEmTii5F6nuN5InlTwGXAP8iX+7bGvgtYA9ASapaTiKXqsV84HpyQ87bgHWAk4FPAksWWJckVTcnkUt1bQiwJ3ArcCewIXmu1BrkjuevF1aZJKkDA5RUjTYHrgL+AewNnE6+c+/z5LlTkqRCGaCkavY+4ELgGXJ4+imwFvmy3kMF1iVJdc4AJdWClYBvkkefziHPkdoI2A24Ge/ck6QBVlGAiohjIuKpiJgVEfdFxLaL2f/YiHgsIt6JiH9ExCGVlSvVuaWBzwH/BH4OvATsCGwKXAFU12LlkjRolR2gIuJA8vrzZwATgduBaRGxehf7Hw1MBU4jrwj2daApIj5cWcmSWAL4OHA/cCN5MeOPke/c+wHwTnGlSVI9KLuNQUTcBdyfUjq63bbHgKtTSlM62f8OYHpK6Yvttp0HbJpS2gbbGEh942/kFghXAquVvt6fgm7wlaQBU/1tDCJiOLAJ+Xfe9m4EturisAZgVodt7wCbRcQSnR3Q2tpKS0vLQg9JizER+AV5WZj3Ax8FPkgOVpKkPlXuJbwVgKHkmRftvQSs3MUxfwA+HRGbRLYpcDj5IsQKnR0wdepUGhsbFzzGjRtXZplSHXsfcA35b96r5F95jgReLrIoSRpcKr0Lr+Nlt+hkW5tvAtPIrQHnAL8DLi29Nq+zA6ZMmUJzc/OCx7PPPlthmVId2wV4EDgf+DU5WJ0FzC6yKEkaHMoNUK+QQ0/H0aYxLDoqBUBK6Z2U0uHAUuSeyqsDTwNvlt5vEQ0NDYwaNWqhh6QKDAOOA54ADgG+BGwAXIutDySpF8oKUCml2cB9wM4dXtoZuGMxx85JKT2XUppHvl/oupSSy6VKA2E08H3yiNQawF7kHlL/V2BNklTDKrmEdw55TtPhEbF+RJxLHlW6ECAipkbEZW07R8Q6EfGJiHhfRGwWEZeTfwf+cl98A5LKMIE8N+p3wJPkyeafA14rsihJqj1lB6iU0hXAicCpwAPAdsAeKaVnSruMJQeqNkOBL5B/9/0jeV35rVJKT1datKReCPII1CPkDm2XkudHNWEjTknqobL7QPUD+0BJRXoJ+ArwE2A8uU3uTkUWJEllqf4+UJIGoZWAi4B7geXIMxo/Ql4uRpLUKQOUpOwD5EWKryBfnB8PnALYx1aSFmGAkvSuIHcw/zvwNfK6eu8DLqaLrm2SVJ+qOkA1NTUxfvx4Jk2aVHQpUn0ZQQ5Qj5Mv6X0a2Iy8dLgkyUnkknrgr8AJwD3AgcB3WfheW0kqjpPIJVWpLcmLMV0K3AqsC3wdmFlgTZJUIAOUpJ4ZAnyKfFnv88B3yEHq57gsjKS6Y4CSVJ5lgDPJy8BsDnwC2Jp8eU+S6oQBSlJl1gJ+DdxMvpS3GXAoMKPAmiRpgBigJPXO9sD95NUwfw+sQ14iZlaRRUlS/zJASeq9ocBngSeAI8krZY4HfoPzoyQNSgYoSX1nWeAc4GFgfWA/YEfgoQJrkqR+YICS1PfWI1/Oux54AZgIHAX8p8iiJKnvGKAk9Z/dyaNRZwOXk5eFOReYXWRRktR7VR2gXMpFGgSWAE4kz486CDgZeD8wrcCaJKmXXMpF0sB6iByo/kweoTqHfMlPkirjUi6S6sD7gZvIPaT+DmxI7mz+epFFSVJ5DFCSBl4A+5K7mX8T+B/y/KgLgXkF1iVJPWSAklScJYEvkedHfRg4GvgA+fKeJFUxA5Sk4o0FLgHuBkYCO5B7SD1ZZFGS1DUDlKTqMQmYDvycHKbWB84A5hZZlCQtygAlqboE8HHyBPOTyMvCbFX6syRVCQOUpOo0krwo8R1AM7mb+bnA/CKLkqTMACWpum0O/I28FMxJwPY4N0pS4ao6QNmJXBIAS5FHn/4M/JvcS+pCetCGV5L6h53IJdWWN4EvAj8CdgEuAsYVWpGkYtmJXJIWaxny6NM04FFyJ/PLcDRK0oAyQEmqTbsBDwN7AZ8C9gFeKrQiSXXEACWpdi1HHn36DfluvQnAVYVWJKlOGKAk1b59yJfzPggcQO4j9VqhFUka5AxQkgaHFcmjTz8nz4+aAPy+0IokDWIGKEmDR1sX80fJjTf3BI4AWoosStJgZICSNPisQh59ugi4knyn3k2FViRpkDFASRqcgjz69DCwNrATcDwws8iiJA0WVR2g7EQuqdfWAP4EnA9cDGxMvmNPknrBTuSS6sfj5J5RdwEnA6cDSxZakaTesxO5JPWrdYC/AN8G/hvYBLiv0Iok1SgDlKT6MhQ4hRycGoAtgNOAOQXWJKnmGKAk1acNyJfyvgJ8ixykHim0Ikk1xAAlqX4tQR59uhOYRb6k911gXoE1SaoJFQWoiDgmIp6KiFkRcV9EbLuY/Q+OiAcj4u2ImBERl0TE6MpKlqQ+tin5kt4JwJeAbYEnCq1IUpUrO0BFxIHAecAZ5F6/twPTImL1Lvbfhrzc58XkxRUOACaRW9xJUnVYkjz6dDvwMrAR8H1gfpFFSapWlYxAnQRcnFK6KKX0WErpROBZ4Ogu9t8CeDqldH5K6amU0l+AH5F/55Ok6rI18CBwOPA5cgPOZwqtSFIVKitARcRw8iyBGzu8dCOwVReH3QGsFhF7RLYSsD/dLPPZ2tpKS0vLQg9JGjAjgR+QG3D+k7wUzMX0oGudpHpR7gjUCuSbgF/qsP0lYOXODkgp3QEcDFwBzAZeBN4gL6rQqalTp9LY2LjgMW7cuDLLlKQ+sCN5KZgDgE+TFyd+odCKJFWJSu/C6/h7WHSyLb8QMZ68iMLp5NGr3YA1gQu7evMpU6bQ3Ny84PHss89WWKYk9VIjefTpWuB+cvuDX+BolFTnhpW5/yvkG3w7jjaNYdFRqTZTgOkppe+V/vxQRMwEbo+Ir3a2lExDQwMNDQ1lliZJ/WhPcp+o48hj6r8BfgisWGRRkopS1ghUSmk2+WbfnTu8tDNdL8+5FIvex9LWZaWQ9WskqSKjgV8CVwK3kEejri6wHkmFqeQS3jnApyPi8IhYPyLOBVandEkuIqZGxGXt9r8W2Dcijo6ItSJia/IlvbtTSs4mkFR7DiCPRm0B7AMcQp7ZKalulHsJj5TSFaUmmKcCY8n/jOyRUmq70XcsOVC17X9pRCxDHvg+m/zPzM3Af/WudEkq0Mrk0afLyO0ObibPldq1wJokDZjobA7SAFtsAS0tLTQ2NtLc3MyoUaMGoiZJ6rlngSOAPwKfBb4HLFNoRVI9KWQ6kGvhSVJvjQP+QJ5U/jNyF/NbC61IUj8zQElSXwjgKHIX89WA7YHPA+8UWZSk/mKAkqS+tDbwZ+As8ojUROCuQiuS1A+qOkA1NTUxfvx4Jk2aVHQpktRzQ8mrhv4NGEVe6Oor5LUYJA0KTiKXpP40F/gO8A3yWgxXkudMSeorTiKXpEFnGHn06S/A88AHyIsUS6ppBihJGgibkdfS+wCwC/AtFl2jQVLNMEBJ0kBZAbie3Ib4VODDwGuFViSpQgYoSRpIQ4HTyEHqTvKI1L1FFiSpEgYoSSrCbuRLemOArYEf04NbaiRVCwOUJBXlPcDt5GVgPgscBrxdaEWSesgAJUlFagAuAP6X3OJgS+CJQiuS1AMGKEmqBp8A7gZmAZsCvy22HEndq+oAZSdySXVlA+AeYGdgX+AUciNOSVXHTuSSVG0ScC45QG0NXA6MLbQiqZrZiVySRP5xcBJ5UeInyK0Obiu0IkkdGKAkqVptS251sB6wA3AWtjqQqoQBSpKq2crAH4GTgS8C+wHNhVYkCQOUJFW/YcC3gauBm8l36T1UZEGSDFCSVCs+AtwHjAS2AC4rthypnhmgJKmWrA38FfgY8ClyB/NZhVYk1SUDlCTVmhHAT4CLgJ8C2wBPF1mQVH8MUJJUq44A7gBeJbc6uL7YcqR6UtUByk7kkrQYHyC3Otga+BBwKjCv0IqkumAnckkaDOaT79T7GrAj8HNgxUIrkgaKncglSRUaAnwZuBF4gDwydWeRBUmDmwFKkgaTHcmX9MYB2wE/wO7lUj8wQEnSYLMacAtwDHA88HHgrSILkgYfA5QkDUbDgfOAK4DrgM2Ax4osSBpcDFCSNJh9FLi79PUk4MoCa5EGEQOUJA1265ND1F7AgcCJwOwiC5JqnwFKkurB0uTWBt8HLgAmA88VWZBU2wxQklQvAjgOuA14ltzq4KZCK5JqlgFKkurNFuRWBxsDuwBnkhtxSuqxqg5QLuUiSf1kRWAa8FXgK8BHgNcLrUiqKS7lIkn17nrgE8CywFXkS3tS7XApF0lSAfYgX9JbHtgKuAi7l0uLYYCSJMEawF+AQ4EjgSOAdwqsR6pyBihJUrYkcCHwU+ByYEvgn4VWJFUtA5QkaWGHAHcCM4FNgd8VW45UjSoKUBFxTEQ8FRGzIuK+iNi2m30vjYjUyePRysuWJPWr9wP3AjsAewP/BcwtsiCpupQdoCLiQPISlWcAE4HbgWkRsXoXh5wAjG33GAe8BvyqgnolSQOlEfg1cBZwNrAz8GKhFUlVo+w2BhFxF3B/SunodtseA65OKU3pwfF7A78B1kwpPYNtDCSp+t1GXkcvyAsSb1NsOVI71d/GICKGA5sAN3Z46Ubyza89cQTwp1J46lRraystLS0LPSRJBdqO3OrgfeR19H5YaDVS4cq9hLcCMBR4qcP2l4CVF3dwRIwFdid3GenS1KlTaWxsXPAYN25cmWVKkvrcWPLaeccCxwAnAfMKrUgqTKV34XW87BadbOvMocAbwNXd7TRlyhSam5sXPJ599tlKapQk9bVhwH8D3y8970u+W0+qM8PK3P8V8u8bHUebxrDoqNRCIiKAw4H/TSnN7m7fhoYGGhoayixNkjRgjgPWBD5Gvrx3LbBKoRVJA6qsEahS8LmPfC9GezsDdyzm8A8C7wUuLuczJUlV6kPk7uUvAZsDDxZbjjSQKrmEdw7w6Yg4PCLWj4hzgdXJ/WuJiKkRcVknxx0B3JVSeqTyciVJVWUj4C5gRfKdedcXW440UMoOUCmlK4ATgVOBB8iDt3u0u6tuLDlQLRARjcB+OPokSYPPquQ2B9sDHwaaii1HGghl94HqB/aBkqTBYB5wMrnV8gnk5ptDiyxIdaKQPlDlTiKXJKlzQ4Fzyb2ijgeeBH4BLF1kUVL/cDFhSVLfOga4DvgzeZLH88WWI/UHA5Qkqe/tDkwnN7/ZnDxjVhpEqjpANTU1MX78eCZNmlR0KZKkcr2ffIfeyuQ79H5fbDlSX3ISuSSpf80EPgFcQ55gfnyh1Wjwqf7FhCVJKttI4Crg88DnSg/X0FON8y48SVL/GwqcRV6P4jjgKeCXeIeeapYjUJKkgXMUeS7UrcC2wHPFliNVygAlSRpYu5Lv0HuNfIfe34otR6qEAUqSNPA2JN+htwp5JOraYsuRymWAkiQVY2XypbxdgL2B8wutRiqLAUqSVJylyHfonUReP+94YG6hFUk94l14kqRiDQG+R15D7xjyGnqXA8sUWZTUvaoegbITuSTVkc8A1wN/wTv0VPXsRC5Jqi6PAHsCs8mTyzcpthxVPTuRS5LEBsCdwGrAdsDvii1H6owBSpJUfVYGbgF2B/YBzqUH1yukgWOAkiRVp6WAK4Evku/SOxbv0FPV8C48SVL1GgJ8h7yG3tHkNfSuAJwOq4I5AiVJqn5HAjcAdwDbAP8uthzJACVJqg07kQNUC3kNvXuLLUf1zQAlSaodE8hr6L2HfIfe1YVWozpmgJIk1ZaVgD8DHwL2Bc7GO/Q04Ko6QNmJXJLUqRHkyeT/BZxMnmDuHXoaQHYilyTVtouBo4AdyW0P/DFRb+xELklS2Y4g36F3J7A18Eyx5ag+GKAkSbVvR/Idem+R79C7p9hyNPgZoCRJg8N48h16awIfBH5TbDka3AxQkqTBYwxwM7AXsD9wFt6hp37hUi6SpMFlBPALYG3yOnpPAD8AliiyKA02BihJ0uAzBDiDvIbeZ8hr6P0KaCyyKA0mXsKTJA1ehwE3kieVbw08XWg1GkQMUJKkwW174K/AO+Q79O4uthwNDlUdoOxELknqE+uR+0S9l3yH3q+LLUe1z07kkqT6MYt8We9y4DvkSeaF9LFWHyrkDDqJXJJUP5YEfk4eifov8h16F+AdeiqbAUqSVF+GAN8kh6gjgRfJa+iNKLIo1ZqqngMlSVK/+RRwLXATsDvQUmw5qi0GKElS/doV+CPwAHk9vVcKrUY1xAAlSapvWwO3AP8GtgOeL7Qa1YiKAlREHBMRT0XErIi4LyK2Xcz+DRFxRkQ8ExGtEfGviDi8spIlSepjGwO3A28B2wD/LLQa1YCyA1REHAicR26SP5H8v9y0iFi9m8OuJA+OHgGsCxwE/L3cz5Ykqd+sA0wHGsgh6qFiy1F1K7sPVETcBdyfUjq63bbHgKtTSlM62X83cseNtVJKr3XylvaBkiRVj5eB3cjr510PbFlsOVqsQvpAlTUCFRHDgU3IKwu1dyOwVReH7QXcC5wSEc9HxOMRcVZEdHnDaGtrKy0tLQs9JEkaEGOAPwMbADuRJ5lLHZR7CW8FYCjwUoftLwErd3HMWuTB0A2AfYATgf2Bpq4+ZOrUqTQ2Ni54jBs3rswyJUnqhUbgD+RlX/YEflNsOao+ld6F1/GyW3Syrf1nJODglNLdKaXrgZOAQ7sahZoyZQrNzc0LHs8++2yFZUqSVKGlgKvJv/ofAFxaZDGqNuV2In8FmMeio01jWHRUqs0M4PmUUnO7bY+RQ9dqnR3Q0NBAQ0NDmaVJktTHhpOXfmkkr6HXDJxQaEWqEmWNQKWUZgP3ATt3eGln4I4uDpsOrBIRS7fbtg4wH3iunM+XJGnADQUuJK+ddyJwGj24/UmDXSWX8M4BPh0Rh0fE+hFxLrA6+X8vImJqRFzWbv9fAK8Cl0TE+IjYDvge8JOU0ju9rF+SpP4XwLeBqcA3yEFqfpEFqWhlLyacUroiIkYDpwJjgUeAPVJKz5R2GUsOVG37vxUROwPfJ9+N9yq5L9RXe1m7JEkD60vAssAx5Mt5F1HBT1INBmX3geoH9oGSJNWWXwKHkO/Q+yWwZLHl1Lnq7wMlSZLI62lcDdxADlFvFVqNCmCAkiSpEh8iB6i7yQ03O1trQ4OWAUqSpEp9kNy1/J+lr2cUW44GTlUHqKamJsaPH8+kSZOKLkWSpM5tAtwOvE5ed+OpYsvRwHASuSRJfeFpclfEt8nr540vtJp64iRySZJq1hrkkagVgO2AewqtRv3MACVJUl9ZGbiFvN7GDqWvNSgZoCRJ6kvLkS/hbQXsBlxbbDnqHwYoSZL62kjgGnKPqH2AnxVbjvqeAUqSpP7QAFxO7lj+SaCp2HLUt1zBR5Kk/jKMvF7essBxwBvAlynovjH1JQOUJEn9aQhwNnlu1FfJ/aK+hyGqxhmgJEnqbwF8jTwS9TnySNSPgKHFlaTeqeoA1dTURFNTE/PmzSu6FEmSeu94oBE4HGgmTy5vKLQiVchO5JIkDbTfAR8Ftgd+Tb5rT5WyE7kkSXXhI8A04C/AruRLeqopBihJkoqwA3Az8Bh5JOqlYstReQxQkiQVZTPgVnJ42hb4d7HlqOcMUJIkFWkD8iLEc4Gtgb8XW456xgAlSVLR1ibPhxpFHom6v9hytHgGKEmSqsEqwG3AmuQ5UbcXW466Z4CSJKlajAZuAjYBdgGuL7Ycdc0AJUlSNVmGHJx2Jbc7uKLYctS5qg5QTU1NjB8/nkmTJhVdiiRJA2dJ4CrgoNLjx8WWo0XZiVySpGo1HzgB+AHwHeCUYsupUoV0Iq/qtfAkSaprQ4DzgeWA/wJeB86koMig9gxQkiRVswBOB5YFvkBe9qWJKp+EM/gZoCRJqgUnkUPUkUAz8FNgiSILqm8GKEmSasXh5GabHwdagF8BIwqtqG45AChJUi3ZH7gW+DOwGzlIacAZoCRJqjW7An8EHiR3Lf9PseXUIwOUJEm1aCvgVuA5YLvSswaMAUqSpFq1EXkR4reBbYB/FVtOPanqAGUnckmSFuN9wHRgOPly3pPFllMv7EQuSdJg8DwwGZhNvrS3RpHFDKhC2opW9QiUJEnqoVXJd+YtQQ5SzxRazaBngJIkabBYjRyihpIv5/272HIGMwOUJEmDyThyiIIcop4tsJZBzAAlSdJgszo5RM0nhyhbHPS5igJURBwTEU9FxKyIuC8itu1m38kRkTp5rFd52ZIkqVvvIYeoOcAO5Enm6jNlB6iIOBA4DzgDmAjcDkyLiNUXc+i6wNh2jyfK/WxJklSGNcghahY5RL1QaDWDSiUjUCcBF6eULkopPZZSOpF8hfXoxRz3ckrpxXaPeRV8tiRJKsda5BD1NjlEzSi2nMGirAAVEcOBTYAbO7x0I7mpfHf+FhEzIuKmiNi+ux1bW1tpaWlZ6CFJkiq0NjlEvUkOUS8VW85gUO4I1ArkmyM7/qd/CVi5i2NmAJ8B9gP2Bf4B3BQR23X1IVOnTqWxsXHBY9y4cWWWKUmSFvJe4BagGUNUHyirE3lErEKehrZVSumv7bZ/BfhkSqlHE8Mj4logpZT2opNO5K2trbS2ti74c0tLC+PGjbMTuSRJvfUPcqPN0cDNwJhCq+kLNdGJ/BVgHouONo2hvCx7J3n1nk41NDQwatSohR6SJKkPrEu+nPcKsBPwn2LLqVVlBaiU0mzgPmDnDi/tDNxRxltNxGlskiQVYz1yiHqJHKJeKbacWjSsgmPOAf43Iu4F/kqe37Q6cCFAREwFVk0pHVL684nA08Cj5LWiP0GeD7VfL2uXJEmVWp8coiaTQ9RN5Mt66pGyA1RK6YqIGA2cSu7n9AiwR0qpbdnCseRA1WY4cBZ5mcN3yEHqQyml63tTuCRJ6qXx5HlQ25OvJf0JWL7QimpGWZPI+8liC2hpaaGxsdFJ5JIk9YeHySHqPeQQtVyx5ZSpJiaRS5KkwWZD8iW8p8kjUW8UWUxtMEBJkiTYiByingR2IfeLUpcMUJIkKduYfAnvCWBXwIVAulTVAaqpqYnx48czadKkokuRJKk+fIAcov4B7IYhqgtOIpckSYu6hzwfagJwA7BMseV0w0nkkiSpSkwCbqTUrAh4q9hyqo0BSpIkdW4z4A/AgxiiOjBASZKkrm1BDlF/A/YEZhZbTrUwQEmSpO5tSZ4HdS/wYeDtYsupBgYoSZK0eFsD04C7MURhgJIkST21LXA9cCfwEfIKt3XKACVJknpuO+D3wHRgb2BWodUUxgAlSZLKMxm4Drgd2Ie6DFEGKEmSVL4dgGuBW4B9gdZCqxlwVR2gXMpFkqQqtiNwDXAzsB91FaJcykWSJPXOH8iTyncBrgKGD+inu5SLJEmqQbsCV5OD1EeB2YVWMyAMUJIkqfd2A35L7hX1MWBOseX0NwOUJEnqG3sAvybfoXcQgzpEGaAkSVLf2ZM8D+oa4OMM2hBlgJIkSX1rL+BK8ryoTwBzC62mXxigJElS39sbuIJ8Se+TDLoQZYCSJEn9Y1/gcuBXwKeAecWW05cMUJIkqf/sD/ySPBp1GIMmRA0ruoDuNDU10dTUxLx5g+S/tiRJ9egAYD55UnkAPwGGFlpRr9mJXJIkDYxfkOdDHQJcTF9dByukE3lVj0BJkqRB5OPkYZNDyCNQP6ZmJxMZoCRJ0sA5mHw571PksaMfUZMhygAlSZIG1ifJIeow8kjUBdRciDJASZKkgdfW1uAIcnhqoqDZTJUxQEmSpGIcTh6JOpIcor5PzYQoA5QkSSrOp8kh6rPkEPXf1ESIMkBJkqRifYYcoo4mh6hzqfoQZYCSJEnFO4ocoo4lh6izqeoQVdVz3puamhg/fjyTJk0quhRJktTfjiHPgzoXuKHgWhbDTuSSJKm6/BXYgp6OQNmJXJIkiS2LLmDxqvoSniRJUjUyQEmSJJXJACVJklSmaphEvlgRMQpoBhpTSi1F1yNJkgaHiNgCaEgp3VrWcTUSoAJYBngz1ULBkiSpJkTEi8CtKaUDyzmuJu7CK4UmR54kSVKfiYj1gZWAO8o91jlQkiSpXm1Ter6t3AMNUJIkqV5tS55j/WC5BxqgJElS3YiIj0XEvRFxL/AJYA5wd2nb5j1+H+dkS5KkehMR44B/A19KKX2n3OMdgZIkSfVo+9JzWe0L2higJElSPZoMzATureRgA5QkSapHk4HpKaW5lRxsgJIkSXWlNP9pTSq8fAcGKEmSVH8Wmf8UEZ8vrXzSIwYoSZJUbyYB84B7ACJiLWCNcpaLK3Qpl3Zr3EmSJFWq3LVyXwPeSCnNjohRwFTgqHI+sNA+UKWimwsrQJIkDQaNKaUer5kbEcsDVwCvkhtpfjOl9Hg5H1h0gOrpCNQywHPAasCb/VjS3cBm/fj+A/k5A/EZnpfq/AzPS3V+huelOj/D81Kdn1HueSl3BKrXCr2EV/pmF5sY283perOchFmuiJjfn+8/kJ8zQJ/R9qXnpbo+o+1Lz0t1fUbbl56X6vqMti89L9X1GW1f9ut56Q0nkS+saRB9zkB9LwPB81KdPC/VyfNSnTwvg0xNrIXXbq5UWdc41b88L9XJ81KdPC/VyfNSnWrhvNTKCFQr8I3Ss6qH56U6eV6qk+elOnleqlPVn5eaGIGSJEmqJrUyAiVJklQ1DFCSJEllMkBJkiSVyQAlSZJUJgOUJElSmaomQEXEMRHxVETMioj7ImLbxez/wdJ+syLiyYgoaxFA9Uw55yUi9o2IP0bEfyKiJSL+GhG7DmS99aLcvy/tjts6IuZGxAP9XGLdqeDfsIaIOCMinomI1oj4V0QcPlD11osKzsvBEfFgRLwdETMi4pKIGD1Q9daDiNguIq6NiBciIkXE3j04pup+5ldFgIqIA4HzgDOAicDtwLSIWL2L/dcEri/tNxE4Ezg/IvYbkILrRLnnBdgO+COwB7AJ8Gfg2oiY2P/V1o8KzkvbcY3AZcBN/V1jvanwnFwJ7AgcAawLHAT8vX8rrS8V/GzZhvx35GJgAnAAMAm4aCDqrSMjgQeB43qyc7X+zK+KPlARcRdwf0rp6HbbHgOuTilN6WT/7wB7pZTWb7ftQmCjlNKWA1FzPSj3vHTxHo8CV6SUTu+nMutOpeclIi4HngDmAXunlDbu71rrRQX/hu0GXA6slVJ6beAqrS8VnJeTgaNTSmu323Y8cEpKadxA1FxvIiIB+6SUru5mn6r8mV/4CFREDCePVtzY4aUbga26OGzLTvb/A7BpRCzRtxXWpwrPS8f3GEJeUdsfEH2k0vMSEYcBa5M7+6oPVXhO9gLuBU6JiOcj4vGIOCsiRvRjqXWlwvNyB7BaROwR2UrA/sDv+69S9UBV/swfVtQHt7MCMBR4qcP2l4CVuzhm5S72H1Z6vxl9WWCdquS8dPQF8lDtlX1YV70r+7xExPuAbwPbppTmtlvlXH2jkr8rawHbALOAfUrvcQGwPOA8qL5R9nlJKd0REQcDVwBLkn+mXAMc3491avGq8md+4SNQ7XS8lhidbFvc/p1tV++Ue17yThEHAacBB6aUXu6Huupdj85LRAwFfgF8PaX0+EAUVsfK+bsypPTawSmlu1NK1wMnAYc6CtXnenxeImI8cD5wOnn0ajdgTeDC/ixQPVJ1P/OrYQTqFfKcjI6/EYxh0cTZ5sUu9p8LvNqn1dWvSs4LsGDi5sXAASmlP/VPeXWr3POyDLApMDEiflDaNgSIiJgL7JJSurm/iq0TlfxdmQE8n1JqbrftMfIPhdXIc9XUO5WclynA9JTS90p/figiZgK3R8RXU0pe3ShGVf7ML3wEKqU0G7gP2LnDSzuTr0d35q+d7L8LcG9KaU7fVlifKjwvbSNPlwIfTyk5b6CPVXBeWoANgY3bPS4E/lH6+q5+KbSOVPh3ZTqwSkQs3W7bOsB84Lk+L7IOVXheliKfg/bmlZ699l2c6vyZn1Iq/AEcCMwmX/tfHzgXeAt4T+n1qcBl7fZfE5gJnFPa//DS8fsV/b0MpkcF5+UgYA5wDPm3hbZHY9Hfy2B6lHteOjn+NOCBor+PwfSo4O/K0sCzwK+A8eQWII8D/1P09zKYHhWcl0NL/4YdTZ6ntjVwD3BX0d/LYHqU/v/fuPRIwOdLX6/exXmpyp/5hf+HbPcf6BjgaaCV/FvDdu1euxS4pcP+HwTuL+3/FHBU0d/DYHyUc16AW0p/GTo+Li36+xhsj3L/vnQ41gBVBecEWI/cN+3tUpg6GxhR9Pcx2B4VnJfjgUdL5+UF4GfAqkV/H4PpAUzu7mdFrfzMr4o+UJIkSbWk8DlQkiRJtcYAJUmSVCYDlCRJUpkMUJIkSWUyQEmSJJXJACVJklQmA5QkSVKZDFCSJEllMkBJkiSVyQAlSZJUJgOUJElSmf4frhu4xEX5zOIAAAAASUVORK5CYII=\n", "text/plain": [ "Graphics object consisting of 1 graphics primitive" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "#Plot first component of solution versus time, as simple example.\n", "yp = yv4.column(0);\n", "p = line(list(zip(tv4,yp)),rgbcolor=[1,0,1],legend_label='Implicit Euler')\n", "p.set_legend_options(loc='upper right');\n", "p.axes_labels(['$t$','$y0(t)$'])\n", "show(p)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "SageMath 9.2", "language": "sage", "name": "sagemath" }, "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.7" } }, "nbformat": 4, "nbformat_minor": 4 }