Skip to content
Permalink
master
Switch branches/tags

Name already in use

A tag already exists with the provided branch name. Many Git commands accept both tag and branch names, so creating this branch may cause unexpected behavior. Are you sure you want to create this branch?
Go to file
 
 
Cannot retrieve contributors at this time
{
"cells": [
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"from numpy import *\n",
"import matplotlib.pyplot as plt\n",
"from pretty_plots import setdefaults"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [],
"source": [
"setdefaults()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Numerically integrate Euler-Bernoulli equation\n",
"$EI\\frac{d^4w}{dx^4}=0$\n",
"\n",
"Boundary Conditions:\n",
"\n",
"1. $w(0)=0$ \n",
"\n",
"2. $w'(0)=0$\n",
"\n",
"3. $w''(L)=0$\n",
"\n",
"4. $EIw'''(L) = -F$\n",
"\n",
"Using constants\n",
"- L=400 mm\n",
"- E=200e3 MPa\n",
"- t=3 mm\n",
"- w=12 mm"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [],
"source": [
"def int_N(N,F):\n",
" \"\"\"numerical integration setup and solution \n",
" for EIw''''=0\n",
" BCs \n",
" w(0)=w'(0)=0\n",
" w''(L)=0\n",
" w'''(L)=F/(EI)\n",
" inputs: N=number of nodes (1-N), node 0 is x=0mm\n",
" F=applied force at x=L\n",
" outputs: w=solution of w(x) at each node\n",
" K=finite difference matrix w''''\n",
" vF=applied forces/EI\n",
" \"\"\"\n",
" L=300; E=200e3; t=3; w=12\n",
" I=w*t**3/12\n",
" dx=L/N\n",
" d1=-4*ones(N-1)\n",
" d1[-1]=-2\n",
" d2=ones(N-2)\n",
" d=6*ones(N)\n",
" d[0]=7\n",
" d[-1]=1\n",
" d[-2]=5\n",
" K=diag(d)\n",
" K+=diag(d1,1)+diag(d1,-1)\n",
" K+=diag(d2,2)+diag(d2,-2)\n",
" vF=zeros(N)\n",
" vF[-1]=-dx**3*F/E/I\n",
" w=linalg.solve(K,vF)\n",
" return w,K,vF"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Demonstrate qualitative convergence between 4, 24, and 44 nodes. The beam shape looks similar between each solution**"
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAakAAAEeCAYAAAAw4+qWAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xt809X9+PHXSXovLWkphQIKBBTxgvaCIuKGWtAN3URb1G1O52zr1Ilzs8Vtun51U6ibG9t3363F6XS6Ca2XL8q+U8rmb264OVrUeZmXFq+0UKDp/Zbk/P74JOkt6b1Jmryfj0ce0E8+SU5P27xzznl/3kdprRFCCCGCkSnQDRBCCCF8kSAlhBAiaEmQEkIIEbQkSAkhhAhaEqSEEEIELQlSQgghgpYEKSGEEEFLgpQQQoigJUFKCCFE0IoIdAOmupSUFL1gwYJAN0MIIaaUqqqqI1rrmcOdJ0FqnBYsWMC+ffsC3QwhhJhSlFIfjuQ8me4TQggRtCRICSGECFphH6SUUhlKqXKlVKPrVq6Uygh0u4QQQoR5kFJKZQNVQA5wzHXLAaqUUjmBbJsQQogwDlJKKQuw2/XlGq31Iq31ImCN61i56xwhhBABErZBCsh3/Vuita50H3T9v2zAORNuV+0u1lasZdkjy1hbsZZdtbsm66WEEGLKCucUdPeIabuX+0oxAtQaoGSiX3hX7S7u/Pud9Dh7AKhrq+MHe3+A/cO9fPHVZ9hlb2TrjGTqzYrZ8WlszNjIOus6z2O3Vm+lvq2e2fGz+90nhBChRoXr9vFKqRrAqrVWPu7XQK1rCtCnrKwsPdrrpNZWrKWurc7rfZFOJz1KgeptlhnF55NOIyE2hYr6l+h2BTeAGBVJcUs36xo+YdfMeWxNslDf0ywBTAgR1JRSVVrrrGHPC+MgpQGGCVJe71dK5eOaCjz++OMzP/xwRNekeSx7ZBmaiev3OKeTC9ra+VN8PD2m3ubGmGMoXlksgUoIEXRGGqTCeboPwDbMfV4TJ7TWZbjWrbKyskYdbWbFz6K+rd7bE/cbQY1Uu8nEswnTBh3vdHRy795iktqb+Eg5eejtR41pwshENjbaWNfwCUyfBxfcBcs2jPp1hRBisslIagwjqb7GMt23q3YXxXuL6XR0eo7FaM1dDUf5ebKF+ojBnx3iMOHQTrpGH8O8inJqio4eY0NrG7sSLWydfZxMEwoh/GakI6lwzu6rnaBzRm2ddR3FK4tJi09DoUiLT6N4wXou6dbcesxGjNPZ7/wYcwx3nXsv//WZzcSYY/rdF+nUfK61lbgBjxlOt0lxz8wZrD5uDt9LTqCupwmNpq6tjuK9xZJtKIQICuE83VcLWJVSGVrr6r539Kk4MdR04Liss64bPFpJPo11e+6GI8d8ZvcBvdl9kYlsrP+Ydc02dsV3UpySTKep93NHBIoMYnlFt4OPEdhRL6M29zTh4kgL73Yf4xf7fyHZhEKIgAjn6b58jFTzEq110YD7tgCFQIFr/cmnsUz3TajXd8Ceu6HJd3afz2zCkayBafoFuBink+J2xbpzZR1LCDF2kt03Au51JyDTPZpyjaKqYPj1KAiCIDUCXtfAnE6+d/QYP0uyeB1NDSXe6eSnDU0czrqaXx6rklGWEGLUJLtvZAowRlNVSil31YnsPveFBHfg6DdN2GhjXWsHkbEzKE6IpFP3XnsViYkl2swbdHsdabWZTOTPSoKPnvXc717L6vt6QggxXmE9kgJPkdktgHsdqhoo6lsqaShTYSQ1HF9VLNaWr6Gu3Uuq/BASzTE8+fnHqbK9J5UxhBA+yXSfn4RCkPLF2zRhlFOT2dnJy7ExQ65nKeh3uXKM1hQvWM+61fdMXoOFEFOGpKCLcfOkykdOR2lNWo+du48cpexQA2mOoVPeB3706VSKn9Y8aSR6CCHECMlIapxCeSTVT58sQqbPY1f6eoo/+VO/UVYkJhY4Ne8pp/dRltas7Oxh7vzP8lJrLYfaD8tUoBBhSqb7/CRsgpQXvtaysrefx6HOIyN+HqkxKET4kSDlJ+EcpHzZVbuL4r8W0dlnNGXSGucQa1jTTFE89bnHSEtZ6o8mCiECTFLQRcCss66Dj/7B1ponqTebmG13sLHRRoYd1s5J8fqYVmc3Fz6Xy2IdQYO202SC2U7YaJVkCyHCmSROiEmxbvU9vLC8mNcb4YVP6lgXMYO0dVtJi0/z+RitFO+ZHNjMCq0UdWbFDw48za4X7/Rjy4UQwURGUmLyLNswqHTSxmnxg9LaI5WJBQ54TzkGJVx0mYxAlTTjBFac+hVMSj5XCRFOZE1qnGRNavR8JVws++2p6CHWrWZoRU9ENC2OLskKFGKKkzUpEbS8VoAHZto1hyN9B6mjSoNrBFbXVscP/n6X5/mEEKFJ5k5E0JjffOagvbSinJpT2xTKy4i/y9lN8Uvf5dXdm+Cnp0KxxfhXLhgWImRIkBJB4wvn/ZATGrJI7XGitCa1x8mi+nRe/ug+tI8NsTpxcvXBXXx+mp3PHD+HZUmw9l/FkmwhRIiQ6T4RNC5Nnwvcy/3Pv0OLrYMESyyf+8xCZn1k48+tFkxRvveg/Dgq0vP/uggzxR88DbUrZCpQiClOEifGSRIn/GPrP7bzm/+UoFW355hymsnsamNfTLTXMkwJpmh251YSH2PxZ1OFECMgiRMipGxccQXWlHi2/PMBmnoacPZY6Gq4kHs6H+bz1sEFbQFanF2s/cO5fGXW2aQsXM22N38rW4cIMcXISGqcZCTlfy2dPfz3n9/nob8f4HP6Jd62PkN95DCft7TuN9qSeoFCBJZs1SFCVkJMJHd8fikvfOuztC+5jKSGLKKd/T9sRaBI6rudyIDpwE5HJ1urt/qjuUKIcZAgJaashSnxPHjNcm6+/AGi267F2W1Ba3B2W2j5dAMnOH7Od6evMkZRXtS1HqTp7mRJWxciiMl03zjJdF9w6HE4efwfH/LA7ndp7rR7jptNipSTSmjXR70+Lt7pZEV7J2/ERHM4wszs+DRZrxLCD2S6T4SVSLOJa89ZyIu3n8fVK+Zjcs3uOZyao59kgzPS6+PaTCb2TIvjUIQZjVHJonhvMbtqd/mv8UIInyRIiZCSHB/FPZeeyh83nsvKRTMAsDen01F3Gc5uC2hI67FzVVMz1u4er89hrFf9zJ/NFkL4INN94yTTfcFLa80Lbx3ih7ve4uNjHZ7jf4u6hXmmIziAMxYc52Ore7j/lDzWZt0sldeFmASyM6+fSJAKfp09Dh76+wH++8/v097t4Aumv7E58kHiVDdr582hboj09TlE0BkVR2N3i1xfJcQEkjUpIVxiIs3cuHoxf/nOai7PmMdO5yo29VzPJ84UvnmsaVD6et9swIPYOdbdjEbLepUQASAjqXGSkdTU8+rHNv7r2TfZ/5FRCzAicT/RM5/HFGljRswsblqSy8f/eYqHuz71OhU4y26nstkMF9w1aFNHIcTIyHSfn0iQmpq01vzvqwe57//e5lBzV7/7vnjGHDZ97iQufOZMr+WW0JqNjU1c3d5D9CU/l0AlxBhIkPITCVJTW1uXnV//vxpK/1pLt723QkVspBnLkhJaHQ0+H2uxO8BkoslkkvUqIUZJ1qSEGIH46Ai+vXYJe277LOtOS/Mc7+hx0PDR+YOur+q7+aItwozNpGS9SohJJEFKCOC45Dh++eUMnshfwdK0RKD/9VVKa9J67Pyw4SjfP3IMk5cZiE5HJz/d94C/my5ESJMgJUQfK6wzeO6bq7h3/Wkkx0dhb06nrWYTq9+9hGc+PsIX2tq5oqXV+1oVcKj9EI//uQi7w/uFwkKI0ZE1qXGSNanQ1dTRw8/3vMcjez/A7tR8wfQ3CiN2MEcdZe3x8zgU4X1Le4BUbcIRncCx7mZZrxLCC0mc8BMJUqGvpqGVHz73Fn95pzeJIiJxP7FzngLVO2JSWqO9Va8AYpxOitsV686VtHUhQBInhJgwi2ZO4+GvncnD1y7HOjMecK1XHeytBzgzZhb3nPldbo0/qV9yhVunycTWaAc8e4tsCyLEKMhIapxkJBVeuu1OHn35A7bueY+WPluCRJgUV589n1svOJFzyzN8Xl/1eN0hlsXMgm+94a8mCxGUZLrPTyRIhaejrV38ZPe7/OGVj/rtqZgUF0ncos00271fX6W0ZnlHJx+mzOdwe4OsV4mwJdN9QkyiGdOiuXf9aTz3zVWctTDZc7yxvYdDH56P0gP2r3JFMq0Ur8TFcqj9sFxfJcQISJASYhxOmTOdJ/JX8KsvZzDXEgsY61XtBy8juifWc33V7cca+Wx7h9fnMPav2urPZgsxZfjeo0AIMSJKKT53WhrnnZTKgy/V8su/1NDRnM6R5nS+YPobRZE7mKPa+DIWzojz/hx1bXXYHT1EmL3vICxEuJI1qXGSNSkxUH1TJ1v+9B+e3v9pv+OpCdFELbyPpp7DXh93ElEUn3sfp1jX+qOZQgSUJE74iQQp4Uv1R43817Nv8drHNs+xiMT9xM15Gq26vT9Ia+KVmXbtYLbdwcYus1xbJUKSJE4IEWAZxyfx9DdW8pPc05mZEA2416vWG9dXAakxqVwUPad340WlaMOJVoq6yAiK4zS7Km+Xa6tE2JKR1DjJSEqMRGuXnf/5y/s8+NIBuh29W4LERZm56bzFXDj7Ha745ya6vRSsmG23s7vZLNdWiZAiIykhgsi06AgKLzqJyts+y0WnzPYcb+92cP/z73Dtzli6fZRUqjeb+Vv3EX81VYigIiOpcZKRlBiLve8f4e7n3uI/9S2eY/GLNmOKsvl8TGbkDD6JjOBw+2G5CFhMeTKSGoZSqlwppYe45Qe6jSJ0rVycwnPfXMU9l55KUpyRdt7VcCEmp7n/iX0+RFb1HOVQ+yG5CFiElbANUkCG61+bj9uxALVLhIkIs4mrV8znL99ZzbUrF6BbM2iry+l/EXBTD2t7zF4fLxcBi3AQzhfzWoFarfWiQDdEhDdLXBTFXziFL591PHc/N4OX3ksHoBm4C1icEgszN3p9bF1bnf8aKkQAhOVISillcf23NqANEaKPE2Yl8Oh1Z/Kba7JYMKO3NMX7Rzo8Keve/OipHNq7mv3RRCH8LlxHUu7FuuqAtkKIAZRSXLB0FueeMJPf7j3Az/e8T2uXna6GC4lJewplGrwt/RMt7/DC789Bac0xE8x2wkbretatvicA34EQEyssR1IYU32AJ4Gi0ZUsUaWUKhzuwUqpfKXUPqXUvoYG71syCDEeUREm8j+ziL98ZzVXLj8OR0s6nXXGJotag9k+jRPsvSnrx0xw1KyMi4DNiuIDT7PrxTsD+B0IMTHCMgVdKbUF6BuMqgELvcGrWmudOZLnkhR04Q9vfNrEfz37Jv/6oNFzzISD9XMeY3fiW163rU9zaF64Ti4AFsFJUtCH5s7sq9BaK611piuBYhFGwMpwBTIhgsKpc6ezo+BsfnFVOnOmxwDgxMyTB6/x+Zg6Ezi10+f9QkwFYRmktNZrXMEpd8DxWsB9bNhpPyH8SSnFJafPYc+3V3Nr9gnERBp/vjPtPmZDlOLGxz/D4cYaP7ZSiIk1Jaf7lFKj/auzjXT6rs/zW4FFrsDlk0z3iUA5aOtg8//9h/oPHuDA7Go6Td4/c8Y6ISZ6GraeNqlUIYLGSKf7pmR2nx+vbZILekXQmmOJ5edXpbNw05UsB47OrKYhQpFi16R1mfl3vFFNvcMEHT2tAJ5KFYAEKjElTMkgNR6ua6QaMUZXST7ut7ru911ITYggMccSyyu2K6HpSsC4CLgWWBX/LP8+7m84ByRVuCtVSJASU0HYrUm5Ak8lYFFKlfa9zxWg9ri+LPJ324QYi9svXEJs5ODSSX9ruwQn3iur17XV4XQ6JrtpQoxb2I2kXHKBKiBfKbUB2Ack0z/rryxQjRNiNC5NnwvA/c+/w0FbB6mJ0ZiV4mBTJ84ei8/K6jc8ciY/snUws+kgTJ8HF8gOwCL4TMnEiYniunB3Db0VKPYBpVrripE+hyROiGDU0e3g9orX+NMHf/RZqQIg3ukkyqmxmU3MdjjZuOhyqVQh/EKukxoBrXWJKx09yXVbM5oAJUSwio0y84ur0vn2yqvo6lOpgp7pnOlMRLk+nLaZTDRGmI1KFRFmij94Wrb/EEElrEdSE0FGUiLY/eWdw9zyh/20dNo9x56wfI2C2amDkioA0uLTeCHnBX82UYQhGUkJIQA4b0kq/3vTOSyaGe85Nq89AV8fT2X7DxFMJEgJEQasM6fxzE3nkL00FYAS+wZm2X1n9xVXfJGO7lZ/NU8In2S6b5xkuk9MJU6n5meV7/LzP7/PmdOfoHb2frpM3tPUU7UJHWvhSGejVKoQEy6kK04IIcbGZFLctnYJJ89J5LYdZrrrlhA983lUpI0EcxJWRxuvq24ADisndBpFV6RShQgUCVJChKGLTk1jQUo8+Y9G81GNsV19K3Da0hl8P2kbPzz2T5BKFSIIyJqUEGHqpNmJ7Lz5HFYtTvEce+Htozz4/lcHBSi3ekmqEH4mQUqIMGaJi+K3X1vO9asWeo69d7gVeizeH6A1f/z5iVBsgZ+eCq/v8FNLRbiSICVEmIswm/j+xSfzwIbTiYow3hI6Dl+IdkYOOlcrRdH0aH44w0JX88fw7C0SqMSkkiAlhADgsox5VNxwNmnTY7A3p9PpqlQBkBKZyAx778XA2xMTWDH/OJbNTWHtvrulSoWYNBKkhBAey+ZZ2HnzKrLmJ2FvTqetZhMtb28m7shmnv2knrWtbZ5z7UoZ5ZTMiuK9xRKoxKQYc5BSSp2vlLpeKXWfUupXSqnvuL4+YyIbKITwr5kJ0fw+bwVfOut4z7E3DzbTopP5ccNREh2DLwJ2Z/4JMdFGlYKulLoeKKB3SwtvKUBaKWUDdmNUFP/L+JoohPC3qAgT964/jVPmJFK88016HJrN3RvYHPkgLT62qa9rq+No88fMSDzOz60VoWzYkZRSKtE1WnIAZUAS8CRwA8a+TJmu2xrX1/dj7NW0AahUSr2nlFo/Se0XQkyiL581n9/nrSBlWhQ7navY1HM9M+2+q9RseHIdr70n035i4gxZFkkpdRnwoOvLMoyR0YERP7lS2Rgjr8uBF4ANWuvmsTc3+EhZJBEODto6KPhdFf/+tImIxP2D96jS2nNtVYTWfD4lg3911lPfVi8llYRXE1UF/UGgSGudrLXeNJoABaC1rtRa52LsevsBcMdoHi+ECA5zLLGU33A2l6XP7Z/5p2FG9CyunXkWiQ4nYCRU7Dy6n7q2OjTaU1JJEivEWEiB2XGSkZQIJ1prfvO3A9z7x7dxut463OtXKyxvcduL3+ZtHyvdsk+V6Cvg+0kppRIn67mFEIGhlOL6c608et1ZTI81Lvbttjv5TvlrPPjmXB7Ked7nY+tbD0qVCjFqowpSriSI00dw3uXAqKYGhRBTx6oTUnj25lUsmZXgOfbbvR+Qt/0jZsXO9vqYaU4nukmqVIjRGe1IahFQrZS619udSqkFSqnngR2Aj+JfQohQcPyMOJ66cSWfO7U3KL1ce5TmujVEmaIHnd9iNvOt1BTa7J2w525/NlVMYaMNUouBV4FNSql3+46qlFLfAWowUtH3uM4VQoSw+OgIfvmlDL695kTPscN1p9Bx8DIskakorYlyOj337YmP4+J5aZyf6GDZI8tYW7FWEirEkMaUOKGUygd+DWiM1PQsjGulbECe1vrJiWxkMJPECSEMlW8d4tbtr9La1Vvj77XE24jrrueBZAuPTfe+TB1jjqF4ZbGkqIeZSU2c0FqXYUz9NQH5GBUoqoAF4RSghBC9sk+exTM3rWRhSrzn2J2tl+FU0RQds3Hf4SPG9VQDSEklMZQxBSlXfb4XMNadDmAEqwygTLL6hAhfi1MTeOamc1i9ZCYAO52ruL3r69Srmaxr6/BaRw2gvq3ef40UU8pos/sSlVK/whg1LcK40HcxYAWewiiFdEAp9fUJb6kQYkqYHhvJb65ZzjdWLwKMQLWiYyvL9HYs0bO8PsasNR8fes2fzRRTxGhHUgcwyhztBxZpre8H0FrbXJUlNries0wp9cqEtlQIMWWYTYqii07iF1elExNpvM20dNmpO3AeEWpw5p9dwZV//DJ733jM300VQW5UiRNKKSdQqLX+8RDnTMcop3SZ1to8/iYGN0mcEGJobx5sIv/RKj61dQAQkbgfy9xKujnGdFMULfZOHCbXRKDWJEbE0uLokpp/IW6yEicWDRWgALTWTa5R1dpRPrcQIgSdMmc6O28+hxXWZADszekceft25jT+N79f91d+t2wjqa66fyhFs6NTav4Jj+GqoCdOZNXyiX6+YDDakZTWmpaWFpqbm2lvb8fhZQM5EbrMZjNxcXEkJiaSkJCAUr5SCUJPj8PJj3a9zW/3fuA5ZomL5JdfyuDE6Bou3HM9PV66I82heSHrLli2wX+NFZNupCOp4TY9tCmlNgObxxNclFLnA1uASsK4ErrWmsOHD9PW1kZycjKzZ8/GbDaH1RtVONNa43A4aG1t5ciRI3R0dJCamho2P/9Is4niL5zCyXMS+f7Tb9DtcGJr7+GrD73Cdz+/FLtSGJde9ldvwiilBBKowtBw032LMabtGpVSTyilzhvpE7tKJH1HKfUexi69lVrrsA1QAC0tLbS1tTF//nwsFgsRERFh8wYljOKsERERWCwW5s+fT1tbGy0tLYFult9tyDqOJwpWkJpgJFA4nJp7nnuLaJXs9XwN/DFKSSmlMDWixAlXhYlCjFRzDVQD+zDKINlcp1mAGa5zsl1fK6ACI1U9JAvOjma675NPPmHatGlYLFLWUIDNZqO1tZV58+YFuikBcai5kxseq2L/R8ZbSETifuLmPI1W3V7PP6+tnf+kLpKNFEPESKf7Rpvdl42xRfwFGMHIl0qM0VOZ1rppxC8wBY0mSL377rtYrVYiIoabZRXhwG63U1tby4knnjj8ySGqy+7gzmfeYMe+TwAjUM1P/QMNEYoUhwMNHPHx9yLllKa2iVqT6kdrXYkRgNyp5lbXLRmoBY5prfePvrnhweFwYDaHfFa+GCGz2Rz2iTPREWa2XL6MU+ZM5+7n3sLenM4prW1sjnyQONVNi1IUpqbwt7jYQY91l1OSIBXaxvyR3jVC2u+6iRGSNSjhJr8LBqUU16xcwImzErjp99XsbFsFPVAYsYO5pqP8oiOK9FjwVlNJyimFvtGWRbpsshoihAhvZy+awc6bz+HktER2OlexqvvnLOx8nC/F/YbUOO8bKU6PGDzCEqFltBfzViilHEqpV5RS33YVmhVCiAkxLymOJ7+xkouXpXmOvfLBMZoOZnvdSNFmb+dnz16L0xne06ahbLRBag/GoDsLuB+oUkoddaWnr5cK6GKsKioqUEqhlKK6unrIczMzM1FKYbPZhjxvrEpKSsjMzJyU5xbDi40y84ur0im66CTcM6JH6k+l/dP1TI9MRQERfRK+fnOsiqseP5c15dmykWIIGlWQ0lqv0VqbMDY4LMKofK4wCss+iXE91b+UUvfJKEuMVV5eXsBeu7KykqKiomEDpZhcSim+sXoRD1+7nIQYY+m8w3YGn7x+G1ekbOfFdc/wGWeU5/y3nC3Utx/qLaf01yJ2vXhnoJovJtBYNz3cr7W+X2udq7VOxti2owAjUGViXFMlVVfFmFRXV1NSUuL317XZbOTm5vr9dYVvq5eksvPmVSxOneY5tu2lA9y88yh3r9/NNRHet/7oVIqtNU/C6zv81VQxScYUpPpyjZguB3Iwrp8CY3Q17uceD6VUvuu6LjGFZGcbP7KioiJqa2v9+tq5ubmTNoUoxm5hSjxP37iS7KW9Aeml945w2bbXuOT8p3w+rt5skioVIWDUgUQpdYar3NHzSikHxgaIJcAaoBEow7jgN2lCWzq6NlqAUlc7hjovQylVrpRqdN3KlVIZ/mml8MZqtbJlyxYACgoK/Pa6JSUlVFZWUlpaitU61HXqIhASYiIpuzqTWy44wXPsw6PtrP/VXixRqV4fM83pRDd94q8mikky2hR0d1DaghGUDgDbcAUlrfVirfUNWusnA1xpYttwJ7hGWVUYI8BjrlsORjJIzuQ2z3+e2f8p52z+Mws37eKczX/mmf2fBrpJwyosLMRqtVJZWUlFRcWkv151dTVFRUVkZ2eTn58/6a8nxsZkUty25kR+/ZUM4qKMi+Lbux3UHTiPaOfgyjktZjPfn5VKT0+Xv5sqJtBoR1KK3kvqtgD5QRKUUErlKKVKlVKNGMFmqHMtGGWbANZorRdprRdhBF6Actc5U9oz+z/ljqf+zae2DjTwqa2DO57695QIVOXl5YCRRDHZU3C5ublYLBbPa4rgdtGpaTx94zkcnxwHGPtTWevTmW13oLQmyun0nLszLprLfn8O2TsukMy/KWq0FSc2YRSPzcLI7it0XTVfjVEu6V8Y1c4DsWfUHcBIp+rcH5dLXKWeAKPsk1KqzHV/PsY0pt8s2DT5fzwdPQ5u3f4qt25/dUKf94PNE1uaJiMjg/z8fMrKysjLy5u0AFJQUEBtbS27d++Wwr9TyJLZCey8+Ry++Yf9vPTeEV5pupIvtMzjsehykvQR7pmZwjPTjAt9P6ALOg4DeDZSBKSc0hQx2tp9JbjeuJVSCzFGHmuAdIygpV331QK7tdY3Tmhrh26b58IW13TdUO9q7hHTdi/3lWIEqDX4OUiJ/kpLS9mxYwcVFRVUVlZ6kiq8KSkp4ejRo8M+5/Lly8nJMQbaFRUVlJWVUVhYOORzi+BkiYvi4WuXs+VP/2HbSwfY6VzFzo5VJMZEsPW805j/1i1sbXtn0OOk5t/UMp7afQcwkiTKwBO0ijDe4BdhFJ71W5AaJSuA1nrQxTBa62rX6NDn6rlr65J8gOOPP36SmigAtm3bRm5uLgUFBdTU1Pg877777hvRtGBOTg45OTnYbDby8vL6JWqIqSfCbOJ7607m5DmJbHry33TZnTR32vn67/ZTeNEWaLvU6+Ok5t/UMa49I1w77q7BmAIcONW2ZzzPPclGkr7l8xyttSc4Z2VljXyvk2FM9JSZe02qo6e3ZEwOT+aAAAAeP0lEQVRspJn7LjuNS9PnTuhrTRZ3UKmoqKCoqMhnQGlsbBzV8+7YsQObzUZycjJr1qzpd5879d193J1UIYLX+vR5LJ6ZQP7v9lHX1IlTw+b/+w8pS2fQxeARdrwpEq21FPmdAkab3Xe+q5rEv1yZfrsxRk+ZGNXQ78dIRDBprddOfHMn1FAfu0PiYplL0+dy32WnMdcSiwLmWmKnVIBy27bNSNYsKSmZ8EoQtbW1VFZW9ru5ub+Wa6emhtPmTWfnzatYvqD36hfbp9koHTXo3FZnN98rv5junk5/NlGMwWiz+yoxqklkAk0YFSbc6edZWutNWutgHkH1NdQqecisoF+aPpe/bzqfA5vX8fdN50+5AAVgsVgoLS0FfGf7JSUleWr/DXVzV5TIz89Ha+315r5Oyv21ew1LBL+ZCdE8fv0KvnyWMQ1vb06n/eB6sBuBK6rPvMezHR9R8MR5NLXK1F8wG+103x7gBaBiIreDV0r5XmzwztY3UWIMahl+ys+/5Q7EkPLz8yktLfU5krrjjjtGnDghQltUhIkfrT+NU+ZM5wc736CnOZ2W5nQiTIrCi+ZR88FNPKmNDzr7nK18tiIbJzDbCRut61m3+p7AfgOin9Fm960Z/qzRc12j5E+1gFUplTEweaJPxQmZ4wky5eXlLFrk/VelsLDQz60Rwe5LZx3PibOmccNj1Rxp7cLu1Nz5x4+5ImMrG033ejL/HK51qTozFB94GkACVRAJaH29AHKnp1/h5T73sVI/tUWMkNVqlWAkRiVrQTLPfvMcls2b7jm2vbqOPx7+FhaHc9D5nSbF1tqn/dlEMYywDFKu7DwwLkb2ZCW6/l844BzhBzk5OWitPWtPvmzZssWzVjRZF9/W1NSg9YQlbYoAS5sey46Cs7msz3ps9Uc2mkzeM/vqTMjPP4iEZZBycVcvrVJK7VZK7cao5df3PiFECIiJNPOTDadz58Un445NM+0+ApFSbH7mChwOu/8aKHwK2yDlGimtwSjplO26VWOk0MsoSogQo5Ti66sW8uh1Z2GJi2RGQwYxzsFTfgC/b36bW7evob2rxc+tFAON62LeYKW1rqC3EO5Q51VipNMLIcLEqhNS2HnTKs7/iZ0M4OjMahoiFCl2jcUB78UYn91f7DnCqidWYkcxO342GzM2SimlAAjJICWEEEM5fkYcDqfmlaYroelKAJqBSLq5YPb9/D3JGEH1ALi3pJfCtAERttN9QojwNscSO+hYD1HsafgeiSpy0H3uwrTCvyRICSHC0u0XLiE20jzoeI8Dmp09Xh9T11Y32c0SA0iQEkKEpYG1LS2xvaMnZ4+Pyxu05tHKb/ungQKQNSkhRBi7NH1uv3qWf3qjjlv+8CpdDRcSk/YUyjRgRKUU93/6AnW/PJHbGw5jmj4PLrgLlm3wc8vDh4ykhBDC5aJT03j4a8uJ7syis+4ynN0WtIakyGTm9+54w2PTojlr/lyWJcHafxWz68U7A9foECcjKSGE6OOcxSn8Pm8F1z5sorEmHYBus4miS+ewZ18Ou+NiAOg0GZ/x6yLMFH/wNNSukMy/SSAjKSGEGOD04yyU37CStOlGQOp2OPnmU59QcqiBeC8XAHcqJZl/k0SClBBCeLE4dRoV31iJdWY8AE4N9c4ZtPvYzbdeMv8mhQQpIYTwYa4llvKCszltrlFFvcS+gVl2h/eTtWbvG4/5sXXhQYKUEEIMYca0aH6fdxZnW2ew07mKpIYsop2Di9Nqpbhp32aefXlLAFoZuiRIiaBQUVHh2eLd1+67bpmZmSilvG4jHwglJSVkZkoJyFCWEBPJw19bztqTZ/FK05XY6q7E2W1cS5USOZ1E1zqVXSm+++5jrPxdFsseWcbairXsqt0VyKZPeRKkRNDJy8vz6+tVVFSQmZlJUlISSikyMzMpKSkZ0WMrKyspKioaNrCKqS8m0sz/fDmD3Mx52JvTaavZRMvbm0lr+zGPnfs/LO6z9UeLswvdp+afBKqxkyAlgk51dfWIg8R45ebmkpubS3V1NcnJyWRkZFBdXU1RUZHPrerdbDYbubm5fmmnCA4RZhMlOcvIO3eh59jemqNs/JOZrWv+QJSXLaqk5t/4SJASQSU7OxuAoqIiamtrJ/W1KisrqaiowGKxUFVVRU1NDVVVVTQ2NpKdnU1tbS0FBb73v8zNzQ2aKUfhP0opvvv5pRRetMRz7PVPmri2/Cg9PjP/6v3VvJAjQUoEFavVypYtxsLzUAFiIri3qt+yZQsZGRme4xaLhfLycgDKyrzvf1lSUkJlZSWlpaVYrdZJbacIPkopbly9mHvXn4Y7LtU2tIHde80/s9bUH33Hjy0MHRKkQt3rO+Cnp0Kxxfj39R2BbtGwCgsLsVqtnpHOZHGvI7lHb31ZLBZP8Bk4WnJPB2ZnZ5Ofnz9p7RPB70tnHc8vv5RBlNl4K20/tBacg7f5sCv4yjPref9HKVPm7zBYSJAKZa/vgGdvgaaPAW38++wtU+IPxD2SycvLm7QptfLycqqqqryOhGw2m2e60WLp/+k4Nze332hLhLfPn5bGQ9cuJy7KjL05nY66y9A9FkBhMUWjtLFQdSjCzFfTZlHdeXjK/B0GA6ndF0yKp0/+a/R0wFN5xm0iFTdN6NNlZGSQn59PWVkZeXl5kxIQ+k7x9WWz2bjgggsAY1TXV0FBAbW1tezevXtQ8BLha9UJ7np/r2BrTqe1OZ0os4nvXZVO4u5VfMsSTbvJRIvZxDVpqQCk7bubjdPipd7fMCRIiaBVWlrKjh07qKiooLKy0uu0nFtJSQlHjx4d9jmXL19OTk6Oz/srKys9gSg7O9uzPgZGqnpZWRmFhYVDtkWEpzOOs1BecDZX/+YV6ps76XY4ufHxKmqij/BwRwTXzZ5Fm9mEexGrzoxsST8CEqREUNu2bRu5ubkUFBRQU1Pj87z77rtvRNOCOTk5XoOUO5OvsrISMEZQfQOUzWYjLy+vX2KHEAOdMCuBim8YgerAkTacGj51zuDk7iPEaydtA1ZY3OnpEqR8kyAVTCZ4ysyzJtXT0XssMhYu+fmU2aTNHVQqKiooKiryGSAaGxvH/BolJSUUFRV5Xm/Lli2D1ql27NiBzWYjOTmZNWvW9LvPvXblPu5OqhDhaV5SHOU3nM01D73CmwebKbFvYHPkgzSYB29VD8aW9A6HHbNZ3o69kV4JZe5AtOduaPoEpuguotu2baOiooKSkhKuuOKKCX3ugoICysrKsFqtlJeX+1yncqutrfV5/ZZ7FDbZqfMi+KVMi+YP+Su4/pF97DywCnpgpn0nhyO9X0d1244L2XL5TmKi4v3c0uAn2X2hbtkG+NYbUGwz/p1iAQqM7Dr3NU2+sv3cJY2Gu/WtEFFWVkZZWRnZ2dnU1NQMGaDy8/PRWnu9uUdd7q+HWvMS4SMxJpJHrzuT7KWz2OlcxYeHr0J7SU8H+HP3YQq2X0BzW4OfWxn8ZCQlpoT8/HxKS0t91si74447Rpw44eaeOpRUcjFZYiLN/PorGRQ++TpPVUMnED3zeUyRNmbFz+bEjnZe0s0AVDvb+Gz5+TiUYnb8bDZmbJS1KiRIiSmkvLzcZz29ganiw+l7HdTChQuHPHc8611CRJhN/DjndCyxUTz0d7A3G1vSW4+z8KNrMtn5/Ff5cevbgHHRL30K04Jk/sl0n5gyrFbrqIORL33XlWw225A3IcbLZFLcefFSvrP2RM+x1z62kVv2T9Ze8FsspuhBj5HCtAaltZeyvWLEsrKy9L59+0Z07ttvv83SpUsnuUViKpHfifDz2D8+5M7/fQP3W+9cSywtabei8fJerDX//mDqJj0NRSlVpbXOGu48GUkJIYQffWXFfH5xVTqRZiPT71NbB9pHYVqA5+Jjp1RJs4kmQUoIIfzs4mVzePCa5cRGGtdOdRxa6z3zTynuSE3hd4kJxvWOe+72c0sDT4KUEEIEwGdPnMnjeWcxPTYSe3M6nX0K06ba7aTa7Z5zS2YkcXNqCmsTHGG3Lb0EKSGECJCM45PYUXA2qQnR2JvTaX1/E+3vbGbnEc1Tn9aR3tnpOff/xcdRFxkRdtvSS5ASQogAWjI7gSe/sZIFM+IAcDg1dzStJ45oSusbWN3W7vVx4ZL9J0FKCCEC7LjkOMpvWMnStEQAdjpX8e3O6+iOnMUDh4+CjyzscNiWXoKUEEIEgZkJ0TyRv4LlC5IAI1Cd3vwAPzj9JdLi07w/JjbFn00MCAlSQggRJKbHRvLodWdx/kmpnmO//+dHJHZ+kRhzzKDzu9sa+LBuvz+b6HcSpIQQIojERpkpvTqT9elzPcf2vbmQ1O6vMDvONaJyTf/ZTPDVP32Vtw5UBqKpfiFBSgghgkyk2cRPck/n2pULPMfefHcx8Yd/wF9zqvjlgsuIcToBOGaCq1+8lc/+YVVIpqdLkBJCiCBkMil+cMnJ3Lamt97fqx/b2FD6Miemf5dtp9xIgitQdZsUx7qbQjI9XYKUEEIEKaUUt1xwAvd88RSUa7/E9w63cvmv9jLdei2PnHkXJi+Zf52OTra+WBQSZZQkSAkhRJC7+uwF/OyKM4gw9db7y/31XrosF6GV991+602ERL2/kA1SSql8pVR2oNshhBAT4YtnzOXBa7KIiTTeto+0dnNV2T9Iip7l9fw4rXGGQL2/kAxSSikLUArkDnFOuVJKD3HL91+LRUVFhWeLd1+777plZmailJq0vZ5KSkrIzMyctPOFGKvVS1J5/PqzSIwx9qtt6bJz+IPziPSyH1WbycT3U2bQ0/SJv5s5oUJ1Z95tIzgnw/Wvr3e6YxPUFjFKeXl5VFVVBeS1KysrKSoqmrTzhRivzPnJbC84m68+9AoNLV102M7AoTXzU5/gkBmitabTZIw/nk2I593oaJrKsznUfnhKbksfMiMppVSOUqpUKdUI5IzgIVagVmud5ONWMclNFj5UV1dTUlLi99e12Wzk5vocfI/7fCEmytK0RJ68YSXHJxv1/rqb0llacwlVHzXwjw8/4fLmVs+570RFUN9+aMpm/oVMkALuAPIB37uHubimAwFqhzxR+F12trGMWFRU1G+Ld3/Izc0d1RTiaM8XYiIdPyOOihvO5qTZCUBvvb+26NncdbSRvObQKEwbMkFKa52ptVZaa8UQa1Eu7i2Lh178EH5ntVrZsmULAAUFBX573ZKSEiorKyktLcVqtU74+UJMhtTEGLbnn03W/N56f8uaHuCuM17iphvf9fm4qVSYNmSC1Ch53lVcCRSNrmSJKqVU4XAPdmUO7lNK7WtoaJjclo7TrtpdrK1YO6WuRC8sLMRqtVJZWUlFxeTPulZXV1NUVER2djb5+cPny4z2fCEm0/S4SH739bNYvWSm59hj//iIjdtf6y2jNMCM6GEnnIJGuAapRa5/CzHWr2pdtwxgi1JqyFV7rXWZ1jpLa501c+bMoU4NqF21uyjeW0xdW92Um48uLy8HjCSKyZ5Sy83NxWKxeF5zos8XYrLFRpnZ9tUsvnjGHM+x516vI7p1ndfCtC0dx6j6z9P+bOKYhWp233DcmX0VWmvP1KBSygqUAxlKqS1aa7+mbZ32yGmT/hqdjk42vbSJTS9tmtDn/fc1/57Q58vIyCA/P5+ysjLy8vImLSAUFBRQW1vL7t27sViG/3Q52vOF8JdIs4mfbjgDS2wkj7z8IQBvvLOYxQuvYnry/3GovR6lNVopukyKgn/cyU86j/LZM64PcMuHFpZBSmu9xsfxWqVULlCDMcqS3OIAKi0tZceOHVRUVFBZWelJqvCmpKSEo0ePDvucy5cvJyfHSP6sqKigrKyMwsLCIZ/bbbTnC+FvJpOi+AunkBQfxc8q3wPg/QNLWNKeReXXz8T2QQUF++7jqNlEl1J889WfMf3NbTTZO4I2PT0ogpRSqmaUD7FprSfl6klXoKoFrEopq9ZaMgADaNu2beTm5lJQUEBNje9fk/vuu29E04I5OTnk5ORgs9nIy8vrl6gxlNGeL0SgKKW4NftELLGRFD/7FgDvHGoh59d7+d116/nduTPI++t3+DTChFYKm93IAqxrq6P4r0Xw0T9Yt/qeQH4L/QRFkNJaLxr+rIDw6wW9Ez1l5l6T6nR0eo7FmGMoXlkcdJ+WfHEHlYqKCoqKinwGiMbGxlE9744dO7DZbCQnJ7NmTf+BtTv13X3cnQ4/mvNlpCUC7dpzFmKJi+I75a9hd2o+PtZBzq9f5tHrzuF3ax9i7Z6vYx9Q969TKbbWPMm65NNg2YYAtby/oAhS/uS6RqoRYzSW5ON+q+v+KX0RjDsQba3eSn1bfdAO54ezbds2KioqKCkp4YorrpjQ566trfV5PVZlpbGRXN9U+NGeL0QgXZo+l8TYCL7xWDVddidHWru4ouxlHrp2OQ5lAgZXUK83m4x6f0ESpMIuu88VeCoBi1KqtO99rgC1x/VlSKxHrbOu44WcF3j9mtd5IeeFKRegACwWC6Wlxo/KV7ZfUlKSp/bfUDd3hYj8/Hy01l5v7uue3F/n5OSM+nwhgsX5J83isevPIsFd76/TztW/+SeWKO+ZydFaB1W9v7AbSbnkAlVAvlJqA7APSKZ/1l9ZoBonBsvPz6e0tNRn8dk77rhjxIkTQoSb5QuS2Z5v1Ps70tpFZ4+TQx+cR+LsJ+g0DZjyM5m4OW0WP+1sIi5meoBa3Cssg5RrNLXIdeHuGnorUFQCpVK3LziVl5ezaJH35cvCwmGvwRYirJ08J5GKG87mK7/5J580dtDVlM5C3sGWWkW92USs1rS7CtPujYkiZ/t5dMcmcbijIaBLBSE53ae1rnCVSBpycUBrXaK1XtOnqOwaCVDBy2q1SjASYhwWpMTz5DdWsmSWUe/vlaYrOanmi/ytroeXP/yUG5t66/19TA+HOg4HvBCA0l62HhYjl5WVpfft2zeic99++22WLl06yS0SU4n8TohAsLV3c91v/0X1R73ru9ecPZ8fXHwy25/9KvfaXgUvO/6mxafxQs4LE9IGpVSV1jpruPNCciQlhBDCN0tcFI9dfxafObE3eeKRlz/k1h2vkXPJo14DFASmMK0EKSGECENxURE8+NUsLl7WW4R252sHyXt0n8/CtMlRif5qnocEKSGECFNRESa2XpnOV1Yc7zn24jsN0Pg5or0UprV12Xhx/4P+bKIEKSGECGdmk+KeL57KLecv9hx7r/ZEYpuuIDV2NgowuXIXHEpx62s/Y+fezX5rnwQpIYQIc0opblu7hLsuPtlz7OOPl9LzwXd57uJ/8tw5P2ae3QkYgep77z3OOQ+fyrLfnsrah05l14t3TlrbJEgJIYQA4LpVC3lgw+mYXRf4fnSsnct/vZe2hJU8uvYhTrD3ZoM3mxRaKerMiuIDT09aoJIgJYQQwuOyjHmUfiWT6AgjPDS0dLHh1y/zkfMEHr5kB5HOwZctdZoUW2snZxNFCVJ+JtelCTf5XRDBKvvkWTx63ZkkRBtFiZo77Xz5wX9S3ZiC3Xt2OvWTFE0kSPmR2WzG4XAEuhkiSDgcDsxmc6CbIYRXZ1ln8If8FaRMiwKgs8dJ3iP7mOXjLWy2c3LaIUHKj+Li4mhtbQ10M0SQaG1tJS4uLtDNEMKnU+dOp/yGlcy1xAJgd2qSG9KJGTDlF+PUbLSun5Q2SJDyo8TERI4dOyajKYHD4eDYsWMkJvr/4kghRmOhq97fCanTAPin7UoW1qeT2uNEaU1qj5Orp62dtN18pXbfOI2mdp/WmsOHD9PW1kZycjLTpk3DbDajfJQgEaFFa43D4aC1tZVjx44RHx9Pamqq/PzFlNDY1s3XfvsvXv148H5usZFm7rvsNC5Nnzvi5xtp7T4JUuM0miAFxhtVS0sLzc3NtLe3y6gqzJjNZuLi4khMTCQhIUEClJhS2rrsZNyzmy774AWouZZY/r7p/BE/10iDVFjuJxVISikSExNlmkcIMeXER0fQ7SVAARy0dUzKa8qalBBCiBGb40qiGOnx8ZIgJYQQYsRuv3AJsZH9L52IjTRz+4VLJuX1ZLpPCCHEiLmTI+5//h0O2jqYY4nl9guXjCppYjQkSAkhhBiVS9PnTlpQGkim+4QQQgQtCVJCCCGClgQpIYQQQUuClBBCiKAlQUoIIUTQkrJI46SUagA+HOKUFOCIn5oTSqTfxkb6bWyk38ZmPP02X2s9c7iTJEhNMqXUvpHUpxL9Sb+NjfTb2Ei/jY0/+k2m+4QQQgQtCVJCCCGClgSpyVcW6AZMUdJvYyP9NjbSb2Mz6f0ma1JCCCGCloykhBBCBC0JUkIIIYKWBCkxqZRS+Uqp7EC3QwgxMkqpQqVUVaDb4SZBahIopTKUUuVKqUbXrVwplRHodvmbUsoClAK5w5w3qv4K9f5VSuUopapc35t2/b9wiPOl/+h9c3X1WaPr/zlDnC/9NoDrA+UWYML6Ydz9prWW2wTegGxAu241rpv765xAt8/PfVHu+r5LJ6q/Qr1/+/SZ+/ur6vu19J/Pfuv7fVQN6LdBv3/Sb1770AI0ur8vH+f4vd8C3jGhdHP9kN0/gGwfPyhLoNs5yX2QgzF6amSIN4mx9Feo92+f76MRyBjwfe8e2JfSf572F7ravnvAcWuf38OB/Rn2/ealH3f3+X60l/sD0m8B75hQuvX5Y9ni5b5S132FgW7nJPdB30+wwwWpUfVXqPcvvaOofC/3ef7gpf98/s5Zh/gdG3M/hGq/+einfFyjnSHO8Wu/yZrUxFrj+ne7l/tKB5wTkrTWmVprpbVWDLMWxej7K9T71z1PXznwDq21DagFz1ofSP+5HQMqtda1Q5wzo8//pd/6cK0PbcHow6Euzg1Iv0UMd4IYFSuA1rp64B1a62qllOccAYy+v0K9f3MBvL3ZugKT+/u3uQ5L/wFaa69vdK4+K3B92feNUvqtv3LAxvAfKgPSbxKkJtZIflGn8i/zRBttf4V0/3r7YwbPm+0e15clfe6S/hvA1VflGN+HFePNt2BA30q/uSilSjHavqbPhx9fAtJvMt038Yb6QQ/3SxCORttfYdW/rpTgKoypwEqtddGAU6T/+kvGeONL7nPM4uW8sO83V3p+PlCitR40xeyD3/tNgtTE8/YHMZL7wtVo+yss+lcpZVVK7cbIuLJivJF4m9aS/utDa12rtV6ktU4CFmGs421RSm0ZcGpY95trxLkNqPXywWcofu83CVITa6iF29GcEy5G219h0b+uC3drMFJ1K4BFPt5IpP+G4Frbc6+z5Pe5S/oNNuAKEkqp3X1vuKbg+hxzV4wJSL9JkJpY7uyrQVdT9zk2JacGJslo+yvk+9e1RrAF43vN1FrnDpG1Fvb95x5xehkpAf2SUPp+ag/7fuvDivFhqO/Nzf21u+8C0m8SpCZWuevfK7zc5z5W6uW+cDXa/grp/lVK5WN84q90TVl5TaToQ/rPSD/Ppv9IyUMp5V6Y7xvow77ftNZl7ktFBt5w9VWfYxWuhwWm3wJ9EVmo3ei9gLXvFe4Z7uOBbp+f+yKHIS7mHUt/hXL/0lsyZsTVC6T/+lVK2DLguIXeC30LB9wX9v02RH96vZg3UP0W8A4JtRvGJzr3D8a98O3+elAlgVC+jTBIjaq/QrV/6V9CpnGom/TfoO+pb/mjRgbX7tvt5TFh329D9OdQQcrv/RbwDgnFG8b0Q98/kir61K4Kl9tIgtRY+isU+5c+ny6Hu0n/ef2e3BX33aPRRtcbos8iptJvPr9Pn0EqEP0m28cLIYQIWpI4IYQQImhJkBJCCBG0JEgJIYQIWhKkhBBCBC0JUkIIIYKWBCkhhBBBS4KUEEKIoCVBSgghRNCSICWEECJoSZASQggRtCRICRECXPsqaaVUo2vXVW/31QSqfUKMlQQpIUKANjb3K8IotDpwA0D3vj65CDHFSJASIkRorUuAaiDfvfOpayPFDKBED7+JohBBR6qgCxFCXDvR1mDsrrrG/X+t9aKANkyIMZKRlBAhpM+0nxVj3x6QaT4xhclISogQ5EqSsAJlWuuCQLdHiLGSkZQQIcaV3Zfs+jIrkG0RYrwkSAkRerZhZPlVAxlKqcIAt0eIMZPpPiFCiFIqByPlvALIAw5gBKxFrvUqIaYUCVJChAjXNN8B15cLtda2PkGrUmu9JnCtE2JsZLpPiNDhnuYr0lrbALTWFUAlkO0KWEJMKTKSEiIEKKWygd1AtdY6c8B97munbLhGWAFoohBjIkFKCCFE0JLpPiGEEEFLgpQQQoigJUFKCCFE0JIgJYQQImhJkBJCCBG0JEgJIYQIWhKkhBBCBC0JUkIIIYKWBCkhhBBBS4KUEEKIoPX/AY3LKonhGxOlAAAAAElFTkSuQmCC\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"for N in range(4,45,20):\n",
" w,K,F=int_N(N,10)\n",
"\n",
" plt.plot(linspace(400/N,400,N),w,'o-', label='N=%i'%N)\n",
"plt.xlabel('x')\n",
"plt.ylabel('w(x)')\n",
"plt.legend()\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Demonstrate quantitative convergence of tip displacement under given load.**\n",
"\n",
"The relative error is calculated based upon finite-difference calculations\n",
"\n",
"$e_{rel}=\\frac{\\delta_{new}-\\delta_{old}}{\\delta_{new}}$\n",
"\n",
"The absolute error is compared to the analytical solution of $EIw''''=0$ for a cantilever beam\n",
"\n",
"$\\delta_{an}=\\frac{FL^3}{3EI}$\n",
"\n",
"$e_{abs}=\\frac{\\delta_{new}-\\delta_{an}}{\\delta_{an}}$"
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"<matplotlib.legend.Legend at 0x7f28ebf3c908>"
]
},
"execution_count": 21,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAZ0AAAEeCAYAAABR+8jUAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3X2MG+l9H/Dvb1cv1uokcVeSbd3FuruhYseFGzvcPacpGrg+cZ0/EqBJTEpFgCIFYu+6fxQokGB5Mpo2QIsopNMABYrUS6VAkBZIJNIJmqStHVIuEiQoUu/y4hT2OU53pDvfuYnvtMvTne4k7cuvf8wz5HA45JDc4fv3AwyWnHk48/DR6vnt8zYjqgoiIqJBmBl2BoiIaHow6BAR0cAw6BAR0cAw6BAR0cAw6BAR0cAw6BAR0cAw6BAR0cAw6BAR0cAw6BAR0cAcGXYGRs25c+f0mWeeGXY2iIjGyubm5huqej4sHYOOzzPPPIONjY1hZ4OIaKyIyMudpGP3GhERDQyDDhERDczEBR0RSYhIQUR2zFYQkcSw80VERBMWdEQkCWATQArAttlSADZFJDXMvBER0QQFHRGJASiZt8uqGlfVOIBls69g0hAR0ZBMTNABsGJ+5lS17O40r/O+NERENASTNGXabdHcDDi2DifgLAPIDSxHPdrdP8DjvQPs7Sv2Dg6wd6DOtu+83jebKqBwfrrcffXXgKq7ByZt/TPqpvOk6TtVQA8gODAZavwCAgVUoeaYtEpj0rk/xZcGqJ8LtfSop/PsF18aBJaGttjvnisofTPp6mm9AWlbfT5gv0SQ36D90vIrdH6Olt8j6HpdfY9RFH0+u/s96tzTH/l7OHVmoS/nBiYr6FgAoKoV/wFVrYhILc0wvft4H3/5ahVff7WKl++9g//35kP8zZsPcf/hLh482sODR/t4vH/Q07mPYRcxvI2YvI2TeIgT8ghzeIQ5PMSceX0CjzAnD3EMeziKPRzDHo5gH0el/v6o2XfM7HPeH0CgmMUBZsw2iwPMiJr3zrHGNFpL5+4/Ir19NyIajG+f+H2cSnyib+efuKDTSxoRWYHpert48WKUear55nfv44t/vIWvfONv8Giv+4r3CbyDD8qr+IB8D0/JPTwpb+BJuYfzUsW8vI0Y3sYT8rAPOSciis4kBR0AqIYcC5xIoKp5mHGfpaWlSNus+weKXyv9Fb74xzb2Dzo79dEZxQ8eeQV/f+YbWMS38P14BU/he1Fma+icto+zofbT7YQwx8Tz2hx30njTm9ci2DtyEm994Hm8e/F57B87bVJ7z9uONO0J/ExzsiEZjYyMRufWaJTFZFDMHo3he9/7HhYWFnDkSPQhYtKCTrvZaQOfuaaq+Pzv/h/c3PhOw37r/EksPT2PH3j/aTwZO4ELZ96D+RNHcWb7RZz8VhGzL/1XyLvb3V9w5ghwYh44sQAcfwI4OgccO2l+zgFHTzrvj80Bs8eB2WPA7FHz8xgwe8Tz+igwc9Tz+ggwMwvIDCDm54z3ddCxVscFMxJtRXFwcICXX34Zx48fx/vPnsWxY8cgEV+DaNKpKh4/fox79+7hO9/5Dp5++mnMzEQ732ySgo6N8C42exAZcX2p8lpDwPnhZxfwiz/xd/CRp87UE6kCW18FvvIrwKv/u/XJZo4A5z4InPt+4MwHzPZ9wKkLwNyCsx0/DUxpRbuzs4MjR47gwoULDDZEPRIRHD9+HBcuXMCrr76KnZ0dnD17NtJrTFzQEZGEfzKB544E7brfIvXu433kvvyt2vuf/qGn8Kvpj2JmxlMhvvU3wB/8C+Db/6P5BCffC1ifAJ75UeCpRSfgHDk2gJyPp7fffhsLCwsMOEQREBHEYjEGnRAFAEkAVwH4Z7BdNT/XB5WZ3//6a/jeW48AAO87fRz/9qc+0hhwXt0Efvsq8OD1+r7ZY8APXgE++jPAxR9xuqioIw8fPsTc3Nyws0E0Mebm5vDd73438vNOTNBR1byIrANYE5GbbmvHtHLW3DSDys/vvfha7fVn/oGFuWOeor77p8B/SQF779b3JX4W+MSa02VGXTs4OIi875loms3MzODgIPolDhMTdIxVOK2ZTRFx70qQ9BwbiJ0Hj/Hnd5yJADMC/KOPPVk/+Pq3gd/5mXrAObEApH/T6UqjQ2HXGlF0+vX/aaKCjmnt2ACyqAebCoCM99Y4/bb58k5tofTf/b4Y3nv6Pc6bvcfArX8CPHzTef/E+4B/+t+Bc5cGlTUioqGaqKAD1O61tjjMPHzt5fp0548/M18/8Gf/HnjdTC44Ogf8zE0GHCKaKuwE74Ovf6c+SW7xaXMPo+p3gD/5Qj3R5X8FPPlDA84ZEdFwMej0gf36g9rrH3j/KefFn38R2Hdms+HCx4CP84bXNF4ymQxEBOXy4XqqbdtGtdq4eiGqc9PoY9CJ2INHe7Wp0kdmBN83fwJ49BZQ+c/1RJ/8vLNCn2jK2LaNeDyOz372s8POCg0Jg07E7t6rt3I+sDCHI7MzwF/8NvDITB44ewm4tNzi00TTaXV1FaVSCUtLS8POCvUZg07E7r7xTu31s+dOOi/+0vOInx/+HBd90lD5u7ZGgWVZSCaTiMXG/+G+7co3yrIfxX/HTrD2i9hr1XrQubgw50wgeG3D2TFzBPjIp4eUM5pGuVwOIgLbtlGpVBCPxzE/P9+QplKpIJ1OIx6PQ0SwuLiIfL67ddS5XA6Li4uYn5+HiCAejyOTyTRUjMvLy4jH4wCAYrEIEUEmkwEAlMtl5HI52LZze8R0Og0RaZmP1dXVwONRfJduztOufDsp+2KxiOXlZczPz2N+fh7Ly8soFotdXWfcMOhE7N6Dx7XX508dB/76j+oHn/lR58acRANm2zYuX74MAEgmk7X9brAoFouIxWJIJBKoVCpYXV3F8nJn3cCLi4vIZDKwbRtLS0tIJpOwbRu5XK52TcAJFCsrzgQay7KQzWZr1yiVSrVzuGkBoFAoBF7z1q1bAIArV65E+l16PU+r8m13LJ1OI51Oo1wuw7IsWJaFcrmMdDpdC8bdXGdcTNw6nWG793Y96Jw9eQy488f1gx/8sSHkiJ554b8NOws9u/srPx7JedLpNK5du4a1tbXaPtu2kclkYFkWSqUSLMu5SXu1WsXly5drrQ/vZ/zy+TwqlQqSySRKpVLDsXg8jkqlAtu2YVkWUqkUEokE8vk8EolE2/O6XW3lchnVarWh283d5+2Oi+K7HOY8QeXb7lixWESxWGy6jm3bWFxcRC6Xw9WrV5FIJELPNW7Y0onYvbcf1V4vzB0B7vxJ/aD1DweeHyLAaVn4Kypva8Kt9AAgFovVWhjXr19ve96FhQWkUilks9mmY6lUCgBqrZduua0Y/zRqN29u/r2vD/NdDnOeoPJtd8xtyfivY1kWbty40ZCm0+uMCwadiG17uteePPgu8O6O82buLHD+B4aUK5p2QV0xGxsbte4jP7e7p1qtth2wTqVSKBQKTeeoVCqHXnOTTqcBADdv3mzYf+vWLcRisVpQA6L5Loc5T7uurqBjtm23vI77vTY2Njo617hh91rE3vB0r73/rW/WDzy1NLUPWBu2qLqoxtlzzz3XtM+tOMNu7Li9vd12Vlm1WsWtW7ewubmJjY0NVCr+J4v0xu0+8w6su11r7tiQNw9ANN+ll/MElW+rY27Lr930cMuyAluI7a4zLhh0InbvQb177fT2X9YPPNX8Fw3RsLiVaywWw7Vr19qmXVhoPfmlXC7XBtfdac+rq6tYWlrCzZs3kcvlDpXPK1euIJ/Po1wuI5lMBnatRfVdojpPlPzjWZOAQSdC7zzew8Nd5/kTx47M4Ogb36gf5H3WaITEYrFaZXaYMQLvGIi3uwto7hbrRTqdRj6fR6FQQDKZxK1bt2BZVkO3VFTfJarzhHHHcIK6z1xu99ukBRyAYzqR2nlnt/Z6Ye4Y5PW/qh9874eHkCOi1paWllCtVgO7w6rVam3dSDvemWlBxw7L7WK7detWrWvN28pxRfFdojxPGHdsKOg6bnfipN6dgUEnQu882qu9fvLo28A795w3R08Cp/lEUBot7oyzy5cvNwWIdDqNarVaG8xvJRaLwbbtps/n8/la5Rk0eN/NavorV640BJugABfFd4nyPJ1exz2ny7bt2n3pWq3VGXcMOhF65/F+7fWHZj3PFj//Qd76hkaOu1amWq0iHo8jHo/XVse7Yyj+AXs/93g8Hq8tdnTvRuAGh0wmU1vN746FuIsgg1bf+7mVvG3bSCQSDVOMo/wuUZ4nTCqVQiqVgm3bmJ+fx+LiIhYXFxGPx2sTJSZhploQ1oQRevC43tJ5VjxB59yHhpAbonDZbBalUgnJZBLb29vY2NiAZVlYX19vWuzZ6vPr6+uwLAvFYrG2UPTOnTu1cRjbtrG5uQnAaRmtra3VZqV10gXnXQQa1LUW1XeJ+jxhCoUC1tfXa2Vk23ZtssT6+npk1xk1ou5zlQkAsLS0pO0G+Nopf/Nv8Znfcj77H973B/iJN3/bOfCJF4BPtp8NQ4fz0ksv4cMf5rgZUZS6+X8lIpuqGjoQxZZOhN7ZrXevve/gb+sHYheHkBsiotHDoBMh70SC83sMOkREfgw6EXrgmUhwlkGHiKgJg06E3jUTCY5hF6d233B2yixw+qkh5oqIaHQw6ETIbemch2cNwhPvBWZ54wciIoBBJ1LvukFH3qzvfOK9Q8oNEdHoYdCJ0AMzkeBcQ9B535ByQ0Q0ehh0IvROraXj6V47yZYOEZGLQSdC75iJBOfB7jUioiAMOhFyJxKc45gOEVEgBp0IPdwN6F5j0CEiqmHQidDjPecBbgvyVn3nyfNDyg0R0ehh0InQ7r4TdM7gQX3neybvyX9ERL1i0InQ3oFzx+4z4gk6Jxh0iIhcDDoR2t1jS4fGXzweh4gMOxsAnAfAiQjK5fKhzmPbdldPK6X+YdCJ0O6B4ij2MCePnB0yCxw/NdxMEU0527YRj8drj4Gm4WLQidDu/oGvlXMGGJG/GImIRgGDToT29hVn5O36Do7nEFGH2nX/Rdk1OOxuRgadCD1uaukw6NDoyOVyWFxcxPz8PEQE8XgcmUymbSWUyWRqYzyLi4vIZDKB6SqVCtLpdC1tPB7H6upqy3MXi0UsLy9jfn4e8/PzWF5eRrFY7Oh7rK6uQkRQqVSajpXLZYhILZ/Ly8uIx+O1a3qPtcr74uIi8vl8R3nx6/RcuVwOIgLbtlGpVBCPxzE/Px96zNVp+XVyroFTVW6ebXFxUXv17At/qD977d+o/uvTzvZbP9nzuag73/zmN4edhZGWSCQUgMZiMU0mk5pMJhWAAtBEItGQ1rIsBVBLE4vFap8HoJZl6c7OTi19qVRqOFcymdRYLFZL65dKpRrSe8+9trbWkHZtbU0BaKlUqu1bWVlRALq5udl0bjcv7nkKhUItvWVZms1mG86VzWZb5iWZTHZVxt2cy01bKpU0FoupZVm1NO2OdVt+YecK083/KwAb2kEdywe9RGT/QHGgvunSbOmMhl86M+wc9O6X3gxPEyKfz6NSqSCZTKJUKjUci8fjqFQqsG0blmU1HCuXyygUCkilUgCcbpnLly+jUqkgk8lgfX0dAGoth1KphGQyWfv88vIyyuUyisVi7RzFYhHFYhGWZaFUKtWuads2FhcXkcvlcPXqVSQSiUN/bwBIpVJIJBLI5/NIJBJYW1urHbNtG5lMpikv7vcsl8vI5XINn2ml13Ol02lcu3Yt8BpBx3otv3bXGTR2r0XEXRh6Wt6p73zPGFd2NDEWFhaQSqWQzWabjrnBwLbtpmNra2u14wAQi8Vw+/ZtAE4gc7vO3G4ub8ABgGw2i2w22xDM3ABVKBQa9luWhRs3bjSk6bfV1dXAvMRiMRQKBQDA9evX+3ouy7JaBoKgY72WX7vrDBqDTkTcoHMSD+s7OV2aRkAqlUKhUGj667dSqbRd/+JWpF6xWKwWXNxA5Z7Xbdm43JaF97q2bSMWiwW2ZNwAt7Gx0elXO5SNjY2WebEsC5ZloVqtdjTw3uu5/IE67Fiv5dfuOoPG7rWI7O47dyOYE0/QOXZySLmhBhF0UY27arWKW7duYXNzExsbG4GD8H7+7jZXIpFAuVyGbdtIJBIoFAq1gFMulxGLxbC0tIR0Oo0rV64gFnO6md0gtbS01PaaQa2ufnADQNhC2O3t7dp3iPpczz33XMu0/mOHKb921xk0Bp2I7JmWzhwe1Xcy6NAIKJfLWF5eBuBUSslkEqurq1haWsLNmzeRy+UOdX7LsrC1tVUbA3KDT7lcRiaTwe3bt7seo6lWq6EV/WG4QSIWi+HatWtt0y4sLAzsXFHpd/kdBoNORB7Xgg5bOjRavOMN3jEaALh582bLzwVNLgDqYzj+QJJMJhu63rLZLPL5PNLpNLa2tmrnatd95nYf9VphdtpK8l7jsGMdUZ6rnUGU3yBwTCcie6Z77aR4WjpHGXRo+Nzg4Q847rFW3NlpXtVqtdaF5nblxONxpNPphnSWZdU+772GO7YR1L3nrjNp133ktb293bTPHbTvxNLSUsu8VKvV2vqXQZ+rnSjLb1gYdCKyy+41GlGxWAy2bTcFmHw+X6uoggbLc7lcw4JDd/ovgFo3kht4isVi06QE97PeFpE7gy6dTjdc07bt2r3RwmavuYs9/UExn8+3nRjh/45uXi5fvtxUNm7+/MG0lSjP1cl1DlN+Q9fJYp5p2npdHPqN197UpzN/qH/6L3+kvjj0/97u6VzUPS4Obc1dYAlAU6mUplIptSxLY7FYbaGhZVm6vr6uqs7iUMuyaosOva/RZqGjmzaZTNYWmCJgEWe7xY0rKyuBefcu6Nza2mq4XiqVqp3DzYt3keTOzk7D9y8UCoFl4+bdXdja7eLQbs7l5tObl06OdVt+YecK04/FoUOv5Edt6zXo/MUrO/p05g+18ouJetB55c97Ohd1j0GnvfX19VogsCxLV1ZWancVcO884FZY6+vrtQC0trZWqzgTiYRms9nA8xcKBU0kEg13IlhZWdGtra2W+XErZfcuCUEVY1DQUVXd3NxsqNQTiUQtjf+uA/7v4f8OpVKpIS+JRKL2/bvV6bkOE3RUOy+/UQw64qSdPCKyAsBW1a4exLG0tKS9rBPYfHkbn/6P/wtfObaGD8286uz83J8B7/9I1+ei7r300kv48Ic/POxsEE2Ubv5ficimqoYOKPU8e01EngdgAYgDiAHYAlCFE+3+otfzRkFEYgDWAeQBHO7pTx16vGfW6XBMh4iopa6Cjoh8BsAqAHdkMGgllIpIFUAJwLqq/s/DZbEnNwZ9wb0DM5GAi0OJiFoKDToichrANQBrcIKMDeBLcILKtnkPAAtwWjzPwQlKVwCkRcQGsKaqvxd57hvzmQKwbK478EnqnL1GRBSubdARkZ8G8Bvm7RfgtFzuhJzzS57PJ+G0jL4kIn8E4Iqq3j9Eftu5hnoLbOB29xUzOMAJeVzfeeTEsLJDRDSSwtbp/AaAjKouqOoLHQScBqpaVtU0nFbQXTiBoS9UdVFVRVUFQFcT4kVkRUQ2RGTj9ddf7+n6u/sHOAHfwtAZLoMiIvJq29JR1UhuFKSqVQCfi+Jc/aCqeTiTDrC0tNTTdL69ffXdAmcukrwREU2SSP4UF5HTInJdRL4iIr9sxoGmyuP9AxyX3foOdq0RETU5dNARkTMA7gDIADgL4AUA9rQFnt39AxyHN+gcH15mptSkrjkjGoZ+/X+KoqVzA85ssbhZGHQJzhhOPoJzj429fWXQGaLZ2Vns7u6GJySijuzu7mJ2djby80bxaIMkgLw7yUBVbRG5AaD5lrYtiMhWl9esqupil5/pm+0Hj/Hv/uivEGfQGZpTp07h/v37OHfu3LCzQjQR7t+/j1Onon/6cdiU6efh3Ermbsh55n3vu2qXqWq8m/SjZvvBY9x/uIdjslffOcugM0gLCwt45ZVXAACnT5/G0aNHQ5/iSESNVBW7u7u4f/8+dnZ2cPHixcivEdbSiQMoi8gXAbzQYo1NAcBnROS6qn5dRH4IwAqc29BMleMNa3QYdAbp+PHjuHjxIra3t3H37l3s7+8PO0tEY2l2dhanTp3CxYsXcfx49PVY2JTpGyJyB8AXAayKSFZVP+9Ltgani61i7j5gwblLwYg/1CE683NH8c+fv4T49mvAt8xOBp2BO378OC5cuIALFy4MOytE1ELoRAKzwPMSgH8G4HMick9Efs5z/E3TPXYNwItwWkSX+njngZFz9onj+PlPfQg/+RHPeAKDDhFRk45nr6lq3iwWzQK4ISJ/LSKf9BzPqeoVVf1CPzI6FvY8dyTgmA4RUZOup0yrag7OlOivArgtIl8Wkacjz9khqGrR3BJndaAX3vcEHbZ0iIia9LROR1WrpkK/BGAWzmLQX5+2BaFN9hh0iIja6SjoiMjzpjttX0S+7XarqaqtqssAfgzAxwHsiMgv9zG/o60h6LxnePkgIhpRoUFHRC7DeXbOm3BucfMWnGnU3vGcsrkbwVUETDaYGg1jOseGlw8iohHVSUsnC2eB6JKqfsHcCeCu2d/AjKW4kw2m6jY4AHxjOmzpEBH5dRJ0EgDKvn0lAC1vQ+OZbDBd9jyPNjjClg4RkV8nQcdd8OnlLgBtSVXf7DVTY2vPe0cCtnSIiPw6CTp5AMvmOTkfE5HrcO5AMHW3uQnlbelwTIeIqEnoXaZVNScicTiTCDIABM5dpX+135kbO/ts6RARtdPRow1UdVVEfgXO+E7FfYwB+TSM6XCdDhGRX9ijDU6791AzgeZQwcZ7vom0x7tMExG1EzamUzVjOYe604BZXPo1ODcFnVwNYzoMOkREfmFB5xKAT8G508DveBeEhhGRZ0TkF0Tkr+FMsS6r6mQHnX22dIiI2gl7no4NYElEVuA8N+eKiCiACoANAFsAqiZ5DMBZONOpk+a9ACgC+NRUjAM1jOlwIgERkV+nEwnyAPIikgSQBnAZQLs7OJfhtG7yU7Vex9vSme2oaImIpkpXNaOqlmHuTiAiZ+C0aiw4dx+wAWyr6otRZ3JsHHgekTxzdHj5ICIaUT3/OW5aMC+ajQBgf7f+epZBh4jIr6fn6VALB3v11zPsXiMi8mPQidKBp6XDoENE1IRBJ0reMR12rxERNWHQidI+WzpERO0w6ESpYUyHLR0iIj8GnSg1BJ3Z4eWDiGhEMehEyRt0OKZDRNSk46AjImdEZF9EvtzPDI01jukQEbXVcdAxi0HvAoj3LTfjTNU3ZZotHSIiv26719IAzorIz/cjM2NNDzxvBJhhzyURkV+3fUDbcAJPQUT+MZz7sH2tVWJV/d1D5G288BY4REShug06NgCF88iCRbNpQDox+6dnChdvgUNEFKrb2vEFBAcZ4ngOEVGobh9tkOtXRsZew2MNpqeBR0TUjUhGu0XkdBTnGWsc0yEiCtVT0BGRj4nIV0TknojsA9gxr78sIp+MOI/jgWM6REShug46InIdwCaAZQDzAO6YbR7ApwCUReTXo8zkWOBjDYiIQnUVdETk0wAyAN4EkFbVGVW9ZLYZAFcA3AewKiI/FX12R1jDmA6DDhFRkG5bOqtwZq8lVPVL/oOqWoQzjVpM2unBMR0iolDdBp0lABVVvdsqgaraACoAnjtEvsYPH2tARBSK92qJSsOYDqdMExEF6TbobAJIiMjTrRKIyBkACQAbh8nY2OGjqomIQnUbdNbhjNeUROSj/oMi8jE4wUYBFA6fvTHCxxoQEYXq9o4ERRH5EoBPA6iISBX1Fs0SgBicoFRQ1d+INKejjut0iIhCdT2mo6ppOFOj78JZm7OMxjU7aVW9GmEexwPX6RARheqpdjRTo4sAICLPmn13IszX+OGYDhFRqEM/rlpV70x9wAE4pkNE1AE+rjoqHNMhIgrFx1VHhUGHiCgUH1cdFd4Gh4goFB9XHRW2dIiIQk3U46pFJAXgGgALzpqhCoCbA3niKadMExGFmpjHVYtIAUDKvLXNloBz255VVe3vBAg+2oCIKNShp0yPAhFJwgk4VQCLqhpX1UU4C1bLACwRWe9rJjimQ0QUalKmTLvP7smoasXdqapVOBMfAGClrzngmA4RUahJmTKdMD/L/gMm8NgAICKxvuWAjzYgIgo1KVOm0+Z6tv+ACTSWOV4N+rCIrMC0hC5evNhbDg4O6q/Z0iEiCjQRU6a9XWoNmXACzm3ztuUkCFXNA8gDwNLSUm+z89QTdIQtHSKiIBM1ZdrLTC5Yh9PKKatqpq8XVM/sNeEDWYmIgkzMlGmXiFhwgk3S7Mr1PeAAvpYOgw4RUZBDDT6IyGkACwC2VfX+Ic6z1eVHqmZKtP88awCy5m0Rzmy2pnGevvAGnRkGHSKiID0FHRH5BTgr/93ZYBkAvyoiXwFwS1X/Uzfni2LhplmHswJn3Cndapynbw7YvUZEFKbr2tEElizqTwoVz+GPA8gPegGpZ/ZZ2SwMHWzAATiRgIioA10FHRH5LJxHU78IwFLVS74kzwL4KoBlEfm5aLLYEXfMJt02VT9xTIeIKFS33WurcGavPR80hmPWwSyLyDaAzwHoqputF951OADuiEjLtKo637eMNIzpsKVDRBSk26CTALDZwaSBDThreAbB8rzu3x0HwnBMh4goVC+LQzthhSeJhhm/ad28GVhGOKZDRBSm2z/JX4TzqICnWyUQkWfhBJ2Nw2Rs7DQsDh1+DCQiGkXdBp3rcFoVJRH5qGe/AoCIfAxOsFHU18tMB47pEBGF6iromK6sFwBcAlARkXtwAsznzetNOFOpc6r61agzO9I4pkNEFKrr2tHcCucSnKnRYrZ58/M2nIeoXYsyk2NBPbek45gOEVGgnu5IYG4ts+y+F5Ez5iFv04s3/CQiChVJ7Tj1AQfg4lAiog6wdowKJxIQEYVi0IkKJxIQEYVi7RgVdq8REYVi7RgVTiQgIgrF2jEq3inTHNMhIgrEoBMVjukQEYVi7RgV3vCTiCgUg05UOKZDRBSKtWNUGtbpsFiJiIKwdowKx3SIiEKxdowKx3SIiEIx6ESFi0OJiEKxdowK771GRBSKQScqHNMhIgo9EPaCAAAM6ElEQVTF2jEq7F4jIgrF2jEqnEhARBSKQScqDYtDZXj5ICIaYQw6UeFEAiKiUAw6UTngmA4RURjWjlHhmA4RUSgGnajwhp9ERKFYO0aFYzpERKEYdKLCxaFERKFYO0aFi0OJiEKxdowKgw4RUSjWjlHhmA4RUSgGnahwTIeIKBRrx6iwe42IKBRrx6hwcSgRUSgGnahwcSgRUSjWjlFpmEjAYiUiCsLaMSq84ScRUSjWjlHhmA4RUSgGnahwTIeIKBRrx6hwcSgRUSgGnahwcSgRUSjWjlHhmA4RUSgGnSioAtD6e5GhZYWIaJQx6ETBfwscBh0iokAMOlHgeA4RUUdYQ0aBN/skIurIRNWQIrImIpsioiKyY16n+n5hTiIgIurIkWFnICoisgXAMm8r5mcCQEFE8qq62reLc2EoEVFHJqKGFJE1OAGnrKqiqouqugggDqAKYEVEEn3LABeGEhF1ZCKCDoCr5mdDa0ZVbQDXzdtk367eMJGAM9eIiFqZlKCzDaeVY7dJc7ZvV1fvGh22dIiIWpmIMR1VXQ7aLyIx1Fs/N1t9XkRWAKwAwMWLF3vIAMd0iIg6MXE1pIjERKRkJhbsAFgAsKqqlVafUdW8qi6p6tL58+e7vyjHdIiIOjJxQQdOkLHMT1esr1fk4lAioo5MXA2pqraqxlV1Hs7sNRtAVkSy/bso1+kQEXViJMZ0TFdYN6pmSnRbqmqLSBrAFpwxm0wv+QvFMR0ioo6MRNBR1XivnxURC8A6gIqqNgUVE3iAfnaxNYzpMOgQEbUyCTXkNpw1OCtBB01QApxutv444L3XiIg6MfY1pKpWAZQBxPzjNmbKdMG8Xe9fJjimQ0TUibEPOsYqnNvdrHlu9LkJZ8p0As7C0Vzfrs4xHSKijozEmM5hmXGbZwFk4XS1JeAEoTKAdVUt9jUDT7wP+PFfc1o87znT10sREY2ziQg6QK2brX93km5nbgF47ueGcmkionHCviAiIhoYBh0iIhoYBh0iIhoYBh0iIhoYBh0iIhoYBh0iIhoYBh0iIhoYUe+jlgki8jqAl7v82DkAb/QhO+OK5dGI5dGMZdJoEsrjaVUNfQomg04ERGRDVZeGnY9RwfJoxPJoxjJpNE3lwe41IiIaGAYdIiIaGAadaOSHnYERw/JoxPJoxjJpNDXlwTEdIiIaGLZ0iIhoYBh0iIhoYBh0qCURWRGR5LDzQTSORGTNPMGYPBh0eiQiCREpmMdj75jXiWHnKyoiEgOwDiAdkq6rchjXchORlHkM+o6IqHm91ib9xJeLW6ma8nAfE59qk37iy8Rl/ljLwnmKcas0U1MeDVSVW5cbnEdiq9m2zOa+Tw07fxF9x4L5PutRlcO4lpunLNx8b3rfT2O5+PK46SuTpt+ZaSgTT95jAHbc/LZIMzXl0fRdhp2BcdvML5T7j51s8UsRG3Y+e/xuKTitmx20qUB6KYdxLTdP/nYAJHzfp+Qvo2koFwBrJl8l337L87vjL6uJLhNfOZQ8+dSA41NVHk3ff9gZGLfN8x8uG3Bs3RxbG3Y+e/xu3r9Ww4JOV+UwruWGeitnJeBYrTKYpnLx/J5YbX4vev6O41gmAXlfgWmNtEkz8eURtHFMp3vL5ufNgGPrvjRjRVUXVVVUVRAyloPuy2Fcy83tMy/7D6hqFYAN1MbAgOkol20AZVW126Q563k9DWUCM76ShVM27RZ7TkV5tDTsqDduG1r89eI5rgjo5x+3DU5XW7uWTlflMK7lBifoJFocC2rpTEW5tCkPd6zB2702FWVi8r0D09XV6ntMS3m02o70FKmmmxVRmnHXbTmMZbmpaiVov2nZ3DZvc55DU1EuLlMOBTh5tABUAaz6ym3iy0RE1uHkaVmdFnA7E18e7bB7rTftfqnCfuEmSbflMBHlZqbDbsJpBZVVNeNLMk3lsgCnwlvw7IsFpJvYMjHTxFcA5FS1qRu2hYktjzAMOr0J+k/VybFJ0205jHW5iYglIiU4s5MsOJVMUF/61JSLqtqqGlfVeQBxOGNcWRHJ+pJOZJmYlt4NAHbAHx/tTGR5dIJBp3vtBk+7STPuui2HsS43sxB0C8401SKAeItKZqrKxUudiQXuBJQVz6FJLpMrMJW+iJS8G0yXl2efe3ePSS6PUAw63XNnKzWtBPbsG6vmbo+6LYexLTfTX5+F8x0WVTWtrWduTXS5uK29gJYMgFrgARr/+p7oMjEsOH+QeDeX+94tk2koj5YYdLpXMD+vBhxz960HHJs03ZbDWJabiKzA+au9bLqRAicWeEx6uWzDqUBXgg6KiDug7Q3KE1smqppXs8zAv8GUgWdf0XxsYsujI8OePjeOG+oLJ73TQhPwTZ8d5w0hU6Z7KYdxLDfUpwB3vOJ70ssF9RX3Wd/+GOoLR9d8xya6TNr97kzj70jbchl2BsZxg/NXnvtL4A4su++bVq6P49Zh0OmqHMat3NB4+5GddtuUlYv3djc7aL73WingMxNdJi3KqV3QmbryqH2XYWdgXDc4XQze/2ib8NwXady3ToJOL+UwTuUGz1+SYds0lYvJr3sXcrcluGMqwpY3n5z0MgnIf9iizqkqD3fj46qJiGhgOJGAiIgGhkGHiIgGhkGHiIgGhkGHiIgGhkGHiIgGhkGHiIgGhkGHiIgGhkGHyBCRLRGZiIVrIpIQkU0RUXPD0pFl8rg57HzQYPDJoUSTyX2aZxnOynWikcCgQzSZLAAVDX7IHNHQsHuNaIjMkyf7dc7tqM9NdFgMOtR3IrJm+u0TZiuZ9zvmdcKXft1NH3CupDmWbXH+lGcsY8uXbt0dtzFpUm3ynPWlDXxomUmbEJGCL33T82Y8+bTMZ7bg3CizI+a7lUy5uWWX8qVZ95zTLau2Yzrd/vt0kx9f+qzn32bTPI21Xb46KtcWabfMv/dYPcp5Kgz7jqPcJn8DsAbzjBXzcwvOmIN7h2IFYHnSr8P37BDPsSR8z3LxnNf9nP/W7wU44xo7ntdBzyfZ8nzee9t+N+0WfM/V8Vzbvetvy1v8e9Imzbm3/GnalGGhzXWyvvLJevK7hpA7EXf779NNfkzaGJrvjuw+GsEt681DlGvSl7bkOf/WsH//ufl+34adAW6Tv/kqEP/DvdzKy1tx9hp0/EHEXxnFAq7hPY+3kk159nsrzXXPfstTSVst0q8F5HPHXw4h5ZdqcR3vc20Svus3Vc4R/vt0mx+3rP3/BinPdTd95+mmXN19SV/e3YDW8nEL3Aa/DT0D3CZ/81RqmwHH3MBQ8OzruaUTkL6pEmxz3S3/uT3HvA90i5l9pTb5dCvOnU7KIaT8ttpcx624SwF57TbodPrv03F+gsrNl37df+0eylUR/EyjhPluTefhNryNYzo0SDcD9tkRnj9oavA2AKhqpYvrNo2BqGoVzvRjwKn4AGAJQDXg3FBV21wjFjCuUPanD2G1uU7Rk5fD6vTfp5v8uGVVNGXoFzTe1G25VgDAjCklPWkrqpoLOg8ND4MODVKUASZIJLO1TMUWxK283Io0Bqfy06DNk27Bd56vdZoXEXHPsdEmmW3yclih/z495GfJs6/Ta3ZbrmlzniQAd2JDSURWOJFg9HCdDlEPPJVZFcD1kOQDmbosIrEWrYmhMGXU9ruralVE/J8BuihX80dC3LRy0nCCj7tlReQyWzujg0GHxo0VnuRwRMRq0dpxpw5XTGVZBQBVzfUrL6pqm0q5XfeZ293V94DTbX5ExPbsa+Kfjn2YclXVMkzXpWmRZQCswJkMEe/mXNQ/7F6jUebvlgKcv2T7bdW/w/wFnoRTmboV6QacbqCg9UQxt5sngvy4YxhB13HXxbTr7opaN/lxyyrVoqurqazRRbmaNU9bIlLwplNVW1Xdc/f9DxXqHIMOjaIt87OhQjILA5PNySO35l3kaCrL2+att8snY37e9ox1uApwxiYKODz3OgVvxW2uecO8bbl4tQ86zo9pfbktltu+9Ck4LZFW5w8tV/MHgAUnqDX8bnj+Ddm1NkIYdGgUuTOgUu5fseLchXgd9QqpX2w4lZS7ut1dyJgAUPZ2+ZhxghycSnDLpC+JyA6c4FhW1fxhM2RmhBVh1sGYlfmbMItVAeRN19JA9JCf63DKNOFJ7y5ALcM3maCHcnV/J0qetO75AeCz0X17OiwGHRo55q/XRTgV0gKctR8AsGwq/QyctRxRy8JZo7MIp9JbgFNRVgBkNODmmaqaAbDsyesSnEp0NSh9r1Q1DaflV4ZT2bt3kE57upEGppv8qGrVU6Zu8FkAkDNllIVv6nQ35Wp+J9Lm3Auot4bzAOKcRDBaRJ1FVERERH3Hlg4REQ0Mgw4REQ0Mgw4REQ0Mgw4REQ0Mgw4REQ0Mgw4REQ0Mgw4REQ0Mgw4REQ0Mgw4REQ3M/wel9xXqdTRc0QAAAABJRU5ErkJggg==\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"Dold=0\n",
"nn=arange(4,450)\n",
"rel_err=zeros(len(nn))\n",
"abs_err=zeros(len(nn))\n",
"for i,N in enumerate(nn):\n",
" w,K,F=int_N(N,10)\n",
" D=w[-1]\n",
" L=300; E=200e3; t=3; w=12\n",
" I=w*t**3/12\n",
" D_an=-10*L**3/3/E/I\n",
" rel_err[i]=abs(D-Dold)/D\n",
" abs_err[i]=abs(D-D_an)/D_an \n",
" #print(D,D_an)\n",
" Dold=D\n",
"\n",
"rel_err[0]=rel_err[1]\n",
"plt.plot(nn,rel_err*100,label='relative error')\n",
"plt.plot(nn,abs_err*100,label='absolute error')\n",
"plt.xlabel('number of nodes')\n",
"plt.ylabel('error (\\%)')\n",
"plt.legend()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Show log-log convergence rate of finite difference method**"
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"<matplotlib.legend.Legend at 0x7f28ebacf780>"
]
},
"execution_count": 22,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAa0AAAEeCAYAAAA5CErsAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3XdcXNed///XAYQk1AZ1yaoz6sWSAZe42wJX2XEBKS6/3U0xJBvvJptkIXIc20kc2+A4391sGjjJJt9vlFiC2I4tywXcW2wDkqxui1G3JEtCoy4EzPn9cWZgGjAM0+fzfDx4SLr3zL2HYXTf3HPOPUdprRFCCCESQVqsKyCEEEIES0JLCCFEwpDQEkIIkTAktIQQQiQMCS0hhBAJQ0JLCCFEwpDQEkIIkTAktIQQQiQMCS0hhBAJIyPWFUg2I0eO1FOmTIl1NYQQIqE0NDQc0lqP6qmchFaYTZkyhfr6+lhXQwghEopSamcw5aR5UAghRMKQ0BJCCJEwJLSEEEIkDAktIYQQCUNCSwghRMKQ0BJCCJEwJLTiQctxqPsRHPwk1jURQoi4JqEVDzY+C+/8HH51PvwuH+r/F84cjXWthBAi7khoxYO1f+n8+56PYNW34Wcz4G9fg6bXwemMXd2EECKOyIwY8eCSb0HWcPjkJXC2mW1tZ2B9tfkaOgEW3gEL74Th1tjWVQghYkhprWNdh6SSl5enQ57G6eQhE1JrlsOB9YHLTLoYzrsL5twC/QeHXlEhhIgjSqkGrXVej+UktMKrT6Hlad8602z48Uo43ey/v98gmPNFE2CTLoY0aekVQiQuCa0YCVtoubW1mGbDtX+BT2tBt/uXyZ4CC+40TYiWSeE7txBCRImEVoyEPbQ8Hd8PH68wzYeHtgYooGDq5bDwLph9E2RmRaYeQggRZhJaMRLR0HLTGvY2wtrlsKEm8PD4zCEw71ZYeDdMvACUimydhBCiDyS0YiQqoeWp9QxsWWWaD5teAwL8PEdMMyMPF9wBQ8dHr25CCBEkCa0YiXpoeTq6Fz5+yjQfNjf571dpYLvaBNjMG6HfgOjXUQghApDQipGYhpab1rD7A1fz4TNw9rh/mQEWmF9oAmx8jjQfCiFiSkIrRuIitDydPQmbnzcBtv2twGVGzTZD589dCoNHR7d+QgiBhFbMhBJaL3y8j3ebDoV8zqED+lGYO4Fpo3t42PjITlj3VxNgjl3++1U6TL/GBNj0ayEjM+Q6CSFEb0hoxUgoofWTVZv4/Tvb+3TefumKkstt3Hv1NAb0S+++sNMJO9814bXp79B6yr9M1giYv8QE2Nj5faqbEEL0JNjQkmkUkkRru+aXr2/j2v96i7c+Odh94bQ0mHoZ3Ppb+N4ncPMvYdIXvMucOgwf/AZ+e6n5+sdv4eThyH0DQggRBLnTCrNQ7rQadjazeV+AwRJB0MAzjXto3OXw2r743HE8sHgOo4f2YoTg4SYzdH7dX+HYXv/9af1g5nXm2a9p+ZAu8y0LIcJDmgdjJBYDMZxOzVMf7eaxFzdz7Exbx/Yh/TMovW4md144mfS0XowOdLbD9jfN0Pktq8yM874GjYYFS02AjZ4Vhu9CCJHKJLRiJJajBw8eb+GR1Zt5Zo33XdKCiRZ+ess85p0zrPcHPe2AjU+bANvbxfc1Psf0fc27HQZmh1BzIUSqk9CKkXgY8v7utkPc/+wGth862bEtTcGXL5nKfxTMYHD/EJv1Pt8C6/4C656CEwf896f3h1k3mgCzXgVpPQwIEUIIFwmtGImH0AI409rOb95o4jdvNHG2vXPl43HDBvDgTXO5du4YVKgPFLe3QdOrZvThltXgbPUvM2Q8LPiSmbx35LQQvwshRKqQ0IqReAktN/vBE9z/7Abea/Ie+Zc/ezQP3TyXCdl9nAn+VLNZuHLtcrMGWCATLzIzb8y9FQYM7dv5hBBJSUIrRuIttAC01jy7di8Pr9rM4ZNnO7YP7JfOt/On85VLp9IvPQxPP+zf4Fq4cgWcCvCwdMZAs3DlwjthymWycKUQooOEVozEY2i5OU6dpfylLfz1w91e22eNHcJPb51P7uQwDaJoOwufvuJauPJlcLb5l7FM6ly4MntKeM4rhEhYEloxEs+h5daws5n7nt7A1gOdz4YpBXdcMImya2cxLKtf+E524iCsX2lGH36+MXCZKZeZvq85N0PmoPCdWwiRMCS0YiQRQgugtd3J79/Zzn/VfcKZ1s6BGiMHZ3L/jXP44sLxoQ/UCERr2LfW1Xy4Es44/MtkDoa5t5hnvyZdJDPPC5FCJLRiJFFCy2138ykefG4jr2353Gv7pdNG8pNb5jF1ZATufNpaYOtqE2Db6kA7/csMt3YuXDlsQvjrIISIKxJaMZJooQVmoMbLG/fz0HOb2H+sc/aLzPQ0vn6FlX+9KohJeEN1bF/nwpWHPw1QQIH1SjjvbvMMWL+BkamHECKmJLRiJBFDy+1ESxtPvLKVP723A6fHx2JC9kAeumku+XPGRO7kWsOeelj7Z9jwNLQc8y/TfxjMu80E2Dm50nwoRBKR0IqRRA4ttw17j/KDZ9azbs9Rr+2LZo3mwZvmMmlEH5/t6knradi8ygSY/U3MtMA+Rs50NR9+CYaMjWx9hBARJ6EVI8kQWgDtTs2Kj3ZT8fIWHKc6Z7zon5HGv145jZIrrJFrMvTk2G2mjVq7HI4EWHNMpZkZ5xfeBTOvh4z+ka+TECLsJLRiJFlCy6355Fkef9n/2a5Jw7P40c1zuWrW6OhURGvY9b7p+9r4DLSe9C8zMBvmF5kAG7dAmg+FSCASWjGSbKHltmbXEX749w1s2Ovd13TNnDH8cPEcJg6PcJOhp5YTsPk5E2A73wlcZsw803w4fwkMHhW9ugkhQiKhFSPJGlpgmgz/8uEuHn9pi9e6XQP6pXHvVdO453Ir/TOiPLN7s93VfPgXOLrbf39aBky/1sw8P/0aSA/jg9NCiLCR0OojpVQ1UKm1ruvN65I5tNwOn2ih/KUtrKzf47V96shBPHTzXK6YEYM7G6cTdrxlwmvTc9B22r9M1kg4d6kJsDFzo19HIUSXJLT6QCmVD9QCBRJaXWvY2cwPn93Ipn3eTYbXzxvL/YvncI4lRs9UnTkKG581gzd2fxC4zLiFpu9rfiFkDY9u/YQQfiS0QqCUKgSeBCyuTRJaPWhrd7L8g1387JWtHPdoMhzYL51/WzSNr11qJTMjhrO5H/rU3H2tewqOf+a/Pz0TZt5gAsx2NaSHuECmEKJPEja0lFLFgL23YRGmc1sAKzAcudPqlYPHW3jsxS38rdG7ydA6ahA/vnkel04fGaOauTjboel118KVL0B7i3+ZwWM7F64cNSP6dRQihSVkaLlC4whQpbUu6aZcDrAMyHdtqgMe1Vo3hrEuGgmtXvtwezMP/H0DW/Yf99p+47njuP/G2YwbFgfTMJ0+AutrzB3YZ118ZCacb8Jr3m0wYFh06ydECkrU0KoGCukmtDz6mwDsrj+trj+LtNY1YaqLhFaI2tqd/N/3d/Lz2k840dLZZJiVaRad/PIlYVp0MhwObDJ3Xx+vgJMH/fdnDIDZN5kAm3qFLFwpRIQkTGi5+pEKgCV09iUFDC2POzHwCBSfIMvWWgdY96LX9ZLQ6qPPj53h0Re38MyavV7bp48ezI++OJeLbTFuMvTU3mpmnF/zZ/jkpcALVw6dYBatXHinmYVeCBE2iRRaDUCOz+auQqsUKAcqtNZlPvsqgWKgTGtd4dpWDOT2UIXaQHdnElrh8w/7YR74+wY+OXDCa/vic8fxg3hpMvR08hCsrzYPLx9YH7jMpIvN0Pk5t0D/wdGtnxBJKGFCy5PrrquarkOrFtOPlevbf+Xq52oA6rTWBWGoi4RWGLW2O/njuzv4r7pPOHm2vWP7wH7p3Hv1NL522dToP5gcjH0fu5oPV8LpZv/9/QbBnC+aAJt0sTQfChGiZA2tJsCqtQ44qZwraOxaa1sY6iKhFQH7j57hkdWbeW6d9/DzySOyeGDxHBbNjuDyJ33RdtY0G65dDp/Wgm73L5M9BRbcaZoQLZOiXkUhElmyhpYG6CG0utzfy7oEHVquZshigEmTJuXu3Lmzr6dPev+wH+ah5zb6jTK8auYoHrhpbmRWTA6X4wfMwI21y+HglgAFFEy93AzemH0TZEZxXkYhElQyh5ZDa53dxeuPAJZwhFao5E4reO4Hk594ZavXXIaZ6Wl89bKp3HvVNAb1j+OHfbU2Q+bXLIcNNWYmDl+ZQ2DerbDwbph4gcw8L0QXkjm0onKnFSoJrd47fKKFn72ylac+2o3nx3Hs0AHcd+Nsbjp3HCreL/atZ2DrCybAml4j4MKVI6a5Fq68A4aOj3oVhYhnyRpaUevTCpWEVug+3uPgwec2smaX9xMLF04dzkM3z2X2uKExqlkvHd0LHz9lAqy5yX+/SgPrVWbwxswbod+A6NdRiDgTbGgl2lAnO3SMFPTisa3Pz2iJ2Dh3goW/ff1ifla0gJGDO1cg/mB7Mzf+4m0e/PsGjnqsohy3hp0Dl30X/q0BvvIy5PyTaSZ0005oehVqvgJPzIBV34G9DRBHv0AKEa8SLbSqXX8uDbDPva0ySnUREZCWpijMncBr37uCr106lYw0c1Pt1PCn93dy1RNv8NcPd9HuTIALvFIw6SK4+X/ge1vh1iozQMPTmaNQ/3t48mr49Rfg3V+YgR5CiIASqnnQVcZd4Y5ntTye0YppfxZI82C4fXrgOA89v5F3tx322j7/nGH86ItzyZkUcExOfDuy07Vw5XJwBBhpqtLNgpUL74QZ10FGZvTrKESUJWWflqtMMZ13U+7h6O6Jc0u01lWRrWX3JLTCT2vNSxv28/ALm9nr8F7c8facCZRdP5PRQxKwX8jphJ3vuhaufBZaT/mXyRoB85eYABt3bvTrKESUJG1oucrlY6ZzcvdjNWKmb4r6cia+JLQi5/TZdn7zZhO/fbOJs23Oju1D+mfwrfzp/PPFU+JnIt7eajnuWrjyL7DrvcBlxs43Q+fnF8GgEdGtnxARlpChlQwktCJvd/MpfrJqE69s8u77mTZ6MA/dNDf2a3f11eEm18KVf4Vje/33p/WDmdeZAJuWLwtXiqQgoRUjElrR89YnB3no+Y3YD5702n7d3LHcv3g2E7ITfCYKZztsf9MMnd+yCtrO+JcZNBoWLDUBNnpW9OsoRJhIaMWIhFZ0nW1z8sf3tvPfdZ96TcTbPyONb1xp4+tX2BjQLw4n4u2t0w7Y+LQJsL1dfL7G55hnv+bdDgMTcICKSGkSWjEioRUbnx87w2MvbuFpn7W7JmQP5P4bZ3Pt3LHxP6tGsA5uNSMP1z0FJwIMj0/vD7NuNAFmvQrSkiC0RdKT0IoRCa3Yqt/RzIPPbWTjZ8e8tl9sG8GDN81l5tghXbwyAbW3mSmj1v4Ztr4I7Wf9ywwZDwu+ZCbvHTkt+nUUIkgSWjEioRV77U7NUx/t4vGXt+LwmEEjTcHdF03mOwUzsGQl2bNPp5phfY0JsH3rApeZeJEZOj/3VhiQIFNiiZQhoRUjElrxw3HqLP+n9hP+/IH3DBqWrH58t2AGd1wwiYxEHSLfnf0bzOjDj1fAqUP++zMGmoUrF94JUy6ThStFXJDQihEJrfizdf9xfvT8Rt5r8p5VY9bYITxw0xwutiX4EPmutJ2FT18xAfbpy+Bs8y8zbJJZtHLhnWYRSyFiJOKhpZS6GrACNsACNGEmq63XWq8N6aBJQEIrPmmteXnjAX66ehO7m71n1bh+3ljuu2E2E4cn+BD57pw4COtXmtGHn28MXGbypWbwxpwvQmYcL8IpklJEQksp9TWghM6ZKAINx9KY8KoFKrXWrwd9giQgoRXfzrS28/t3tvPL17ZxutV7iHzJ5Va+fqWNrMwkflhXa9PntXY5rK+G00f8y2QOhjm3mACb9AVZuFJERdhCSyk1FFgGlGJCyo6ZNqkWaHb9G2A45o7rfEyo5WMCzA6Uaq2fCek7STASWolh/9EzPPbiZp5d+5nX9nHDBvD962dx84LxyTNEvittLWbU4drlsK3OLJniK3uqGXm48A4YNiH6dRQpIyyhpZS6Dfid659VmDun7b2oRD7mzux24BVgidb6WPevSmwSWomlYWczDz23ifV7j3ptP39KNg/eNJd55wyLUc2i7Ng+M3Bj7XI49EmAAgqsV5oAm70Y+g2McgVFsgtXaDVjJqJ9so+VsQCPAUe01sv6cqx4J6GVeJxOTU3DHipe3sKhE53POikFS/Mm8r1rZ3otSpnUtIY99Sa8NvwNWgL8jtl/GMy7zQTYhDxpPhRhIaMHY0RCK3EdO9PK/7z6Kf/77g7aPIbIDxmQwbcWJfgs8qFoPQ2bV5kAs7+Bae33MXKmGXl47lIYOi7aNRRJJKqh5dHv5V6M8bFkbwbsioRW4ms6eIKHV23i9a0HvbbbRg3ih4vncOXM0TGqWQw5dncuXHkkQA+BSjMzzi+8C2ZeDxkpcmcqwiZqoaWUGoYZbJGNGaCRgxmgYU3F4JLQSh6vb/mcn6zahP2Q9yzyi2aN5v7Fc5g6MgWHhWsNu943Q+c3PgOtJ/3LDMw2a34tvAvGLZDmQxGUaIbWSsxAi2la6+1KKSuwDViptf5Snw6egCS0ksvZNid/em8H//3qp5xo6Xw4t1+64iuXTOXeq6cxZEC/GNYwhlpOwObnTIDtfCdwmdFzzdD5+Utg8Kjo1k8klGiGVjOwQmv9DY9tlUCh1jrllleV0EpOB4+38PjLW6hu2IPnf5lRQ/pTeu1Mbs+ZQFpaCt9RNG83i1au/Ssc3eW/Py0Dpl9rAmz6NZCeokEvuhSu0YNXA3at9Y5uyjQDr3jeVSmlfgsUSWiJZPPxHgc/en4TDTu9H8pdMGEYD948l5xJKb6OldMJO942fV+bnoO20/5lskaagRvn3QVj5ka/jiIuhSu07gEqgd8C3w/UR+W6q/oakKO1XqeUOg8zGKPS8+4rVUhoJT+tNX9f+xmPvriZA8davPbddt45lF0/izFDB8SodnHkzFHY+KwJsN0fBC4zbqHp+5pfCFnDo1s/EVfCOSNGPia0pgLlWuv7fPYPwwzAmIIZkGEFtmNCTAZiiKR1sqWN37zRRNXbds62dc4mkZWZzr1XT+Mrl0xNjlWTw+HQp2bi3nVPwfHP/PenZ8LMG0yA2a6G9CSeSksEFPY+LaVUMeYBYY2Zlun3PvtLgTzgI631472vcnKQ0Eo9uw6f4pHVm3lp436v7ZOGZ3H/jbMpmDMm+aeECpazHZpeN3dfW16A9hb/MoPHdi5cOWpG9OsoYiJiAzFc4fQYZlb34lSbELcnElqp691th/jR8xv55MAJr+2XTR/JDxfPYcaYJFo1ORxOHzGzbqxZDp81Bi4z4Xzz8PK822FAikyplaIiOnrQNS1TOXAPZk7BEq31zl4fKAlJaKW2tnYnyz/Yxc9rP+Ho6c5Vk9PTFHddOIn/yJ9B9qAkWzU5HD7fbO6+1q2Ak5/7788YALNvMgE29QpIk2bXZBOVIe+uZ7IqgatdfwYcrJFKJLQEwJGTZ/l57Scs/2AnHjNCMWxgP76dP527L5qcWlNCBau9Fba9Cmv/DFtfAmerf5mhE8ys8wvugBG26NdRRERYQ8s19L0SM8iiCXNn9brH/nxMk+F5BBiskUoktISnLfuP8ePnN/mtmjxt9GDuv3F2ak4JFayTh82aX2v/DPvXBy4z6eLOhSv7S/NrIgvn6MFFmCbANcAK4EvAQiDftz9LKVWIWcIk4GCNVCChJXxprXll0wEeWb2ZnYdPee27auYo7l88B9uowTGqXYLY97EZffjxCjjd7L+/3yATXAvvhMmXQJrcxSaacIZWPTBMaz3dY1sTcFhrfUEXrykFHtVap1zDs4SW6EpLWzv/++4OfvnaNq8poTLSFP/0hSl8a9F0hmXJTBHdajsLn7xkAuzTV0C3+5exTDYjDxd8CbInR7+OIiThDC0nPg8Ku2a8uKe7UFJKDdNaH+1qf7KS0BI9OXi8hZ+9vJWVDbu9poTKzurHd66ZyR3nTyRD+rt6dvxA58KVB7cELjP1clh4txnEkZkV3fqJXglnaG0DmrTW13psewWY6nn3JQwJLRGsDXuP8uNVm/hwu3dz18wxQ3jgpjlcMm1kjGqWYLQ2Q+bXLIcNNWYmDl+ZQ2DerSbAJl4gM8/HoXCGlvu5rMeAlcBSoAzTZ/WzMNQ1qUhoid7QWvPihv389IXN7HV4z9N3zZwx3HfDbKak4hIooWo9A1tfMAHW9BoBF64cMc30fS24A4aOj3oVRWDhHj1YiXkmSwMKqNJaf73PtUxCEloiFGda2/nd23Z+/UYTp8529tNkpqfx5UumpPYSKKE6uhc+fsoEWHOT/36VBtarzOjDmTdCP5kvMpYiMY3TVMwCj41a6wBLlwqQ0BJ9c+DYGcpf2sLTjXu9to8cnMn3rplJUd5E0lN5CZRQaA27PzRD5zc8A2eP+5cZMAzmFZoAG58jzYcxEK5Z3oeG82HhcB8vHkloiXBYu9vBj5/fSOMuh9f2ueOH8sDiOVxoTblVf8Lj7EnYvMoE2Pa3ApcZNds0H567FIaMiW79Uli4QsuJqz+rL2Hjeji5HKjTWi8L9TiJQEJLhIvWmufWfcZjL25h39EzXvtunD+O718/i4nDZURcyI7sNLPOr10OjgCz0Kl0s2DlwjthxnWQIdNvRVK4QsuKGXxxHlCNGfoe1AS5SqkpQCFQgplJoyLZAwsktET4nT7bTuVbTfz2zSbOtHYugZKZkUbxZVa+caWNQf1lKY+QOZ2w6z3T97XpWWg95V8mawTMX2ICbNy50a9jCgj3QIxioBQTPhqzflY9Zkond/uFBRjhKpPv+rcCaoCyVOkHk9ASkfKZ4zSPvbiF59Z5r0c1ekh/Sq+bxW3nnUOa9Hf1Tctx2PR3E2C73gtcZux8M3R+fhEMkmbacInIhLmuOQaLgEWYcOpKHVCLGWWYMA8Yu2avf5LO0K3DBG4X6yb4k9ASkdaws5kfPb+Jj/d4/9daMGEYD9w0h9zJsgJwWBxugnV/hbV/hWN7/Pen9YOZ15nZN6YVyMKVfRTxWd5dKxZbXV/DMasWN2ut14R0wDiglGrAhFUZ5vspxyxsOVVr7ejutW4SWiIanE7N02v2UvHSFj4/7r2Q4hcXjqfsulmMtwyMUe2SjLMdtr9ppo7a/Dy0nfEvM2g0LFhqAmz07OjXMQlEZWmSZOLqv2sCcj3vrJRSGijSWtcEcxwJLRFNJ1va+PUb23jy7e2cbevs7xrQL42vX2Gj5HIbAzNTbgrQyDntgI1PmwDb81HgMuNzzND5ebfDwOzo1i+BJWxoufrP7Frruiif14JZibnCZ9sRoCDY+khoiVjY3XyKx17cwgvr93ltHz9sAGXXz+LmBeNR8uxReB3c2rlw5Yn9/vvT+8OsG02AWa+ShSt7kJCh5RESVVrrkm7K5QDLMH1PYPqeHu1N31MP9chxHXsppsmzINjXSmiJWPrAfpgfr9rExs+8n1DJnZzNA4vnsGCiJUY1S2LtbWbKqLXLYetqaD/rX2bIeDPr/MK7YOS06NcxASRqaFVjhsl3GVquwSC1rn/aXX+6B4UE3YzXQz1K6Ryq322A+pLQErHW7tTUNOzm8Ze3cuiE9wX0tpxzKLtuFmOGypRFEXGqGdbXmADbtzZwmYkXmvCaeysMGBrd+sWxhAkt18KRBcASzCAI6CIoPO7EwKPJzifIsoMdNBFE3axAg6s+ZcG8RkJLxIvjZ1r55Wvb+MO722lt7/x/npWZzjevmsZXL53KgH7SZBUx+zd0Llx56pD//oyBMOdmE2BTLkv5hSsTKbQaMHMaeuoqtEoxI/oqfEPENalvMWaIeoVrWzGQ20MVarXWNa4mwTytdZXPcasBS7BNhBJaIt7sOHSSn67eTO2mA17bJ2QP5L4bZnP9vLHS3xVJ7a1mwco1y+HTl8HZ5l9m2CRYeId5eDl7StSrGA8SJrQ8ue66quk6tGoxfU25vv1XrtBpwEwVFXQflMfr3XdrXndqrlCtD7aJUEJLxKt3tx3ix89vYusB7wljL5g6nAcWz2HeOcNiVLMUcuIgrF9pAuzzjYHLTL7UDN6Y80XITJ1laZI1tJoAq9Y64K+FruHpdq21LcTzH8HM9FEONGP6tYoJEJJdkdAS8ayt3clTH+3miVe2cuRUa8d2pWBp3kS+e81MRg3pH8MapgitYd860/e1vhpOH/EvkzkY5txiAmzSF5J+5vlILE0yDHMhr9VaX9fH+nV1jp5CSwP0EFpd7g/i/FZMYLlHJdYTxIwYrmbIYoBJkybl7twZYPJNIeLI0dOt/OLVT/nTeztoc3ZeAwb3z+CbV03jy5dMkf6uaGlrga0vmgDbVgfa6V8me6rp+1p4BwybEP06RkGkpnFqApxa6+l9qVw3xw8mtBxa64BP7LnulCyhhlY4yJ2WSCRNB0/w8KpNvL71oNf2icMHct/1s7lO+rui69g+M3Bj7XI49EmAAgqsV5oAm70Y+iXPrCeRCq0czDNRP9VaP9GH+nV1/JjeaYWDhJZIRG9s/ZyHX9jMts9PeG2X/q4Y0Rr2NsCaP8OGp6ElwBSu/YfBvNtMgE3IS/jmw0iF1hTAhgmWJkyAdTGXCWitnw764MS+TyscJLREomprd/KXD3fx89pPcPj0dxXmTOA/r53JaHm+K/paT8OWF0yA2d/ALLThY+TMzoUrh46Ldg3DIlKh5cS8Y56hEegACtBa6141iodp9GCj1rqnYe4RI6ElEt3RU6384jX//i55visOOHbDx0+Z57+a7f77VRpMyzcBNvMGyEicQTWRCq1SAodUQFrrx4M+OEGFVjFQSeDntMoxa36V+D5rFU0SWiJZ2A+e4JHVm6nb/LnX9nMsA/n+9bNYfO446e+KFa1h1/um72vjs3D2hH+Zgdlmza+Fd8K4hXHffJiUQ95dZdwV7rjb8rjLiml/FkhoieTz9qcHeXjVZr/nu/ImZ/NDmc8w9lpOwObnzN3XjrcDlxk91wydn78EBo+Kbv2CFNXPZXSfAAAgAElEQVTQUkoN1Vof67lkj8cJJrTcd1tg+tSgc4h6TO+yQEJLJKe2dicr6nfzxCuf0HzSfz7D0mtnMXaY9HfFXPP2zoUrj+7y35+WAdOvNQE2/RpI7xf9OnYhoqGllFpI5wKJ7l+zHJhBGeVa69d7fVCCCy1XuXzX+d3TPzVinqeK6nImgUhoiWR29HQrv3p9G//rM5/hwH7pfONKG/dcZpX1u+KB02nuutYuh03PQdtp/zJZI83AjfPugjFzo19HHxELLaXUo5i+I3cznO9M6xqo1Fr/a68OnCQktEQq2HHoJI+s3swrPvMZyvpdcejMMdj4jAmw3R8ELjNuASy8G+YXQtbw6NbPJVIDMW7H3Ak5gHu01n/z2V8IPAkMBQq11s/0qtZJQEJLpJL3mg7xk1Wb2bzPu3fgvEkWfrh4DjmTZOXeuHJom2vhyqfg+Gf++9MzYeb1JsBsV0N6RtSqFqnQegVYBNi01ju6KGMFtgGvRGq6p3gmoSVSTbtTU12/m5+94r9+1y0Lx1N63SzGW5Jn5oak4GwH++tm4t4tL0B7i3+ZwWPNwpUL7oBRMyM++jBSodUMNGmtz++hXD0wVWs9IuiDJwkJLZGqjp9p5VevN/GHd7Zztr1z/rwB/dIovtzG16+wkpUZvd/cRZBOH4ENfzMB9lkX06xmTzWzzi+8C0bNiEg1JLRiREJLpLpdh0/x6IubeXHDfq/tY4b2p+y6Wdyy8BzS0qS/Ky59vtnVfLgCTn4euMyECyDn/zPPgIVx7sNIhVYtcDVmKqWAU5m7ZoM/gpkN/tqgD54kJLSEMP5hP8xPVm1i42fe/V0LJgzjgZvmkDs5Nh3+IgjtrWbG+bV/gabXAj+83H8ozFoM82+HqVf2uf8rUqFVCKzE9FkVaa3X+exfiBmoYcU8M/W7XtU6CUhoCdGp3an5W+MeHn95KwePe/ebLD53HN+/fhYTsrNiVDsRlLYWE1xrl5slVAKtvDx0AlxYDDn/DANDe9g8kkPeq4HbMUPbHZg1p6DzmS0FVGutl/bqwElCQksIfyda2vjNG9t48u3tnG3r7O/qn5HGPZdZ+caVNgb1l/6uuHfykLn7qv8DHNnuv7/fILj3w5DW/Ir0w8WFmId7p/rssmMe8v2b/6tSg4SWEF3b3XyK8pe2sOrjfV7bRw/pz39eO5PbcyZIf1ci0Br2NsLaP8Omv8Opw2b7hAvga7UhHTJq0zgppaYCaK0DxG7qkdASomcf7WjmJ6s28fEe73Wi5p0zlAcWz+WCqdLflTBaz8D6lfD+r+HKMph7a0iHCXtouQZYNGMGWKTc81fBktASIjhOp+aZNXupeHkLB45593fdMH8sy66fzcTh0t+VMLQ2X2lpIb082NAK+uha66PADswikEII0SdpaYrbcyfw2nev5N+vnkb/jM7L0er1+1n0xJuUv7SF42dauzmKiBtKhRxYvdHbMxQBI5RS341EZYQQqWdQ/wy+c81MXvvelXxx4fiO7WfbnfzmjSau+tmbPPXhLtqd8bOMkoid3g55n4K506oGmjBLg3zUVXmt9dN9q17ikeZBIfqmYecRfrJqE2t3O7y2zx43lAcWz+ELtpSbsyAlROo5LSdmqLvn8J5AB1CA1lqn3BoFElpC9J3TqXlu3Wc89uIW9h8747Xv2rljuO+G2UweMShGtROREGxo9fbBiO8TOKSEECJs0tIUt5x3DtfMHUPVW3Yq37RzurUdgJc3HuC1LZ/zlUum8s2rpzF0QPwsZCgiLywrF4tOcqclRPjtO3qax1/aytNr9nptHzEok+9cM4MvnT+JdHm+K6GFffSgUmqYUqpdKfVS36omhBC9M27YQH6+dCHPfvMSciZ1ThN0+ORZfvDMBm78xdu8u+1QDGsookWGvAshEsbCiRb+9o2L+Z87zuMcjzW6tuw/zl2/+4Cv/ame7YdOxrCGItJkyLsQIqEopbhpwXhe/e4VfO+aGWRldo73qtt8gGv+z5v8ZNUmjp6S57uSkQx5DzPp0xIiug4cO8PjL2+lpmGP1/bsrH58p2AGd1wwiYz0yD/0KvpGhrzHSG9Dq62tjebmZo4ePUpbW4Ap/4UQPcrIyOCUzuTxN/fyTtMRr33TRw/m/sVzuGLGqBjVTgQjUqFVSi+GvGutHw/64EmiN6HldDrZuXMn/fv3Z8SIEWRmZqKUjIASoje01pw9e5bDhw/T0tLC5uP9eeTFLew5ctqr3FUzR/GDG+cwbfTgGNVUdCdqs7wLb70JrcOHD3Pq1CkmTJggYSVEH2mt2bNnD1lZWQwaauEP727nV69t4+TZ9o4yGWmKuy+azLfzp2PJyoxhbYWvsA957+IkQ5VSU5RSQ/tynFR14sQJLBaLBJYQYaCUwmKxcPLkSQb0S+dfr5zG6/95JUvzJuL+L9bm1PzxvR1c8fgb/PHd7bS2O7s/qIg7IYWWUup7SqnDwBHMgIxi1/aXlVJfDWP9ktqZM2fIypKlF4QIl6ysLE6f7mwWHD1kAOWF5/L8vZdyoccaXUdPt/LQ85u47r/e4vUtnyMtTomj16GllHoZs2pxNrAd70EZFwBV8gBycJxOJ2lRmMpfiFSRlpaG0+l/9zTvnGE8VXwRv707l0kea3Q1HTzJl//4Ef/0hw/Zuv94NKsqQtSrK6ZS6h6gAFgDWLXW03yKTAVeAwrkjis40jQoRPh09/9JKcV188ZS+53LWXb9LAb375x69e1PD3H9f7/FD55Zz+ETLV0eQ8Reb3/NL8GMHrxaa73Dd6fW2qG1LgCOAl/ve/WEECK8+mekU3KFjTf+80ruvHAS7ikLnRqWf7CLKx9/g8o3m2hpa+/+QCImehtaOUCj1vpYD+XqAWtoVRJCiMgbObg/j9w6n9XfuozLpo/s2H68pY1HX9xCwc/f4sX1+6S/K870NrTsQZaTwBJCJIRZY4fyf79yAX/4lzysozrX6NrVfIpvLG9kadU/2LD3aAxrKDz1NrTWADlKqcldFVBKTcWElsxlJGKurKwMpRR1dXV9Oo7dbsfh8F5JN1zHFrGnlOLqWWN4+duX89BNc7Bkda7R9eH2Zm765Tt8r3odB3wWpBTR19vQehQzWrBWKbXAY7sGUEotxISVxowwFCLh2e12bDYb99xzT6yrIiKsX3oa/3LJVN743pV85ZKpZLg6vLSGmoY9XPWzN/ifVz/lTKv0d8VKr0JLa92IWb14GtDoelZLA/e5/t6AGQpfobV+LdyVFSKelJSUUFtbS15ejw/xiwRjycrkgZvm8PJ/XE7+7NEd20+dbeeJ2k+4+mdv8Pe1e6W/KwZ6/ZCQ1roCE1qvYe66FCaoFPAqkKu1XhbOSgrh2zQXD6xWK/n5+Vgslp4Lx7nu3t9wvvfx+HPsjm3UYH73z+fz569eyMwxQzq2f3b0DN96ai23/eY9Gncd6eYIItxCerJVa23XWhdorYdrrdOAbNffr9FarwlzHUWKqaioQCmF3W6nsbERm81Gdna2V5nGxkaKioqw2WwopcjNzaWqqqrX58nNzSU7OxulFDabjbKyMq8La0FBATabWfe0pqYGpRRlZWUA1NXVUVFRgd1uxicVFRWhlOqyHiUlJQH3h+N76c1xunt/g3nva2pqKCgoIDs7m+zsbAoKCqipqenVeRLNpdNH8sK/X8pPb53HiEGdcxau2eXgtl+/x7//dQ17jpyKYQ1TiNZavsL4lZubq4O1adOmoMumkvLycg3o2tpabbFYtNVq1fn5+X77AZ2Tk6NzcnI6/u1ZTmutS0tLO47lyf0ai8Wi8/PzdX5+vtcx3aqrq3VxcbEGtNVq1eXl5R3H8j12bW1twDq4WSwWDegjR46E9L0E854Fc5zu3t+e3vvCwsIuz1NaWhr0eSIp0v+vjp4+qx9ZvUlPv2+1nly2quNr+g9W60dXb9bHz7RG9PzJCqjXQVxjOx8JF3FnyvdfiHUVQrbjsRv7fIyioiKWLVtGaWlpxza73U5ZWRlWq5Xa2lqsVvN0hcPhYNGiRR13P56v8VVVVUVjYyP5+fnU1tZ67bPZbDQ2NmK327FarRQWFpKTk0NVVRU5OTndHtfdVFhXV4fD4fBqNnRv82xODMf30pfjBHp/u9tXU1NDTU2N33nsdju5ublUVFSwdOlScnJygj5PIho6oB/Lrp/NXRdM5tEXN/Pihv0AnG1z8ts3m3i6cQ9l183i1vPOIS1NZrwJN5n4zoNSqlgppQN85ce6bqnIarX6XehKSkoAqK6u7rhoAlgsFqqrqwF49NFHuz3u8OHDKSwspLzcf4BrYWEhQEeTX28tWbIEwG8YvLtu7vp7/r0v30tfjhPo/e1un7tZ1Pc8VquVJ5980qtMsOdJZJNGZPGbu3NZWfIFFkzs/AXl8+MtfLd6Hbf+5j3WSH9X2EloebNhHqAu8vmSZ85iID/f/3eF+vp6LBaL32/zYC6OVqsVh8PRbYd/YWEh1dXVfsdobGzs8zNXRUVFAKxYscJr+8qVK7FYLB2hCOH5XvpynEDvb3f77HZ7l+dxf1+B1pLr7jzJ4IKpw3nmGxfz8yULGD2kf8f2dbsd3Prr9/jX5Q3sOHQyhjVMLnHXPKiUKgbsWutYPLFpBeq01v69yjEQjia2RHb++ef7bXNfeHuaaLi5ubnbUX0Oh4OVK1fS0NBAfX09jY2Nfausi7v5z3NggrtpsLi42K8OEJ7vJZTjBHp/u9rnvvPsbni/1WoNeIfa3XmSRVqa4racCVwzdyy/fn0bv3t7O2dda3WtXr+f2k0H+KcvTOHfF01n2MB+PRxNdCeuQkspZQEqgSqgy9BSSuUAywD3r3B1wKPaPEfWF1bgI9c5rFrr0NqIRES4L84Wi4Vly7p/qmL48OFd7qurq6OgoADoHLZeUlJCXl4eK1asoKKiok/1XLJkCVVVVdTV1ZGfnx+waTBc30u4jhNOvv15qWRw/wxKr5vF0vMnUvHSVl5Yvw+A1nbN79/Zzt8a93Dd3LF855oZjB4yIMa1TUxxFVrAkz0VcPUvuXvP3aFSCBQqpYr6eJeUA1iUUuWucwGUafNsmogxi8XScTHsSx+JZx+QZ3Md+DfrhaKoqIiqqiqqq6vJz89n5cqVWK1Wr2a1cH0v4TpOT9x9WIGa/9zczYepGlieJo8YxK/uyuGe3Q5++sImPtph+rYcp1p56qPdrF6/j1/flculHhP1iuDEvE9LKVWolKpUSh3BhE93ZS10BlaB1tqmtbZh1vgCqHaVCaUe7tc1Yh6WzgbKgHKlVLf1EtGTl5eHw+EI2JzncDg6nhvqjufIwED7+srdRLhy5cqOpkHPuyy3cHwv4TxOT9x9Y4HO424OldlBvC2caGFlyRf45Z3ncY5lYMf2Y2fauPv3H7Dkt++zeV9Pi2YITzEPLUwzXzEQTNi4OwUqPPu8XH+v8inTK9qsBaa01kWuvztcd1h1mHXERBxwj/hbtGiRX8AUFRXhcDg6BkN0xWKxYLfb/V5fVVXVcfENNPihN7M5LFmyxCusAgVkOL6XcB4n2PO4j+lmt9s75mUMNHow1SmlWHzueN4qvYoHb5pDZkbnZffDHc0s/p93eOi5jXwuk/EGJeahpbXOdYWFwozU647718VAbTiVPmXcQ9gre/jq6S6qEYhOZ4DokftZKYfDgc1mw2azdczO4O5D8h3w4Mu932azUVRU1DGLRFlZWUe4lJWVdcwm4e4Lqquro6ioKODsD77cIWG328nJyfEaIh7O7yWcx+lJYWEhhYWF2O12srOzyc3NJTc3F5vN1jHQJNlHCvZFepriy5dM5YV/u5TLpo/sWHyy3an543s7uOCRV/mX//2Qo6dbY1vROBfz0OolK3RM3OvFY5vVY1uV1rqkh68aMH1lSqkjSinfq0s+MuQ9rpSXl1NbW0t+fj7Nzc3U19djtVqprKz0e1i4q9dXVlZitVqpqanpeNB4+/btHf1QdrudhoYGwNyZlZaWdowKDKYJ0fMh4kBNg+H6XsJ9nJ5UV1dTWVnZ8R7Z7faOwSaVlZU9H0AwfcwQ/t9XL+Slb1/OtNGDvfa9sfUgl5a/xrNrZDLerqh4emNcdz3VQJXW2u9/ulJKA7juygK9vtv9QZy/CXBglmBxYO78igFbdyMJXcP0iwEmTZqUu3PnzqDOt3nzZmbPnh1KVYUQXUik/1ftTs1rWz7n639uoN3pfS3Oykznj1++gAumpkZDj1KqQWvdY6doot1pgQmTUPYFIxczIrEcE55Weggs6Lijy9Na540aNaqPVRBCpIr0NEXBnDE03l9A3mTvyYRPnW1nSeX7TPn+C6zdnViz40dSIoZWdwM2+jTW1jX4osg1KjFbm5ns5VktIUREDcvqR803LmbdA9dw4/xxfvtv+dW7fHN5I80nz8agdvEl0UIrmACRkBFCJKRhWf341V05/OVrFzLvnKFe+15Yv4+cn9Qy5fsvsLs5dZdBScjQcs2I4cVjm9xHCyES2sXTRvL8vWaUYSCXVbzOHVX/4EAKDpNPtNCqdv25NMA+9zYZwiSESHhKKf7fVy9k04+vZfa4oX7737cf5sJHXuVP7+3gTGt7DGoYGwkVWlpr9wPEpZ53W66/l/qUEUKIhJeVmcGL37qM8tvnB9z/4HMbmfXDl3h2zd4o1yw2Eiq0XNxD4RuUUrVKqVqgwWefEEIklaXnT+L9ZVdzxwWTAu7/9oq13PW7f3D8TCsnWtqiXLvoibcJc3ukta5SSrmHpbsfv2/ETGwbi+VMhBAiKsYNG8ijt81nV/NJ3t122G//u9sOM/+hVzr+fe9V0/hOwYykWkE5rkLLNTtFj++uK5xyI18jIYSIP7+8I4e6zQeYOnIQhb99v+tyr2/jl69vY/OPr2NgZnoUaxg5idg8KIQQKS17UCZFeRPJmzKcTx6+nmvmjOm2/OwHXmLX4VMcOXmWGfe/yA3//TYHj7dEqbbhFVd3WkIIIXonMyONqn/Ko92psd23ustylz/+esffN+07xkPPb+RXd/o9PRT35E5LCCGSQHqaYsdjN3Ljuf4zagTywsf70Fon3MS8ElpCCJFEfnVnDh/9ILglYqYuW83UZavJe7iONz85SGu7k//3/g5W1u/G6YzPMJPmQSGESDKjhvRnx2M3surjz7j3L2t6LH/oRAv//IcPOccykL2O0wCMHJzJ1bO67ytza2lrp39GdAZ6yJ2WSAo2mw2l4mNYb1lZGUop6ur69gSG3W7v1WrJQvhafO54PrhvEWOG9g+qvDuwAL7yx3pe23IgYPPhoROdgzieX/cZC39Uy5eq3o/K3ZmElhBxyG63Y7PZOpaxFyJUY4YO4N2yq0N67Vf+WM/UZav55l8acTo1H+1oprRmHXkP1/Htp9bgOHWWf/vrGk63tvMPezOrN+wLc+39SfOgEEIkuYz0ND76QT6Xlr9GS5uz169/4eN9vPCxdyA9u/Yznl37mde2zzzu1CJF7rSEECHrrvkynE2b0kzad6OG9GdFyRd4+JZ5sa5Kn0hoibhWUVFBbm4u2dnZKKWw2WyUlZV1exErKyvr6OPKzc2lrKwsYLnGxkaKioo6ytpsNkpKSro8dk1NDQUFBWRnZ5OdnU1BQQE1NTVBfR8lJSUopWhsbPTbV1dXh1Kqo54FBQXYbLaOc3ru66ruubm5VFWFNld0sMeqqKhAKYXdbqexsRGbzUZ2dnaP+9yCff+COZYIzcKJFu6+aDLDBvaLyPFVzxMa9ZmElohb7sCx2+3k5eWRn5+P3W6noqKCRYsWBXxNQUEBFRUVNDc3k5OTQ2NjIxUVFdhsNq8wqqurIzc3l5qaGiwWC/n5+TQ3N1NVVUVurv8MYUVFRRQVFVFXV4fVasVqtVJXV0dRUVGXoRiqkpISiouLAbBarZSXl1NQUNCx3x3k7rq7v8+SkhKvcsEI5Vh2u73j/c/Pzw9qXyjvX3fnEX3zRNGCjr//y8VTwnbcaIyFkj6tePbQsFjXIHQPHe3Ty6uqqmhsbCQ/P5/a2lqvfTabjcbGRux2O1ar1WtfXV0d1dXVFBYWAqZZadGiRTQ2NlJWVkZlpVluzX2hrK2t9bogFhQUUFdXR01NTccxampqqKmpwWq1Ultb23FOu91Obm4uFRUVLF26lJyc8MwuUFhYSE5ODlVVVeTk5FBaWtqxz263U1ZW5lcX9/dZV1dHRUWF12u6EuqxioqKWLZsWcBzBNoX6vvX3XlE31w9azT/tXQhLW3t3JYzgYusI/j6nxt6fmEckDstEZeGDx9OYWEh5eXlfvvcYWK32/32lZaWduwHsFgsvPrqq4AJQvfdlruZzvc3+PLycsrLy73C0B1w1dXVXtutVitPPvmkV5lIKykpCVgXi8VCdbVZI/XRRx+N6LGsVmuXQRJoX6jvX3fnEX2Tlqa45bxzWHr+JPqlp3HdvLGsKL6Iuy8KvOxJPJHQEnGpsLCQ6upqv9++Gxsbu33+yX0h9uRu/oPOoHMf131n5ea+s/E8r91u72g6C1RPgPr6+mC/tT6pr6/vsi7uZjeHwxHUwIVQj9VdU12gfaG+f9IkGF0XWkfw8C3z+cUd54V8jGg8KynNg/Gsj01sic7hcLBy5UoaGhqor68POIjBl29zoVtOTg51dXXY7XZycnKorq7uCKy6ujosFgt5eXkUFRWxZMkSLBYL0BlyeXl53Z4z0F1fJLgDpKeLQ3Nzc8f3EO5jnX/++V2W9d3Xl/evu/OIyLl5wXje+fQgK+v39Pq10Xi8X0JLxKW6urqOgQBWq5X8/HxKSkrIy8tjxYoVVFRU9On4VquVpqamjj4wd3jV1dVRVlbGq6++2us+KofD0WNQ9IU7ZCwWC8uWLeu27PDhw6N2rHCJ9PsngtcvPX4b4SS0RFzy7G/x7KMCWLFiRZevCzQ4Azr7sHyDKD8/36vpsLy8nKqqKoqKimhqauo4VnfNf+7mr1AvuMHepXmeo699PeE8Vnei8f6J8BsxKDOk151ubQ9zTfzFb5yKlOYOH9/Acu/rint0oCeHw9HRBOhuirLZbBQVFXmVs1qtHa/3PIe7bydQ86T7OaPumr88NTc3+21zD3oIRl5eXpd1cTgcHc8/RftY3Qnn+yei42uXB25m78l/v/ppmGviT0JLxCWLxYLdbvcLqKqqqo4LXaDBBhUVFV4PrLqHbwMdzWDu4KqpqfEb1OF+recdmXsEY1FRkdc57XZ7x9yAPY0edD8s7BuqVVVV3Q4s8f0e3XVZtGiR33vjrp9vGHclnMcK5jx9ef9EdA0dENrDx2dDmCKqtyS0RFxyP1zrviNyz9hQVlbWcfdVVlbmNXOD1WolJyeno6x7Jg33816ezWDuC6l79gn3n+6LtHsoNpgRboWFhdjtdrKzs8nNzSU3N7fjgeXi4uIeR7p5PvPlPk9ubi4lJSUBh/W7+5HcD+B6hmlpaSkOhwObzdZR9+zsbOrq6sjPz+9473oSzmP19L339f0T0fcF64hevyY9TWbEECmqvLycyspKrFYrNTU1HcGzfft2qqurO2bHaGgwD0SWlZVRVlZGQ0MDpaWlNDc309jYSE5ODuXl5X4PKJeWlnYMqW9ubu642ykuLqapqcmv76u6uprKysqO89rtdvLz8zu298RqtdLQ0NAx84Y7hGprayktLfWb9cJisVBaWorFYqGmpsbrTsj9/biPVV9f39G06ft9BvM+h+tY3enr+yei71d35fB44bl8t2AGmRkmKob0734YRFHuhIjXSyXaUsvxLi8vTwf7zM7mzZuZPXt2hGskRGqR/1eRsefIKSxZmcx78OUuyzQ9ckPId1tKqQatdY+dmzJ6UAghRI8mZGf1WCYKrYPSPCiEECI8ojEjhoSWEEKIoN171bSYnl9CSwghRNC+e82MmJ5fQksIIUTQotEE2B0JLSGEEL3y67v85+X845ejM8GxhFaMySMHQoSP/H+Kjhvmj2PGmMFe266cOToq55bQiqG0tDSczshPeyJEqnA6naSlyWUtGrIyY/PElPx0Y2jAgAGcOnUq1tUQImmcOnWKgQMHxroaKSFWXVsSWjE0ePBgHA6HNGkIEQZaaxwOB4MGDYp1VUQESWjFUHZ2Nm1tbezbt4+WlhYJLyFCoLWmpaWFffv20dbWRnZ2dqyrlBJiNYZQpnGKobS0NCZOnEhzczO7du2ira0t1lUSIiFlZGQwbNgwRo8eLX1aURKrPi0JrRjLyMhg9OjRjB4dnZE3QggRDj/64lwWPfEmABWF50btvBJaQgghes02ajCv/MflHDrREtLaW6GS0BJCCBGSGWOGMGPMkKieUxp/hRBCJAwJLSGEEAlDQksIIUTCkNASQgiRMCS0hBBCJAwJLSGEEAlDQksIIUTCUDLfXXgppQ4CO0N8+TDgaBReF2z5YMp1V6a7fSOBQ0HUIZ6E+vOJ1Xn6cpxIfKbCUaar/fJ5is55InmNmqy1HtXjkbTW8hUnX0BVNF4XbPlgynVXpod99bF+v6P184nVefpynEh8psJRpqv98nmKznmidY3q7kuaB+PL81F6XbDlgynXXZlQv594Fa3vJ1zn6ctxIvGZCkeZZPpMJdrnqS/HClsdpHlQxIRSql5rnRfreojkIJ+n1CF3WiJWqmJdAZFU5POUIiS0RExorYO6yCilqpVS+ZGuj0hsPX2elFLlSqkmpZRWSjUopQqjVTcRXhJaIm65wkouLqJPlFLVQClQDuQC9YD8MpSgpE9LxB3Xb8FPAhbXpgKtdV0MqyQSlFLKAhwBSjzvxpRStQBa64JY1U2ERu60RFgopYrD+JtrHbAIkAtKigrj52k40Ij5THlyuPaJBCOhJfrM9dtsJVDUQ7kcVx/VEddXtVIqx7ec1tqhtW6Uu6vUFM7Pk9barrXO1VrbfY6fj3+QiQQgoSXC4cmeCrh+a27A9FE1u74KAekUF74i9nlyvcXhz4oAAAcwSURBVG47YNdal4WnuiKaJLRESJRShUqpSqXUEXoYLOH6zbbW9c8CrbVNa22js/mv2lVGpKhIf56UUhZXP1YtsFJrnRvmb0FEiYSWCNUyoJjOwRLdKXb9WeHZ5Of6e5VPGZGaIvZ5UkpZMXdXwwGb1rokLDUWMSGhJULi6idQWmtFD30PdP4GvCLAvkqfMiIFRfjzVAvU+fZticSUEesKiJRgBdBaN/ru0Fo3KqU6yggRhKA/T64+LCtQGaCvyyGDfRKPhJaIhmACSUJLBKs3nyf3n+UByjRiHjYWCURCS0SLo4d9AfsyXM1FQvgK6vPkeqBY5iVMItKnJaKluw52GTkoeks+TylKQktEQzCd39JBLoIln6cUJqElosEOZgYD3x0e27pr7hHCk3yeUpiEloiGatefSwPsc2+rDLBPiEDk85TCJLRExHnMrl3q+dux6++lPmWE6JZ8nlKbhJaIFvcsBA1KqVrXlDoNPvuECJZ8nlKUDHkXUaG1rlJK2THPy7iXnGgEyuQBT9Fb8nlKXbIIpBBCiIQhzYNCCCEShoSWEEKIhCGhJYQQImFIaAkhhEgYElpCCCEShoSWEEKIhCGhJYQQImFIaAkRJkqpJqVUUjz4qJTKUUo1KKW0Uiqu5/Fz1bGh55IiGciMGEKIQKoxq/7W0Tk9khAxJ6ElhAjECjRqrQtiXREhPEnzoBAJTCkV9lV6PY7ZHO5jC9FXEloi7imlSl39Fjmur1rXv4+4/p7jU77SXT7AsfJd+8q7OH6hR19Ok0+5Sne/latMYTd1LvcpW95N2RylVLVP+eJu3ger6zVNwJEg3kL36wtd79cRj/eu0KdMpccx3e9Vt31avf359KY+PuXLPX42DUqp0h7qFdT72kXZJtfPO+y/FIg+0lrLl3zF9RdmjSTt8WcTps+lyfVvDVg9yle6tuUEOFa+a195gOO7X1fr+nIfuxrTr3PE4+/a9xwe9XG/9ohP2SbA0sX3pl1lPcvXdlE233XsJt8y3byH1d2cp9zn/Sn3qG8pkB/On09v6uMqa/HZ7/5ZeL7XDX14X/N9ytZ6HL8p1p9/+fL5vMW6AvIlXz19+VyASn32uS9+nhfeUEPLN4R8L2aWAOfwPI7nRbrQY7vnRbfSY7vV4yJv7aJ8aYB6HvF9H3p4/wq7OI/V4+Kc43N+v4t7GH8+va2P+732/RkUepy3wec4vXlf3dvyferuDsTCYN4H+YrOV8wrIF/y1dOXx0WxIcA+d7BUe2wL+U4rQHm/i2g3523yPbbHPncQaPeF1+OiGKie7gvvkWDehx7ev6ZuzuO+8NcGqGtvQyvYn0/Q9Qn0vvmUr/Q9dwjvqwZ0gLLulZD9jiNfsfuSPi2RSFYE2GYP4/EDDe1uBtBaN/bivH59QFprB2b4OJgLJ0Ae4AhwbLTWdtc5LAH6VXq7yKG1m/PUeNSlr4L9+fSmPu73qsb1HvoK1N/W2/e1EcDVp5bvUbZRa10R6DgidiS0RCIJZ0AFEpbRcq4LYyDui5/7QmzBXDx1oC+PcsN9jvNRsHVRSrmPUd9NMburLn3V488nhPrkeWwL9py9fV+LXMfJB9wDQ2qVUsUyECP+yHNaQsSAx8XQATzaQ/GoDD1XSlm6uJuJCdd71O33rrV2KKV8XwO9eF9dv2TYXHdZRZjwcn+VK6UWyd1W/JDQEqnG2nORvlFKWbu423IP/W50XWwdAFrrikjVRWttd13Uu2v+czfXRTywelsfpZTdY5sf3+H0fXlftdZ1uJpeXXeEZUAxZjCJrTfHEpEjzYMimfk2q4H5TTrSSnw3uO4A8jEXY/eFuB7TjBXoeTKLu5kqDPVx9+EEOo/7uajumuvCrTf1cb9XhV001fm91/TifXU989aklKr2LKe1tmut3ceO+C86IngSWiIZNbn+9LqguR4szfcvHnalng/Jui62r7r+6dlkVeb681WPvh63akzfTDV95z5PteeF33XOJ13/7PLh5wgIuj6uuz/3HdOrPuULMXdCXR2/x/fV9QuEFROKXp8Nj5+hNA3GEQktkYzcI9AK3b9FKzMLeCWdF7RIsWMucu7ZFdwPwuYAdZ5NVq5+kgrMRbTJVb5WKXUEE651WuuqvlbINSKvBtdzUK6ZIRpwPewMVLmaxqIihPo8inlPczzKux9grsNnMEYI76v7M1HrUdZ9fIB7wvfdi76S0BJJx/Xbcy7mgjYc8+wPQIErNMowz/KEWznmGa1czEVzOOZC2wiU6QCTz2qty4ACj7rmYS7CJYHKh0prXYS586zDhIV7Bvcij2awqOlNfbTWDo/31B1ew4EK13tUjs/Q9968r67PRJHr2MPpvBuvAmwyCCO+KG0eohNCCCHintxpCSGESBgSWkIIIRKGhJYQQoiEIaElhBAiYUhoCSGESBgSWkIIIRKGhJYQQoiEIaElhBAiYUhoCSGESBj/P/E6oye5yKweAAAAAElFTkSuQmCC\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"plt.loglog(nn,-rel_err*100,label='relative error')\n",
"plt.loglog(nn,-abs_err*100,label='absolute error')\n",
"plt.xlabel('number of nodes')\n",
"plt.ylabel('error (\\%)')\n",
"plt.legend()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Around 300 nodes, the relative error is subject to numerical noise**"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.6.5"
}
},
"nbformat": 4,
"nbformat_minor": 2
}