Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
TB Lecture 27
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Lecture 27: Power/inverse iteration\n",
"\n",
"## Rayleigh quotient\n",
"\n",
"Suppose $A$ is $m\\times m$. We can define a function called the **Rayleigh quotient**,\n",
"\n",
"$$ r(v) = \\frac{v^*Av}{v^*v}, $$\n",
"\n",
"for all $v\\in\\mathbb{C}^m$. As shown in the text, the stationary points of $r$ are the eigenvectors of $A$. More specifically, if $v$ is close to an eigenvector $x$ with eigenvalue $\\lambda$, then \n",
"\n",
"$$ |r(v)-\\lambda| = O\\bigl( \\|v-x\\|_2 \\bigr) $$\n",
"\n",
"in general, and \n",
"\n",
"$$ |r(v)-\\lambda| = O\\bigl( \\|v-x\\|_2^2 \\bigr) $$\n",
"\n",
"if $A$ is hermitian. Thus we can turn an eigenvector estimate into an eigenvalue estimate. \n",
"\n",
"For example, let"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"A =\n",
"\n",
" 2 1 1\n",
" 1 3 1\n",
" 1 1 4\n"
]
}
],
"source": [
"A = [2 1 1; 1 3 1; 1 1 4]\n",
"[V,D] = eig(A);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We'll color the unit sphere according to the value of the Rayleigh quotient at each point on it."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAjAAAAGkCAIAAACgjIjwAAAACXBIWXMAABcSAAAXEgFnn9JSAAAA\nB3RJTUUH4AoUFDQeFVGscwAAACR0RVh0U29mdHdhcmUATUFUTEFCLCBUaGUgTWF0aFdvcmtzLCBJ\nbmMuPFjdGAAAACJ0RVh0Q3JlYXRpb24gVGltZQAyMC1PY3QtMjAxNiAxNjo1MjozMAvpoZQAACAA\nSURBVHic7L15mBxV2ff/PaequmfLZIeECSTshF0UvNhcohgBgRd48IkLeqG4gEtAfXyQRXkFHn2J\nPj8XNIqAEkBW2QQiyhqIEJYYQliC2WYyk0xmJsvs3V1V5/z+ON1nqmvrquqepWfO5+orqa6urU9X\n93fu+3zPfQjnHAqFQqFQjDZ0tC9AoVAoFApACZJCoVAoxghKkBQKhUIxJlCCpFAoFIoxgRIkhUKh\nUIwJlCApFAqFYkygBEmhUCgUYwJ9tC8gEuvWrXvuued8X6qpqZk/f/4xxxzT0NAwTGf/29/+tmHD\nBgBHHnnkRz7ykRHYEcBDDz3U1tYG4IQTTjjhhBNi7Ts2Kac1xjLj9X2VTzXew2Pq06zGBiwXXg38\n9re/DX8XU6dO/eMf/zhMZ7/gggvEWb761a+OzI6c81NPPVXse80118Tdd9T53e9+t2TJkiVLlrz5\n5ptyZTmtUeaph5WReV9jHN9mr8Z72PfTHPmbSlCNDVgm1REhlWT37t0XXXRRTU3NokWLRvtaFLjh\nhhu2bt0KYPr06UceeeQEOfVEZnw3+/h+d2OK6hOkH/7wh3PnzhXL27Zte+KJJ1566SXxdPHixWNK\nkK655pqvfvWrAJqamkb7WkYf1RqKsY+6S0eX6hOkc84557jjjpNPr7766hNPPPHll18G0NHRsWnT\npgMOOMC5PWOso6Njx44d6XR61qxZU6ZMqez12La9cePG2bNnT5o0yfXSUUcdddRRRwXt2NbWRgjZ\nZ599Ip6oq6uru7t7//33pzSGFaWjo2NwcFBKeKVobm7mnM+dO5cQEnGX8NYIaUYvHR0dvb29cZsi\nMcN9CwmitEB3d/f27dsty5o9e/b06dPjniLWpzYcjRzlHo51J5SD7xsMv0ujE+uTamlpSafTe++9\nd8nDVrwB29ra+vv7DznkkMRHqDCjnTOMhLMP6fXXX3e9euONN8pX169fL9d3dnZedtlle+21l3yV\nEPKBD3zg/vvv55ybpnnkkUc2NTU1NTX97Gc/cx7w97//vVh/3HHHMcaCOglefvnlU089ta6uDgCl\n9Kijjrr++usZY3KDr3/96+I4//3f/y1XMsauu+66OXPmiGPOmTPnvvvuu/XWW8WW8hTO9PEbb7xx\nzDHHiB+R2traiy66qK+vr2Sj3XDDDQceeKA4yIwZM2688cYVK1aIs3z0ox8V25x//vm+LXDiiSeK\n9XfddZdz/bvvvnv66adPnjxZHHbSpEmnnXaaM7F+8sknNzU1aZomNpg6dWpTU9M555wT1BpRmlFe\n5MMPP/zkk08efPDB4uCupgg5tYtFixaJA1588cWuFhPrL7roIrEm/BYSeG+PuK1asgU453/5y19O\nPPFE5zf38MMPX7p0qffdeSn5qUVvZC8hzR73Ho7SDi6iN3XEN+i6S6PfVJKSn5TzSm6++eZZs2aJ\nzfbff/8rr7zSsiy5ZcUb0Hnq1atXH3/88QBmz55dzkdQWcaDIMk03d577y1Xbt++ff78+c7fEect\nsmzZMs751772NfH06KOPdh7wwx/+sFh/+eWX84B+zp/+9KfyNnXyH//xHwMDA2Ib3x0//elPu3bR\ndX3BggVi+fzzzxebyXvxjDPOmDp1qmuX8847L6S5GGOf+cxnvNcmj3n44Ye7zuLqNd1///3F+ptv\nvlmuvP3222tqaryHNQxDft+kBDo58cQTy2lGeZGf/vSnU6lUUFOEnNrF0qVLxauTJk3KZrNy/RFH\nHCHW/+lPf4p4C/m+r1itGqUF7rjjDnl212VcddVVIXdCxE8teiN7CWn2WPdwlHbwEr2pI75B16cZ\n/aYSRPmk5JXIBSenn356JpNxbVmpBpQH/PGPfyyjHylIyT6CylLdgrRt27Y//OEP9fX14qWvfe1r\n8qWrr75arJwxY8azzz6bzWa3bdsmlea0007jnMvOJzhCq46ODvmprFmzhvv94qxcuVLcbel0+qqr\nrnrooYcuv/zymTNnyg9bbObd8a9//as84wUXXPDAAw8sXbrU+avnFSQAkyZNuvTSS7/97W9Lazul\ntKWlJai5HnjgAbnvueeee/fddy9duvTwww+XKxMIUnNzs2znBQsW/OlPf1q2bNnChQvFmnQ6/e67\n73LOn3322ccee2zGjBli/eLFix977LGVK1eW04wRmyLk1C46Ojp0PZ+s/tvf/iZWCrMvgFQqtWfP\nnoi3kO/7it6qEVvgAx/4gFhz1VVXZbPZzs7OJUuWyKvt6ekJuhMifmrl3G8hzR79mBHbwUsCQQq/\nGNenGf2miv5JOa/koIMO+tWvfnX//fd/9rOflSvlHwoVb0B5QHlXzJ079/jjjy/nI6gs1SdIxAEc\nfP/737dtW+7y1a9+9eSTTz755JNvuukmufKmm24SG8tY6tBDDxVrbrjhBrHmlltuEWuOPfZYscb7\ni/O+971PrHHmXh555BGxcubMmeIPipCfqrPPPlvu2NXVJbslvII0ffp08fvIOX/22Wfl+12+fHlQ\nc8k+tvPPP1+G2zt27JC3VwJB+vznPy/WnHTSSblcTqy0LOu0007zvqN9991XrLztttvkysTNGKsp\nfE/t5ROf+ITY7JJLLhFrfv7zn4s1Z511llgT8RYqR5CitIBlWVI+f/GLX4htGGNXXnnlpZdeeuml\nl77zzjtBbzP6p1bO/cYDmj36MSPeCV4SCFL4xfjG8RFvqoiflLySxsbGbdu2yd2/9KUvifXz5s0T\nP2UVb0Cnwp111lmdnZ1y48QfQWWpPlMD95tRcNasWQsWLHB29P3+9793btDR0fHaa6/97ne/E08Z\nY2Lhi1/84pVXXgng/vvvFwsPPvigeOkLX/iC7wV0d3f/61//AkAp3XfffYWfAsD06dOnTp26e/fu\nzs7Op59++lOf+pR33zVr1ogFmS0UO55//vm33nqr7+nOOOMMmf0//vjjNU2zbRvAjh07fLfv7+9f\nvXq1WL788sulbO+1116f/exnf/nLX/ruVRIZTV566aWGYYhlTdO++c1v/uMf/3BuEJEEzRi3KUJY\ntGjR3//+dwCPPvrob37zG0KI/O7953/+p1iIeAslJnoLHHzwwe+88w6Ayy677OGHH164cOFJJ530\nwx/+MJ1Oh58iwadWwUaOeMxyvlAVv5hy0DQt1id1wQUXzJ49Wz797ne/e9tttwHYsmVLW1ubVMGS\n15ygAWfMmHHvvffW1taKpyP8EYRQfYJ0wQUXiL/0OeevvvqqyOC1t7efccYZy5cvl3/5Ali7du0f\n/vCHl19+ef369b29vb5Hu/DCC6+++mrG2Jo1azZu3Dhz5synn34agK7rn/vc53x3WbdunVhgjH3o\nQx/y3Wb79u3elZ2dnfIyXLY331S1wHnL1tfXz5kzp7m5GUAul/PdftOmTXLZNWbi6KOPDjpLOAMD\nA5s3bxbLzgQjAJkJ7Ozs7OzslEFYSRI0Y9ymCOHcc8/9+te/nsvl2traXn/99f3333/lypUAampq\nzj77bLlZlFsoMdFb4A9/+MOZZ57Z3d0N4LnnnhNVS+rq6s4+++xrr71WRvkukn1qFWzkiMdM/IUa\njospk1iflOtDOeywwwzDME0TwKZNm5yCVPEGPP3006UaJTvCMFF9gnTFFVc4bd/vvPPOUUcdZds2\nY+zGG2+UgnTLLbd85StfEct1dXUnnHDCEUccUVtb6yr6MGfOnI997GPir8X7779/3rx52WwWwOmn\nn+70VjkRGwCglB500EG+2/j++dzY2EgIERGeuGUl/f39Qe/X9bdVSQ+u83slEwgCb19uEOJbIbEs\nS74jV7en86lrr3ASNGPcpghhypQpCxcuFF16jzzyyEEHHST+5Dz99NNlZ2/EWyg6rvaJ3gInn3zy\nli1bbr/99uXLl69YsWJwcBDAwMDAPffc8+STT7700ku+mpTsU6tgI0c8ZuIvVAght+JwvEFJrE/K\n9aFQSjVNE1fu+kGoeAO6RlkNx0eQjOoTJBfz58//1Kc+JfItMtLM5XLf+ta3xPIVV1xx7bXXik/0\njjvuEL8mzv6nL37xi0KQHnjgATmG6Ytf/GLQGaUXi3P++uuvR6+hJ0axiD80nnrqqZNOOkm+JC6g\nIhx22GFS9jZs2HDMMcfIlzZu3OjaWKZxZJc+AM55R0eHc7PGxsZ9991XDFZ3HVPuOGXKlOhjqlBG\nM1aKRYsWOQVJrJT5uli3kIuIrRqrBaZMmbJ48eLFixdns9mVK1c+8MADN998s23bu3fvvv322//n\nf/7Hu8twfGrDQTl3QsSmHkmif1LOawawbdu2TCYjlkNSJl4SNKBst8RHGCbGQ7Vv6Rrq7+8Xkce6\ndevkR/uDH/xA/n0hI1NnR9S5557b2NgI4PXXX3/00UcBTJs27ayzzgo63d577y1GEXHOxfYCxtji\nxYvPO++88847r7293XdfmQ76+c9//uKLLwKwbftHP/rRK6+8Ev99+1NfXy/zgdKgASCbzd55552u\njWUU+Oqrr8qV//jHP7zpi/e///1iwdXXJU8hN3AS8rtQTjNGoeRP0tlnny2yFm+++ebjjz8OoK6u\nTmbJY91CLiK2asQWeOONNw488MADDzxQGPzS6fSCBQt++9vfnnPOOWL7nTt3Bl1Jsk+tHBIoQTl3\nQvQbuCKEv7u4n9S9997b09Mjn4oOJACNjY3SlBGF8r9Kw/1ljE7VR0gAnH/i7dmzZ/LkydKmCeDh\nhx/+whe+YNv2P/7xDzkAReRnBHV1dRdccIH4xorQddGiReHZrZ/+9KfCv3TJJZesWLHiggsu2LNn\nz7Jly8Rn+dGPflQOdnPx/e9//9Zbb7Usq6en59RTT50zZ87u3bv7+/tTqVQFv0Jf+tKXfvjDHwK4\n6aabampqzjvvvF27dv3v//6v6y8yADKBsGHDhgULFlx00UXPPPPMvffe6z3m9ddf//jjj5umuXz5\n8gsvvPCiiy6ilN5555333XcfAErpT37yE7mxzHotXbq0r6/voIMO8o04EzdjCBFPDaChoeHMM88U\nFnnR+Geeeaa0w8a6hVxEb9UoLTBjxoyurq6enp5NmzYtWbLksssuMwxj/fr1q1atEgc5+eSTg64k\n1qdWDtGb3ZfEd0L0pi6HiO/uiCOOiPVJdXV1LViw4Kqrrmpqalq+fPmPf/xjsX7x4sXRs+uC8r9K\nw/FlTMIIOPnKJ3xg7JtvvilffeGFF8RKOSAAwMyZM8UPjQxUdV13HuGFF15wtsmqVaucr0Yc3yrY\na6+93n777ZAdn376aVcpkY9//OM33HCDWPbavqMMrnSRy+V8S0DKrmxp+966daurnwlAU1OTDOGd\nZ1myZElQwv1HP/qR8wKkgVUQMjA2SjPGagrfUwfhHLAFwFl/gUe+hbzvK1arRmkB6TUHUFdXt88+\n+8iE4fve977u7u6Q9xjxUyvnfgtq9ljHjNIOXqI3dcSL8b1Lo99UUT4peSUHH3ywN/F7xBFH7N69\nO8GHUs5XqZyPoLKMh5Td/PnzpWNEeHkBPPDAAx/96EfFcmdn58DAwOc//3n5p4plWWvXrpVHOOWU\nU2TsP3/+/ChTj9x7771//OMfnX2DlNJzzz13xYoVLvOMiwULFqxbt+7OO++87LLLvv/97z/44IOP\nP/64/HPbaX1JjGEYzz33nNO23tDQcNtttzm95oI5c+YsX77c+bfPMcccs3LlSt/ad9/73vdWrFhx\n7LHHOseiH3nkkc8888y1117r3PLnP//5ueeeG+WvvMTNGET0UwM444wzZLpcBEzOV2PdQk5itWqU\nFvjGN75x9913i5EiAwMDYvBKfX39l770pb///e8i4RxE9E+tHGI1uy/J7oRYTZ2Y6O8u1ie1aNGi\nu+66S45BTKfTF1xwwUsvvZSsWGL5X6WKfxkTkO/9Hq80Nzdv2rRpypQpBx98cEhPXTabnTlzpvD1\n/uQnP7niiiuin6Krq+utt96qr68/4IADpk2bFr7xP//5TzHj1rx580QhKcGHPvQhEaVdeeWVMloq\nn927d69bt27mzJmHHHIIpfRXv/rV4sWLARx++OFvvfWW3Ixz3tLSsnnz5vnz50ep8Njf3//2228z\nxo444oiQVrUsa9euXbZtT5s2reSImVjNWJJYpw4n4i3kJW6rRmmB1tbW1tbWvr6+vffe+8ADDxQ1\nxyIS8VMrh4o0e4I7IW5TJyPWuwv5pOSX/ZprrhEVEN59993e3t5jjz22HEWXlP9VquyXMR4jE4iN\nWS655JL/+q//kgmuVCrV1tY2fKf76U9/Kk7U0NDw3nvviZWvvPKKSDsQQtauXTt8Z5ejYmXKTqFQ\njDATcNq96IwHU0M5CEemfHrJJZcMqwv2c5/73C9/+cvt27f39fWJmdcppW+88YZlWQAuvvjiipS+\nVygUimpkPPQhlQ+l9Nhjj/3BD37ws5/9bFhPNGfOnBUrVpxyyikAbNtevXr1a6+9ZppmXV3d0qVL\nb7755mE9u0KhUIxlxnkfUkk6Ojq6u7v33nvv8J7hivP222+vXr16+/btDQ0Nc+fO/chHPhKrPyAZ\nu3btEsWv0um0axpDhUIxMjQ3Nw8MDACYMWNG9FJbE4SJLkgKhUKhGCOolJ1CoVAoxgRKkBQKhUIx\nJpjoLjvFCJPJZNasWSMKZ+nFAPCdbFuhUEwQlCApRg7Lstrb23t7e2fNmmVZljC7ZzIZy4FQJpdE\nOUVLoVCMV9Q3XDFCWJbV2to6Y8aM1tZWp7q4qgYIlZL6JOQKgKi97VQmb4ClUCiqGvU1VowQXV1d\nU6ZMKVm0RoqN9yWnViE0tPLqVsXfjkKhqDjqi6oYCdrb23VdT1Y1UuLSKm9opdKACkVVo76NimFH\nTO3lnGFoOEiQBgwPrZTDQqEYYZQgKYYXoUYjNLtXABVJA0KFVgrFMKO+VIphZM+ePZZlCZP32KRk\nGhDKYaFQjBTqm6MYLvr6+vr6+sayGpVEOSwUipFEfTcUw0Imk+nq6hrdTN2wohwWCkXFUV8AReXJ\nZDLt7e2zZs2asL4A5bBQKBKgBElRYSzLmuBqFI5yWCgUQaj7WFFJrEI5BqVGCaiIw6Knp2fatGnK\nYaGoRtTNqqgkEcsxKBIQMbTasWOHVYxyWCiqBXU7KipGRcoxKBLg1Kr6+nrXGGTlsFBUC+qeU1SG\nkSnHoEiAclgoqgUlSIoKMBbKMSgSoBwWijGFunUU5TL2yzEoElARhwX8AiyFIgh1fyjKYhyUY1Ak\noCKhlXJYKFyoO0CRnHFfjkGRgJKhlUoDKoJQH7MiIaocgyIBymGhCEEJkiIJlirHoKg0ymGhUJ+W\nIjaWKsegGFmUw2KCoD4SRWxUOQbFmEI5LMYNqtEV8VDlGBRVhHJYVBeqZRUxUOUYFOOJZA4L27bn\nzp078lc7EVCCpIiKKsegmDiEpAG3bNky0lczYaCjfQGK6kCUY1BqpFAohg8lSIrSqHIMCoViBFCC\npCiBKMeg+o0UCsVwowRJEYYqx6BQKEYMJUiKQFQ5BoVCMZIoQVL4o8oxKBSKEUYJksIfVY5BoVCM\nMEqQFD6ocgwKhWLkUYKkcKPKMSgUilFBCZKiCFWOQaFQjBZKkBRDqHIMCoViFFGCpMijyjEoFIrR\nRQmSAlDlGBQKxRhACZICmUzmvffeUwNgFQrF6KIEaaIjyjGk02mlRgqFYnRRgjShkeUY0un0aF+L\nQqGY6ChBmtCocgwKhWLsoARp4qLKMSgUijGFmsJ8gqLKMZRPc3+WEN4yMLh1YKB1sF8nbOvAYGdm\nj821NDU7s91ZO5XWchScAIRwcJJjRq2WYSBTUtN1Yu9dO8lk2j61DfvU1jPQObV1Bzfoe9U0jvY7\nUyhGByVIExFVjiEuWwcGNveZK3fu1Kn1ctfOHYN9u82dGrENwjQwg9g6YRphGpgGphGugetgaXDK\nOCGcghNwAtQDYACQM1uywKY+MJDNIIwTBmKDMpAes352Q02Dsdfsurpjp+6zT239PjWT5tarTj7F\n+EcJ0oRDlGNQA2DDae7Pruza2TI4+MquHa92ddXomRSxjaEHm2XYOphOmI78ghAhjXANTAOn4l/C\nCTgFp4DQJAIuTkEADvGEcICDMBAGYmt7bJNaZnPrANncqZugNid1NTNqjb1n1k46buqswxr0w6eo\nj08xDlGCNLFQ5RhCaO7P3tncvmp3+5rd2zUwg1pChGanrRS1DWLrYAaxDTAdtkFYQZDyD40wDVwD\n18AoOAXXCjoksnb53F1Bk4jj1EKQeGGBEWIDDMQm1AKxQU2zJ5drbu/XHuzUB7leXzPjkClzD5+8\n7141k0+ZMXW0WkyhqCxKkCYQohyDytQ5aR7I3LF5R9vgwEPbNhrUThErRc0pui2kKB8VQfxbUKO8\nJuXVSAPTkY+KtLwI5f+lAHWER84gSaoRAQAugiQUYiYOwkEYkI+ZQCwQmxAL1CTUBM2Ze7Z3bNnc\noWe4cVt65gGT958/Zb8PTJ49d1JqlBpSoagASpAmCplMRs1HLhHB0H2tzdszvWlqpYg11TANwlLE\nKuiQZRBmwE651YgJNdLygVFehPR8mm5Ik5yCRADKuVONhA6R4jgJ+QxecbQEwkhenGxQIU4mqEVo\nDjRHaNbs2dHZsrnT+AvXe+mBJzVM/1DDpJNnTB/xRlUoykUJ0oRAlGNQaiR06P+9vTFtZFLUSlFr\nii7yclaKFqIiWEJ7UvnYiBX+HcrO6QUFErGRDk7B9SEd4pQ7ln2ydkNSRDj3vVROhjSJ83yoJAMm\nu5DNs5CXpRxohugZrGnuTS1du+oP6b33mXTI/5mz/wnT9hrJFlYoykEJ0vhHlmOYsGrU3J+9Y8uO\nu7du3Z7pNag1ucZKEcugdorYhZDIMvJqZKeGAiOmw07JBB2YDp4CE6JSC0uISg0sIS0UPMWZ6B/S\nwFLCTgcQIMVtgBRHSCJsAjz9SQKhUxyEk3zAxAr/MkJsEMYLqTxQi9C8MoFmSCZDtAGzt2dn801d\naZqe9cG9jzlp2sxjps0eoeZWKJKiBGn8M5HLMVy3buuLnbtf2t1haKZBrUm6bRDLoPm8nCNBZ6dg\np4mtgdUQS2TkaolFgRqYQj9qYROgHpY4cl1hoZZb8nS1cCxz23kltdx0XZtI9IllDiLddw6IECxH\nxxIYCBcecVKIlgixORHRkuhkykHLgmaINkgG+82+tVu3vtxiWKl9zp176n/su1/ZjapQDBdKkMY5\nE7McQ3Nf7o4tO37y7gadWrpm1RkyLzekQGli6SQfANURk4DXk5yIVBpIDkA9TDgUKESKgnTIqVVp\n2N4wyLmBCz60QCAFiYg4iTBAyBIDbE5sUJsQG8Qi1OLEBDULypSBNkiyA0Tvswb+tmHr37fOev9e\n7zt9n0P2q6uN06IKxUigBGk8MwHLMTT35a5f13JXa6tOrRrdMqhlUCufkZO9RMQSHu40seqIScFR\nECGnFKEgQlGkSOpQDbedI43SGNInCp4uDpvkLuFvasjjwGXf0pAHj4lQief/tUBtktckCyQHLUe0\nLLRBaNNJpi/Xv3Zr62Pvznz/3KPmT9nvP/dVAwAUYwglSOOWiVaO4bo3tt+5tW1rpk+jVo1u6dQW\nHUUGsVLE1omdJhYlvJbmgyEUi1AsKZI6RBxy4gx3auAfJyFYfmqCoyUJB3LQSHEGjzkN4oTaILbo\nXuLUJNSCnQM1QbPQskSbDG0qyc1sGOze2f5EZ+19mw/4zP4nqDyeYoygBGl8MqHKMVz3xvY7mrdv\nzfZqmp3SRFQkdMjSKUsRi4LXULOO5gB4pQihObooUZF8KUiHvCIUIj81CHwpA90ZY+VDJU4Yyefx\nbJ63PNigFiE2p1ZBmXJ514PWQLRB6I0wp9JsL1v3xIaNf21pOm7m+xcfMj/ovArFyKAEaRwylssx\nZDKZHTt2tLa2AtAdOJ9GP9p1b2y7o6W9ZbBf06yUbunE1qmVIrZGWYqYIiQSOhSkRlECo1hSRBxi\n49IhrwiFaI8vQdtnuYahaIkyEBvMzo+opRYhFqcWqMlpjtg50DpY9dAmEb0R5hSS7TYH3mhr+2Ln\na8fMPP47hxwW65IUigqiBGm8MZbLMYjhUI2NjbNmzbIsy7IsAJlMxrkMwKVPLsUSh7pu7bY7NnW0\nZPuoZhu6pRFbp/l4KE0sQlgtNVHQoShqFB4YRZEiZxeRU4pcOuQVlbiy5CIDXfZUMRCLi2iJ2Jwy\nQmyI4kP5Kg8Wz8tSDbQcrDpoDUSfBLOR5HrMgbfatl3Y+fpxM99/uZIlxWigBGlcMZbLMcjhUK2t\nrU5pcfnRhTJ55Uqg6/or3eYP1vdtyfUSamm6rRFbo8ygZopYlPAamgNQETWKFRjV+kVFNUW9SiVk\nSZCmuVINWUSWpVyHyhS+1DZowelArYIyCVlKc1oDOwtaA60Wdh2semI1wmokZrc5+GZr62c6/vWN\nwxaeMmNarItRKMpECdL4YYyXY4g4HEoGQ+Kpc/vmvtyXXmx+fvdOQi2q2Rq1NWprhOnUqqEmAJca\nEfBaagpLtxhXhGA1Ck/TDZ8UhSuQoWV911t2StjBXbtbXHeeOgudAzaIzYXZgQpfuEWoSWiKUxN2\nlmtpotXAroVVT6wGmJOoudt6b+natuf2+uBFBx7XVFsfcoUKRQVRgjROGOPlGMofDnXdmu3Xv7UV\nmk00i1JGqa0Ty9AsAu6rRqLyQn3BtjBiahQkRSV1KEh7fNE1Hxmz7JROLL0w3CnHjHThpFnoDGCg\nlvCFgwqngwlqEJbmdhZaGloadhp6HbHqYU6iufbO56/t3vi+vT5w6cFHRr82hSIxSpDGCWO5HEOZ\nw6FWbO+7+OXNzYMD0G1CbE2zCLhOrJQmdMhHjQAYZCheCc/UDYcaRZeiIB2iRozcHTNTCFCpwhbi\nH2KC2qA2pxaoAWqCpkBMaAZYGrbQpBpYNbDriFVP7D3mpn+1tn+uY+3xM993xbTUIAAAIABJREFU\nmbLhKYYZJUjjgbFcjqGc4VAtPbkv/7N5RVc3NAuaTQjTNAuATkuoUfR+I1k9IZYahafpgtTIKUVe\nHYqlQOH7Cn2SpzDttDh1lqXSsDlsE5oOYnOqQygTM0Bz0AywlHgQluZ2DbHrYNURa6f571dbt12b\n233R/u9Xc9cqhg8lSFXPWC7HUM5wqDv+vevil7ZAs6HZohyBUCODmpVSIxR0qEw18k3TRZSiIB3i\nhrv2XRSIabiOycyUV/wISwHIQdPBhCNcB9XBDVCDMJ2zFLdThKVgp2GniV0Hq55a2zqfv6Z701Ez\nTvivQw9JcG0KRUmUIFU3Y7kcQ+LhUM29ua+uaHluZzc0Bs0COKG2VCOxTTI1cuFyeAsSq5FvYOSU\nIo1a1JFI9EpRSRFiKbc3j5q6syir6wjENORZfJSJpQDYoCaozqkOZkDTOTUIyxJN59zgLEWYAZYm\ndhp2HbF2mhtXNPe+vbv1xmM/MLOmMfxqFYq4KEGqYsZyOYbEw6FWbO877W//hs4KagSXGonwSBKu\nRi58u47yxykeb+RcieISDEikRk4xiCRFBMwoPT7Jdxua012HDVGmHDM0MBOaDmZypkHTQXXODMIM\nwgzODOQfacJqYNemtnUM7ln8evtx04//3mEHl7xChSI6SpCqlTFejiGZAf361e3Xrd0OjYFKNbK8\nauQMj8IJSdYJXMk6gWvmCK+LIb++EmrkkiKuMa6xkHfkr0Bm0RfZGUgJcfIqk1OWOIjGmcWpBqqB\n65zq0MTsGDphOmd5fQIzCEuDpYndaW78V9v2m7SF31QGPEXlUIJUlYz9cgxx1ai5N/eVFS3Pd/aC\ncEdslP9hFfWzy+k6chE9WRc+3shFuBqFS5E3HYcA+YmymZQoeVinMnllybRTGqEa1zTOLVANTOOa\nDk3jXCdcI1znXC+ESkKW6jTr1a3Lz29/53+PO2tu/VgcbKCoOpQgVR9VUY4h1rXduX73xS80c52D\ncOi2U42kra7MC0ucrHM9DUnWRVej0lJEwPSy3rJTooQ4OZVJXgAVG5spQ8sBoHZaSJEGTQPXOMsv\ngGmEaZzpnGuUG9zWCUvBNijbYW24cs2Dx0w//vsqfacoGyVIVcb4KMfg5IbXdlz3xnZuMADQCyVK\nHWoUlKxLHB4JXOGRIChZl386/GoUFBLZ0UIlieZI4sljOpXJGTBJWdK1rMYpZQbljBZkiXKuQaOF\nUEljTCc6BdcJM8BThHWZG19u67xJ+/g3Dz4i1kUqFC7oaF/AeKa7u3vBggV33nlnpQ44zsoxtPTk\nFv5143VvbM8/F2qklfjlpT5TfZfAGR5pxbvHCo9c+Dq8JSXViKWsIjXyczHYhiUevhcQgtzRuS8z\nLPGQZxcXwA0ThFMjR40cIczQsgaxUtRMw0rDqoFVw61abtVyq45bdbAauNnIzUaem8qzMzE4G/3z\ntD0rtv79rBUPNvcmH02lUKgIaRi55ZZb2tra+vv7K3XA8VSOoaUnt/CxTZsHMwC4XujG9yTrvOFR\nipqIFh75IoKeioRHXmR4FEWNnDv6SlHQWTjlPDShR00dfGjCdOehROQkTicDJprTxQFl35IO6ICJ\nNOEa5ZyCU8gFrhFOOQoLXCNcBzco225Z31nzyP/3/k+p+dEVyVCCVGFs2966dWtzc/PDDz/8xBNP\nVPDI46kcwwtt/Qsf28h0BqlGeqBTIDEh5jqBb+iTODzyqpGT6GoUJEUsFXWorOuANGe4Du6SJZnE\nc1oehN+B2Gmaz9o5NIlzSvILGuEUXBf/Urbd2vz91Y8cPfWYS5tmy5ruYzOgV4xBlCBVmE2bNn3q\nU5+q+GHHUzmGp3albly3URTtKVKjCOFReO9RSfxlyW9lgvDIi28VhrhqFF2HgpBHkMpUUpZcoRK1\n05RzwoUmgaLwIJxyTgAKTgjXwCnh7eaWtV3Zn/XNv+zAA5zzhsAxo5WQqARTMirGN+pWqDCzZ8/+\n1a9+JZbXrl17yy23lH/M8VSO4etPtS7bVgMCpoeNtklGyXydIEq+LpxywiMn4WpUQor8bOJucq7x\nSUXKFCZLAFyhEgNhnHBOwAmEOHECTghE+o6iEDYRvoPZq3v7lu2a/J1DDhVndE5zZVmWmIlRyZXC\nhfqwK0xDQ8PChQvFckW+S+OpHMMZD256vqNPqlGU8EinCVN5rrnJvSTI1yULj4KSdQnVKIoO+W7s\nECdxcF9ZKupYKoRKzEzpNMt5moBTDqlJBCDghDiWAUpACajGn9v6d53Y3z74cHimuXLhlSvnGiVX\nEwf1cY5pxlM5hjMe3PTCtn5ulN7SCSUMwfm6iPhm6gSx8nVBhIdHkrLUiHAY5XWzybMXlMlXlpyh\nErU0py/cAEw7bRATTCcABQgHoBGAEBAhVASEc0I4IZxQPNPyVNtg//87+vjwSyspV3IGYadcRZ/w\nXlEtqA9sTPDKK6/8+te/bmpqmjNnTlNTU1NTE4Cjjz563JRjOOOhTS9s67cNDviFR0kpvwNJUGa+\nzotveBREaTUqGRWlA86S9dN/cbQAWdIsHTwfKjHdRnGvkjgcQZpwTjgATjgIQABCeN5RTwAuVoFo\neKfz5Y89t+PpjyTvWA2RlvAJ71EcSDmXlc9ibKIEaUzQ1NR03nnnAVi1alVra2tbgTfeeGMMfnPi\nDoc66/5NL+wI9b578nUJKKcDyZcgf135OMOj5GoUJEJB27jEKUCWbN0vVCpO3xlaFjYIAK4RALww\ntIsAYmZ1AsILk01RNNubP/7cY0+VoUlBhEx4j1JypTKBYxDV7mOCpqamc889F4D4V3DooYeOQTVC\nzOFQZ92/6bkd/QCc4VFJfAt7R6ekLPkSrkAlO5CCatYJfEvVebaJpkZRpChkL6cyeWRpKFRy9CoF\naRLLZ+3y/wKF8AgAAecihQdGyVZ7y0l/f/qfn/hYkitPSnS5iuWzGMm3MNFQjauIR6zhUGfdv+n5\n7f3eeiA8jsUuem1vF3EdDclIB1+Yb74uJDxy41WjZFLke5CcY/CsQ5acmoRCqJTXJFPnhkks3alJ\nnDAxP3o+NkKRJoGTfM8SBcc7H/w7XfWJj1bgLVSCKD4LqVLOjquenp4RvMyJhRIkRQxiDYc66/5N\nL7bmXQwiPHITvwOpllbiF7kShHQgeUkSHiVSI1Lj3oZnAmwk4viuaKmgSXD2KklNMoYc4bRQKdy0\n05qY6Q9AoRdJhEcM4IRwDkbACWxKGX/ne280/uyY95d8I6OOECrfNMCWLVtG+momDEqQFFGJNRxK\nqFGk48bvNIplsQvB12LnKnZX0mI3QgSrkVeEgl71ESdxWClLAaGSb5eS1CQAWZZKwUahO4kXpIiD\ncAIGwkBs0m9rWNv56uLVxi+POzre21dMDFRxVUUkxHCouGpkFf8AhnQgkbInmAgixFzni4GwdGKs\nwCicorKnJcOjAMLVyLuxeLhfSJtFglc4u7wqeamyMCsAbpiit8zQsmmao+AG7BRYittpbqe5VQtb\n1GOt52YjctOQ24sMzNF61u16acn6f0e/bMXEQUVIitLEGg5148qOqLHRxCPihHtu/MKjWFLku687\nYEqbRaFScfouyObgipPSsDlE7o5w2PmUHYiImSwQm/SblP615Xmb0ysOOzDxW1CMS5QgDSMf+9jH\n1q9fP9pXUS6xyjGs3Np/48oOuDJDpOhZiKMhlue7Uom7MYcrPKq0GrkOUiRLLk0yNWF8EOk7X00C\nMORxALIsVQPOAcYJBzjASH6ZgdiE2CAWpSboIy0rGNOuPHxe+W9EMW5QgqQII1Y5hpVb+8+5e7N3\nva3HnsFIMWK4ZcnZqyRqQzi6lHw0qbjkHYAsSxngHLboQ2KwGSGcE5sQm+c1yaQ0p9EXO9es2DHt\nQ3s3jsLbVoxJVB+SIpBY5Ri2dufOuWcz4I6Hho9MtDJE/TH/6spAS3Q5eYYGIZGoMhy3nnd4eKTV\nmCGPqMcM7lLyWtXz5YWMHAACnqa5FCwDTPYn1XC7BnYtt+tgNXCrEeZ0ZGbRfmq1/t93Vqg5/RQS\nJUgKf+KWY/jm8jZhTxuxeIjxcqVvMH6GIFNqF2bmLdAo+/J8KalG4buHyFIUTRIITXIZHABQI6dr\nOQBpmkvDSsEWNoc07DTsWm7VcrsWZgPMRuSmkcFZtI9arRe++nT4NSsmDkqQFP7EKsdwzj2bV7YU\njAwFPRqeH+QJRvSRsISXVCNJkCyRGrMosPNokst359UkOCpWpGGlwVLcNjhLczsFVsOtWm7XcWsS\nzMkwp5HB2bQPdtu3//VG1LepGNcoQVL4EKscw40rO4bUyEHklFVsBgoDMcPpd5srkpAhZWXwYhNr\ndgkHWjr2jr6yRNJWUajkUcRwTXIawQEQ8BRsqUlpsBpu18ISmjQFuWl0YJbW+3LHv/7n7S1xr18x\n/lCCpHATqxzDyq39S57vdI4l1azkkZFtl0iIDbIKaAyAQeJzouwIa0/liB4bRdzXX5M845PyTx2a\nhOLOpDQsCq6BCVlKczsN0aUkNSk7gwzsrfc+17F2RWd34nehGB8oQVIUEascw9bu3P+5YwsAYo+t\n9FwfjxRCuWBx/BglO5MSkhsF46tWYxLqjmejaFKIwUF2JtXASsM2ODPAUrDTsNPcruV2Lax6WJNg\nTiGZmXSQWduvevOlrQODw/EGFdWCEiTFELHKMQD41mPbomymmUU/9MTyues4K/FDnImWpispRQOR\nhcQVRWX8gqqKQHOVCfvKuga/VGF0TQrvTAKQhpXizOAsxVkKLJ333VkNMBthTiUDe9F+brd947W3\nKv/eFNWDEiRFnriz0974QsfK5pGoyJCJlqYrKUUh/u/BMpJ12WhKKRAThFeWkvm6Qw7qu/46/6ls\nVzw/pEOJ835eTYLD/i6CJOT7k5jBbWFwSHO7BlYttxtgNcKcRgdm0IGWgXevf6s52WUoxgFqYKwC\nALLZbKzZaVe29C95odO5hliEJzN823pQfdWcbUSZEqmfp8Inje2HHlTRbhC6t8RqBpprFvPiV3Xf\ninamnQ6awpzm9CgFv4eJ5ma27PZ8M1734/wVZs3GFc9bp318IGsOjUvVaky7uJ4QqTGLxsyKAbOO\nwkK+4Z1zhlnTTqdpThQIzHDdIMQGS8G2ObFALWKZoJNAs6DT6MCArj/V8daHZkz70N6TKtwKimpA\nRUgKZDKZzs7O6PORA1iyorP0Rh6oX7LOF4tH+lOppN0uxGjn62vwJa7RTg5FImbZubgIXUp20AQT\nDoQmSVkCINTIu6WP7y40cSfMlBETdzWw0tx2Ju7E+KS6QmfSdDpomu3fefOVku9IMS5RgjTREeUY\nZs6cGV2N8sm6yOoSfO7AH3oePIhJGO1KSlEyX4NLpcJFq0xfgxYlfVe5wVzNzay5eaiKoK8aCWLl\n7obKxZKhpyGJOwAGmHikYKc4SyNvcGhEbioZnK71M2vb116p+iKQigQoQZrQyHIM6XQ64i4rm/uX\nPLsTNgUQpEnSxaDH75UIcX4Pq68hSjdSkK+h/G6kosTXaBjtXLg0KZK7Qfd33MERJNXASsHWOTM4\nM7htiCCJ22lu18KuhzWJ5KaRwan6wDt9m1e091b8fSnGOEqQJjSxyjEIlqzoKnru0CQSPAKpTKNd\nSV9DP08hVIqyoAjwNfgWEAqvaBcUG5l2Xte9WTs6ejJzy621FThKhHHOJd0NADQw8RBGcFHyLg27\nhlu1sBpgNtLsNDqYNdsXr1lTgctWVBVKkCYuscoxCG58oXPl5sBUT3KCo6KcXUKKInYjWX63um9G\nzndlUDdSdgSydtmiFgickjyYC79gvLeh4cIv+O+YNnrSRo/PtbmCJGcliFJDZYeu1qeekK2D6Zzp\n4PnBSTyfuKuDXQdrEnKTSWaaNsDRfP1bLZHfpWI8oARpghKrHINga7e55Nkunxf8wh1XSCSJ7msw\nS0VF5XQjVWo0kowaRjdrV9LXMHcuveaH6Wt+GDUxm5goQRIF18GELBmcGaK8EFiK23Ww6mBNItkp\nNDOJZh9rX//idh+xVIxXlCBNRGKVY5B865HSw2C9WbsS3Uh+voaS3UghBYT6SuXuoo9GGrNZuwRB\nEoC5c+nV16Tf29Awd27Cb32sekIS354kA0wTQRKYDjlg1hZxUj3MSSQ7WRscNHf8z/pNya5WUY0o\nQZpwxC3HIFjZ3L9y8wDnhFmUeaOcCHHPyHQjhW3g72XQ4elGipW1GwmvXba0AkUxfwOYO5f+4+k6\nZ6gUEjmVUyLPNVm7K0gCoIHrYFohSDLA0rBTnNXCqoM9iWQbSXaqNvj2rq0rOlSQNFFQgjSxiFuO\nQfKtR7Yzm/JCzTofTSoTv6hIZO1CupFE1i4kdydCJd/RSNG9dhFrCI1k1i5ZkCQQoZIUoauvSV99\nTXmpPL8xSYKgMUnCbieCJD0fJNkGF7m7fOKukWQn08z0uj3fevXNsi5PUT0oQZpAZDKZrq6uWP1G\nghtXdDZ3uSsXBGmSzNoFmb/d3Ugxs3Y5ZiA0azfADcTP2kXx2mXcOb3KZO1GOEgSJEvflZiZwlNr\n1Yk3SCJ50x3XC0FSCqwwMsmuJ+Ykkm3UB/ux7fo3W2Ndp6JKUYI0UchkMtHnI3fh72VwadJIZe1C\nanKLrJ3vBiJICsnaha8MD5KkMg1vkBTBbhdLk0T6Lvr20fG1NghchRt0Lhx3Qx6HFOw0ROLObCC5\nSSQ7RR98sH3DcFynYqyhBGlCIMoxJFOjbz26jdnB94nf6BRvkBSJRObviFk7/32TjpCVDHuQNGyJ\nO0HJCMk7LYU/pSr1yWkp8ocFL1RchQYuHnmDgyjAClYLu46YjSQ7Scv05jq/8vLGSFeiqGaUII1/\nZDmGBGq0tTv359fDBswXaVVwkJQ4a+dr/i7HayeCpEpl7RybJbE2lA6SXFQ6cVcS32kphgieTzac\nFB3aLMVtnXOdMw3c4EznLAWWgl0Duw5WPTEbSHayPvh6X0tzX1gJXcU4QAnS+CdBOQbJNx9uFzGQ\nMNf5WuxiGRyGz2sngqQQr12ItcGToNP8Vvpn7VxSJLN25QRJZbobKqhJPoeKODW9X2k7gdPaIOIk\nmg+SRITEdeSrr9bAqiPmJJJroNm+XNf1aprz8Y4SpHFOgnIMkpUt/S9uHIRHcsIUqPBSyayd/wjZ\nmF67koVWS4ZK7gOOjSApDE+QNNya5KKoZEMwgaXtCnomrQ0ACKBzpnOuSdMdZGeSWUfMeppt0DKv\nd7et2N5XmbehGJMoQRrPJCjH4OTGZ7sQID8xJKpApUbIiiAppNBqSF27YbI2jHSQFFmThkOWyuy4\non73QZpbIkLSwHXODc514bgDS8OuI2YDMetprtfsumNrpEmKFVWKEqRxS7JyDJKWPeaLGzMhShP4\nUqkgqUTWztYRkLVLZm0oGSS5rA2D+SReWNWGkkFSFLud1CRnkFRZTUJ5oVLCfSPk9LxZOxQn7oTp\nLl1w3NUSs4Fm62ju1T3bX2gbiXmKFaOCEqTxSbJyDE6WPN/JWQmPnFOTEgRJI2Nt8GW4gyRJSJDk\nJFLirpQmBTF86bsiZA2h4tFIvrP2BUCEtUHjXAyYzY9Mgl0HsxZmHc3tzu1a1tI+TO9AMeooQRqH\nJC7H4OTPq3ukINkWzT88/m9/HfIESS4SWxtEkJTM2uBbtWEEgqSKJe5KwTNGSJwUV5aGT8a8I2QB\n1HArzW0NnIogqTA4SSTuaolVR8wGmq3Xsqt6t73YqoKk8YkSpPFG4nIMTr716DZmaihI0dALHHaC\nJF70qg2CYGuDLyHWhtimhmEIksKpbOIOoX08QpZKKs0wdT6FILN2KCTuKLgu+pNEjTvYtcSqJWY9\nze3M7r69uWMkL08xYihBGleUU47ByZ9fDytnaUdxNMQNkkQ0Vpy1E0FSXGtD3CApv1cZQVLJwg0J\nEnfDpEkCrzLJNeFSFMvRYIdWEnJ2IxW/Qii4xof6k/LDkohVQ6w6mq3XMv/s3R79MhRVhBKk8UM5\n5Ric3LOmW4ZHQdsEvRQlSHKRD5Jc/VWV83/HHSQba0xSUHW7KNGSb+KO+nWkuamEJgkixkyRKHty\nDQpew+0abunC3cAL0yaBC02qhVlLrFpqdub23LFxZwWuWTHGUII0TiinHIOLP6/uQagaCewAR8PQ\nckCQJLN2/hI1GkFSSE/S0FOiwxEklW8BdzKkSXyoTcI6k7ya5Oe7K7+2kOuA7lWRjRVBOLuRUnyo\ngG+hMykvS/lhScROE7OO5upo9ob1qpLQOEQJ0jihnHIMTlY2D7y4IeNSI9Oipp8+BWlSEOGl7Xz9\n3y6GKUhyH8ovSJIMh7vBSaTEHfwikjJCpZKUf5xwo50ncYdCT5IMkuw0sWuJVUPNGprryu1U1obx\nhxKk8UA55Rhc3L26mzmGkTilKEiWvEQPkgSJJ6QYlSAp6Kk3SApJ3IU77iquSWXKCc/6vZEywiNp\ntAsiDVvjTOP5OMkYqt1gp2HVabmUZl63rjXxBSjGJkqQqp4yyzG4+PNrvbJeqq/8uFaWGSQlGCRb\nwSDJt3BDSHU7SUR3g8QbJDkZAU1CGbLEM4YzkTgcOMOjGm7Xcgv5CAkUvFBSiBuwU8SqIVaaWGlq\ntfDdw3pVipFHCVJ1U2Y5Bhf3rOmWY49CgqEomhS3J2nEgiQOgmjpuyizm0d3N8TuTCpTk0JlKaIy\nhW2ZDbwGeXklJx50zo3kS8ECns/aiXGyaWKniVVLcx3Z3df/Sw2SHVcoQapiyi/H4OKnT+9ijCBU\njQQhmlSScLvdcARJzup2/dyVqQsMkoaOWWwBj+tuiNuZ5NQkYg/pXAlNihwqCaQyOVXHd+VIUsOL\nepIIQAvuhnxPErEN2DXETBMrrZl3blGCNK5QglStVKQcg5OVLf3Nu03OScSOoqDNEgdJbsqz22WZ\njnjpu6gW8KHzlkrcRTF/h2sSYa4MZ2gRhzihkot4CpSo98i3LEU4HNAAWgiVdDAddorYKWKnqFVD\nzdZMz4p2Vf97/KAEqSqpSDkGF3ev7mFm0f1gWdT1cO3i1CTfIMmrSYLwwg2J7XbO6nZ28b3tWwK8\n4u6GBI47lNIkV5k7lyaVTt8hhixFwnsov3xdAlxGuxrYtdzSOKOFwg2iJykF2wBLEytNrFTN4PPt\nYRNIKqoLJUjVR6XKMbh4cUPGGR555cd3pa8mhbgb4pUAjxAkhVS38y0BHt0C7nI3DF/izklcTYJv\n+s7008WKyFKoGiWDEuZdWVvI3Yk5zimXmsQMwgxip4htUDulmWpCivGEEqQqo1LlGFw8ttls3jMU\nrfiqkXzJ9apv7i4kceciQXW7uCXAk1nAJZVK3EU0OAQNTgryOMCrSZwE6oSQpbgqIvaKo2fyCksW\nMjdIiRn/CLijJ4lrYAZsg9gpYqWp2ZrrVRNSjBvK/etm3PPUU0899thjGzZsaGhoOOSQQ7785S/P\nnTs3fJcnnnji2Wefda2cPn36FVdcUebFVLAcg4vX2sBMTUhLiBo5roTqus8ftrZFNb/1TohFuM4B\naCaxjaFBT66nxKJcZ7A06EMD+DnTCXX/fuVsI6WZGWbU0CFxG2RGLTUHWKqOuoe89PFUA8mJf+XK\nfhj1MPuh18MCIBYGoNc58kiD0GudT4ku/pDPQKuBPfSUaDXcBpAhek1+A13ko+RClqVEkQLTTguz\nGTNTcnQOMQ0xgJTmdOf839TU5eQOmqk7i8XRnMFSxdouVCfl93PPSZG6pINnhQgXoVyoLvrhldtI\nc0c53A0iSNK5nSJWitgpmnu+s/fUpvqSB1GMfZQghXHDDTcsW7YsnU4ffvjhPT0999133yOPPPKb\n3/zmlFNOCdnr8ccff+aZZ+rri74h++yzT/nXU6lyDF5e28Z56WnVinBqkmlRo7AsNYlZlBYvwKLQ\nGRyaJNBNWI6fKWpR5lQ1oUm2Dm3ot9W2dU2zTGYYDhESmpRhqRqPCAll6uepeocIBShTXpPyO0Kv\ngyXFRmjSINFquQ2HJglCNCmEkpqEQpkDlNIkADFkSZIslRdfjSJSw21vbxxxBEmimJBOmE4tnVrL\nWrZffWzFvKaKUUQJUiCrVq1atmzZfvvt96c//ampqQnAk08+efnll1911VVPPvlkSIyyfv364447\n7q677qrs9VSwHIOXrbuoECRneMSK+3Wo4ZasIE0aOoJXk4qJFCRJbB2a5QqShCaJIMl1cN8gSWiS\nS4TEUxEkDW1ZrEwe7clrkkAESd7NHBsEBkkopUkoDpVcmoTiuto+oRKiyVJ0SmX8So5A8h2DFQL3\nCZJsndgpYqep1ZLpfbG1/5Q5KkiqelQfUiB/+ctfAHzve98TagRg4cKFn/zkJ9vb21euXBm018DA\nQGtr6/z58yt7MZUtx+Di7jXd4MQqnnyPeazYzCTelU4Bk51JYcOSRmOcbHR3g2tY0kCo485VCLyc\nziRE6E8KGjMLvy4l/3jFd7hSXPyOEBQeyev09Xz72jp8qYWw2wn/N9NIPkgyiG0Q29DM53co8/d4\nQAlSIKtWrdI07cMf/rBz5cc+9jEAL7/8ctBe7733Huf8sMMOq+CVVLYcg5e7/9XDREZNOhGCq6CG\nvOTE67iroAU8yjjZ6O6GkGFJMljzddxJXJrkdYFH1CRJ+ZqEkrKUQJkC9nKdpWR4FItabolIVAMn\nhWJCGpgGrhOmE6YRWyP2io7uCp5UMVqolJ0/pml2dHTMmTPHlZo78MADAbS2tgbtuH79egCzZs1a\nsmTJunXramtrDz300AsvvDBxcCPKMVRwAKyXFzYOMj50J5SUHGYSZ/ouMHHHAQIEJ+7KdzeIniS5\nMrq7oeKJu+idSVEMDoicuwMQkr5DUMdSvr2Kv/6+Cb1SuhWiRl5/na+BMDpa3m4HDVwjXOe2yN0Z\n1NpsqQhpPKAiJH96enoYY5MnT3atF2u6uwP/HBOCdNllly1btmzXrl1kgVJoAAAgAElEQVQvv/zy\n7373u9NPP/2f//xngsuoeDkGL/es6WYmZYyI8ChiAOTazD9xZydP3An8E3cRLOCVTdy5hsq6atwF\nucCTxUm+A2YRHCfBL1SKES05kZGT8xFKRCNDghoNvhBw8dCKgySd2FszfX97d8eePXv6+voymUwm\nk7GsCnWYKUYQFSH5Y5omAF13t49YI171Zf369YSQiy+++KKLLkqn07Zt33zzzb/4xS+uuOKKJ554\nIsgg99BDD/36178G0NTUJLqsPvjBDwLo6uoavkyd4MXmARZay5mY4H4/O1HipCiOO8EwuRuE404E\nSXJlmY47+VS6wEWc5HKBo4w4CYXCo1HiJBTbHAoNqHtnEA+LlmLiK0Xh4ZGE+3X+RYHkrQ2FkUmE\na4SJB6X2i7sHPzBDsxwA0Au4lis+cEJREZQg+aNpGvyER6wRr/ryi1/8wjRNafLWNO2SSy7ZtGnT\no48++uSTT55//vm+e5177rknnHBCW1tbW1sbgLa2tlWrVgGo+ABYLy9uHpAVvp1xD3G8dbnsUqYQ\nTZJ4NUmm8vIHL07cSQu4eOpvAS8wAok7oUlSmXxd4ENnCdYkL+X77lBsB4cjfYeADB4cWpJMmYKi\noqCuo6G6fIUgjzMKP0dDpHnf8z1JIkjiOmyNMJ3aGrHva++5/sR5zo2FJklxEmGTxKlSKEiUU7QU\nI49qd39qa2sBDA4OutaLNeJVX2bOnOld+YlPfOLRRx/997//HXJGGRtJHnrooRH4O27LDs64u/gC\nCfilEuudsuTSJImvCzy/i51kWJJvkDT0LPKwJKFJLmVyhUexOpOCRiZJTcpfTOjIpFiaBMe8q95Q\nCdFkCcXSUlKcwrNzbptfZF+DyFVGkSIBKTwoGAUTBgcdTCN2S9bdjSQlRzx15SeccmVZViaTcT4N\nkauIl6pIgGpcfxoaGiZNmrR161bbtp3x0JYtWwDMnj07aEfGGCGEkKIkmPgm7Ny5c7guNyn3vNHN\nGBEjkGR4FKRGElcSz6lJI5G4k0HSsCXu5NPw8g0lDQ7RKzgggibBMdFqSKiEgAweAmSpcJCEdgOf\n/qqi8q/uV6OMQAqqUZs/puhJItA4p4TrhFHCNMqoZq1s6T95v6ijkVxy5UIGVWJZdkplMpmenp6I\np1DERZkaApk/f34ul3vrrbecK1evXg3g8MMP991ly5Yt8+fP/853vuNaL5wOBx100PBcaRlw8OIw\npqQa+W7mzPX5Vh4qWXe1nDnOoxddLTlbEkJd4JKSBoeS5cB9p6gIH5+E0JJ3XqeDN0wRfgfNqsyf\nof7uiQA1KtNfN3ScoQipUOCOiIfNwSs4GkkEQw0NDQ0NDVOmTJkxY8asWbNmzZo1b968xsbGSp1F\n4UIJUiCf+MQnACxZskSu2b59+5///GdN0xYsWCDW7Nq168knn3zyySdF39LcuXOnTZv2zDPPvPvu\nu3KvPXv23Hrrrbqun3baaSP7DkrzwsYM4yRK8TovUTQpZGqloNmSXMOS5NMoM1NUxHHneirrrrom\nTAofLTt09jimO8TXJJf7LoosgReUKemAoaB9gzJ1XjWKPiS26DhDC1xqEgGjYDphlPIXtinzd3Wj\nUnaBLFq06J577nnllVfOO++8M888c8eOHY8//vjAwMBXvvIV6Vn497///e1vfxvAK6+8MnnyZELI\ntdde++1vf3vRokWf+cxnDj300O3bt991112dnZ2LFy8+4IADRvUN+fDCxkExRayQE5fGNGSLMjx9\nae+wktK5u+pK3IXUXXU8LWFwSGa6Q2FCIGfuDoDvECUEpO/g6VWCX8eSo22LFT04p+fd2HuKojV+\nVm9Xvs47KDhDIv0iEQDghIBwTgFKOCWcwt6cdXf6KqoLJUiBGIZx++23//jHP37qqadE4q6+vv67\n3/3ul7/85ZC9Fi5cuHTp0p/97Ge33XYbAELIfvvtd9NNN43B8AhA806Lc/98vUuN5BqXLIVokiRE\nk4Icd4JIRVcLJHbcRXGB+3YmDR0wvunO1+Pg9d0hwA6OUk4HFPcqIVSWJMlipnA1KhkeCTEO8jUM\n+kkUKWgSLQxLomCE8C0Dg819ubkN8QrlKcYOJHaRZ8VIceihh4rOp2HinrV7vnLnzpypMZswkzjD\nI68auXDJklOTnIIkDQ7ScScnp6CeBRkkCU2SUZHQJPlUaFLecScFSbMAyCBJaJJw3MkgSWiSdNyJ\nOEloknTcCU2S4VFD8VMRJ0kpEgtyZJI0OAhNkqVX5XppuhtaU9jG6buTUuScQTXt8ARKTYIjTso3\nneHuA2QBBVXDlSkivjm6IDWS4ZEUJKfFLv+vmPaQaHL+w3wvHdEGoQ8SfRBaP4x+6L1IdfN0N0vv\n4nXddu0uq747Vz+QnfTUqccM91QUW7ZsmTdv3rCeYsKi+pAmLi27LM4JZ+6yCyXVyLuNU8zCO5MS\n1LjzfeqeUjZmjTtJ9M4kSXgFB4nX4BClPwl+dRwQ2qUU0qsEv46l/HpTlw/vq+GE7BhXjeJT8IIi\nP5ybyv4kwghhL7SpbqQqRgnSxKVlt8kYOI+tRr5bVkaTxKHKKAQexXEXpe5qn5/1zmVwGNo3suku\nXJN8awsF2RyCnA4IlqWgEj5OcfIqTfirrlM4ryHo8px483Vez3eWBP1M8YLpjlNwQgDCXlSzx1Yz\nSpAmLi/8O8s5ERWtI7q9XUTUJEnp+SlCHXcSt+MuQo27KC7wkFrgkpKmu3I0CRGsdxFDJfjJEoID\nJvdmceIn7zGDTN7h4ZHL0ZAr/DoxRzcjd/wLAKRQu4EwQjghsGptKKoWJUgTly19Oc7h7EJ0CczU\nwUHXw3uQKJoUYiuvTOJO4k7cEZRygScYmSQ1qfA0qhG8IpqEUqFSdFmKKE4h+B7BdTpvsk4SXqDB\nqUNO5B1L8v8OxUmEsC0DmUiXrhiTKEGaoLT05HhG45xw2/9r7ys/vrIUokmSchJ3QoSCEnfuziQA\nQ4m7/MoEnUlRNMnVmYTgaZMiapLvEKWgLqWgUAmRZUkQV5nCxSxIjZy4wiOnnUHgddYNeqRLyJLU\npLzjjvDmwczWbvcE9opqQQnSBOWfLYVUe+EPTqeu+KqR81XXBkGaFLczqbBpBSanSNCZFMvgMBya\nhAg2h4ihEoJlqaQyhT9894Wf5jnP7k3W+YZHrg4kX883zz+cwwVAAEI4AQdhLT0VKGeuGBWUIE1U\nbMI4EY4GV0wTrkZBmwW5IWJ1JkWcVbac8g0hnUmSKAYHyQhrEkqFSj6yZLkVKFyZ4uIrRb5q5CUo\na+eaDx6OJiIFNeKeOIkQDsKbe1WEVK0oQZqgvLhlkDOUOQgtRJMq1ZnkStxJErvAJeGdSclMdy73\nXVxNKmm9ixIqwStLnHijJYFUprjiFLKj60TOKwkKj+TwIwS751HcS8dBuFQigJChiKm5TwlStaIE\naaLCwUHgEaSI4VHQ9tE1aVRc4K7OpJAyd5Ig0124JnkHJ0XRJKCEzQGRQyUEJ/GCTNhOjSn58D2C\n9+DhahR91gkvvPAfd6wg4CB8S78SpGpFCdIEpXm35QyPghJuDYO9rod3m4iaJH85EncmJXOBh3Qm\nScI7k3wNDqi0JkW0OcCjSXFlCX7KUSa+B/RVIy9B9eu87YZC23IQXvjXgxKkKkYJ0gSluccE3AGS\nU1qC5Md3fRRNYlZYZ5KLke9MSqZJ+VfL06SsVJ1SXUq+6Tt4Qo1YslSOOAXt7jqd82KCwiNnvs7b\ngeSkHzoF45wwDuSVKY/oQ9rSnw3ZXTGWUYI0QSFZKjzfCYfEJtOk4M6kuC5wSazOpFiaJIhlBEci\nTWIgEbuUEDlUQrAs+VoMnOLkKzMlN3CewnUZ3gseWggt7x04vxQM0YfEAIA4q40QqF+1KkZ9dBOU\nLT0m5wEjD/30xncb12ZRfHdJOpNKucCjj0yS+BocJL4GB0kCTfL67gbza7RkXUohoVJJWUJAwOQl\nVgjlK0UR1cgbHvkaviXS+e2XteMtu1XKrlpRgjRRMfPGWYmUkyhqJImoSYkNDnnidCY59grrTAo3\nOEQx3SGyJsHhBS/H5hASKiWWpYjiFETQEXxHv8KTqYNf8Tong8UNKGEgzKNJBBwgnKBljxqKVJUo\nQZqIbN1jchBeyL5HL6jqS7gmSeJWcBCU25k09jQJZdgcEBwqwdsrEyBLvv4CpziFqFSUbbyBUcAY\nWB+LfLESu4XK6aqXKbt8qMTFH1gE4DzFWpTzuzpRgjQRaenO5V2yHrzhUSrb63z4HjBEk8I7k0Ko\nzMgkiZ/BQRKlqhBGXJMQ3KUUMVTyXYNQZXISRX68h3Wd3fdpxi8VKfEMRfLRMw4wER5x4o6TcpT4\nOWgUYx8lSBORPXv2cEQaFetVoCBZSqxJI1fmrgyDw8hrUpmhkq8shShT0tmJig7iPaPvU68aBZnr\nvCXs+qGLthWBESu47ETNevFnFjeYc95hRRWhBGnCsWfPHsu2XHX8qZ86BcVDCJAll82hspqUJ2aZ\nu+HQJIFXk1xPE2tSOaFSuCwFrRQ4xSlEpUpu4yuHzouUF+99m/KpK1/n6kDq4ykCMKlJ+SBpqMBd\nS49K2VUlSpAmFn19fX19fYP6ZNf6yZkMAMKZXBOiRs5tfGVJLgdpkiSKJuW3TNSZNEQi012UwUkI\nLnYXV5NKhkryCl2hUklZihgweYkiUa7Delc6ry14X39rg3yn/R5N4pwwma/jZMg0SvxHzCrGPkqQ\nJhCZTKarq2vGjBnwDIkV1GeSzLaZQJMiGhwEiTuTyjc4SEZGk+CXvnNu4ErfRZclhAZM0fXJl6Dd\n3Vk7v3yjK1nnDI9883XOpxYIA2G8YLfjQbUbFFWDEqSJQiaTaW9vnzVrVk1NTcmNo4RH4dtH16QQ\ng0OZnUmOzSpmusOwaVKC9B38ZMn5xoNkKbjGtu59xNrAuVnQZYQn67xrBgKuVuTrmGOErO9miipC\nCdKEwLKsIjUiIJX+8kbXJElZBodko2UlAaa7SmlSLI9DSJcSIjgdEJrBg58sIU7WLlYI5aNhLkUM\nHXSFgPDI0W5DketQB5LsQxKTx3KAK3GqSpQgjX8sy2ptbZ0xY0aU2MgfzrRcr/NBmO3dKjyuGhp4\nOzwGh5KjZaMYHCRlahI84uTUJPHzKhcQuUsJnlApPIPnTeKFKFPirF1gOOUJjELUSKzJEvcvkjc8\nkk0qFIiBClkS/UlB4xkUVYESpPFPV1fXlClTGhoaIm7v0hUt16uZ7r4lag1oOR/5cdkcgoo+JNak\n/JYeg4OvJpU2OEQw3aFsTZIL/TDK6VLyDZUQmsFDcMAU5C/wTcoJpQl5yX0Qb6owoEvMtYaByHft\neiNeU8MgNwopu/yYJADwm1RFUS0oQRrntLe367o+ZcoU50pOomY0fFXH+WqQLMnlSpnu8hsEGByG\nLilWBYcCJTVJ4NUkQXRNQqIupaBQKSSDV1KWUEqZ3BtHC558pSiiGjmTdc7eI998HQAOYnPCOGGc\n5hN3YgZkKG9DtaIEaTzT3t4OQNjq/HF8bzW/LFwUfGUpliZFqThefjnwKKY7RBuchOACrOVrUvRQ\nCcEZPATLUogyRdenoN29p/O9Hr8ZyktMOSEWXKO+RL5OWhvEPOaqA6l6UYI0bhFqNGvWLO9L+05J\nSV+D/FuyNjvg2iw8PCq5cTJNqojBoRzTnSS6JsXN3SFAkyK671Aqg+eVpYjKVHhVj/XwO4L74CE5\nRjiNG8HhkcSpSXnbNy8MRRKmBmDfyfFmZFeMEZQgjU/27NljWZavGglIodhKBTPuyTQpZedH4w5r\nBYdA010EIzhKaRLi9yfBz+YAv/QdApwOCM7gwSNL3g3ymxWUKbzkdkR8DxXU4+W9Hu/AI0HBAzKk\nMX08JVqYc7DCaKRCvWBCuOpGqlaUII1DRDmGOXPmBG3A8+FRYcZNBzRp4k6QQJNS1lB5iOEzOCCy\nEXyENWnAI0VB6bugXqUEsuRbrjSZOIXsFXIxRdnIYp93SHjkytcNsJQFzQZlPG+0y6fsGPabnLw0\nn2IUUYI03nCWYwhiXn0KhPsORdLNfOIuKF9Hcr3iEXTwBJo0MgYHSbgRHCOlSQXvMomYvkNArxKi\nyVJ0ZUKxzIQ/fHf3SRIGjJoaDM7sOZFN5O5DEqYGEM4pA4WKjaocJUjjiojlGDjAU/mx7dE7gF06\nFCJL5WhSeFUhQZjBIaCqUJARfPg0SdYFl3Ofu0YmIbRLCTFDJZSSJfgFTHAoU/gkrSUJOk7IGF7f\nriMUh0deO0NxH5KIjagcHgtOCFOmhmpFCdL4wV2OIQQOUjA1BM9jXkSQ9gTJUrghwnd8UpSqQomn\nqEB8Tco/S6pJqGiXEkJDpQSy5FvM1Ckq4UIVZZsQ11+RKaM4L+dUowFo8AuP+nle7zkH49QulGwQ\nFjvCsO8UZWqoSpQgjRNilWPYb6rBDUbyf1CWPnhIgk5u4N3GZQd3jbeVlcVHzHQ3dGGJBsxiRDSp\nZPrON1SCnyw5qx6ExC5BZbadRJGf8MMGpel8u44k8q8l3/BogKVyXMu77PL+OgJOiKl+1qoV9cmN\nE+KWY5g3KUWKK9rl9HL7gUuGSk5NclYWr6Am5YlsBA8cMDvimpQ4VArqWGIgrmjJu33+4gsSElGf\nvITs7r7CYDWS9Bel7NyxTn9RN5IoHUTF2FhwCo55yvNdtShBGg/4lmMIZ+40nRRNaYacUWKC6igE\nhUpyOYoRfOholTPdScagJoWk7xAhVEJoxxI8STygRL+RS5+iPLwH8R7flabzqpFv15FsFldsJBYY\nCmUaCik7wsm8WmWxq1aUIFU9pcsx+LGfSLKT0oUoS+brouwSS5OSTeUXYroLN4KjEpoUXu8OAZoU\nPX0XUZaGlqH79i15CyJUytTgOppzZZZQ38AIAWokcYVHLk0aZIbJdDs/NpbmazQwMrdRRUjVihKk\n6iakHEM4+002KOFk2FyyldKkKFWFJOVoUmGzhJqEUjVYUaxJcdN3CDDgITSDh+CAybdUj1Oc7Mg/\nDiGSJk7ECn/1OC8m49h+wK1DQ8k690QexbZvDkctO1BhsSO2ctlVK0qQqpiS5RhC2HeqTuhQ9a/d\ntbUVvzyvJhE2JC9BmkR5XiIrY3BA6cFJJYo4IIkmybnPveWFEqTvEoRKJWUJocoEIEdokOkuisHB\nJ09YfIXyLyHX0OCgWkGyQWQDAuCArKwqwiPCySn71vtekmLsowSpWilZjiGc/SanKIEw2mX0oR+O\nXHpS6Z2zve5HAC5NolbG+dRXkyZnhraprCYFGcExDJqEpF1KiB8qhctSeB6vsKXmfCApQQfxJA99\nwjuXGjnDI0m/Q5MGmZFhKc4J45SBOj3fXEVIVYsSpKokSjmGcPadbBDKKc1HSOV+g4NlyaVJrvFJ\nvppUeSO4vJixrUm+5uaSoRI8GbwBd47OP2AKqh3n0ievxpTcIOgsg8VpOm/ZJPlOURwY+YZHAHJc\nE74GgIJTMEJtcvI8FSFVK0qQqo+I5RjC2a8xRTROCCcU/kWEggiZFjZAlqJrku/gpKFXh2dwUonC\nQhgJTSqZvgsKlXwzeCj83IcHTHBoRpA4FW0cOYTy67jSfQMjeIr4eX3e+Vf9U3bEmbIjnBCO/aZW\nwJ2hGBXUJzcsPPXUU4899tiGDRsaGhoOOeSQL3/5y3Pnzq3IkWOUYyjFqQfWPLXOAuEVLraS7YUn\n70dyvTw1tFLL9dqOp4QzTiiA+kx/X21+/dTBQdGz1ZC1+tL5G5WY4AYAMJNQgwOwLKrrDIBpUUNn\nAGyLajoDwCxK9aHKrbAodAaAWITrHIBmEtvgAHQTljH0lFqU5bekXBzB0qDbAGDr0CzOdEItALat\na5oFwGSGQU0AOdtIaSaADDNqqFhI1dAcgEFm1FITwABL1dEcCr+tYrmfp+pJDgVNaigsiwXnsvjJ\nrkdeSPuh18OSywDkUxR++usca6Qw1PIiN6NLQmrh43UMIkTPXBLoK0XwqFF4eCTydQCG5uUTLjuG\n/fVyvxeKUURFSJXnhhtu+MY3vvHMM880NDT09PTcd999Z5999osvvlj+kWOVYygJ0ZlGOaEAgU0I\nB6QYBBISHrk248y1LiROMnJDg2QraAQf2jJ4cFLUwkKofJwUkr5DtFApJIPnSoJ583goVV/VGTyV\nfPgfwXPwoMCopBr5hkcZZlhcY5zyvBoRwsmp+9b5XoyiKlCCVGFWrVq1bNmy/fbbb/ny5ffcc88T\nTzzxy1/+0jTNq666KpPJlN4/lLjlGMI5af86QjkhnBMwSmxa0Zsh1+9VrxBNqmz1VUlJIzhGQ5Pi\npu9ClkvKUngeL39JkcsClSToOK7zeitQ+KpR/iWeNyjmT+GYP54PjUCiYIQwsp8ahFTNKEGqMH/5\ny18AfO9732tqahJrFi5c+MlPfrK9vX3lypXlHDlBOYZwTppXTwin9P9v7+zjo6ruhP8759w7eZsY\nogiEhISKgGC13bZCW6m4oQW11hZEF634qKyPXV23q7I2gqz70Gr1o59VQddqsa2o1LKVWCtYurzo\nU/UxYyufSmlFECMwQCAoCXmZzNx7z/PHmbm5c9/mzlvuvPy+n3zo5MzMvSepyTe/3/md3wGgoBIy\nKOfhJzkXThoezFX31fSb3bk4Se8Lbj373Lpt1n1JKWWlg+1jcF5YchpxMlN8bpn2V7W+zHoj235I\nTqdLWA9AEg8iWiCiyVFVTjT5pmIBiaoET0IqaopYSH/7298OHTrk9Oyvf/3rhx56aCTnI+jo6GCM\nzZ492zg4Z84cAHj77bczvmxm7Rjcaa6V4ZSYqPxWKdESpQ2KnNOkR9ZOykvRXZrN7sDZSWAJlWKa\nnJP0nZcCPNsMHtiFR9YRMAjDSU46afVXBbskoVNloNVG1mSdbXgEACpner4ONEJVOL8FS+yKmCIW\n0k033XTxxRdv3LjR9tlXX331pz/96QhPKRaLHT16tLGx0bTGM2nSJAA4ePBgZpfNuB2DO82nBKQI\npZSLHbIqjf8K1igDAGPRQbaUn5MgF+k78LCqJD4d4HrWTrY1ExiwNRMkyymln6y4vN1licu7jYzh\nkX6p4Xwdp0QjLEomYGfVYqaIhQQAkUjk9ttvf+SRRzgviHMie3t7NU2rq6szjYuRnp6eDK6ZTTuG\nlMyaXCnqGjgFlRDwUteQGTl10vBFcu0kqgGMlJO8p+9SZvA0IOY9pKkCJqdBHZNgUn5Yr+AUpRkn\nqX8tpge2DA5/9+SoKsc0OZ6v04hYQPru9JwltBFfKG4hnX766VOmTHniiSduueWWgYEBv6cDsVgM\nACTJ/MMpRsSztoTD4VAoFA6HTeNZtmNIyZfPrKKM00SElI8GQsOkcpKRgF05n23RXW6dpLdBGwEn\ngef0HTiESmCnKKeSh8SIjYF0c7j4yTsOuUFmUpGLjTyGRwDAExES0YCq5POjoLu7W/zURCIRRUmj\nch0pBIp7H1JdXd2vfvWrpUuXbt269corr/zJT36Sv9/dXmCMgZ14xIh41pZwONzW1hYOh0UpRGNj\no3jQ3d2dp9hIcP6EGsp6KOVAQaG0Esy12gDAA7XuDb95tM80QgIOdYCWLUrG/UmmzUk6wcGT7puT\ndNLanDS8Syl5f5JpcxJAGvuTAMDjFiUAsO5SAgD3jUoAYN2rlPJT074lSN4AVGPZb2TrJOvLnF7p\n/hqjIJ0aJoGdjUzhUfxtifo60AjV4PwJNZJEdBWJulYpgfGxPoIUFEX/f0l1dfXjjz/+8MMPP/nk\nkwsXLly1atWMGTP8mkxVVRUADFqSS2Kkyjn+mDFjxrZt2wBABEnhcDgcDre3t+dkA6wLzaNkIJxR\nTijnlOh5z2hFrW2MYsVqIzGYvZMCQyf1xnruTtJ3y0LWTopPI59OAgCnnbOQvHkWEirS/QSG/bOQ\nkZYg2UxgEIatcqwv84jt69OykZHB5E7qABDTZFWVRH2d2IFEFZjWVG96ozCToii6ohQDTpZCXflF\nKXzTCSG333775MmTly9ffv31169YsWLRokW+zCQYDNbW1h44cEBVVWM81NnZCQANDQ0pr6BHSADQ\n1taWVxsBQEtNxQWTK7ftUggFjYJCCQAMVlRXDQ0AgBqoZekfhiQQorLXUp6dpJOBk5yaOEAWTgIA\nxhSRu3Pp5gAApoYO4CFUgjS1BB7MlHg2vXyXi7Gs3YBcuvaZ9sCamtclhUcAIPbQcUJVcs1ZZhtB\nIluu28W0h083E7i6CgDETyKKKt+Uzjf3W9/6VktLy80333zPPfd88MEHLgs2eWXatGmhUGjXrl3n\nnnuuPvjuu+8CwPTp032Zkjtfnlz52t9OUso1ShRGT1ZU1A4N5erijqGSXXshHS9O0rE6SQ+SbHFz\nEgA4NxaC9J0EADlJ34EhVIIstATezJR4Nge/HFxUZPvY1kaDw7u44lcT4ZGwUbyiQSWmo4G9YBSM\n1VVgMJaeBlQUpbe3N90bIR4p7qIGE+eee+6LL7549tlnP//88x0dHb7MYe7cuQDw4IMP6iOHDx9e\nt24dY6y1tdWXKbkza0KQShpjnFBQKFES/Ro8nUPB1dQvifbZpvVMNQ4jVgjuhocGrF5rHLyVOXis\nvoPkQMFU7GCqd3AqeXAaAUOxuLUIIgOcLmWdm2nc3UbiKWN4xDkVLeyIRpgK5zflpoOJQLiqsrIy\nGAyOGjVq9OjR48aNGzduXFNT0ymnnJLDGyFGSkpIADB27Nh169Zdcsklfk1g0aJFZ555ZigUWrBg\nwdNPP33fffctXLhwYGDghhtuGD9+vF+zcqGlNgB1EcY0wjhQUO16CA33RTVZyvPZOQXlpMyaOLg7\nKa4lhaVVegcequ+cCvDS1ZLtiG2ZtclP7pZK+Uqnu+uP9S9B/9LAkqkDw/cnpsmqygBAbxdENEJV\nWPR583YLpOggBbKDJwOi0SgABAL2Gxd+/etfd3V13XLLLSM7KepNkm4AACAASURBVACA7u7ulStX\nbtmyRVVVAKipqfne9763ZMkSlyo7W6ZOnbp79+78zDGJS3+9740/aYODkhKh8hA/ZSgaHBqqifQJ\nDQgxDNsi2SL2pnHAPn2XLDmeXGhnrLvTgzZj4k4vVTcW3enrSXruTkok5eTEA5Z4QC0P9AIHnjjX\nVE1cR0lcWR/RDOk+PnyFROzIEqak8QeMDa/NyIm1IpG+E1QmBvVKB0ik7wTVhnHrpzUk6dNg8qe2\nI+7jmWGrOseOfJbeqdbwSF89ipczaBJoDBRGYkwaYtdPrF99yQhV2HZ2dk6cOHFk7lVuFPEakpOK\nBAsXLhyxmZgYPXr0qlWr/Lp7BpzfXPPWjh5KORE/4JTm6Y8U+yWl5PUk00EVtqQsBM/slAqnQnBw\nWk+SOJB4nJRBmQN4qL4DQ6UD2K0qWT91X1sCh/UksFNIuopy2dbq0qTOo43EU8bgEjQKGhXh0fmN\n+dnQjYwspZayQzLg6mn1VNYkxkXWLsaoyih4XEZKk3Rzd5l1BDeSj9xdvC+4QtLdNus9fedlVcm0\nsJQyiWe7nuSUuDO9wMuHyxVMI9ZJprSR/s2Jh0cqEweWE5UwBb7ajKdOlAIoJARaagNfmgCMaZRy\nTkFhNCLJeggikmaOy0jpMzJOMp5SkXMnQaatHMBhScna0AE8rCpBmloCBzOBBzmlhYv/jHMzTtv4\nwNRBFRLflnhxXZxEQ1WNUIU012MLu1IAhYQAAHy+kVCm6Vm7aJrLXeniZeUpXScND+bZSbkuvYu/\n16n6zilUMmnJxVJCSx7NBA6BkfVl3l9sVZExMNKrGCwdGeyTdVyTQJWM5QxXn3Wq0/SQ4gKFhAAA\nfKtFprImMY0wzinEGBUHyIqsXTxIknKZFbFxkmtviCyL7oxk6STIoBzcrfSOuafvwCFUgmQtcSDu\n0RLYBUzgOTbynqazvsVpAtaufbY2Mibrhq+uUdAIVQlTYNG52FO1REAhIQAAjVV01uRKJnEWX0Yi\nvdYmETTxW9VYg+DUIsgDKZ2Up0JwyMZJiXqP9JwEaSwppRUqQZpJPHAImATeZWOLy9utKkrLRubw\nCADi4RGh/ez8FlxAKhFQSEicL0+qZpLKqEYY1xjRj0fKR2mDTrpOMuLeETzpIs4nzEK6TjJ0BEhv\nixKkXlJKN1RKV0tOZrKVE+SiqMF6fevE9PnrX5fpKxXfjaTwSGGEE6IRqsHVn8PwqHRAISFxrvrs\nKCJpksRpImvXU1npWNqQO9JyklNvvZRFdy4bZiFNJ0Ga22bB85ISpBkqQZpaAgczgQc5pYXtpVxC\nNycbGZJ1DIzhEY+HR0yBOy84PfsJIwUCCgmJ0xIMfHVaBWUaY5wwiDLKiXOQlKOsnSN5KAT3xUne\nl5Ss6TtwDZWcFpbAQUsezQTJcvIiqpSvNN1oUJONKnK3UeI7QIa1rTC9nOFrY2vxiNhSAoWExFEU\nZW4LsIAaL21gEGVUobS/ska8YESDJABQk3q8FqmTwPOSEjin72xDJXBdWAKLlsAuYAKDmWzlZCRd\nUVkvKz7licLCQcvkY5zpXylYk3WJB0QjRAWqwlXnYL6upEAhIXG6u7u/Nj5IJC1e2sAgymh/IMCJ\nw38kOQ2SbJykuLUJyKwQ3Iitk3RSOClR2uDupDTKHFzTd2AJlTwuLIGzlmz1Y5RTSj+l9XZrkGQK\njMQXonIKdjZKCo8UmsjXEayvKzGKuHUQkkOOHTs2YcKEz00YM2P6yf+3Q5WYpjLCGYkyGmMMko/s\nS3mGrOCdv96lP66Q6z/TuHAo+unoUV90er1NYyHXrkIpT6nQuwpB+qf5gXtvIZUCQMr2QsZB44kV\noHe90zsMgeO5FWDoMwSG9nembkNgaIKn/4o3nvsnHhh74hkNYWqIZ31BZtheYdAu2WibqUuykSrF\nFc6BaEBVuNruACSkqMEICYEjR44AgDgr/QdfHcMCiiRpVGTtJDooy/2VcU+Yjxh3DZLOnXyn/ngo\n9ulH4V8fOrbVfSZ5LbpLd8MsuMdJ4Km9EHgpc3BN34ElVPK4sATeAiZBNoGRy9Wst7am6cCuisGM\nnqxTKNEIUQlTYdE5KKRSA4VU7pw4cQIATj89Xqp0wbjar55ZRSVNYhqVuMbIEKNRSYJEaYP3laQK\nud7kpKHYp+FjW1K8jWvmkRwV3cHIOgm8lzmAY/rOdlUJUmXwPGrJi5y8WCrli603Ms7HaiOnZJ34\nX6oRpsLXGoK4/aj0QCGVNX19fX19fSI20ln4xToqa0lBkiTpQZIZ5yApfGxL94k/mXJ0x0+86z4l\nHhtI60tIq8ABAEjivBWPTRwg106ClLuUwC1U0rT4y6yhkhctRRxsYWsmI7aWcneV7ZVNgZEnGxnD\nI5UQlVAFrpqO4VEJgkIqXyKRSHd39+jRo03js5uq4NR+JmuMaUQESRJTKTV3EkoVJB06tvXQsa3d\nJ/5kHPQSJOW8gwMYnFQTHT7l1mMTB8jISR7LHMAlfWcXKnFO3DN4Vi2layYvfnLC5QqmmxpVlNpG\nChPfJaoRqgJVySLcD1uKoJDKlEgkcuTIkXHjxlVaWgQ1V1fd/pUmKmuyxBnjnMGQxE4GKhyvlWa5\n3fET754c2Of+mrw6yb0QHHLkJLArvYOUS0rgNVTyriWwC5jA2UwCk588flivY72LU2CUbCPDN0FP\n1ilUrB5hOUOpUsQnxpY8eToxVlGUI0eOjBo1Khgclof1EMybXv/4xe2DsUE5MsSUKCVRqB1SRg0O\nBgdPUk2VYgMuJ8kadXJyYN/7nT+1TmP86XMaT/+6+1Rt9GbZomsM1GyPlwWHE2Yh1SGzYHfOLLgf\nNQs2p82C4XhZsDtzFpyOnQXDybMwfPisQD+CFpJPoY1PmJozksZDaXUqLS+ze03qY/qcrGZ4QZKr\njL7UVaqaat9NsZFCmUKkKFwz+dTVlzamnFL+wBNj8wdGSGVHd3d3MBg02siWZV9sIJUqk1V9MSnG\n6Imqqr6qWo0ycD0nSReJk40gR3GSiXTXk3IYKtl2Yk25quQpWrJL4oFrwASWmAnswiZwyOlZXhNI\n+eH8XvP1TdNIaSMxSBUKHKgKNEb+bdYYl9kiRQ0KqbwQFd6jRqXOv7cEA5d/qY4FVEnSmMSJBBGZ\nagQ0QkCvuJNrwHkxycVGkCgEz+RryF3uDlLtmQXPTgLP3RzAOYNnX+wAbkm8zLTkYqaUfkqJ03Ws\nKvJiI/37wBRCVbj67PrmUdgrqGTBjbFlhHG/kRfu+vKYF/94gkWZpGqqSjRGhhgbCASYVlMT6Y8F\nauRof9IbKmp1W5BA8P2/JtnorIk36n7ykq8T2OyWtWDaMGvEdsMspLNnFgA4B9HVz33bbNJjheq5\nO+POWUhk6mw3z4Jh/ywAEIUOp++Ek0QGT/y+TqTpjBtpIXkvrUD/vW/M4+lusE3lZekkEyb/GR2p\nWmM+i42oQlmMUBWoQhadjatHpQwKqVwQ+4282wgAWoKBBTPqNrzZJ6lUUQiXqKKSQY3JcgAAaiL9\nAKAGaln05HDvBoOTZnx+tZeTYTMhuYODCWMHB/DmJB1bJ6kq1ReTsnQSJDd0gMSSkouTwLiqZOrs\nAGloCSy9HgRGW9jKKWNsd7na2mg4zlPNv5H0qJFqcPX0+vMn4t6jUgZTdmWB7X4jLyybOYYHVCar\nAVmjjHMJhiQ2KMtiGcm+CjzXjcDtreY5cecF28UkSNWDFbzk7lyr78B5ScnTqhLYLCxZk3i2ebyY\nXQykJ/Rs03opcXm76abGWdnYyJSsixGmAlXInbPwpIkSByOk0kfsN8rARgDQUlOxYGZd+xsnmUIl\nlXKNaxqJaFSSZVJVa9vJ1AQJBHWjnDf9xxnMATJK3KXsdAcOiTtI1e8O0oqTYDhUso2TwHP6DpxC\nJRjugxf/diVHS+AaMIFdSZ4gAydZsZrPJk0H9jYSyTpRznD1WfUT6nLQ1ggpZDBCKnFc9ht5ZPl5\np2uVKpNVWVKZxAmDmEQjsqQR4rhVNjmflpcDk8Cm4o7E+m1fCFkU3UH2cRKkUXoHOQqVwBItgUPA\nBIbwxTZsygCnC5runtJGYoypQGPkzguwuK70QSGVMmLLUTY2AoCWmspHvt1IKlQmaxLTqMw5g4jE\neiviW2W9OCl7vCTuTH3wTIm7bJykk72TwLmhg231HViKwrPUEjibCZJd4l1Rqd5CrCryYiMWI4QD\nVeGxi5vwIL5yAIVUsggbjR49OhsbCf6+sbJhnMQCqqT3E5IgIjP9jPORWUzygvfFpOw3J0EunAQe\nKsLBLlTKQEtpmUnHqijv0kpcf3huZhU5rBsBAOEAAEyBidUVeO5RmYBCKlk8boD1QlN19WPfHM9F\n4k6Ob0uKSTTGKCT3QRA4OSlbLXHVZtCauEuWEFWHGw3YHlEhyGxzEuTfSU6hEpgyeJBaS+AQMIHB\nTO5y8oLTpcy3Nj5rnHbiy6QKkRSgKll9iZ99GZCRBIVUmnjfAOuRC8bWfml6BQ2okqTKskYlziWI\nyOyTqioA6K+sMQZJkJ/1JB4bTP0iCyT5KPSUx8ua3263mATenOTe9Q5cnZRhqATmX+5xLXkLmOLv\nMBglpaW8vNLmXnY2MiXrAIBocH5jzfktNU53R0oMFFIJku4GWI88eWGLEuCsQpUkTZI0KnFNgiGZ\n9VRWimPOTU4aJt81DqmCpFwl7iBNJ4FLJ1aHMoeMQyW3DF58Wo4Bk4ucht9tZyn3cMr+yqY0XSJT\nZ7WRFAOikscvbnKfGFJKoJBKjQw2wHqkJRi4Y1692JYUr7iTuCKRiMROBgJui0kwUnV33jAl7kbS\nSfkLlcBJSx4CJoF3ObnjeBHTfZMDI7DaSIMffHnMBGwUVE6gkEqKjDfAeuTaiaPPO6uaVihM1oST\nuARRmUZkKUbZSDjJep6sIM0gKYPFpJw4CdJP30GqUCnFwhK4BkwOIY5RTqaPlC+zuZxVRQ42SlwX\nAOBrjcE7v4al3uUFCql0cDpwL4c01wV+8o1GpUoTFXfimHOxmNRbWaESIpwUC9SANyelq6V0z5M1\n4rEK3ISPTnJK34HFWClCpfiN7LQEKcxkJY1AyvbKhjlYbRQPjxSgKrnzK2ijsgOFVCJkvwHWIy01\nFXe0nsYr4hV3YjGJSxCRaW9lJQD0VdWKJSXw4CRIU0tpCYyk00DIKXFnvmb+neQ9fZdeBi9+L7s8\nXnyuUrpySu8iyfd1tFEMCIcffHkM1jKUISikUiAnG2C9c+2UU780rZpUqCygypLGJI1IXJVIRGYn\nqqoAQO/gAN6cBDlZVXI9KkmQfeLOjOF4y1w5CTyn76zPOmkpPTNBsldMH+6vscWiIhcbAcDE6oo7\n8dCjsgSFVPTkcAOsR5qrK5+cMz5WrdGAKsWdxIkEikSGJKp3h4smi8eLk0ggSGTHds4ZSCvjIAk8\nLiYpST7I1kme03cuoRLYaQmcAiZIZSYr3qMouysb52CyUWKQrL4YNx6VKSikoieHG2C901JTufpb\n47VKVSwmybJGJU0UOHBCPq2qMhU46HESlxN5mIpa+/ZChMbNlNAPkatysKk2gXuQlOVikom0nQSO\nS0rgIVTKSktg8Id3Oblfx/XWVhvFk3VfwWRd+YJCKm5yvgHWO4snnfbFaVW8QhU7k2RZYzLnEvRX\nMAAwOQn0/UmEpgyVdOIeIln9frQJkpJL9Twm7iD9AgfIhZO8h0rWF4CrlhzNBMlycrGUt5eZVaRS\nWxsBwNcagpisK2dQSEVMnjbAeufpv28cM14igfhikl7g0FchceLsJGv6LledWD0sIwEAc24KDq6J\nOxMenWTExUkel5QgVahkfQE4aAm8mMmIF0XZXd80E33hzWSjidWB31z1GU8zQUoUFFKxkr8NsN6Z\nUF315LyGWJVGA4rYLStJGpU5l6E/ICmUmLqvgpOTIEfdwTO9iHvijvLh0gWXAgcvRXfg7CTzp6mc\nlG6oBAktuZvJq5xcsb2U6dYmG7EYeewSbMpQ7qCQipJ8b4D1ztfG1N0xp16p1miFygJaciE4izJq\ndJJxfxKIJaU8hUrJWLN2aR0pWxeJGD91OaIi905yTt+BXajkRUvgHDAlbkTT9VPKt5huZ7IRANx5\n/pjzJ+DSUbmDQio+RmADbFqsOKfpipk1WpUocIgXglM57qQYo3rujhNq7XdnEyrlR0vuuAdJHhN3\nJnLgJEgvfQcOGTwXLbmYKXFHmvLD5e3WwMhqowsagneej0tHCB5hnootW7a88sore/fuDQaDU6ZM\nWbJkSUtLi/tbNm3atH37dtPgaaed1tbWlv18RmwDbFr8n78b81Z43+FOlfH41pwoAOeUAxkEVhNV\nhk8KHzwZragNDMXPFxeRCg/UmoMY4SRva0JxnOvFc47xvHNwOPJc4HTwORjOPgfTkeeWT/VD0OO3\nMxyFDgknqYabil/0SnITOOvLdHRhaJJDZ6b0sfWc0Yv6ceoXNAR/swiXjhAAFJI7995779q1aysq\nKqZPn97b27t+/frf/OY3jz/++KxZs1zetXHjxm3bttXUJOUfxo8fn/18RngDrHcaqoIPTyY39Sq9\nx4Fx4EBkAPHXPAfSD1JwSPm0qqp+cLCvqlZ3EgCogVrdSWDNremhUkozpR9UsehJU2PywNBJ496p\n4OBJ41FP9YODulbB4iQjI+MkAMhYS5DKTJCpnFziLdtOSBOrA4/hriMkAQrJkY6OjrVr1zY3N//i\nF79obGwEgM2bN992223Lly/fvHmzixJ27979hS984fnnn8/tfEZ+A6x3du7cOY7B9kXTzn5md6UG\nDAC4BBzE38EcSJ8HJ4GTliDZN0MO0soDhGt6GySwOCnplYYgCSxOMpKVkwBcQiUAYDFiMo2tlsA1\nYBLYqsVkKaIRojluwzLdyzQlAAAOj13SNKEukPIKSJmAa0iOvPjiiwCwdOlSYSMAmDdv3kUXXXTk\nyJE333zT6V0DAwMHDx6cNm1azufjywZYL+zfvz8SiZxzzjktNRVPXDxuKMhJhcYqVEnWZEmVJC7q\n7voqpL4K6dNEbyEAiFbUmnbOCsyrSibEIpP+kU9qIq4F4skVdx4XkyDVepJbOTikWFICu1UlsKt3\n0F/stMJki77sJD5S2sh6ceM0XrniDCxkQIygkBzp6OhgjM2ePds4OGfOHAB4++23nd71wQcfcM7P\nOuus3E7Gxw2w7nR1dXV1dZ133nni08WfGfODC0cNVWui050UEKdUaFTiXAYgMBBgRieBXUU4WAvw\nssP7pVw2yQpctsqaMHVw8O4km0ulcpJ7UbiOk5YgfTOlxPZqxrvfMBEm8ONdXV09PT09PT2R5FJG\npDzBlJ09sVjs6NGjTU1NpvzYpEmTAODgwYNOb9y9ezcAjBs37sEHH/zLX/5SVVU1derUxYsXZ1MU\n5/sGWCd6enr27NlzzjnnGAfv/uwErrGHXjseAEgceiP+pRohmkIGZQZgWIxJTt+BoSZbF0laLely\ni2klyR33xJ1xMQnSyd3ZjCQvKYFD+g7sknJOSTzju2zf6AUXpRltdPX0+js/LwsPHT16NBKJDA0N\nVVRUAEBlZaV4UFdXp39agGlqJB+gkOzp7e3VNE38SBgRIz09PU5vFEL613/916GhoYkTJx44cGD7\n9u3r1q179NFHv/rVr2Ywk0LYAGtLT0/Pzp07zznnHOt3acW54wHgwdeOV3BgBACAAMQIKEA1IAoh\nA4RVR1XTkhIkAhSTlsBlbckDOYy0bHGpuAPXxSR3RFTkvcwB7CodwG5VSeCuJbBTi+11PAZVRhvN\naqp5fF4TAIwdO9b4GhEkDQ0NCT/19PQMDQ2JcZOujJay/ueHFC8oJHtisRgASJL5+yNGxLO27N69\nmxDyj//4j9dff31FRYWqqk899dQjjzzS1ta2adMmpxWgUCgUDocBoDGBGBcbYJuaCm4HeyQScbKR\n4IbJdR+fjK7/08kKAgyAJNrFqIRqhKhA+khSmQMkQiXCNTnaD8nFDoIMAqZ82MiltCEl3oMkW1KW\nOUA6oRJ40JL1OuliShLOaqr57RVn2L5SCKaystL6H5XRVQCAoVWpgkKCpUuXij/EdFauXMkYAzvx\niBHxrC2PPPJILBbTi7wZY//0T/+0b9++l19+efPmzZdffrntuxobGzds2CCcFAqFxAgAdHd3F2Bs\nFIlE/vjHP06ePNnlj9PGqpoVn5c+7ot2vD8UAM6AAQFCeIyAAkQjhBPSF//PrwoA9FBJbJ41hkpg\naaxg0gxxfTZ7Umbt0gqS0nJS6twdeErfgQctgWczpcH778Ck8/TPXGzkjslVHkMrkcnQIyqwC7CQ\nggKFBFu3bh0YSDoY+6677jrllFMAYNCyP1+MVDn/gXz66adbB+fOnfvyyy/v2bPH6V2NjY3333+/\n/qkwU2trawFuOQKAPXv2NDc3m34pWGmpqXh2dsNiONzxwVAAVBZfSQJCqEKoSoATQlToq7APlcBQ\nYmBN4hnJzEBqFt7KJkiCkXISWNJ3kKraO4dmkmKgfPhO3+NLRv3ne2IkYxulJCehldFSqCu/QCHB\njh07bMdra2sPHDigqqoxHurs7ASAhoYGp6tpmkYIISQpuSEydcePH/c4JT1lV4A/FTt37qyoqGhu\nbvby4nGVtc/OhpsCH27/CwsQlREAwgkAASCEqgQ0QsDgJDCESpAowPOopbTIxkZeyHIlKSdOAodQ\nCZwXlnSyMZN4r/bJob7Hl+iD+bOROylDq3TTgNu2bXv44YdbWlrED+nMmTMBYMaMGfrPLJINKCRH\npk2bFgqFdu3ade655+qD7777LgBMnz7d9i2dnZ3z5s275JJLHn74YeO4qHQ488wz8znfkWDnzp0A\nMGXKFO9vGVdZ++RXJv1v+HD7LlZBVAacECCEE8IJoWJJiRPoIxJwsIZK4KAlyMJM+bCRS+MGQVqJ\nOxgRJ4GHUjrjCpC7nJJe+eE7RhudP6Hmtwt9sFFKKisr000DfvDBB11dXbfddpt4WUdHBySWfkd6\n9qUICsmRuXPnhkKhBx988NlnnxUjhw8fXrduHWOstbVVjHzyySfvvPMOALS2tsqy3NLScuqpp27b\ntu3999/XtyKdOHHi6aefliTpG9/4hi9fSK4QNjIVeXthXGXt6i9NX3vasQfeOF5BgBIAAhIBQoAY\nlpSIXagEDlqCZK94lJMq1wBx2/ETzd1OW1OQlBKrk0xk4ySwS9+BZy0JnPYwmYi+85uBX67QPx38\n/RPnfXksLLzH05sLBts0YDgc/tnPfvbss8/OmDFDjMyfP9+f+ZUohPNMalLLgVgs9p3vfGfv3r1n\nn332N7/5za6uro0bN3Z3d994441Lly4Vr+no6Lj22msBIBQKif9wN2/e/C//8i9VVVVXXXXV1KlT\nDx8+/Pzzzx87duz73//+zTffnNYEpk6dKkKrEaCzs3PixIkuL9i/f79xA2wGfHwyurbz6P1vflLR\nT1iE8hjTYlRVqKJSRaGaRrhKuApEAUhsODXuQjW23065fRUAgGtAaPxfzzgJyVrUYF1DskZIJiFZ\ns3ZSsjysQrIW3ZkNZDcCYC69i8/Hzkk6mW08MhHZ/ERk8xM2ty7+3zPhcLi1tfX+++9HCeUPFJIb\n3d3dK1eu3LJli6qqAFBTU/O9731vyZIl+qqSVUgAsG3btoceeujDDz8EAEJIc3Pzv/3bv2UQHhWO\nkLq6uvbv35+NjXR+tOvAfW9+UtFPpQiFGFVjTFOoqhBFpapKNJVwlYgePEzjVTEVnLUEHs3kGe82\ngoyEBBYnSRZt+OskyEJL2ieHBl64W9n7R9tnS+D3zOLFi2fOnPnP//zPfk+klEEhFS4FIiSXDbCZ\ncd9fP/rhmz2BASYPEhJjWoxqigiViKpSTSWaQiDRGk7vzeOiJcjaTLFADXcOpPInJMhRkOQ0OGJO\n0j451Puji1xeUOy/ZxYvXmwqhUXyAfayQ9zIuY0AYNn0z/x83unRGnWohvMKlVaoTFalgCrLmixr\nkqyxACeMcxm4FG/JCgCfVlXpDuirqjVJQvRp1bu1eke8JV0beSX9X8IxS1M7a5s7c2u7NLH2vjOR\nblM75cN3nGw0ceJEznkJ2AgA0EYjAAoJcSRlO4aMWdTcuPvqsxrGQySoaZUqEU6SVUlSZUmTJJVJ\nnMqc0HiQ4aIlqzCMchIf4uh0SByj7t1eTjbyugkpZ61KU2NvKbdzylM36vZ4a2NBnQmxTaKoETbS\nK5vyQU9PT2tr63PPPZe/WxQLKCTEHi/tGLKhubpyy8VnfvaMTweCWqxK4xUqDagsoDFZi4dKkkbl\neJtwFy2Bs5l0RPcHoR+XYMhEVrGRA5rlt7ziIdzxGCSl66SU5LD5d5Hy2GOPhcPhvNoIANasWRMO\nh/v73c46KRNQSIg9HtsxZENzddXvLjx3+az6SK0areZapUYDCguoTNZkSZUlTZY0JmksED9RyVZL\ntmbK0iXuV8imR4MXrFm7bHE4tShlkARpOunCCy/UH//85z/3/sbCpL29fcOGDdu2bcvHxVVV7ezs\nfP3112+77bannnoqH7coRnAfEmJDWu0YsmFUZf3dZ9dLVLmno0cdYBWMsCgwqmmUEYVTSqlKVA1U\nlYAMAKBBvN5BOEmUPOiGMBY+GI3CNLVqKKk7lC1eNOZkI/ddsTokVUGBd2w2ITm+lAC1v6/Ttlmd\nlAUO+m7Z7du3v/7666+99pr49LrrrgOA66+/3tMMC49QKNTW1pa/2Gjfvn2XXnppni5evKCQEDMZ\ntGPIkrZpZ3y3ZXDO/3xwoIsHIiBHCKUcKCUKJ5RSlVNKVZVoGgEZgJu1BK5mAgCVsuzzby6BkUcb\nAQDJcw4sDUtljSIBEGg+JfBXgO3bt1944YWvv/668QXCScVIKBRavHixcQNszmloaFi1apV4/N57\n761ZsyZPNyousOy7cPGl7DvjdgzZc3iw7+l9x3/U8YkUbLFypgAAGCdJREFUkQIRYFEKMcJVqilU\nU4mmxZ0kPsCwHqOfHU44r4mqpsvWWzrkpot7js7JRrZtGqxl34xxQsyDXoq/Iev67/g8M4qQRGB0\n9fRRT7ZOcHl7MSI2wObVRia2bt16880333777TfddNPI3LFgwQgJGWb//v2RSCQnG2AzoKEqePfZ\nwTljgxe++pHK5IAE8hAhMc4oJwqlGqeUaipRNaJpRFUJFUm8GBG/+kkMOCGmmAksOvHuJy9rRVna\nCACsNiocbG0kVNR8SuCJrzdeMN7+fK/iRW/HMGI2QoygkJA4XV1dWTYHyglfGX3ahwurfvHh8R+G\nulUmyVHOopwyzhVCqMjgEU0jjBJVI5zHgyRdS5AImKxmEuSqJMElTZdWCzsr7h3tcotLeGS1kb5c\ndH5z9e8um5S/WflIW1vbrbfeis2B/AKFhAAAnDx58uDBg75k6qw0VVXf/dlqSYo+sONwpK9WilA5\nqrIYJQonKucqISoVfcLFWR+cE0K4kJOLmcAipwxIuVzkZCPb8MjaOsg7aeTr0oQzrlkK/RKBkXzV\nZ0fd/YWCOzQyJ4h2DLltDmR7/md9fX0Ob1FKoJAQ6Onp2bNnz7Rp006ePDk4OCgZAP/OZGo7a9I1\nLQ2/2Nf9w3e61SFJjmhSlFCFE4UQyjklRKNUI1QlGhcLS1xVqfi9zzXgqtlMYKcTvWOeE/0Bxr2V\nIrgERt5tZBseuZ9rnhqHt9uGR06BESdwfkv15ktLMzCCvLVjsD3/E4XkBAqp3BHtGGbPnl1XV6co\niqIoYlAxIMxkUpRRWnmiqar67rObzzoFlr37yf7jXI1IUpRLUUJUThRCVM41QinVNMI10DRCaTxO\nUhRKEoXOxpgJDHISqJR4r5GzxT1B53QuX5Y2yr6cwYuNdBU1j5J/0tp0QUOprRjp5K8dg9P5n4gt\nKKSyRm/HcNpppwGAbhdxxK2OsJTuJ6ErSBxiZjSTlExOJrlwQvOs0+t/tu/4D0PHlQCThpgc5SxG\niMKJSjjjRCVcI1QjVCNc5PFkVTeTrgRRlWf1B/F2zI+JlAtFLkfE+msj7yoCAE3id80YU6o5OoFo\nx5CnDbBIWqCQyhqP7Rh02VifMroKXEMrq7e8z3NcZe2y6bWLJ4555qNjP/zTcXWISVEqCS2pnKiE\nq5RrJG4mLrTENU70X/1GM0FyC58saxBMuB9V7rRo5KON7KMiAC7xWU01T104oSUYsJ1zaZDXdgxI\nuqCQypectGMwucoaWuUwDTihuvrus1uu/cyYm995/38+UpUAk4aoFBvWEjDCVcI1g5Y4sW60M8kJ\nEmtO2Xwf3D0EaaoI8mAjjyoCAE3iE+qlpy5oLuEcnSDf7RiQdEEhlSkj047BaBePaUD30KqysrK5\nuuqV2X+3/owDP3rv0O5jshJlLEblKKcKoSohKhCNcnH+LCdEI1wDzoksa+IwiJhCjW4QvU0JBeLQ\nXCfr74BjPUIBqkiVeHOdfM20UXd/vsF2bqXECLRjQNIFOzUULvnr1OBjOwYvWNOAYFCX7qoepf8X\nBz/91SEePslolLEYlWLAYkBUQlQCGgEtHi1xsWOJA+dEM0RCtm1MvbTfdiFlJbfLNiOnajqvNvKQ\nozOqSKOgMQAATkCTeHOd/N2zRq0oAxWBH+0YEC9ghFR2+NuOwQsp04Di36AS/EHwtOsn9v7h04EH\n9g58/ImqxCiLURYDSeFUIUQTAVPCTJxwTigF8YBzCMia8e8x4ads9ga5kCsP2Y+nGRUNLxQxzilM\nqJOumVa/4nNloSLAdgwFDEZIhUs+IqSurq79+/cXso0y41ikd/uxnmc/Ovb7jzmJMRqjLEaYAlQF\nYSbQADQCPB4wiWhJhE0C65lDkMVJEJLkcOSDAZfdRTlM0CV5SIqfGcgpaJRzypvr5GvOGrXic+NT\nTbakWLx48cyZM3O7ARbJCSikwiXnQsrHeeSFxpvdnzzz0bFnPuiFqEQVShXCFGEmPVoSKSrCOQEN\n9FRe/N8Etn7KEvf9rcSuyyrkaKHIWLDAKXBJba4N3P2FsYsnjU4979JCtGPA88gLExRS4ZJbIZWD\njXQODg4823nkP3d29g7VkRgjCmEKpQowFahKiAaGmAmAE7OW4oOO13fXlffGCi6dfnLoIT01Fw+J\nptZfe2Z9S40/DTj8ZQTOI0eyAYVUuORQSGIDbJnYSKcrcnL70Z5fdh753QENYhJRKVUoVQlVgGlA\nVKAqAQ5EaCkhJz1UMib0BPZnhHvGvdEcYZp9iyImjn9KwiUvB4nUHCfAadxDnPKWU8k1Z45ecU5j\nRnMvBdBGhQ8WNZQ+ejuGsrIRAIytrF3UXLuouSk8OPDsx4e3HR58PTyoKoyolKqUKoRonKpAVU41\nQjSuy4lwAJE/E0k8DhwIcGAB1ZTZyxhCuVuheZp7WjkZrpeLGyhuI95cJy+eUv+/zjy1ubocQyId\nbMdQFGCEVLjkKkLauXNnXV3dCJxHXvgcHBx4rvPwL3b3fNSvgsKISkWNOFUJ1YBqQOIfwkZA4uk7\nAIOcAMyRU7Z4bjenScM2FBKKT0cPhgiHeGpOWjz51GsnnVqeqTkT7e3tq1evRhsVPiikwiUnQhLt\nGEbyPPKi4MDAwLr9h7YeHnj98CCoDFRKVEo0QlRCNUKMcuJAOCFCGbqlxGMAwwPDIAEAbs2zuWPb\nYs5oIADgBFQGAAAUOAFOkiTECW+pkxdPHnXtpNEt5R0PGcENsEUECqlwyV5IBb4BthA4NNj/Rvcn\n//doz/Mfdg/EqkFloBGiUaISwuOFeYTrZop/6MGTUA8AmGMm46dWM5HhBzzxKSdcxDrxCxAAAI3G\nA6B4GCSOJKTiARcPGmvgq6dXnjWm5u6zMQg2gzYqLlBIhUuWQkIbpcvBgYH/d7x7S9fAWwfe/2Cw\nATQGGgEROXECGhAu5BRfZxqWE8QNRJLDJqOJjHFUfITER7jxX5J4StcPSbiK6v/yllPkC8dXja+g\nM6ur/+4Um35LMIKnhBQs2I6h6EAhFS7ZCGn//v2FcB558RIe7H/7ePdrR3s+6OXvHPukP1YtdpOK\n3bVEE1UPBOyiJTIc5liumySkuLVMQRKQeBUFxOMhDoSfMYo0BwOzxtSdUSt9t8W+pYKpPaCp91K+\nTwkpNPR2DHgeeRGBQipcMhZSqbZj8JFDg337BwY7+4fe6D6xp1fr7FG6hnoHlSrg8cUc0JJKHoh1\nYckwAGLUoCKIq0hsgIIz6qGxunr2uCAQfuGY2lmjT81y/h7bA2Z5SkhBge0YihEUUuGSmZDKagOs\nvxyJnNw/MPhR/xDn9KP+oY/7hyjR9vdHBoYoB3IiNtCvViTn7YiwUwWLjqsKylRprKo+IyhrQAjh\nXzt9FNfoBWN8+H/NGk75flhwlmA7hiIFhVS4ZCAktJHvRCKRI0eOjBs3TvzuLnac0oApTwnxcc64\nAbZ4Kdy/cZB0iUQiaCN/KTEbQY4OC4YRDK3QRkUNCqlEKNt2DIWDoiglZiN3PJ4SYjqAMa8VFtiO\nodhBIZUIe/bsaW5uHjt2rN8TKV+6u7tHjRpVJjZKSU5Cq7QqLNrb2zds2IA2KmpQSKWAaMeAzYF8\n5MiRIwAwatQovydSBKQMrdJNA7733nsbNmxob2+fMWNGe3t7Y2MjAODeo2IEhVT0iA2w2BzIR4SN\nxo0b5/dESgFjJOQxDfjnP/+5vb391ltvPXjwYEdHRzgcDofDt956K+5AKjpQSMUNtmPwnRMnTgDa\naESwTQOGw+GHHnoI2zGUBrk/FhMZMfbv3x+JRNBGPtLX19fX14c28gu9HQPaqDRAIRUrXV1d2BzI\nXyKRSHd39+jRZXcKeOHQ1taGqblSAoVUlPT09OzZswfXjXyk9LYcFR2iHQM2ByolUEjFB7Zj8J1y\n23JUgIgNsNgcqMTAooYiA9sxFAK45chfctKOYcuWLa+88srevXuDweCUKVOWLFnS0tKSowkiGYIR\nUh7p6elpbW197rnncnVBbMdQCBw5ckSSJNxy5BeiHUOWNrr33ntvueWWbdu2BYPB3t7e9evXX3bZ\nZW+88UauJolkBgopj6xZsyYcDvf39+fqgtiOwXfEliMsZPCLnLRj6OjoWLt2bXNz86uvvvrCCy9s\n2rTp0UcfjcViy5cvF52NEL/AlF2OUVX1wIEDH3/88UsvvbRp06YcXhnbMfhOd3c34JYj/wiFQm1t\nbdk3Tn3xxRcBYOnSpaKnAwDMmzfvoosu2rhx45tvvjlnzpxsJ4pkCgopx+zbt+/SSy/N+WWxHYPv\n9PX1RSKRpqYmvydSpoRCocWLF+dkA2xHRwdjbPbs2cbBOXPmbNy48e2330Yh+QgKKcc0NDSsWrVK\nPH7vvffWrFmT/TWxHYPviC1HGBv5RTgczpWNYrHY0aNHm5qaTDUpkyZNAoCDBw9meX0kG1BIOSYY\nDM6bN088zsm5L6IdA26A9RHccuQvuW3H0Nvbq2matSxIjPT09GR/CyRjsKihoMF2DL6DW458J7ft\nGGKxGNj9sShGxLOIX2CEVBCEQqENGzYAQFNTU2Njo77WumfPHszU+YtoDoQ28ouct2NgjIGdeMSI\neBbxCxRShixdunRoaMg4snLlyvr6+syu1tjYOHPmTADo6Og4ePBgOBwOhUIA8O///u/i2cbGRqOr\nsJXkyCC2HJkOQUBGjHy0Y6iqqgKAwcFB07gYEc8ifoFCypCtW7cODAwYR+66665shCQyEuLfxYsX\nz5gxQ/wcisNdwOAqgTCTrijhs8YEWX1tCADgliO/yUk7BivBYLC2tvbAgQOqqhrjoc7OTgBoaGjI\n7e2QtEAhZciOHTvydOW2tjYw/BzqdjHl0IWldD+Jc8kAQIRWRjPpoRW6yju45chfRDuGPJ1HPm3a\ntFAotGvXrnPPPVcffPfddwFg+vTp+bgj4hEUUmHR3t4eCoW8/BzqEZL1KaOrwDW0MqYBUVc6uOXI\nX3LSjsGFuXPnhkKhBx98UP+z7/Dhw+vWrWOMtba25ummiBdQSAVEKBRavXp19hlzk6usoRWmAV3A\nLUf+kqt2DC4sWrTohRdeCIVCCxYs+OY3v9nV1bVx48aBgYEbb7xx/Pjx+bsvkhIUUgHR2Ng4Amdf\nGu3iMQ3oHlqVUoUFbjnylxy2Y3BBluVnnnlm5cqVW7Zs2bVrFwDU1NTccccdS5Ysyd9NES+gkAoI\n3wORnKQBoWhDK9xy5C85bMeQktGjR69atSoWi3V2dtbU1DQ0NBBC8n1TJCWEc+73HJCixxRaCVdB\nUVVYKIrS3d0dDAaxyNsX9HYMeB55OYNCQvKLNbQCg7oKp8JCbDnCIm+/WLx48cyZM/E88jIHhYT4\nianCAhxcBXlOA4otR1jI4BeiHQOeR46gkJACxSkNmPMKi+7ubkVR0EZ+kacNsEgxgkJCig/vaUBI\nFVr19fWdOHECtxz5BdoIMYJCQkoNjxUW4XB45syZFRUVU6dOnTRpUkFVWJQJjz32WF43wCJFBwoJ\nKSN0V7W1tenJvQKssCgH2tvbV69ejTZCjKCQkLLDKU1UIBUW5cDIbIBFig4UElJeZLZoMWIVFuUA\n2ghxAoU0cvT09MyfP/+GG2645ppr/J5LmZKPRYscVliUA2IDLNoIsQVbB40ca9asCYfD/f39fk+k\nrMn5okVjcr8lPCXEBb0dA9oIsQUjpPyiquqBAwc+/vjjl156adOmTQBw++2333TTTX7PC/GfYulh\nkUOwHQPiDkZI+WXfvn2XXnqp37NACpGUoZXRVSXQyla0Y0AbIS5ghJRf+vr63nzzTfH4vffeW7Nm\nDUZISJYUY4UFboBFvIARUn4JBoPz5s0TjyUJv9tIDjCFVkasacBCCK3QRohH8FckgpQOBVhh8dhj\nj4XDYdwAi3gBhYQg5UJOQqu0Kiza29uxORDiHRRSbli6dOnQ0JBxZOXKlfX19X7NB0HSIucVFgCw\nYcOG9vb2W2+9NRQKYZ034gUUUm7YunXrwMCAceSuu+5CISGlgTES8pIG1O01f/78jo6ODRs2CHUt\nWLAAq+wQF1BIuWHHjh1+TwFBfMCaBgyHw21tbbfeeuuMGTNM44VclY4UAigkBEFySVtbW2NjoymQ\nAoe1KwQxgkJCip4tW7a88sore/fuDQaDU6ZMWbJkSUtLi9+TKlNEhTceRo5kBvV7AgiSFffee+8t\nt9yybdu2YDDY29u7fv36yy677I033vB7XuVIe3s74H4jJAtQSEgR09HRsXbt2ubm5ldfffWFF17Y\ntGnTo48+GovFli9fHolE/J5d2TF//ny0EZINKKSRY86cObt378a+QTnkxRdfBIClS5fq6xPz5s27\n6KKLjhw5ondsQoqOnp6e1tbW5557zu+JICMNCgkpYjo6Ohhjs2fPNg7OmTMHAN5++22fJoVkCx7U\nUrZgUQNSrMRisaNHjzY1NVVWVhrHJ02aBAAHDx70aV5IJlgPakHKEBQSUqz09vZqmlZXV2caFyM9\nPT1+TArJEDyoBQEUElK8xGIxsOuhLkbEs0ix0NDQsGrVKvFYHNTi73wQX0AhIcUKYwzsxCNGxLNI\nsYAHtSCARQ1I8VJVVQUAg4ODpnExIp5FEKSIQCEhxUowGKytrT1w4ICqqsbxzs5OAGhoaPBnWgiC\nZAoKCSlipk2bFo1Gd+3aZRx89913AWD69Ok+TQpBkAxBISFFzNy5cwHgwQcf1EcOHz68bt06xlhr\na6t/80IQJBNQSEgRs2jRojPPPDMUCi1YsODpp5++7777Fi5cODAwcMMNN4wfP97v2ZUa2EAByTdY\nzYIUMbIsP/PMMytXrtyyZYtI3NXU1Nxxxx1Llizxe2olCDZQQPINCgkpbkaPHr1q1apYLNbZ2VlT\nU9PQ0EAI8XtSpQM2UEBGEhQSUgrIsjx58mS/Z1GCYAMFZCRBISEI4gg2UEBGEhQSgiCO+NJAQRzU\nMjL38k5nZ+eSJUtisdh//Md/GGs4f/nLXz7xxBMtLS0/+9nPZFn2cYYlAFbZIQiCpGbixInXXHNN\nV1fXihUrPv30UzG4Z8+eH//4x59++umyZcvQRtmDQkIQBPHEddddN3v27O7u7nvuuQcAotHoHXfc\nMTQ0tGzZsmnTpvk9u1IAU3YIgiCeIITcf//93/72tzdv3vzb3/52586du3fvvvjii6+66iq/p1Yi\nYISEIAjilVNPPfWhhx6ilK5YsWLt2rXNzc0/+tGP/J5U6YBCQhAESYOZM2dee+21g4ODnPNly5YF\ng0G/Z1Q6oJAQBEHSIBKJvPHGG+Lx008/rWmav/MpJVBICIIgafDAAw/s3bv38ssvnz59+jvvvPPk\nk0/6PaPSAYWEIAjile3bt69bt278+PHLly+///77ZVlevXr1n//8Z7/nVSKgkBAEQTzR3d29bNky\nQsh9991XU1MzderUW265RVXVpUuXYs/ZnEA4537PAUEQpNDhnN94441/+MMfrr76arEPCQBUVb3y\nyiv/8pe/fOc733nggQf8nWEJgBESgiBIap599tk//OEPEyZMuPPOO/VBxtgDDzwQCAReeumlV155\nxcfplQYYISEIgiAFAUZICIIgSEGAQkLKiHvvvfeCCy74r//6L9M453zBggWzZ88+ePCgLxNDEARQ\nSEhZMWfOnK6urv/+7/82Zap37Nixa9eucePGNTU1+TU3BEFQSEgZMXPmzJaWlkOHDv3xj380jovD\nuS+//HKf5oUgCAAKCSkrCCFXXnklALz88sv6oKZpr776anV19SWXXOLf1BAEQSEhZcb8+fMlSfrd\n734Xi8XESCgU6u7unjt3LnbJRBB/QSEh5cVpp5329a9/vbe397XXXhMjGzduBICFCxf6OS0EQVBI\nSBnyD//wDwDw29/+FgAURfn973/f0tLypS99ye95IUi5g0JCyo6vfOUrzc3N27dvP3ny5FtvvXXi\nxIkFCxYQQvyeF4KUOygkpOwghFxxxRXRaHTz5s0bN25kjM2fP9/vSSEIgq2DkLKku7t79uzZ55xz\nzp49e774xS8+9dRTfs8IQRCMkJCyZPTo0XPmzNmxY0dfXx9uP0KQAgGFhJQpYkNSfX19a2ur33NB\nEAQAhYSULeJEtcsuu0yWZb/ngiAIAAoJKVteeOEFALjiiiv8ngiCIHFQSEh5cfz48Vgstn79+rfe\nemvWrFmTJ0/2e0YIgsTBKjukvPj+97+/efNmznl1dfX69etRSAhSOEh+TwBBRpTzzjuvp6dn7Nix\n1113HdoIQQqK/w+jOR0kqfZZsgAAAABJRU5ErkJggg==\n",
"text/plain": [
"<IPython.core.display.Image object>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"[X,Y,Z] = sphere(100);\n",
"v = [X(:) Y(:) Z(:)]';\n",
"Av = A*v;\n",
"Ray = sum( conj(v).*Av, 1 );\n",
"surf(X,Y,Z,reshape(Ray,size(X)))\n",
"hold on\n",
"plot3(V(1,:),V(2,:),V(3,:),'k*')\n",
"\n",
"title('Rayleigh quotient values on the unit sphere')\n",
"shading interp, axis square\n",
"view(107,20)\n",
"xlabel('x'), ylabel('y'), zlabel('z')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The eigenvectors (their projections onto the sphere) are shown with black markers. The max of $r(v)$ is in the middle of the yellow spot. It's the third eigenvector of $A$ as computed above, and close to $(0.4,0.5,0.75)$. "
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"evec_err =\n",
"\n",
" 0.0216\n"
]
}
],
"source": [
"x = V(:,3); \n",
"v = [0.4;0.5;0.75];\n",
"evec_err = norm(x-v)"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"eval_est =\n",
"\n",
" 5.2134\n",
"\n",
"\n",
"eval_err =\n",
"\n",
" 9.5213e-04\n"
]
}
],
"source": [
"eval_est = v'*A*v/(v'*v)\n",
"eval_err = norm(eval_est-D(3,3))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Power iteration\n",
"\n",
"Now that we can turn an eigenvector estimate into an (even better, in the hermitian case) eigenvalue estimate, we turn to finding such eigenvectors. Clearly, if $A=XDX^{-1}$ and $k$ is a positive integer, then $A^k=XD^kX^{-1}$. That is, $A$ and $A^k$ share the same eigenvectors, but the eigenvalues of $A^k$ are raised to the $k$th power. Note that if there is an eigenvalue $\\lambda_1$ such that $|\\lambda_1|>|\\lambda_j|$ for all $j>1$, the ratio $|\\lambda_1/\\lambda_j|^k\\to 0$ as $k\\to\\infty$, and in that case $A^kv$ is dominated by the eigenvector $x_1$ that belongs with $\\lambda_1$, for practically *any* vector $v$. This leads us to the **power iteration**."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAjAAAAGkCAIAAACgjIjwAAAACXBIWXMAABcSAAAXEgFnn9JSAAAA\nB3RJTUUH4AoUFDQfYlac5QAAACR0RVh0U29mdHdhcmUATUFUTEFCLCBUaGUgTWF0aFdvcmtzLCBJ\nbmMuPFjdGAAAACJ0RVh0Q3JlYXRpb24gVGltZQAyMC1PY3QtMjAxNiAxNjo1MjozMXzukQIAACAA\nSURBVHic7d1/VFR1/sfxD40oxphG0u4IAQoxwdofS3q2VssUcjxpndV0D+0xf4A/TnZqI93NXfec\n1G+aR+vs4mY/PJBCaayblQa4FOuPs7mBGraWpYY4LhhmhoIKxDTw/eP2ne80DMOvO/d+7p3n4y/4\n3GF4O+c6Lz53Pu/7Cevo6BAAAOjtOr0LAABACAIJACAJAgkAIAUCCQAgBQIJACAFAgkAIAUCCQAg\nBQIJACAFAgkAIAUCCQAgBQIJACAFAgkAIAUCCQAgBQIJACAFAgkAIAUCCQAgBQIJACAFAgkAIAUC\nCQAgBQIJACCFAXoX0AuHDx/Oy8s7ffp0VFTU+PHjFy9ePGjQIL2LAgCowzAzpH/961+PPPJIc3Pz\nvHnzxo4dm5eXt2TJko6ODr3rAgCoI8wo7+kPPPBAeHj4W2+9dd111wkh3nvvvWXLlr322mvjxo3T\nuzQAgAqMMUO6ePHiqVOnHnjgASWNhBBTpkwZOHDghx9+qG9hAAC1GCOQamtrhRAJCQmekfDw8JiY\nmHPnzulWEwBAVcYIpObmZiFEZGSk9+D1119/7do1nSoCAKjMGIGkfNDl83FXR0eHUT4AAwB0S9Nl\n342NjdOnT8/Kypo9e3bno+Xl5cXFxdXV1VarNTk5OTs7Oz4+Xjk0ePBgIURLS4v341taWmw2mwZl\nAwA0oOkMKS8v79y5c36vs61Zs+axxx7bu3ev1WptamrasWPHgw8+6FmzoASP0+n0PN7tdtfV1RFI\nAGAaQQ8kt9vtdDoPHDiQk5OzefNmv4+prKwsLCyMi4vbs2dPUVFRaWlpbm6uy+VasWJFa2urEGLE\niBE2m+2DDz7w/Mj+/ftdLteYMWOCXT8AQBtBD6SamhqHw7Fo0aLS0tKuHrNz504hxLJly2JiYpQR\nh8MxZcqU8+fPHzx4UBmZM2fOkSNH1qxZc+rUqeLi4mefffaWW26ZOHFisOsHAGgj6J8h2Wy2jRs3\nKl8fO3YsLy+v82MqKystFsuECRO8B9PT00tKSioqKtLT04UQc+fObWxs3LJlS2FhoRDi9ttvf+65\n5yIiIgL8arvdrto/AwDM4uTJk3qX4F/QA8lqtTocjh9+2QA/v87lcl24cCE2NtYnXRITE4UQdXV1\nyrcWiyUnJ+eJJ544d+7cjTfeOGTIkJ78dmlfdwOx2+28jP3Hy9h/vIaqkPkvdf1vrtrU1NTe3j50\n6FCfcWWksbHRe9BiscTFxWlXHABAK/r3IblcLuFv8qSMKEf14mxo1fG3A0BI0T+QLBaL8Bc8yohy\nVC+r3j8zcs2/95++rGMNABAi9A8kv02vnhHlqC72n7689XC9s6F14ktVE186GrKxxFV7VfAy9h+v\noenpH0hWq3XIkCG1tbVut9t7XGmD1bH1dVXZGc/X+09fmvhS1fyiL7iIBwBBon8gCSFSUlLa2tqO\nHz/uPVhVVSWESE1N1ako4bzkO2nberh+4stVq8rOEEsAoDopAmny5MlCiA0bNnhG6uvrt2/fbrFY\nJk2apFdV+x5NWzl5ZELUjxajOxtaV75/RoklvQoDAFOSIpAyMzOTkpIOHTo0Y8aM/Pz8tWvXzpw5\ns7m5OSsra8SIEXpVlRAV8YxjpBJLPoeUWBq55t/EEgCoRf8+JCFEeHh4QUHB6tWry8vLlQt3kZGR\nS5cuzc7O1ru0H2Jp7lhbweH6le//KH6UWNp6pH7eGNszDt/QAgD0SphUWwq5XC6n0xkZGWmz2cLC\nwvr5bKr3dTsbWjvHkiIhKoJYAiA/mW94IVcgqStIr3vgWHpm8sh5Y9kUA4CkZA4kKT5DMhblIt6Z\nFb/sHDzOhtb5RV/QSwsAfUAg9VFCVMSWzJQzK355b+KNPoeUXlqalgCgVwikfkmIiti35OdbMlN8\nVocLmpYAoJf4DEkdzobW/acvrXrfT/yw3gGAPGT+DIlAUpOy3mHrkXq/sbTv0bTOEykA0JLMgcQl\nOzUF7qXl/g4AEAAzpGDpanU4V/AA6IgZUijyrA73fze8l46y2AEAvBFIwaV8dNT5Ct7+05e4ggcA\n3gikoPNMlXw6ljx3aKWLFgAEgaSZrjqWlC5a2pUAgEDS1LyxNr9X8NhjCQAIJK11ewWPqRKA0EQg\n6UO5FR7tSgDgQR+SzmhXAqAlmd8YmSHpLHC7ErcMBxA6CCQpdNWu5LlluC5VAYCWCCRZsNgBQIgj\nkOQSqF2JqRIAUyOQZOS3XYk7OwAwNwJJUgEWOyh3dtCrMAAIEgJJal0tdlCmSsQSADMhkGQXeLED\n21gAMA0CyRiUxQ5sYwHAxAgkI1GmSl0tdmCqBMDQCCSDUa7gsS4cgPkQSIbEunAA5kMgGRXrwgGY\nDIFkbKwLB2AaBJLhBV4Xzv3CARgFgWQSXe34x/3CARgFG/SZTYAd//Y9mubzgROAUCPzGyMzJLNh\nXTgAgyKQzIl14QAMh0AyrQDrwucXfc5UCYBsCCST87sunKkSAAkRSObXbQst68IByIBAChUBWmhZ\n7ABABgRSCGFrJQAyI5BCDlsrAZATgRSi2FoJgGwIpNBFCy0AqRBIoY4WWgCSIJBACy0AKRBI+EHg\nFlpiCUCwEUj4fwGmSmytBCDYCCT46qqFlq2VAAQVgQQ/aKEFoD0CCV2ihRaAlggkdIMWWgDaIJDQ\nPVpoAWiAQEJP0UILIKgIJPQCLbQAgodAQq/RQgsgGAgk9EXgdeG00ALoAwIJfZcQFbElM4UWWgCq\nIJDQL7TQAlALgQQV0EILoP8IJKiGFloA/UEgQU200ALoMwIJ6qOFFkAfEEgIClpoAfRWWEdHh941\n9NSBAwfefPPN6urqm2666b777pszZ87AgQMDPN5ut588eVKz8uCXs6G14HD9yvd9EyghKmLeGNsz\nDt91EACCSuY3RsPMkN58881FixYNHjw4Ozv79ttvz83NfeKJJwyUpiGLFloAPTRA7wJ6pKWl5fnn\nn3/ooYfWrl2rjCQmJq5cufL48eOjR4/Wtzb0hNJC23mqtPVw/f7Tl5gqARBGmSHV1dVdvXp12rRp\nnpG0tDQhhNPp1K0m9BIttAACM0YgxcTEvPXWWz//+c89I59++qkQYuRI/qw2GFpoAXTFSIsaPL78\n8su5c+cmJSUVFhYGeJjMn90hwGKHfY+m+azNA6AWmd8YpfsM6cqVK8XFxZ5vo6OjMzIyPN92dHT8\n/e9/X7du3ahRo3Jzc/UoEOpQruDFR0Wsev+M98U6pYWWT5WAEBSsGVJjY+P06dOzsrJmz57d+Wh5\neXlxcXF1dbXVak1OTs7Ozo6Pj1cOnT179le/+pXnkWlpafn5+crXX3311R//+McjR44sXLhwyZIl\n4eHhgWuQ+Q8BeASYKm3JTL03cZguVQFmJfMbY7BmSHl5eefOnbt27VrnQ2vWrCksLBw0aFBqampT\nU9OOHTt27dq1adOm8ePHCyHi4+OPHj3a+ae+/PLLefPmjRgxYvfu3aNGjQpS2dCeMlWaO9Y28eUq\nn6nS/KLPmSoBoUPNRQ1ut9vpdB44cCAnJ2fz5s1+H1NZWVlYWBgXF7dnz56ioqLS0tLc3FyXy7Vi\nxYrW1i4XWbW3ty9btiwpKWnbtm2kkSmxCy0ANWdINTU13iuz/dq5c6cQYtmyZTExMcqIw+GYMmVK\nSUnJwYMH09PT/f7URx99dOLEid///veHDh3yHh89evSwYVzSMQnPVGl+0Rf7T1/yjCux5LzU+szk\nkSx2AExMzUCy2WwbN25Uvj527FheXl7nx1RWVloslgkTJngPpqenl5SUVFRUdBVI//nPf4QQ69ev\n9xnPz89XLvR1xW63e76W9rIpvNFCC6jO+51QZmoGktVqdTgcPzzvAD/P7HK5Lly4EBsbGxHxo79z\nExMThRB1dXVdPfOSJUuWLFnSh5IIISMKPFXaf/py5+0tAATg/U4oczhp2hjb1NTU3t4+dOhQn3Fl\npLGxUctiIDlaaIFQo2kguVwu4W/ypIwoRwFv7EILhA5NA8lisQh/waOMKEcBH+xCC4QITQNp8ODB\nQoiWlhafcWVEOQr4xS60gOlpGkhWq3XIkCG1tbVut9t7XLlpt81m07IYGA670ALmpvXdvlNSUtra\n2o4fP+49WFVVJYRITU3VuBgYES20gFlpHUiTJ08WQmzYsMEzUl9fv337dovFMmnSJI2LgUGxCy1g\nSloHUmZmZlJS0qFDh2bMmJGfn7927dqZM2c2NzdnZWWNGDFC42JgaEoLbed14VsP17PYATAirQMp\nPDy8oKDA4XCcOHFi/fr1BQUFLS0tS5cuzcnJ0bgSmAC70AJmotsGfS6Xy+l0RkZG2my2sLCwYPwK\nme+yDtWtKjvjdw8L7jYEeJP5jdGQO8b2kMyvO4KBXWiBbsn8xqj1JTsgeGihBQyNQILZ0EILGBSB\nBBOihRYwIgIJpkULLWAsBBLMjBZawEAIJJgfLbSAIRBICAm00ALyI5AQQtiFFpAZgYSQwy60gJwI\nJIQiWmgBCRFICF200AJSIZAQ0mihBeRBIAG00AJSIJAAIWihBSRAIAH/jxZaQEcEEvAjtNACeiGQ\nAD9ooQW0RyABXaKFFtASgQQEQgstoBkCCegeLbSABggkoEdooQWCjUACeoEWWiB4CCSgd2ihBYKE\nQAL6ghZaQHUEEtBHtNAC6iKQgH6hhRZQC4EEqIAWWqD/CCRAHbTQAv1EIAFqooUW6DMCCVAZLbRA\n3xBIQFDQQgv0FoEEBAsttECvEEhAcNFCC/QQgQQEHS20QE8QSIBGaKEFAiOQAE3RQgt0hUACtEYL\nLeAXgQTogxZawAeBBOiGFlrAG4EE6IwWWkBBIAH6o4UWEAQSIA9aaBHiCCRAIrTQIpQRSIB0aKFF\naCKQAEnRQotQQyAB8qKFFiGFQAJkRwstQgSBBBgALbQIBQQSYBi00MLcCCTASGihhYkRSIDx0EIL\nUyKQAEOihRbmQyABBkYLLcyEQAIMjxZamAOBBJgBLbQwAQIJMA9aaGFoBBJgKrTQwrgIJMCEaKGF\nERFIgDnRQgvDIZAAM6OFFgZCIAEmRwstjMKQgVRWVnbXXXedP39e70IAw6CFFvIzXiCdP3/+T3/6\nU0NDQ3t7u961AAZDCy1kZrBAam9v/93vfnfnnXfqXQhgVLTQQloGC6S8vLywsLA5c+boXQhgbLTQ\nQkJGCqTPPvvstddeW7du3XXXGalsQE600EI2hnlnb2lpWbp06dNPPz1ixAi9awHMgxZayGOA3gX4\nunLlSnFxsefb6OjojIwMIcSzzz576623Tp8+Xb/SAHNSpkpzx9rmF32x//Qlz7gSS85Lrc9MHukz\niwKCIViB1NjYOH369KysrNmzZ3c+Wl5eXlxcXF1dbbVak5OTs7Oz4+PjlUMNDQ3r16/3PDItLS0j\nI2Pv3r1lZWXbtm375ptvhBCXLl0SQnz77bfXX3/9sGHDgvRPAEKK0kJbcLh+5fs/mhVtPVy///Sl\neWNszzh8l4wD6gpWIOXl5Z07d+7atWudD61Zs6awsHDQoEGpqalNTU07duzYtWvXpk2bxo8fL4SI\nj48/evSoz4988sknV65cefDBB70HZ86cOXbs2DfeeCNI/wQg1ASeKu0/fbnz2jxARWEdHR1qPZfb\n7a6trT179uy7775bWloqhHjqqacWL17s/ZjKyso5c+bExcVt3bo1JiZGCFFWVpaTkxMdHV1WVhYR\n4f9cv3jxYkNDg+fb48ePL1++PD8/Pz4+/pZbbumqHrvdfvLkSXX+bUCIWVV2xmeqJIRIiIpgqmR0\nMr8xqjlDqqmpmTZtWuDH7Ny5UwixbNkyJY2EEA6HY8qUKSUlJQcPHkxPT/f7U8OHDx8+fLjn2ytX\nrgghRo0axQIHIEiUqZLPFTxlqrT1SP2+R9OYKkF1aq6ys9lsG//PggUL/D6msrLSYrFMmDDBe1DJ\noYqKChWLAdBPtNBCY2rOkKxWq8Ph+OF5B/h5ZpfLdeHChdjYWJ9Lc4mJiUKIurq6Hv6iO+64o4dT\nTrvd7vla2lkqILN5Y233Jt7Y1VRpS2bqvYksLJKd9zuhzDRd9t3U1NTe3j506FCfcWWksbFR9d9I\nCAH951nsMPHlKu/73SkttHyqJD/vd0KZw0nTxliXyyX8TZ6UEeUoADnRQotg0zSQLBaL8Bc8yohy\nFIC02IUWQaVpIA0ePFgI0dLS4jOujChHAUiOXWgRJJoGktVqHTJkSG1trdvt9h53Op1CCJvNpmUx\nAPqMXWgRDFrfXDUlJaWtre348ePeg1VVVUKI1NRUjYsB0B/sQgt1aR1IkydPFkJs2LDBM1JfX799\n+3aLxTJp0iSNiwHQf+xCC7VoHUiZmZlJSUmHDh2aMWNGfn7+2rVrZ86c2dzcnJWVxW0XAIOihRaq\n0Hr7ifDw8IKCgtWrV5eXlysX7iIjI5cuXZqdna1xJQDURQst+knNm6v2isvlcjqdkZGRNpstLCws\nGL9C5nsIAiamTIx8LtZxY1ZJyPzGqFsgaUDm1x0wN2dDa+etlQSxJAGZ3xgNs4U5AAOhhRZ9QCAB\nCBZaaNErBBKAIKKFFj1HIAEIOlpo0RMEEgCN0EKLwAgkANqhhRYBEEgAtDZvrC3A1kr7T1/WqzDo\ni0ACoAPPYofOU6X5RZ8zVQpNBBIA3bALLbwRSAD0RAstPAgkAPqjhRaCQAIgCVpoQSABkAgttKGM\nQAIgHVpoQxOBBEBGtNCGIAIJgLxooQ0pBBIAqdFCGzoIJAAGQAttKCCQABgDLbSmRyABMBJaaE2M\nQAJgMLTQmhWBBMCQaKE1HwIJgIHRQmsmBBIAY6OF1jQIJABmQAutCRBIAEyCFlqjI5AAmAottMZF\nIAEwG1poDYpAAmBOtNAaDoEEwLRooTUWAgmAydFCaxQEEoCQQAut/AgkAKGCFlrJEUgAQgsttNIi\nkACEHFpo5UQgAQhRtNDKJqyjo0PvGoLFbrefPHlS7yoAGICzobXgcP3K931DKCEqYt4Y2zMO3xV6\nxiXzGyMzJACgY0kKBBIA/ICOJX0RSADwI3Qs6YVAAgBfdCzpgkACAP/oWNIYgQQAXaJjSUsEEgB0\ng44lbRBIANA9Nv3TAIEEAD3Fpn9BRSABQC/QQhs8BBIA9BottMFAIAFAH9FCqy4CCQD6jhZaFRFI\nANBftNCqgkACABXQQtt/BBIAqIYW2v4gkABATbTQ9hmBBADqo4W2DwgkAAgKWmh7i0ACgCCihbbn\nCCQACDpaaHuCQAIALdBC2y0CCQC0QwttAAQSAGiKFtquhHV0dOhdQ09dvHhx27ZtFRUVHR0dU6ZM\nmTZt2vDhwwM83m63nzx5UrPyAKBXnA2tBYfrV77vm0AJURHzxtiecfiug1CFzG+MhgmkhoaGOXPm\nfP/99/fff397e/s777wTGxu7bdu2AD8i8+sOAApnQ+v8oi/2n77kMz5vrO2ZySN9ZlH9J/Mb4wC9\nC+ip3Nzc1tbWnTt3Dh06VAhxzz33PPzwwx9//PEdd9yhd2kA0HdKC23nqdLWw/X7T18K3lRJQsaY\nIbW1tY0bN27x4sULFizwDLrd7rCwsOuu6/JjMJn/EAAAH11Nle5NvLHz2rw+k/mN0RiLGmpqapqa\nmsaOHVtcXPyHP/zhscce27x5c1tbW4A0AgBjoYXWGG/oFy9eFEK8/vrr69evHzhwYERExF//+tdZ\ns2a1tLToXRoAqCmUW2il+wzpypUrxcXFnm+jo6MzMjKuXr0qhKiqqiouLr7hhhuEEJWVlXPmzCks\nLFy8eLFutQJAECjrwuOjIla9f8Y7gZQWWhN/qhSsQGpsbJw+fXpWVtbs2bM7Hy0vLy8uLq6urrZa\nrcnJydnZ2fHx8cqhhoaG9evXex6ZlpaWkZFhtVqFELNmzVLSSAjxi1/8wm63Hz58mEACYErzxtru\nTbzRZ7GDMlXaeqR+S2bqvYnDdCwvGIIVSHl5eefOnbt27VrnQ2vWrCksLBw0aFBqampTU9OOHTt2\n7dq1adOm8ePHCyHi4+OPHj3q8yMxMTFCCGV9nccNN9zg9/kBwByUqdLcsbaJL1f5TJXmF31uvqmS\nmp8hud1up9N54MCBnJyczZs3+31MZWVlYWFhXFzcnj17ioqKSktLc3NzXS7XihUrWlu7vDY6cuTI\n+Pj4yspKz0hTU9MXX3xx2223qVg/AEgodHahVTOQampqHA7HokWLSktLu3rMzp07hRDLli1TJj1C\nCIfDMWXKlPPnzx88eDDAky9evLisrCw3N/fs2bOfffbZ448/3tbWNn/+fBXrBwA5hcgutGpesrPZ\nbBs3blS+PnbsWF5eXufHVFZWWiyWCRMmeA+mp6eXlJRUVFSkp6d39eQPPfTQ1atX33jjjZdeeiks\nLOzWW2/dtm1bXFxc4JLsdrvna2mX3gNAT/S5hdb7nVBmagaS1Wp1OBw/PO8AP8/scrkuXLgQGxsb\nEfGjDq/ExEQhRF1dXeDnnzt37ty5cy9cuBAZGRkZGdmTkgghAGbi+VTJp4VWmSrtP33Zbwut9zuh\nzOGkaR9SU1NTe3u7z9oE8X+rFRobG3vyJDfffHMP0wgATMmsLbSaBpLL5RL+Jk/KiHIUANAT5muh\n1TSQLBaL8Bc8yohyFADQQybbhVbTQBo8eLAQovP9fpQR5SgAoFdMswutpoFktVqHDBlSW1vrdru9\nx51OpxDCZrNpWQwAmIY5dqHV+uaqKSkpbW1tx48f9x6sqqoSQqSmpmpcDACYidFbaLUOpMmTJwsh\nNmzY4Bmpr6/fvn27xWKZNGmSxsUAgMl020Lruv4mvWrrltZ3+87MzCwqKjp06NCMGTOmTp369ddf\nl5SUNDc3L1y4cMSIERoXAwCmFKCFNvzOHGdDq+o7o6tC60AKDw8vKChYvXp1eXm5cuEuMjJy6dKl\n2dnZGlcCACbWVQvtTV+WJETN0rGwAHTbwtzlcjmdzsjISJvNFhYWFoxfIfNOvQCgmVVlZ5Sp0r2J\nN36VmyntG6NugaQBAgkAFM6G1oLD9ROSblx8/y+kfWOUbsdYAIDqlCt4elfRDa1X2QEA4BeBBACQ\nAoEEAJACgQQAkAKBBACQAoEEAJACgQQAkAKBBACQAoEEAJACgQQAkAKBBACQAoEEAJACgQQAkAKB\nBACQAoEEAJACgQQAkAKBBACQAoEEAJACgQQAkAKBBACQAoEEAJACgQQAkAKBBACQAoEEAJACgQQA\nkAKBBACQAoEEAJACgQQAkAKBBACQAoEEAJACgQQAkAKBBACQAoEEAJACgQQAkAKBBACQAoEEAJAC\ngQQAkAKBBACQAoEEAJACgQQAkAKBBACQAoEEAJACgQQAkAKBBACQAoEEAJACgQQAkAKBBACQAoEE\nAJACgQQAkAKBBACQAoEEAJACgQQAkAKBBACQAoEEAJACgQQAkAKBBACQAoEEAJACgQQAkAKBBACQ\nAoEEAJACgQQAkMIAvQvohQ8//HDLli1nzpyJjo6+6667Fi9ePHjwYL2LAgCowzAzpD179mRnZw8c\nOHDhwoWpqamvv/76woUL29vb9a4LAKAOw8yQCgoKRo8e/dJLL4WFhQkhoqOjc3NzT5w4kZqaqndp\nAAAVGCaQwsPDIyMjlTQSQkRFRQkhrr/+el2LAgCoxjCX7DIzMysqKt55552WlpZDhw4VFBTcdddd\nCQkJetcFAFCHYWZIU6dO/fzzz5cvX758+XIhxM033/y3v/1N76IAAKqRLpCuXLlSXFzs+TY6Ojoj\nI0MI8eKLL27ZsmXWrFnjxo27cOFCfn7+r3/967y8vNjYWP2KBQCoJliB1NjYOH369KysrNmzZ3c+\nWl5eXlxcXF1dbbVak5OTs7Oz4+PjlUMNDQ3r16/3PDItLS0jI6O5ufnVV1+dO3fu008/rYxPnDhx\n6tSpu3fvXrJkSZD+CVDY7faTJ0/qXYXh8TL2H6+h6QUrkPLy8s6dO3ft2rXOh9asWVNYWDho0KDU\n1NSmpqYdO3bs2rVr06ZN48ePF0LEx8cfPXrU50c+/fTTtra2cePGeUbi4uISExM/+eSTINUPANCY\nmosa3G630+k8cOBATk7O5s2b/T6msrKysLAwLi5uz549RUVFpaWlubm5LpdrxYoVra2tXT1zTEyM\nEKK8vNwzcu3atbq6ultuuUXF+gEAOlJzhlRTUzNt2rTAj9m5c6cQYtmyZUrGCCEcDseUKVNKSkoO\nHjyYnp7u96diY2PHjx+/a9eum266yeFwtLa2vvDCC999992sWbNUrB8AoCM1A8lms23cuFH5+tix\nY3l5eZ0fU1lZabFYJkyY4D2Ynp5eUlJSUVHRVSAJIV544YX169e//PLLL774ohAiMTHxlVdeue22\n2wKXZLfbe/3PQCe8jKrgZew/XkNzUzOQrFarw+H44XkH+Hlml8t14cKF2NjYiIgI7/HExEQhRF1d\nXYAnHzZs2Nq1a1etWvXVV18NHTp02LBh3dbD558AYCCaLvtuampqb28fOnSoz7gy0tjY2O0zhIeH\ne9bjAQDMRNM7NbhcLuFv8qSMKEcBAKFJ00CyWCzCX/AoI8pRAEBo0jSQlO2LWlpafMaVETY3AoBQ\npmkgWa3WIUOG1NbWut1u73Gn0ymEsNlsWhYDAJCK1nf7TklJaWtrO378uPdgVVWVEIKdjQAglGkd\nSJMnTxZCbNiwwTNSX1+/fft2i8UyadIkjYsBAMhD67t9Z2ZmFhUVHTp0aMaMGVOnTv36669LSkqa\nm5sXLlw4YsQIjYsBAMhD60AKDw8vKChYvXp1eXm5cuEuMjJy6dKl2dnZGlcCAJBKWEdHhy6/2OVy\nOZ3OyMhIm83m2ZgcABCydAskAAC8SbdjbP8F2P0PPVRaWrpv3z6fwZtuuknZPx4B9HlrSngL8DJy\ncvbEu+++e+DAgVOnTt1www2jR49esGDBT37yE5/HSHg2mi2QAu/+hx4qKSnZuNIOTAAABG5JREFU\nu3dvZGSk9yCrTnqiz1tTwluAl5GTMzC32/3444//85//tFqtP/vZz86fP19YWLhjx478/PwxY8Z4\nHibp2dhhIhUVFcnJyRkZGXV1dcrIP/7xj5SUlHvuuaelpUXf2owlPT39N7/5jd5VGMb3339/5syZ\n/fv3P/nkk8nJycnJya+88orPYzg5u9WTl7GDk7M727ZtU6Y7nvNqx44dycnJd999d1tbmzIi7dmo\ndR9SUHW1+9/58+cPHjyoa2lG0tzcXFdXl5KSonchhlFTU+NwOBYtWlRaWtrVYzg5u9WTl5GTs1tb\nt24dMGDA2rVrPbv8zJo16+677/76669PnTqljEh7NpoqkLra/U8IUVFRoVNRxnPq1KmOjo5uNz+E\nh7I1pWLBggV+H8PJ2a2evIycnIF1dHR89dVXdrv95ptv9h4fOXKk8NpzTtqz0TyfIfVn9z94UzY2\n/OlPf7phw4bPPvts8ODBdrv9kUceGT58uN6lSSqoW1OGjm5fRsHJ2R23271q1arO6xeqq6uFEAkJ\nCULus9E8gdT/3f+gUP7PP/nkk999911CQkJtbe2+ffu2b9+em5v7y1/+Uu/qDImTUy2cnIENGDDg\noYce8hn897///dFHHyUlJSUlJQm5z0bzXLJj9z+1nDx5MiwsbMGCBUeOHHnvvfc+/vjjJ598sqmp\nafny5VevXtW7OkPi5FQLJ2dvffDBB0uWLBk0aND//M//eO9IJ+fZaJ4ZErv/qeUvf/mLy+XyrKO1\nWCyPPvpoTU3N7t27y8rKOv/9hW5xcqqFk7Pnmpqannvuubfffnv48OG5ublpaWnKuMxno3lmSOz+\np5bo6OjOXR3Kbdq//PJLPSoyPE5OtXBy9tDevXvvv//+t99+e+rUqcXFxd4dSDKfjeaZIXnv/ucd\n8uz+11vt7e1hYWE+Nxi0Wq1CiG+//VanooyNk1MtnJw9sWXLlnXr1sXExBQUFNx5550+R2U+G80z\nQxLs/qcGp9OZkpLy1FNP+YwrHyYrH4qiDzg5+4+Tsyf27t27bt26MWPG7N69u3MaKaQ9G00VSOz+\n13/x8fFRUVF79+49ceKEZ/Dy5cv5+fkDBgy47777dKzN0Dg5+4+Ts1tut3vdunVDhgx59dVXlYmj\nX9Kejea5ZCfY/U8NYWFhK1eufOKJJzIzMx9++GG73V5fX79t27Zvvvnmt7/97ahRo/Qu0Kg4OfuP\nk7Nb1dXVZ8+ejYuL+/Of/9z56Pz582NjY4XEZ6OpAond/1ThcDhefvnl559//rXXXhNChIWFxcXF\nvfjii/wF2h+cnKrg5Azs2LFjQoj//ve/b7zxRuej06ZNUwJJ2rPRnPshsfufKi5fvlxfXx8XF+dz\nZ2X0ByenKjg5VSHb2WjOQAIAGI6pFjUAAIyLQAIASIFAAgBIgUACAEiBQAIASIFAAgBIgUACAEiB\nQAIASIFAAgBIgUACAEiBQAIASIFAAgBIgUACAEiBQAIASIFAAgBIgUACAEjhfwGglvAROb9gFQAA\nAABJRU5ErkJggg==\n",
"text/plain": [
"<IPython.core.display.Image object>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"v = rand(3,1);\n",
"for k = 1:20\n",
" v = v/norm(v);\n",
" evec_err(k) = min(norm(v-x),norm(v+x));\n",
" v = A*v;\n",
"end\n",
"semilogy(evec_err) \n",
" "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Since $A$ is symmetric, we should expect around 12 accurate digits in an eigenvalue estimate derived from a Rayleigh quotient with our final $v$."
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"ans =\n",
"\n",
" 1.8918e-13\n"
]
}
],
"source": [
"abs( v'*A*v/(v'*v) - D(3,3) )"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now how much would you pay? But wait--there's more!"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Inverse iteration\n",
"\n",
"If $\\lambda_J$ is an eigenvalue of $A$, then $(\\lambda_J-\\mu)^{-1}$ is an eigenvalue of $(A-\\mu I)^{-1}$, with the same eigenvector. If $\\mu$ is closer to $\\lambda_J$ than to any other $\\lambda_j$, a power iteration on $(A-\\mu I)^{-1}$ should converge very quickly. This amounts to repeatedly solving a linear system with the matrix $(A-\\mu I)$."
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAjAAAAGkCAIAAACgjIjwAAAACXBIWXMAABcSAAAXEgFnn9JSAAAA\nB3RJTUUH4AoUFDQg1DCx2AAAACR0RVh0U29mdHdhcmUATUFUTEFCLCBUaGUgTWF0aFdvcmtzLCBJ\nbmMuPFjdGAAAACJ0RVh0Q3JlYXRpb24gVGltZQAyMC1PY3QtMjAxNiAxNjo1MjozMuXnwLgAACAA\nSURBVHic7d1/cFRVnvfxk+40IUkHBIQ1v1GchIQfO6KUDqIIYUhqEWtA4EELUIlACSuCxK0peWpn\nYIbRgXHXuOOMUlGHuGI2kqlFk7A7RAyjQBIcWHBTD4khttOJgfBjSCDdwaaT54+LTdvd6XSSvn3P\n7X6//rCS053uM7funA/33nO+J6q3t1cAAKA1g9YdAABACAIJACAJAgkAIAUCCQAgBQIJACAFAgkA\nIAUCCQAgBQIJACAFAgkAIAUCCQAgBQIJACAFAgkAIAUCCQAgBQIJACAFAgkAIAUCCQAgBQIJACAF\nAgkAIAUCCQAgBQIJACCFaK07MADHjh0rKio6c+bM6NGjZ86cuXbt2piYGK07BQAIDt1cIX366acr\nVqyw2WxPPvnk9OnTi4qK1q1b19vbq3W/AADBEaWXMX3BggUmk2nv3r0Gg0EI8dFHHxUUFLz99tv3\n33+/1l0DAASBPq6QLly40NjYuGDBAiWNhBB5eXnDhg377LPPtO0YACBY9BFIVqtVCDF+/HhXi8lk\nSk5Obm1t1axPAICg0kcg2Ww2IUR8fLx7Y1xcXFdXl0Y9AgAEmT4CSXnQ5fG4q7e3Vy8PwAAA/Qrp\ntO+Ojo6FCxeuWrVq+fLl3q9WVVWVl5c3NTWZzeaMjIz8/Pz09HTlpdjYWCGE3W53f7/dbk9MTAxB\ntwEAIRDSK6SioqLW1laf99m2b9++fv36gwcPms3mzs7O0tLSRx55xDVnQQkei8Xier/T6WxpaSGQ\nACBsqB5ITqfTYrEcOnRo06ZNu3bt8vme2tra4uLitLS0/fv3l5SUVFZWFhYWOhyOLVu2dHd3CyGS\nkpISExMPHDjg+pPq6mqHw3HPPfeo3X8AQGiofsuuubn54Ycf9v+esrIyIURBQUFycrLSkpubm5eX\nV1FRcfjw4ZycHCHEypUrf/3rX2/fvn3JkiWNjY2vvPJKamrq7NmzVe189ZnLQojxo4aPHz1c1S8C\nAKgeSImJia+99pry86lTp4qKirzfU1tbazQaZ82a5d6Yk5NTUVFRU1OjBNITTzzR0dHxzjvvFBcX\nCyGmTJny0ksvDR/uLycyMzOH0nNH3Jiv5vzS9avJdlEIYbJfjL7xwyWT/aLSGP3dDwAgv4aGBq27\n4JvqgWQ2m3Nzc298WbSPr3M4HO3t7SkpKR7pMmHCBCFES0uL8qvRaNy0adOGDRtaW1tHjRqVkJAQ\nyLcP5bhXn7k8+3fHb/YzbsyN/47x/X7lKmr8qNjvfhiePnr4+NGxQufXWJmZmdKevjrCYRw6jmFQ\nDPFf6qrSvrhqZ2dnT0/PyJEjPdqVlo6ODvdGo9GYlpYWmo5ZLtn7f9P33t99479nfL/BlUnhF1oA\nMHTaB5LD4RC+Lp6UFuVVTYwfHTt+9PAeW8tfu28NygcqiSUGFVpKfwS5BSB8aR9IRqNR+AoepUV5\nVRMPTbjl1EP/5rS3XMp6xxCbYvlbt3LNdOjMZeG6Hvqb3RUzQRFIaAm324NuP9/MrfGjhgu3bAMA\nXdA+kHwuenW1KK9q4srhxxwXa4QQo//fUyNmvD9+QoqYcIsQ4snpnoufvgunG4n19aVuy9+6VUqs\n732j8skB5JYrnGZNuEUM8GKLu/ZBwWEcOo5h2NM+kMxmc0JCgtVqdTqd7tdDyjJYrZa+2hteVdJI\nCNFja+k88tgtcz/t6803xv3Rw5XE8uZKLPHdoylVL7M8v9ftYusPx9q83xaU3AKAIdI+kIQQWVlZ\ndXV19fX1U6dOdTUeP35cCJGdna1Jl2JSF1+/UOueSV0nXoi/a+fgPu1mYgnR72WWCG1oCXILgByk\nCKR58+bV1dXt3Lnz3XffVVra2tr27NljNBrnzJmjSZcMcSnxd+3sPPJYj+3GvPNr1r1CiEFnUr8C\nCS3Rx+1BIXFuMSkDQICkCKRly5aVlJTU1dUtWrRo/vz5586dq6iosNlsq1evTkpK0qpXhriUETPe\n98gkQ1xybOZGrbokArg9qJAqt7wxmRCANykCyWQy7d69e9u2bVVVVfX19UKI+Pj4zZs35+fna9sx\n59XeETPev1z1gKvF3lAohNA2kwIheW4FazIhoQWEkyipthRyOBwWiyU+Pj4xMTEqKmqInzb0dd0t\nP1sshIiZYDDecszVaIhLMf9wZ/St9w2xezqi4fVWILhJCARO5oIXcgVScA39uDcuvnHDMO7vo+P+\n/uYMwAjMpH7dvOiRMrf83yRk5RYih8yBJMUtOzmd/e3N+3K2k9eFEK5M6rG1XP2fF0bMeN8Ql6JN\n5+Rzc8T3e5/QezKhe25Z/mYXbtkWRAO9ScjFFhB6XCH1yXV5pDCao2ImGD2uk/wsTsJQeNwkFCGc\nBN8vnmxB12S+QiKQ+uRot9rrj3ZWl9rqjygtRnOUeUa06babuxoa4lLip+4xjUsdal8xKDI/3OJi\nC3IikLQRrOPuaLdeLH2ls7pUCGE0R42cZzKYb064cHYmOS9PHzn7/8RO+tHQvwtBp+GK4355X2yx\n1hhqI5C0Edzj7mi3dlaXdn5S2mNv8cgk20mn7eR109jUEbOXjlm6OVjfiFDiYgsRgkDShhrHXbmP\nd/Xz0uGZf3FvVzJJCKHE0oiHlnIfL/zo4mKL0IJ/BJI2VD3u1y/UdB55zPVrz9Ve20ln9xmn8qtp\nbGrspB9xHy8CyTz9va/pGIRWRCGQtKH2cb9m3dt14gXXrz1Xe68cvu441+P+Hu7jwSfuEEIrBJI2\nQnDc7Q2vKsWEFD1Xezv+5HBe9TykXDBhEPpds0VoYRAIJG2E5rgHmEkKnjAh6KS92OIOoZwIJG2E\n5rgrWyW5dk4SQvRc7R2WuMN9AZMHLpgQYrpYaMzFVmgQSNoI2XH3zqSY1MXxd+10X8DkExdMkIf/\n6RiEVtggkLQRyuOubHPu2jlJCBGb+ZyyS4V3xQcPygVT3KQZI2YvDU1vgUHjDqHeEUjaCPFx77G1\nuO+cJNwyScEFEyIEC7ZkRiBpI/THvd9MElwwAUII7hBqh0DShibH3WNxkp+dk5RaRBdLX+nro7hg\nAnR0h1AvdQgJJG1oddw9JoIb4lL87JzU7wWTEMI0NnXM0s1cMAE+6XQOoVbbQhJI2tDwuA8okxSB\nXDAxWRwYHFYZuxBI2tDwuPfYWuwNhdese10tAe7mF+AFE7fyADVEwh1CAkkb2h73vhYnBfjn/V4w\nCSHiJs0Y8dBSbuUBoSR5aIm+7xAqjQSSNjQ/7n4WJwVIuWCy1R/xP1mcW3mAVPq6Q6h5aAkhbj/4\nf5tP9HkDRlsEkroCmQgeCOWCSQmnvt7DrTxAXzS52Moof0bzgbEvBJLqPDLJEJcSm/lcTOriwX1a\ngHMfWMYEhIeghxaBpA1JAkkMZHFSgAKc+8CtPCASBH6HcPzo4cOKn5JkYPRGIIXIICaCB4JbeQAC\nZLnUPX70cKkGRg8EUuh4Z1IgE8EDFEihPNO4VGblARFOtoHRHYEUOt4TwU1j7ku4//0gfgW38gD4\nJ9vA6I5ACinvieADWpwUOOVWXucnpY7z1r7eo9zKi5s0g2QCIoeEA6MLgRRqQ1+cNCD93soT3yXT\nmKWbVeoDAHnIOTAqCCQNBGtxUuACv5XHfHEgvEk7MAoCSStBnwgeIB4yARFO5oGRQNKMShPBAxT4\nQybmiwPhROaBkUDSkraZJIRwtFuvn2/p+OQ/+n3IxIZMQHiQeWAkkLQ06F0qgo6HTECEkHlgJJA0\nNsRdKoIuwNIPPGQCdErmgZFA0l6IJ4IHKJANmXjIBOiOzAMjgSSF0E8ED1AgGzIJHjIB+iHzwEgg\nySK4u1QEHQ+ZgPAg88BIIElEq8VJA8JDJkDXZB4YCSS5aD4RPHCBP2SiXB4gD5kHRgJJOqruUhF0\ngT9kolweIAOZB0YCSTreE8ElzyQFD5kAXZB5YCSQZBSyXSrUwEMmQGYyD4wEkqTkXJw0IDxkAiQk\n88BIIMlL2sVJAzKgh0yssQXUJvPASCBJ7fqFms4jj7l+lW1x0oAE+JDJNC51xENLecgEqETmgZFA\nkp0uFicNCA+ZAA3JPDASSDqgo8VJAxLgnkwkExBEMg+MBJI+hGsmKez1RwPZk4mHTMDQyTwwEkj6\nIM/OSeqhkCsQAjIPjASSbsi2c5J6WGMLqEfmgZFA0hNdL5gdBB4yAUEn88BIIOlMGCyYHQTW2ALB\nIvPASCDpT3gsmB0E1tgCQyfzwEgg6VL4LU4aENbYAoMm88BIIOmV90TwiMokBWtsgYGSeWAkkHQs\nvBcnDQjTH4AAyTwwEkj6pq/d/NTmaLdeP9/CGlvAD5kHRgJJ33S6m5/a2McW6IvMAyOBpHuRtjhp\nQFhjC3iQeWAkkMJBZC5OGhCmPwAKmQdGAilMROzipIFijS0inMwDo24C6fr16+vXr+/p6XG1GAyG\nN99808+fyHzc1eC9m19M6qNkkk+ssUXEknlg1E0gNTQ0PPLII48//vjw4cOVFoPB8MILL/j5E5mP\nu0oifMHsILDGFpFG5oFRN4G0b9++X/ziF8eOHYuKigrwT2Q+7uphcdLgBP6QiekP0DWZB8ZorTsQ\nqNOnT2dlZQWeRhFLuUfnyiRlvgOZ1C/TuFRl/refNbaO81ZH9Y1HUEx/AIJON1dITz31VEJCgsFg\nOHHixIgRI2bNmrV+/frY2Fg/fyLzPwRUFQm7+YVA4PvYMv0BOiLzwKibQPrRj37U1dW1bNmyzMzM\nL774orS0dOrUqXv27DEYDH39iczHXW0smA0W1tgizMg8MEoXSFeuXCkvL3f9Onbs2Llz5zqdzvfe\ne2/atGmTJ09W2svKyl588cXCwsK8vLy+Pkrm4x4CLJgNLtbYIjzIPDCqFUgdHR0LFy5ctWrV8uXL\nvV+tqqoqLy9vamoym80ZGRn5+fnp6enKS19//fVPfvIT1zunTZv21ltveX+C0+m86667VqxY4Wei\nnczHPTRYMKsG1thC12QeGNWa1FBUVNTa2trV1eX90vbt24uLi2NiYrKzszs7O0tLS/ft2/f666/P\nnDlTCJGenn7ixAmPP7lw4UJzc/Pdd99tNBqVFqPRGB0d7fPz4aJMsXNfMKtMdiCThsJj+oPPNbau\n6Q88ZAIC1+cDmEFwOp0Wi+XQoUObNm3atWuXz/fU1tYWFxenpaXt37+/pKSksrKysLDQ4XBs2bKl\nu7u7r0+2WCwrVqw4evSoq+XkyZNdXV1ZWVlB7H9Y8n50ZG8otDe8qlV/womSTLf/rjZ1a9mIh3zf\no3Oct14sfcX6s0e/eubei6WvONr73B0DQDADqbm5OTc3d82aNZWVlX29p6ysTAhRUFCQnJystOTm\n5ubl5Z09e/bw4cN9/dXUqVMzMjK2bNly4MCBtra2P//5zwUFBenp6QsXLgxi/8OVIS7F49HRNWvZ\n9Qs1fb0fA2Ialxo76Ue3/eOrt/+u9rb1r8ZNmuHzbUoyfbXu3q+eubfzE3/zI4CIFcxnSFevXnWF\nyqlTp4qKip5//vm1a9e6v2fWrFnnz58/fvy4q+CCEKKiouL5559fuXLlli1b+vrwb775Ztu2bdXV\n1b29vSaT6cEHH9y6devYsWP99CczM9P9V2lvm4YGO8yGDNMfIBu9DIbBfIZkNptzc3NvfG60j092\nOBzt7e0pKSnuaSSEmDBhghCipaXF+09ckpKS3njjDbvd3t7enpyc7PPzvUl73EPPe8Hs1f95gQWz\najCNSzWNSx0xeylrbCEJ95HQI5ykEsxbdv3q7Ozs6ekZOXKkR7vS0tHR0e8nxMbGpqenB5hG8BCb\nuTE28znXr95z8BBcNx4y/b729t/V+nnI1Fld6nrIZK8/6vNtQCQIaSA5HA7h6+JJaVFehapiUheb\nxty8Tadkkob9iRCmcamuh0xMfwD6EtJAUiZteweP0uKa0g31KBMcPDLJYyMlqES5jxfg9IeWny1u\n+dlipj8gooT03pdSes5ut3u0Ky3+C9MhWJRMcr9Zp9QZoohDyHg8ZPK5xtZx3uo4b7XVH+EhEyJH\nSAPJbDYnJCRYrVan0+l+PWSxWIQQiYmJoexMJFMWzLpn0jXrXkNcMgtmQ2xAJcap/oCwF9JbdkKI\nrKysb7/9tr6+3r3x+PHjQojs7OwQdyaSKZnk3sKCWQ0pyZSyda//Nbbu0x94yITwE+pAmjdvnhBi\n586bd4fa2tr27NljNBrnzJkT4s5EOIo4yMZjja3/6Q/KGlufhYsAnQr1/Olly5aVlJTU1dUtWrRo\n/vz5586dq6iosNlsq1evTkpKCnFnoDxPct/1/Jq1zDTmPhbMasv1kGnM0s1+1tgqydT5SSlrbBEe\nQh1IJpNp9+7d27Ztq6qqUm7cxcfHb968OT8/P8Q9gSImdbGyoZ/yq7JgliIOkgh0+gNrbBEWNNsP\nyeFwWCyW+Pj4xMRElTYml7nKumy8CwtRxEFOfkqMu1BiHH7IPDBKt0FfEMl83CVEJunIgPaxHfHQ\nUtO41JD1DZKTeWAkkHADu57rUSCFXIUQcZNmjHhoKQ+ZIOQeGAkk3EQm6Vfg+9gy/SHCyTwwEkj4\nHu+KqzGpiynioCN+1ti6sMY2ksk8MBJI8OSdSbGZz1HEQXcc7daLpa8E8pCJ6Q8RReaBkUCCD94V\nV8kknRrQ9AeljhHCm8wDI4EE365fqPHYmYJM0rUA97E1jUtl+kN4k3lgJJDQp2vWve5FHAxxKTGp\nj5JJehf49AceMoUlmQdGAgn+eC9OoohD2GD6Q2SSeWAkkNAPMim8Odqt18+3dHzyH6yxjRAyD4wE\nEvpHEYdIEOAaW9PY1DFLN/OQSb9kHhgJJPSvx9ZyzbqXTIoQAU5/YI2tTsk8MBJICAhFHCIQD5nC\nkswDI4GEQJFJEYsS4+FE5oGRQMIAUFgoklFiPDzIPDASSBgYCguBNba6JvPASCBhwCgsBAUlxvVI\n5oGRQMJgkElwx/QHHZF5YCSQMEgUFoI3e/3RANfYMv1BKzIPjAQSBs87kyjiAEGJcbnJPDASSBgS\nCgvBD9bYSkjmgZFAwlBRWAj9osS4PGQeGAkkDBWFhRA41thqTuaBkUBCEFDEAQNCiXENyTwwEkgI\nDjIJg0CJ8dCTeWAkkBA0FBbCoDH9IWRkHhgJJAQTmYQhYo2t2mQeGAkkBBlFHBAUjnbrxdJXWGMb\ndDIPjAQSgo9MQrBQYjzoZB4YCSSowqOIgyCTMDSUGA8WmQdGAglqodgd1MAa2yGSeWAkkKAiCgtB\nPUx/GByZB0YCCeoik6Aq1tgOlMwDI4EE1VHsDiFAifEAyTwwEkhQHcXuEEqssfVP5oGRQEIoUFgI\nocf0B59kHhgJJIQImQStUGLcncwDI4GE0KGwEDTEGluFzAMjgYSQIpOguQBLjMdNmhGWa2xlHhgJ\nJISadyZRxAGaCPwhUzhNf5B5YCSQoAGK3UEqEbXGVuaBkUCCNq5fqOk88ph7C5kEzUVCiXGZB0YC\nCZqh2B3kFN5rbGUeGAkkaInCQpBZWJYYl3lgJJCgMTIJ8gunNbYyD4wEErRHsTvoRRhMf5B5YCSQ\nIAUyCTqi6xLjMg+MBBKk0GNrsTcUXrPudbWQSZBfgGtsTWNTxyzdLMlDJpkHRgIJsqDYHfRLRyXG\nZR4YCSRIxLuIA5kEfZH/IZPMAyOBBLlQ7A7hQdoS4zIPjAQSpEMmIWxIWGJc5oGRQIKMKHaHMCPP\nGluZB0YCCZIikxCWNC8xLvPASCBBXh7F7gSZhDCi1fQHmQdGAglSowArwp69/miAa2yDMv1B5oGR\nQILsKHaHSBCyEuMyD4wEEnSATELkUHuNrcwDI4EEfaDYHSKNSiXGZR4YCSToBpmEyBTcNbYyD4wE\nEnSDAqyIZMEqMS7zwEggQU8owAoMscS4zAOjpIF06tSptWvXfvzxx3Fxca7GY8eOFRUVnTlzZvTo\n0TNnzly7dm1MTIyfD5H5uGPQKMAKKAY3/UHmgVHGQDp37tyzzz578uTJ48ePx8fHK42ffvrp6tWr\np0+fnpub29bW9u67706fPr2oqCgqKqqvz5H5uGMoKHYHuBvQGtsfLnpS2oFRrkD64IMPiouLm5qa\nenp6hBDugbRgwQKTybR3716DwSCE+OijjwoKCt5+++3777+/r08jkMIYmQR4c7RbL5a+4v8h08d/\nMz3z8dch69KAGLTuwPdMmTJlzZo1O3bsWLx4sXv7hQsXGhsbFyxYoKSRECIvL2/YsGGfffaZFt2E\n9pTpDO4t16x77Q2vatUfQAamcam3/eOrt/+u9rb1r454yPcSpY//ZgpxrwIXrXUHvmfixIkTJ04U\nQnR3d+/de3MyldVqFUKMHz/e1WIymZKTk1tbW0PeR8hCeXTkXoBVmRROYSFEONO4VNO4GwUdvB8y\nfdFl1LBv/skVSH2x2WxCCNftO0VcXFxXV5dGPYIUlOukziOPuVrIJMDFlUyuNbbRY1PEF/u17lef\ntAmkK1eulJeXu34dO3bs3Llz/bxfedDl8birt7dXqgdg0ET0rffF37XTvQArmQR4MI1LvVn+7t/C\nNJA6OjoWLly4atWq5cuXe79aVVVVXl7e1NRkNpszMjLy8/PT09OVly5durRjxw7XO6dNm+Y/kGJj\nY4UQdrvdvdFutycmJg6l/wgPMamLlTWzrpZr1jJDXEpM6mI/fwVANkMKpKKiotbWVp/3zbZv315c\nXBwTE5Odnd3Z2VlaWrpv377XX3995syZQoj09PQTJ04E/kVK8FgsFleL0+lsaWnxM8UOEUW5HnJl\nkpJPxtgUCrACOjLgWXZOp9NisRw6dGjTpk27du3y+Z7a2tri4uK0tLT9+/eXlJRUVlYWFhY6HI4t\nW7Z0d3cPopdJSUmJiYkHDhxwtVRXVzscjnvuuWcQn4awFJu5MTbzOdevPbaWq//zgvu8cACSG3Ag\nNTc35+bmrlmzprKysq/3lJWVCSEKCgqSk5OVltzc3Ly8vLNnzx4+fHhwHV25cuXnn3++ffv2xsbG\n8vLyX/7yl6mpqbNnzx7cpyEseWeSx1olADIb8C27xMTE1157Tfn51KlTRUVF3u+pra01Go2zZs1y\nb8zJyamoqKipqcnJyRlER5944omOjo533nmnuLhYCDFlypSXXnpp+PDhg/gohDHluZH7vbvOI49R\ngBXQhQEHktlszs3NvfHH0T7+3OFwtLe3p6SkeKTFhAkThBAtLQH9c3XJkiVLlixxbzEajZs2bdqw\nYUNra+uoUaMSEhIC+ZzMzEzXz1RtiATKXIbrF2pdBViVTKLYHSKZ+0gos+BP++7s7Ozp6Rk5cqRH\nu9LS0dExlA83Go1paWmBv58QikCGuBRlIrh7Jl2ueoBMQsRyHwllDqfglw5yOBzC18WT0qK8CqhK\nyST323Q9tpYrhx/z8ycANBf8QDIajcJX8CgtyquA2rz37nNcrHFfPwtANsEPJJ+LWF0tyqtACHhn\n0jXrXjIJkFbwA8lsNickJFitVqfT6d6uLGultgJCiaLggI6osv1EVlbWt99+W19f7954/PhxIUR2\ndrYa3wj0xXs/WXtDIZkESEiVQJo3b54QYufOm7ultbW17dmzx2g0zpkzR41vBPxQ5ji4t5BJgIRU\nqfa9bNmykpKSurq6RYsWzZ8//9y5cxUVFTabbfXq1UlJSWp8I+CfsmDW/QHSNWuZoCg4IBNVAslk\nMu3evXvbtm1VVVXKjbv4+PjNmzfn5+er8XVAIDyKgvfYWq5Zy0xj7qMAKyCJKFW3FHI4HBaLJT4+\nPjExMSoqSr0v8ikzM5OFsfBgb3jVfaMKQ1yK+Yc7ySREDpkHRlWeIbmYTKYf/OAHSUlJoU8jwCeK\nggPSUjeQAAnFpC6mKDggIQIJEUcpwOq+nyyZBMiAQEIkMsSlxGY+Zxpz89GRkkkadgkAgYQIpSxO\n8siky1UPaNglIMIRSIhcFAUHpEIgIaJRFByQB4GESOezKDiFhYDQI5AAH0XBKXYHhB6BBAhBUXBA\nAgQScANFwQFtEUjATTGpiz0y6Zq1jEwCQoNAAr7Hu7DQNWvZ9Qs1GnYJiBAEEuDJZwFWMglQG4EE\n+EBRcCD0CCTAN4qCAyFGIAG+KUXBvQuwkkmASggkoE8+C7BSFBxQCYEE+OOzACtFwQE1EEhAP7yL\n3VEUHFADgQT0j6LgQAgQSEBAKAoOqI1AAgJFUXBAVQQSMAAUBQfUQyABA0NRcEAlBBIwYBQFB9RA\nIAGDQVFwIOgIJGCQKAoOBBeBBAweRcGBICKQgCHxziQKsAKDQyABQ8VGFUBQEEjAULFRBRAUBBIQ\nBGxUAQwdgQQEBxtVAENEIAFB43OjCoqCAwEikIBg8lkUnEwCAkEgAUHmXRScjSqAQBBIQPAZ4lJG\nP/KVe707e0PhpQ9vJ5YAPwgkQC0+a7Bes+7Vqj+A5AgkQEXea2btDYXUuwN8IpAAdVHvDggQgQSo\njnp3QCAIJCAUfNYW0rA/gIQIJCAUfNYWoo4D4I5AAkLEZ20h1swCLgQSEDo+6ziwOAlQEEhASHnX\ncbA3FJJJgCCQgNAzxKXcMvdT9xYyCRAEEqAJ5XmSewtFHAACCdAGRRwADwQSoBmKOADuCCRASxRx\nAFwIJEBjMamLY1IXu36liAMiFoEEaMwQlxKb+RxFHAACCdAeRRwAQSABkqCIA0AgAbKgiAMiHIEE\nSIQiDohkBBIgF4o4IGIRSIB0KOKAyEQgATLyWcSBe3cIbwQSICnvTLI3FF6ueoA6DghX0Vp3wLdT\np06tXbv2448/jouLU1quX7++fv36np4e13sMBsObb76pUQeBUIhJXdxja3V/gKTUcYhJfTQ2c6OG\nHQPUIGMgnTt37pe//OWlS5d6e3tdjWfOnKmurn788ceHDx+utBgMXN4hzClFHAxxyfaGQlfjd4+U\naj3W0gJ6J1cgffDBB8XFxU1NTe5XQorTp08nJCT88z//c1RUlCZ9AzRhiEuJAGBlcwAAD+hJREFU\nzdwYk7rYo+iq42INl0oIM3IF0pQpU9asWSOEqKmp2bv3e/NcT58+nZWVRRohMilrZq9Z93pfKgkh\nyCSEB7kCaeLEiRMnThRCdHd3ewfSqFGjNm7ceOLEiREjRsyaNWv9+vWxsbEa9RQINeVSyTTmPo89\nk+wNhdesZVwqIQzo5jHM6dOnq6urx40bt2HDhrvvvvvtt99+6qmnvO/sAeEt+tb7Rsx43332nfju\nUolJ4dA7ba6Qrly5Ul5e7vp17Nixc+fO9fN+p9P5zDPPTJs2bfLkyUKIRx99dMqUKS+++OKf/vSn\nvLw81bsLyES5VBJCuN++E99dKpl/uDP61vv6+FNAakMKpI6OjoULF65atWr58uXer1ZVVZWXlzc1\nNZnN5oyMjPz8/PT0dOWlS5cu7dixw/XOadOm+Q8ko9G4cuVK95af/OQnW7du/eKLLwgkRCZlpoP3\nU6XOI4/FZj7H7Tvo0ZACqaioqLW1taury/ul7du3FxcXx8TEZGdnd3Z2lpaW7tu37/XXX585c6YQ\nIj09/cSJE4F/0YULF5qbm++++26j0ai0GI3G6Ohon18NRAgulRBmBvwMyel0WiyWQ4cObdq0adeu\nXT7fU1tbW1xcnJaWtn///pKSksrKysLCQofDsWXLlu7u7kH00mKxrFix4ujRo66WkydPdnV1ZWVl\nDeLTgHASm7nxlrmfum84Kyg1BH0acCA1Nzfn5uauWbOmsrKyr/eUlZUJIQoKCpKTk5WW3NzcvLy8\ns2fPHj58eBC9nDp1akZGxpYtWw4cONDW1vbnP/+5oKAgPT194cKFg/g0IMwY4lIS7vc904FSQ9CR\nAd+yS0xMfO2115SfT506VVRU5P2e2tpao9E4a9Ys98acnJyKioqampqcnJyBfumwYcPefPPNbdu2\nPfvss729vSaT6cEHH9y6deuwYcP8/2FmZqbr54aGhoF+L6AjylOlrhMvOC7erAtOqSGI74+EMhtw\nIJnN5tzc3Bt/HO3jzx0OR3t7e0pKiqvGj2LChAlCiJaWgP6xtmTJkiVLlri3JCUlvfHGG3a7vb29\nPTk52edXeyOEEFGUvZR8rp+9Zi0zjbkv+tZ7hRDG2BRDXAplhyKH+0goczgFf9p3Z2dnT0/PyJEj\nPdqVlo6OjqF8eGxsrGuqHgBvfZUa6rG1XLPtda/TqgSSEk5CiOhb7zXGKj8wFQLaCH4gORwO4evi\nSWlRXgWgKp+lhjwocdVjaxEXhRCir6wyxCUb4lKUrAJUFfxAUmZmeweP0uKatw1AVX1dKgXCI6sQ\nXvp59K6h4AeSUl/Obrd7tCstVJ8DQkm5VHJcrLl+oVYI0WNrcdpbmHcHOQU/kMxmc0JCgtVqdTqd\n7tdDFotFCJGYmBj0bwTghyEuJSZucUzqYvfGHluLEk5CCLIKklClll1WVlZdXV19ff3UqVNdjceP\nHxdCZGdnq/GNAAZEmWWn/P/fPatcN+tcWaW0kFUIAVUCad68eXV1dTt37nz33XeVlra2tj179hiN\nxjlz5qjxjQCCQpnO4DOrECZeiKRp30KIZcuWlZSU1NXVLVq0aP78+efOnauoqLDZbKtXr05KSlLj\nGwEAeqdKIJlMpt27d2/btq2qqqq+vl4IER8fv3nz5vz8fDW+DgAQBqJ6e3vV+3SHw2GxWOLj4xMT\nE0O/+3hmZiaVGgDAncwDo7ob9JlMph/84AeqfgUAIDzoZgtzAEB4I5AAAFIgkAAAUiCQAABSIJAA\nAFIgkAAAUiCQAABSIJAAAFIgkAAAUiCQAABSIJAAAFIgkAAAUiCQAABSIJAAAFIgkAAAUiCQAABS\nIJAAAFIgkAAAUiCQAABSIJAAAFIgkAAAUiCQAABSIJAAAFIgkAAAUiCQAABSIJAAAFIgkAAAUiCQ\nAABSIJAAAFIgkAAAUiCQAABSIJAAAFIgkAAAUiCQAABSIJAAAFIgkAAAUiCQAABSIJAAAFIgkAAA\nUiCQAABSIJAAAFIgkAAAUiCQAABSIJAAAFIgkAAAUiCQAABSIJAAAFIgkAAAUiCQAABSIJAAAFIg\nkAAAUiCQAABSIJAAAFIgkAAAUiCQAABSIJAAAFIgkAAAUiCQAABSIJAAAFIgkAAAUiCQAABSIJAA\nAFKI1roDng4dOvT+++83NTWNGTPmxz/+8cqVK4cNG6a8dOzYsaKiojNnzowePXrmzJlr166NiYnR\ntrcAgGCR6wrp/fffX7NmTWxsbH5+/pQpUwoLCzds2NDb2yuE+PTTT1esWGGz2Z588snp06cXFRWt\nW7dOeQkAEAYkukKy2+2/+c1vHn300V/96ldKy4QJE37+85/X19dPnjx5x44d2dnZu3fvNhgMQoiJ\nEycWFBQcOXLk/vvv17TXAIDgkOgKqaWl5erVqw8//LCrZdq0aUIIi8Vy4cKFxsbGBQsWKGkkhMjL\nyxs2bNhnn32mTV8BAMEm0RVScnLy3r1777zzTlfLF198IYS4/fbbrVarEGL8+PGul0wmU3Jycmtr\na8i7CQBQhURXSHFxcVOmTImNjVV+/fLLL//lX/7l3nvvnTRpks1mE0LEx8d7vL+rq0uDjgIAVKDN\nFdKVK1fKy8tdv44dO3bu3LmuX3t7ez/44IOXX375jjvuKCwsVFpc/3V/G5MaACBsDCmQOjo6Fi5c\nuGrVquXLl3u/WlVVVV5e3tTUZDabMzIy8vPz09PTlZcuXbq0Y8cO1zunTZvmCqRvvvnmxRdf/Pzz\nz1evXr1u3TqTySSEUC6b7Ha7++fb7fbExMSh9B+ByMzMbGho0LoXusdhHDqOYdgbUiAVFRW1trb6\nvG+2ffv24uLimJiY7Ozszs7O0tLSffv2vf766zNnzhRCpKennzhxwvuvvvzyyyeffDIpKenDDz+8\n4447XO1K8FgsFleL0+lsaWlhih0AhI0BP0NyOp0Wi+XQoUObNm3atWuXz/fU1tYWFxenpaXt37+/\npKSksrKysLDQ4XBs2bKlu7u7r0/u6ekpKCi4884733vvPfc0EkIkJSUlJiYeOHDA1VJdXe1wOO65\n556B9h8AIKcBXyE1Nze7z8z2qaysTAhRUFCQnJystOTm5ubl5VVUVBw+fDgnJ8fnXx09evT06dP/\n9E//VFdX594+efLkW265ZeXKlb/+9a+3b9++ZMmSxsbGV155JTU1dfbs2QPtPwBATgMOpMTExNde\ne035+dSpU0VFRd7vqa2tNRqNs2bNcm/MycmpqKioqanpK5BOnjwphHB/tqR46623Zs6c+cQTT3R0\ndLzzzjvFxcVCiClTprz00kvDhw/339vMzMzA/mfBHw5jUHAYh45jGN4GHEhmszk3N/fGH0f7+HOH\nw9He3p6SkuKRFhMmTBBCtLS09PXJ69atW7duXV+vGo3GTZs2bdiwobW1ddSoUQkJCf12leefAKAj\nwZ/23dnZ2dPTM3LkSI92paWjo2MoH240GtPS0obyCQAAOQV/YazD4RC+Lp6UFuVVAAA8BD+QjEaj\n8BU8SovyKgAAHoIfSD4XsbpaXJWBAABwF/xAMpvNCQkJVqvV6XS6tyvLWqmtAADwSZXiqllZWd9+\n+219fb174/Hjx4UQ2dnZanwjAEDvVAmkefPmCSF27tzpamlra9uzZ4/RaJwzZ44a3wgA0DtVqn0v\nW7aspKSkrq5u0aJF8+fPP3fuXEVFhc1mW716dVJSkhrfCADQO1UCyWQy7d69e9u2bVVVVcqNu/j4\n+M2bN+fn56vxdQCAMBCl6pZCDofDYrHEx8cnJiZGRUWp90UAAL1TN5AAAAiQNjvGqsrPxoAIUGVl\n5SeffOLROGbMmJ/+9Kea9EdHBr1rJdz5OYycnIH4z//8z0OHDjU2No4YMWLy5MlPP/303/3d33m8\nR8KzMdwCyf/GgAhQRUXFwYMH4+Pj3RuZkBKIQe9aCXd+DiMnp39Op/PZZ5/9+OOPzWbzpEmTzp49\nW1xcXFpa+tZbb7lvICfp2dgbRmpqajIyMubOndvS0qK0/Nd//VdWVtaDDz5ot9u17Zu+5OTkPP74\n41r3QjeuX7/+1VdfVVdXb9y4MSMjIyMj44033vB4DydnvwI5jL2cnP157733lMsd13lVWlqakZHx\nwAMPfPvtt0qLtGejKuuQtNLXxoBnz549fPiwpl3TE5vN1tLSkpWVpXVHdKO5uTk3N3fNmjWVlZV9\nvYeTs1+BHEZOzn794Q9/iI6O/tWvfuXaAGjJkiUPPPDAuXPnGhsblRZpz8awCqS+NgYUQtTU1GjU\nKf1pbGzs7e2dOHGi1h3RDWXXSsXTTz/t8z2cnP0K5DBycvrX29v7zTffZGZmjhs3zr399ttvF27b\n0Ul7NobPM6RBbwwID8rGhrfddtvOnTv/93//NzY2NjMzc8WKFbfeeqvWXZOUertWRpR+D6Pg5OyP\n0+ncunWr9/yFpqYmIcT48eOF3Gdj+ASSqhsDRhTl//MbN268du3a+PHjrVbrJ598smfPnsLCwhkz\nZmjdO13i5AwWTk7/oqOjH330UY/GI0eOHD169M4777zzzjuF3Gdj+NyyY2PAYGloaIiKinr66ac/\n//zzjz766C9/+cvGjRs7Ozt/+tOfXr16Veve6RInZ7Bwcg7UgQMH1q1bFxMT84tf/MJ9szo5z8bw\nuUJiY8BgefXVVx0Oh2serdFofOaZZ5qbmz/88MP//u//9v73F/rFyRksnJyB6+zsfOmll/74xz/e\neuuthYWF06ZNU9plPhvD5wqJjQGDZezYsd6rOpQK7l9++aUWPdI9Ts5g4eQM0MGDB//hH/7hj3/8\n4/z588vLy91XIMl8NobPFZL7xoDuIc/GgAPV09MTFRXlUXvQbDYLIS5evKhRp/SNkzNYODkD8c47\n77z88svJycm7d+++7777PF6V+WwMnyskwcaAwWCxWLKysp5//nmPduVhsvJQFIPAyTl0nJyBOHjw\n4Msvv3zPPfd8+OGH3mmkkPZsDKtAYmPAoUtPTx89evTBgwdPnz7tarx8+fJbb70VHR394x//WMO+\n6Ron59BxcvbL6XS+/PLLCQkJb775pnLh6JO0Z2P43LITbAwYDFFRUT//+c83bNiwbNmyxx57LDMz\ns62t7b333jt//vxzzz13xx13aN1BveLkHDpOzn41NTV9/fXXaWlp//qv/+r96lNPPZWSkiIkPhvD\nKpDYGDAocnNzf//73//mN795++23hRBRUVFpaWm//e1v+RfoUHByBgUnp3+nTp0SQvz1r3/993//\nd+9XH374YSWQpD0bw3M/JDYGDIrLly+3tbWlpaV5VFbGUHByBgUnZ1DIdjaGZyABAHQnrCY1AAD0\ni0ACAEiBQAIASIFAAgBIgUACAEiBQAIASIFAAgBIgUACAEiBQAIASIFAAgBIgUACAEiBQAIASIFA\nAgBIgUACAEiBQAIASIFAAgBI4f8DO8Ej/RcfVosAAAAASUVORK5CYII=\n",
"text/plain": [
"<IPython.core.display.Image object>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"for mu = [4 4.5 5.1]\n",
" [L,U] = lu(A-mu*eye(3));\n",
" v = rand(3,1);\n",
" for k = 1:20\n",
" v = v/norm(v);\n",
" evec_err(k) = min(norm(v-x),norm(v+x));\n",
" v = U\\(L\\v);\n",
" end\n",
" semilogy(evec_err)\n",
" hold on\n",
"end"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Rayleigh quotient iteration\n",
"\n",
"The coup de grâce: a better shift means a better eigenvector estimate, which means a better eigenvalue estimate, which is an even better shift, etc. This feedback loop causes *cubic* convergence in the hermitian case, which is rare--and fast! It's so fast that we pretty well need extended precision to see it."
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"lambda =\n",
" \n",
" 5.21431974337753518741549770084858048890791963721949943433138\n",
" 2.46081112718911088347412409730147999190011289045787329828077\n",
" 1.32486912943335392911037820184993951919196747232262726738785\n"
]
}
],
"source": [
"digits(60)\n",
"A = vpa(A);\n",
"lambda = eig(A)"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"error in eigenvalue = -0.21431974337753517678706316473835613578557968139648\n",
"error in eigenvalue = -0.00120498927917453164254668607213716313708573579788\n",
"error in eigenvalue = -0.00000000019350335636704886234483316195391670519466\n",
"error in eigenvalue = -0.00000000000000000000000000000086272041087647542775\n",
"error in eigenvalue = 0.00000000000000000000000000000000000000000000000000\n",
"[\bWarning: The system is inconsistent. Solution does not exist.]\b \n",
"[\b> In symengine\n",
" In sym/privBinaryOp (line 937)\n",
" In \\ (line 335)\n",
" In pymat_eval (line 31)\n",
" In matlabserver (line 24)]\b\n"
]
}
],
"source": [
"v = vpa(ones(3,1));\n",
"for k = 1:5\n",
" v = v/norm(v);\n",
" mu = v'*A*v; \n",
" fprintf('error in eigenvalue = %.50f\\n',min(mu-lambda))\n",
" v = (A-mu*eye(3))\\v; \n",
"end"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Count 'em: The number of leading zeros *at least triples* with each iteration, and in 5 iterations we get 60 digits of accuracy. 😎"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Matlab",
"language": "matlab",
"name": "matlab"
},
"language_info": {
"codemirror_mode": "octave",
"file_extension": ".m",
"help_links": [
{
"text": "MetaKernel Magics",
"url": "https://github.com/calysto/metakernel/blob/master/metakernel/magics/README.md"
}
],
"mimetype": "text/x-octave",
"name": "matlab",
"version": "0.11.0"
}
},
"nbformat": 4,
"nbformat_minor": 1
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment