Skip to content

Instantly share code, notes, and snippets.

@BrianChevalier
Created March 21, 2017 06:44
Show Gist options
  • Save BrianChevalier/25473cdc2295088353c151d6e95559e2 to your computer and use it in GitHub Desktop.
Save BrianChevalier/25473cdc2295088353c151d6e95559e2 to your computer and use it in GitHub Desktop.
Notebook
Display the source blob
Display the rendered blob
Raw
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Theory and Background"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Problem Statement"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The vertical stress under the corner of a rectangular footing subjected to a uniform load of intensity q is given by the solution of Boussinesq’s equation,\n",
"\n",
"\\begin{equation}\n",
"\\sigma = \\frac{q}{4\\pi}\\left[ \n",
"\\frac{2mn\\sqrt{m^2+n^2+1}m^2+n^2+2}{m^2+n^2+1+m^2n^2+1}\n",
"+\\sin^{-1}\\left( \\frac{2mn\\sqrt{m^2+n^2+1}}{m^2+n^2+1+m^2n^2}\\right)\n",
"\\right]\n",
"\\label{Equation}\n",
"\\end{equation}\n",
"\n",
"\n",
"Where m and n are dimensionless ratios with $m=a/z$ and $n=b/z$ , where $a=4.6$ and $b = 14$ are the dimensions of the footing and z is the depth below the footing. If $a = 4.6$ and $b = 14$, use a third-order interpolating polynomial to compute the stress, $\\sigma$, at a depth $z = 10m$ below the corner of a footing that subjected to a total load of 100 tones. Note that you will need to use Boussinesq’s equation to generate enough data points over different depths to create an interpolated polynomial."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Solution Strategy\n",
"To solve this problem, Newton's Divided Difference Interpolation is used. Interpolation creates a continuous function, $f(x)$, of order $n$, that passes through $n+1$ points."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"To create a general equation to find the produce an interpolation we will look at the formula for a second order interpolation which is shown in the following equation:\n",
"\n",
"\\begin{align}\n",
"f_2(x)=b_0+b_1(x-x_0)+b_2(x-x_0)(x-x_1)\n",
"\\end{align}\n",
"\n",
"where:\n",
"\n",
"\\begin{equation}\n",
"b_0=f(x_0)\n",
"\\end{equation}\n",
"\n",
"\\begin{equation}\n",
"b_1= \\frac{f(x_1)-f(x_0)}{x_1-x_0}\n",
"\\end{equation}\n",
"\n",
"\\begin{equation}\n",
"b_2 = \\frac{\\frac{f(x_2)-f(x_1)}{x_2-x_1}-\\frac{f(x_1)-f(x_0)}{x_1-x_0}}{x_2-x_0}\n",
"\\end{equation}"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"There is a pattern to this polynomial which is given by:\n",
"\n",
"\\begin{align}\n",
"f_n(x)=b_0+b_1(x-x_0)+...+b_n(x-x_0)(x-x_1)...(x-x_{n-1})\n",
"\\end{align}\n",
"\n",
"where:\n",
"\n",
"\\begin{align}\n",
"b_0&=f[x_0]\\\\\n",
"b_1&=f[x_1, x_0]\\\\\n",
"b_2&=f[x_2, x_1, x_0]\\\\\n",
"\\vdots\\\\\n",
"b_{n-1}&=f[x_{n-1},x_{n-2},...,x_0]\\\\\n",
"b_n&=f[x_n,x_{n-1},...,x_0]\n",
"\\end{align}\n",
"\n",
"and finally the $n+1$ divided difference is given by:\n",
"\n",
"\\begin{align}\n",
"b_{n+1}&=f[x_{n+1},...,x_0]\\\\\n",
"&=\\frac{f[x_{n+1},...,x_1]-f[x_{n},...,x_0]}\n",
"{x_{n+1}-x_0}\n",
"\\end{align}\n",
"\n",
"This same process is followed to solve this problem."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Problem Solving Process"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The first step to solve this problem is to generate data so that a cubic interpolation function can be produced. A MATLAB function is defined to take in some inputs including a, b, z, and q. The code is as follows:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"``` matlab\n",
"function [sigma] = func(a, b, z, q)\n",
" \n",
" %define m and n\n",
" m = a./z;\n",
" n = b./z;\n",
" \n",
" %define some variables to break up the equation\n",
" var1 = q/(4*pi);\n",
" var2 = 2.*m.*n.*sqrt(m.^2+m.^2+1).*m.^2+n.^2+2;\n",
" var3 = m.^2+n.^2+1+m.^2.*n.^2.*m.^2+n.^2+1;\n",
" var4 = asin((2.*m.*n.*sqrt(m.^2+n.^2+1))./...\n",
" (m.^2+n.^2+1+m.^2.*n.^2));\n",
" \n",
" %compute sigma based on Boussinesq’s equation\n",
" sigma = var1.*(var2./var3+var4);\n",
" \n",
"end %function\n",
"```"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Next, the constants must be defined. The variables a, and b are given in the problem statement. The variable q is given in metric tons which is a mass and not a force so it must be converted:\n",
"\n",
"\\begin{align}\n",
"q&=\\frac{F}{A}\\\\\n",
"&=\\frac{(100Mg)(9.81m/s^2)}{(4.61m)(14m)}\\\\\n",
"&=\\frac{981kN}{64.4 m^2}\\\\\n",
"&=15.2329kPa\n",
"\\end{align}\n",
"\n",
"The code to initialize these variables is as follows:"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"a = 4.6;\n",
"b = 14;\n",
"q = 981/(a*b);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Next the function is plotted to visualize the the function. Note that that the variable z in the scatter command has a mius sign in front. This is to show that z is a depth into the ground."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAjAAAAGkCAIAAACgjIjwAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAA\nB3RJTUUH4QMVBRAeh61SIwAAACR0RVh0U29mdHdhcmUATUFUTEFCLCBUaGUgTWF0aFdvcmtzLCBJ\nbmMuPFjdGAAAACJ0RVh0Q3JlYXRpb24gVGltZQAyMC1NYXItMjAxNyAyMjoxNjozMJVvuZsAACAA\nSURBVHic7d1/XFR1vsfx7wjq6AbXQSF+pPIjQ3NXb9KokSaOStbqkntTcJXEHvvwei1l7bFSN72x\ntBVd9aqYeuve3QWvabh5t9SHmtYS/mqJH+6tNEAZGWUUC2RggnFQYO4fx51YQJJfc74Mr+ejP873\nO2fmfL6gvf1+z5lzNA6HQwAAoLZ+ahcAAIAQBBIAQBIEEgBACgQSAEAKBBIAQAoEEgBACgQSAEAK\nBBIAQAoEEgBACgQSAEAKBBIAQAoEEiCj69evnzhx4tKlS61fampqOnHixFdffeX6qoAeRSABMlq5\ncuW0adMKCwtbv9SvX7/nnntuypQpdXV1ri8M6DkEEiCd8+fPZ2Zm+vr6zpw5s80dFixYYLVa3333\nXRcXBvQoAgmQzv79+x0Ox/z58z09PdvcYcGCBUKI7du3u7YuoGcRSIA4fPjwpFYuXrzYfJ+lS5dO\nmjRp7969SvOdd96ZNGnSv//7vwshDh48OHny5HvuucfLy2vq1KmffvppF+s5cuSIEOKnP/3pnQoL\nDw8PCwv76quvvv766y4eC5BH2//+AvqUhoYGm82mbJtMptraWiHEzZs3m+8TFhaWkZGxa9eu2NhY\nIcQf/vCH3Nzc5OTkixcvzps3T6vVJiUlVVVVvfXWW0888YTRaAwKCvrB43799dd5eXkNDQ1K8623\n3lq/fv2sWbMKCgqEEMHBwSUlJXcqLDg42Gg0FhUVPfjgg931cwBU5gDwN1lZWQMGDBBCrFy5ssVL\n58+fF0Jotdra2tqrV69qNBqdTnfz5s19+/YJIYKCgkwmk8Ph+J//+Z/f/va3X3/99Q8eKzExUQjx\n2muv/dM//ZMQYuLEiT/5yU9OnDhRVlam/N2sqalpp7AlS5YIIVJTU7tv9IDKmCEBt509e3bevHk3\nb96cN2/eli1bWrw6atSoCRMmnDlz5ujRo5WVlQ6HY968ef3794+MjPT3979y5UpISMiDDz44bdq0\nhISEMWPGtH+stWvXpqWlffvtt76+vkIIjUbTr1+/L7/8UgjxxRdfCCG8vLy8vb3bKey+++4TQpSU\nlHTrzwBQE+eQACGEMJvNTzzxRE1NTWRk5J49e/r1a+OvhrJYd+DAgf3794u/XVkQEBDw9ddfv/XW\nWz/96U8vXbq0Y8eOyZMnHz16tJ1jlZWVvfHGG+vWrVPSSOG8hvtHP/qREMJutzscjnYKcy70AW6D\nQAJETU3Nk08+aTabH3jggQMHDmi12jZ3UxLo4MGDWVlZQ4cOnTFjhhAiKytrw4YNDzzwwMGDBysr\nKx999NGmpqaPPvrIbrd/+OGHH374YWVlZYvP+eyzz4QQTz31lNJUFgNXrFihNENCQvr373/r1q1v\nv/22ncIuX74shGgeaUBvx5IdINasWaPc+GDkyJEvvfSS0jl//vzo6OjmuwUHB0+aNOnzzz8XQixe\nvFi5JruxsTE1NfXdd9/dvHnz4MGDq6qqhBATJkywWCzz5s0TQmRlZU2fPr355wwdOrR5Mzw8PDw8\nfPny5UrTw8MjNDS0uLjYbDa/8847dypMOdU0bty4bv9pAKpR+yQWoL6nn3669V+NDRs2tN5z8+bN\nyqsff/yxs3Pjxo0+Pj5K/+DBg1944YWmpqarV68qPVlZWa0/RwgRGBh46NAhIcQvf/nLFq/GxcUJ\nIXbt2nWnwhobG/38/IQQpaWl3fqTANSkcTgc3Z9yQB9z69atsrKyW7dujRgxYtCgQXfzlps3b373\n3XctZkuKU6dOTZ069Yknnjh8+HCb7/3kk09mzZo1adKknJycLtUNyIQlO6Ab9O/fPzQ0tENvGTBg\nQJtpJISYMmXKxIkTjx07du3aNX9//9Y77Ny5UwixatWqTpQKSIuLGrrKZrOdP3/earWqXQjcym9+\n85tZs2b93//9X+uXmpqabty4kZCQsHDhQtcXBvQcluy6ZN++fevXrw8KCrp48eKCBQvWrl2rdkVA\nr5dttBwvqVa2R/poE/QB6tYDl2HJrvNsNtsrr7yye/fuhx566Nq1azNnznzyyScfeughtesCerGM\nvPKdeeVRYTqlmXKsVAhBJvURBFLnaTQaT0/Pe++9Vwjh7e3dv39/5eYuADot5VhpetwYZyBNu3/I\n0sxCAqmPYMmuS9599919+/ZFR0efPHly1KhRKSkpGo1G7aKA3irlaKnJYk+P+7sbLy3NLBRCtOiE\nW2KG1AFWq9V548vAwEBvb+/Lly/X1NRcvXq1vr7ebDZXVVW1vm4qPj4+NzdXCLF48eL4+HhXF93d\nrFar8x5rvZ07jUW4xXB+l1O2YaavyWRqPpbHh4u0zy0mk0nV0rpEo9GMHDlS7Sp6AQKpA/Lz81NT\nU5XtxMREPz+/Q4cOffTRR15eXg6HY8mSJXv37nXeAMYpNze3uLjY5cX2FJPJFBwcrHYV3cOdxiLc\nYji/fXLgUWN13KPBzcey9PCZ2WMDevXQenWauhKB1AEGg8FgMDibu3fvDgwM9PLyEkJoNJqIiAij\n0ahedUCvFxWmSzlWmm20BHvc7sk2WkwWe/LjIarWBRchkDpv/Pjxv/3tb0+fPv3oo49WVVUdOXLk\nmWeeUbsooBcL9tEmR4cszSyMuX+QrtghhMg2Wjh71HcQSJ334x//eO3atYmJiVqt1mazxcTEKLcg\nA9BpCfqAYB/th/kmpblEH+C84g5uj0Dqkvj4+EWLFlVWVup0uv79+6tdDuAOosJ0wR41vfqkETqH\nQOqqfv36KfddBgB0BfeyAwBIgUACAEiBQAIASIFAAgBIgUACAEiBQAIASIFAAgBIge8hAehBzrvd\nu6WJEyfu2rVL7SrcB4EEoAe52d3uWwgPD1e7BLfCkh0AQAoEEgBACgQSAEAKBBIAfM9sNtvt9uY9\njY2NpaWlatXTp3BRAwDZZRstx0uqle2RPtoEfUAPHaiuri42NjY7O7t5p4eHxy9/+cv09PQRI0b0\n0HGhYIYEQGoZeeUpR7+foKQcK83IK++hY7355pvPPPNM62ebrVq1at26dT10UDgxQwIgtZRjpelx\nY5zPjZ12/5ClmYU9MUlqaGh46623vv322x07dmzbts3Zv2XLlp/97GfPP/98RUWFr69vtx8XTgQS\nAHmlHC2NCtM1f4q50lyaWZgeN6aLH97U1FRUVHT9+nUhRFBQUFNTk4+Pz4ABA+bPnz9t2jQhRHZ2\n9ubNmydPnqzRaEaMGHHq1Kl58+Z18aBoB4EEQF4Z+eWtg2eJ3r/5Il7n/OUvf3n99dfvvffe999/\nf/78+ZGRkQEBAcOHDxdC+Pr6+vr6VlRUbNy48f333/f29hZCDB8+vKSkpIsHRfsIJADySo4O2Zl3\nrfkMSfxt2tSVj71+/forr7xy5MgRT09PT0/PX/ziF9OmTXv77bedK3IOh2PJkiUrV658+OGHlR5/\nf/+ampquHBQ/iIsaAMgrKkyXbbRkGy3OnmyjxWSxJz8e0pWPfe+99xISEjw9PYUQNptNySEvL69v\nvvlG2WHDhg0ajWb16tXOt1RWVirzJ/QcZkgA5BXso02ODlmaWZjw8O2rGLKNlq6fPbp8+fLo0aOF\nEBaL5fr16w8++KAQIiwsrKysTAiRk5OzY8eO/Px8jUbT/C3cua6nEUgApJagDwj20Tq/h7REH9DF\n9TohxPLly3fu3CmE+PTTT9PT05XOhx9+uLa21m63p6SkVFRUjBs3Tun/9a9//cILL5hMpjFjuhqE\naB+BBEB2LS6067rQ0NDk5GS73T5z5kxnp6en5+LFi3ft2nXkyJEW+3/88cePPfbYvffe2401oDXO\nIQHoi/r16zd48OAWnUlJSf/93//tcDha9G/dujUlJcVVpfVdBBIA3BYYGPgf//EftbW1zTsbGxtX\nrVoVFhamVlV9B0t2APC9qVOntujx8PCYNWuWKsX0NcyQAABSYIYEoAdNnDjRja+WnjhxotoluBUC\nCUAP2rVrVyfeZTKZgoODu7sWyI4lOwCAFAgkAIAUWLL7ATabzWw2+/v7K3f8beHmzZv19fXO5j33\n3NP8XiMAgLtHILVn375969evDwoKunjx4oIFC9auXdtih7S0tL1792q1WqX5wQcf8PwuAOgcAumO\nbDbbK6+8snv37oceeujatWszZ8588sknH3rooeb7FBUVbdmyZcqUKWoVCQBug3NId6TRaDw9PZW7\nV3l7e/fv33/AgAEt9ikqKho9erTZbLbb7WrUCADugxnSHQ0aNCgpKWnFihXR0dEnT56cO3euco96\np+vXr9fU1CQkJNhstoqKimXLlq1cubLNj1K+h7F48eL4+HhXlN6TzGaz2iV0G3cai3Cv4bjTWIQQ\nVqtV7RJ6BwLpe1arVXkaihAiMDDQ29v78uXLNTU1V69era+vN5vNVVVVQ4cOde5fU1Mze/bspKQk\nPz+/c+fOLVq0KCIiIjIysvUnFxcXu2gMLuFOXxBxp7EI9xqOO43FZDKpXULvQCB9Lz8/PzU1VdlO\nTEz08/M7dOjQRx995OXlpTzPeO/evStWrHDuHxoaunHjRmV77Nix0dHROTk5bQYSIJtso8X5hKGR\nPtoEfYC69QCCQGrOYDAYDAZnc/fu3YGBgV5eXkIIjUYTERFhNBqb75+bm1tdXR0dHa00m5qaGhsb\nXVkw0DkZeeU788qdTxhKOVYqhCCToDouarij8ePHf/XVV6dPnxZCVFVVHTlyJCIiQgixZ8+egoIC\nIYTyZEllsbukpCQrK2vGjBnq1gzcjZRjpcmPhzj/S48bo2QSoC5mSHf04x//eO3atYmJiVqt1maz\nxcTExMXFCSG2bt0aGxsbERHx2GOPxcbGxsTE+Pj41NXVJSUlTZgwQe2qgR+QcrS0xQNYlebSzML0\nOB7RDTURSO2Jj49ftGhRZWWlTqfr37+/0pmTk+PcYdWqVStWrLBYLHwfFr1FRn556+BZovdPOcok\nCSpjye4H9OvXz8/Pz5lGrXl6epJG6EWSo0N25l1r0alMm1SpB3AikIC+JSpMl220ZBstzp5so8Vk\nsSc/HqJiVYBgyQ7oa4J9tMnRIUszCxMevn1ZXbbRwtkjyIBAAvqcBH1AsI/W+T2kJfoA1usgAwIJ\n6ItaXGgHyIBzSAAAKRBIAAApEEgAACkQSAAAKRBIAAApEEgAACkQSAAAKRBIAAApEEgAACkQSAAA\nKRBIAAApEEgAACkQSAAAKRBIAAApEEgAACkQSAAAKRBIAAApEEgAACkQSAAAKRBIAAApEEgAACkQ\nSAAAKRBIAAApEEgAACkQSAAAKRBIAAApEEgAACkQSF1ls9kuXLhgs9nULgQAejcCqUt27NgRFRWV\nlJT02GOP7d69W+1yAKAX81S7gF7s9OnTu3btOnz48LBhwy5evDh//vzJkyeHhYWpXRcA9ErMkDqv\noKDgkUceGTZsmBAiNDR03LhxJ0+eVLsoAOitmCF13pAhQ65evapsOxyOa9euOZsthIeHCyEWL14c\nHx/vuvp6htlsVruEbuNOYxHuNRx3GosQwmq1ql1C70AgdYDVai0rK1O2AwMDZ8+evWnTpk2bNk2f\nPv3AgQN1dXX19fVtvrG4uNiFZfa44OBgtUvoNu40FuFew3GnsZhMJrVL6B0IpA7Iz89PTU1VthMT\nE+fMmfPee++98847mzZtmjVrlsFg8PHxUbdCAOi9CKQOMBgMBoPB2ayurm5qatqyZYvSjIuLc4MV\nOQBQCxc1dJ6Hh8ezzz5bWFgohDh16tSFCxeaxxUAoEOYIXWel5fXmjVrEhMT6+vrm5qaduzYMWjQ\nILWLAoDeikDqkqeffvrpp5+2WCw6nU7tWgCgd2PJrhuQRgDQdQQSAEAKLNkBLpVttBwvqVa2R/po\nE/QB6tYDyIMZEuA6GXnlKUdLnc2UY6UZeeUq1gNIhRkS4Dopx0rT48ZEhd0+6Tjt/iFLMwuZJAEK\nZkiAi6QcLY0K0znTSAihNJdmFqpYFSAPAglwkYz88iV6/xadS/T+pqobqtQDyIZAAlwkOTpkZ961\nFp3KtEmVegDZEEiAi0SF6bKNlmyjxdmTbbSYLPbkx0NUrAqQBxc1AC4S7KNNjg5ZmlmY8PDtqxiy\njZb0uDHqVgXIg0ACXCdBHxDso3V+D2mJPoD1OsCJQAJcqsWFdgCcOIcEAJACgQQAkAKBBACQAoEE\nAJACgQQAkAKBBACQAoEEAJACgQQAkAKBBACQAoEEAJACgQQAkAKBBACQAoEEAJACgQQAkAKBBACQ\nAoEEAJACgQQAkAKBBACQAoEEAJACgQQAkIKn2gXIzmazXblyJSgoaPDgwa1fvXnzZn19vbN5zz33\naDQaF1YHAO6DQGrPjh07MjIygoKCysrKVq9evWjRohY7pKWl7d27V6vVKs0PPvjA19fX5WUCgDsg\nkO7o9OnTu3btOnz48LBhwy5evDh//vzJkyeHhYU136eoqGjLli1TpkxRq0gAcBucQ7qjgoKCRx55\nZNiwYUKI0NDQcePGnTx5ssU+RUVFo0ePNpvNdrtdjRoBwH0wQ7qjIUOGXL16Vdl2OBzXrl1zNhXX\nr1+vqalJSEiw2WwVFRXLli1buXJlmx8VHh4uhFi8eHF8fHxPl93TzGaz2iV0G3cai3Cv4bjTWIQQ\nVqtV7RJ6BwLpe1artaysTNkODAycPXv2pk2bNm3aNH369AMHDtTV1TW/fkEIUVNTM3v27KSkJD8/\nv3Pnzi1atCgiIiIyMrL1JxcXF7tiAK4SHBysdgndxp3GItxrOO40FpPJpHYJvQNLdt/Lz8//1d+c\nPn3az8/vvffeu3z58qZNm0JCQgwGw9ChQ5vvHxoaunHjRj8/PyHE2LFjo6Ojc3JyVKodAHo9Zkjf\nMxgMBoPB2ayurm5qatqyZYvSjIuLa7HglpubW11dHR0drTSbmpoaGxtdVi0AuBlmSHfk4eHx7LPP\nFhYWCiFOnTp14cIFJa727NlTUFAghLDb7SkpKcpid0lJSVZW1owZM9StGQB6L2ZId+Tl5bVmzZrE\nxMT6+vqmpqYdO3YMGjRICLF169bY2NiIiIjHHnssNjY2JibGx8enrq4uKSlpwoQJalcNAL2VxuFw\nqF2D7CwWi06nu9OrDQ0NFoulne/DhoeHu9NFDSaTyW3ONrvTWIR7DcedxiLcbjg9hyW7H9ZOGgkh\nPD09uTsDAHQdgQQAkAKBBACQAoEEAJACgQQAkAKBBACQAoEEAJACgQQAkAJ3akDflXPlxs7iUmV7\npI82QR+gbj1AH8cMCX1URl552ucWZzPlWGlGXrmK9QBghoQ+KuVYaWqULu7REKU57f4hSzMLmSQB\nKmKGhL4o5WhpVJhuctAgZ09UmC4qTLc0s1DFqoA+jkBCX5SRX75E79+ic4ne31R1Q5V6AAgCCX1T\ncnTIzrxrLTqVaZMq9QAQBBL6pqgwXbbRknPl+/lQttFistiTHw9RsSqgj+OiBvRFwT7a5OiQNYcv\nFNdqlZ5soyU9boy6VQF9HIGEPipBH6C9WV1ce7u5RB/Aeh2gLgIJfdfkoEFxPMcTkAbnkAAAUiCQ\nAABSIJAAAFIgkAAAUiCQAABSIJAAAFIgkAAAUiCQAABSIJAAAFIgkAAAUiCQAABSIJAAAFIgkAAA\nUiCQAABSIJA6pra2tqysrHmP3W6/cOGC1WpVqyQAcA9uEkjbtm07fPhw855vvvlmxowZPXGg7du3\nO5sHDx6cOnXqr3/96+nTp2/YsKHbDwcAfYebBFJOTs7q1as3b97scDiUnsbGRrPZ3I2HSEtLW7hw\nYXp6urOnurp63bp127Zt279//0cffbRv377PPvusG48IAH2KmwSSEGLZsmVHjhx57rnnbDZbT3x+\nZGTkc8899/Of/9zZk5ub6+/vP2nSJCGEr6+vwWA4ceJETxwaAPoC93mEeVhY2L59+1atWrVgwYK3\n3367X79uzlq9Xi+EOHv2rMlkUnq++eabgIAA5w7+/v6XLl1q873h4eFCiMWLF8fHx3dvVa7XvfNO\ndbnTWIR7DcedxiKE4BzzXXKfQBJCeHt7//73v3/99deffvrpF198sSsfZbVanRcvBAYG6nS61vs0\nNDR4eHg4m56enrdu3Wrz04qLi7tSjGyCg4PVLqHbuNNYhHsNx53G4vxXLNrnVoEkhPDw8HjllVdG\njRq1bt26rnxOfn5+amqqsp2YmDhnzpzW+wwcONButzubN27cGDhwYFcOCgB9mZsE0vLly0eMGOFs\nLly4MDQ0dPfu3Z3+QIPBYDAY2t8nKCio+T98SktLR40a1ekjAkAf5yYXNUyZMqV5IAkhJk2atHXr\n1h49qF6vr6+vz8zMFEKcO3fu+PHj06dP79EjAoAbc5MZkioGDx68fv36l19+eevWrd99993zzz8/\nfvx4tYsCgN6KQOqY5cuXN28aDIa//OUvlZWVOp3O05MfJgB0Hv8P7SqNRuPr66t2FQDQ67nJOSQA\nQG9HIAEApEAgAQCkQCABAKRAIAEApEAgAQCkQCABAKRAIAEApEAgAQCkQCABAKRAIAEApEAgAQCk\nwM1VIZ1so+V4SbWyPdJHm6APULceAK7BDAlyycgrTzla6mymHCvNyCtXsR4ALsMMCXJJOVaaHjcm\nKkynNKfdP2RpZiGTJKAvYIYEiaQcLY0K0znTSAihNJdmFqpYFQDXIJAgkYz88iV6/xadS/T+pqob\nqtQDwJUIJEgkOTpkZ961Fp3KtEmVegC4EoEEiUSF6bKNlmyjxdmTbbSYLPbkx0NUrAqAa3BRAyQS\n7KNNjg5ZmlmY8PDtqxiyjZb0uDHqVgXANQgkyCVBHxDso3V+D2mJPoD1OqCPIJAgnRYX2gHoIziH\nBACQAoEEAJACgQQAkAKBBACQAoEEAJACgQQAkAKBBACQAoEEAJACX4z9O7W1tRaLZfjw4c4eu91e\nVlZ27733ent7t97/5s2b9fX1zuY999yj0WhcUSgAuB0C6e9s27aturr6zTffVJoHDx589dVXAwMD\nzWZzXFzcmjVrWuyflpa2d+9erVarND/44ANfX1+XVgwA7oJAui0tLS0nJ+fMmTPz5s1Teqqrq9et\nW/df//VfkyZNqqiomDNnzqOPPhoZGdn8XUVFRVu2bJkyZYoaJQOAW+Ec0m2RkZHPPffcz3/+c2dP\nbm6uv7//pEmThBC+vr4Gg+HEiRMt3lVUVDR69Giz2Wy3211aLgC4HWZIt+n1eiHE2bNnTSaT0vPN\nN98EBAQ4d/D397906VLzt1y/fr2mpiYhIcFms1VUVCxbtmzlypUuLBkA3EofDSSr1VpWVqZsBwYG\n6nRt3Fu6oaHBw8PD2fT09Lx161bzHWpqambPnp2UlOTn53fu3LlFixZFRES0WNNThIeHCyEWL14c\nHx/fncNQg9lsVruEbuNOYxHuNRx3GosQwmq1ql1C79BHAyk/Pz81NVXZTkxMnDNnTut9Bg4c2Hwh\n7saNGwMHDmy+Q2ho6MaNG5XtsWPHRkdH5+TktBlIxcXF3Va6BIKDg9Uuodu401iEew3HncbiXHdB\n+/poIBkMBoPB0P4+QUFBzf8YlZaWjho1qvkOubm51dXV0dHRSrOpqamxsbG7KwWAvoKLGu5Ir9fX\n19dnZmYKIc6dO3f8+PHp06crL+3Zs6egoMBut6ekpChrCyUlJVlZWTNmzFCzYgDozfroDOluDB48\neP369S+//PLWrVu/++67559/fvz48cpLW7dujY2NXb16dWxsbExMjI+PT11dXVJS0oQJE9StGQB6\nL43D4VC7Bqk5HI7KykqdTufp2XZ4NzQ0WCyWdr4PGx4e7k7nkEwmk9ss7rvTWIR7DcedxiLcbjg9\nhxnSD9BoNO3ffMHT05O7MwBA13EOCQAgBQIJACAFAgkAIAUCCQAgBQIJACAFAgkAIAUCCQAgBQIJ\nACAFAgkAIAUCCQAgBQIJACAFAgkAIAUCCQAgBQIJACAFAgkAIAUCCQAgBQIJACAFAgkAIAUCCQAg\nBU+1C0BvlW20HC+pVrZH+mgT9AHq1gOgt2OGhM7IyCtPOVrqbKYcK83IK1exHgBugBkSOiPlWGl6\n3JioMJ3SnHb/kKWZhUySAHQFMyR0WMrR0qgwnTONhBBKc2lmoYpVAejtCCR0WEZ++RK9f4vOJXp/\nU9UNVeoB4B4IJHRYcnTIzrxrLTqVaZMq9QBwDwQSOiwqTJdttGQbLc6ebKPFZLEnPx6iYlUAejsu\nakCHBftok6NDlmYWJjx8+yqGbKMlPW6MulUB6O0IJHRGgj4g2Efr/B7SEn0A63UAuohAQie1uNAO\nALqIc0gAACkQSAAAKRBIHVNbW1tWVnY3nQCADiGQOmbbtm3bt2+/m04AQIcQSHcrLS1t4cKF6enp\nP9gJAOgErrK7W5GRkREREYcOHXI4HO13AgA6gUC6W3q9Xghx9uxZk8nUficAoBMIpLZZrVbndQqB\ngYE6XZe+cBMeHi6EWLx4cXx8fDcUpyqz2ax2Cd3GncYi3Gs47jQWIYTValW7hN6BQGpbfn5+amqq\nsp2YmDhnzpyufFpxcXF3FCWL4OBgtUvoNu40FuFew3GnsbCCcpcIpLYZDAaDwaB2FQDQh3CVHQBA\nCgQSAEAKLNl1zPLly++yEwDQIcyQAABSIJAAAFIgkAAAUiCQAABSIJAAAFIgkAAAUiCQAABSIJAA\nAFIgkAAAUiCQAABSIJAAAFIgkAAAUiCQAABSIJAAAFIgkAAAUiCQAABSIJAAAFIgkAAAUiCQAABS\nIJAAAFIgkAAAUiCQAABSIJAAAFIgkAAAUiCQAABSIJAAAFIgkAAAUvBUuwC4TrbRcrykWtke6aNN\n0AeoWw8ANMcMqa/IyCtPOVrqbKYcK83IK1exHgBogRlSX5FyrDQ9bkxUmE5pTrt/yNLMQiZJAOTB\nDKlPSDlaGhWmc6aREEJpLs0sVLEqAGiOQOoTMvLLl+j9W3Qu0fubqm6oUg8AnZ+mJAAADadJREFU\ntEYg/Z3a2tqysrK76VTcvHnzu2YcDkfP19gZydEhO/OutehUpk2q1AMArRFIf2fbtm3bt2+/m05F\nWlra9OnTn/ibysrKnq+xM6LCdNlGS7bR4uzJNlpMFnvy4yEqVgUAzXFRw21paWk5OTlnzpyZN29e\n+53NFRUVbdmyZcqUKa4qs5OCfbTJ0SFLMwsTHr59FUO20ZIeN0bdqgCgOQLptsjIyIiIiEOHDjVf\ndmuzs7mioqLRo0ebzeZhw4ZptVpXFdsZCfqAYB+t83tIS/QBrNcBkAqBdJterxdCnD171mQytd/p\ndP369ZqamoSEBJvNVlFRsWzZspUrV7qo3E5pcaEdAEiljwaS1Wp1XqcQGBio03Xmf9M1NTWzZ89O\nSkry8/M7d+7cokWLIiIiIiMjW+8ZHh4uhFi8eHF8fHxXypaB2WxWu4Ru405jEe41HHcaixDCarWq\nXULv0EcDKT8/PzU1VdlOTEycM2dOJz4kNDR048aNyvbYsWOjo6NzcnLaDKTi4uJOlyqh4OBgtUvo\nNu40FuFew3GnsbS5xILW+mggGQwGg8HQxQ/Jzc2trq6Ojo5Wmk1NTY2NjV0uDQD6KC777ow9e/YU\nFBTY7faUlBRlbaGkpCQrK2vGjBlqlwYAvVUfnSF10datW2NjY1evXh0bGxsTE+Pj41NXV5eUlDRh\nwgS1SwOA3koj7c0FeouGhgaLxeLr63unHcLDw93pHJLJZHKbxX13Gotwr+G401iE2w2n57Bk11We\nnp7tpBEA4C4RSAAAKRBIAAApEEgAACkQSAAAKRBIAAApEEgAACkQSAAAKRBIAAApEEgAACkQSAAA\nKRBIAAApEEgAACkQSAAAKRBIAAApEEgAACkQSAAAKRBIAAApEEgAACkQSAAAKRBIAAApEEgAACkQ\nSAAAKRBIAAApEEgAACkQSAAAKRBIAAApEEgAACkQSAAAKRBIAAApEEgAACkQSAAAKRBIHVNbW1tW\nVta8x263X7hwwWq1qlUSALgHAqljtm3btn37dmdzz549U6dOffHFFw0GQ1paWs8dN9toSTlaqvyX\nkVfecwf6Qbt27VLx6N3LncYi3Gs47jQW4XbD6TkE0t1KS0tbuHBhenq6s+fChQsbNmzYt2/fn/70\npwMHDuzatevMmTM9ceiMvPKUo6XOZsoxNTPp3XffVevQ3c6dxiLcazjuNBbhdsPpOZ5qF9BrREZG\nRkREHDp0yOFwKD1FRUWRkZEjR44UQgQGBo4YMeLSpUsTJkzo9kOnHCtNjxsTFaZTmtPuH7I0szBB\nH9DtBwIAFTFDult6vX7KlClK/Cjmzp3rXL4zmUxGo3H8+PHdftyUo6VRYTpnGgkhlObSzMJuPxYA\nqIgZUtusVqvz4oXAwECdTtfOzgUFBatXr16+fHloaGjrVydOnBgeHt7pSkoNr937xf+EJ59v3mkb\n+sD1B376WfJTnf7YrujKcGTjTmMR7jUcdxrLxIkT1S6hdyCQ2pafn5+amqpsJyYmzpkzp83dbt26\ntXHjxsOHD//bv/1bdHR0m/t08XxmRl75ceOP0+PGNO+cvuPML8N0yRnLu/LJACAVAqltBoPBYDD8\n4G4rV6709PQ8dOiQt7d3D1USFaZLOVaabbQ4V+2yjRaTxZ78eEgPHREAVEEgdd4nn3xiNpv379/v\n4eHRc0cJ9tEmR4cszSxMePj2VQzZRkuLCRMAuAECqfNOnz7d4kKGTZs23WnhrisS9AHBPtrjJdVK\nc4k+oPk1DgDgHjTOi5gBAFARl30DAKRAIAEApEAgdbPWd19V3Lx587tmestK6Z2GoygvL7fZbK6s\np4vaHE5DQ8N3f89ut6tSXoe086ux2Wznz5/vXTf8bX84Fy5c6F1/0py4+XKHcFFDN9u2bVt1dfWb\nb77Zoj8tLW3v3r1arVZpfvDBB76+vi6vrsPuNJxLly4tX768sbGxrq5u+vTpr732mirldVSbw/nk\nk09eeuklZ/PmzZvz5s17/fXXXV5dx9zpV7Nv377169cHBQVdvHhxwYIFa9euVaW8jrrTcHbs2JGR\nkREUFFRWVrZ69epFixapUl7n7NmzZ/PmzcOHD798+XJ8fHxiYqLaFUnPgW6yZcuWuLi4Bx544MUX\nX2z96rPPPnvy5EnXV9Vp7Qynqalp1qxZBw4ccDgcN27cmDFjRl5enho1dkD7vx2nioqKmJiYixcv\nuqywTmhnLHV1dWPGjDlz5ozD4SgvLx87dqyyLbN2hnPq1KnJkydXVFQ4HA6j0ThhwoSSkhI1auyM\n8+fP/+M//qPJZHI4HFeuXImIiCgoKFC7KNkxQ+o2re++2lxRUdHo0aPNZvOwYcOc8ySZtTOcM2fO\naDSauXPnOhwOrVZ77NgxjUajSpF3r/3fjqKxsXH58uUvvPBCSIjUXzpuZywajcbT0/Pee+8VQnh7\ne/fv33/AgAFq1NgB7QynoKDgkUceGTZsmBAiNDR03LhxJ0+eDAsLU6PMDnPZzZfdCYHUbfR6vRDi\n7NmzJpOpxUvXr1+vqalJSEiw2WwVFRXLli1buXKlCiV2RDvDOX/+fHh4+G9+85vDhw/369dv4cKF\nq1atUqHEjmhnOE5//OMfvby8HnvsMdeV1SntjGXQoEFJSUkrVqyIjo4+efLk3LlzH3zwQRVK7Ih2\nhjNkyJCrV68q2w6H49q1a86m/ObOnTt37lxlu+duvuxmuKjBFWpqambPnv2HP/whKysrMzPz97//\n/WeffaZ2UZ1XXV2dlZXl7+9/8uTJ9PT0P/7xjwcOHFC7qK6y2+1bt26V/x8K7WtsbLx8+XJNTc3V\nq1fr6+vNZnNVVZXaRXXe7Nmzi4qKNm3a9Ne//vXVV1+tq6urr69Xu6gOKygoeOaZZ+5082U0RyC5\nQmho6MaNG/38/IQQY8eOjY6OzsnJUbuoztNqtf/wD//wz//8zwMHDhwzZsz8+fM//fRTtYvqqqys\nrMGDB/f2FZWCgoJDhw4dOHDgtdde+9///d+Ghoa9e/eqXVTn+fn5vffee5cvX960aVNISIjBYBg6\ndKjaRXXArVu3UlNTf/WrX61bt+5f/uVf1C6nF2DJzhVyc3Orq6uddxVqampqbGxUt6SuCA4OHjBg\ngPO80YABA3r1cBQHDx6cPXu22lV01YULFwIDA728vIQQGo0mIiLCaDSqXVTnVVdXNzU1bdmyRWnG\nxcXFx8erW1KHuODmy26GGVLP2rNnT0FBgd1uT0lJMZvNQoiSkpKsrKwZM2aoXVpnKMN55JFHbDbb\nkSNHhBBVVVUffvhhVFSU2qV1hjIcZfvLL7/s1Uv8yljGjx//1VdfnT59WghRVVV15MiRiIgItUvr\nDGU4Hh4ezz77bGFhoRDi1KlTFy5cuJt78EtCuflyWloaaXT3mCH1rK1bt8bGxq5evTo2NjYmJsbH\nx6euri4pKamXLg0pw4mIiHj77bf/9V//df369Tab7amnnnrqKXUeFdhFzuHU1NRUVlb2lsu32uT8\nk7Z27drExEStVmuz2WJiYuLi4tQurTOcv5o1a9YkJibW19c3NTXt2LFj0KBBapd2t1x282V3ws1V\nXaehocFisfSK78Peperqai8vrx59+gY6oampqbKyUqfT9e/fX+1auofFYmn/qc1wDwQSAEAKnEMC\nAEiBQAIASIFAAgBIgUACAEiBQAIASIFAAgBIgUACetyZM2dadzY0NHzxxReuLwaQFoEE9KysrKzt\n27e37vf09Fy7dm1ubq7rSwLkRCABPWvbtm13uiXoL37xi9YP7Qb6LAIJ6EGffvrpd999d6eH/sXE\nxJhMpi+//NLFVQFy4uaqQA/avXv33Llz+/Xrl5WV5XyMgmLq1Klr1qyZMWPG/v37x40bp1aFgDwI\nJKB7HD9+vKKiwtkcO3bsmDFjiouLf/aznwkhRo8e7XwcbVFR0fbt23/1q18JIUJCQvbv369KwYBs\nCCSgq2pra5cvXx4cHDx8+PCdO3fOnTt31KhR99xzj9Vq/fbbb4OCgoQQgYGBgYGBQoiqqqo33nhj\n1apVyqN97rvvPpPJZLPZBg8erPIwALURSEBXvfzyy/Pnz4+JiRFCaLXaqqqqp59+Wgjx17/+VQih\nBJKioaFh5cqV48aNW7FihdKjvHr9+nUCCSCQgC4pLi4uLS1V0kgI0dDQ0NDQ4NwWQjR/vvurr75a\nW1v7u9/9ztmjPAl+4MCBrqsYkBVX2QFdUlRU9JOf/MTZPHHixLRp05TtoUOHCiGUR9cLIXbv3v3x\nxx//53/+Z/PHnl65csXb29vPz8+FJQOSIpCALhk4cKCn5+2Vhi+++EKr1U6cOFFpjhw50tvbWwmk\nzz///I033khKSnI4HFeuXLly5co333wjhDCbzb360elAN2LJDuiS6OjoP//5z5988klNTU1RUdHm\nzZudL3l4eDzyyCMXL14UQuzevbuhoeGll15yvnrffff9+c9/vnjx4oMPPqhC3YB8eIQ50A2qq6t/\n9KMf9e/fv0V/Tk7OCy+8kJ2dPWDAgDbfNXPmzD/96U8jRoxwSZmA1FiyA7rBkCFDWqeREGLy5MnD\nhw8/cOBAm+9677335syZQxoBCgIJ6Flr1659//33W/c3NDQcOnTIef03AJbsAABSYIYEAJACgQQA\nkAKBBACQAoEEAJDC/wNjBNWM1EnbSAAAAABJRU5ErkJggg==\n",
"text/plain": [
"<IPython.core.display.Image object>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"z = linspace(8, 12, 8); %Set input array for function file\n",
"\n",
"fig1=figure(1); clf; grid on; axis square; hold on;\n",
"p = scatter(func(a, b, z, q), -z); hold on;\n",
"xlabel('\\sigma(z)'); ylabel('z'); title('z vs. \\sigma(z)');\n",
"legend('\\sigma(z)');"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This plots shows the section of the function where the interpolation will occur. Since a third order interpolation function must be used four points must be selected. The points *closest* to the point $z=10$ must be selected."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"z = z(1,3:6); %redefine z to be vals. closest to z=10\n",
"T = func(a, b, z, q); %find the corresponding outputs"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Therefore the following data will be used:\n",
"\n",
"z = [9.1429, 9.7143, 10.2857, 10.8571]\n",
"\n",
"T = [1.7597, 1.7023, 1.6522, 1.6084]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Calculate $b_0$:"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Therefore b0 is 1.75969\n"
]
}
],
"source": [
"b0 = T(1);\n",
"fprintf('Therefore b0 is %.5f\\n', b0);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Calculate $b_1$:"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"collapsed": false,
"scrolled": true
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Therefore b1 is -0.10041\n"
]
}
],
"source": [
"b1 = (T(2)-T(1))/(z(2)-z(1));\n",
"fprintf('Therefore b1 is %.5f\\n', b1);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Calculate $b_2$:"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"collapsed": false,
"slideshow": {
"slide_type": "slide"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Therefore b2 is 0.01116\n"
]
}
],
"source": [
"%T[z2, z1]\n",
"var1 = (T(3)-T(2))/(z(3)-z(2));\n",
"\n",
"%T[z1, z0]\n",
"var2 = (T(2)-T(1))/(z(2)-z(1));\n",
"\n",
"b2 = (var1-var2)/(z(3)-z(1));\n",
"fprintf('Therefore b2 is %.5f\\n', b2);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Calculate $b_3$:"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Therefore b3 is -0.00094\n"
]
}
],
"source": [
"%T[z3 ,z2, z1]\n",
"var1 = (T(4)-T(3))/(z(4)-z(3));\n",
"var2 = (T(3)-T(2))/(z(3)-z(2));\n",
"\n",
"temp1 = (var1-var2)/(z(4)-z(2));\n",
"\n",
"%T[z2 ,z1, z0]\n",
"var1 = (T(3)-T(2))/(z(3)-z(2));\n",
"var2 = (T(2)-T(1))/(z(2)-z(1));\n",
"\n",
"temp2 = (var1-var2)/(z(3)-z(1));\n",
"\n",
"b3 = (temp1-temp2)/(z(4)-z(1));\n",
"fprintf('Therefore b3 is %.5f\\n', b3);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Next, plug in the solved values into the interpolation function:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\\begin{align*}\n",
"f_n(x)&=b_0+b_1(x-x_0)+...+b_n(x-x_0)(x-x_1)...(x-x_{n-1})\\\\\n",
"f_3(x)&=b_0 + b_1(x-x_0) + b_2(x-x_0)(x-x_1) + b_3(x-x_0)(x-x_1)(x-x_2)\\\\\n",
"f_3(x)&=1.75969 - 0.10041(x-9.1429) \\\\\n",
" & +0.01116(x-9.1429)(x-9.7143)\\\\\n",
" & - 0.00094(x-9.1429)(x-9.7143)(x-10.2857)\n",
"\\end{align*}"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This can be implemented in a MATLAB anonymous function to be plotted and to solve the problem."
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"f_3 = @(x) 1.75969 - 0.10041.*(x-9.1429) + ...\n",
"0.01116.*(x-9.1429).*(x-9.7143)- ...\n",
"0.00094.*(x-9.1429).*(x-9.7143).*(x-10.2857);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Results"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Problem Solution\n",
"To find the stress at depth $z = 10$ the interpolation function is evaluated at this point with the following:"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Therefore the stress at z = 10 is 1.67643\n"
]
}
],
"source": [
"stress = f_3(10);\n",
"fprintf('Therefore the stress at z = 10 is %.5f\\n', stress);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Next, the interpolation function can be plotted on top of the data to visualize the interpolation and to verify that it is accurate. "
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {
"collapsed": false,
"scrolled": false
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAjAAAAGkCAIAAACgjIjwAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAA\nB3RJTUUH4QMVBRAgRsxPiAAAACR0RVh0U29mdHdhcmUATUFUTEFCLCBUaGUgTWF0aFdvcmtzLCBJ\nbmMuPFjdGAAAACJ0RVh0Q3JlYXRpb24gVGltZQAyMC1NYXItMjAxNyAyMjoxNjozMnth2LcAACAA\nSURBVHic7d17XBNX3j/wEwkQsaARUQwqCD/Loru4ShGk0pJUXW2h6gqIK1Rtuy1VFKnXLW6RZ9vS\nakVApLW7W3WteKndl4UHtNZFBUGqwD69UAUKRBORmxAjxmAg+f0xNkYu4ZbMTJLP+6+ZyUnOlxD5\neGZOznA0Gg0BAABg2jCmCwAAACAEgQQAACyBQAIAAFZAIAEAACsgkAAAgBUQSAAAwAoIJAAAYAUE\nEgAAsAICCQAAWAGBBAAArIBAAgAAVkAgAbDRnTt38vPzb9y40f0htVqdn5//448/0l8VgFEhkADY\naN26dc8///y1a9e6PzRs2LC1a9fOmTPn/v379BcGYDwIJADWqaysPHbsmJOT09y5c3tsEB4eLpfL\nv/jiC5oLAzAqBBIA63z99dcajSYsLIzL5fbYIDw8nBCyb98+eusCMC4EEgDJzc3166ampka3zerV\nq/38/I4fP07t7t+/38/P76OPPiKEZGdn+/v7P/XUU/b29oGBgefPnx9iPadPnyaEvPTSS70V5unp\n6eHh8eOPP/78889D7AuAPXr+/xeAReno6FAoFNS2WCxua2sjhDx8+FC3jYeHx8GDBw8fPrxs2TJC\nyOeff37lypWEhISampolS5bweLwtW7a0tLTs3bt34cKF1dXVLi4uffb7888/X716taOjg9rdu3fv\nzp07582bV1paSghxc3P75ZdfeivMzc2turr6+vXrU6dONdT7AMAwDQD8Ki8vz8bGhhCybt26Lg9V\nVlYSQng8XltbW11dHYfD4fP5Dx8+PHnyJCHExcVFLBZrNJp//etff/vb337++ec++4qNjSWEvPfe\ne0uXLiWEzJo163e/+11+fr5EIqH+bd69e1dPYStXriSEJCUlGe6nB2AYRkgAj/z0009Llix5+PDh\nkiVLUlJSujw6ZcqUmTNnlpWVffPNN83NzRqNZsmSJdbW1gEBAc7Ozrdu3Zo8efLUqVOff/75VatW\neXl56e8rPj4+NTW1sbHRycmJEMLhcIYNG/bDDz8QQr7//ntCiL29vYODg57CJkyYQAj55ZdfDPoe\nADAJ15AACCFEKpUuXLjw7t27AQEBmZmZw4b18E+DOlmXlZX19ddfk19nFowfP/7nn3/eu3fvSy+9\ndOPGjYyMDH9//2+++UZPXxKJ5IMPPti+fTuVRhTtHO4RI0YQQpRKpUaj0VOY9kQfgNlAIAGQu3fv\nvvjii1Kp9Omnn87KyuLxeD02oxIoOzs7Ly/P0dHxhRdeIITk5eXt2rXr6aefzs7Obm5ufvbZZ9Vq\n9ZkzZ5RK5alTp06dOtXc3NzldYqKigghixcvpnapk4Fr1qyhdidPnmxtba1SqRobG/UUdvPmTUKI\nbqQBmDqcsgMgmzdvphY+cHV13bZtG3UwLCxs/vz5us3c3Nz8/Py+++47QkhkZCQ1J7uzszMpKemL\nL77Ys2ePnZ1dS0sLIWTmzJmtra1LliwhhOTl5QmFQt3XcXR01N319PT09PSMjo6mdq2srNzd3Ssq\nKqRS6f79+3srjLrU5O3tbfB3A4AxTF/EAmBeaGho938au3bt6t5yz5491KPffvut9uDHH388evRo\n6ridnd3bb7+tVqvr6uqoI3l5ed1fhxAiEAhycnIIIa+//nqXRyMiIgghhw8f7q2wzs7OsWPHEkJq\na2sN+k4AMImj0WgMn3IAFkalUkkkEpVKNWnSpOHDh/fnKQ8fPrx3716X0RLl0qVLgYGBCxcuzM3N\n7fG5586dmzdvnp+fX3Fx8ZDqBmATnLIDMABra2t3d/cBPcXGxqbHNCKEzJkzZ9asWWfPnq2vr3d2\ndu7e4NChQ4SQ9evXD6JUANbCpIahUigUVVVV2m8vAhjEjh075s2b93//93/dH1Kr1Q8ePFi1atXy\n5cvpLwzAeHDKbkgyMjIOHjzo4uIikUji4uJWrFjBdEUAAKYKp+wGr7Cw8PDhw7m5uWPGjKmpqQkL\nC/P39/fw8GC6LgAAk4RTdoNXWlo6e/bsMWPGEELc3d29vb0LCgqYLgoAwFQhkAZv1KhR2qm9Go2m\nvr5euwsAAAOFU3YDIJfLtQtfCgSCBQsWJCcnJycnC4XCrKys+/fvt7e3d39WVFTUlStXCCGRkZFR\nUVG0VkwIIUQul2tXRWOEhRdg4T8+CiCEcDgcV1dXBgswFZjUMAB5eXlJSUnUdmxsbHBw8LVr1/bv\n33/nzp158+bV1NSMHj26+0xcT0/PiooK2ot9TCwWu7m5oQDL7B0FoAATghHSAIhEIpFIpN2VyWRq\ntVq7+nJERAQjAyAAAPOAa0iDZ2Vl9eqrr167do0QcunSpaqqKt24AgCAAcEIafDs7e03b94cGxvb\n3t6uVqszMjL6uWYMAAB0h0AaktDQ0NDQ0NbWVj6fz3QtAACmDafsDABpBAAwdAgkAABgBQQSAACw\nAgIJAABYAYEEAACsgEACAABWQCABAAArIJAAAIAV8MVYADAi7Wr3ZmnWrFmHDx9mugrzgUACACO6\ncuUKs6vdG5WnpyfTJZgVnLIDAABWQCABAAArIJAAAIAVEEgAAI9JpVKlUql7pLOzs7a2lql6LAom\nNQAA212obr34i4zadh3NW+U73kgd3b9/f9myZRcuXNA9aGVl9frrrx84cGDSpEkDejVt2a0y2e+b\nbI1XttnACAkAWO3g1duJ3zweoCSerT149baR+vrwww9feeUVa2vrLsfXr1+/ffv2Ab0UnWWbDYyQ\nAIDVEs/WHojwCvJ4dNex5//fqNXHrhljtNHR0bF3797GxsaMjIz09HTt8ZSUlJdffjkmJqapqcnJ\nyWkQZYvF4sXPuBmpbHOCERIAsFfiN7VBHnxtGhFCqN3Vx64N/cXVavXPP/9cUFBQUFBQU1MjFotH\njx5tY2MTFhb25Zdffvnll2vXrn348KG/vz+Hw5k0adKlS5fYULYZQyABAHsdLLm90te5y8GVvs7i\nlgdDfOXLly+//PLLu3fvfumllw4ePHj+/PnKysqJEycSQpycnKZNmzZ27NiPP/742LFjDg4OhJCJ\nEyf+8ssvAy1b1SipT99gwLLNG07ZAQB7JcyffOhqve5Qg/w6/hjKy965c+fdd989ffo0l8vlcrl/\n+tOfnn/++U8//VR7Rk6j0axcuXLdunXPPPMMdcTZ2fnu3bsDKnuW8po0IXRC4kmlgco2exghAQB7\nBXnwL1S3Xqhu1R65UN0qblUm/GHyUF726NGjq1at4nK5hBCFQkHlkL29fUNDA9Vg165dHA4nLi5O\n+5Tm5mZq/NTPsp/6dl9DetyExJN20wIMVbbZwwgJANjLbTQvYf7k1ceurXrm0XSAC9WtByK8hviy\nN2/e/M1vfkMIaW1tvXPnztSpUwkhHh4eEomEEFJcXJyRkVFSUsLhcHSf0s+V61SNEqt9G5bZKUXO\ne1ZJxxNpbatM9v2dlqGXbfYwQgIAVlvlO173T/lK3/FDP/EVHR1dUFBw7ty5jz/++MCBA9TBZ555\npq2tTalUJiYmNjU1eXt7CwQCgUCQnJxMCBGLxV5efSeKoryodo2f3bSAOXtzDF622cMIqQ8KheLW\nrVsuLi52dnY9NlAqlRKJZNy4cdSVTwAwuC4z1obO3d09ISFBqVTOnTtXe5DL5UZGRh4+fPj06dNd\n2n/77bfPPffcuHHj9L/snRO75edPaE/TacsWi8Vubpjw3TeMkPTJyMgICgrasmXLc889d+TIke4N\nMjMzAwMDt27dKhKJUlNT6a8QAAZn2LBh3f+XuWXLlr///e8ajabL8bS0tMTERD2vpmqUSBKWKsqL\nJn/yHZVGMAgYIfWqsLDw8OHDubm5Y8aMqampCQsL8/f39/Dw0DaoqqratWvXqVOnXF1d6+rqXn75\n5cDAwJkzZzJYMwAMhUAg2L17d1tbm729vfZgZ2fn+vXrdf/td6EoL5ImhDqGb3QM30hLmWYLI6Re\nlZaWzp49e8yYMYQQd3d3b2/vgoIC3QbXr18PCAhwdXUlhAgEgkmTJt24cYOZWgHAQAIDA3XTiBBi\nZWU1b9683trfObGbmk2HNBo6jJB6NWrUqLq6Ompbo9HU19drdykhISEhISHUtlgsrq6unj59Ot1V\nAgCjqNN0TFdhJhBIj8nlcmrSJyFEIBAsWLAgOTk5OTlZKBRmZWXdv3+/vb29xyeWlpbGxcVFR0e7\nu7v32ICaLRoZGRkVFWWk4vWQSqX0d4oCWNI74wV4e3ub8X2+vb29O1fvEovF+pvJ5XJayjF5CKTH\nSkpKkpKSqO3Y2Njg4OCjR4/u378/OTl53rx5IpFo9OjRXZ6iUqk+/vjj3Nzcv/71r/Pnz+/tlSsq\nKoxYdz+4ubmhAIvtndkCvvzyS7FYzOw7YPAC5OdP3Dmxe1zMnn7OX+gzsYCCQHpMJBKJRCLtrkwm\nU6vVKSkp1G5ERET38c26deu4XG5OTg7mfANYCEnCUkIITtMZAyY19MrKyurVV1+9du0aIeTSpUtV\nVVVUXGVmZpaWlhJCzp07J5VKU1NTkUYAloCa2203LWBi4ldM12KeMELqlb29/ebNm2NjY9vb29Vq\ndUZGxvDhwwkhaWlpy5Yt8/HxKSws7DKRITk5Wc+JOwAwXYryoob0uP6fpoNBQCDpExoaGhoa2tra\nyuc//pZ4cXExtZGQkJCQkMBQaQBAH+0SDNZj+7u+KgwCAqlvumkEABZF1Sip37fB2mkiLhrRAIEE\nANAzagkG7dp0YGwIJACAHnRZKRVogEACAHgCdZqOYG437TDtGwDgMe0NjTC3m34YIQEAPILTdMxC\nIAEA4DQdK+CUHQBYOpymYwmMkADAouE0HXsgkADAQuE0HdvglB0AWCJVo0SaEIrTdKyCERIAWJyB\n3tAI6IFAAgDLghsasRYCCQAsBXXRyG5agGP4RqZrgR4gkADAIuCGRuyHQAIA84cbGpkEBBIAmDNV\no0T1yVrVxCm4aMR+mPYNAGaLWoLBav5rzjEpTNcCfcMICQDMk/Y0XeMIAdO1QL8gkADADD0xt1ss\nZrYY6CcEEgCYFcztNl0IJAAwH4ryImlCKFZKNVEIJAAwE1i329QhkADAHGBBIDOAad99UCgUVVVV\nCoVCf7P6+vrm5mZ6SgIAXapGiSRhKdbtNgMIJH0yMjKCgoK2bNny3HPPHTlypLdm9+7dW758eXZ2\nNp21AQD59ZtGjuEbMYXBDOCUXa8KCwsPHz6cm5s7ZsyYmpqasLAwf39/Dw+P7i0TExOfeuop+isE\nsHC4aGRmEEi9Ki0tnT179pgxYwgh7u7u3t7eBQUF3QMpOzubz+f7+PgwUSOA5cJFI/ODQOrVqFGj\n6urqqG2NRlNfX6/d1aqrqzt06NAXX3zx4Ycf6nkpT09PQkhkZGRUVJSRqtVDKpXS3ykKYEnvZlmA\npuV2x/H3hnnMtJr/mrgfX3pl/B2Qy+XMFmAqEEiPyeVyiURCbQsEggULFiQnJycnJwuFwqysrPv3\n77e3t+u2V6vV27Zti4+P5/F4+l+5oqLCWEX3j5ubGwqw2N7NrABFeZH0gwF/04jZd6A/qQkEgaSr\npKQkKSmJ2o6NjQ0ODj569Oj+/fuTk5PnzZsnEolGjx6t2/7QoUN8Pn/EiBGVlZWtra12dnZ1dXUC\nAVbNAjAWXDQybwikx0QikUgk0u7KZDK1Wp2S8miR4IiIiC4n3GQyWW1t7aZNmwgh9fX1NjY2KpUq\nPj6ezpoBLAcuGpk9BFKvrKysXn311YMHD3p5eV26dKmqqoqKq8zMTE9PTx8fn7i4uLi4OKrxjh07\nXF1dV69ezWjJAOYJy9NZCARSr+zt7Tdv3hwbG9ve3q5WqzMyMoYPH04ISUtLW7ZsGabVAdADtx63\nHAgkfUJDQ0NDQ1tbW/l8vvZgcXFx95Y7duygrywAi4Fbj1sUBFLfdNMIAOhBnaYjuGhkSbB0EACw\njqpRIk0IxfJ0lgYjJABgF1w0slgIJABgEVw0smQIJABgBVw0AlxDAgDm4aIREIyQAIBxuGgEFAQS\nADAJF41AC4EEAMzARSPoAteQAIABuGgE3WGEBAB0w0Uj6BECCQBohYtG0BsEEgDQBBeNQD9cQwIA\nOqgaJbVr/HDRCPTACAkAjI66aIRbj4N+CCQAMK7Os/9s+O9ZnKaDPiGQAMCIJAlL1UqlB9II+gGB\nBABGQU1hsJsW0DlrKdO1gGnApAYAMDxFeVHtGj/H8I2O4RuZrgVMBgIJAAwMUxhgcHDKDgAMifre\nK6YwwCAgkADAYCQJSwm+9wqDhUACAAPQTmHARSMYNFxD6oNCoaiqqlIoFHra3L59W38DAPOGKQxg\nEBgh6ZORkXHw4EEXFxeJRBIXF7dixYouDW7cuBEdHd3Z2Xn//n2hUPjee+8xUicAgzCFAQwFgdSr\nwsLCw4cP5+bmjhkzpqamJiwszN/f38PDQ9tAo9H8+c9/XrduXUhIiFKpDA4OLikpeeaZZxisGYBm\nmMIABoRA6lVpaens2bPHjBlDCHF3d/f29i4oKNANpLKyMg6HExISotFoeDze2bNnORwOc/UC0A1T\nGMCwEEi9GjVqVF1dHbWt0Wjq6+u1u5TKykpPT88dO3bk5uYOGzZs+fLl69ev7/GlPD09CSGRkZFR\nUVHGLrs7qVRKf6cogCW9G6kATcvtjuPvDfOYaTX/NbFYTH8BA8J4AXK5nNkCTAUC6TG5XC6RSKht\ngUCwYMGC5OTk5ORkoVCYlZV1//799vZ23fYymSwvLy8mJqagoKCmpub11193c3NbtGhR91euqKig\n4wfonZubGwqw2N4NXoCqUVK7aanz2hQHYTgjBQwCswX0mdlAwSy7x0pKSjb8qrCwcOzYsUePHr15\n82ZycvLkyZNFIpGjo6Nuex6PN3LkyDfffNPW1tbLyyssLOz8+fNMFQ9AD0V5kTQhdELiyf6nEUA/\nYYT0mEgkEolE2l2ZTKZWq1NSUqjdiIiILifc3NzcbGxstNeNbGxsOjs7aasWgH64+zgYFUZIvbKy\nsnr11VevXbtGCLl06VJVVRUVV5mZmaWlpYSQ2bNnKxSK06dPE0JaWlpOnToVFBTEaMkARiRJWKoo\nL5r8yXdIIzASjJB6ZW9vv3nz5tjY2Pb2drVanZGRMXz4cEJIWlrasmXLfHx8eDzep59++pe//GXn\nzp0KhWLx4sWLFy9mumoAw8MqDEAPBJI+oaGhoaGhra2tfD5fe7C4uFi7PWPGjDNnzshkMnt7eysr\nKyZqBDAuVaOkdo3fgKYwAAwOAqlvumnUo1GjRtFTCQDNsAoD0AmBBAA9wxQGoBkCCQB6gFUYgH4I\nJAB4AqYwAFMw7RsAHqOmMIwMWoY0AvphhAQAj2AKAzALgQQAhGAKA7AAAgkAMIUBWAGBBGDRMIUB\n2AOBBGApLlS3XvxFRm27juat8h2PVRiAVRBIABbh4NXbh67eDvJ4tOxI4tnaB+WXF3z3EaYwAHsg\nkAAsQuLZ2gMRXtpACq7NlJ85PiEtB1MYgD3wPSQA85f4TW2QB1+bRpKEpWNvl/1r8fE38tqYLQxA\nFwIJwPwdLLm90teZEKJqlEgSllo7TZyY+NVKX2dxywOmSwN4DKfsAMxfwvzJh67WP2vfVr9vw8ig\nZdQUBmrYxHRpAI9hhARg/oI8+Iryoto1fo7hG6k0ulDdKm5VJvxhMtOlATyGERKA+Rt7u+zD5s+2\nef1tqnQ8kdYSQi5Utx6I8GK6LoAnIJAAzJz6am5D3qEJiSffuveU9ntIK33H43wdsA0CCcCc1adv\n6JRUTfnkO0JI0FiCEAI2QyABmC1qhTrrt/YxXQhAv2BSA4AZoqZ3200LmJj4FdO1APQXRkgA5kbV\nKJEmhGon1AGYCgQSgFlRlBdJE0KxQh2YIpyy64NCoaiqqlIoFHoaVFZWyuVyOqsC6BFu+QomDSMk\nfTIyMg4ePOji4iKRSOLi4lasWNGlwcmTJ3fu3Oni4lJTUxMeHh4fH89InQCEEPn5E3dO7MZN9sB0\nIZB6VVhYePjw4dzc3DFjxtTU1ISFhfn7+3t4eGgbKBSKd99998iRIzNmzKivr587d+6LL744Y8YM\nBmsGi1WfvkHVJEEagUnDKbtelZaWzp49e8yYMYQQd3d3b2/vgoIC3QYcDofL5Y4bN44Q4uDgYG1t\nbWNjw0ytYNmo6d2YUAemDiOkXo0aNaquro7a1mg09fX12l3K8OHDt2zZsmbNmvnz5xcUFISEhEyd\nOpWJSsFy4QbkYE4QSI/J5XKJREJtCwSCBQsWJCcnJycnC4XCrKys+/fvt7e367bv7Oy8efPm3bt3\n6+rq2tvbpVJpS0uLo6Nj91f29PQkhERGRkZFRdHwg3QhlUrp7xQF0NC7puV2xycxVvNfu+f74j2x\nmP4C+gkFYNJTPyGQHispKUlKSqK2Y2Njg4ODjx49un///uTk5Hnz5olEotGjR+u2Ly0tzcnJOXPm\njL29vUajWbly5fHjx9esWdP9lSsqKuj4AXrn5uaGAsysd0V5kfSD/k7vtvD3n/ECxL3/dwF0IZAe\nE4lEIpFIuyuTydRqdUpKCrUbERHRZXxTVVUlEAjs7e0JIRwOx8fHp7q6ms6CwWJhejeYJUxq6JWV\nldWrr7567do1QsilS5eqqqqouMrMzCwtLSWETJ8+/ccffywsLCSEtLS0nD592sfHh9mawRLcObG7\nIT1u8iffIY3AzGCE1Ct7e/vNmzfHxsa2t7er1eqMjIzhw4cTQtLS0pYtW+bj4/Pb3/42Pj4+NjaW\nx+MpFIpFixZFREQwXTWYOUzvBjOGQNInNDQ0NDS0tbWVz3+8aH9xcbF2OyoqasWKFc3NzXw+39ra\nmokawYJIEpZaO03E9G4wVwikvummUXfDhg0bO3YsbcWAZcL0brAECCQAtlM1SmrX+GEKA5g9TGoA\nYDVFeRHSCCwEAgmAvTC9GywKTtkBsNSdE7vl509gQh1YDgQSABthejdYIJyyA2AdrN4NlgkjJAAW\nwfRusGQIJAC2wPRusHA4ZQfACoryImlCf1fvBjBLGCEBME9+/sSdE7vHxexBGoElQyABMAzTuwEo\nCCQAJmF6N4AWriEBMAbTuwF0IZAAmEHdS8I5JoXpQgDYAqfsAOhGfdloZNAyB2E407UAsAgCCYBW\nVBo5hm/EhDqALnDKDoA+1FdfkUYAPUIgAdAEX30F0A+BBEAH6s5G+OorgB64hgRgdNRCDBMST1qP\nnch0LQDshUACMC711dy7P/0HX30F6BMCCcCI6tM3dEqq3D/KYboQABOAa0gD09bWJpFIdI8olcqq\nqiq5XM5UScBakoSlqiaJ9Vv7mC4EwDQgkAYmPT19377Hf1+ys7MDAwM3bdokFAp37drFYGHANtRC\nDFgWCKD/zCSQ0tPTc3NzdY80NDS88MILBuwiNTV1+fLlBw4c0B6RyWTbt29PT0//+uuvz5w5c/Lk\nyaKiIgP2CCZK1SiRJCy1mxaAZYEABsRMAqm4uDguLm7Pnj0ajYY60tnZKZVKDdhFQEDA2rVr//jH\nP2qPXLlyxdnZ2c/PjxDi5OQkEony8/MN2COYIu2yQLgHOcBAmUkgEULeeOON06dPr127VqFQGOP1\nfX1958yZ4+rqqj3S0NAwfvx47a6zs3NjY6MxugZToV2IAYvUAQyC+cyy8/DwOHny5Pr168PDwz/9\n9NNhw4aUtXK5XDt5QSAQ8Pn87m06OjqsrKy0u1wuV6VS9fhqnp6ehJDIyMioqKihVDU4hh0pooDe\nqKvLOj6J4b6V3jhCQMRimnvXAwUwXgAmPfWT+QQSIcTBweGf//zn+++/HxoaunXr1qG8VElJSVJS\nErUdGxsbHBzcvY2tra1SqdTuPnjwwNbWtsdXq6ioGEoxQ+fm5oYCjPr6ivKihpMf9bYskNn/+ChA\nP7HOf1BAD7MKJEKIlZXVu+++O2XKlO3btw/ldUQikUgk0t/GxcVF93NWW1s7ZcqUoXQKJgrLAgEY\nhJlcQ4qOjp45c6Z2d/ny5Z9//vkf/vAHo3bq6+vb3t5+7NgxQkh5efnFixeFQqFRewQWkp8/0ZAe\nN/mT75BGAENkJiOkOXPmdDni5+dHzX8zHjs7u507d77zzjtpaWn37t2LiYmZPn26UXsEtrlzYrei\nvAjLAgEYhJkEEm2io6N1d0Ui0eXLl5ubm/l8PpeLN9Oy1KdvUDVJ8NVXAEPB39Ch4nA4Tk5OTFcB\ndJMkLCWEII0ADMhMriEB0AnLAgEYA0ZIAAOgXYgBX30FMDiMkAD6C2kEYFQYIQH0C7UsUG9ffQWA\nocMICaBvivIipBGAsSGQAPpALcSANAIwNpyyA9BHfv7EnRO7JySetB47kelaAMwcAgmgV1QaYSEG\nAHogkAB6Ri3EgDQCoA0CCaAHWIgBgH6Y1ADQFRZiAGAEAgngCVQaOcekMF0IgMVBIAE8hjQCYBCu\nIQE8IklYimWBABiEERIAUTVKkEYAjEMggaXDkqkALIFTdmDRqDRyDN+IZYEAGIcRElguagFvpBEA\nSyCQwELhdhIAbINAAkuENAJgIQQSWBxFeZE0IRRpBMA2CCSwLNTNjcbF7EEaAbANAgksCNIIgM0Q\nSE9oa2uTSCS6R5RKZVVVlVwu7+0pfTYAlkAaAbAcAukJ6enp+/bt0+5mZ2cHBgZu2rRJKBTu2rWr\ne/vMzMzAwMCtW7eKRKLU1FQaK4WBodJo8iffIY0AWAtfjH0kNTW1uLi4rKxsyZIl1BGZTLZ9+/bP\nPvvMz8+vqakpODj42WefDQh4/Oesqqpq165dp06dcnV1raure/nllwMDA2fOnMnQTwC9Ul/Nbcg7\nhFvtAbAcRkiPBAQErF279o9//KP2yJUrV5ydnf38/AghTk5OIpEoPz9f9ynXr18PCAhwdXUlhAgE\ngkmTJt24cYPmsqFP8vMnOktykEYA7IcR0iO+vr6EkJ9++kksFlNHGhoaxo8fVnCxqAAAGOlJREFU\nr23g7OzcJW9CQkJCQkKobbFYXF1dPX369B5f3NPTkxASGRkZFRVlhNr7IJVK6e+UJQWor+Z2luQ0\nhPzF+tdfK/0s+f1HARRcY+4nCw0kuVyunbwgEAj4fH73Nh0dHVZWVtpdLperUql6fLXS0tK4uLjo\n6Gh3d/ceG1RUVAy55CFxc3OzwALq0zeomiTuH+VYi8XMvgOW+f6jAC0xc/8fMi0WGkglJSVJSUnU\ndmxsbHBwcPc2tra2SqVSu/vgwQNbW9subVQq1ccff5ybm/vXv/51/vz5xisYBopKI9yGHMCEWGgg\niUQikUikv42Li4vu/2tqa2unTJnSpc26deu4XG5OTo6Dg4PBi4RBQxoBmCJMauiVr69ve3v7sWPH\nCCHl5eUXL14UCoXUQ5mZmaWlpefOnZNKpampqUgjVpEkLCWEII0ATI6FjpD6w87ObufOne+8805a\nWtq9e/diYmK0cxbS0tKWLVsml8u7TGRITk7GiTtmSRKWWjtNdI5JYboQABgwBNIToqOjdXdFItHl\ny5ebm5v5fD6X+/i9Ki4upjYSEhJorQ/0QhoBmDQEUh84HI6TkxPTVUDfJAlL7aYFOIZvZLoQABgk\nBBKYPOo25CODljkIw5muBQAGD5MawLQhjQDMBkZIYMKoNHIM34glUwHMAAIJTMaF6taLv8iobdfR\nvBWuHUgjAHOCU3ZgGg5evZ34Ta129x+5V2vX+CGNAMwJRkhgGhLP1h6I8Ary4BNCVI2S5X+P/qvX\n344jjQDMCEZIYAISv6kN8uBTaaQoL6pd4zch8aTdtIDVx64xXRoAGAwCCUzAwZLbK32dya83fqXS\naKWvs7jlAdOlAYDB4JQdmICE+ZMPXa2fpbzWkB43LmYPdd2IGjYxXRoAGAxGSGACgjz41NhIm0YX\nqlvFrcqEP0xmujQAMBiMkMAEjL1d9mHzZ9vGvDFVOp5IawkhF6pbD0R4MV0XABgSAgnYTjs2eovn\npf0e0krf8ThfB2BmEEjAarpn6oIIQQgBmDFcQwL26nLdCADMGwIJWAppBGBpEEjARkgjAAuEQALW\nQRoBWCYEErAL0gjAYiGQgEWQRgCWDIEEbIE0ArBwCCRgBaQRACCQgHlIIwAgCKQu2traJBKJ7hGl\nUllVVSWXy/U/sb6+vrm52ZilmS2kEQBQEEhPSE9P37dvn3Y3Ozs7MDBw06ZNQqFw165dvT3r3r17\ny5cvz87OpqVGs4I0AgAtrGX3SGpqanFxcVlZ2ZIlS6gjMpls+/btn332mZ+fX1NTU3Bw8LPPPhsQ\n0MPfzcTExKeeeorees0B0ggAdCGQHgkICPDx8cnJydFoNNSRK1euODs7+/n5EUKcnJxEIlF+fn73\nQMrOzubz+T4+PnRXbOKQRgDQBU7ZPeLr6ztnzhxXV1ftkYaGhvHjx2t3nZ2dGxsbuzyrrq7u0KFD\nGzdupKlKc4E0AoDuLHSEJJfLtZMXBAIBn9/DTQ06OjqsrKy0u1wuV6VS6TZQq9Xbtm2Lj4/n8Xj6\nu/P09CSEREZGRkVFDbX0gZNKpfR3qqcAdXVZ57H3rSLiG0cIiFhMfwE0Y9v7jwLo1+esKKBYaCCV\nlJQkJSVR27GxscHBwd3b2NraKpVK7e6DBw9sbW11Gxw6dIjP548YMaKysrK1tdXOzq6urk4gEHR/\nqYqKCoOWP2Bubm4sKUBRXtRw8iOXDXtpHhsx+w6w5/1HAYwQ0/IfLzNgoYEkEolEIpH+Ni4uLrof\no9ra2ilTpug2kMlktbW1mzZtIoTU19fb2NioVKr4+Hgj1GsmcKYOAPSw0EDqD19f3/b29mPHjkVE\nRJSXl1+8ePGNN96gHsrMzPT09IyLi4uLi6OO7Nixw9XVdfXq1czVy3ZIIwDQD5MaemVnZ7dz586U\nlJSAgICIiIiYmJjp06dTD6WlpeXn5zNbnmlBGgFAnzBCekJ0dLTurkgkunz5cnNzM5/P53Ifv1fF\nxcVdnrhjxw4ayjNR6uqyhpMfIY0AQD8EUh84HI6TkxPTVZgwRXlR57H36Z/FAAAmB6fswIioM3VW\nEfFIIwDoE0ZIYCza60aNI3qYCg8A0AVGSGAUmMUAAAOFQALDQxoBwCAgkMDAkEYAMDgIJDAkpBEA\nDBoCCQwGaQQAQ4FAAsNAGgHAEGHaNxiAqlEiTQidkHgSaQQAg4YREgyVqlFSu8YPaQQAQ4RAgiFB\nGgGAoSCQYPBUjZL6fRuQRgBgEAgkGCQqjRzDNyKNAMAgEEgwGEgjADA4BBIMRv2+DSODliGNAMCA\nEEgwYJKEpSODljkIw5kuBADMCgIJBgZpBABGgkCCAZAkLLWbFoA0AgBjQCBBf0kSllo7TXQM38h0\nIQBgnhBI0C9UGjnHpDBdCACYLQQS9K0+fQMhBGkEAEaFQII+1KdvUDVJJiZ+xXQhAGDmEEhPaGtr\nk0gkukeUSmVVVZVcLtfzrNu3bysUCiOXxgykEQDQBoH0hPT09H379ml3s7OzAwMDN23aJBQKd+3a\n1b39jRs3Fi5cuHLlynnz5m3fvp3GSukgP38CaQQAtMH9kB5JTU0tLi4uKytbsmQJdUQmk23fvv2z\nzz7z8/NramoKDg5+9tlnAwIer02g0Wj+/Oc/r1u3LiQkRKlUBgcHl5SUPPPMMwz9BAYmP3/i7oXj\nSCMAoA0C6ZGAgAAfH5+cnByNRkMduXLlirOzs5+fHyHEyclJJBLl5+frBlJZWRmHwwkJCdFoNDwe\n7+zZsxwOh5nqDQ1pBAD0wym7R3x9fefMmePq6qo90tDQMH78eO2us7NzY2Oj7lMqKys9PT137Njh\n5+fn7++/d+9e+so1JkV50Z0Tu5FGAEAzCx0hyeVy7eQFgUDA5/O7t+no6LCystLucrlclUql20Am\nk+Xl5cXExBQUFNTU1Lz++utubm6LFi3q/lKenp6EkMjIyKioKEP+GP0jlUr731hdXdZ57H3r+K/E\nYjEjBRgDswVY+I+PAggh+mdFgZaFBlJJSUlSUhK1HRsbGxwc3L2Nra2tUqnU7j548MDW1la3AY/H\nGzly5JtvvsnhcLy8vMLCws6fP99jIFVUVBi0/AFzc3PrTzNFeVHDyY+m/L2UqQKMh9kCLPzHRwEG\n/O+debPQQBKJRCKRSH8bFxcX3Y9RbW3tlClTdBu4ubnZ2NhorxvZ2Nh0dnYaulL6KMqLGtLjJiSe\nZLoQALBQuIbUK19f3/b29mPHjhFCysvLL168KBQKqYcyMzNLS0tnz56tUChOnz5NCGlpaTl16lRQ\nUBCDBQ8FlUbjYvZYj53IdC0AYKEQSL2ys7PbuXNnSkpKQEBARERETEzM9OnTqYfS0tLy8/N5PN6n\nn36ampoqFAoXLlwoFAoXL17MbM2Do00j3HAPABhkoafsehMdHa27KxKJLl++3NzczOfzudzH71Vx\ncTG1MWPGjDNnzshkMnt7e90ZECZE1SiRJoROSDyJNAIAZiGQ+sDhcJycnPS3GTVqFD3FGJyqUVK7\nxg9pBABsgFN2lgtpBACsgkCyUKpGSf2+DUgjAGAPBJIlotLIMXwj0ggA2AOBZHGQRgDATggki4M0\nAgB2QiBZFknC0pFBy5BGAMBCCCQLIklYajctwEEYznQhAAA9QCBZCknCUmuniY7hG5kuBACgZwgk\ni1CfvoEQ4hyTwnQhAAC9QiCZP/XVXFWTBDfcAwCWw9JBZk5RXtRZkuP+UQ7ThQAA9AGBZM6oZbyt\ntx5nuhAAgL7hlJ3Zwg33AMC0IJDMk6K8SJoQihvuAYAJQSCZIdziCABMEQLJ3GAZbwAwUQgkc4Ol\n6gDARCGQzAqWqgMA04VAMh9Yqg4ATBoCyUxgqToAMHUIJHOApeoAwAwgkEzenRO7sVQdAJgBBJJp\nk58/oSgvQhoBgBlAID2hra1NIpHoHlEqlVVVVXK5vLenKBSKyspKPQ2MR1FedOfEbqQRAJgHBNIT\n0tPT9+3bp93Nzs4ODAzctGmTUCjctWtX9/YnT54MCgraunVrYGDg+++/T2Olj5aqm/zJd3R2CgBg\nPFjt+5HU1NTi4uKysrIlS5ZQR2Qy2fbt2z/77DM/P7+mpqbg4OBnn302IODxV3wUCsW777575MiR\nGTNm1NfXz50798UXX5wxYwYN1VJpNC5mDw19AQDQA4H0SEBAgI+PT05OjkajoY5cuXLF2dnZz8+P\nEOLk5CQSifLz83UDicPhcLnccePGEUIcHBysra1tbGxoKBVL1QGAWUIgPeLr60sI+emnn8RiMXWk\noaFh/Pjx2gbOzs43btzQfcrw4cO3bNmyZs2a+fPnFxQUhISETJ06tccX9/T0JIRERkZGRUUNsU5N\ny+2OT2K4b6U3jhCQX0vVTyqVDrHTIbLwAiz8x0cBhBBGrjGbIgsNJLlcrp28IBAI+Hx+9zYdHR1W\nVlbaXS6Xq1KpdBt0dnbevHnz7t27dXV17e3tUqm0paXF0dGx+0tVVFQYqnLJgY3jNuwd6NjIzc3N\nUAUMjoUXYOE/PgoQ9+//jmChgVRSUpKUlERtx8bGBgcHd29ja2urVCq1uw8ePLC1tdVtUFpampOT\nc+bMGXt7e41Gs3LlyuPHj69Zs8Z4ZWOpOgAwYxYaSCKRSCQS6W/j4uKi+/+a2traKVOm6DaoqqoS\nCAT29vaEEA6H4+PjU11dbYRiH8FSdQBg3jDtu1e+vr7t7e3Hjh0jhJSXl1+8eFEoFFIPZWZmlpaW\nTp8+/ccffywsLCSEtLS0nD592sfHx0jFYKk6ADB7FjpC6g87O7udO3e+8847aWlp9+7di4mJmT59\nOvVQWlrasmXL4uLi4uPjY2NjeTyeQqFYtGhRRESEMSrBUnUAYAkQSE+Ijo7W3RWJRJcvX25ububz\n+Vzu4/equLiY2oiKilqxYgXVwNra2iA1XKhuvfiLjNp2Hc0Lqc3EUnUAYAlwyq4PHA7HyclJN426\nGDZs2NixYw2VRgev3k78pla7+8O//4Wl6gDAQmCExC6JZ2sPRHgFefAJtRyD7N/PTdhT2+fTAABM\nH0ZILJL4TW2QB/9xGqXHTf7kuyAP/upj15guDQDA6BBILHKw5PZKX2dq+0H5ZWqpupW+zuKWB4zW\nBQBABwQSiyTMn3zoaj217Ri+kfoCLDVsYrQuAAA6IJBYJMiDf6G69UJ1q/bIhepWcasy4Q+TGawK\nAIAemNTAIm6jeQnzJ68+dm3VM48Wdb1Q3XogwovZqgAA6IFAYpdVvuPdRvO030Na6Tse5+sAwEIg\nkFhHO9EOAMCi4BoSAACwAgIJAABYAYEEAACsgEACAABWQCABAAArIJAAAIAVEEgAAMAKCCQAAGAF\nBBIAALACAgkAAFgBgQQAAKyAQAIAAFZAIAEAACsgkAAAgBUQSAPT1tYmkUj6cxAAAAYEgTQw6enp\n+/bt689B9jh8+DAKsNjeUQAKMCEIpP5KTU1dvnz5gQMH+jzINl988QUKsNjeUQAKMCG4Y2x/BQQE\n+Pj45OTkaDQa/QcBAGAQEEj95evrSwj56aefxGKx/oMAADAICKSeyeVy7TwFgUDA5/MH/VKzZs3y\n9PQ0UF2DhAKYLcDCf3wUMGvWLAZ7NyEIpJ6VlJQkJSVR27GxscHBwYN+KVzPBADoDwRSz0QikUgk\nYroKAAALgll2AADACggkAABgBQ7mKwMAABtghAQAAKyAQAIAAFZAIBme/rVWb9++rVAo6Oyxx06V\nSmVVVZVcLmeqAEp9fX1zczNTBRjwdzGIAhQKRWVlpaF+BX0W09HRce9JSqXS2J3q75eeD2GfP7hB\nPoSD690Yfw1MGqZ9G156erpMJvvwww+7HL9x40Z0dHRnZ+f9+/eFQuF7771n7B576zQzM3PPnj0T\nJ068efNmVFRUbGwszQVQ7t27t3z58ldeeWX16tU0F2Dw38VACzh58uTOnTtdXFxqamrCw8Pj4+OH\nWECfxZw7d27btm3a3YcPHy5ZsuT99983aqd6+qXtQ6j/BzfUh3CgvRvvr4Fp04DhpKSkREREPP30\n01u3bu3ykFqtnjdvXlZWlkajefDgwQsvvHD16lWj9thbp5WVlb///e/FYrFGo7l165aPj09paSmd\nBWgf3bhxY3Bw8Oeffz7o3gdXgGF/F4Mo4P79+15eXmVlZRqN5vbt29OmTaO2h05/MVpNTU2LFi2q\nqamhs1Pdfun8EPZYgPbI0D+Eg+jdSH8NzABGSIakZ63VsrIyDocTEhKi0Wh4PN7Zs2c5HI5Re+yt\n0//93/8NCAhwdXUlhAgEgkmTJt24cWPmzJm0FUA9lJ2dzefzfXx8BtfvUAow7O9iEAUolUoulztu\n3DhCiIODg7W1tY2NzaAL6H8xlM7Ozujo6Lfffnvy5Mm0ddql3+zsbNo+hD0WQB0xyIdwEL2XlpYa\n46+BGcA1JEPy9fWdM2cO9c+si8rKSk9Pzx07dvj5+fn7++/du9fYPfbWaUhIiPbuTWKxuLq6evr0\n6XQWQAipq6s7dOjQxo0bB93vUAow7O9iEAUMHz58y5Yta9asycjIeO2110JCQqZOnTqUGvpZDOXE\niRP29vbPPfecQXrsZ6dd+qXzQ9hjAcRwH8JB9G6kvwZmAIFEE5lMlpeX5+zsXFBQcODAgRMnTmRl\nZTHbaWlp6SuvvBIdHe3u7k5nAWq1etu2bfHx8Twez0j96i+Azt9Fj311dnbevHnz7t27dXV17e3t\nUqm0paXFSAV0oVQq09LS1q1bR093ffZLw4ewxwLo/BB2752RvwYmAYFEEx6PN3LkyDfffNPW1tbL\nyyssLOz8+fNMdapSqZKSkjZs2LB9+/a33nqL5gIOHTrE5/NHjBhRWVnZ2tra1NRUV1dHZwF0/i56\n7Ku0tDQnJycrK+u999776quvOjo6jh8/bqQCusjLy7Ozsxv0yTED9kvbh7DHAuj8EHbvnZG/BiYB\n15Bo4ubmZmNjoz1TbGNj09nZyVSn69at43K5OTk5Dg4O9Bcgk8lqa2s3bdpECKmvr7exsVGpVIad\nZqa/ADp/Fz32VVVVJRAI7O3tCSEcDsfHx6e6utpIBXSRnZ29YMECevrS3y9tH8IeC6DzQ9i9d0b+\nGpgEjJCMLjMzs7S0dPbs2QqF4vTp04SQlpaWU6dOBQUFGbVHQkiPnZ47d04qlaamphrvD4H+AuLi\n4rJ+9eKLL7722msG/0OgvwAafhf6C5g+ffqPP/5YWFhIHTx9+vTQJ3f0pxhCyA8//DCUqzWD67R7\nv3R+CHsswNgfQv290/nXwLQgkIwuLS0tPz+fx+N9+umnqampQqFw4cKFQqFw8eLFRu2RENJjp4WF\nhdQ15N/+6uzZs3QWYNi+BlEADVXpL+C3v/1tfHx8bGzsnDlz5s6dO3v27IiICMMW0GMxd+/ebW5u\n9vDwMF5f3TvtsV86P4Q9FmBs+ntn6t8F+2FxVbrJZDJ7e3srKyuz75T9BdBZVfe+1Gp1c3Mzn8+3\ntramoQBgIcb/XbANAgkAAFgBp+wAAIAVEEgAAMAKCCQAAGAFBBIAALACAgkAAFgBgQQAAKyAQAIw\nurKysu4HOzo6vv/+e/qLAWAtBBKAceXl5WlvtaCLy+XGx8dfuXKF/pIA2AmBBGBc6enpUVFRPT70\npz/9qce7ngNYJgQSgBGdP3/+3r17vd0Nb9GiRWKx+IcffqC5KgB2wu0nAIzoyJEjISEhw4YNy8vL\nS0lJ0X0oMDBw8+bNL7zwwtdff+3t7c1UhQDsgUACMIyLFy82NTVpd6dNm+bl5VVRUfHyyy8TQn7z\nm99obxh6/fr1ffv2bdiwgRAyefLkr7/+mpGCAdgGgQQwVG1tbdHR0W5ubhMnTjx06FBISMiUKVOe\neuopuVze2Njo4uJCCBEIBAKBgBDS0tLywQcfrF+/XiQSEUImTJggFosVCoWdnR3DPwYA0xBIAEP1\nzjvvhIWFLVq0iBDC4/FaWlpCQ0MJIf/9738JIVQgUTo6OtatW+ft7b1mzRrqCPXonTt3EEgACCSA\nIamoqKitraXSiBDS0dHR0dGh3SaE6N6d+n/+53/a2tr+8Y9/aI9Q97G2tbWlr2IAtsIsO4AhuX79\n+u9+9zvtbn5+/vPPP09tOzo6EkKkUim1e+TIkW+//faTTz4ZPny4tv2tW7ccHBzGjh1LY8kALIVA\nAhgSW1tbLvfRmYbvv/+ex+PNmjWL2nV1dXVwcKAC6bvvvvvggw+2bNmi0Whu3bp169athoYGQohU\nKqXz1toAbIZTdgBDMn/+/P/85z/nzp27e/fu9evX9+zZo33Iyspq9uzZNTU1hJAjR450dHRs27ZN\n++iECRP+85//1NTUTJ06lYG6AdgHtzAHMACZTDZixAhra+sux4uLi99+++0LFy7Y2Nj0+Ky5c+f+\n+9//njRpEi1lArAaTtkBGMCoUaO6pxEhxN/ff+LEiVlZWT0+6+jRo8HBwUgjAAoCCcC44uPjv/zy\ny+7HOzo6cnJytPO/AQCn7AAAgBUwQgIAAFZAIAEAACsgkAAAgBUQSAAAwAr/H4uFbSxO0TybAAAA\nAElFTkSuQmCC\n",
"text/plain": [
"<IPython.core.display.Image object>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"fig1 = figure(1); clf; grid on; axis square; hold on;\n",
"p = scatter(func(a, b, z, q), -z);\n",
"xlabel('\\sigma(z)'); ylabel('z'); title('z vs. \\sigma(z)');\n",
"legend('\\sigma(z)');\n",
"\n",
"x = linspace(z(1), z(4), 10); %Set input array for function file\n",
"p1 = plot(f_3(x), -x);"
]
},
{
"cell_type": "markdown",
"metadata": {
"collapsed": true
},
"source": [
"## Discussion\n",
"Most of the MATLAB work in this problem is hard coded. This script would be more useful if the code is taken and generalized. Newton's divided difference method is recursive which mean that the computations from the previous step are used in the next step. This fact could be taken advantage of to create a recursive script that could generate an interpolation function given the order of the function. This nice thing about this script, however, is that any function can be inputed for the `func`. The code would interpolate a third order function as long as the appropriate data is user selected."
]
}
],
"metadata": {
"celltoolbar": "Hide code",
"kernelspec": {
"display_name": "Matlab",
"language": "matlab",
"name": "matlab"
},
"language_info": {
"codemirror_mode": "octave",
"file_extension": ".m",
"help_links": [
{
"text": "MetaKernel Magics",
"url": "https://github.com/calysto/metakernel/blob/master/metakernel/magics/README.md"
}
],
"mimetype": "text/x-matlab",
"name": "matlab",
"version": "0.14.2"
},
"latex_metadata": {
"author": "Julius C. F. Schulz",
"title": "Amazing Article"
}
},
"nbformat": 4,
"nbformat_minor": 2
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment