Skip to content

Instantly share code, notes, and snippets.

@vuddameri
Created November 6, 2021 07:37
Show Gist options
  • Select an option

  • Save vuddameri/75212bfab7d98a9c75861243a9f8f272 to your computer and use it in GitHub Desktop.

Select an option

Save vuddameri/75212bfab7d98a9c75861243a9f8f272 to your computer and use it in GitHub Desktop.
Thomas Algorithm.ipynb
Display the source blob
Display the rendered blob
Raw
{
"nbformat": 4,
"nbformat_minor": 5,
"metadata": {
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.8.10"
},
"colab": {
"name": "Thomas Algorithm.ipynb",
"provenance": [],
"include_colab_link": true
}
},
"cells": [
{
"cell_type": "markdown",
"metadata": {
"id": "view-in-github",
"colab_type": "text"
},
"source": [
"<a href=\"https://colab.research.google.com/gist/vuddameri/75212bfab7d98a9c75861243a9f8f272/thomas-algorithm.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "e1314e17"
},
"source": [
"<h1><font color = 'blue'><center>Tridiagonal Matrix Algorithm (Thomas Algorithm)</center></font></h1>"
],
"id": "e1314e17"
},
{
"cell_type": "markdown",
"metadata": {
"id": "d50d8c69"
},
"source": [
"<p > Sparse matrices are commonly encountered when finite difference and finite element methods are used to solve differential equations. <font color='red'>A sparse matrix is a matrix where a majority of elements have zero values.</font> The presence of a large number of zeros creates both challenges and opportunities when solving matrices. The primary challenge arises when inverting a matrix using traditional methods. However, for sparsity also helps reduce the data storage requirements and lead to simpler approaches to solve systems of equations involving sparse matrices.</p> "
],
"id": "d50d8c69"
},
{
"cell_type": "markdown",
"metadata": {
"id": "1022efea"
},
"source": [
"<p> <font color='red'>A tridiagonal matrix only contains non-zero elements along the principal diagonal as well as along the left and right columns immediately adjecent to the diagonal. </font> If the coefficient matrix is a tridiagonal matrix, then Gaussian elimination can be readily applied in an efficient way. <font color= 'blue'> The tridiagonal matrix algorithm (TMDA) or the Thomas Algorithm </font> is an efficient application of Gaussian-elimination to solving system of equations whose diagonal terms are sparse. "
],
"id": "1022efea"
},
{
"cell_type": "markdown",
"metadata": {
"id": "ea0eec40"
},
"source": [
"Consider a system of linear equations represented using the matrix form below:\n",
"\n",
"$\\begin{equation}\n",
"[A] [x] = [d]\n",
"\\end{equation}$\n",
"\n",
"If A is a tridiagonal matrix then it will be of the following form:"
],
"id": "ea0eec40"
},
{
"cell_type": "markdown",
"metadata": {
"id": "da0d2a40"
},
"source": [
"![tridiag1.png](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAUEAAACOCAIAAAAgpS2zAAAAA3NCSVQICAjb4U/gAAAAGHRFWHRTb2Z0d2FyZQBtYXRlLXNjcmVlbnNob3TIlvBKAAAY3UlEQVR4nO2dwY/bxhWHh2sfCsMOuUBqwICbpdykh6ZdUS7apBeL3EOPu1R6ycmiAvTQtLCof8Cici1gUbf2YnF9bWNRBdpbwmEuBVLAHLmHoj2EVBGgBzfWMNi2F3fZw2sIdXe9q12Rkrj7vsOCoqSZ4Yo/zsybN+8JSZIQBEEKy9qyG4AgyFyghhGk2KCGEaTYoIYRpNighhGk2KCGEaTYoIYRpNighhGk2Fye8XOHXUEEQci6MUghOXBv4I2xYGbthyuVytoU3/nOd3JtFlIUkiR5/fXXp++Nt99+e9mNuljM2g8TQt5+++0f/ehHcHz9+vUD7yZJQimtVCqSJGXWOmTlEQTh3r17z58/h5e///3vl9ueC8gpNPzDH/6w0+kc+ZbjOJZljcfjyWSSUcOQhRIEwe7ubhiGmqaZpnmq7zabzfT42bNnT548ybp1yHFkY9MyDKNer5fLZeyEiwhjrFarmaY5HA7DMDQMY9ktQk5BZnZp3/dVVc2qNGSRtNttXddlWSaEdDqd3d3dIAiW3ShkVrLRcJIkvu9Xq9UkSXLdzJh3+RcT3/fTARQcUEqX2SDkNGSjYfjJd3d3LcvSNC2nOwC6C8uyFEXhnOdRxQUkSZI4jg88GfHfWyBOYdM6BkqpKIqO40iS1O/3G41GGIbwFuc8k0myqqqyLA+Hw8Fg8MEHHwRBoGna/MUi8MA9sKjr+35+NX766ad7e3vpy2vXrn3/+9+fs8wXL1588skn02du3rz5rW99a85ij+Ef//jH06dPp89sbm6++uqr85T597///c9//vP0me9973uiKJ7wtWQ2yuVyq9V62bvVavXevXtw3G63odh+v28YRjrAnof79++Lovj8+fMkSSaTycOHD5Mk8TzPNE1JkuYv/yLz2WefEULu37+fniGE3L1792yl/fSnP33rrbeO/8ybb7554NY/W13THF4Qef/99+cv9hg+/PDDAzX++te/nrPMX/7ylwfK9DzvxG9l0w/7vl+v1+E4iqLXXnstSRIwk1iWNWfhSZIMh8Nqtbq+vk4IkSSp0WgQQhRFUVXVtu05y7/ggClrlpMZ8u677/7kJz+B42vXrs1f4NWrVz/66KP05c9//vO83cWg/N/97ndf+9rX4Mzm5ubLPswYo5TGcQw93MvY3t5+44034PjZs2fvvvvu2trJs90MNBxFESEkHdlSSnVdFwQhq3UmQRBGo9HOzs6B87iOlQmCIGxsbBw4mbeGv/GNb2xtbWVY4OXLl6cLvHbtWpKz7RPKv3PnztWrV0/8MGPMtm1Zlo/X8I0bN27cuAHH4/GYELK/v39i4RnYtOAOgFF7v98nhLzMFWTOKlL6/T48OJBMqNfro9EIjmFVSdf1pbbovAFL7tVqNY/CM+iHBUG4f/9+o9FQFIVSGgRB5j2k67q6rpdKpSiKwjBMFzORTGi1Wqqq9no9sHo8fPgQpi1IVnDOx+NxTg4U2cyHO50OSHf+2e+RKIoSRZHneeVyWVEU3BmTLZIkMcYGg0EQBI8fPy6VSsttTxRFjDFFUWRZ9jxPkqRKpbLcJp0NzjljjBACJreV1jAh5Mj/crbLjIcXk3AZM0Nqtdqym0AIIa7rep5Xq9VgtKWq6u7u7oMHDwq3lMgYMwyj3W5TSh89elStVnPqe3KMAeA4DmNMVdWcOmdKqW3b7XbbsiycHp8POOeO4/R6PVVVJUna2NhQFAW6spRC+JDB7pEHDx7UarVer5ckyfRkOO2fMyGzfvgwebvOq6qKHtrnD1gsTJLE9/1ms6mqamphdhzHdV3O+erL2LIsURRh7DCZTOI4hns1iiLHcRzHMQxDUZRM6sJYPMgKIUkSWCtBpQee0YZhnHZf5FJIPRpg8AzXAooFj4lsLbKoYWSF4JxblsU5Hw6H5XIZzONFdOOJ4zgVKqUUriWnSSVqGFkhut1up9PhnLuuCyuUQRBkOHVcDIIgVKvVtP2PHj1SVXUwGBz2pcmEHOfDCHJaarWa7/umafb7/d3dXVVVy+VyEfth27YtyxqNRqIo9vt9MLu6rptHXahhZIUANyE4Ltxi0jSKokwrNle/NxxLI4WBMeb7fhzHvu8X1zUANj+Mx+Os5gjYDyNFolqtwm7WZTdkLrrdboalFUPDsCYeRRE47iy7OchyyGpBdblkfgMXYywN7juNRqPoD2AEyZxiaFiWZdg/jJ0wghygGBomhFBK8/MaR5DiMpOGkyT55z//+cUXXxz/sSAIcopLnHwV+9bzvMX4ykJFxTV+Losvvvhib28PpzyLZCYNv3jxQhCEY36YIAhgi6/jOJqmZb7bgXM+Go2GwyHsIs51l5zneXAtjLFCL1EuHrhDBEH4z3/+s+y2XCBmtUtfuXLlZTEywzDc2tr68MMPt7a2arWaIAgPHjyAt0B7oijOaVGEvtdxHFBXqVQKgqBSqTDG4jjOMEcMPCCCICiVSrquQ4RduAqSWyyV84QoileuXFl2Ky4WGcyHp7dZwT5esDyB/3qSJJZlzbndBCbD8CBIdQX+a2EYZhgyvtFo7OzsQCAL27bhcuAq+v3+iuySR5Bp5tUwbLPa2dkBa5PneWmv6zgOuLzatt3r9eaRme/75XIZjsG7RVGUTqcjSZJhGLIsZ+JSGwTBeDxO3eJkWZYkiVIKVwGbV1d/5yqygsDKaE4elxn0w3Ecp0NlEFscx7DL2fO89GPzDHc55+mq0nA43N7eXl9fh0n4/IVP10IOORLAjpP0JQbERc6Aoij5paGbV8MQnRjubMbYcDjUNG0wGCRJIkkS6MGyrOPD6p5IOlr2PI8x1uv14CQEc+OcZ2JFK5VKoihCXF9CSJrFEx4ftm3X6/Xz4SqELBiw46xubFrXdU3TpJTKstztdiEWdrppAwxRc86HLcsyDGM0GsEEOI2cyjm3bZtSmkn3KMsybBOrVCqTyWRjYyMdojPGGGOO48xfC3IBCcMwDceTORloGPaLJUkCU+LpLhFmj6ZpOo6j6/qZlQaB0dIqAM65aZq2bUuSZNt2JlFaarXaYcMVqNdxnCMDxCCnYn9///PPP3/y5Mkxn3n8+PGf/vSnzc1NTdPg9z3Vj7u/vx+G4eXL+e4F+Mtf/kIIOX4VLQzD3d1dQogoiqddnXn27Bkh5LPPPjvxfsvsOg97UDHGGo2GLMu7u7tRFM0/3D1QhWEY4/G4VqtxzpvN5pyFvwzOuaZpiqJomgar0zlVdEHY39+/efPmMamJ2u025DQBs6XjOLVaTRRFSLI1C0mSlEql7373uwce+tly69YtcqyG+/1+r9ejlH788cc//vGPT+tlCEu5r7322omfzPFZBXlM0+PMy59+Nuc6TZ22aWF+iVwJw3A8Hqe5fsBJfnrlstfrSZJUr9dX/IfwPO+999578uSJJEnvvPOOKIrpZNh1Xd/3ZVmu1+uZ6CJHDcuynOs/ejFjWkmScPC8MEqlEhgdwK+m2+1Ox5c2DINSCikgGGOrLONer7e5uZl6NMRxDI8haH+32zVNs9FoTHcPZ6Ywex6QCwXMWaafnhDEI83hssoL9eA0oarqdGxa6IchiSkhRNf1rMJroYaRFQJ81CF83J07d0AD9Xo9SRJFUeAv5B9b/UW+NG0VbNeJ49g0TVjBIYRwzrMKc4kaRlaIwWAAvrS+74OAHz9+LMvytDXItu1ms7nKGhYE4e7du+Co2+/3wVN4MBikObRhTRT7YeQc0mq1yuVyt9tljO3s7Oi6HsfxdDprx3EgVcKKZ9iCHEtgtXJddzQajcfjdBucZVmwZpZJXcWIp4VcECRJGg6HcNxqtVqt1vS7MAo1DAOc81bZpiVJEngTAtNdLjgsgeNTJn7+2A8jxSCKolqt1mg0BEGoVCoFdVx3HKfT6VQqFUEQsjLLYT+MFANZliETd6ExDCPzCBnYDyNIsUENI0ixQQ0jSLFBDSNIsUENI0ixQQ0jSLFBDSNIsUENI0ixQQ0jSO6YplmpVHJyDkUNI0ju2LYtimJOwSRQw8uh2Wym+0uRiwBjLKfYtKjhpfGy/FXI+SMIgpWOTXuesCzL930ItZtrRdMb05ADpNEODcOwbTtJEgiFt9xWnQEIAEC+SqWQ03wYNfx/QEqKc7A/prgYhtFsNiuViqqquq57nnf79m1VVWePTbsiGIYxmUyGw2EYhrdu3UpzkmUOavgglNIlBrLMNSTy6hOGYblcrlQqhJAoiprNpiAI9XodAvdTSqMoiqJIluXMd/BlS7/f393dhbhCsixPx6ZNE6GYpplJRCGcD/8fSZJ88skny8oz3O1219bWLnIY+lKpBLE7INC0pmmQ1BIG0mn4i0ajseKxeGzbLpfLMHiGLNnQMVBKe72eruuyLGeSmYQUqB8eDAaUUkiP9vTp05wSO8Dsxff90WgEKVoWOQ2TZblaraK9mhBCKU2Tm0DHS77qwaIoEkVx/rnlYDCA1PBHctoEMdMkSfL06dP79+/DS0joW6lUOOeKosCtm2EsoWJoOE17keZSzE/DoijC877dbjebTciXsxiOzPZ0ofA8b2trKwgC13XTcWaz2XRdVxAEXdcdx4E7Yf66YIibvsx2CiOKYprZ77e//a2qqmEYQiJuXdc7nQ6l9ALF4gnD8L333hsMBtAlpjImOUyQfN+fTlD66NEjx3GGwyHnnDEGVpb5a0FeBmOsXC5zzuM45pxHUWTbtq7rqcDg8QoJH+YcIsGsOw8EQWg2mxBW2nGcKIog+RtMEyCbVBiGrVYrTWY0DwXQMPi4QGzeJEl8359ee3BdV9f19fV1VVXnH5xMFx5FEWSsgkCKsizXarWlaDgIglKpVMTFldOSToYppZ7n9fv91BmGcz4cDuEJG0WR4zhzzieDIICk1kciCMI8ps1OpwMLS6ZpwnGpVJJlmTEWRZGu64ZhaJrW7Xbn/1kLoOHRaKQoCjyJfd8nhKiqCn2vZVkZTpDAmFQul+Gl7/uapgmCYFmWJEmc86V4ZQRBcPv27bt37y5yVL9E0ni0mqalAZkJIYwxwzCq1aosy3Ecz2/RjaIozeR0JHMuT0wvhqUX5bouY0zXdbidziZgz/OmI3sWQMNk6r/pum65XIYpq+u6hmFkO0EiX6XYAPtZt9slhJimaVkWpTSrwPynolQqbW9vF251NHNUVW2327Ztj8fjdrs9//rfUqwPcC91Oh2Y85+hhCAItra2dnZ20q8XQMPNZrPX60GyKcaYIAi9Xq9er8O7GU6QSqUS5NRYX183TdPzvNQsAYH5W61WOi1fGNNh0y84lmUtuwnzAvnu5ymhUqk8ePBgeoRSgPXhWq3WbrchXQ2l9N69e6IoQupwGF4qijIajSDn5Zy4rhtFURiGruuCzYNzDn6RkO3q+NEXgiyAVqs1PZUoQD9M/n9qlI4qGWPNZhMmSCSjNOKSJB3ID+K6Liwygf3jIhiWkGJRDA0fiaqqkLEmqwnSkRiGEUWRZVmMsX6/v8rp9pCLSYE1TBY1QToH0zDkHFOA+TCCIMeAGkaQYoMaRpBigxpGkGKDGkaQYoMaRpBigxpGkGKDGkaQYoMaRladKIqyCnmxRDjnOe17Qw0jq45pmpqmFV3G3W63VqvlEfCw2L6WyEXAsqxqtbrEgMGZAHuV8wgAhBpGVh1FUc7BVpP8rgLH0ghSbFDDCFJsUMMIUmxQw8iq0+12NU1b8eQsJzIYDCqVCmRgyhbUMLLqjMfjIAjyziabNxAK95iI1mcGNfw/BoPBYDBYWHVBEGQSxO8iYNs257zoaaharVaSJLi2lCONRiOO44U97FutFuSFuciZSpFMQA3/j8FgkMc452XYth2GIQoYmR/U8P+YDrq9AM6H38IZWFtb+/zzz58+fZpfFYIghGF46dKlXB+Rf/3rXwkhL168yKl86FFmseThfBhZKGtrazdv3rx9+/bsX+GcB0FwqipKpVJ+aQ2BN998kxBy5cqV2b9yqqu4fv06IeT1118/8ZOo4eNgjFmWtbAxNmPMNM1FDukLQb1ev3379qkEsIJYlpXTVaCGj6Pb7XY6nTzW9I5kMBj0er08trYAnudBLricys+JWq22vb1ddLu0qqo7Ozt5XAVq+Dh6vd6TJ0/yHpWltFqtjz/+OI9kfJ7nQdpextiCZ/7zYxjGcDgsepYcVVVd183jKtCmdRySJC1MwFBdHgLzPK9Wq0EicsMwzjys4JxHUXQxTXGrDPbD559Go5GO4mzbPpvTomVZ0JNk3DhkbmbV8L///e8vv/wy16YUhSiKNE1bWBJTznmtVjuzeIIgGI/Huq7DS0mSzjacsywrLeQYvvzyy3/9619nKB85MzONpS9durS3t1c4W0hOBEFAKYUZ5gKqg2TIr7zyyiwSOgz8ageayjl3HIdzLkkS/DVNc/6mCoLAOd/b27t06dL8pSEzMpOG19bWvv71r3/zm9/MuzWFAOaWC5sWViqVeexqlUpFFMXRaARj6TAMO52OoiiqqqqqyhiTZVkQhEw0TAi5desW5xz9zxYJ2rTOwoLtOvPY1SRJ8jyv1WoNh8MkSTY2NmzbliTJcRzDMGRZZozt7OzAh4+MO1f0QFbnHtTw+adSqRwWJ6UUBueO4+i6TimFnnnxzUPmBO3SFxQQLSGEMSZJ0okmOkrpeDwej8e+7y+ifcjMYD98QXEcBwzUjuNEUTTLfLherxNCir4X//yBGr6gpMNmWZZlWZ7988gZCIIgjuMZ/9WnBcfSCJIvlNJGo6FpWk6O96hhBMkXVVVhqpLTWAY1jCC5QymtVqs5LZujhpGVYzKZTCYTOF7Yxs88SJIEDP6g4ZxqQZsWslo0m01CyGg00nWdMZYkyfr6um3by27XqWm326PRaGNjA9bk8jMKooaRFaLf7+/s7GxtbQ0Gg3feeScMQ8MwxuPxstt1atrt9u7ubhAE6+vruU6GCWoYWSnW19e3trYIIYyxcrksyzJ4hi67XacjDMMPPvjg4cOH6+vrhJDJZHLnzp38fMhz0TC4AaDj+yrDOYcB6kq5WKZ7s4bDIUwgixhyADaKwn81SRJKqWEY+VWXvU2r3+/rum5ZlqIoRY9jdo7hnHPOO51Oaj1aBcIwjKJoMpmMRiPQAOccnL0ZYxDAwHEcyPyw3KYeA+d8Y2MD3DmiKPrb3/6mqipcheM4juPAtWQVUCFjDXue1+l0HMfpdDr9fl/TtJW6RZAUWZaho1udThhSmViWNRwOyVe7tbrdLgzrXNeNosiyLMMwYP/zclt7DJIkwY5OQggMdjRN63a7sNOz0WhAAFPLsjKpLmMN27atqipMAyqVShzHGL1lZaGUlstl+LFWAUEQFEVZX1/3PK9er3c6nVartbGxATHGDMOglMJ9D2JYbmuPodFoEEJarVa9XpdlWRTFVqtlGAZou1wuw9A6q0yOWc6HkyTxfR/WBoBqter7PlwSciRBEMA6yuItN77vV6vVwWAQx/FSGnAYSiljDObAYRhKkpQ+YiDkCEiXMaaqKkQgWWJrX4YkSZTSMAwh7kKj0UhjCTqOAwMf13V1XY+iaP6HUcb9cBzHB84UPW1sfnDODcNwXXdjY8MwDAhks7DaYc7p+74gCJPJZPEj6mSK6fOpEatUKk2PEdLNknDguu6Bf9fLCsyVYypNQ0lLkpRe1PRV6Lp+eEZwhqvIUsNoiD4VzWYTrEqapimKsrGxscheBUws0Bu0Wq3RaJRfbPrD/OIXv1j7ihktz2nQL0hVdWAPEOd8bYpPP/00p5Yf4JVXXkkr/c1vfnPi51VVBdu7YRgwK55+91e/+lVa2uzR5DMeS2dY2vkmCIJHjx49fvwYXqYDyIVBKd3c3AQZLHj54OHDh3t7e+nLa9euzfKtdKRwZAS/q1evfvTRR9Nnbt68OVcrT+LOnTsHatzc3DzxW+ki05FJ87a3t994443pM7OEYcpSw4IgiKI4fYZzXsT1vQUApr50OZRS2m63F9kAmFKmtRNCFmYl+sEPfpB5mZcvXwbnkIXx6quvZl7jjRs3bty4cdpvZTwf3tnZmfaMG41G+bl6F51yuQyzD8ZYHMeqqsI6hG3bEODKtm1YR8mjdkEQUtH6vg8JjSBQnuu6uVaNZEvGGoYFAFgT9jxPFMU8sgedA9L1Q5gSE0IqlQpjDFZNYLOLaZqQISmnBoC50fO8KIp6vd50YK1cq0ayJWNfS03T2u12rVZTFIVSSildTev/0oH1BlgzbLfbSZKkzm2QiiENWJeTxdi2bdM0IVWq53lgAZZlmXMOI/zFT9GRs5G9v3Sj0Wg0GkmSoJn6eKD7BaY9YSilIB4wGuekJQgxfeBkuvLhuq6qqrB6iU/hFecUGv7DH/6Q2l2uX7/+s5/97JgPo4DPTBr5GSQE7sELqxo0fKqqe73e8+fP4fiPf/wjJmpZNMlslMvl6W99+9vfnvGLyGkJgiA99jxvkVWHYXjaqvf392/dujV9b7z11ls5NQ85EiGZbVH38Mewp0UAvDeWy38B21kzcWm6ucAAAAAASUVORK5CYII=)\n"
],
"id": "da0d2a40"
},
{
"cell_type": "markdown",
"metadata": {
"id": "a981f017"
},
"source": [
"Just like Gaussian Elimination, the TMDA proceeds with the <font color='blue'> forward sweep </font> which results in a last equation having just one unknown leading to a solution of $x_n$. This is followed by a <font color='red'> reverse sweep (backward substituion) </font> where the estimated $x_{i+1}$ value is used to compute the unknown $x_i$. "
],
"id": "a981f017"
},
{
"cell_type": "markdown",
"metadata": {
"id": "c39fcbfa"
},
"source": [
"<h3> <font color='blue'> Forward Sweep </font> </h3>"
],
"id": "c39fcbfa"
},
{
"cell_type": "markdown",
"metadata": {
"id": "2a8a5f12"
},
"source": [
"The forward sweep perform row operations wherein the pivoting row is modified and subtracted from other equations to eliminate all but one unknown $x_n$. These operations can be succinctly summarized as follows:"
],
"id": "2a8a5f12"
},
{
"cell_type": "markdown",
"metadata": {
"id": "a832d3b5"
},
"source": [
"![forwardsweep.png](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAUEAAADsCAIAAABlmYPiAAAAA3NCSVQICAjb4U/gAAAAGHRFWHRTb2Z0d2FyZQBtYXRlLXNjcmVlbnNob3TIlvBKAAAgAElEQVR4nO2df4wj6ZnXX8+sVkd2Z1zOBe2GXOLykORymzAuz6K7KHcZV/VeDoKgXT0cSqKIdnmIAgnCPxYhTkHbLg9IBxdolxed4CLiqp7TBQS37fKKPy5k1i4PcEouol1eSEQ2S1c5kJDcSe3XyWT32M24+OPZqXXc3e7q7ir/KD+fP0btdvmtp8f++n3e532e5404jkMQBFlaLszbAARBzsVKaNiyrI2NjXlbgSCBEH4Nq6p65cqV69evz9sQBAmEkGtYVdWbN29ub2+XSqV524IggRAJcUzLsqxr164lk0nDMOZtC4IERZjn4UqlQiktl8vzNgRBAiS087DjOBcuXIhGo4PBIBKJzNscBAmK0M7D4D9zHIcCRsLNQ/M2YIFoNBq6rkciEYZhZFlmGGbeFiHIyaCG3yCbzQ6HQ13X2+322tqaKIo8z8/bKAQ5mdBq+FQutKqqt2/fPjg4IITEYrF8Pg8C1nXdMAxFUQIyEkHOT2g17D1W5ziOoiiZTCYWixFCOI6r1WrwFM/zHMcFZSKC+EFoY1qn4sUXX0wmkxO/tG271+uxLDsPixDEK6Gdh0/lS0ej0XGtqqr65JNP1ut127Ydx8GFMbLIhFbD4Eun0+kTr4xEItVqtVarXblyxXEcXdcJIYIgFItF9KWRxSe0Gn7ppZcIIaZperk4l8ulUql2u80wTLlchoWxrus8z+MOE7LghFbD7373uwkhh1e5x8Fx3MSUq+u6KIqyLMuy7Lt5y4tt25RSf90T27Zt28Y1y9nAmNaxcBxnmqYoivM2ZLGQJEmSJF+Gsm1blmXXCfJlzBUktPMwrIdt2z7zCMVi0TdrQoSmaX4NxTBMOp1mWdbHMVeQ0GoY6Pf78zYhbPi42cYwjCAIYa26mRnh96XxI+IXsLJA92TRCPk8jPgFpVSW5UajceHChSOTyVVVHV+5HN6fnx4axPKyMxN+DU//cHzve98bDAYex3niiSd8Mmr56Ha7hUJhZ2eHEHJkUDqXy53H5UF36cyEX8OO40yR8Re/+MVGo+FlnIsXL37961/3z64lQxAEQkgmk8lms8ftmeNcOhdCq2H38zT9g/XMM88888wzM7Fo6bEs6/nnn2+1Wkc+22g0er3elJfjNntAhFbDrm82fR4+G7qu9/v9QqHg77ALjq7r8Xic5/lsNgtO9Tgsy0aj0eNee9q3oFqtJhIJ3Jz3Qmg1HKhfpyhKp9O5evUqeJgrAqWUZVld14/cXkqlUmcYs91uv/jii4SQu3fvapqWTqcTiQQh5Omnn2YYBjXshdBq2CUIMWuaBqURK0WpVGIYhlJaqVT8GpNSevXqVdc/Hw6H8MPe3p6iKEH4UOEjtBoO1JdmWdZxnFUrh2AYxvdW+8edoZNKpYbDIQrYC+HP8fDyOTAMQxAEnuc9dpOnlHY6nbN5j4gXDMPwXq+y4oR2Hnal62Ue5nnetu1cLhePx70MbpompvgGir8ee7gJrYZdPPpjlmXF43EIqJwIVskFDUazvBNaX/revXuEkLe97W0er+/1eqhMZBkJrYZfeeUVQsgjjzzi5WLHcQzDiMfjuVxOkqTzVCwiyIwJrYYff/xx4tmRNk1zOBz2+31VVTOZzErt+iLLTmg1DHisH4ZJGMJUDMPYtg1tJUzTxAxBZMEJrYZPtbXYbDbdDpiUUkKIu/eLi2RkwQltXPpU5zyYppnJZOChYRjRaBT2flmWXbVEDmTpCPk87LFxDMMwbk1ss9mEY8dlWeZ53mN3WwSZF6HVMMzDXjQciUR4nu90OoQQSZLS6XSpVDJNU5IkSinOw8iCMzcNm6YJrdgDws2e94KiKNFoVJIkt6qO4zjbtjmOw/OWkAVnbuthXdenlJuen263SwjxeIsjs/k1TZMkCf713TwE8Yu5zcOGYQSqDVjHnieqzDCMaZo4DyMLTmRevcgikYhlWQEpxHGcWCzGcZzHOiQEWV5m50t3u91ms0kIEUWRUprJZIKb4hRFGQ6Hq9YrB1lNZqThQqHQ7/er1Wqz2RQEoVwuT3ekVVWdkmLFMMyUTuXdbvfWrVtbW1vH1ZcjSJiYhS8Np/uapgkdkhzHgQl5Ct1ud7zt8+Gkq+NSmimlgiCsr69j9SmyIgSu4cFg8Na3vnV7e9v3Ni5HAqd+K4oyg3shyCIQuC8NezynPa622+1C3vKRMAxzXB8cRVEkSVIUBY8FQlaEGa2Hx/d4bNuGaJYoiqIoHrkwNk1zShHvFA1zHKeqKpymN5uZH0HmjBMwBwcH8Xi83W47jjMajba3t7PZbNA33draIoTATREk3MwipgXHvcP3hSRJEI7SNC24JuDQzRxPl0dWgfnkeGiaxnFcKpUK7u7ZbPb27dt7e3vYQRYJN/PJtRRF0bZtt2Q3CKBDJeZpIaFnPhpmGEbXdSjuC+gW0GE8uPERZEGYW82DruscxwV3ahFULHnspzUDTNNUFEWWZWyaifjL3GoPi8VioKVLkNq1IIIpl8u9Xq9arfZ6vVQqNZ6ChiDnZG4aXp1+kdVqdWdnB1JNS6WSexwMHAJ62uwXBJkgtD3xFuTIvMFgcOvWrXw+Dz19xtcOWJyM+EJoNTyXPbPDmKZJKT2yFcHqeCJIoIS2J95CMZFqSiktFovoRSO+ENp5+Ec/+hFZgL0ljuPi8Xin0+F53nEcTdN6vR7LssViMbiYPLJShFbDo9GIEHL9+vX5mhGLxQzDKBQKhmFYliVJUi6XIw/qq+ZrGxIOQqvhS5cukbEjV+YIy7KHex4oimIYhq7reNAuRC68xyAdD4fCrxShXQ/DJ2PuvvRxQKbaMi6JC4WCj2aXy2VBEHieZxim0WhMv9g0TVEUi8Uiz/MbGxsL++aemW63K8tyLBY7VUQ2tPMw0Ov15m3C0ei6zjDMIrgJpyWRSPh1jlyj0RgOh5DTXq1Wb9y4sbu7O6ULGqgX4vmRSMRxnDDFFIrFouM4/X7/1N9N8yl5DJ6vfOUrhJDr16+PRqN524IczebmJiFkb2/PcRyIX6yvr0+5Ph6PR6NR9+d4PD4DI2cMlL6f6kMbWl/aBddOfkEprdVqPpZkb2xsrK+vQ4UZvE3TpyDDMCB51rKsfr+fzWb9ssQ7lFI36b3b7SqKMvca9ZD70oiPSJIkSdLa2lq32z28JJ7eAg0Oppv4JTRjgp9BFdO9dJZlKaW6rlcqlfX19bn0WoKN/VQqlc/nIVVWFEVVVefYCDn8GnYwjOkHjUZjfX0dkr2PPMXKtu3p57xO12exWEwmkyd2FKaUWpYVj8fn8p52u914PM5xHOTeQUcalmVN00QNB8j0N1uW5eeee25mxsyFz33uc5/4xCfOOUgqlWIYplAoJJNJ8H4n2NjYOPPnuFqtEm8NG1iWLZVKpVKJ4zhBEKBr6sxIJBKlUqlarSaTSRDwYDDo9XpwYPUEkiQdVzYXi8VODMJ7J/wanj4Pf+Yznzn/53vBeeyxx84/CMuyjuM0m03Q25E4U3dEjnsXVFU1TRMizJZlHfkFATQaDUEQIJjP8zwszo87LSAI4NaGYaTTafgNmH2ki6Gq6mysCq2GPfpajz32mC8f8WWh2+2eucGYpmnD4VAUxSPbd6uqalnWlJffunXr8C91Xe90OnDms+M4uVwOZmPLsmKx2PjeW6PRuHHjRjab1TSN/HT0i1IKh0VPDG4YxmF1HR4ZRhgMBoe/PkzTnBjWcZy7d++6f36z2dzc3IzFYrIsTxSxzM7bDyI+vgi4e0vzNmSBaLVahJBCoXC2l5fL5UwmMxgM8vm8L/Z0u12GYXieT6fTcECHO3I0GoWZ32Vvby8ajUK/4cFgwLLs1atX4SnYo2q1WuPX1+t1Qsj29vb4L/f398lRO1jpdJphmIkdHdjmmRgW/g8PDg5cO+v1+vb2dr1eP9P/wSRw0N9gMPD+kvDvLTmLUYS4CCQSiUwmc+b0JuhkKMuyX2dZaZp29epVUM5oNLp8+bLrI0iSBB9Q9+JUKqVpmizLPM+nUql0Ot3pdOCpjY2NZDI5sZwWBCGdTk942olEYnNz83B+ayaT2dzcnJg5BUHIZDITw/b7/UwmE4vF4GGxWOx0OpRSSIM/D6qq8jxvmub169czmcwpEml8+fJYQHAePpLRaLS5uTlvKzxxnJ2j0ehwCkSr1fJrJpwYdmImX0BwHiaEEEqpYRiKoixI/63D2LbtS4awpml+ZUoGymAwOC4RNRKJHF5qapoWxO7Ozs7O4p+AG1oNw9t85Pt9GNM0NU0rlUoL261O1/V2u33+/Op+v39+r28G1Go17x67bdtuvNpHbNvmeX7x+yXN55yHGdBqtZ566qn19fUTzzoGqtVqpVIZDAaLmRAiiiLDMBCSRZBxQru39PLLL5MHnQC8AJt+iylgQohhGHhmBXIkofWl3/ve95IHHXlOxHEcaJfTbrcDlYplWYZheFzWgjGUUuiAuYzFxsgMCK2GT4VpmsPhsNlsUkobjQbs0flONpsFd12SpFgsNiV+1mg0OI6zbdswDEEQDMPAI9GRY5lzXDwwYC8+nU57uXh7ezsajcLGOpSS7e/vO44zGAwsy/LFnnw+7xpTLpfdOli4S7fbHbc8Go2CAel0Oh6Pn2rHH1k1QrsePtXKFhbDENh0xpr46LruHp48wZFp7gDLshOxX8uynn322d3dXXhomqabcAsPDcNwXeWNjY1sNgt5f3BK8zK2+0Bmx7y/RILiVPNwPB7f2tqCn/P5PPHQSGE0lYmLIWtvPEHvuMwBMHsivw9BpoDrYUII6ff7buZDp9PJ5/ORSESW5Smnq0amcvj6ZDIJCXqw9uZ5XlEUQgic4TqxNj4yfCXL8ontoyilOzs7uAW1UoRWw6fypZPJJGgV6sUqlQqUvJim6Uvm1njNOkg3lUq122048wHKbuBZlmWj0ehwOISH3W7XTe7lOO7EFCtN02DY89uMLAuhzfHodrvXrl1Lp9Ne9opM04QeK7Ztq6oKE6Zt26IoTu9N4R1o/hSJRCRJKhaLLMsWCgVBECilcF/3ykajUavVEomE4zjxeLxUKjEMAxd4yRmC1Tse5rQ6hFbDjuNcuHCBZdn9/X2Pc7Lz090CQNUsywaaYAxJ2pIkTdn+NU0TvmW8bCyjhleNMPvSm5ubJzZ5mnjJ+EOo//ZrHj4OcNqnOwtQ4O5OwoZhdA4RqJHIIhNaDRNCZFmORqNHto/wgqZpMPv5a9UEiqLAMYhTrhFFUdM09xqGYQ4HJwM1EllkQutLA6qq3rx5s16vL0WxzhQYhoFt5BNPWkNfetW4GO43O5VKQVjoscceO3MfqUXgq1/96ve//31Jkn7mZ35mymXQvhxSrJeiThg5PyGfh4Fut1upVMJ0Ng+CuKyEhhEkxIQ5poUgqwBqGEGWG9Qwgiw3qGEEWW5Qwwiy3KCGEWS5QQ0jyHKDGkaQ5QY1jCDLDWoYQZYb1DCCLDeoYQRZbmah4fv/619n3v+x+v/2evQRgiDewXkYQZYb1DCCLDcnavj+d/7gn/6dX//VX0o+8fN/4Zd+TfrHzZdfhSd+8ke3fuX9n/rCl3/7s+vXn0wmf/Gv5tXevQe1yKM/+a//8tN/+Ref+IXkL//13/i33/xxgH8Bgqw2J8/Df/ra5V/+9D/bef7OC7u17Nv+yz/67LN7rz147rU/VP/9QzfV1n/ba1U/+PLn/8G/+u8/IYSQ+/bv5j/7xT/+lcrv3/mDnfw7v7L9uy/fD/BvQJBV5kQNX3zv+mf/5q89+b74n/u593zwk/Lf+/AP/tOXv/GTB0+yf+Pvf+ov/uxD5OLPfvDjf+093/3a1747IuT+N//D7+29ffOf/MZHn3jH29+d/sxv3nz/CJuFIEgwnHjuoXPvm7//z39L/XLX+pN7r993HBK5tHHwIMJ88R3sz12EHy9cil4iP/rhjxxCXnn5pe8+krz23jfGvvBY6sl3XvxaQH8Agqw4J83Df/rV3/rbt74e/7v/5st/9M1vf9v6n7+zcckZvekYX7jw0ycoOA/+ufjQQ28+89D4AwRB/OQEDd//vy++ePDz4t/66Psfv/Twhcj9/re+/cqJbvEjf/4977j37Ze+/8Zs7dz79kv/B9fDCBIMJ2j44p9NJP6M/YfG/quEjH74P9Tf/L1vnbyyvfj+X//41Zdu/4v/+J3/R8j9P/7Pyu90XsX1MIIEw0m+9KNP/cPPf+z+zic+9EH+Vz/x+e/+lU999PLJbvHFK9lnnxXv/fbHPvQh4S99+t8xH//YL5y47kYQ5Exgf2kEWW4wTwtBlhvUMIIsN6hhBFluUMMIstyghhFkuUENI8hygxpGkOUGNYwgyw1qGEGWG9Qwgiw3qGEEWW5Qwwiy3KCGEWS5QQ0jyHKzQhqmlHa73XlbgSA+syoabjQaiUSi0+nM2xAE8ZmV6K+haVoul6vX67lcbt62IIjPhL+PR7vdXltby+fztVrN+6ts29Y0jRAiiiLHcUEZhyDnJvy+dC6Xi0ajsiyf6lWUUtu2K5VK6L/jkGUn5BpWVbXf74uiGIvFTvVCjuOSyWQ0GsVJGFlwQq5hwzAIIZlM5gyv7XQ6PM9HItjeHlloQq7hfr9PCDnVJNxutxVFsW3bMIx0Oh2YaQjiDyHXMOBxTUspFQTBMIxkMimK4nA4FAQhaNsQ5JyshIY9IkmS4ziVSkUQhEwmE41Gk8kkPEUp1XV9vuYhyJGsxP6wlzVtt9ttNpu7u7vw0Lbt8cUwChhZWFZCw14AlfI8TwhxHKfT6RQKBfdZSZLmZBeCnAD60m+STCYh+mXbdr/f53leURRCiKZpoihSSudtIIIcAWr4DViWBc+ZUnrjxg1CSCqVMk3TMAyO42zbRg0ji0nIfel79+55vDKXy9m2LUlSNBpVVVWWZZ7nq9UqKJkQwrJsgIYiyFkJc7604zjve9/7Xnrppbt37374wx8+8zjFYpHjOJZlYbWMIAtFyH3pt7/97YSQ119//TyDmKbJMAz60oexLMv1U3xhMBhA9MGvAVeBkPvSwDl9DU3TKKWYOH2Yfr//rne9y69VRrFYdOMOjuNglqtHwqxhx3F8WSngSvg4eJ73cX0BuwDlchlbNZyKkPvSAH6jBwGklM/bCiTU8zASHO12e2dnh1LabrcrlcrhC6Z7QPit6iMroeEQx97nhSzLhmGkUql4PH74Wdu26/X6lJevra1NccLx/ToVK6Hh6d/6lNJPfvKT3keLRqNf+tKXzm3UcrOzs2Pbdq/Xq1arh59lWfbWrVuzt2o1WQkNT/9ef8tb3lIsFr2P9vDDD5/boqWHZdlyuRyPx48rz0RfemashIan8/DDD3/kIx+ZtxXLx87OTj6ft2272+1ubGyMP2Wa5njFyGEqlQomzPgFavhciKLY6XQsy2IYZt62zBTTNPv9/o0bN1RVPVzUxXHcmfeHhsPhD3/4w/GJulwu37p1q91uo+yPZCX2loIjmUzG4/EVPD4CegYqisKybCKR8GXMYrHI87xpmtFodHznOZFIpNPpdrvty13CR5jn4UgkEvS6q1KprNoM7OJjiiUAOR6HkSQpHo/7frvQEOZ5+Ax5Wu122zAM27a9v8Q0TWy7FTSGYWAS9XGEWcMuHmdj27abzaYgCKqqehxZUZSJcA7iO6ZpDodDvzz28BHm2sPRaCQIwt27d+/cufPUU095eYllWVeuXGm1Wji1IstCmOdhd/r1viqGBGCMfyJLRJg1TAj58Y9/TE6Tuwd94WHPM0i7EMQ3wqzhSCQCvXguXPD6Z3Y6neFwWKvVdF3H2RhZCsK8t0QIefzxx7/1rW959KUty+r3+8ViEVIvI5FIo9GAkJVt21hFjCwmYZ6HTwsshiFJENxv2JM0TRM3NpCFBTX8Jp1OJ51Ow6Q9HtziOA4TDJCFBTX8JrZtuwcd6roej8dBvcViETWMLCyo4TdJJpPQkG0wGNy+fbtarcZiMTjDBZvOIAtLyGNap6JUKomiWKlU2u327u6uIAiUUmgujRpGFpY5aLhcLvf7fZZlZVmWJIlhmOOS3f3CY0iZZVnTNB3HKZfL8BuGYeCslgCNQ5DzMQdfWhCEnZ2daDRKCNF1Pbg9Gzi+MBqNHtnz6TgmNqIg217TNJ+NQxCfmIOGoVgPKkWHw2Fwx4JC+dH4McJnABqXn6pZD4LMkjlouN1uR6NRWGRms9ng6m8hHJXNZs8zCMMwsiyvbJEwsvjMrm6p3W53Op2rV6/evn3bcZxmsymKIrRuCOJ2juMkEgmGYXBbCAk3M4ppZbPZSCSiaVqxWGw2m9vb2/D7KQKmlE6PdYmiOCXaVCwWoYP5WU1GkOVgFhqGE3QsyyKEZLPZWq0G0gVfdwrXr18ffzixrJ3i36qqahiGYRhYOI6EnllouFarwTxMCIGOZ6lU6sRXMQyztrZ2htsZhvH000/v7e2hgJFVIHANt1qt4XCYyWTgoWEYHhfAlNLpRbyJROLIfSme57e2tkRRbDabWGyEhJ4ZrYevXLkCP3Q6nUKh0Gg0CCEbGxtQEnRkDzpK6fTsqFgsdtxTpVLJtm1BEPb29qZchiBhwAmYg4ODaDTabrdHo9Hm5iYhpN1uFwqFg4ODoO8bj8d5nj/zCKPRaG9vbzAY+GjVeQB7LMtyHKfRaHS73XlbhCwEge8Px2KxdrutKArU1jcaDVmWOY6LxWK2bVcqFSgzCOK+2WzWMIyzhaZ1XU+lUteuXVuQpjzgVly7dg1cmGKxeGJEEFkR5tnXslgsQmxZluUgxocmlZlM5mwfd13XNzY2RqPRghzwZdt2IpGA1UEikbAsC1f7CJlv7SFMJsG1yIC4dKfTOdv3FPTHWxABE0LAoUilUmAYChgB5l8/HOhnMZ1OU0r7/b73l8AagxBiGEYymQzIMPcuHi8mD9qMgGHBJZkjS8c864c1TZMkSVGUgHxpF+8d7VRVhVKqSCTS6/XcIkR/ce9CCPnOd76jqupx+SqGYRSLxXQ63ev1+v0+BAVN08Q6KuRN5hRLcxzH6Xa75XI50PgqTFztdtvLxfV6PR6PQ8A8n88TQoIIno/fxe2/dyTgPO/t7cGrCCGtVgt+77tVyPIyz3mY47jFKa8fDAZPP/301taWu5+cTCbdnxVFOW72a7Vax62ZGYaZyEiDu2xvb8PIlNLxfFJd1zVNcyNwkiRtbm6OjwDpMdj4GhkHe/G8gaZplFI3wOYuPoHp9cOO55Ut3MUVIVRfus+Kouga0O12+/3+eH7bQgXYkMUBNfwGw+EwGo3CGnUwGPR6PVmWITNU1/V+vw9+72FOldQ9HA7j8Tjcxbbtfr/P8zzcxTTNTqfjFlRDySRMwo7jwFOQ0MayLKVU13UvkS0IOni3EFk65h+XDhTvMyS0JYCJDnxmnuer1apt27ZtNxoNX9riwdcE3AV8ZkEQFEWBxNJut+tWO4POYXvMNM1+vy8IQrVahb+IUnriIckQD8vlcuc3G1lkwnx2KSHkAx/4wDe+8Y12u33iGhJc3GKxCDWStVqtVCqNRqNcLscwDMdxpzpYfMpdOI6TZbnb7TIMA3e5fPkyKA3aWbsxalEUE4mE4zjRaBQORh6NRrVa7VR3jERC/hYj84xLB839+/effPJJ8iDNwwutVgsixgcHB27AvFqtlstlSFT2hVarBWnY+/v77l1UVc1msxN32dvb29/fdxxnMBi44ehqtZrJZDxGp8P9FiPODPKl58sjjzxCCHn99dc9Xi8IAkSMY7GYGzOHlaePW7KCIMBkm0gkxu8iiuJEWmgqlQJ3mmEYcCVs22YYhmGYgPLMkaUjzI7WaDQSBOHu3bt37tx56qmnzjyOruu2bUMrbB/NmwAaj7g55FMAhxx8e3jVxAU8z7trh0B9aUopxMz8+p+BpmuEkGw2e2ILh8FgACsLlmVXOm43b0cgQO7fvw+7r3fu3Jm3LX6iqmqhUFBV1cvFgb7FL7zwAnmQhXJ+CoXC1tZWq9WCBJvd3d0pF+/v7yeTSVhrJJNJjuN8sWEZCbkvHUoopYlEwktculKpEEIqlUpAZ82sra3t7+976ax0IoPBYGdnJ5VKCYJQq9XS6fTNmzenXK9pWq/Xi0QiiUQik8mYprmy/Q9XYn84ZKkRcA7jiSlu4FEHlPLt4lfTMtu2KaWqqkKWSzqd7nQ67XZbEIQjr9/Y2Oj1euM+/Mr2AF+JedgJ3Zp/7jmqlFJJkiA1xZcBU6nU3t7ezs4OPBwOh2RqTRvHcVA3Qim9fft2Npv1xR04FZTSbDYL/wnVarVQKPA8P4dylHk78wES1vXwIpDP5/f397PZbDabPfxst9tNT2X6xthgMGAYJp/Pn2jG9evXGYYRRfGsf8e5KBQKlmWl0+l4PA57hNvb2yzLztiMlfClp/Paa6/duXNn3lbMiEcffXSia/fZYBgmkUjoul6tVg8/y3Hc9NXp9NVNsVhMJpOwmJ8OBLEzmQwc/TNjd5phGDgrE9pLkQcZu7O0gazIeng6r7zyyhe+8IV5WzEj3vnOd/qi4UqloqrqcDg8rg3LmWMQkE96qiBcsVhcW1urVqteZO8jkG833npZ1/Xg+kYcR5g17PFjxDAM9pc7A81mc3NzE3obTixcbdtWVXXKawVBODL7VVVV0zRhVayqaiqVOm7lXygU+v0+vHHwRsOcPGMMw4jH4xDYg1KZcrns5uHMyIgZ++6zBNfD09na2pJl+Wyvhazy3d3der1er9cPXzCaypFjtlotcMLb7TbUcsGV7XZbFMXxJsGj0Sgej0ejUXi4u7tLCNne3nYcZzAYHJmIurW1BReM02g0stnsRPthy7IkSTrcmmJ9fb3RaEz8MpPJbG5uws/b29tgUj6fh4zdzc3Nwy/xnTBreDQaQQ0wavhIYOUG+din5eDgAGahTCbjizGWZU1MXK5EocR6QoH1eh1iY/v7+xzHJZNJkCL0EgoLyKkAAARgSURBVE6n0xPjR6NRd0CX9fV1cihHBQ70m7gdDHs4XhWNRt1km+3tbY7jFEXZ2tpyHGdvb48Qsr6+fur/i1MSZl+ahHFXyUdM0wTH7wx7vLFYrNVq7e3t+bWVwjDMc889N3EL+EFRFJ7nYeZ3yeVyLMtC5Gxzc7NUKsHvOY7b29s7fGJmt9s9vLba2dmxLGtiU6pUKqVSqQlXH4YtlUqO44yPoyiKm+ZZKpUgygV72rBbNotDv4L+kpgjri/9wgsvzNsWP/GxA1mhUDjOs10oVFWdnno5zmAw8LIvdQZct3mhwByPn0LXdZ7nE4lEQMmJ5wcOoPBlKDiDcvGT2KBpycbGhsfrZVkOomm5LMuLWVkRZl/6DJ9OURQppblcbu6JUMcBjbV8GYpSOuPNmDNzKo9dFMUg2gaOV4MtFCGfh3/wgx8QQl599VXvL+l2u8lkcmGTb31sEL+Yn8jDnHafJqC/a2H/u8Ks4UgkAj0AHn30Ue+vmuhouVBAG63F9OiQeRFmX5oQcunSpVNdTynt9XrJZLJQKAyHw2KxGIRT3e12IY0hnU73+/3pjW+r1aplWZFIBKKdKGBkgjDPw2cA9iqghDWTyXiPo3jHNM2NjY1sNqsoys7OjrspchhKqSAIlNJarRaJRDY2Niil0wWPrCCo4Z+i0+lcvXoV5jqGYWzbhs19XdcPbzkCU4L+hy8eDAYgYIgtsyw77rfruj4uUSgJqlQqkUiEYZjt7W1JkvC4Q2SCkPvSp8UwDLfoHCpjQTMsyx7nVEMo+8ineJ6fCPy6rbnc241reLwvlGVZzz//vJt1HPSxcsjyghp+E8dxxs86hHT2WCwGEj1uAmw2m95vAVuyMBSsvd3bwV3cbwrYoD4yuuY9pd77gY/I8rISvrT3HA/yIMWPUgrTIKUUNvd96REP7eZh4xpUKggC5ANrmibL8pHnPADuHmmxWDwxBaVYLM6npwQyc1ZCwx6JRCL5fN4wDNu2BUGo1+uCINi2XSwWKaW+TGiiKNq2bVlWo9GABTZMqrZti6I4fsiDIAhuHXy73R4PR0Mn6uk3ghzj8xuMLD4h1/BpC7JrtVo0Gm00GvV63T0/BergfLEnlUo1Gg1N0xiGMQyjXq9D90mWZW3b5jhu/JvCNM3Lly9DWYKbWz+DI9eRJWMWSdnzA47ePmcGfCaTgePOfTLqaLLZbKPRgC4Wx6GqqmVZh2vojqRcLgdtM7IIhDymJUmSLMu3b9+WZdmtZTstLMsahhFEGv3EXSb2lg4DB526lhw5IeMsvWqE+awWoN1ur62tbW1tLUt+/3R4nlcUBepUp18JYkZJh56Qr4cJIRCdunXr1vQOT8sC9FL3cshDv9/v9/tz6TKFzJKQ+9IAnCGcy+Uikciy5xvDgeMnOvbQvpwQMhgMZmIXMjfC70u7DAYD27ZnfxoAggTKCmkYQUJJ+NfDCBJuUMMIstyghhFkuUENI8hygxpGkOUGNYwgy83/B7xF+13+k7dqAAAAAElFTkSuQmCC)"
],
"id": "a832d3b5"
},
{
"cell_type": "markdown",
"metadata": {
"id": "acbd761b"
},
"source": [
"Where $c'_i$ and $d'_i$ are modified coefficients for c and d in the original matrix. Notice that b and a coefficients are not modified in this process (as they are eliminated)."
],
"id": "acbd761b"
},
{
"cell_type": "markdown",
"metadata": {
"id": "633537c9"
},
"source": [
"<h3> <font color='blue'> Reverse Sweep (Back-substitution) </font> </h3>"
],
"id": "633537c9"
},
{
"cell_type": "markdown",
"metadata": {
"id": "09da384c"
},
"source": [
"The completion of the forward sweep results in a value of $x_n$. This value can be used to compute the value of $x_{n-1}$ which in turn can be used to obtain $x_{n-2}$ and so on using the transformed coefficients as follows: "
],
"id": "09da384c"
},
{
"cell_type": "markdown",
"metadata": {
"id": "27678e23"
},
"source": [
"![backsubstitution.png](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAT4AAABCCAIAAABBzS1pAAAAA3NCSVQICAjb4U/gAAAAGHRFWHRTb2Z0d2FyZQBtYXRlLXNjcmVlbnNob3TIlvBKAAALi0lEQVR4nO2dv56q1hbHt+fe4lQBu1QBzwtsOFWqiOcBIjNVqqCpkmbU9BHn9GfE02fEPMCILzBAHkCgSzXsdKca8AXct1g3fLyoMzjqONxZ30oRNpsNv7X2H1yrwjknCIKUjTenrgCCIE8BpYsgpQSliyClBKWLIKUEpYsgpQSliyClBKWLIKUEpYsgpQSliyClBKWLIKUEpYsgpQSliyClpEzS5Zx7nkcISdN0MpmcujoIckpKI91ut1utVlutFiHEsqxOp3PqGiHIKfn3qStQFMuywjCUZZkQ4nkeaBhBXi2Vsvxfl3P+5s2b6+vrdrtdqVTiOAYZI8jrZIcOcyZyzvlzCh5OB6NcTdM8z6vX66hb5JVTqMOcpmmr1aKU+r7f6XRmsxljzLIsRVGKn+kBtVcqlW3n7Xa7hBDG2GKxkCSpVqtNJhPsLSNIIel2u93xeFytVjVNG41Gw+Hw/fv3rusWl65pmuA2gZxWBUGYzWbrR2ma1mw2Ly8v0zRVFKVer8NGTdMKnhdB/m/hj5Ekieu68FmSpKurqyRJxuPxowfuydXVFSHk/v6ec75cLgVBuL6+PvZJEaQs7DBNFcfxu3fv5vO5qqrHNCb/RVEUURTBVwdB8P79+7u7u1qt9gynRpCXT6EOcxAEqqp6nicIAujW8zxZlmVZtiwLPsxmM0qprusbSzBN0/f9beWLoug4Tm5jFEX9fh8+u64LA13HcTRNsyxL13XGWBRFhmHglBXyGnnUL19fXxNCbm9vm81mvV7nnC+Xyx9//JFz7rqu67qEkOl0yjkXBOGBcpbb2bg/pbTf73POkyRRFMUwjLu7u0+fPpmmOZ1OJUkKgsB13WazmR0Sx7GmadivRl4Djy8OybIsCILv+5RSxthkMjk7O+v1evCTKIrgbNM0fbicynY27g+O2rKsbrdLKY2i6PLy8vz8vNVqhWHY7XYVRWGMrbrcNE2DILBtu6jdQpDyUkTf8/n87u4OPt/e3q7+NBwOTdPknI/HY8Mw4jg+oF1JkmQ+n2d1yAqnlAZBwDlvNpvT6TRJkuyQ5XL522+/HbAOCPIyKfRKhqqq2fxQo9FY/cnzPFiq8TxP1/XDejxRFLMpMVVVMwfLGIN1Kc/zFEVZHSenafr3338fsA4I8jL512Aw2Of4L1++/PDDD4QQURT/+uuvn3/++e3bt4ep2hYYY19//fW3335LCHn79u2XL19W39D45ZdfLi8vRVE8ah0Q5OSU5h3mgqRpirpFXgP/b9JFkFdCaf6viyDIKihdBCklKF0EKSUoXQQpJShdBCklKF0EKSUoXQQpJShdBCklKF0EKSUoXQQpJShdBCklKF0EKSUoXQQpJS9duowxz/M8z0vTNE3Tdrv9aCSdffA8D+PjIKWgBNIdDAaNRkMQBNDVUaXb7XZXQ72XAsZYo9EIw/DUFXmcJElGo5FpmiesA2Ps7OxMUZRGo3FCM805D8Pw7Ozsyc/bS8/0p2kaxLKpVCqe5zWbzeOFbk3TNIoiy7KOVP6RYIwJgvDCI9pyzhuNhiiKYRgahnGqaiRJkqXyGI/H7XY7juPLy8tnrkYYhp1Op1Kp+L5/cXHxxFJOGxqrCFlUV0opRI09EhAj9njlI8vlsl6vw908Cbe3t4SQLNwvpZQQsi2c8PNUJhensTg7dJhd102ShBASx3EQBE80FYWJ4xiGuFEUKYoCY91tIdqfDOfcdV3GGCEEYsSu7xMEAVxvmqYvrTvNGDvsvYDJBfgMt+CAhRdntamzJ2H/YkVRlCQpC4G0U767l0ahDnMYhoPBQNf1Xq+n67ogCJPJxDTNs7OzR4/lDwbQ2RiEGRILQhODlhqNxjHSYUOo506n02q1Wq2W53nrSRggxWEURYvFIkkSSqllWeu7nYQgCEajUZqmmqZtNDq7Nj5jDMqxbZtzDpGue73eM1jqHN1uN01Ty7IopYvFQhAEyHexZ+AxVVXBTANhGEqStC0S+AunaKY/13UrlYpt257nDYfDX3/9NbOC4A83jrWy9Ajb+Omnn9YP1HVdkqTRaEQIiaKIECKKoqIoh03w1+l0oigC035+fk4ptW0792SYpmkYRqPRmE6n5+fnSZKA5cp2yMVw30bB3XZlMBg4jqOqahzHG08KmSu28eHDh1yTTqdT0zQnk4lt20EQ1Go113U/f/58pPpvYzweN5tNCJpvmqaqqpzzjx8/BkGQiyW8D67rRlEEvdYyUki6lmWBZQrDEJrSdd3srtu2Del21w9sNBq7tvV4PPZ9fz6fw1fGGOTmfPjRYYw9MFtIKc11EKbT6efPn7OzrF7OKqqqQv0h+YMoiqvyZozpup6b2s095eDJjzQx3m63IfHScDhc/1WW5Y8fP+5UoKqqqqq22+1OpwORt8FHrTf+rg2+E7VaTdM027YNw8hyXMH29Z1N09zmNiVJarfbG39ijLXb7el0ekBb8MwUki4MCYIgWCwW8IivPugbu2oZu/bZHMeRJAnOmCRJFEVF1hJkWf7uu++2/fru3bvcFkh9lg11tvnzbGjted66BZFleX1Jxrbt1dDWkA0YehAHR9f1wWAgSdK252/Xxtc0Ddo8swWe533//ffrx+7a4DuhaRrn3Pf9rBqO41BKN5pv2HljOdVqdeN2GJFNp1PosJQ0fWQh6YJxymX6gySalmVtc7mEENu2x+PxAyVPJpPc/VgsFoqiwFMFtlbTNMdxQEXg9DZOVn348KHItQBRFFFKN1prXdezoazrupBXKbMgaZqGYahpWhiGtm23Wq3TTnWAa2KMMcZyBogx9vAyDMxc5DZCm4Mt4JzPZrPhcMgYg+zkq3vu1OC7EoZh5ic4547jtFqtjZe5q9uEyc7hcAhPcq/Xm06npRzuPjoHDeOl+Xxer9ch09/9/T3k14OhbJYB6CAYhgFpBCHBH6X0/v4etnDOcxmGnszV1ZWiKNnX33//HTIncc6HwyF8WC6XhBDDMMD6wDX2+32YzTdNM8u3tMr6Fs75ejt3Oh1RFPdM0QSzR3Ec9/v9LCnUnlxcXMBd5v+sXiRJAmkWD1L++uJQvV6XZTl3W6+urrKFOhjXzOfzfr+f5aB6MoZhNJtN0zQHgwFkooPtFxcX63fk5uaGEHJzc7O68e7uThTFTqeTK5lSuvpQAQ/c6PXFISi54OLZ415X0zRK6Wg0ajabs9lsMBjAhDNZ6T0e0PMMBoNWqwUTjIZhWJbV6/UgsSCkGjpIcoNer+f7frvdliQpDMNmswlXtLo+VKlU6vV6tVoNguDi4mIwGMiyTCkFM9/tdmE8Rv4Z7cNRnudlHebVnnNuDCyKYpqm4/F4n/cBwLSNRqNvvvnmUL2+KIoyV1yr1QRB6PV69Xp9//I1TVssFmEY+r7/559/wlw9+WfckZuCylWDUvrHH3989dVXe6Zlt217MpkQQmazGWyBcRAhpFqtcs6DIMh1A1cXk7I9BUFYna0EZFlen9HYeKPDMOz1ejDaOj8/VxRlOBxCZ1MQhKJPeBF9c86zZevc+nWn0xkOhwf0utvOG8fxeDw+7CsTuey+4EU3uk2+duFxHFNK4zjOuYttXtd13dzG29vbg6QCPtUbBQfk06dP+7vTgjyQ2/nm5ubJL0g8wKFudI6ir2Rkg4HcqCAMQ13Xj7dwn50ORpiHzSeUy+4LVnPbKXIXbts2ZDZ8uEppmvq+TwiJomh1RZEQ4jjOQZa7SjlO+1+iKNrTnRbngdzOvu8fY+biUDc6x75/P4D3Ex6eZD4IoJODv5WROwVMhxTZGS65yM6cc5gRWN0YhmGlUinp3OZhsW17n5WkQ+E4jizL2yaln8zxbnSZ0oXBoMi27SNZCliAbbVa+7j3Z357AXm1vPQ//a0Ca1HHc7wguT3TfKJukeehTF4XQZCMMnldBEEyULoIUkpQughSSlC6CFJKULoIUkr+A5V19QgpZ5C2AAAAAElFTkSuQmCC)"
],
"id": "27678e23"
},
{
"cell_type": "markdown",
"metadata": {
"id": "85b402c5"
},
"source": [
"<h3> <font color='blue'> Python Implementation </font> </h3>"
],
"id": "85b402c5"
},
{
"cell_type": "code",
"metadata": {
"id": "5e7017e7"
},
"source": [
"# Python Function to Implement Thomas Algorithm\n",
"# Venki Uddameri\n",
"\n",
"# Step 1: Import Libraries\n",
"import numpy as np\n",
"\n",
"# Function for TMDA Algorithm\n",
"def thomas(a,b,c,d):\n",
" \"\"\" A is the tridiagnonal coefficient matrix and d is the RHS matrix\"\"\"\n",
" N = len(a)\n",
" cp = np.zeros(N,dtype='float64') # store tranformed c or c'\n",
" dp = np.zeros(N,dtype='float64') # store transformed d or d'\n",
" X = np.zeros(N,dtype='float64') # store unknown coefficients\n",
" \n",
" # Perform Forward Sweep\n",
" # Equation 1 indexed as 0 in python\n",
" cp[0] = c[0]/b[0] \n",
" dp[0] = d[0]/b[0]\n",
" # Equation 2, ..., N (indexed 1 - N-1 in Python)\n",
" for i in np.arange(1,(N),1):\n",
" dnum = b[i] - a[i]*cp[i-1]\n",
" cp[i] = c[i]/dnum\n",
" dp[i] = (d[i]-a[i]*dp[i-1])/dnum\n",
" \n",
" # Perform Back Substitution\n",
" X[(N-1)] = dp[N-1] # Obtain last xn \n",
"\n",
" for i in np.arange((N-2),-1,-1): # use x[i+1] to obtain x[i]\n",
" X[i] = (dp[i]) - (cp[i])*(X[i+1])\n",
" \n",
" return(X)"
],
"id": "5e7017e7",
"execution_count": 1,
"outputs": []
},
{
"cell_type": "markdown",
"metadata": {
"id": "ccc366f2"
},
"source": [
"<h3> <font color='blue'> Illustrative Example </font> </h3>"
],
"id": "ccc366f2"
},
{
"cell_type": "markdown",
"metadata": {
"id": "46422ba0"
},
"source": [
"Use the thomas function to solve the system of equations shown in the matrix form below:"
],
"id": "46422ba0"
},
{
"cell_type": "markdown",
"metadata": {
"id": "b3caecd8"
},
"source": [
"![example.png](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAbUAAADXCAIAAAASg1hvAAAAA3NCSVQICAjb4U/gAAAAGHRFWHRTb2Z0d2FyZQBtYXRlLXNjcmVlbnNob3TIlvBKAAAgAElEQVR4nO3dbWzbeJ4f8C93b4dMMRsxs5MRM+1GHMyDmR30zCxQmOmLmmmvFYMCtXIvanlbwPId0CgLLKxsX1huC4yzL07KvRgrfTFW+mIsF+hYeXNWDti1XBRrem47lotDTA/aMXPXG9OL7Zjaa2NqX4ypfTD7QnImiSk/yKJEOf/PK0OUyZ/48CP5f6QcxwFBEASxz9c6HQBBEIRPkfxIEAThjuRHgiAIdyQ/EgRBuCP5kSAIwh3JjwRBEO5IfiQIgnBH8iNBEIQ7kh8JgiDckfxIEAThjuRHgiAIdyQ/EgRBuCP5kSAIwh3JjwRBEO5+z+P127Z94HKGYTyOgCBa5bCz+ains22WCrlcvlB6RF2dXJxUvLsGnNXMtZGcxYlyNBaPyjy52o7J8dJ6qufgrdODc54GQBAtszM7QB98Ovem1g9dy9b8WDhEB4TB1Nz6dhuidhxnZ+OnH97oC9KsGJ9d32nPNpu1s7E4PTYc7usR5NRyMyvY/Xx6sLent39wdGr+5DvY6+dHAKHwaFRscN9iJcH7AAiiFRgxmhhrdL5aWv7ewqGrMHIRaaTIDue1XIRvaXAHYfirI9lSZDh1/dr3rmjG/GJS8uNz5M767K2RW7lHb0RTqXy22Wdd6o1YXottrc7eHY+L42eup2an4iLbdFAnTbAHWk/1AOHpNt0mCaJjtqb6D31+XJ/soxEcnNtqV1DP2f3yk7FemuLji/57iPzy4fv9AbBScrF1e+fLhx8MhMAI8cWmMxCpnyGItrCL6fQKE85kIlyHIqDOSO9Nxl7bzN3Omx0KwZ2zszw+9MOPL8Tn5lNy6/bOmcs35+bf79u4N3SzYDW3CpIfCaIdbDVfKAeUeMeyIwBQjDwcvVgtzRT9lCAdq3An91fc8OSk3PyLsDtKSEwm3infv53Vm/r/rsqPtlHMxBWRZxmGYVheVGLpgn5IjWKXMku5hCLG1U7H8SyzlEtEJIHjOI7jeFGOpQt6kzfmjuhk/IamVSDKHS75oyD2iy9VtSWto2E8w7GKM8VfhaI3Dq/H316dTUblvQMoSNFkXjv4AFKUeOPG72PtfqG5BNmyl303LSx/3FmdGgi5xE+HBqb9Xid3LNurs6PhEA0AgeH5Tkfzle3Fsb4AEBqYWq4dz63FVH8AdGiwO/a/x/EfWv44PxhAYLDjB3R395PREEXJU50qBN1nd/enN4IUE5k75Ch8+fD9/iAoVk7ViygfP/xgIASKUyaXD8ow9fUrs82koS7JjxvT4QAQ6B0cTU1OTU2mRgf7gl+1tQgNz5+GKqDt1dmxwXB/X0+w/rt8lB83pvoDAN2bWn36LN5evBECOlnlcFSex39YftyZHaARGm2qyUpL7X7+fh9FiUdoidQmu48/CoPiEwfumt0v/mwwCDxXt7T75cOxXoBilamNA/738/d7KUpo6hd3RX7cmg4HAn2p5+4SG3M3evdyZM/Y6gm30Xnb66u19lo7i8NBf+XH9ck+GkBgcO75Q7k11U8DCAw0dXdulzbEf1h+3J4N0744TXe/+KDfX/nx8w8OS9i7X3wYDgAUo0w/dyPbffxnAwGAYuTGGfJICbiBbih/1HMZTc4WktKzhbd8JFvI9Ncy5KNC3kcFKs1hBVFgAYDhhE6W4e9jFyfSK1UgoMSU58vPuUhMooHKg4m0bw+Aj+JnfdDu0H9d1mwADNuwasZxSnduL1RA0fLg87VbFKsMKwE49tLt28VGVREMw9CH931y1QX5US8UkEhH3VIGH0soAQCAYRins6Km86xCtlAGQEuKS+UCJykCADzKZdU2B3ZEPorfb4nJRxrvGsdW7+U3AQpCeH/9NsVIYRGAUy7MNEqQJ7gjdEF+tLlYJt6g2wIjSHztD/jurnhK2GpBrQIAJ7p2Q+BFIQAA5UK+1NbAjqjb43/ROdDuF8sAEBAF3uULnCCGADiV4v2GT5BN64L8KMbijZtF7WVFnnRU9IimlioAAJ53fetnuPrnZa1ktC2qo+v2+F94uloqA6DA9/AuiynwPRwAOBVtpbk2PAfogvx4IEMzAKA3qvCdDeS0snRtEwBAc5z7XWovv0DXWn56nlyXxm+bWjGbjNUa+/LRp/q77KzPxmWePXPhasY/4XrHsfW1+u9kOfdieZbjarUQhtby1qxdnh+1gloB6HAiRh4fvWHq9UuTYRpULjAsUzs9K6bpv8bi3Rm/oRY0WxgYHuDtSrW6eX8iV8sRO6upq1dv5kubFbtcWjqkafTpYJlmFQBAsw2K0J4c2arV8gPY1fnRVrOFTdC9E+mYryp8TxNrr9qvYSH3VwuspqoIvdWd8fNKPBaRr8jxdKwHAB4VVBPOVn7k+58mSpZlba3O/3R1Otrq3nh+ZNt7Oe+wI+jAKrc6P7ZhfDPPmPmJ3Cbdl8olxE6HQgAAfJNfmuS3+CkIESV059EmdLW0iuw9cWo+ygPgROXFeSI4esVr649f9+ZHq5BILrEDs/lko8El2+nYrav81wrNHcuxwCZwwMn3ZMEBbdg6prvjpyDKUuDuZqVSjP9hOTm/6Idzvb3YJ8UiNmz3ZFm/+CiwwVbvnm59vzbzsXiRTxXzUb7ToQAws8qZ4+H9NvREIxxfL/xueAewbLtWQBTgG9SAdFK3x8/IikgD+NWZSKpRK7dTjeX5WhPnqm25H0F779DSnHsThRPoyudHW0tHEka0oPri0REAK8XGxqTj/Acn8x7F0mKsKAtYWAMqlmkBLgnkSaG4IPrw+u32+MFKioClNVjWi1Absw8FoV/EzBIA0zQbHMB6DQ4nCK2+wXVhfjRy0UhOzJUyLR8srmmMGEuf1jJQQRIDWKsAhu5+epq6AQAISfW2+v7S7fFbhm4CcCqlJR3KaT3LDsCJUg+WHjkwN01g3y2s/jlABaRwy29w3fZ+bRZiSprNqNl9PWkJj0jRWh9OQ9fc3m8MTa8AQFCJ+PPi7eb4HSMXS5hSfwDAI1U1Oh1PB1AQB5QQAJS1vbZaz7B0zQBABeSB1g+u2VX50VLjStKaKObcx2DWi0WjzRG9CBg5HgkCqJbU0v4EY5VUDQBC0Zjsj8KO53Vv/M763ViWz+ZzMYmGA22htP8V295YzMSuDOVO78s3BWk42gM40F12gGNrC1oVoIKRG5HWPzN1T360SkklZsQLDWpkrFIynseL0+ihjRg5meyjgXIhty/BWMVcsQoEwhPJY5W/tlN3xu/sLI9/fyGay8gsKysi4FRLD+rxm4VMToNtqLOTiZGhW/95xYTPWia1EkWJo2PhAJyKOlN8LkE6tjpTLINi+sbGvJhGvEvKH20tE1EytjJhF7OZ4r6lpq4W8rpcnPDbE0BzLKN+FvjlpBcS2YmcNL6WT2YnSgn+yed2aWJioYrgYC7n6xb6XRO/VcoXDE5S5NDWvaF/Vxmbq1VZ87Lcg5VHleJM0Qr3r4+PPBCnp8FAHkr0c48K/+Q/dTpuj1EXYlOpmcvfX3pw+44WSX1VK+vo2dv5MsXKk7MJbyrXmhgz8uhaND7u+nQ4eOgP6fHPiJ8ntJ7qrf+m3tRGp4N5Ymt+tDcAumd4uj7r+vbyZDgIBPsnV7thfgWP4z/S+LiHzP/qLN6onec0zXDK1NPfXR6tTS5CB4KcmPhqvtK96RKeHze2sd3HH4Z9Nj7u+31HmfDh8Sd/0h8AxUX2plN4/NmHgz00WClx8Pytu7s/GaQpLr7YRHD+f7+2inFlZKF8yLfo/kTUj60zjsMoZpKJWESUxtfqn6xNSFI0nkim852fhoxTMpq+nFGsbETkeZ7nhFiBjc2u6mrCJ62sDtYN8YuRaG+AptlL16cW555p7SiNpgZ7aDrwxpWR2fnWT/PXDc5dGVc3Hv6XYboQqx3AS9fu2UpmeX3Zw/3h//drVskaTrbTUbQDryTSSqeDOBAnxTOFeKejaJ7f42eVjGZlXBe9MZTXh9ocjv+cuzyUzg+l27dB/z8/EgRBdAbJjwTRFruA0+kYavwTie95nR8dYNfjTRCEHxycdSqfPKzir9TP2xePOweV//HQwV8tdTySp/g3S3ibHz9XDeDhzyqeboQgOu+XH38K/O+f/Z9Gy+nzr34NZ89/s50xuaFAX/BHJHUOfvHxZw4+/YtfereNv1yp4m8/bmaCSm/rZ755/mXg1Qv04d8kiK5Gc68A9t8922g58zb/DTDvnm9nTAdEUv6ODyKpO/vaOeAbnIdZ4hz3Dfzy9Vea+E9v8+P5d88D/Jt+aT5BEF4JvHMRsC+1P+vUBvfqePOvplE4K7xO4czbAe+28WbP34Hx1sUm/pPUzxBElzK1+dz47ULZcbTcRK5QIlPAt5z/2z8SBOGKE6+NiNdG2tgc8IVDnh8JgiDckfxIEAThjuRHgiAIdyQ/EgRBuCP5kSAIwh3JjwRBEO5IfiQIgnBH8iNBEIQ7kh8JgiDckfxIEAThjuRHgiAIdyQ/EgRBuCP5kSAIwh3JjwRBEO5IfiSItvDPrFj+iaTuRZ1/hiAIAL6an2v5YdVH83O1Zf6Zn1WanH+G5EeCIAh3JD8SBEG4O6350TbVbEyW0808U3vGNoqZuCLyLMMwDMuLSixd0E/nnCFmKZdQxLja6TieZZZyiYgkcBzHcRwvyrF0Qbc6HRRxdNurs8movHcABSmazGueHsBTmB/NUjYu8/zVmzNLhtnpYJ6wtWxEeOParXsLa5uVarVarWyuLcyMXxeFSO5U5UhLyycUnr8ycnfBsH30wyw1KQlXRgqI5XTTNE2zlFGM9HVRjJ6u/X9q7axOype++6/ulZWsZpqmaa7ODtv5oe9eupYpeZYjT9X8XGYpl87kNcPQVsrVTgfzDCMXkW8uoHdwNCoJLCy9VMgXVspVoLr5YESJc1pOYTsd5ElZWj6dzpUMQ1/Z9NfuB4xsJHJnxe5NFfJxsTbfMCcnCwVDvHpvRAar5SNch0MkDuBszY1c++HSL/n4T+eTcu0Anrscn53bvPLdOz+8NsSszsd5TzbspfVUDxCe3vZ0I09sry9vbDuO42zPDtQmGw+NrrZn0wfbmg4HAn2p5Wf3w8bcjd69OdF7xnwR6Ilsr6+ubzuO4+wsDgcBAIHh+U4HVbM+2UcDCAzOPX8qbk310wACA7MnPEm3pvqBvsmNRsuNH3wLYMI/PtlWTm53d+MH3wJ1Rul4JHW7ux9HXgbFfm/tgO988WE4AFCMMr313JLHfzYQAChGntpovIUfvQnq69+ZbCK60/R+zQoSzwIAy3F+ehjTcxlNzhaS0rNB8ZFsIdNfy5CPCnlflZQ2gxVEgQUAhhN89SxmFyfSK1UgoMT2PaRzkZhEA5UHE/4qqia+4jilO7cXKqBoefD5p3yKVYaVABx76fbtogfFJKcpP/qUXiggkY66pQw+llACAADDIJO7e8QqZAtlALSkSMy+pZykCADwKJdV2xwYcSSOrd7LbwIUhLC877mHYqSwCMApF2Y8SJAkP3rO5mKZuOC+jBEkvvYHmP3XLtECtlpQqwDAiaLbWwUvCgEAKBfypbYGRhyJA+1+sQwAAVHgXb7ACWIIgFMp3m99gvQ6P/q351DbiLH4/tvenr2syEsNMihxQppaqgAAeN71rZ/h6p+XtZJxkg35rtfeKaGrpTIACnwP77KYAt/DAYBT0Vb0Vm/c2/z4uWoAD39W8XQjXc3QDADojSp8ZwM5rSxd2wQA0I0KpffyI3TtBJfXLz/+FPhs6ReNltPnX/0azp7/5uFrso1SIZOIyiLHMueEePGrtivbq1NRiWfPvHE9bzQbJwX61aNG4geOra/VDwvLuRdrsxxXK8U3tJa3ZvU2P37z/MvAqxfow7/5gtIKagWgw4kYeXz0hqnX28AyDOtegsGwTO0ErZhm85cXzb0CnHvtbKPlzNv8N3Dh3fOHrkgvFg1GGhwOs3alaj3KpfMmADjbi4mrV8eL2mbF3iwtGE0HWo/kO4dH0i5nXzsHvMI1yBKWadYai9FsgyKoJ0e2ajU4gEHuJdCvv9JEbN7mx/Pvngf4N0nJmjtbzRY2QfdOpGO+qvA9Tay9NupMowLerxZYJ2jPHnjnIvD6pYb58ciESDyqSFeU5HvRIIBqqahagH5v6EfVtG5Z1sby/E8Xp+QTb8cvKJwVXqdw8e2A+3Lb3st5hx1BB1bZPT9++60zOPvWxSaiI/UznWPmJ3KbdN9ELiF2OhQCAOCj/j4UpLAcAFAtqaVSMv4gMj2lcADDS4osvFBPHEf/sa0/fn7pP2Mf96c1vJl0C6uQSC6xA7P5pOiDX3Jq9z/LscAmcMDF82QBw/qo2SzFyGGRvr9ULedjMTmtpvlOR9QZ7JNiERu2e7Ksn7wU2GCrT0p/PD+aWeXM8fB+G/rgmMx8LF7kU8V8lO90KDjV+5/j6wVbDe8Alm3XCrgCvK+6FVCsFBYAoMLFXuACGJbna2/eVdtyP4L23qGlOfcmCifgj+dHVoqNjUnH+Q9O5j2KpQ1sLR1JGNGC6otHR5zq/c+KsoCFNaBimRbgkgCfFOoLos8qyQRZDmFtE42qHV4EFIR+ETNLAEzTbHAA6zU4nCC0+gbnj/zIiLH0C1MGZ+SikZyYK2Uat4pst9O8/wVJDGCtAhi6++Vl6gYAICTV2+r7hqHX2qtslkoG/BZc23Ci1IOlRw7MTRPYdwurfw5Qgfrzdiv54/36xWEWYkqazajZ7h+up1tI0VofTkPX3N7PDE2vAEBQifjqFuGs3x25y8i9AKCp3o3g5XcUxAElBABlba+t1jMsXTMAUAF5wKX/6AmR/NhGlhpXktZEMec+lpZeLBptjuhFwMjxSK2ljFranyCtkqoBQCgak/1R2AEAzs7y+B89UHKzmVgvgKpadAkd2F7O1sZbZnnx+q1TOtYyBWk42gM40Bf23yYcW1vQqgAVjNyItP6Zg+THdrFKSSVmxAsNamSsUjKexwtbCu8lRk4m+2igXMjtyzJWMVesAoHwRPJY5a+ecrYXbv5wc3Q2KTK8LIcAVNSF+vBCRj5d6zzj7Cwnr40UIMcS8ZgEvZD5w6tDBf8MCN06FCWOjoUDcCrqTPG5BOnY6kyxDIrpGxtTPLi/+aP8scVsy9rbjT65o9paJqJkbGXCLmYzxX1LTV0t5HW5OOGfJ5iTsIz67vfJ3oeQyE7kpPG1fDI7UUrwTz63SxMTC1UEB3O5jlcQm2q+aPJyRAqs3hm6c2FyvjbikxCRg3dnypuFmVKqL6TevLkUnY0CcMz7/6mcVIu1V5HR94bHr/zDP31wO6dHkj6rZWoB6kJsKjVz+ftLD27f0SKpr2o1HT17O1+mWHlyNuHNz25izMija+/4uHt2lkdDtR9Ht33brtanw8FDD0RPar3TcbbIeqq3/pt6UxudDuaJrfnR3gDonuHp2ii+zvbyZDgIBPsnV3dasPpDxsfdnv2Dl/D3Gx7jnbn6kM40zfDRuadGgd2ZH6y1b6GDQU5O7cW6PT81+/TGdj9/vw8UF188LNDdxx/9AU39vn/Ott3P3/8HFPWPprYO+d7jT/6kPwCKi0zWB5p+/NmHgz00WCmxeOBVvrv7k+GzVPDwXePiVOXH9bnU2OhwuOepnkp0qH/wxtjY5PxGm2LYZ3v+Rujw2xTd33j8426xMT85Njo80Pv07g/2Dd4YHUvNrrcgBZ3c1vLU6EBfTygUCoWCwZ7+wbHZ1RadnCcdP3xrdrAnQNNBKT73/DrWp8Ihmg6I0cnlxsHufvFBP0VJjQPY+143jh/+xOOHH40N9u8dwFDfwOjU8mF59UTjh5+q92shkkxHOh3E81glazjZTkfRDrySSCudDuJAnBTPFOKdjsIVF83rUfdFQrxoHBa0Y+srOt0/EeFbHZivnLs8lM4Ppdu3wVOVHwnixeRYxXsLFxLzMb7TkZwyJD8SRJdzdkq3/7Q8NpvzSXesU4S07yGIruZsL4z/6ZnUXOLFGtWnPUh+JIgutqPdG19SplL+6ax6qpD8SBDdytFzt2bYZEr5qvGmXiqdxjbinULKHwmiHX7T6nnqnK25oWiWjcWK2XrzCNvSF1b4ybmDewLRAODUx3Tzg+pvdgF8w8NZWF4CgN81858kPxJEG9j0WeBl/u+1aHXO9uLNq0P3H1Vxa+WrTylG/nDj4H4kFNjXL34Nf834Z67Fl77xEvW1v/tWg/kVWuHit1/G5tnfazC87kHI+zVBtMGWbSPwDyPvtGh11LmrWd3e1xJ6Z/HwfpJM/z/6Nn6LX/tlRCD7/zk7+PY//QPPOnhSkP7Fm/jd13/bRMEDyY8E4T1D08oQJB9UMVMQeoWXqvpSy+eKbooDfUn7NS30edlrnBOEs9BXXMe3OxjJjwThOaOQX0FvxBeTnFOsHBadzeL9UqcjAeBAu1/8OS0Neln/TjHSgPRSRX2gHjtBkvxIEB6zihOZFbrfL5OcU3x0NHx2M3c7f4wXTvs4jrxWx8zfzv5VIDLqPiJqq1BsZDT6Wjl/O3vch2aSHwnCSzur2WhsxurPZDs+hNoe6kJ0avKfOQs3r6ePMi65szUlU8eave1rF26qR1jv9vL49Zt/fmZgatKDkW2fQZ1RUlP/MlAav36zeKxSSFJ/TRAesU01m4gnC7Yypebj/nh4rKHeGJlbtIeu35Ivr0xMTSYU/oCCUYq90vLp2+yN+cmbN28vnRn6aHEq2oYbB3Xh+vTyRzevjVy7rI9NTU1EjlgU3MSYP0fXmfEfCaLtnh3fbGtxanSgNxTqHRidWtzwxdhubnY25qduhHuCwd6B0XYNAbgxP3kj3BMM9Q2OTR8+NlnLffHJ9OhAbzDYE76Rmjt81D2SHwmiBZ7Nj4vDwUDvYGp2ccP/5/72+k8+uNHPsUo7rtPdLz8aCAT7hifnljt329jZWPxwNBxixUNHy3S8Ln90gFb3GyAIP9oFnrS5lnOGllGsfFzk+KvxnOaXtobP2lqejEkcJ99eERLz63OxNnThps4MFQx1os/IRgXu0rVkm+cU216djss8KyYWmFh+dfmpqTYa8TRR/zj8EvCtPzY83QhBdN7aUAB4OfIXz328/cn74SDo0OBc+18lD/Tlww/CQbBi/EhDu+8c1+Gr3P3iJ2P9ATDiIbMjtMzu5x8NhmhGiE4fY8x4b/Pjx//6FaDnP/js1CCIltNvvQF8+9+4TBKwuzE9EATdN+mbGV+c3cc/vRECIyYPmLDhqW9/8UH/8R7TjjIPjuM4zpefvd8foFhl2vMMsfvlw7FemuJj88dLxt7WX59/9zzAv9n5TgME4a3AOxcB+9LZ/UsoPpZN54WRiWQhVvC6IcuROHru1r2fC8lCSjpKOF7UX9ecuZSYfu/+pX87fluNTnk5/bhj5sczn54b/umUcrwDQNr3EITXuGgikriWzxWtSLTzCdLRCzOfUn0fjIpH/IczYix91O8eE8XHbsjjf1yYKU3K3iVIxyzMqL/uuTV67E2Q9uEE4TlGUmS6WlK1TgcCwDHVhU/ROyD7ork6xcphEWWtZHi3DccuLZV+HZQHjp/kSX4kCO+xgsCjbBidr8h2YKzpoAWB73QkdbwonIWxZni4CUPXf42mhsAg+ZEg2oDjOcCyOp8fYZumBY7jfFMrwHEcKpbp2a5xYG6aCHBcE0UbJD8SRBswDAMcZ+gG79g2wDCdLwetocAEGcCGh7vGtm0wbDM3BJIfCaIdGHg4gcBx2ADAeDha9/ExgJfp8QRIfiQIgnDXXfnRLOUSEUngOI7jOF6UY+mC7oMSHQ+YpVxCEeNqp+N4hm0UM3FF5FmGYRiWF5VYus39w9rGl/ufqLO3FrOxq1fT3jcH6J78aKlJSbgyUkAsp5umaZqljGKkr4tiNHeqrlFLyycUnr8ycnfB8EVxVY2tZSPCG9du3VtY26xUq9VqZXNtYWb8uihEyP4n2mZreSou82/845szS0bZ+811S/twIxuJ3Fmxe1OFfFysFbRycrJQMMSr90ZksFre2xGI28HS8ul0rmQY+sqmf2bfBAAYuYh8cwG9g6NRSWBh6aVCvrBSrgLVzQcjSpzTcsfsmOBDHu//p8ev6Cz/RFJ3pFFstpZzd+7mNcPQVsrtuzw86u9Y06rxzdYn+2gAgcG551e1NdVPAwgMzPp/IKlDbK+vrm87juPsLA4HAQCB4flOB+U4juNsTYcDgb7Uc511N+Zu9O5VOfSMrXYottY52f5/dnyz/YwffAtgwj8+eZwns7u78YNvgTqjdDySut3djyMvg2K/59J1/Wnb68sb247j7D7+aIAGAIpPHO2k29390Zugvv6dySai64b3a7s4kV6pAgEltu8hhYvEJBqoPJhoQ2GEt1hBFFgAYDjBV8/Cei6jydlC8rnOunwkW8j01zLko0K+23e/x/ufPv/q13D2/DdbvuJjokC/6o9Injj72jngFe6Q6n1WkHgWANjjN2QMci+Bfv2VJmLrgvxoFbKFMgBaUqT9TZg4SREA4FEuq7Y5sBeEXiggkXYdA5+PJZRaQxHDMEhZ3UGYt/lv4MK75zsdx14k3/FBJAAACmeF1ylcfNvDFkfffusMzr51sYn/9H9+tNWCWgUAThTd7hu8KAQAoFzI+2HCytPH5mKZRpOnMILE1/4A45v+GATRKv7Pj5paqgAAeN71rYfh6p9728f9xSXG4o0nJ97Lirzkp+mnupylF3PpeETiWYbhlKfmJN3ZmEsqAsucuzJeIo/r7eD7/Gjp2iYAgG5U7LCXH6Frx53dljgpQzMAoDeq8J0N5PQw1ULJ5vsHBwSmUq2WFyYyKgDA2Zi+fmUkW3pUqVY0VTvWNKVEk3zfvsfU6ycC06gDJcMyNFAFKqZpAV3fzKSbaAW1AtDhRIw8PrYKJ8diANDPj87MfH8JZbWoQQ4Vb37vzyOqWRQsTdUgkPtRW/j/+XGvjS7TqIDrqwUWac/bVraaLWyC7jITnRUAABGWSURBVJ1Ix3xV4X46UJwc7gWAR6WSlr95h0pNxwUG4ERZEcn+bg/fPz8eC8mP7WTmJ3KbdF8ql/BocOkXnSDLIaxtojQRsYfyauNSYMIrvs+PLMcCm8ABye/JAobt2Bl07JGrGj4OdwurkEgusQOz+aTog19yGvc/BbFfCtzdrFRtKZlyadtGeM7379ccX2842vAKsGy71t0owDczAmYrmFnlzPHwXT70gZmPxYt8qpiP8p0OBad2/1OMrIg0ANuyyJtRR/j/+VGUBSysoT7CsEsCfDLysCB2qpKAlTya4M2fbC0dSRjRguqLR0ec2v3v2LpmVgFUS6pmx7yc4I9w5/v8CEESA1irAIZuuuZHUzcAACGp3la5/RjvJnjzHyMXjeTEXCnjn/KwU7n/ne3izZuaEA4+WihXSgs65FP3E/3P9+/XgBSt9WEzdM3tJcPQ9AoABJUIOX88ZxZiSprNqNnuH67H35yN+yP/rjqRzyfkAIBHqmrs/9LO+mxCETmWYTnx2njB5RvEyXRBfmTkeCQIoFpSXToNWPVJM0NR8v7hNUuNK0lrophzH0tOLxaNNkd0Wu2s371579LUdIRjpIhEA9AW1FoxkqVmsiUbgLOeHRlfk5I5dXE+pdiLd/7wWteP0eI7XZAfwcjJZB8NlAu5fQnSKuaKVSAQnkgeq/yJOC6rlFRiRrzQoEbGKiXjeZBmec2z9EIuV9Qs7GzkR/5oJTZbq7Jm5YgIoKrOFEzsrGZG7jKyxADWwhIzkU9HZVG4fDU+PR2/iEcL2ukcTb9z/F/+CABCIjuRk8bX8snsRCnBP/ncLk1MLFQRHMzlTlMLZWtvmmTfVFraWiaiZGxlwi5mM8V9S01dLeR1uThxOp7gO7L/9Wz0+vgaMELT7OWJxcUnAyZxcqQXK2u2OiJw4xeuTS1OCwDAyrHYkzZKFLgQR4fqA7SdbvZXs+R6f3y8zo9HGhn4CMRksWAq0btJJcYWMjGBBaxSJhq5+yjYP9noha9L6Wqp1uW8omsGfNCRTM9FlFtLFeD++K37jb7UkyqckiZ6Te7/E47KzcvR/qBest+4NjE9nXi6YYAQT9/IR3PGuctDk9OTT071pxtw7mj3//xMYva90/8O5djaUqnWns/UdQuSt3cEb/Pj56oB/L+fVRA7+a/glIymR7PpdDYiTtiAbTOCHJtdTUZdhz3rQkYxky1qulp4sFb/ZG1CkrSIxPNSLBEVOpN9rGJcGVk4bKoPuj8R7fYu2Cfa/7/8+FPgd0u/wFOvN8fCSEnVTLouYpWsZmUb/ePOxvztkZF7letTTc3w3DX0Qjqn6lqxsFC7ezl2MS7KRUXiBTme8Ogxwtv8+M3zLwOvXmjZxL+cFM8U4q1am9/wSiKtdDqIfVglazgNL87T5ET7n+ZeAX772tlWBnQElpbP5oq6DUvLfu+yZq4uJ7r9NtWIEEmmI+3eqLf1M+ffPQ/wb57q2xpBAAi8cxF4/VK78yMrRpOZXKG08fD9/rP2yp27apsDOOW6of6aIIiDnbmcmBrrRdkwSA12K3VH/TVBEIfgBYEOUBx5WWsl8vxIEKeCuWmwkeFT0obAL0h+JIgu5GwVEtFYulCfNnJnY2b8weXpFOn22Vrk/ZoguhFDw1TTUSErRaOKwHFien6KJw+PLUbyI0F0IerctYxqZDodxmlH3q8JgiDckfxIEAThjuRHgiAIdyQ/EgTRWa0axab1SH4kiLY44fg+LeSfSAA4+MXHnzn49C9+6d02/vJnFfztx80MHkzyI0G0gf3XP/8NzP/1t52OYy+Sz3wQSd3Z184Br3AtG8VmvyD3EujXX2niP0n7HoJoA+Zt/htg3j3f6TjqkZS/44NIAAAUzgqvUzjzdsC7bXz7rTP4m7cuNvGf5PmRIAjCHcmPBEEQ7kh+JAiCcEfyI0EQhDuSHwmCINyR/EgQBOGO5EeCIAh3JD8SBEG4I/mRIAjCHcmPBEEQ7jzOj7/5HfCb33i7DYLwgd/8FvjtQaf67xy/DFLzO78Nl7PrdUC/c+D8tpl/9DY/fv7ffw6srVRatT6zlEtEJIHjOI7jeFGOpQt6N8332+3xH51ZyiUUMa52Og5XtlFMR6VI1mjhOn+5/D+B9aW/abS88snqr/HX6uct3GRTHFQ+Xq3ir5c6Hkmdg198vH7U8Xu2V2eTUXnvAhKkaDKvHeEC+suf/Qr/95PTPX6PpSYl4cpIAbGcbpqmaZYyipG+LorRnG53Orgj6Pb4j8rS8gmF56+M3F0wbL/9MNsoZqIi98a18fsrpt+CIw60szopX/ruv7pXVrKaaZqmuTo7bOeHvnvpWqbk3UOG46Ufh18CvvXHxsnXtDHVHwDo3tTqzlOfbi/eCAEIDs5tnXwTnur2+I9ie3V2bDDc39cTrJ9bgeH5Tsf0xM7G4uSNgf7+vtDeOFp9kxstXP/aUAD4O//8vzVabvzgWwAT/nELN9mU3d2N+LdAnVE6Hknd7u7HkZdBsd9bO/BbX/zZYBCg+PjiUxfQ7pcPx3oBilWmNg7awo94UF//zmQT0XXH86OeiSaWKghEJuLi01NYsvJEsp9G+X4snvfze2q3x38klgkxmS+qJS0/HDz8621mWUwkU1DVkprp63QsxDE5W7mRkftligm/95781AVEnRHH3hsIONbCrZGWlpY80Q350S5OpFeqQECJ7Zv+nIvEJBqoPJhIN1O80BbdHv8RsYIosADAcALX6WD24USpNjs0x3k5EivReo5TunN7oQKKlgcjz51ZFKsMKwE49tLt20UPCky8zY/fPP8ycP7Cyc5Hq5AtlAHQkiLtn/+ckxQBAB7lsuqJNuOZbo//lNl/CFqC5l4BAq9905u1v8gcW72X3wQoCGH5+ecLUIwUFgE45cKMBwnS2/x4/t3zQOjNE52StlpQqwDAieK+vQOAF4UAAJQL+dJJtuOVbo+fOJLAOxeB19452+k4Th0H2v1iGQACosC7fIETxBAAp1K83/oE6f/3a00t1doH8bzrWxvD1T8vayWjbVEdXbfHTxAdpaulMgAKfA/vspgC38MBgFPRVvRWb9z3+dHStU0AAM1xbo9fwF5+ga61fPecXLfHT7SbbZQKmURUFjmWOSfEi19V3G2vTkUlnj3zxvW80bn42sux9bX6ZcFy7sXa7F6JsqG1vDWx7/OjqZu1PxiGdX9PZ1imtnsqpum/SuBuj59oM71YNBhpcDjM2pWq9SiXzpsA4GwvJq5eHS9qmxV7s7RgdDjK9rFMswoAoFmmwQW0d2VVrZZfQL7Pj9ZeG2Omwd55aoHlu/bI3R8/0WZCJB5VpCtK8r1oEEC1VFQtQL839KNqWrcsa2N5/qeLU3Knw2wb297LeYddQQ6s8guXH4+l2/NLt8dPtAwFKSwHAFRLaqmUjD+ITE8pHMDwkiILHlXC+9PRf2zrrx/fz3/NciywCRzw458sYFj3Er5O6pL47eOeWg1v5kRrUIwcFun7S9VyPhaT02qa73REncE+KZayYbsny/rJS4ENtvqk9P3zI8fXC18bXsGWbdcKKAJ8gxqQTuqK+M2scuZ4eJ8OPXGKUKwUFgCgwsXSMf+1uG8TlucDAICqbblfQfbepUVz7k1ETsD/z4+iLGBhDahYpgW4JJAnhbKCKLQ3tqPoivhZKTY2Jh3nPziZ9ygW4glBlkNY20Trqx26BwWhX8TMEgDTNBtcQPUaHE4QWv2A4fv8CEESA1irAIbuvntM3QAAhCSJb2tkR9MN8TNiLC12aNtEQ4Zea6+yWSoZ8OXJ3Q6cKPVg6ZEDc9ME9j1C1D8HqED9ebuVfP9+DUhRJQAAhq65PV8bml4BgKAS8ecl3u3xEx3hrN8ducvIvQCgqR6O4OVzFMQBJQQAZW2vrdwzLF0zAFABecCl/+4JdUF+ZOR4pNbSQS3tTzBWSdUAIBSNyf6sMej2+In2c3aWx//ogZKbzcR6AVTVosup88z3F+OX3riptie69qIgDUd7AAf6wv7bhGNrC1oVoIKRG5HWF993QX4EIyeTfTRQLuT2nSVWMVesAoHwRPJY5Wft1O3xE23mbC/c/OHm6GxSZHhZDgGoqAv14Z2MfHpf5xlnW701cu9RtzYOO/vaOeCVxqMqUZQ4OhYOwKmoM8XnEqRjqzPFMiimb2xMafh8EXzt66Bff6WJ2LohPwJCIjvRS6OcTz47yptdmphYqCI4mMv5uoKv2+M/Jsuon8V+vGIty6qV5rc3OPr8+a/h7PnG4/uYaj6XLxk2tpfvDN25MDkd5QAKQkQOAtgszJRsZ6sYv7kkKfwz/+hsL92+ZwpHHnKTfvU8DoykvSicFV6ncPHtwAHfuRCbSvUHHOvB7TvPFFI5evZ2vkyx/ZOziQPKHr/9xu/h7FsXmwmviTF1j2491QOEp7dbsa6t+dHeAOie4en12vq2lyfDQSDYP/nMmNx+1e3xH916qrd+cvWmNjodzPO25wbrmSQwON/C/b411Q/0ptYbbnc2TB+wfGduoPb8RNMMH316OPmd+cFa5qCDQU5OPX+u7D6eHx2e+uyT0RDFxRePEOju4w/DFCU2jrTddj9/v4+i5KnDhtB//Mmf9AdAcZHJ5doF9PizDwd7aLBSYvHADLO7+5NB+oh753nd8fwIAJyS0fTljGJlIyLP8zwnxApsbHZVVxNiNxTcdXv8hzOKmWQiFhGl8bX6J2sTkhSNJ5LpfMen2LG1fDoRj8pi9H659knlflRUYolkOqu6Ffu3FyNFIz0Bmg5eHpldnH16FFhGmUiHQzQduHQ1OTeXfPZccbYXbs8IyRs+bNjWcueujKsbD//LMF2I1S6gS9fu2UpmeX15cv+wkK3i//Y9T+OkeKYQ73QUzev2+A/GK4m00ukgGmHEaDIT7XAQjVo4A+Cieb1BeEK8aLieNM72wviM8N6sAGf5WGEc58ttcrSYzl0eSueH0t6s3E135UeC6FIM0OL+wc72g/EZ8b1ZAXCO838+7eTvYVh2LUE285ZG8iNBtAPDALbVsmaMzvaD8fu17HhMlmUBZ3zUfZ5hGAfllg+98xXbsqto7geT/EgQbcCwHAvLbNDH9Lic7Qe37mIgZWuaBsBxdNOGvWVomsZwgsA1zgUOrE0L9AUfDVXAciyNNdO7QmDTNClwoWZaiJD8SBDtIIgCqoZmIHryblIOjFJJvTdz5d4zHz8YufyAEpIP11MHbMLSDRPc9Zb3VG4eIwg8/lzXDTzXcqk1HNtYMxC43tSkmt1Tf00Q3YyTlV48KhZbMYcGRV1O68+0YHm4175n98DkCMdSF0q/Dsqt76ncNApCWH4N2oI3zQgcu/Sg9CtaCjfV+ZDkR4JoCyEa66PXshm1g/Ujjp69U/hVKHrDT31ZKcg3Yu/YC3fvejAFvGPm7+R/HoyM7pt6/khIfiSI9uDjmUSPeS8WK3SmvaWzvXhr6Pb/ODc49V7rB3I4CYoSxyaHL2qZoZtqa2tpdlbT12/9VyqcSjWXHkl+JIh2YaR0MT/MFK5LSqalL5MUdTlj7G4dNCvN9ur0yJVr/7F89YPF6SZzhYeoc9em5j8M70xfuzKU1VqTI3c25uJXrv77R5f/ZH42xje7lib63BzdeqoHAN1YYHDe0wAIomV25gYDjU9l4MD+hU+tZn1ubKAnEAiFb0zOLW943Ld0Z2t1fmp0oDcQCIVHp1db0tXXM48ffjQaDgUCPQOjU/OrW03umu31xemxwb5gICjFppYP67Z4MMpxjtW49HjMYnqiYBzwBUZKZGL+KSomiANouUT2wHHGhGg6cdS+bpauFtXSCuRU3MO3XUfPj+dMQZYVWTyg1Y+v2Kamqqq6wkYnY03U9RuFdNbgJElWJP7Ev9jb/EgQBNG9SPkjQRCEO5IfCYIg3JH8SBAE4Y7kR4IgCHckPxIEQbgj+ZEgCMIdyY8EQRDuSH4kCIJwR/IjQRCEO5IfCYIg3JH8SBAE4Y7kR4IgCHckPxIEQbgj+ZEgCMLd/weOw3E83BR9yQAAAABJRU5ErkJggg==)"
],
"id": "b3caecd8"
},
{
"cell_type": "code",
"metadata": {
"id": "3ebaa2dc"
},
"source": [
"# Main Program\n",
"#step 1: Create Coefficient and LHS matrices store (a,b,c,d) seperately\n",
"a = [0,-1,-1,-1]\n",
"b = [2,2,2,1]\n",
"c = [-1,-1,-1,0]\n",
"d = [0,0,1,0]\n",
"\n",
"# Call thomas function\n",
"x = thomas(a,b,c,d)\n"
],
"id": "3ebaa2dc",
"execution_count": 2,
"outputs": []
},
{
"cell_type": "code",
"metadata": {
"id": "4e06a951",
"outputId": "784b1f12-4062-49ba-8319-e23c24821203"
},
"source": [
"# Print x\n",
"print('The values of x are: ',x)"
],
"id": "4e06a951",
"execution_count": null,
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"The values of x are: [1. 2. 3. 3.]\n"
]
}
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "c65d59bf"
},
"source": [
"<h3> <font color='blue'> Closing Remarks </font> </h3>"
],
"id": "c65d59bf"
},
{
"cell_type": "markdown",
"metadata": {
"id": "b8333d71"
},
"source": [
"Tridiagnonal matrix algorithm (TMDA) or the Thomas Algorithm is an efficient way to solve systems of equations whose coefficient matrix is tridiagonal. This approach comes in handy when solving $2^{nd}$ Order ODEs (Boundary Value Problems) as well as 1D parabolic partial differential equations such those depicting the flow of heat, mass and momentum in natural and engineered systems. Understanding this algorithm is therefore beneficial to engineers. "
],
"id": "b8333d71"
}
]
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment