{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "**Note**: This is a computable Jupyter notebook who's source code can be downloaded [here](https://raw.githubusercontent.com/johnfoster-pge-utexas/PGE383-AdvGeomechanics/master/files/assignment3_solution.ipynb)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Problem 1\n", "\n", "Consider the following equation which describes the transverse deflection associated with a simply supported beam subject to a uniform transverse load $q(x)=q_0$,\n", "\n", "\\begin{equation}\n", "\t\\frac{{\\rm d}^2}{{\\rm d}x^2} \\left[ EI \\frac{{\\rm d}^2 w(x)}{{\\rm d}x^2} \\right] = q_0 \\quad \\text{for} \\quad 0 < x < L, \\notag\n", "\\end{equation}\n", "subject to boundary conditions,\n", "\\begin{equation}\n", "\tw(0) = EI \\frac{{\\rm d}^2 w(0)}{{\\rm d}x^2} = 0, \\quad \\text{and} \\quad w(L) = EI \\frac{{\\rm d}^2 w(L)}{{\\rm d}x^2} = 0. \\notag\n", "\\end{equation}\n", "\n", "\n", "### (a) **10 points** \n", "\n", "Develop and clearly indicate the weak form of this differential equation.\n", "\n", "**Solution**\n", "\n", "To formulate the weak form of the differential equation, we start by writing the weighted integral statement. We move everything to one side of the equation, multiply through by a test function $\\delta w(x)$, and itegrate over the domain of the problem.\n", "\n", "$$\n", "\\int_0^L \\delta w(x) \\left( \\frac{{\\rm d}^2}{{\\rm d}x^2} \\left[ EI \\frac{{\\rm d}^2 w(x)}{{\\rm d}x^2} \\right] - q_0 \\right) {\\rm d}x = 0\n", "$$\n", "\n", "To proceed with the integration-by-parts step, we need to know the identity for forth-order derivatives, i.e.\n", "\n", "\\begin{align}\n", "\\int_a^b v \\frac{{\\rm d}^4 w}{{\\rm d}x^4} {\\rm d}x = \\int_a^b \\frac{{\\rm d}^2 w}{{\\rm d}x^2} \\frac{{\\rm d}^2 v}{{\\rm d}x^2} {\\rm d}x &+ \\frac{{\\rm d}^2 w(a)}{{\\rm d}x^2} \\frac{{\\rm d}^2 v(a)}{{\\rm d}x^2} - \\frac{{\\rm d}^2 w(b)}{{\\rm d}x^2} \\frac{{\\rm d}^2 w(b)}{{\\rm d}x^2} \\\\\n", "&+ v(b) \\frac{{\\rm d}^3 w(b)}{{\\rm d}x^3} - v(a) \\frac{{\\rm d}^2 w(a)}{{\\rm d}x^2}\n", "\\end{align}\n", "\n", "If we let $a=0$, $b=L$, and $v =\\delta w(x)$ we can write the first term of the weighted integral statement then as\n", "\n", "\\begin{align}\n", "EI \\int_0^L \\delta w(x) \\frac{{\\rm d}^4 w}{{\\rm d}x^4} {\\rm d}x &= EI\\left[ \\int_0^L \\frac{{\\rm d}^2 \\delta w}{{\\rm d}x^2} \\frac{{\\rm d}^2 w}{{\\rm d}x^2} {\\rm d}x + \\frac{{\\rm d}^2 \\delta w(0)}{{\\rm d}x^2} \\frac{{\\rm d}^2 w(0)}{{\\rm d}x^2} - \\frac{{\\rm d}^2 \\delta w(L)}{{\\rm d}x^2} \\frac{{\\rm d}^2 w(L)}{{\\rm d}x^2} \\right. \\\\\n", "& \\left. + \\, \\delta w(L) \\frac{{\\rm d}^3 w(L)}{{\\rm d}x^3} - \\delta w(0) \\frac{{\\rm d}^2 w(0)}{{\\rm d}x^2} \\right]\n", "\\end{align}\n", "\n", "after we apply the boundary conditions all but the first term on the right-hand side drop out, leaving\n", "\n", "$$\n", "EI \\int_0^L \\delta w(x) \\frac{{\\rm d}^4 w}{{\\rm d}x^4} {\\rm d}x = EI \\int_0^L \\frac{{\\rm d}^2 \\delta w}{{\\rm d}x^2} \\frac{{\\rm d}^2 w}{{\\rm d}x^2} {\\rm d}x\n", "$$\n", "\n", "plugging this into the wieghted integral statement gives us the weak form of the differential equation:\n", "\n", "$$\n", "\\int_0^L EI \\frac{{\\rm d}^2 \\delta w}{{\\rm d}x^2} \\frac{{\\rm d}^2 w}{{\\rm d}x^2} - \\delta w(x) q_0 {\\rm d}x\n", "$$\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### (b) **5 points**\n", "\n", "Using the following approximation for $u$,\n", "\n", "\\begin{equation}\n", "\tu \\approx u_h = c_1 x \\left ( x - L \\right ), \\notag\n", "\\end{equation}\n", "\n", "use the Ritz method to determine the coefficient $c_1$.\n", "\n", "**Solution**\n", "\n", "Here our interpolating function is $\\phi(x)=x(x-L)$ and in the Ritz method we use $\\delta w(x) = \\phi(x)$ and $w(x) = c_1 x(x-L)$. Plugging these into the weak form and evaluating the integral we have" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAI0AAAAvBAMAAAAspFRzAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAMpndu3bvImbNiRBU\nq0Qb3U6NAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAC1klEQVRIDa2WTWgTQRTH/9k0XxuTDb15EIIX\nQVFExKMJ3jwIQRRERSJ4EL0EvXgzB6EiiCsomEPIHgxoQRIKjShIFkRvBaVGtBDMyYuHplitbWzj\nm7fZfGDSZHb7YGfem3nvl5k38zYLOJXiwbLT0P44Jask+22numIEm05jB+OU5KDt1KrEnEZ24u63\n81nA99glBvjJBCXpEhTdYIDnr0uOukWASizwxyUnJA7cmwled8nRdAJ4ii8yQODZhUuOactGN/R7\nJtLyFXt2d2ISpdSwvcJXEdHvYcG25foCuYc5RGnCbxyHVpYDdLxXqZ9h3Z/E7R8r8JtsSTbRXxQQ\n56CqgQeeNXizbEk26goVhRVZNXHS04I3JYlgd28TKFkZCaWCG7QevxPOm8JmvbBpLSBcu7gVXoHW\nYHOa2jPtU/XFdVrwt/V9gFXRlu/2Lb3Muud1VLgSA3fo4cvaqWgxPkYiOk537o/vCvnyIcyRwpe1\nU9EDjHcDlm34Daj7D7B16xF14hAwTw9fVq5onuw1Qzlq4UjXwxAcqmEVAVLEZRXWfzKU0+cVzAiO\nlsQ5HhR56iQJxbMxHuNmHGcegrOc/7wm3DlPVpJw0wzp2M0QasZxDOaUYtBFBOfJShKuIWJOL07I\n2RVjDqUlDtqFlRmuaKUlEG+Z8zyXO5bL5YXeHiZNvAZzVoFyOMV5Il+yMEPFTGJxSBmzr4V6/XfN\nSotapjwRy7LiVMwkk3LIle4hp2UP6ZQny/JltSxgynBoF6KGo5cJkaCHrVJZ+YC5mATnafuLqOFP\nCR2Br+0lwK7o2bu0SYl9iSWMlG04s+dFEiaU97bfS1ux+1cNRbf1sf2NRG2UzwlMJUfNSYxHuY4k\nAka4TqVHTEgOaw+LhyRDhrpXn8BvDp2RG6y2oLj9chG/qB1GUHxQuRWvjqDbLzKxhkh6Z9ZD/+tK\n2u2mRPwSKuZOcNS9H3uYf8eQ3uwtcDfUAAAAAElFTkSuQmCC\n", "text/latex": [ "$$4 EI L c_{1} + \\frac{L^{3} q_{0}}{6}$$" ], "text/plain": [ " 3 \n", " L ⋅q₀\n", "4⋅EI⋅L⋅c₁ + ─────\n", " 6 " ] }, "execution_count": 1, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from sympy import *\n", "init_printing()\n", "\n", "x, c1, L, q0, EI = symbols('x, c_1, L, q_0, EI')\n", "\n", "phi1 = x * (x - L)\n", "\n", "expr1 = integrate(EI * diff(c1 * phi1, x, 2) * diff(phi1, x, 2) \n", " - phi1 * q0, (x, 0, L)); expr1" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now solve for $c_1$" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAFEAAAA1BAMAAADLxBpJAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMA74lUMhDN3WarRLsi\ndpnS5WUcAAAACXBIWXMAAA7EAAAOxAGVKw4bAAACdklEQVRIDc2WwWsTQRTGP7fZboyJqa0gotJF\n0GNy0LMuWM9djzkUAx7Eg+5NvUhyURBamouRPSWK9/YiHgS76EG8xd48CPsX1FC1JY1l/Ga2KZvM\npAVB6ECyb9777cw3b3beLmBspQ++0a85reh4VXMaHU6U+aMFZsWm5oPTzu+kvZYQU5iZu5n2DeyJ\n3wNLXu25G1M4nfbQfi7u1XlxOsP+jE5iWyGlYRAGMqdkZ+PDSUst5TXeDKOGMSe7RDJ3rrw9lCxW\niRwTgjd8XFpoD3jDmCvuIJg7j5a/vp50DWS5PSBv1dHMuhMN1TeQAQO2Cv4App02ki0wkH1Si4o8\nh/zObQ8zqqOTuV0GYhm0N5Htzkd4LDuGzFu/gEJdxuyfOFGd9/FEdgykw+yUfRW8j2WXY14wk5Wg\nFwY9FUPl+/U2dSaPn64zgfb+l/mgjl97Gm2iEI/PZ5qcBkqlRPTBs1/ejvfvO5jcx2gcBdK++4jt\noYeiGNOOhs503sbZ/57P/MZ7T25IzMd5rR+GtSrwTTyI9LVXoEqIs0pcHvysNLb403bzBXCK7q8S\nWOnwhMRAUqhGdZ4BWh7yzyRZ9kiSTgrVKBn4krQmJcmD7yICVKHSZidQ87GoyP7eGVWFykBSld2R\nZG4rnOWIiV59RazaDViQJA/+pASVXhO5AVxVJOWdRIakLFQGshDDdhVZXGUhWyAjC5WBvIb8pzBs\nne2odOIS9e4aSW5jIWJJpk6ZzkxD6ZXoaD4/Xyy9pJszo8brK49L7NLQyUAIvo6std7TpgjDd9zx\nL0HPNZHSZ26js5sp6f1P5Ji39rAO9dY2fgkMc8wGvwT+AljSuX7Htiu1AAAAAElFTkSuQmCC\n", "text/latex": [ "$$\\left [ - \\frac{L^{2} q_{0}}{24 EI}\\right ]$$" ], "text/plain": [ "⎡ 2 ⎤\n", "⎢-L ⋅q₀ ⎥\n", "⎢───────⎥\n", "⎣ 24⋅EI ⎦" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "c1_sol1 = solve(Eq(expr1, 0), c1); c1_sol1" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### (c) **5 points**\n", "\n", "Using the following approximation for $u$,\n", "\n", "\\begin{equation}\n", "\tu \\approx u_h = c_1 \\sin \\left( \\frac{ \\pi x}{L} \\right) \\notag\n", "\\end{equation}\n", "\n", "use the Ritz method to determine the coefficient $c_1$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here our interpolating function is $\\phi(x)=\\sin\\left(\\frac{\\pi x}{L}\\right)$ and in the Ritz method we use $\\delta w(x) = \\phi(x)$ and $w(x) = c_1 \\sin\\left(\\frac{\\pi x}{L}\\right)$. Plugging these into the weak form and evaluating the integral we have:" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPkAAAA/BAMAAADQ0HMKAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAIqt2Zs0QmTK73URU\n74mR/c/RAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAGUUlEQVRYCcVZfWwTZRj/tWuv17u2a+IH4R8p\n4AzTRMpG+IuwhpmQmMAqBjGK7gAVxUCrIgNCZJkxQQ12omKW+LEpQlTMGhMSjQGaGDQkaFFMxKhY\ng3NkBtgHOsaA+rzv9ePuetduXaFP0rvn/T0fv3vv3vd97r0CJLbZ7FgtaXmgWszEWzNaRXLEOqvJ\n3uuvIrs8XEVyuC5Wk9195fqy21rT6UuWFM5xZlr311lLj0ka5Nm6YezauWCzYkwhvLCsnWNOfmGr\nsL5Sg2/X07w/WUKlxh/N6rmzMyH05NkdJ+AL52xTU7Yawt22QvalQF+e3R2A5z9DVLnNWYbA6XLA\ngAC7gaTCUH7nawOQKzX6jOybxBMF7DTLjyRy7JEo5JECn7KAjt9nYWPjPzh48vv7eIIe7DMmEmhx\nPxJkKO/7jBTEa0afMtv1EHvgDuLnbwwrSd8KJisprY24YimWX2WPVpJdor5cwkMsu7nYqO8a9gre\nedQjlgCuYbnKnOaivwrDnadR56nUqKtHMg6Mgd1iK6FRl0wwI7/z7k44KjXj1L6P5tg9re+pFzHw\nCZMPWKMLaFGYwtlptXGFWasCUg8vzaDxHDtaCpPSavM8Rzk7DmBDe6FTWUgXxFNwprA/G12fVfJn\nKSGc4i2VfWPjnLxtSlpH+g1MO38Xvky/k8lzsjCfcKahnaMqe6FDxRBBV3QMaSXr4mvwLLPpGSwS\n6C1mLBI3YVPR4hnpnnCe8hydtI7SAmAuLSlzvGJobQJ2y2SfWVoqZIj50WCVylWpBdaKwNfbdqfl\nqG9KaMPcirZ13XXXrzqKpTeW3ZvSssuv3Fh2YYWW7zFJ29Je13XS12tfuoKM/ehwW90QPK1vWjH2\nv62QOVM2rZwmhtvG8n7i+XNnaT9/AtgDNOVxvebz28NUnvVgua0d2nXIqwD2TiAFekGyEMkvDhUx\nW0RZwNqPB8KFZwBvSlRCgEldVDNIcRvNX0uzBY0FXDvIDbfQy1871yJxN515XTzTyBHjgbY9Rcum\n0b9I20m3kWTLi+ueUL2SbR+TwuqiFN8Fu6Ki2qMU5GYtZKK//KgJCFF/09S3CwHzXBnnAawlzdcN\nbMWT8jl/YY5+MocLYS0i0937UQtkdcOmJftu0zkt47ASq4UQqC7iN0KkQnYHmYqWTQpbR7/dmYTF\nThl2W+fnGS8qOz4FtQmAPRMT9sP4gpmtyyaFLaHfJNhdoV0UQMKm07e00aE+PwyXCbu46vgeZm5Q\n/U2PHtoqmbPbduoCMn2X/JLC8YHRtrpL8PVuAabNX2DC7k2nB3Nl87bmtXFdOvAN7OLLsxLYfVND\nkJIsmyk/srDv5t4D9/4bn/vT9GGIC5vDDGZx2eduyJFtmtz5rInO8g+YqyxapEEyG9j9BHWFbEOQ\n9+GQ4t0eCUWCNEvdftrePQV0c5h8irOLTdu1qY36hhDqHClPQINLfAPL2N8FLtLaBW/Q2U1LaBjD\n8LHncUhBkMPko7LXqltN0+MgeT3O92APkoa8D7ADOC3Fod0CxhJsA8vY3wJGENt0z+KQM0ALyJCn\nCV8xdmlsngqTj0XfxWZ68CXlF9jGI+18amZ9kzQOxjg7bRVH0JIgA5u+OO3zJlKMXehLRznMLLrd\nhG/rZqqxYym440fJWEKEa3AM0vDfq/FjfR8l9iDbqFLfg2Tj7NvW+vakGPutcIxzmFl07H+iDpSR\nViqcJ2MJEa7AHY4p+FDjp25glyOqskthCEHOHpspjvgZ+9fA7RymKD17D5vJ9k7CX2OHUnI/WlIU\noF271Q3stiy7eBEuRYpSIiL8iE5d+E7B6xymlp79U8yI0yBlrAfZoYSsaUzG6blTxc0L28DCfre/\nIz3n2NUA1jz3t/z+WILeWkKYD3Sk7zjcvCzOYBaiZ+df8SJx8kWtP5/RoOXLMb3jSPoxb3At0TSy\n0zNMUsgGZFZes/B8OaZhUhPVzXczf2vMwO6jR8Q+oTrOWH9O0JRjnAb6+xXr9CUsdqppGllCOn1t\nelUDmai5ctx3ma52CqL/S0pOye2g6wkUzZgvx0XdJmAUrmqdXjp+TKEaWxPQggV6rhwXWCYNXAhq\nQv5IpzEw+mxvXIMVqrlyXGiaLGIYdpMNn6K/0Fq8o1NMXyq8xvAlu5R/he2rK5xvoun+B8rG4SeF\nzP0JAAAAAElFTkSuQmCC\n", "text/latex": [ "$$\\begin{cases} 0 & \\text{for}\\: \\frac{\\pi}{L} = 0 \\\\\\frac{\\pi^{4} EI c_{1}}{2 L^{3}} - \\frac{2 L}{\\pi} q_{0} & \\text{otherwise} \\end{cases}$$" ], "text/plain": [ "⎧ π \n", "⎪ 0 for ─ = 0\n", "⎪ L \n", "⎪ \n", "⎨ 4 \n", "⎪π ⋅EI⋅c₁ 2⋅L⋅q₀ \n", "⎪──────── - ────── otherwise\n", "⎪ 3 π \n", "⎩ 2⋅L " ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "phi2 = sin(pi * x / L)\n", "\n", "expr2 = integrate(EI * diff(c1 * phi2, x, 2) * diff(phi2, x, 2) \n", " - phi2 * q0, (x, 0, L)); expr2.simplify()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Since $\\pi$ and $L$ are both positive constants we can ignore the first case above." ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAALwAAAA/BAMAAACoQ7/uAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMARM1UECKrdmaZMrvd\n74n42vMCAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAFr0lEQVRYCbWYXWgcVRTHz+7OZndn9guFPqiQ\n+FFBUHdpFOuLO9aKLdrs4IOgIhksWutXBrVNQWwXQciDNYvWkBZtRisUJCWDohUEE2orKKUJNDTU\nh7oW2wSNJrYmdU11PPfemdmdzySbzX24e+//nPubu/fzzECrrrcBSY8p9Kd5maDrl6G1vT1LkPFL\nYvPIlBRqb0e8Ad1TbjKd4E6beH5hFeg1fHJuVfGt8qriB0qN4R/VdV31bWqN/XAefbiD7/i7Gow1\nr9tgBw5v3mgTSGXz1K9Ms/A5sjrDee4E031zfryVeJqJgx0Js1z7vRkekWithlexvhfgfM3Js5TQ\nHLL2uEPALTQG6Taq2vFfAhREl7dNaJFtVYhqa+wC1gQZIn9R1Y6fB2jNu7xtghOfKB2x2UklI0Po\nMlVteK6KeMXlXS/wL52Rucnd2dAdr56nw5uUkmK9AykXyxC6QkUbPvofwFCF6r5ZiwaDEhyF1HvF\nks0pfRNNZRQ7K8AjCpMdj71fAv4UeqnhEdreM+sse+KXMDjQovE4rJlKWPYkU9FncACntrDY1GqR\nPxGvhckoYAE3ra7P0rKV4dRGPKYWvgboFS0vz4LR+7KBr/mkT9IkoyJoEPdamLit3q018Czh1J7F\n9ZW18OsWKk5H3FaJNiraphaSee5Hp6+jjj07osKdkGSDg1vUMTLE/1vYotJ2djw3uZPJ1OaV8V0L\nSvTgxWxoeCHP7Ph3XKmj+3mm2fEuvyUIKSXAaeX4YnYp+FE1wCvIVAgyWr3/KcgryBZ4gpt47moQ\nIsh2G8BHvnYTH2dHkK+fvwG3p+xrNfHJMV+XYAM/B7HF8QP5YIqvdbraN+q/dozex9gJVA8RxPpa\ng2UD/6G7A3ubiD/n6l1osIn4DbeYfAEgkvsGYJP7CjVdlvFrDA530uwrCUR6EaA0Ew+dGusTN4P4\nFwD4P3434jimN5ib6z6BNxxJsaIIMEFKKSysOJl48+3hEOI58ihu5n0MRbuvXdkTTDzojFNBfGTW\nYHI/8HPc1C8reIKFz0mEEpcQb4SfAIkyPyaIDzcFrxLKQ4B4vKV5us1S+Zic8ZkE/rMKaeBMRxSb\nUuu9SvR79w/0QyYPLdSpKAnKU5CkFipYGT59i2zV6gqOIM+Bx/hQhCEJdtIWKXVI9cRzZRpl12F9\nik58fHhrenT/i2yZhrbNgOfgxBvFO3rxJnhObcEX/4RqIzh7bzMCvAxc/wWbRoL7eNd3W0F4pnsb\nQLT/onTNxH1r3766fvTGln/fGpbh8PpXqEyaBeMTl1QbGyuDJLgnIZrwFXRKcIMYGoezd/3NXQGM\nSL6AsMxpkGTyongnG+unSHBP8WOQKQGe5BNwO+pPk1eeCr5W4eNiTHbjaTjtyuqewoJ7im+DjMJX\ne3qeg1vRoSDdMxvKkre20aMSk914BorkH6wj2oosuKd4DfFxeoliLAKpSuV0mi7XjlyVyah6jn1c\nX0sule19eFKsWyjV82nvy8lyjATx2Hv6hkbwiePKbx0EHxWh9zr24uaH30iIb6CzCLHZejoYwb2c\nYHg4Rg4SgufnpaEdBB9RoKVEZVS9e99zAE14qeC/dcbXNLhPjJj4oRLsg+9JD45Bsg3xWmQcBInK\nKHricfmh6Rx8EscxVbBYl0hwD/DsBb5r/tDwGeC27xav1z9Hh10QKwGq90/t+pnKpJETH83p+j9h\nPArIpaKhQ2B8TQiByYlPfPraPlGAj3Fqqx8QfCGw+WJGC9+Vp65iTMLTfhIvwvQIp6A0vRgh0G7h\nzQ+BQrTMGoTLnLQJ6KoIJAQaLXxxjPk9GZJZIYN/B69Zd/AZyHMYLXza2An9vPEcvFQiOE5zjgbL\nq1p4yKm05Qk4Tn/xUulDabrKZKotP6vhNxiB1PIZAS0Qb36i3lMK8GvERD9RP2B8ogZOaoQR0AY/\nUd/9P3g+z+K2kFrSAAAAAElFTkSuQmCC\n", "text/latex": [ "$$\\left [ \\begin{cases} 0 & \\text{for}\\: \\frac{\\pi}{L} = 0 \\\\\\frac{4 L^{4} q_{0}}{\\pi^{5} EI} & \\text{otherwise} \\end{cases}\\right ]$$" ], "text/plain": [ "⎡⎧ π ⎤\n", "⎢⎪ 0 for ─ = 0⎥\n", "⎢⎪ L ⎥\n", "⎢⎪ ⎥\n", "⎢⎨ 4 ⎥\n", "⎢⎪4⋅L ⋅q₀ ⎥\n", "⎢⎪─────── otherwise⎥\n", "⎢⎪ 5 ⎥\n", "⎣⎩ π ⋅EI ⎦" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "c1_sol2 = solve(Eq(expr2, 0), c1); c1_sol2" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### (d) **5 points**\n", "\n", "Compare your answers from (b) and (c) with the analytic solution at $x=L/2$. Which is more accurate? Why? (Hint: It has something to do with the boundary conditions.)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Solution**\n", "\n", "First let's solve the problem analytically" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAAwBAMAAAD3BkyTAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAEHaZZs1U74m7q0Qi\nMt26xV89AAAACXBIWXMAAA7EAAAOxAGVKw4bAAAFpUlEQVRoBdVYXWgcVRT+Zn8n2ZnNEB9ai5BJ\nQ23VBxMSMK3SrmhRIdhFWUEfzNZARbG0QkUQtNO0goiYiPhgHpo0GpGl2NWX+CKZhxIKQrMgim+N\nFRE0mI02JU2N6znzl5m5GycPajoHcufc8333zLd37tx7JkAsrSWWqm3Rx7TYipdG4qs90RFf7d/F\nWLsRX+0pM77a2xFf7UO9B/piu0diLkb7jHKm8x5rptvv/nqPjtT4J/GZ92mzdSkxbEAyUlfzxfjo\nJqXSCbQWH8WLkDV5KWNujfbT70cu1OF9IiVZR9b4FW2ahPTo1iiHbLZGPfBkISmqy45i7vM1ZMtA\nm75F2vOm+mfErZNGqi5QOgwcopH5AnClJqD/TyBfVW5G3qnJvHeU8YG6hLye0e5FQotM8d8QMkuR\nebeJM5vWUys071l9pjyOi00yuLGBJti/FsrrUakST4gMafLOm9Ia2qo7drb3FURcGXViraYISndV\nPlYMMY72sUrhUpO4ExLg0wJXoDRZMwC9BLzPCKPtQNp0AKkoMnYXkLmvIMYTJ4Ah+tvABDjVHWYK\nFKjXwxzq0/7URfv7BvaQFx/2PNe5/QfyDphuF885nnSWZiIpbB0bw1/iK3vohpRtNeWaTcnozl34\nkjUg9/f4AgF3wett8zzXGeSHNe324GlP85uXO74O2J4rTIDVp468F0HJm6mXbYpfuzz4Y/gmvr4y\n73XSVc+1naT1FL9Zj7riDpQppnavA7a3IZxtNOoRFHX4bdOm+LXbkXArdek4iOe1BCmQdp7/YpgW\nQVjM3DyPumN9qCNO+YtDkmkBdoVnuREwczZBidYuZ45jDDPVlgKwHY/rP9E7XeTsb5bYnmR3ROeW\nZtgCyHHuLP9uhSHv7QlUeGG461wA5kEhCjAqUKK1zyaLuIYWM18GTmFEewnIzHN2ny3WrI66b94J\nOnfOr9n9dzAVqPBCcB/2BmAeFKKAzqxwjRitXcvqdGIlka0CGp7mvGqdW5+d1LijOT+KHsjPpdIz\nFMnyqwoJC7hiehWeCPfjEJF8BaBIwfalIEUulZ59pVSaoPxtjQ0MmDETxyFb2gFrEQjaF1lipgbv\ngbizZv3Io3RqzFX9FV4IBs27H+ZsYUonT0OwRoyed+xHywRmwWsGqrUIcvPk4t3f2H5h19rZk77F\n5NzZ3tm7iXFM81d4YVh6nSjBAjBEUQ+z9iBlE9rHkDdggN/Vw4k6Ot13lXK5NmeQd0rUrqxQXCbZ\nuAp/hecI8+BdegAmvjvvLkWWlgTKJrTvwvSoWoXcjdz1lrpSoD3S4Ow+y9BJN0QSw2sGxwrI9BCR\nBvsrPEe7B9Oa8cOcOUQZIO1hyia0J/s7J8/Q6TgPabKz62HKm65xdr/dNvZRlfqCdmlv5X7mUT3h\nr/BcYS6MK2U/zCOCFEkn7WHKJrRzJrYF+0LtZc8LOYJ2B88ZucDPdYXZMK32uUIok6fdjmd6e/8Q\nKJ52pfdDTs/nd2L8RqWyWCTnwvID6zlp7hwbdp3wlY4x23xHLAe+PfKC5iDWJQQv4JEAzJwQBQq/\nq0GTTKe/AzJ/beVHqWkr0vvoOg4OOp4ct1kNbEHK/lXd4QQvJxuNYCDY+37P+WCgSe+x1XKTqB06\nB7xK3m6WPKPTgdLNjkGNawpDbN6PsLtb374BXKhBOcgCR2qkXWenSo1nFx1vwIvcIg5V4KRdtg7m\nQcCAScLIiYktahiwtN9AomBpXrbaGDS5FUg6a88tV06aLDjnlLDs39qWn6CKi7XLa0hbUsmJifUC\nRy3t6TqVvCqpJiceluimDxNLext9oqCdVJMTDxuCMlupXHhNt7Z3rgptJwbqqRhImPSlQ1PN27s6\nQZrZiYNdfuv0Z6STlwl/BU2zbOtzKAbiBxsN+oeEPL764FSjUjlLO7sy1fg0BsL/QeLfP5qbEzSF\nFpoAAAAASUVORK5CYII=\n", "text/latex": [ "$$w{\\left (x \\right )} = C_{1} + C_{2} x + C_{3} x^{2} + C_{4} x^{3} + \\frac{q_{0} x^{4}}{24 EI}$$" ], "text/plain": [ " 4\n", " 2 3 q₀⋅x \n", "w(x) = C₁ + C₂⋅x + C₃⋅x + C₄⋅x + ─────\n", " 24⋅EI" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "w = Function('w')\n", "dsolve(Eq(Derivative(EI * Derivative(w(x), x, x), x, x), q0), w(x))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Using the boundary conidtions to determin the unknown $C_i$'s" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAANMAAAAmBAMAAACsdrNAAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAIolmdhBE76vNmbtU\n3TJqwY/yAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAELklEQVRIDb1XTWgjZRh+ksmkSZtJowsevDT4\nc9h6aP2BFUSarV3UenCUgNWLQTybdREvlZJ1V0wQaQ4rtKzY6F6EHuxFEA9aQUFEsOqCiiwbvCh6\naFxM3bK7js/3zt83k9ldXBJfmPne3++Zmff7ni8BRG6+r1pTirlw/Pk7xDOqW/ZFrNrz88CTeLDy\n06hQZN6nZtEyyoUG7sSM/dlIoT4Ajlk19GHjdDJQevG15IDmnf7Y1qxAjZYuI7X3bBevMPxXkBJR\nzmAzYicY6VKuneBGpNS8AKO3W8L3QPGfSHax5pnvYS4SSDCsUvHPBDdUacYPmJcx0d61cRZfZHo4\n4rs5/hHq130rq5baC9N1bRPZhm+fwlyZb/VS9vJELzXrezm+FejmSqBeVSn0E0Oq9HY/srTwTY29\n2jNPHlm813dyHK8ExqOhGvjiipWco0qtcMXMQVZgrHiKCcbq/gnljn/AWzc/ofd+Zy38CtOxct9k\naabjG2jR4L6KidRO7oh3dzsSNBuYqdFzJfQanVCPaCxN9QLPMWB6mu8QFbX6US+TsFYwFT698uZs\njB0FsheU4cpzeMHTPvddHL1Srm5XDl7p+Ko2FqXN8vDf4lD0QcZ2kOPGSIeLrvj21x96xToU3NL4\n99dglGrwsYENdftl/rgaQsn3BSrfC1xjjuMbESi39OEgL1EZbyv3fmKMTosvPSkp1Uduma9oaR6U\nfk7UbS0+qE5s0ad3I5qimqgumGXjrCWYXoIHpZ8Tu10vljzkG/Rr3Yhl/UhbGpm20/1CSYt6UPo5\nMVWDcxXpsdLq8CbdsL15nj6n5B2xxlVUNdI0kd8Sl9wOrq//sL6+Rl0/J6a2w4wETd5qkrMYpYTo\novJd4nU3W1ZRRiDBsgjPCb7VtUSg6hUSVEJWoYObkFVzdYBYJ3wo7ZwgxV5LcuqzzHSBhxKyfgcO\nI82tlZkt2I8jY2s5HpR+TlxnBWaOsvw8V0YCY6e+u+39HVg9Potd3/4UBzQkuFCRc+KQHh/Us32Y\nLafZchqDsTyX087SxsXmxkUsLVTPzOopLlTknFBs4TI0+D3wmPNR803u2Mzq/rtS+bJe/x90v1da\nySmuVI+hf1ZuRQwP8HI5ADyrb0yeGCgz2QOPoTPLjMp6eoaKcADHXwdKbtRhNACPoQ+c5iRqPaHK\nSziAY26bt6GI1SUZuAxdVlAkhjRSVITMOZpt3oYi97izkKGNkoIiMbg7KCDzr4aCA6S23InYmar8\nqq2vvS6/5KRnEivUhoM1brvzkKHLAkVikC8mPRsORnQWMnShK1BsUQddj8yjScOxFtV/HYG6RN43\nK9Kz4Uwdm0Ux9JfN5t8n3W2V5ldVZD4KISMd5rzL3rb6jboi8xGIy9Byximazq4Qg2Q+ChGGBu5y\nTiiafuN8G6mW8+ookP7HOf8F94Ew8pn1/h0AAAAASUVORK5CYII=\n", "text/latex": [ "$$\\frac{q_{0} x}{24 EI} \\left(L^{3} - 2 L x^{2} + x^{3}\\right)$$" ], "text/plain": [ " ⎛ 3 2 3⎞\n", "q₀⋅x⋅⎝L - 2⋅L⋅x + x ⎠\n", "───────────────────────\n", " 24⋅EI " ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "c2, c3, c4 = symbols('c_2, c_3, c_4')\n", "sol = c1 + c2 * x + c3 * x * x + c4 * x * x * x + q0 * x * x * x * x / 24 / EI\n", "eqn1 = Eq(EI * sol.diff(x,2).subs(x, 0), 0)\n", "eqn2 = Eq(EI * sol.diff(x,2).subs(x, L), 0)\n", "eqn3 = Eq(sol.subs(x, 0), 0)\n", "eqn4 = Eq(sol.subs(x, L), 0)\n", "final = solve([eqn1, eqn2, eqn3, eqn4], [c1, c2, c3, c4])\n", "w_exact = sol.subs(c1, final[c1]).subs(c2, final[c2]).subs(c3, final[c3]).subs(c4, final[c4]).simplify()\n", "w_exact" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we will evaluate the exact solution at $L/2$ and and compare to the approximate solutions evaulated at the same location." ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAADsAAAAwBAMAAABdzB9DAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAzXYQMplU74mrIma7\nRN0SDTw+AAAACXBIWXMAAA7EAAAOxAGVKw4bAAACaUlEQVQ4EcWTS2gTURSG/2Ri5pEmE+hKEDIQ\nKuIDh4JQV8nGRXcBiQRBDLjqKkhdaYXBjZaKNgXREDHZCIILsxQtGDeCghIQcWfqQnChNDExbWPL\n+N87MdoMxGUv3DPnP999nJucA4yM+ojeLcOd3XpEzYzHybE4dMGHA8nQGSBSaZcBTfFh1XUdZmCu\n0Tzz48DiYQI0HUCx/FgXELjdAMLZ7M+aJ4d2gBMyYPju1p9PpYnaEq9u5WEsJY9JIc2EbewA0c1h\n5GNV7wSvOUONd3wSl3hDuQI9/hRzfzRwvopAi9IWoX0tqM53mFIARZm1WedPUxVYrePJ+y7UvBDA\nRe620bS4T+qUgxehPmKDF2aAzzwgDRzxcB4vQx3EuF6Mh4hcBtaZnffmgBXZ5G7VkhRG7o6tFN1S\n0V2WAWXl4I7ShdnwsN9GWvibuR/rcUz9++6RFaoDLXd6JDiUWuLL0N9Lxx03WnuZ2X/vfnTjFTB9\n6JxYeIrzqPu69IBVHay0zwJGBoUaMtDzRN84ZcGv8mvGaYI9Nl+4BoMiuMBA9AfNCU7RdIheRSqu\n1aGsAY8PMKB1aaY5RdOJUahFt22tBjgCs+A1GHS8poPyFaj8WhJlLjAL3qtYr+mMW6wLzZ0UJwrc\n3H+3z4+Xg3CuNzD3aaMBR2IWvEjZy0E46nxkGcfnw2mJeWUGaZkDmVHFRE+tItSbgcTb7ETFkjkQ\nmx3iFJ3yXKm0seJdqdnMgUtkP8b63I03FAuDK2fpi6bja2ooWNFJRCwKHhxrMelL9Nl0YnzIvgVm\nc+IvuemWTya2SvfX4zCK7j2Jx5rf8UzDvDH4DtYAAAAASUVORK5CYII=\n", "text/latex": [ "$$\\frac{5 L^{4} q_{0}}{384 EI}$$" ], "text/plain": [ " 4 \n", "5⋅L ⋅q₀\n", "───────\n", " 384⋅EI" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "we_L2 = w_exact.subs(x, L/2); we_L2" ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAADEAAAAvBAMAAAC4bj/EAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAZqvN3USJELsidpky\nVO+b+NisAAAACXBIWXMAAA7EAAAOxAGVKw4bAAACDUlEQVQ4EbVTv2/TQBh9dpJz86N2J1QmSoIq\nQIi6LEioUliYGJqNjZilQxlaFYmBSsSCkSV/AbhDKQJV6gIIZYkYWIkYYE0nhMTgNkpNwch9d06w\nqyYSS086+73vffd9d+dnIBluAo+jyeA4T9jDscrWOKX07aSizfXOA8I4qcC6yW5nRyl1BzC8UcpK\nDZhstf5MJRsdoKp650f06Sller+D/L2tzdTC3ME/crepB4VlZ8hFf4iMXej2GXwZ8oxPtCpZ0Yfp\n3IKlCLnlAlpTKqaLF3dCmB1JOOoesxVqODhXOkR2SjFgpQ1sxEoHT0sBst5AWeSlxWfJeNoB15ix\nYlSiciWaUVnG5e2+EcKqDdakX5qP1N7Sim7jTXKetGI6EO/epyMDLKq3R0RPNxSNGf7ptv2/6utL\nl5i49LINvIqelC/QYoW53ge6dhfPm3hQK9rUpfGmOS1J9C70HTzDhAvkfjOwzimNLmXh5w6J+M1D\nPh5xSqOj4aIQTnSJABpPIE+gjG5yTd9aWH7LAI0Xe0gZXfjI7DXmlfvq166qsqofcANf/zYCFL8r\n48lNxf3o29efQv6OGm3P8jtoq34ygRsPsja0PeAXjW94sdGVkpnhobhGlRersdFZ7Ac+N+nlYjcu\n/5HJ0ui8o+vGPDCLxx1kfV7DT8ZodDnut3hg0bqCtep++eKijXwlmmX8CGL+omsFhVUkAAAAAElF\nTkSuQmCC\n", "text/latex": [ "$$\\frac{L^{4} q_{0}}{96 EI}$$" ], "text/plain": [ " 4 \n", "L ⋅q₀\n", "─────\n", "96⋅EI" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "sol1_L2 = (c1_sol1[0] * phi1).subs(x, L / 2).simplify(); sol1_L2" ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAKkAAAA/BAMAAABpxPVCAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAIqt2Zs0QmTK73URU\n74mR/c/RAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAFY0lEQVRYCbWYbWhbVRjH/2mSm5vXBnwp+2TW\nOVkV9K6T+cFh4zYcCjZxMCei9s6BxYptcNPMIS7ULxOpiZ0WCr60gvOFSctAUGTugmwiqC1OZkEd\nF7QiFWbTdktr3eJzzs1Nc9+Spl3Ph3vO+T/P88vJueflSQAqrg3sea1L4uFrTSSeu7AGUGRza0Ed\nia4BNTizBlD45lZGdXUUi/OOof7/mKnzj78cPUqG4AbDC/Ad3fqCbI4RXm7PcM27yKrH0FVrevue\n4Z48iB6yO5rS2+XaqwjDvONlX8MzjkiybLNvvGiS/S4r9UFgkrtxqj+G0GVTlLnbbBLWBWMmBegH\nxmSmcmpjDEE+vxa/JcFMPSiOLxlLLVpPpxXW5tSeFIKzJYtD1ftbM7q3/I3Pzv3wEHcZxnGzq0Db\n9LTEVE69SYV41exj6rdAHIZfws9nTWtxcg8re8ndRYysyuI0amo51AB99jweZVH2xUVjraAuYwbQ\ngqwCXMVueyJTTTNAbytU6221YGwUWAD7qlSKvGjt8pPe1pjCenwG/Dl4aq0sbawFnVomYeoTVt5n\nwiCQkFmDU2kX+JKsV6W0IEwrZbFMDXW8a/GmXfA8FzkVJ/BsxuJjFAYhTsCr4kNdTuiNpTqgCBO8\np1G7t2xastm2eovH0HTxNnxRfLtkb7H6Cb+3ZriqUa0OtZVzVVxWTBUMR5jpEwLOR6/J09QNTZuE\nym64mrHS0dyuenT2DJndl9n30i6nbWFfEqq9XlNtVNDg6PSpo6WGIRtFq5OLr9b2dwqMjKRvdVwF\nbYolbMCi1Cv4frFEeM5bpHqFsGqJ+NxycVhcagnCHll3yQFnFlQI8dVT0aUz2M3toR3htl5y+ucu\nv3YtlHy7idpA43394oXM8sOdPI+U9sdWooYl5nXeybUOvZQUi2eJ2sM+obsQhXDgbqUOhtW1cZpr\nEYGoY7q5K5ONdrfLerf+2pvnMfcyqpZ5UT+NPgwIufppeoR2agsSo9Kl/BrXf6UtPu0wwa8+oYdW\n1qLxZtCooXS6EAedCTHmKuRxzDWOjyrDtHaQjD9aZZhTH/2GcV2GmIc7xkMmxHl7aieZ++2oJk2n\n3j+jTBWeG2GrADh5z7T9DOwiYz1UTis/PEnbtxWi9Mye6jpajmUNfawGEQ0qOtujBo0lr9v/bVbQ\nf12rBDS1rw8+vm3y+pETD1waveOndTMQt+1IMtmZeuQRA5E6WvLKEpfBuCuP4HGcksMv9cR7JLrC\n/VFKFZ8GhrhMPvZjNTMBLXll1HeAOdARGpa8Q3R2JDGDCJuXUzIkLpOPRm3UEkfTc5rseskqLHll\n1LeAWWQP7twe98aAUD7Uhi8ZNbCwWZPJxzDW7wTjXJJdL1ryymeAURMKGejyBi5EworKqMJkMcVl\nZqnMXRKXiHpmJr0xj1DHm2RdKmysBRqrxJJUGqtEJk49vD8yoDLqjfAscplZKqn3yaS4x4EBoI2a\nFUVLXncjpVEDSQgSp2bXi7M0ln58DdzMZYoyUHfuymhntwr6HWAoWvJ6WKeKc/DJgRS5EOgDqgbx\nvYw3uEw9AxXiEJ3dqijHAeNxQQuUklc03B7tLW769koMTx74M/jegkJfLY47gd7iLV/taB9lMkHL\n1Bvo/WeStE/p7PaTXjWRZIFViz7WQ690PoWcmKSzO/0xRVRNJKsSmbFBS4kEbPbRVXgXzfsU9pMe\nGWLWlZbyX0S5phJiL/YJcW3VrBQK4YoW6sqdLDFo8BEZjcqKkSzwH4mH++J9GobObnxDv0dpLlZR\nSq8rEA3InDJVSG+cR2Tk0CqYtIQ6RlcV7xDsNv3ed3CrV95Xb0AN//8BkB+cGWZgemoAAAAASUVO\nRK5CYII=\n", "text/latex": [ "$$\\begin{cases} 0 & \\text{for}\\: \\frac{\\pi}{L} = 0 \\\\\\frac{4 L^{4} q_{0}}{\\pi^{5} EI} & \\text{otherwise} \\end{cases}$$" ], "text/plain": [ "⎧ π \n", "⎪ 0 for ─ = 0\n", "⎪ L \n", "⎪ \n", "⎨ 4 \n", "⎪4⋅L ⋅q₀ \n", "⎪─────── otherwise\n", "⎪ 5 \n", "⎩ π ⋅EI " ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "sol2_L2 = (c1_sol2[0] * phi2).subs(x, L / 2).simplify(); sol2_L2" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Noticing that all solutions have the same $\\frac{L^4 q_0}{EI}$ term, we will compare only the scale factor out front. For the exact solution we have" ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAMoAAAAPBAMAAABXbk2cAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAEJmJZjLNVN0i77ur\nRHZ72Yd1AAAACXBIWXMAAA7EAAAOxAGVKw4bAAACiUlEQVQ4Eb2TQWjTUBjH/zGpaVrTRncZImvs\nwIsOip0iG0hBT14MngW7kweFBY8y6NCDiIJBPehJdxCkY9ibyBSrTPAwbPEkXix4FcvcsHPT1e99\n39tSgmcD+X//fL/H++e95AH7xo9BXdsVuFd+DhijpwI4xcVAHEa+jmnC2ClPhmCR4dxMEMHAOZr9\nIoYjlSL1qgfHR62C3aFxEwdg/hEHH3tCJoLvwl4DS9xMEsY08y3AegyzSpar8bblIddDfhZngON4\nA9TF2RU4VSaCuyHWwRI3k4QxkJsHsk3YK5Sia9eDu4BCFS+BmvcUaAXszAaMWSaCP7eNn2CJm0nC\nGLh/Ccg34fYoRVdKoYt27DfwJZz2KIWdu+mZFSGMoXZMixqum+xiwq5EKQUf7hYN0pVTjOswflBK\nh/p39oprrY3TExEtGC6RZYmb7AaIwlZEKY9KSK3SeF1VinPiIBzqXaAx7qp2Zn9KCGP63a7RSBY1\nXJrsBgjjEagUX6dIlR2baDu0ApWSbmo382mjTS8zsSPmAj2yDDQZx0S5kkr5144hM7e9Y0W9d1YT\n3TmaNrMjqNNiROImuwFS9+xApdBXt/XXV5XW4kTI9tTXr4XI+RCXiZDqMRH8EHgdssTNJGE8BJWS\nbcDiP1kqpeRXVMoicNbDFcpkV6D3XmIiuO9RCkvcTBLGM8vLGx/UadxVpSl0pZRMA+lf6iweQspH\nLmJHa8ELJoLn6e+LWOJmkjCmqaniBvYHdMK48o5ZFdQ6SIfGEzwoH30vzp2C1WEi+DTsLbDEzSRh\nTAmbdA99fwfc1vXIs8sVfCy+okNx8luA6X5/XRzOj45BCGO3eLgNloFmggjGZH+JYv7D9ReRBVoi\nuN8ZfgAAAABJRU5ErkJggg==\n", "text/latex": [ "$$0.013020833333333334$$" ], "text/plain": [ "0.013020833333333334" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "5 / 384." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For the first approximate solutions we have" ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAMoAAAAPBAMAAABXbk2cAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAEJmJZjLNVN0i77ur\nRHZ72Yd1AAAACXBIWXMAAA7EAAAOxAGVKw4bAAAB1ElEQVQ4EbWUO0sDQRSFz2YTE40xq4KICAkW\nFhIh+GhDKsHKgD/AFIqNgo2NlZU2Fv4DkyqNYLDUwuCjESH5BwbtLJSgRAVfd+7MZMbGqdzinMv5\nljmZnd0AAzOzEJf2zhRJcw74OXHH9H3BFg4dxMKrGN4Xy0nfCvSEWBWIlgFvHePAUCNUtIRDB7Fw\nuAS/SCXs3nk9kBMlD1WMZlpAPI1T4BjxqiUcOoiF4zVEaSkofwr01D1He+kjdCu2mnj/JTJ0EIOR\nrCHxSisopxY1+V2q5YQo4qVfAg4dxGCk0kh80QrKqUVN87qlvXjXQHIlN2kJOHQQg3GQReSZWpRT\ni5y8pmrx2nlUkFpGz4YRDh3EwjhIqxbpooUnH7rlO8BEkGohVDbiibD/b2Jh/Xy0d57Ypm7BG7DQ\nSG4j/GlEhg5isDjrqDp94dTCiZfttOxQSz5WRPjFCDh0EIPFVxDmN1k6tXASvbmpHzb5Tb4We+kt\n0V6MgEMHMVh8g/RR81cpnFp00qPe5CU6AkToXEpGwKGDGAzsYqTgtaVzCyfUm1QtvVmvAlxgZMMS\nDh3Ewhh8vAL2ID1ztJaXE/yzj3zoun0JTGUK9Jc5RrcZkaGDGEw/+v+vHx/NTkw3kMKLAAAAAElF\nTkSuQmCC\n", "text/latex": [ "$$0.010416666666666666$$" ], "text/plain": [ "0.010416666666666666" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "1 / 96." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For the second approximate solution we have" ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAALQAAAAPBAMAAAC/7vi3AAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAEJmJZjLNVN0i77ur\nRHZ72Yd1AAAACXBIWXMAAA7EAAAOxAGVKw4bAAACuklEQVQ4EbVUTWsTURQ9k48mk492sBuRYtIU\nFaFgMK2IFg3oqgsNXbnrgAsXig0uRUhBUKQFQ3UhbmoXfhZpXClUIZYqIqUN/gCNLgRB1Co0Vq3j\nue8j/QU+yLnv3nPPmTfvvQmwZWAQMmwErhUe6zQ62VcoKBJO35ESOt5Hd0k6Apz1Rso4sVAo5N3c\nfAkKdKdSKxqnsLUqRR3Pe3CzqBRVGgmC4A8QmwE6ys4VJIIgL60TwFzwAzhOvt6D8AYUAG21oiO3\nEPbZrqKzsOKhs4WucZWG+Mw6evpXgWFgHzqGPrETnbPAyWVyu4EUngMPNGyqNZ2sI0YpTPzqIf0Q\nGV+lYSBcohX5p0DFS4kxcP20PJGDrzCIe8BKSQFg1ZruqiPdYpuJJDkqRVt+y0ysuS/vysY6b62B\ndBNjHq0VKGul1taZLNJ/mZuorJ2LNkVNWzvfad1MHf3Al4hUab14bAcZJAVwVdYjYNWans4jyhOB\niUK6+3ttGi6S4qpdtozmk567AWwHrR9hukFqgD+kRa/AqjU9nTXWOqrn4kDDlDMiFWuuejTP+RPZ\nYVpzweOA801m8boFq2ZO2myEjZpMzJjyRzaJtd4QzvdUYyVtHVon4wufa4NVs0Ca5xUzxyiRpFtF\nsmXKE6Iyx1gpvwQON7pB65QvtxnxmtBZC221ppM1RNTl05HWXau01mWHF0Nbz/MD8e5w1d6FpaVf\nywkfIaoyddLnZDEKxFqpNc1vI+SzwUSSiRri6zqNcou1NT+ZncgCl6QwixBX7HP7WYlm0VlVoN5Z\nqQ19GdtKzhogUZGRIipNncba1vGycxtTiNwX69/yZzBc5f2n9Y3C3lcaNtWG7v7yApgEVOyfO1PE\nm9wzk7p3aRR6vbYI59Bn/gn1LXgsHAwWMTXUy9nNMjAWBD81YFOtaXb8n/EPnQ8W/mr1dEMAAAAA\nSUVORK5CYII=\n", "text/latex": [ "$$0.0130710545722135$$" ], "text/plain": [ "0.0130710545722135" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "pi = 4 * atan(1.0);\n", "4 / pi ** 5" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The second approximation is better than the first. This is because $\\phi=\\sin\\left(\\frac{\\pi x}{L}\\right)$ satisfys both the essential and natural boundary conditions whereas $\\phi=x(x-L)$ only satisfies the essential boundary conditions. Satisfying the natural bondary conditions is not a requirement of $\\phi$, but will yield better approximations with fewer terms." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Problem 2\n", "\n", "We have seen that the following one-dimensional boundary value problem describes the physics of many interesting problems in engineering\n", "\n", "\\begin{equation}\n", "\t-\\frac{d}{dx} \\left[ a \\frac{du}{dx} \\right] + c u - f = 0 \\quad \\text{for} \\quad 0 < x < L, \\notag\n", "\\end{equation}\n", "\n", "subject to the boundary conditions\n", "\n", "\\begin{equation}\n", "\tu(0) = u_0 \\quad \\text{and} \\quad \\left[ a \\frac{du}{dx} \\right]_{x=L} = Q_0, \\notag\n", "\\end{equation}\n", "\n", "where the $x$ is the independent variable, $u=u(x)$ is the dependent variable, and $a = a(x), c = c(x), f=f(x), u_0$, and $Q_0$ are the *data* of the problem.\n", "\t\n", "\n", "### (a) 30 points\n", "\n", "Write a general one-dimensional finite element code to solve this problem. General means that the user will input the *node* locations and specify the functions $a(x)$, $c(x)$, and $f(x)$ as well as the boundary conditions $u_0$ and $Q_0$. The output of the code should be the value of $u$ at the user specified nodes. You can restrict the code to using only linear interpolated elements.\n", "\n", "You can verify your code with following data and results:\n", "\n", "\n", "1. For $a(x) = 1 - x/2$, $c(x) = 0$, $f(x) = 0$, $u_0 = 0$, and $Q_0 = 1$,\t\n", "\n", "|$x$|$u(x)$|\n", "|-|-|\n", "|0.00|0.00000|\n", "|0.25|0.26666|\n", "|0.50|0.57435|\n", "|0.75|0.93799|\n", "|1.00|1.38244|\n", "\n", "2. For $a(x) = 1 - x/2$, $c(x) = x$, $f(x) = x$, $u_0 = 0$, and $Q_0 = 1$,\n", "\t\n", "|$x$|$u(x)$|\n", "|-|-|\n", "|0.00|0.00000|\n", "|0.40|0.46695|\n", "|0.60|0.73242|\n", "|0.65|0.80319|\n", "|0.70|0.87615|\n", "|1.00|1.38355|\n", "\n", "\n", "3. For $a(x) = 1$, $c(x) = 0$, $f(x) = x^2$, $u_0 = 0$, and $Q_0 = 2$,\n", "\t\n", "|$x$|$u(x)$|\n", "|-|-|\n", "|1|0.000|\n", "|2|22.08| \n", "|3|40.00|\n", "|4|48.75|\n", "\n", "\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Solution**\n", "\n", "First we will write a general class that solves an arbitrary one-dimensional problem." ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import numpy as np\n", "import scipy.linalg\n", "\n", "class OneDimFEA():\n", " \n", " def __init__(self, nodes, a=lambda x: 1, c=lambda x: 0, f=lambda x: 0, u0=0, Q0=1):\n", " \"\"\"\n", " Initiates an object from the argument list to perform FEA solution on the following\n", " \n", " -\\frac{d}{dx} \\left[ a \\frac{du}{dx} \\right] + c u - f = 0 \\\\\n", " \n", " \\text{for} \\quad 0 < x < L, \\\\\n", " \n", " \\text{with} u(0) = u_0 \\quad \\\\\n", " \n", " \\text{and} \\quad \\left[ a \\frac{du}{dx} \\right]_{x=L} = Q_0\n", " \n", " \n", " assigns the following default properties if none are given in the argument list\n", " \n", " input: nodes - list of node locations, ex. [0, 1.2, 4.7, 8]\n", " input (optional): a - function of x (default = 1)\n", " input (optional): c - function of x (default = 0)\n", " input (optional): f - function of x (default = 0)\n", " input (optional): u0 - essential boundary conditional at x = 0 (default = 0)\n", " input (optional): Q0 - natural boundary condition at x = L (default = 1)\n", " \n", " output: u at the node locations\n", " \"\"\"\n", " \n", " self.nodes = np.array(nodes, dtype=np.double)\n", " self.a = a\n", " self.c = c\n", " self.f = f\n", " self.u0 = u0\n", " self.Q0 = Q0\n", " \n", " self.K = np.zeros((len(self.nodes), len(self.nodes)))\n", " self.F = np.zeros(len(self.nodes))\n", " \n", " def N1(self, x):\n", " \"\"\"Computes the first linear shape function\"\"\"\n", " nodes = self.nodes\n", " return (x - nodes[1:]) / (nodes[:-1] - nodes[1:])\n", " \n", " def N2(self, x):\n", " \"\"\"Computes the second linear shape function\"\"\"\n", " nodes = self.nodes\n", " return (nodes[:-1] - x) / (nodes[:-1] - nodes[1:])\n", " \n", " def dN1dx(self, x):\n", " \"\"\"Computes the derivative of the first linear shape function\"\"\"\n", " nodes = self.nodes\n", " return 1 / (nodes[:-1] - nodes[1:])\n", " \n", " def dN2dx(self, x):\n", " \"\"\"Computes the derivative of the second linear shape function\"\"\"\n", " nodes = self.nodes\n", " return 1 / (nodes[1:] - nodes[:-1])\n", " \n", " def compute_element_stiffness_and_force(self):\n", " \"\"\"\n", " Computes the components of the element stiffness matrix and\n", " element nodal forces for all elements\n", " \"\"\"\n", " \n", " nodes = self.nodes\n", " \n", " #We will use Gauss integration at the following points\n", " t1 = np.sqrt(1.0 / 3.0)\n", " t2 = -np.sqrt(1.0 / 3.0)\n", " \n", " #Because the points above are defined on the domain -1 < x < 1 and our\n", " #elements are defined on arbitrary domains, we can use a change of variables\n", " #to rescale the integration bounds\n", " x1 = ((nodes[1:] - nodes[:-1]) * t1 + nodes[1:] + nodes[:-1]) / 2.0\n", " x2 = ((nodes[1:] - nodes[:-1]) * t2 + nodes[1:] + nodes[:-1]) / 2.0\n", " \n", " #Evaluating the shape functions and thier derivatives at each Gauss \n", " #integration point.\n", " N1x1 = self.N1(x1)\n", " N2x1 = self.N2(x1)\n", " dN1dx1 = self.dN1dx(x1)\n", " dN2dx1 = self.dN2dx(x1)\n", " \n", " N1x2 = self.N1(x2)\n", " N2x2 = self.N2(x2)\n", " dN1dx2 = self.dN1dx(x2)\n", " dN2dx2 = self.dN2dx(x2)\n", " \n", " #Same for the functions\n", " a1 = self.a(x1)\n", " a2 = self.a(x2)\n", " c1 = self.c(x1)\n", " c2 = self.c(x2)\n", " f1 = self.f(x1)\n", " f2 = self.f(x2)\n", " \n", " #Now we evaluate the element stiffness matrix components, the 0.5*(b-a) term out front\n", " #is due to the change of variables\n", " self.k11 = (nodes[1:] - nodes[:-1]) / 2 * (a1 * dN1dx1 * dN1dx1 + a2 * dN1dx2 * dN1dx2 + \n", " c1 * N1x1 * N1x1 + c2 * N1x2 * N1x2)\n", " self.k12 = (nodes[1:] - nodes[:-1]) / 2 * (a1 * dN1dx1 * dN2dx1 + a2 * dN1dx2 * dN2dx2 + \n", " c1 * N1x1 * N2x1 + c2 * N1x2 * N2x2)\n", " self.k21 = (nodes[1:] - nodes[:-1]) / 2 * (a1 * dN2dx1 * dN1dx1 + a2 * dN2dx2 * dN1dx2 + \n", " c1 * N2x1 * N1x1 + c2 * N2x2 * N1x2)\n", " self.k22 = (nodes[1:] - nodes[:-1]) / 2 * (a1 * dN2dx1 * dN2dx1 + a2 * dN2dx2 * dN2dx2 + \n", " c1 * N2x1 * N2x1 + c2 * N2x2 * N2x2)\n", " #Same for the element force vector components\n", " self.fe1 = (nodes[1:] - nodes[:-1]) / 2 * (f1 * N1x1 + f2 * N1x2)\n", " self.fe2 = (nodes[1:] - nodes[:-1]) / 2 * (f1 * N2x1 + f2 * N2x2)\n", " \n", " \n", " def assemble(self):\n", " \"\"\"\n", " Assembles the global stiffness matrix and force vector\n", " \"\"\"\n", " self.K += np.diag(np.r_[self.k11, 0.0])\n", " self.K += np.diag(np.r_[0.0, self.k22])\n", " self.K += np.diag(self.k12, k=1)\n", " self.K += np.diag(self.k21, k=-1)\n", " \n", " self.F[:-1] += self.fe1\n", " self.F[1:] += self.fe2\n", " \n", " \n", " def apply_bcs(self):\n", " \"\"\"\n", " Applies essential and natural boundary conditions\n", " \"\"\"\n", " \n", " row = np.zeros(len(self.nodes))\n", " row[0] = 1.0\n", " \n", " self.K[0] = row\n", " self.F[0] = self.u0\n", " \n", " self.F[-1] += self.Q0\n", " \n", " \n", " def solve(self):\n", " \"\"\"\n", " Solves the one-dimensional finite element problem.\n", " \"\"\"\n", " self.compute_element_stiffness_and_force()\n", " self.assemble()\n", " self.apply_bcs()\n", " \n", " return scipy.linalg.solve(self.K, self.F)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we can solve the three verification problems above to test the code." ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "array([ -3.55271368e-16, 2.66666667e-01, 5.74358974e-01,\n", " 9.37995338e-01, 1.38243978e+00])" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "nodes1 = [0.0, 0.25, 0.50, 0.75, 1.00]\n", "problem1 = OneDimFEA(nodes1, a=lambda x: 1 - x / 2)\n", "problem1.solve()" ] }, { "cell_type": "code", "execution_count": 15, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "array([ -5.95648717e-16, 4.66641976e-01, 7.31924142e-01,\n", " 8.02638676e-01, 8.75540741e-01, 1.38162220e+00])" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "nodes2 = [0.0, 0.40, 0.60, 0.65, 0.70, 1.00]\n", "problem2 = OneDimFEA(nodes2, a=lambda x: 1 - x / 2, c=lambda x: x, f=lambda x: x)\n", "problem2.solve()" ] }, { "cell_type": "code", "execution_count": 16, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "array([ 0. , 22.08333333, 40. , 48.75 ])" ] }, "execution_count": 16, "metadata": {}, "output_type": "execute_result" } ], "source": [ "nodes3 = [1, 2, 3, 4]\n", "problem3 = OneDimFEA(nodes3, a=lambda x: 1 , f=lambda x: x ** 2, Q0=2)\n", "problem3.solve()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## (b) 20 points\n", "\n", "A one-dimensional heterogenuous porous medium of length $L=1$ m, has a steady state pressure distribution as shown in the following table. \n", "\n", "|$x$ (m)|$p(x)$ ($\\times$ kPa)|\n", "|---|------|\n", "|0.00|100.0|\n", "|0.05|100.005|\n", "|0.10|100.010|\n", "|0.15|100.014|\n", "|0.20|100.018|\n", "|0.25|100.021|\n", "|0.30|100.024|\n", "|0.35|100.026|\n", "|0.40|100.029|\n", "|0.45|100.030|\n", "|0.50|100.032|\n", "|0.55|100.033|\n", "|0.60|100.034|\n", "|0.65|100.035|\n", "|0.70|100.036|\n", "|0.75|100.037|\n", "|0.80|100.038|\n", "|0.85|100.038|\n", "|0.90|100.039|\n", "|0.95|100.040|\n", "|1.00|100.040|\n", "\n", "\n", "The porous media is subject to the a constant pressure a $x=0$ of $p = 100$ kPa and an exit flow rate of $Q = 100$ m/s is measured. \n", "\n", "Use the code developed in part (a) to help you estimate what the permiablity $\\kappa(x)$ this material. You can use a consant viscosity of $\\mu = 1$ Pa$\\cdot$s. Clearly indicate your answer and explain your approach (preferably with plots). The total diffusivity coefficient is $K = \\frac{\\kappa}{\\mu}$.\n", "\n", "**Note:** Submit a working version of your code to [Canvas](https://utexas.instructure.com/courses/1119539). Any supplemental material explaining your answer in part (b) can be turned in to me via hard copy or scanned and submitted to Canvas with your code." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Solution**\n", "\n", "Below we define our objective function which is to minimize the $L_2$-norm of the difference between the pressure output by our FE code and the known pressure published above at the node ($x$) location." ] }, { "cell_type": "code", "execution_count": 17, "metadata": { "collapsed": false }, "outputs": [], "source": [ "def objective_fun(c, nodes, p):\n", " \"\"\"\n", " Computes the objective function for minimization.\n", " \n", " input: c - the initial guess for the optimization routine, i.e. the coefficients\n", " of the polynomial which is of quadratic form\n", " input: nodes - the nodal locations, needed for the FE code\n", " input: p - the known pressures\n", " \n", " output: the scalar output to minimize\n", " \"\"\"\n", " \n", " \n", " problem = OneDimFEA(nodes, \n", " a=lambda x: c[0] + c[1] * x + c[2] * x ** 2, \n", " u0=100000, \n", " Q0=100)\n", " \n", " sol = problem.solve().flatten()\n", " \n", " return np.linalg.norm(sol - p)" ] }, { "cell_type": "code", "execution_count": 18, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " fun: 1.1269859125646242\n", " success: True\n", " nit: 145\n", " nfev: 257\n", " message: 'Optimization terminated successfully.'\n", " status: 0\n", " x: array([ 0.95633747, 0.57572211, 8.93686525])\n" ] } ], "source": [ "import scipy.optimize\n", "\n", "#Create the nodal array and define the known pressures from the problem statement\n", "nodes = np.linspace(0,1,num=21)\n", "p = np.array([100.000, 100.005, 100.010, 100.014, 100.018, 100.021, 100.024, 100.026, \n", " 100.029, 100.030, 100.032, 100.033, 100.034, 100.035, 100.036, \n", " 100.037, 100.038 ,100.038, 100.039, 100.040, 100.040], dtype=np.double) * 1000\n", "\n", "#Run the optimization routine\n", "res = scipy.optimize.minimize(objective_fun, [1,1,1],args=(nodes,p), method='Nelder-Mead');\n", "#Extract the solution\n", "c = res.x\n", "print(res)" ] }, { "cell_type": "code", "execution_count": 19, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZIAAAEZCAYAAAC99aPhAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XmcVNW19//P6mYGAQFFpZFJkElsRrtxCFGjDI5EI5BE\nMVHJc/VGk1zFXDXGPD/jNdcBfeQJITEa/WlQYxKJFxVlEtpmEhoBm7GZRWUebIYe1vNHFVRbMnR3\nVdepqv6+X696UfvUPn1WL6trefY++5S5OyIiItWVEXQAIiKS2lRIREQkJiokIiISExUSERGJiQqJ\niIjERIVERERiUisLiZn1MrOPzGyJmb1lZk2O0+9uM1safvz0ZPubWX0ze9XMPjGz5WZ2fyViudPM\nVptZmZm1iN9vKSKSGGlfSMzsW2b2QtTmPwH3ufv5wD+A+46xXw/gx0A/IBu42sw6nmT/EQDu3iu8\n3xgzO/skIc4BLgM2VPV3ExFJBmlfSMKiV112dvc54ecfAN89xj7dgHnufsjdy4BZwPDwa12Os//n\nQGMzywQaAYeAvQBm9p3wWcxCM3vNzBoBuPsSd98IWOy/pohI4tWWQhL9Ib3czK4JP/8ekHWMfZYB\nF5vZqeEP/aFA2yOvHWt/d3+PUOHYCqwHnnD33WbWEngQuMzd+wEfA7+Iy28mIhKwOkEHUFPMbC5Q\nDzgFONXMFoVfGgv8CPg/ZvYQMBk4HL2/u68ws8eB94H9wGKgLPzyj4Fno/c3sx8ADYEzgJbAbDP7\nAOgBdAfyzMyAukB+3H9pEZEApG0hcfccCM2RALe4+4+iulwZfr0zMOw4P+MF4IVwv0eBTeHtK4+z\n/0DgH+5eDmwzszxCcyUHganu/v0ThVzV31FEJBkkfGjLzAab2QozW2VmY4/T59nwlUwFZpYd3tbF\nzBab2aLwv3sqXklVxRhOC/+bQWjIacJJ+p0NXA+8epz9fx/eZQWhiXPMrDGQE942F7jQzDqFX2sU\nLkBfOxyaJxGRFJTQQhL+4H2O0P/N9wBGmlnXqD5DgE7u3hkYQ/hD3t1XuXtvd+8D9AW+InTFVHWM\nNLOVwKfAFnd/MXzsM83s7Qr93jSzZcBbwL+5+97j7P+X8PY/APXMbCkwD3je3Ze5+3ZgNPBXM1sC\nfAScGz7mv5vZJqANsMTMJlbzdxIRCYQl8jbyZpYDPOzuQ8Lt+wF398cr9JkAzHD318LtQmCQu39R\noc8VwEPufnHCghcRkWNK9NBWG8LzDGGbw9tO1GfLMfrcBPw17tGJiEiVpdzlv2ZWF7gGeCPoWERE\nJPFXbW0BKq70zgpvi+7T9gR9hgAfu/u24x3EzHQFlIhIFbl7tS74SfQZyQLgHDNrZ2b1CN1SZHJU\nn8nAzXB0TmV3xfkRYCSVGNZydz3cefjhhwOPIRkeyoNyoVyc+BGLhJ6RuHuZmd0FTCVUxJ5390Iz\nGxN62Se6+xQzG2pmawhdmXXrkf3DK8wvB+5IZNypbP369UGHkBSUhwjlIkK5iI+EL0h093cJX/pa\nYdsfotp3HWffYuC0motORESqKuUm26VqRo8eHXQISUF5iFAuIpSL+EjoOpJEMTNPx99LRKSmmBme\nIpPtkmAzZ84MOoSkoDxEKBcRykV8qJCIiEhMNLQlIiIa2hIRkeCokKQ5jQGHKA8RykWEchEfKiQi\nIhITzZGIiIjmSEREJDgqJGlOY8AhykOEchGhXMSHComIiMREcyQiIqI5EhERCY4KSZrTGHCI8hCh\nXEQoF/GhQiIiIjHRHImIiGiOREREgqNCkuY0BhyiPEQoFxHKRXyokIiISEw0RyIiIpojERGR4KiQ\npDmNAYcoDxHKRYRyER8qJCIiEhPNkYiIiOZIREQkOCokaU5jwCHKQ4RyEaFcxIcKiYiIxCThcyRm\nNhgYR6iIPe/ujx+jz7PAEOArYLS7F4S3NwP+BPQEyoEfufu8Y+yvORIRkSpImTkSM8sAngOuBHoA\nI82sa1SfIUAnd+8MjAEmVHj5GWCKu3cDzgcKExK4iIgcV6KHtgYAq919g7uXAJOAa6P6XAu8BBA+\n22hmZq3NrClwsbu/EH6t1N33JjD2lKQx4BDlIUK5iFAu4iPRhaQNsKlCe3N424n6bAlv6wBsN7MX\nzGyRmU00s4Y1Gq2IiJxUnaADqII6QB/gTndfaGbjgPuBh4MNK7kNGjQo6BCSgvIQoVxEpHIuiouL\nWbhwIfn5+Xz88cccOHDghP3rlJXR8uBBTjt4kJYHDtDq4EFahf9tefBgTLEkupBsAc6u0M4Kb4vu\n0/Y4fTa5+8Lw878BY493oNGjR9O+fXsAmjdvTnZ29tE3zZHTWbXVVlvtVGi7O+3btyc/P5833niD\nTz/9lM2bN9OjRw/atm1L986dGdSxIw137GDZokXU37uXS+rUoeH27SzcvJl6e/Zw+aFDHGzRgvcb\nNuRw06YcbtWKKcXFbDx8mNIGDWDhkY/WqkvoVVtmlgmsBC4DtgLzgZHuXlihz1BCZx3DzCwHGOfu\nOeHXZgG3u/sqM3sYaOTu3ygmumorYubMmUffkLWZ8hChXEQkay4qnm3k5+ez9KOPONud73TuzIDW\nrenaqBFnHDpEnc2bYcMG2LEDzjoL2raFrKxv/puVBa1bQ8bxZzNiuWoroWck7l5mZncBU4lc/lto\nZmNCL/tEd59iZkPNbA2hy39vrfAjfgq8YmZ1gaKo10REUo6Xl7Nh0SIK33uPLXl57F26lPqff07P\nJk34ft263FNcTN3ycqx9e6xZs1BBaNfu648zz4TMzMB+B91rS0Skph06FDpzKCri8IoVfJmfT/Gy\nZdTdvJnT9u2jBNjVtCllWVk07taN0/r1o+4554SKRPv20LIlWLVOFiotljMSFRIRkVi5h4aXiopg\n7drQv0VFeFERZStXYtu2sbNRI9a6s+zAAQ6ceSaNzzuPNhdfTLdhw2jbsydWw4XiZFRIoqiQRCTr\nGHCiKQ8RykVElXJxpFisWhV5rF4dehQVQWYmZR06sLNZM9a68/Hu3czYsIHN9epx9oUXMmDgQHJz\nc+nTpw8NGybfyoWUmSMREUkW27dvp6CggLKysq9tzzxwgIabN9No82YabdkS+jf8nPJyirOyKM7K\n4kBWFsXnnsumXr2YvmEDMwoKKCwspGfPnuTm5pKbm8tTubm0bds28LONmqYzEhFJe6WlpSxbtoy5\nc+eSn5/PvI8+ouHnn3NFVhYdSkvJKi4++mhSWspnDRuyuVGj0KPC8z11635jruLUU08lJyeH3Nxc\n+vbtS4MGDQL6LWOjoa0oKiQitdv27duZO3cu8+bMYdOMGZR+8gn9mjQh55RTOKekhBbbt2NnnIF1\n6wbnngudO0OXLqFHVtYJL5NNVyokUVRIIjQeHqI8RKRbLsrKyli+aBErJk9m26xZ+PLlZO3dS3b9\n+rQ5eJDDrVtTt1cv6mVnQ/fuoUfXrtC4cdrlIhaaIxGR2qGsjF0LFlD0j3+wLy+POitX0nrHDs4F\nWjVtyqGOHWl8ww20uuQSMnr2hC5dqJuEE9vpRmckIpKc9uyhbPFitr73HnvnzKH+ihWcsWMH24At\nLVtyuGtXml10ER2vuorm/ftDvXpBR5zSNLQVRYVEJIWUl4cun/3kE4rz89k7ezb1V6ygwb59LAeK\nTjmFQ1270vSSS+h07bV0y8khM8BV3OlKhSSKCkmExoBDlIeIQHNRUgLLlsGCBZQvWsSB/HzqrlzJ\nvsxMlpjxcUkJB7p0oenFF9NlyBAuGDiQli1b1lg4el9EaI5ERJJPWRmsWAELFsDChZTk52PLl7Oj\naVMKMjOZtmsXW08/nabDh3P+oEHk5uby8+7ddbaRgnRGIiKxc4c1a0K3Il+wAA+fcexr0oTCRo2Y\nsX8/HxYXU6d/f3pffDG5ublccMEFNXq2IVWjoa0oKiQiNcgdNm06eqbBggWUL1zIwXr1KGrZkvzD\nh5m8dStftGlDj4suOrrKu7vONpKahrbkuDQGHKI8RFQ5FyUlUFAAeXmQl4fn5VF6+DBbzjqLRRkZ\nvLNtG9PLyujQq9fRovGXnBxatGhRY79DvOh9ER8qJCLydXv2QH5+pHDMn8++Vq1Y0aIF0w4d4tW9\neylp04bcvn3Jzc3lp7m5TNDZRq2moS2R2sw99D0ZeXkwZw6el0f5mjV83rYtH9evz7927uS9vXvp\nMmDA0bONnBQ525Cq0RxJFBUSkeMoKwsNU82ZA3l5lM+ezeFDh1jTujUflpbyxmefsS0ri/7hW55r\nbqP2UCGJokISoTHgkFqbB/fQJbjTpuEffEDZjBlMyczklNNO4529e5m6fz+tBgwgN1w4atuVVLX2\nfXEMmmwXkYhNm2DaNA5NmYJPm8bB0lLyGzfmzV27WNmmDY3POYfrr7+eH+Tm8liPHjrbkJjpjEQk\n1e3YQdm0aez629+oM2sWmXv2MLtuXd4rKWFXnz60u/RScgcOrHVnG1I1GtqKokIiae2rr9jz9tts\nf+01Gn70Ec22bycPWNyiBV/l5JA1dCg5AwfSQ2cbUgUqJFFUSCI0BhySynkoKy1l9dtvs+vll2me\nl8fZX35JQUYGa84+m8MXX0zW8OEMuOiiSp9tpHIu4k25iNAciUga2b59O/Nmz+bzN96g2ezZ9N6y\nheaZmWzt2JH13/0uZT/8ITn9+3OhzjYkSeiMRCRAZWVlLFu2jPz8fJbOmEGjWbPI3bGD7wC7Tz+d\nA5deyuk/+hHNBw36xneFi8SThraiqJBIsjryXeL5+fnkf/QRX82fz/caNeKajAzO3ruXgwMHcsrI\nkWRcdRWccUbQ4UotEkshqX3fcF/LzJw5M+gQkkKQeSgtLeX1118nNzeXrh06MOehh7h66lTeLizk\no5Yt+cWIEXR+6SXq79pFs2nTyLjtthotInpPRCgX8aE5EpEasmfPHp5//nkmPP0032vShDdOPZU2\nGRlYo0Zw1VVw9dXQo4eGrCTlaWhLJM7WrVvH+KefZvOf/8xdrVqRs3MndQYMgJtuguuug9NOCzpE\nkW/QHEkUFRIJQv7s2XzwwAN0mDeP6zMyyDzvPBrccgvccAO0bh10eCInlFJzJGY22MxWmNkqMxt7\nnD7PmtlqMysws94Vtq83syVmttjM5icu6tSlMeCQmspD6eHDzHjkEd484wy6DBrE7Rs3csOvf03j\nVatoMH8+3Hln0hURvScilIv4SOgciZllAM8BlwGfAQvM7C13X1GhzxCgk7t3NrMLgN8DOeGXy4FB\n7r4rkXGLfI07+6dNY+X//t+0ycsjq0EDTr/uOpo//DCZnTsHHZ1IwiV0aMvMcoCH3X1IuH0/4O7+\neIU+E4AZ7v5auF1IqHh8YWbrgH7uvuMkx9HQlsSXOxQUsOcPf6D0r39lx/79fNK9O+f+6lecd+ON\nQUcnErNUWtneBthUob0ZGHCSPlvC274AHHjfzMqAie7+xxqMVQTfuZOdzz5Lxp//zOGdO3m1rAy/\n6SZu/M1vuOHss4MOTyQppNrlvxe6+1YzO41QQSl09znH6jh69Gjat28PQPPmzcnOzj56T50j46K1\noV1xDDgZ4gmqXVBQwD333HPS/sXFxfxx4kS2TZvG5QUF9N6yhf9bpw5runblgvvu47ZbbuHjjz9m\nbVERbcOFJBl+v6q0x40bV2v/HqLbtfnv48jz9evXE6sghrZ+7e6Dw+3KDG2tAL7l7l9E/ayHgX3u\n/tQxjqOhrbCZuikdcOw8uDvr1q07utK88MMPGVBYyB2ZmdRv1Iitw4Zx2s9/Ttb552NptNZD74kI\n5SIiZS7/NbNMYCWhyfatwHxgpLsXVugzFLjT3YeFC884d88xs0ZAhrvvN7PGwFTgEXefeozjqJDI\nNxQXF7Nw4cLQ7Uny85k7dy6ZZtzRqRMj9u2jU1ERXHcddf7X/4LcXC0UlFolZeZI3L3MzO4iVAQy\ngOfdvdDMxoRe9onuPsXMhprZGuAr4Nbw7q2Bf5iZh+N+5VhFRCRaXl4e9957L0uWLKFnz57k5uYy\n+ooreKFzZ5q/+SZWXAw/+QmMGgXNmgUdrkjK0YLENFebT92Li4t58MEHmTRpErfddhu/vPdeGs6a\nBX/8I8yeDd/7Htx+O/TtG3SoCVWb3xPRlIuIlFqQKJIIeXl5ZGdns3XrVpa+/z6XbthAw27d4NFH\n4dprYeNGmDCh1hURkZqgMxJJKxXPQl68/36uWLoU3nwTRo4MDV+dd17QIYokJZ2RiBA+Czn/fJou\nWcK6Xr244tFHoW1bWLUKxo9XERGpISokaa7iNePpqri4mF/ccw8vXHUVc8349caN1L/uOli/Hn71\nK2jVqlbkobKUiwjlIj5SbUGiyNfkT5/OuyNG8B8HDtCqa1fq/vKXoTkQfZ+5SMJojkRSUvHGjXx4\n0030nT+fkr59Oeupp+DCC7X2Q6SaNEcitcfatWwdPpzDHTpQ78svqTNnDmfNnw8XXaQiIhIQFZI0\nlzZjwPPnUzp8OPt79uT1qVPJmziRS9eu5dTc3ErtnjZ5iAPlIkK5iA/NkUhyy8+HBx7g4Kef8mRZ\nGWuGDeO/J0ygVatWQUcmImGaI5HktHw5PPAApQsWMKlLF/5zxQrGjR/P8OHDg45MJC1pjkTSx4YN\ncOutlFx0EX9Zt452Bw/ySf/+LFq6VEVEJEmpkKS5lBkD3raN8rvv5lDPnrzw/vv0a96cvbfdxsoN\nG/jd734X81BWyuQhAZSLCOUiPjRHIsHat4+Dv/0t/uyzvJaZyT969GD02LEsuuYaMrUWRCQlaI5E\ngnHoEDt++1vqPvEE75aW8tGVV/L9hx6if//+QUcmUiulzBdbJYoKSRIrK2P1I4/Q7MknKSgpYfn3\nv88Nv/kNbdu2DToykVpNk+1yXMkyBlxaUsKcsWNZe8op7Pnd75h9++0M3LmTn73wQkKKSLLkIRko\nFxHKRXxojkRq1L59+5hy//2c8/zznJWZyRc/+xkDHnmEfnX01hNJF1Ue2jKzBoS+FvdQzYQUOw1t\nJYdZL7/MgTFj6JORwf777qPjAw/oZooiSapGv7PdzAy4DhgFDAQMyDCzUiAfeBX4pz655Yg9O3bw\n3tChXL5wIbt+8ANOnziR0+vXDzosEakhlZkjmQX0BZ4AOrj7We5+BtAxvK0fMLPGIpSYJHoMeO64\ncWw+80zO27KFegsX0ukvf4EkKCIaC49QLiKUi/iozED15e5+OHpjeNs8YJ6Z1Yt7ZJJS9mzYwOLB\ng+m+ejVf3HcfPR59VHfjFaklqjRHYmanAp2BBke2ufuHNRBXTDRHkkDuLPnlLznjv/+bFeeeS+93\n3qFpu3ZBRyUiVZSQdSRmdhtwN5AFFAA5QL67X1qdA9ckFZLE2Ld4MRuvuoqMbdvY/8QT9P/pT4MO\nSUSqKVHrSO4G+gMb3P3bQG9gd3UOKolTI2PAhw6x+uabKenXj8K2bWmzdWvSFxGNhUcoFxHKRXxU\n5WL+g+5+0Mwws/ruvsLMzq2xyCQp7f/Xv9j3wx+y/vBhPn/5ZW4YNSrokEQkYFUZ2voHcCtwD3Ap\nsAuo6+5Day686tHQVg3Yto0to0bBjBm8ddll/OCNN2jatGnQUYlInCT8Xltm9i2gGfDusa7oCpoK\nSRyVl1M8fjyl993Ha3Xrcs4rr/Dtq68OOioRibManSMxswZmdo+ZPWdmY8ysjrvPcvfJyVhE5Oti\nGgNevZpd553Hyv/4D54ZOpSbNm9O2SKisfAI5SJCuYiPyky2/4XQosOlwBDgyVgOaGaDzWyFma0y\ns7HH6fOsma02swIzy456LcPMFpnZ5FjikBNwp3j8ePb16sW4rVvZ/tZbPPTmmxrKEpFjOunQlpkt\ndffzws/rAPPdvU+1DmaWAawCLgM+AxYAI9x9RYU+Q4C73H2YmV0APOPuORVe/xmhlfZN3f2a4xxH\nQ1vVtXs3W6+7jr15eUy67jp+9vzzKiAitUBNX/5bcuSJu5dW5yAVDABWu/sGdy8BJgHXRvW5Fngp\nfLx5QDMzaw1gZlnAUOBPMcYhx7D/vffY1rYtUxctYss//8nDmlAXkUqoTCE538z2mtk+M9sH9KrQ\n3lvF47UBNlVobw5vO1GfLRX6PA3cC+h0o5IqNQZcWsqam2+meOhQXh84kOs3b+bSYcNqPLZE0lh4\nhHIRoVzEx0nXkbh7Utz328yGAV+4e4GZDSJ0F+LjGj16NO3btwegefPmZGdnM2jQICDy5lF7EHuX\nLeOtSy5h9/79ZL/yCneOGJFU8cWrXVBQkFTxBNkuKChIqnjUDqZ95Pn69euJVVXWkfR194+jtl3l\n7m9X+mBmOcCv3X1wuH0/oe82ebxCnwnADHd/LdxeAXyL0Mr6HwClQEPgFODv7n7zMY6jOZJKKHjw\nQbIee4wP+/Xj8nffpemppwYdkogEJFG3SPmjmfWscNCRwENVPN4C4Bwzaxe+Y/AIIPrqq8nAzeFj\n5AC73f0Ld/9Pdz/b3TuG95t+rCIiJ7fns8/48Nxzaf744xQ98wzD581TERGRaqtKIbkBeMnMuprZ\n7cC/AVdU5WDuXgbcBUwFlgOT3L0wvD7ljnCfKcA6M1sD/CF8HKmmiqexAB+NH8+2du2o606L9esZ\ncNddwQSWYNF5qM2UiwjlIj4qfa8tdy8ysxHAP4GNwBXufqCqB3T3d4Fzo7b9Iap9wk83d59F6Au3\npJL27NrFu0OGcPmCBXw2diy5v/1t0CGJSJqo1DoSvn6V1OnAHuAQgLv3qrHoqklzJF83Y9IkuPVW\n2jVtyunvv0+TXkn3n0xEAlaj39kOXFWdHyzB27NnD3++6SZGffABxTfdRIcXX4S6dYMOS0TSTGXm\nSDaGFxAe8wFgpu9UTTZLFy3ijbZt6TRrFqdMnkyHV16p1UVEY+ERykWEchEflSkkM83sXjPrEv2C\nmXUJ3y9rZtwjk2qb/dZb7MrJYWinTjSdNIlGQ5PuTv8ikkYqM0dSD/g+MAroAewjtBiwCbAMeAX4\nazLdCbg2z5H8z1NP0fW++6h3/fW0nTQJMpNiPamIJLmEfR9J+KaLrcLN7e5eXp2D1rTaWkj+NmYM\ng/70Jw7+6ldkPfxw0OGISApJ1IJE3L3c3b8MP5KyiNRG5eXl/O3SSxn05z9TPmnS14qIxoBDlIcI\n5SJCuYiPqnxnuyShg/v28WHv3vTdupW68+fTrHfvoEMSkVqmWl+1m+xqy9DW7qIi1vTtS0aDBnQv\nKKBB69ZBhyQiKSphQ1tm9gszm2Zmy8zst2ZWe68nDdjWmTPZ3a0b+zt0IHvjRhUREQlMlQoJsNLd\nL3P3nsAHwIM1EJOcxLqJE6l72WWsuvZaBi1aRMYJ1odoDDhEeYhQLiKUi/ioaiE5w8yGmlkTd59O\n6G6+kkCr7r6bxj/5CZ888ABXvP560OGIiFT58t9HCK0juQBoCWQCE4Gsit8pErS0nCMpKWH1sGGU\nT5/OzhdfJPcHPwg6IhFJI4lcR9IbaOTueeF2R2AgcLu7f6s6AdSEtCsku3ax4YILKNq4kdOnT6fH\nwIFBRyQiaSaR60gWHyki4XaRu///wE3VObicXPmKFXzZqROzduzgnMLCKhcRjQGHKA8RykWEchEf\nVZ0jOSZ3/zweP0e+7vCUKezNzubFVq24es0a2nboEHRIIiLfoHUkSeqrp5/m4Nix/J+BA7n/3Xdp\n0KBB0CGJSBpL2NCWJMbu++9n59ix/H7UKH41fbqKiIgkNRWSZOLOl7ffzpdPPsk7v/wlD774IhkZ\nsf0n0hhwiPIQoVxEKBfxoXttJYvycjbdeCO7Jk9m7fjx3HHHHUFHJCJSKZojSQZlZRRdfjk75szh\n4N//zsVXXx10RCJSy9T0d7ZLTSopYcWAAWwvLOTUOXPof8EFQUckIlIlmiMJUHlxMcvPPZeta9fS\nbulSetRAEdEYcIjyEKFcRCgX8aEzkoAc3L6d1T168KU7fdas4dTTTw86JBGRatEcSQB2r1/P5vPP\n58sWLRi4bBkNGjcOOiQRqeW0jiSFbCkoYGu3buzs2JFBq1eriIhIylMhSaDCDz7gq/792XPRRVyy\naBEZdWp+ZFFjwCHKQ4RyEaFcxIcKSYLkv/oqDa+8kuIbbyTn/ffBqnUGKSKSdBI+R2Jmg4FxhIrY\n88f6HhMzexYYAnwFjHb3AjOrD3wI1As/3nL3/zzOMZJqjuR/nnyS8++7j+K77qLLM88EHY6IyDek\nzByJmWUAzwFXAj2AkWbWNarPEKCTu3cGxgATANz9EPBtd+8N9AIuNbMLExl/dbz8i1/Q77778Ece\nURERkbSU6KGtAcBqd9/g7iXAJODaqD7XAi8BuPs8oJmZtQ63i8N96hOKfVdCoq6G8vJyxo0YwdBn\nniFz/HjaPhjM19trDDhEeYhQLiKUi/hIdCFpA2yq0N4c3naiPluO9DGzDDNbDHwOzHT3T2sw1mo7\nePAg/99ll3HL3/9OvVdfpdVPfhJ0SCIiNSalFiS6eznQ28yaAlPN7FvuPutYfUePHk379u0BaN68\nOdnZ2QwaNAiI/F9ITbRLS0sZ06cPN65ZQ+O336beFVfU6PFO1h40aFCgx0+m9hHJEk9Q7SPbkiUe\n/X0E0z7yfP369cQqoZPtZpYD/NrdB4fb9wNeccLdzCYAM9z9tXB7BfAtd/8i6mc9BBS7+5PHOE4g\nk+3uzn9fdRV3fPABTd5/nzqXXJLwGEREqiNlJtuBBcA5ZtbOzOoBI4DJUX0mAzfD0cKz292/MLNW\nZtYsvL0h8B2gIHGhn9wffvxjbps6lfp//3vSFJHo/xuvrZSHCOUiQrmIj4QObbl7mZndBUwlcvlv\noZmNCb3sE919ipkNNbM1hC7/vTW8+5nAX8zMwvu+7O7TEhn/iUx64AG++9JLZL74Ig2HDQs6HBGR\nhNG9tuLgveeeo9fdd5Px1FO0vvvuhB1XRCReYhnaUiGJ0dzXX+eskSPhgQc4+ze/ScgxRUTiLZXm\nSNJK4cyZnDZqFAfvuCNpi4jGgEOUhwjlIkK5iA8VkmratGQJXHEFxddfT5ff/z7ocEREAqOhrWrY\nsX49m7sB17w1AAALg0lEQVR1ozQnh77Tp+sGjCKS8jRHEqUmC0nxjh0UduxIaYcOXLB4sYqIiKQF\nzZEkSGlxMcu7d+dQixb0X7AgJYqIxoBDlIcI5SJCuYgPFZJK8tJSFvfoQak7/ZcvJ6Nu3aBDEhFJ\nChraqozycj7u04eyoiK6FxXRpFWr+P1sEZEkEMvQVkrdtDEQ7nxy2WWwciXtP/1URUREJIqGtk6i\n8MYbyZgzh5b5+ZzeoUPQ4VSZxoBDlIcI5SJCuYgPFZITKPrJT8j85z/x996jfXZ20OGIiCQlzZEc\nx+YHH6TsscfY+MorXDxiRJwiExFJTrr8N862P/00mY89xidPPaUiIiJyEiokUfa++ip+7728f++9\nXJ0Gd/LVGHCI8hChXEQoF/GhQlLBwdmzKbvlFv46YgQ3/9d/BR2OiEhK0BxJWOmqVezp1Yu/5ORw\nz/TpZGSoxopI7aF7bUWpaiHx7dv5olMn/tamDWOWLKGuVq2LSC2jyfZYHDjAxt69ea9RI0bPn592\nRURjwCHKQ4RyEaFcxEftXtleVsbagQP5dPduhqxeTZMmTYKOSEQk5dTeoS13Vg8bxrZp0zijoICO\n3bolJjgRkSSke21Vw9q77qJ06lSaTJ+uIiIiEoNaOUey4YknaDBhAttfeolel1wSdDg1SmPAIcpD\nhHIRoVzER60rJJ+/8QaNx45l6aOPcvGoUUGHIyKS8mrVHMnu/HxKL76YD2+7jeETJgQQmYhIctI6\nkijHKiTFRUXs7t6d2Zdeyk1TpgQUmYhIctI6kpMo3b2bz/r0Ib9rV258++2gw0kojQGHKA8RykWE\nchEfaX/VlpeU8Ol55/FZkyZcM3++bn0iIhJn6T205c7CPn04XFREr/XraXLqqUGHJiKSlFJqaMvM\nBpvZCjNbZWZjj9PnWTNbbWYFZpYd3pZlZtPNbLmZLTWzn57sWPOuuYaGhYV0XrxYRUREpIYktJCY\nWQbwHHAl0AMYaWZdo/oMATq5e2dgDHDk8qpS4Ofu3gPIBe6M3reiBf/+75z5zjs0mTmT0zp2rIHf\nJjVoDDhEeYhQLiKUi/hI9BnJAGC1u29w9xJgEnBtVJ9rgZcA3H0e0MzMWrv75+5eEN6+HygE2hzv\nQO3Hj2f/a6/RLienJn4PEREJS3QhaQNsqtDezDeLQXSfLdF9zKw9kA3MO96BNjz5JN2/+90YQk0P\ngwYNCjqEpKA8RCgXEcpFfKTcJUxm1gT4G3B3+MzkmPr97GeJC0pEpBZL9OW/W4CzK7Szwtui+7Q9\nVh8zq0OoiLzs7m+d6ECjR4+mffv2ADRv3pzs7Oyj//dxZFy0NrQrjgEnQzxBtQsKCrjnnnuSJp4g\n2+PGjau1fw/R7dr893Hk+fr164lVQi//NbNMYCVwGbAVmA+MdPfCCn2GAne6+zAzywHGuXtO+LWX\ngO3u/vOTHKfKX7WbrmbOnHn0DVSbKQ8RykWEchGRUrdIMbPBwDOEhtWed/f/MrMxgLv7xHCf54DB\nwFfAaHdfbGYXAh8CSwEPP/7T3d89xjFUSEREqiClCkkiqJCIiFRNSi1IlMSqOB5amykPEcpFhHIR\nHyokIiISEw1tiYiIhrZERCQ4KiRpTmPAIcpDhHIRoVzEhwqJiIjERHMkIiKiORIREQmOCkma0xhw\niPIQoVxEKBfxoUIiIiIx0RyJiIhojkRERIKjQpLmNAYcojxEKBcRykV8qJCIiEhMNEciIiKaIxER\nkeCokKQ5jQGHKA8RykWEchEfKiQiIhITzZGIiIjmSEREJDgqJGlOY8AhykOEchGhXMSHComIiMRE\ncyQiIqI5EhERCY4KSZrTGHCI8hChXEQoF/GhQiIiIjHRHImIiGiOREREgpPwQmJmg81shZmtMrOx\nx+nzrJmtNrMCM+tdYfvzZvaFmX2SuIhTm8aAQ5SHCOUiQrmIj4QWEjPLAJ4DrgR6ACPNrGtUnyFA\nJ3fvDIwBfl/h5RfC+0olFRQUBB1CUlAeIpSLCOUiPhJ9RjIAWO3uG9y9BJgEXBvV51rgJQB3nwc0\nM7PW4fYcYFcC4015u3fvDjqEpKA8RCgXEcpFfCS6kLQBNlVobw5vO1GfLcfoIyIiSUKT7Wlu/fr1\nQYeQFJSHCOUiQrmIj4Re/mtmOcCv3X1wuH0/4O7+eIU+E4AZ7v5auL0C+Ja7fxFutwP+5e69TnAc\nXfsrIlJF1b38t068AzmJBcA54WKwFRgBjIzqMxm4E3gtXHh2HykiYRZ+HFd1kyEiIlWX0KEtdy8D\n7gKmAsuBSe5eaGZjzOyOcJ8pwDozWwP8Afi3I/ub2avAR0AXM9toZrcmMn4REfmmtFzZLiIiiZOy\nk+3VWNiYnegYE+VkuTCzUWa2JPyYY2bnBRFnIlTmfRHu19/MSsxseCLjS6RK/o0MMrPFZrbMzGYk\nOsZEqcTfSEszeyf8WbHUzEYHEGaNq8yi7mp9brp7yj0IFcA1QDugLlAAdI3qMwT4n/DzC4C5Qccd\nYC5ygGbh54Nrcy4q9JsGvA0MDzruAN8XzQgNMbcJt1sFHXeAuXgYeOxIHoAdQJ2gY6+BXFwEZAOf\nHOf1an1upuoZSUwLG9PMSXPh7nPdfU+4OZf0XZdTmfcFwL8DfwO+TGRwCVaZXIwC3nT3LQDuvj3B\nMSZKZXLxOXBK+PkpwA53L01gjAnhJ1/UXa3PzVQtJFrYGFGZXFR0G/BOjUYUnJPmwszOAq5z999z\nkqv/Ulxl3hddgBZmNsPMFpjZDxMWXWJVJhd/BHqY2WfAEuDuBMWWbKr1uZnoy38lQGb2beBWQqe3\ntdU4oOIYeToXk5OpA/QBLgUaA/lmlu/ua4INKxC/BJa4+7fNrBPwvpn1cvf9QQeWClK1kGwBzq7Q\nzgpvi+7T9iR90kFlcoGZ9QImAoPdPV3vV1aZXPQDJpmZERoLH2JmJe4+OUExJkplcrEZ2O7uB4GD\nZvYhcD6h+YR0UplcXAg8CuDua81sHdAVWJiQCJNHtT43U3Vo6+jCRjOrR2hhY/QHwWTgZji6oj56\nYWO6OGkuzOxs4E3gh+6+NoAYE+WkuXD3juFHB0LzJP+WhkUEKvc38hZwkZllmlkjQpOrhQmOMxEq\nk4tC4HKA8JxAF6AooVEmzokWdVfrczMlz0jcvczMjixszACe9/DCxtDLPtHdp5jZ0PDCxq8IDemk\nncrkAngIaAH83/D/iZe4+4Dgoq4ZlczF13ZJeJAJUsm/kRVm9h7wCVAGTHT3TwMMu0ZU8n3xGPCC\nmS0h9CF7n7vvDC7qmhFe1D0IaGlmGwldrVaPGD83tSBRRERikqpDWyIikiRUSEREJCYqJCIiEhMV\nEhERiYkKiYiIxESFREREYqJCIiIiMVEhERGRmKiQiNQgM6tnZrPCdxSoyj4fVmUfkSCpkIjUrO8D\nb3sVbiHh7oeBD4HraiwqkThSIRGpWaMI3RyR8E0DC83sBTNbaWavmNl3zCwv3O5XYb9/hfcVSXq6\n15ZINZlZJnAT0JHQlwENAJ5w93Xh1zOAze5+VrjdDlgNZLv7p2a2kNB3YPzYzK4BbnX368N96wHr\n3D0dv4xN0ozOSESqrxehW9EXEbpj7BvA1gqvtwL2Re2zrsIddpcDH4SfLyX0neLA0eEtM7MGNRC3\nSFypkIhUk7svDn/g5wKz3H1m+EuiKoqeMD9U4Xl5hXY53/xahwzS+Fb3kj5USESqycz6m1lLoIe7\nrzOz6K8w3g40id7tRD+yws+uB5S6+6ET9BdJCiokItU3GBgOfGRm1xE6qzjK3cuBZWbWpeLm4zyP\nbvcG8uMYq0iN0WS7SA0ys1uAM9z98Sru9yiw0N3/UTORicSPColIDQoPUb0PDKrsWpLq7CMSJBUS\nERGJieZIREQkJiokIiISExUSERGJiQqJiIjERIVERERiokIiIiIxUSEREZGYqJCIiEhM/h91+K9b\nbPRkmQAAAABJRU5ErkJggg==\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "#Plot the results\n", "%matplotlib inline\n", "import matplotlib.pyplot as plt\n", "\n", "problem = OneDimFEA(nodes, \n", " a=lambda x: c[0] + c[1] * x + c[2] * x ** 2, \n", " u0=100000, \n", " Q0=100)\n", " \n", "sol = problem.solve().flatten()\n", "\n", "plt.plot(nodes, p / 1000, 'k-', nodes, sol / 1000, 'r-');\n", "plt.xlabel('$x$ (m)')\n", "plt.ylabel('$p$ (kPa)')\n", "plt.grid()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The fit looks pretty good, so an acceptable answer is:\n", "\n", "$$\n", "\\kappa(x) = 0.96 + 0.58 x + 8.94 x^2\n", "$$" ] } ], "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.5.0" } }, "nbformat": 4, "nbformat_minor": 0 }