{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "remove_cell"
    ]
   },
   "source": [
    "<div style=\"width: 100%; overflow: hidden;\">\n",
    "    <div style=\"width: 500px; float: left;\"> \n",
    "        <h1>Ingenieurinformatik – Übung</h1>\n",
    "        Lehrstuhl <b>Computational Civil Engineering</b><br>\n",
    "        Kontakt: <a href = \"mailto: cce-inginf@uni-wuppertal.de\">Email senden</a> | Individuelle Kontakte siehe Webseite des Lehrstuhls<br>\n",
    "        Links: \n",
    "        <a href=\"cce.uni-wuppertal.de/inginf\">Vorlesungsskript</a> | \n",
    "        <a href=\"cce.uni-wuppertal.de/\">Webseite des Lehrstuhls</a>\n",
    "    </div>\n",
    "    <div style=\"float:right;\"> \n",
    "        <img src=\"../../logo_cce_combined.png\" style=\"width:150px;\"/>\n",
    "    </div>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Ableitung und Genauigkeit"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In dieser Übungsaufgabe werden Sie 2 Näherungen der 1. Ableitung einer Funktion implementieren und deren Genauigkeit analysieren."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Aufgabenteil A\n",
    "\n",
    "Schreiben Sie die Funktion `ableitung_o1(funk, x, h)`, die für eine beliebige übergebene Funktion die 1. Ableitung 1. Ordnung der Funktion an allen Stellen des Arrays x berechnet und das Ergebnis als Array wieder zurück gibt. Sie können Funktionen wie Variablen übergeben und dann mit dem zugewiesenen Namen, hier `funk`, auf diese zugreifen.\n",
    "\n",
    "Schreiben Sie zusätzlich eine Funktion `ableitung_o2(funk, x, h)`, die die Ableitung 2. Ordnung berechnet."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "loesung",
     "hide-cell"
    ]
   },
   "source": [
    "### Lösungsvorschlag"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "tags": [
     "loesung",
     "hide-cell"
    ]
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "def ableitung_o1(funk, x, h):\n",
    "    dy = []\n",
    "    y = funk(x)\n",
    "    yh_pos = funk(x + h)\n",
    "    for i in range(len(x)):\n",
    "        dy.append((yh_pos[i]-y[i])/h)\n",
    "    return dy\n",
    "\n",
    "def ableitung_o2(funk, x, h):\n",
    "    dy = []\n",
    "    yh_pos = funk(x + h)\n",
    "    yh_neg = funk(x - h)\n",
    "    for i in range(len(x)):\n",
    "        dy.append((yh_pos[i]-yh_neg[i])/2/h)\n",
    "    return dy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Aufgabenteil B\n",
    "\n",
    "Schreiben Sie eine Funktion, die ein Array entegen nimmt, die Funktion\n",
    "$f(x) = x^3 +2x^2$\n",
    "auf jedes Element des Arrays anwendet und alle berechneten Werte wieder als Array zurück gibt.\n",
    "\n",
    "Berechnen Sie analytisch die Ableitung der Funktion $f(x)$ und implementieren Sie diese nach dem gleichen Prinzip."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "loesung",
     "hide-cell"
    ]
   },
   "source": [
    "### Lösungsvorschlag"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "tags": [
     "loesung",
     "hide-cell"
    ]
   },
   "outputs": [],
   "source": [
    "def f(x):\n",
    "    return x**3 + 2*x**2\n",
    "\n",
    "def df(x):\n",
    "    return 3*x**2 + 4*x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Aufgabenteil C\n",
    "\n",
    "Schreiben Sie ein Programm, welches die Genauigkeit der numerischen Lösungen mit der analytischen vergleicht. Werten Sie dazu alle drei Ableitungsfunktionen an jeweils 10 beliebigen Stellen für alle $h \\in \\{1, 10^{-1}, 10^{-2}, ..., 10^{-15}\\}$ aus. Bestimmen Sie die Differenzen zwischen der analytischen Lösung und den numerischen und bilden Sie immer den Mittelwert über die 10 Stellen. Plotten Sie die Differenzen doppelt logaritmisch gegen $h$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "loesung",
     "hide-cell"
    ]
   },
   "source": [
    "### Lösungsvorschlag"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "tags": [
     "loesung",
     "hide-cell"
    ]
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAVMAAAFQCAYAAAAC353aAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOzdd3xT9f7H8den6aQtZZZV9m7D0DIFBEQFFBApQxS3qNctoCKK4gQnonCd158DrwgoG0W2soessgtllFVaZkt3v78/TvCG2tK0TXqS5vt8PPJoT5Jzvu8mzSdnfr+ilELTNE0rGR+zA2iappUFuphqmqY5gS6mmqZpTqCLqaZpmhPoYqppmuYEuphqmqY5gS6mXkJExojIVw4+d5yITC1GG3VEJEVELEVP6HwislNEuonIWBH5LM9jSkQamZXNGURkhYg8VMBj9Wx/o68Dy7lLRH53fkLvUugLrRWNiNwBPAtYgVQgHvgW+FSZeFKvUurtUmjjCBByeVpEVgBTlVIOFXEX5Imy/brC1W2JiAIaK6XinLxcAQ4A6UqpSGcu+zKl1A/AD3ZtuuRvKev0mqkTichIYBLwHlAdqAY8CnQC/E2MprmII2t+JVzG9UA40EBE2pa0Lc2FlFL65oQbEIaxJhpTyPMCgPeBI8Ap4DMgyPZYNyABGAkkAieA++3mvRXYAlwAjgLj7B7rBiTkaesQcKPt93EYa4mXH7sHOAwkA2MLei7gB/wI/IzxheADjMZYW0oGpgOVbM+tByiMLZ63gBwgHUgBJhfwelwthyNt3Wt7LZOAl+yW2wlYD5yzvY6fAP52jyuML7r9wFlgCiB2jz8A7LY9tgiom2fex23zxgN/2O5Ltf2tQ2zP6wNstWVYA7TM8968AGwHMgDfAl6frzHWGn/J+xpirHGPBzYA54E5+b0Xdv+f/7G9FseANwGL7bH7gFW23//xt9g/nuc1aGT7/Rvb67cAuGh73RvaPfdmYK8t47+BlcBDRXitC3yf3OlmeoCycgN6AdkFfSjsnvcRMBeoBIQC84Dxtse62ZbxOkYRuwW4BFS0e7wFRpFpiVGM+9s95lAxBSJtH5TOGAXyfSAr73OBINsH5Bu7D94zwDogAuOL4XPgR9tjeT/AK+w/NPm8FoXlcKStL205W2EUpea2x9sBbQELUN/2YX3Grm0FzAcqAHWA00Av22P9gTigOcYXw8vAmjzzLra9h0F29zWye861GF+I7W0Z7rW9HwF2781WoPblZeTz+pTD+OK8BYjB+MKw/0JYgVEYrUAwxhfe1ALei9m21y8YY013A/CI7bH7sCuW+fwtVzye9zkY/x9nbK+5L0bxn2Z7rIrtbxhge+xp23v8UBFe63zfJ3e7mR6grNyAYcDJPPetwVgrScPYXBOMb3z7b+2OQLzt92625/raPZ4IdCigzY+AiXbzOlpMX8FWlGzT5YDMPM+di7EG8TFXrrHtBnrYTdewfTh88/kAr+DqxbSwHI60FWH3+AbgjgLaegaYZTetgM5209OB0bbffwUetHvMB+NLra7dvDfkWX7eAvQp8Eae5+wFutq9Nw848D912vb3Btj+l263e3wFMMFuOtL2+lm4ciuhGsYXTZDdc4cCy22/30fJi+lXdo/dAuyx/X4PsNbuMcHYqrpcTB15rfN9n9ztpg9AOU8yUEVEfJVS2QBKqesARCQB45+kKkbB2GwcVwCMfy77o9/Jl+e3uYTtoI6ItAcmYKyJ+GN8wGYUI2tNjH9obDkviUhynud0wFg7Hqps/8U2dYFZIpJrd18OxgfW2Tkcaeuk3e/2r1UT4EOgDcZr7gtsztN+vvPa2p0kIh/YPS5ALYxdEtjnLkBd4F4RedLuPn+Mv/mywpZxLzDd9v+QLSK/2O6bVcAyDmO8Z1XyyeIHnLD7v/NxoP2iKOi1zPseK9vnwT5bYa91Qct2K/oAlPOsxfj2v+0qz0nCWPOMUkpVsN3ClFKO/nP8F2ONsbZSKgxjf+vlT0cqRtEAwHZ6UtUClnMCY9P58nODgMp5nvM7xv64pSJiX7yOAr3t8ldQSgUqpY7l047K576i5ChKW3l9CuzBOCpdHhjD/16rwhzF2AS2bzdIKbXG7jmF/W1HgbfyLKOcUupHR5YhIhHADcAwETkpIieBgcAtImJfLGvb/V4HY809KZ8sGUAVuyzl1f/OdihM3v+t6g7OB/98j8V+Gsdea4+gi6mTKKXOAa8B/xaRgSISIiI+ItIaYz8VSqlcjH18E0UkHEBEaolITwebCQXOKKXSRaQdcKfdY/uAQBG5VUT8MPY9BRSwnJlAXxG5TkT8bbn/UWiUUu9iFPCldh/gz4C3RKSuLX9VESnoC+QU0OAqf09hOYrSVl6hGPvqUkSkGfAvB+e73O6LIhJlazdMRAYVMk/ev/VL4FERaS+GYNt7E+pghrsx3tOmQGvbrQnGAcqhds8bJiKRIlIOY1/7TKVUjv2ClFInML4cPxCR8rb/y4Yi0tXBv2UbECUirUUkEGM3kKMWAC1EpL/trIXHMc50uaw4r7Vb0sXUiWzFZwTwPMa+zlMYO/1fwNh/iu33OGCdiFwAlmB8YBzxGPC6iFzE2N843a7t87bHv8I4KJGK8cHLL+dO4ElgGsaaw0Vb3ox8nvsGxsGLJSJSCePUr7nA77Yc6zAOsuRnEjBQRM6KyMfFyFGUtvIahfFlcxGjsP3k4HwopWYB7wDTbO9RLNC7kNnGAd+KyDkRGayU2gQMByZjHIWOw9j36Kh7gX8rpU7a3zCKz712z/seY5/lSSAQeKqA5d2DsZthly3PTIx90I78LfswCvUSjKPqqxz9I5RSScAg4F2MXWGRwCZs73ExX2u3JFfuDtO8kYiEYBzcaKyUivf2HJrriIgPxpf8XUqp5WbncSa9ZuqlRKSviJQTkWCMU5J2YBxh9socmuuISE8RqSAiAfxv3/U6k2M5nS6m3us24Ljt1hjjlCIzNlPcJYfmOh0xLrxIAvpinBudZm4k59Ob+ZqmaU6g10w1TdOcoEyetF+lShVVr169Is2TmppKcHCwawJ5UAadQ+fwhBxmZdi8eXOSUir/87fNvgTLFbfo6GhVVMuXLy/yPM7mDhmU0jny0jmu5A45zMoAbFIF1B29ma9pmuYEuphqmqY5gS6mmqZpTqCLqaZpmhPoYqppmuYEuphqmqY5gS6mmqZpTqCLqaZpmhPoYqppmuYEupgC6w4mk6s7fNE0rQS8vpj+deQsd3yxjt/is8yOommaB/P6YnpN7Qr0tlbnl/1Z7D150ew4mqZ5KK8vpiLCm/2tBPnBiOlbycrJLXwmTdO0PLy+mAJUDgngvqgAdh6/wORlcWbH0TTNA+liahNdzZcB19Ri8vI4diScNzuOpmkeRhdTO6/2jaJqSAAjpm8lPSun8Bk0TdNsdDG1E1bOj3cGtmR/YgoTF+8zO46maR5EF9M8ujapyp3t6/DFnwfZdOiM2XE0TfMQupjmY8wtzYmoGMTIGdu4lJltdhxN05wtJwsunnTqInUxzUdIgC/vDWzFkTOXmPDrHrPjaJrmbItfhU87QWqS0xapi2kBOjSozAOd6vPd2sOs2u+8F1zTNJPF/gzrpoA1BoKrOG2xuphexXM9m9KwajDPz9zGhXR9uammebzEPTDnSVLCo3k7+y5yc53XJ4cuplcR6Gfhg8GtOXUxg9fn7TI7jqZpJZF+AX4aRpZvOW5LHM6ve5I5l+a8lSRdTAvRunYFHuvWkJmbE1i865TZcTRNKw6lYM5jqDMH+Vf646QFhPPfhzpQKdjfaU3oYuqAJ29oTPMa5Xnxlx2cSc00O46maUW15mPYPY+PZBjbfVvw3+EdqF2pnFOb0MXUAf6+Pnw4uBXn0zIZOyfW7DiaphXFwZWoJeNYKh2ZSh/+O7w99aoEO70ZXUwd1LxGeZ65sQkLtp9g3rbjZsfRNM0R54+RM+N+DlODsfyLqcM70Cg81CVN6WJaBI9c34DWtSswdk4siRfSzY6jadrVZGeQ+eMwMtJSeSp3JJ8/2I3mNcq7rDldTIvA1+LDB4NbkZaZw4u/7EDpoU40zW2lzX8B/5N/MSb3X7z6wABaRIS5tD1dTIuoYdUQXujVjKV7EpmxOcHsOJqm5SN141SCtv4f/8ntwx33PUl03Youb9Pti6mIBIvItyLypYjcZXYegPuuq0f7+pV4fd4uEs5eMjuOpml2Ug5vwXfBCNblRtL0rg/o0KByqbRrSjEVka9FJFFEYvPc30tE9opInIiMtt09AJiplBoO9Cv1sPnw8RHeH9QKpRTPz9zu1KsoNE0rvpRzSaR8ewdnVTmyBnxF56bVS61ts9ZMvwF62d8hIhZgCtAbiASGikgkEAEctT3N+T02KwUnY7FkpxVpttqVyvFyn0jWHEhm6vrDTo+laVrRpGVksffTO6mcc5pD3afQpXVUqbZvSjFVSv0B5O0stB0Qp5Q6qJTKBKYBtwEJGAUVXJH3+F/wWScqJ68v8qx3tK1Nt6ZVGb9wD/FJqU6PpmmaY9KzclgwZSTRGevZ3Wo0HbrdWuoZxKwj0iJSD5ivlLLapgcCvZRSD9mm7wbaAy8Ak4F0YJVS6ocClvcw8DBAtWrVoqdNm+ZYEJVLh3XDOR9Ym93XjCvy33E2PZeXV6dRI9iHMe0D8REp8jIuS0lJISQkpNjzO4vOoXO4ew77DNm5ilUb1/PKpQnsDO1McvRIKMHn8Gq6d+++WSnVJt8HlVKm3IB6QKzd9CDgK7vpu4FPirPs6OhoVSSLXlI54yoplZpctPlsZm9JUHVfmK/+vTyuWPNftnz58hLN7yw6x5V0jiu5Q47LGbKyc9SYr+erM6/UVGfei1YqI8Wl7QKbVAF1x52O5icAte2mI4DSudTIGoOPyoY984s1e79WNeltrc7ExfvYc/KCk8NpmpafnFzFCz9tYEj8SwT7QsX7fwJ/518m6ih3KqYbgcYiUl9E/IE7gLml0nKN1qQFVofYX4o1u4jwZn8r5YP8ePanbWRm5zo5oKZp9nKVYvTP22mzawItfeLxH/QFVG5oaiazTo36EVgLNBWRBBF5UCmVDTwBLAJ2A9OVUjtLKRCJ4V0gfiWknC7WIiqHBDBhQAt2n7jApKV6ZFNNcxWlFFN3ZeKz9XuG+i6HLiOhWekfcMrLrKP5Q5VSNZRSfkqpCKXUf2z3L1RKNVFKNVRKvVWamRLDu4DKhV2zi72MGyOrMbhNBJ+uOMBfR846MZ2maZd98Ps+Tifs4y3/b1ENukH3l8yOBLjXZr6pUkPqQtXmxd7Uv2xsn0hqhAUxcroe2VTTnG3m5gSmLt/C10GTsJQPR2K+Bh+L2bEAXUyvZI2BI2vg/LFiLyI00I/3B7UiPimVd/TIpprmNOsPJvPqL5uYVv5jKquzyODvILh0LhV1hC6m9qwDjJ87Z5VoMR0bGiObfqtHNtU0pziUlMqj32/k30Gf0yxzJ7ubPwO1os2OdQVdTO1Vbgg1WsHOkm3qAzzfyxjZ9LmZ2zjvxEG7NM3bnL+UxQPfbGSE+o6u2Wvg5jc5Hd7F7Fj/oItpXtYYOLYZzsSXaDGBfhYmDmlN4sUMXptbOiclaFpZk5WTy79+2EyP8zO5mwXQ/lHo+ITZsfKli2leUbcbP52wdtoyogJPdG/EL1uO8VvsiRIvT9O8iVKKl2fFUiF+AWMsU6F5X+j5tssuFS0pXUzzqlAHarcv8VH9y564oREtaoUxZlYspy9mOGWZmuYNvvjjIAc3L+bjgM+Q2u1gwJduc+Q+P7qY5scaA6diIbHkR+P9LMbIpikZ2YyZpYc60TRHLNp5khmLlvJN0EQslerC0GngF2R2rKvSxTQ/kf1BfJyyqQ/QuFooz/dsyuJdp5iphzrRtKuKPXaet6Yt54fA9ygXFIQMmwnlKpkdq1C6mOYntBrU7WRs6jtpTfKBTvX1UCeaVoiT59N58ps/+NL3XapaUpC7pkPFembHcogupgWxxkDyfji5wymLuzzUSa5SPDdDD3WiaXmlZmTz8DdreSPzPZpwGJ/B30HNa8yO5TBdTAvSvB/4+ELsz05bZO1K5XilbyRrDybz7dpDTluupnm6nFzF0z9u4e6kiXSWbUjfj6DxTWbHKhJdTAsSXBkadHfqpj7A4Da1uaFZOBN+3UNcYorTlqtpnmzCr7uxxv2bQZaV0HU0XHuP2ZGKTBfTq7HGwPkjkLDJaYsUESbEtKCcv4WR07eSnaP7PtW8248bjnB+zdc84/sLtB4G3UYXPpMb0sX0aprdApYAp27qA4SHBvJm/xZsSzjPv1cccOqyNc2TrI5LYvGcqYz3+w+5DXtA34/c9qT8wuhiejWBYcZ+m52zINe5o0zf2rIGt7WuycdL9xN77LxTl61pniAu8SIfT53BFL9JUC0Kn8HfgsXP7FjFpotpYawDIOUkHFnr9EW/3s9K5RB/nv1pK+lZzi3WmubOzqRm8tLXC5jCBPxCq2AZNhMCQs2OVSK6mBamSS/wK+f0TX2AsHJ+vDuwFfsTU/hwsR7qRPMOGdk5jPx2GePTXqOCv8L3nl8gtLrZsUpMF9PC+AdD096waw7kOL8rva5NqnJX+zp8+edB9p7Ra6da2aaUYuyMTTx+cix1LUn43vUTVG1qdiyn0MXUEdYYuJRsDLjnAmNuaU6dSuX4akcGKRl6qBOt7Pp58xG673qJaJ/9WGK+hLodzY7kNLqYOqLRjRAQ5rSepPIKDvDlg0GtSEpTvLVgl0va0DSzZeXkEr9oCr0tG+HmNyGqv9mRnEoXU0f4BhhDye6eD9mu6UavTb1K9K7vx48bjrJ8T6JL2tA0M/28MZ6hmT9zvvI1SMfHzY7jdLqYOsoaAxnnIW6py5q4vbEfzaqH8vzP2zmbmumydjSttGVk53BgyddESBLle77oseeSXo0upo5q0BWCKrnkqP5lfj7CB4Nbce5SJi/PiXVZO5pW2mZsPMzQzJlcrBiJNL7Z7DguoYupoyx+EHkb7F0ImakuayaqZhjP3NiEBdtPMHfbcZe1o2mlJT0rhz1LvqeBz0lCbhpdJtdKQRfTorHGQNYl2LfIpc08cn0DrqlTgbGzYzl5Pt2lbWmaq01bf4hhWTO4FNYYadbX7Dguo4tpUdS9DkKqu3RTH8DX4sOHg1uTmZ3LCz9v10OdaB4rPSuH2OXTaOZzlKAez4NP2S05ZfcvcwUfi3E6x/7FkH7BpU3VrxLMmFuasXLfaf674YhL29I0V5m69hD3ZM0gPbQuEjXA7DgupYtpUVljICfD2HfqYsM61KVL4yq8OX83h5Jct59W01zhUmY221b8TEufeAK7jwKLr9mRXEoX06KKaAthdVy+qQ9G36fvDmyJr0UYNWMbOXqoE82DfL/mEPdkzyAzuCa0vMPsOC6ni2lRiYD1djiwDC6dcXlzNcKCeP22KDYdPsuXfx50eXua5gwpGdlsXDmPtj778O86Anz9zY7kcrqYFoc1BnKzYffcUmmuf+ta9LZW58Pf97H7hGv31WqaM3y75hD3Zs8kK6gqXDPM7DilQhfT4qjeEio3KpVNfTA299/sb6V8kB8jpm8jM1sPdaK5rwvpWaxd+RtdLLH4dXka/ILMjlQqdDEtDhGIGgCHVsHFU6XSZOWQACYMaMHuExeYtFT3faq5r/9bdYj7cmaSHVARou83O06p0cW0uKwxoHKNfk5LyY2R1RjcJoJPVxxg8+Gzpdaupjnq/KUsVq1ayo2WLfhe9zgEhJgdqdToYlpc4c0gPKrUNvUvG9snkhphQYycvpVLmbrvU829/GfVQe7L+Zkcv1BoN9zsOKVKF9OSsA6Ao+vg3NFSazI00I/3B7XiUPIlJvy6p9Ta1bTCnE3NZMXqVfS2bMTS4REIqmB2pFKli2lJWG1XdOycVarNdmxYmQc71+e7tYf5Y9/pUm1b0wry5Z8HuS/3F/ANhA6PmR2n1OliWhKVGkDNa0t9Ux/guZ5NaRQewvMzt3P+kvPHptK0okhOyWDpmnXcZlmDT9sHIbiy2ZFKnS6mJWUdACe2QvKBUm020M/Ch4NbcTolg1fn6r5PNXN98cdB7lOzEYsfXPek2XFM4RHFVET6i8iXIjJHRNyrZ9mo242fO10zPtTVtIyowJM3NGL21uMs3HGi1NvXNIDTFzNYtHYTgyx/4nPtPWVi2ObicHkxFZGvRSRRRGLz3N9LRPaKSJyIjL7aMpRSs5VSw4H7gCEujFt0YRFQp6PLBtsrzOPdG9GiVhgvzdpB4kXd96lW+j5dcYD7mYfFB+j0tNlxTFMaa6bfAL3s7xARCzAF6A1EAkNFJFJEWojI/Dy3cLtZX7bN516sMZC4C06V/siifhYfJg5pRWpmDi/+vEP3faqVqlMX0vlt/Tbu9F2OtLoDKtQ2O5JppDQ+fCJSD5ivlLLapjsC45RSPW3TLwIopcYXML8AE4DFSqklBTznYeBhgGrVqkVPmzatSBlTUlIICSneCcZ+mee4bs39HK47kEP17yrWMkqaYdGhLH7ck8kDVn+uj/ArdoaS5nAmncP9c3y/K4N2x79juO9CNrabQlq5mqWeoTR17959s1KqTb4PKqVcfgPqAbF20wOBr+ym7wYmX2X+p4DNwGfAo4W1Fx0drYpq+fLlRZ7nCt/2U2pSa6Vyc4u9iJJkyMnJVYM/W6Mix/6qjiSnFns5Jc3hTDrHldwtx7Gzl1TbMT+p9NeqKTXzQVMylDZgkyqg7ph1ACq/EbUKXEVWSn2slIpWSj2qlPrMhbmKL2oAnDkIJ7aZ0ryPj/D+oFaIGH2f5uq+TzUXm7I8jrt9fiUgNw26jDQ7junMKqYJgP3OlQjAs4fibN4XfHxNOef0stqVyvFKn0jWx5/h69XxpuXQyr6jZy6xcNNeHvL/3fjfD29udiTTmVVMNwKNRaS+iPgDdwCl0zmoq5SrBA17GFdD5ZrXRd6gNhHc2DycdxftZf+pi6bl0Mq2KcvjGOazmKCcFOgyyuw4bqE0To36EVgLNBWRBBF5UCmVDTwBLAJ2A9OVUjtdncXlrDFw/igkbDQtgogwfkBLQgJ8eXb6VrJydN+nmnMlXspl/uY4Hg34DRrfDDVbmx3JLbi8mCqlhiqlaiil/JRSEUqp/9juX6iUaqKUaqiUesvVOUpF097GdckmbuoDVA0N4K3+VmKPXWDysjhTs2hlz9wDWdxpWUZw9jm9VmrHI66A8hiB5aHxTbBrNuTmmBqld4sa3H5NLSYvj2Pb0XOmZtHKjvikVDYeu8STgb9CvS5Qp73ZkdyGLqbOZo2BlFNweLXZSRjXL4qqIQGMmL6V9Cxzi7tWNny8dD93+K4kNCsJrn/O7DhuRRdTZ2vcE/yCTd/UBwgL8uO9QS05cDqVd3/ba3YczcPFHjvPgq2HeSpgPkS0g/rXmx3Jrehi6mz+5aDZLcZwJjnmd43XpXFV7ulYl69Xx7PmQJLZcTQPdf5SFo/98Bd3l9tA5ZzTxlqp5He6uPfSxdQVrDGQdhYOrjA7CQCjezejfpVgnpuxnQvp5hd4zbPk5ipGztjKyXOpPBe8gIshDYxjA9oVdDF1hYY3QGCYW2zqA5Tz9+WDwa04cT6NN+aVfmcsmmf77I8DLNmdyLctdxJ4IZ7DdQfrtdJ86GLqCr4B0Kwv7J4PWe7RLd61dSryr24NmbE5gd93njQ7juYh1sQl8f6ivdwd6UuHg59Ag24kVelgdiy3pIupq1gHQOZFiMu3kytTPN2jCc1rlGfMrB0kp2SYHUdzcyfPp/Pkj1toUCWYV+UrROVA30l6rbQAupi6Sv2uUK6y22zqA/j7Gn2fXkjLZsws3fepVrDM7Fwe+2Ez6Vk5TG1/FN8Dv0OPV6BiPbOjuS1dTF3F4guR/WHfb5CZanaavzWrXp4RNzdh0c5TzNpyzOw4mpsa/+tu/jpyjol9I6i++hWIaAvtHjY7llvTxdSVrDGQdQn2/mp2kisM79KANnUr8uqcnRw/l2Z2HM3NzNt2nP9bfYgHOtXn5sMTITMF+k0GH4vZ0dyaLqauVKcjhNYwbXyoglh8hA8GtyJHKZ6bqfs+1f4nLvEiL/y8nTZ1KzKm0SGInWlcfx/ezOxobk8XU1fy8TFGL41bDGnudX183crBvHRrc1bHJfP9usNmx9HcQEpGNo98v5ly/hamxDTGd+FICI+Ezs+aHc0j6GLqatYYyMmEvQvNTvIPd7arQ9cmVRn/624Onk4xO45mIqUUo3/eTnxSKh8PvYZqG96GlJPG5r2vv9nxPIIupq5WKxoq1HGro/qXiQjvDmxJgK+FEdO3ka37PvVa36w5xPztJxjVsynXWfbApq+hw2MQEW12NI+hi6mriRhrpweWQ2qy2Wn+oVr5QN7ob2Xr0XN8tvKA2XE0E2w+fIa3FuzmxubVeLRjTZj7pHEKVPeXzI7mURwqpiIyyJH7tAJYY0DlwO45ZifJV79WNenTsgYfLdlP7LHzZsfRSlFSSgaP/fAXNSsE8cHgVvj8McEYGLLfJ0anPZrDHF0zfdHB+7T8VLNC5cZud1Tf3hu3WakU7M/I6dvIzNFH971Bdk4uT/53C+cuZfHpsGsJO7MD1nwC196ju9crhqsWUxHpLSKfALVE5GO72zdAdqkkLAsub+ofWgUXTpidJl8Vg/15Z2BL9p66yKw43bOUN/hw8T7WHkzmzf5WoqqVMzbvg8PhpjfMjuaRClszPQ5sAtKBzXa3uUBP10YrY6wDAGUMaeKmujcNZ2i7OvwWn8XGQ2fMjqO50OJdp/j3igMMbVebQW1qw+pJcCoW+nwIQRXMjueRrlpMlVLblFLfAg2VUt/a3X5RSp0tpYxlQ9WmUK2FW2/qA7x8a3OqBAkjpm8lJUNvfJRFh5NTGTF9K9Za5Xm1bxSc3gcr3zHOiW52q9nxPFZhm/nTbb9uEZHteW+lkK9ssQ6AhA1w1n1Pkg8O8GV4ywASzqbx1oLdZsfRnCw9K4dHp/6Fjwif3hVNoEVg7hPgHwy93zU7nkcrbDP/advPPkDffG5aUVgHGD93zjI3RyGaVLTwcJcG/LjhCMv3JpodR3MSpRQvz45l94kLfMR0JQUAACAASURBVDSkNbUrlYONX8HR9dBrAoSEmx3RoxW2mX/C9vPyqlRj2++JgN6pVlQV6xkn8bvhCfx5PXtTE5pWC+WFmds5m5ppdhzNCX7aeJSZmxN46oZGdG8WDueOwJJx0OhGaDnE7Hgez9HzTIcDM4HPbXdFAO57JMWdWWPg5HZI2m92kqsK9LPw4ZBWnL2Uydg5sWbH0UroUFIqr8zdSZfGVXj6xiagFMx7xniwz0Td4bMTOHqe6eNAJ+ACgFJqP6C3CYoj6nZA3P5AFEBUzTCe7tGY+dtPMHfbcbPjaCXw8bL9+Ah8MKgVFh+BbdPgwFK4cZxxubNWYo4W0wyl1N/beiLiC+gzu4ujfE2oe52xqe8BPd0/2rUhrWtXYOzsWE5dcI/xrLSiOXA6hdlbjjGsfV3CywdCSiIsehFqd4C2D5kdr8xwtJiuFJExQJCI3ATMAOa5LlYZZx0ASXsh0f1HCvW1+PDh4FZkZOfw/MzteqgTD/Tx0v0E+Fp4pGtD445fnzdGf+j3idFNpOYUjr6So4HTwA7gEWChUkr3glBczW8DsXjEgSiABlVDeLF3c1buO81/NxwxO45WBHGJF5m77Tj3XFeXqqEBxoi5O2dB1xegahOz45UpjhbTa5RSXyqlBimlBiqlvhQRfWpUcYVUNa599pBNfYC7O9Slc6MqvLVgN4eT3WdMK+3qPlqynyA/C49c39DooHzBSOPikU5PFz6zViSOFtMvRaTF5QkRGQq87JpIXsIaA2cPwfG/zE7iEB8fo+9Ti48wcvo2cvRQJ25v78mLLNhxgvuuq0elYH9YPBZST8Ntk8HiZ3a8MsfRYjoQ+FZEmttOk3oMuNl1sbxA8z7g4+cRR/Uvq1khiNf6RbHp8Fm+/POg2XG0Qkxauo9gf1+Gd2kAh1bDX9/BdU9AzdZmRyuTHCqmSqmDwB3AzxiF9WallO74siSCKhonS++cBbme08P97dfUoldUdT78fR97Tl4wO45WgN0nLrBwx0nu71SPioFibN6H1YGuo82OVmYVdm3+Drvr8GcClYB6wHp9bb4TWGPgwjHjcj4PISK8dbuV8kG+PPvTNjKzPeeLwJt8tGQfoQG+PNS5Aaz/HE7vht7v6A6fXci3kMf7lEoKb9W0F/gGGgei6nY0O43DKocEMH5AS4Z/t4lJS/fxXE89DLA7iT12nkU7T/F0j8aEZZ+GFeOhcU9o2tvsaGVaYdfmH7Zdi38xn5u+JKakAkKhSU+jj9Mcz+ru7qbIagyKjuDTFQf464jujdGdfLRkP+UDfXmgc334/WXIyTLWSvUloy7l6AGovzDOM90H7Lf9Hi8if4mIHr6wJKwxxhHWQ3+anaTIXukbSY2wIEZO38alTM/6MiirtiecY8nuUwzv0oCwE6uNrZ4uI6BSfbOjlXmOFtPfgFuUUlWUUpWB3sB0jKP6/3ZVOK/Q+GbwD/GYE/jthQb68f6gVsQnpTLh1z1mx9Ew1korlPPjvg41YeFzRk9lnZ4xO5ZXcLSYtlFKLbo8oZT6HbheKbUOCHBJMm/hF2T0br57HpLreWMvdWxYmQc61ee7tYf5c/9ps+N4tS1HzrJsTyLDuzQgdMsXkLQPer8HfoFmR/MKjhbTMyLygojUtd2eB86KiAVw+eFcEQkWkc0iUjYPiFljIP0cFc9uMztJsTzfqykNqwbz3IztnL/keV8IZcVHS/ZTKdif+6J8YeW70KwPNNGng5cWR4vpnfyvD9M5QB3bfRZgcEEzicjXIpIoIrF57u8lIntFJE5EHDnx7QWM3QplU4PuEFiB8ETP228KRt+nE4e05nRKBuPm7TQ7jlfafPgMK/ed5uHrGxC8fKxxmXKv8WbH8iqFnRoFgFIqCXiygIfjrjLrN8Bk4LvLd9jWZqcANwEJwEYRmYtRmPO++w8ALYFdQNndVvH1h+Z9qbJ9BmSlGZv+HqZlRAWe6N6ISUv3c3NkNXq3qGF2JK8ycfF+Kgf7c1+1A7B8LtwwVvdTWsrkal2qichHSqlnRGQe+fRfqpTqV2gDIvWA+Uopq226IzBOKdXTNv2ibVn5fo2KyFtAMBAJpAG3K6X+sWtBRB4GHgaoVq1a9LRp0wqLdoWUlBRCQkKKNI8zVTyzlVbbXyU26gWSql5nWg4o/muRnat4c106SWm5vNk5iAoBJevezez3xFNy7D2Tw/gN6dzZRHj59AhA2Nj2Y5SPa66/d4fXw6wM3bt336yUapPvg0qpAm9AtO1n1/xuV5vXbhn1gFi76YHAV3bTdwOTHVjOfUAfR9qMjo5WRbV8+fIiz+NU2Vkq483aSv10t7k5VMlei/2nLqjGLy1UD36zQeXm5pqWw5ncPccdn69V0W8sVplLJyj1anml4paakqM0mZUB2KQKqDtX3cxXSm22/VxZgmKeV35nDhfaBZFS6hsnZnA/Fl9OV+1ErX2/Q8ZF44R+D9QoPJTnezblzQW7mbEpgcFta5sdqUxbeyCZtQeTebdHGH5rPoTI/tDwBrNjeSVHB9SLF5GDeW/FbDMBsP+ERaCvpgIgMbwLZKfB3t/MjlIiD3SqT/v6lXh9/i6OnrlkdpwySynFxCX7CA8NICZxstHheM+3zY7ltRw+zxRoa7t1AT4GphazzY1AYxGpLyL+GL1RzS3mssqU82HNILSmR57Ab8/HR3h/UCsARs3YRq7u+9Ql1h5IZkP8GcZbj2HZ/yt0ewHCapkdy2s52gVfst3tmFLqI6DQbQkR+RFYCzQVkQQReVAplQ08ASwCdgPTlVL6fBoA8THGh4pbAmmefb177UrleKVPJOvjz/D16niz45Q5Sik+XLyPeuWF7vEfQJWm0P5fZsfyag6dGiUi19pN+mCsqRa6U08pNbSA+xcCCx1p2+tYB8DaycZYPdfebXaaEhnUJoLfd53k3UV76dqkKo2reeZ+YHf05/4kNh0+y5yolfgcOAz3zjdOsdNM4+hm/gd2t/FANFc5WV8rgZrXGtdT7/ScHvgLIiKMH9CSkABfRkzfRlaO7vvUGS7vK21X/hwtD30DLQZB/S5mx/J6jm7md7e73aSUGq6U2uvqcF5JxLi89OBKSPH8a92rhgbwVn8rO46dZ/Kyq13foTlqxb7TbDlylo/C/otY/OHmN82OpOH4mikicquIPC8ir1y+uTKYV7PGgMqB3XPMTuIUvVvU4PZrajF5eRzbE86ZHcejKaX4aPE+hpbfQc3Tq6D7GAitbnYsDcdPjfoMGIJxSakAg4C6Lszl3cIjjQMKHjTYXmHG9YuiakgAz/60lfSsHLPjeKxlexLZl3CKly3fQngUtHvY7EiajaNrptcppe4BziqlXgM6cuW5opozXd7UP7wGzh8zO41ThAX58d6glhw4ncq7v+k9RMVxeV/piyELCU47Abe+DxaHjiFrpcDRYppm+3lJRGoCWYDuutuVrAMAZQxpUkZ0aVyVezrW5evV8aw5kGR2HI+zJTGHS8f3MCxnNrQaCnXN7cNBu5KjxXS+iFQA3sMYwuQQULSeRLSiqdIYqrcsU5v6AKN7N6N+FaPv04vpuu9TR+XmKmbtz+SdclMR/3Jw0+tmR9LycPRo/htKqXNKqZ8x9pU2U0qNdW00DWsMHNsEZw+ZncRpyvn78sHgVpw4n8br83aZHcdj/L7rJC0uraVtzlbkhrEQEm52JC2PohzNv05E7sQ4EHWbiNzjulgaAFG3Gz/L2NrptXUq8q9uDZmxOYHFu06ZHcftZefk8tni7Yzzn4qq3hLaPmh2JC0fjh7N/x54H+jM/67Rz79PP815KtaFiLZlrpgCPN2jCc1rlOfFX7aTnJJhdhy3Nnl5HL2SvyOcM8itH4KPxexIWj6K0tFJJ6XUY0qpJ223p1wZTLOxxsCpHXC6bB0B9/f1YeKQVlxIy+alWbGX+6zV8th46Azrls1luO+vnKh+I9Rua3YkrQCOFtNYQJ8ZbIbI/oCUybXTZtXLM+LmJvy28ySzt5aNU8Cc6fylLCb89zc+9/8IqVSfAw3vNzuSdhWOFtMqwC4RWSQicy/fXBlMsylfA+p1Nq7VL4Nrb8O7NKBN3Yq8Mmcnx8+lFT6Dl1BK8erMdbyd/jbBAT743DWdbD/zh07RCuZoMR0H9Afe5spOT7TSYB1gjIF+Krbw53oYi4/wweBW5OQqnp+5Xfd9avPj+sPcsv9VGluO4zv4W6jc0OxIWiEcPTVqJca5pX623zdinG+qlYbmtxm9qHt4p9EFqVs5mJdubc6quCSmrj9sdhzT7T91kQsLX+Vmy2aj5/yG3c2OpDnA0aP5w4GZwOe2u2oBZefSHHcXXBkadDOKaRnc1Ae4s10dujapytsLd3PwdIrZcUyTnpXDz99+xKM+s0lrMQyf9o+YHUlzkKOb+Y8DnYALAEqp/YA+a7g0WWPg3BE4ttnsJC4hIrwT05IAXwsjpm8j20v7Pv12xi88kzqJc1XbEnTbRKOfBs0jOFpMM5RSmZcnRMQXB0YU1Zyo2a1g8S+TR/Uvqx4WyBv9rWw9eo7P/yjueI2e64/N2+m/dxRpAVWocN803XO+h3G0mK4UkTFAkIjcBMwA5rkulvYPQRWg0U3GUf3csrvW1q9VTW5tWYOJi/dx+IL3dNV3MvkslebdR3lJJ/je6RBcxexIWhE5WkxHA6eBHcDDwAKl1EsuS6XlzzoALp6AI2vNTuJSb95mpWKwP19sz/CKvk9zcnKJ++p+rBzgXO8p+NdqaXYkrRiuWkxF5DYReVwplauU+hKjk5M2wBgRGVgqCbX/adob/MqV2aP6l1UM9ufdgS05lqJ4b1HZuvIrPxumjqVz2nJ2NH2KGu31x8pTFbZm+jxXjmnvjzGYXjdAjytb2vyDoUlPo4/TnGyz07hU96bh3FDHl/+simd1XNnt+/TAqhm0PziFTaE9sA55zew4WgkUVkz9lVJH7aZXKaXOKKWOAMEuzKUVxBoDl5IhfqXZSVxuSFN/GlQJZtSMbZy/VPb6Pk05up0aS55gr09Dmjz8DeLjcCdumhsq7N2raD+hlHrCbrKq8+NohWp0E/iHlomhoAsTYBEmDmnN6YsZjJ1Ttq7+UqlJZHw3iIsqkKxBUykfWt7sSFoJFVZM19tO2L+CiDwCbHBNJO2q/AKheR/YPQ+yy37Xda1qV+CpHo2Zu+04c8pKZyjZmZz+zxBCMpP549pJtIxsbnYizQkKK6bPAveLyHIR+cB2WwHcBzzj6nBaAawxkH4eDiwzO0mpeKxbQ66pU4Gxs2M9vzMUpbgw61nCz2zii4rPMqDvbWYn0pzkqsVUKZWolLoOeAPj2vxDwOtKqY5KKd1FulkadIOgimX+qP5lvhYfJg5uTXauYtSMbR7dGUrW+i8ov3MqX3Mbg+4ficVHX+FUVjja0ckypdQntpt3rA65M4sfNO8HexZC5iWz05SKelWCGdsnkjUHkvl6dbzZcYrn4Aosv73I4pxrqT1wAtXDAs1OpDmRPnzoqawxkJUK+xeZnaTU3NG2Njc2r8a7i/ay9+RFs+MUTfIBsqbdw/7cGqy/ZgI3WWuanUhzMl1MPVW9zhAcXqav1c9LRJgQ04Lygb48PW0LGdkecnVU+nmyfxhCamYu48NeZVRfPXxaWaSLqafysRijl+7/HdIvmJ2m1FQJCeCdmJbsOXmRDxfvMztO4ZRC/TwczhzkqZxneGlYbwL99IB4ZZEupp7MGgPZ6bD3V7OTlKoezasxtF0dvvjjIOsOJpsd5+riliD7FzE+6w569hlE42qhZifSXEQXU08W0RbCanvNUX17L9/anLqVyjFy+jYupLvp1VFKkfrbaxxVVTnZ9B7ubFfH7ESaC+li6sl8fCCqPxxYCpfOmJ2mVAUH+PLhkNacvJDOuDk7zY6Tr3NbZhGcvIP/Bt7B24OiEd3Rc5mmi6mns8ZAbjbsmW92klJ3bZ2KPN69Eb9sOcaC7SfMjnOFzKxszi98jUOqBv3vHUlYkJ/ZkTQX08XU09VoDZUaeOWmPsCTNzSiVUQYL83ewakL6WbH+dsvUydTN/sQZ9qNoGnNioXPoHk8XUw9nYixdhr/B6Qkmp2m1PlZfJg4pDXpWTluc3XUjI2HaBP/OaeD6nNt7wfNjqOVEl1MywJrDKhc2DXH7CSmaFA1hJdujeTP/Ul8t/aQqVm2J5xjw9zPaeRznEq3vGqcwqZ5BV1My4Lw5hAe6bWb+gDD2tehe9OqjP91D/tPmXN1VHJKBo9/t56nfX8hO9yKJUp3YuJN3L6YioiPiLwlIp+IyL1m53FbUQOMsaHOJ5idxBQiwjsDWxIc4Muz07eSmV26gw5m5+Ty5I9buD5tCRHqJL49XjbOttC8hkvfbRH5WkQSRSQ2z/29RGSviMSJyOhCFnMbUAvIAryzUjjCOsD4uXO2uTlMFB4ayPgBLYg9doFJS0v36qj3Fu1l04GTjAmZB7WioUmvUm1fM5+rvzq/Aa74rxIRCzAF6A1EAkNFJFJEWojI/Dy3cKApsFYpNQI97lTBKjc0jux78aY+QM+o6gxuE8GnKw6w6VDpnHs7f/txPv/jIB802k5w2gno/pJxYFDzKqKUa49+ikg9YL5Symqb7giMU0r1tE2/CKCUGl/A/MOATKXUdBH5SSk1pIDnPYwxDDXVqlWLnjZtWpFypqSkEBISUqR5nK2kGWofmUXDg9+wrv1npAfVMC2HsxQ3R1q24pXVRifSr3cKIsi3ZIXtajkSLubyxro0GoZk84t6hvSg6mxt/bZLiqmnvy9lIUP37t03K6Xy76lGKeXSG1APiLWbHgh8ZTd9NzD5KvOXA/4DfAI87kib0dHRqqiWL19e5HmcrcQZzh5R6tXySq18z9wcTlKSHBvjk1X90fPVXV+uU4eTUl2S49ylTNX13WWqzZuL1YVlE43X/uAfJWqrODlKmzvkMCsDsEkVUHfM2EOe31d2gavHSqlLSqkHlVJPKqWmuDCX56tQG2p38Kpu+QrSpl4l3uzfgr+OnOXGiSv58Pe9pGU6r8u+3FzFiJ+2knA2jc8HNyN002So3xXqd3FaG5pnMaOYJgC17aYjgOMm5CibrAMgcSck7jE7ienubF+HZSO70SuqOh8vi+PGD1eycMeJy1s8JfLxsv0s3ZPIK30jufbUDEg9DTe87ITUmqcyo5huBBqLSH0R8QfuAOaakKNsiuwP4uMVQ0E7onpYIB8PvYafHu5AaKAvj/3wF3d9tZ59JTgXdenuU3y0ZD8x10Zw9zUVYfUkYwju2u2cmFzzNK4+NepHYC3QVEQSRORBpVQ28ASwCNgNTFdKuWe3P54otJrRC3/sz+Dig4uepH2Dysx/sjOv3xbFzuMX6D3pT16ft6vI3ffFJ6XyzE9bsdYqz1u3W5F1n0HaWbjhJRcl1zyFrysXrpQaWsD9C4GFrmzbq1ljYN7TcHI71Ghldhq34Wvx4Z6O9ejTsibv/76X/1sTz9xtx3i+VzMGXhuBTyEjhaZmZPPI95vw9RE+GxZNYPYFWDsFmvWBmteU0l+huSt9iUZZ1Lwf+Ph6/TmnBakU7M/bt7dg3hOdqVs5mOdnbmfAp2vYdvRcgfMopXj+5+3EJabwydBriahYDtZMhozz0O3FUkyvuStdTMuicpWg4Q3GUX29qV8ga60wZj7akQ8Ht+LYuTT6/3s1L8zcTlJKxj+e++WfB1mw/QTP92pG58ZVIDUJ1n1qXMZb3WpCes3d6GJaVllj4PxRSNhodhK3JiIMuDaCZSO7MrxLA37+K4Hu76/g/1bHk51jXN+/KzmHCb/u4ZYW1Xnk+gbGjKs/guw0vVaq/U0X07Kq6S1gCdDnnDooNNCPMbc057dnrqd17Qq8Nm8Xt368irnbjvPp1nQaVg3h3YGtjKFHLp6EDV9Bi8FQtYnZ0TU3oYtpWRVYHhrfBDtnQa6HjC/vBhqFh/DdA+34/O5oUjOzeerHLWQr+PzuaEICbMdr//wQcjKh2wvmhtXcikuP5msms8YYY0MdXqOvzCkCEaFnVHW6NqnKD+uPkHP6IA2q2q4DP58Am/8PrrnLGC5G02z0mmlZ1qQn+AXro/rFFOhn4cHO9WlS0a63/D/eNw7qXf+cecE0t6SLaVnmHwxNexvDmeS46djynuTsIdjyPUTfBxXqmJ1GczO6mJZ11hhIOwMHV5qdxPOtfNc4f7fLSLOTaG5IF9OyrlEPCAjT1+qXVFIcbPsR2jwI5YvfV6xWduliWtb5BkDzPrB7HmT/82R0zUErJ4BvIHR+1uwkmpvSxdQbWAdAxgWIW2J2Eo8UnHIYdsyE9o9ASFWz42huShdTb1C/K5SrrI/qF1O9Qz9CQChc95TZUTQ3poupN7D4QeRtsPdXyEw1O41nObGdqklrocNjRp8HmlYAXUy9hTUGsi7Bvt/MTuI5zsTDwufI8g2Bjo+ZnUZzc7qYeos6HSGkur5W3xFJcTDrX/BJNBzfwoGG90NgmNmpNDenLyf1Fj4WiLodNn0N6ed1cchP4h74831j37IlANo/Cp2e4uTmPTQzO5vm9nQx9SbWGFj/KexZCK3zHQTBO53aCX+8Bztng185uO5J6Pik3ZF7PTihVjhdTL1JRBsIq2OseeliCie2GVc17ZkP/qHQZQR0eByCK5udTPNAuph6ExHjnNO1kyE12XuLxrHNsPI92PercXVY19HGOaT6aL1WAvoAlLexxkBuNuz2wtG1j26AqTHw5Q1wdJ0xzv2zO6D7i7qQaiWm10y9TfUWULmRsanf5n6z05SOQ6vhj3fh4Arj4oUbx0Hbh4wT8TXNSXQx9TYixtrpyneN4TdCq5udyHUSNsPiV+DwKggOh5vfhDYPGF0TapqT6c18bxQ1AFBGP6dl1cWT8P3tkBwHvd6BZ7YbR+l1IdVcRBdTbxTeDKpZy/a1+gufg+x0uH8hdHgU/ILMTqSVcbqYeivrADi6Hs4dMTuJ8+2eZxxg6zYaKjc0O43mJXQx9VZRA4yfO2eZm8PZ0s7BglHGgbbrnjQ7jeZFdDH1VpXqQ81ry96m/pJXITUR+n1i9JalaaVEF1NvZo0xrgJKPmB2Euc4tAo2fwMdH4ea15idRvMyuph6s6jbjZ9loSeprDSY+xRUrAfdxpidRvNCuph6s7BaUOe6srGpv/JdOHMA+k4C/3Jmp9G8kC6m3s46AE7vhlO7zE5SfCe2w+pJ0HoYNOhmdhrNS+li6u0i+4P4eO7aaU42zH3SuEz05jfMTqN5MV1MvV1IVah/vVFMlTI7TdGt+zec2Aq3vKs7K9FMpYupZhzVPxtvFCVPcuYgLH8bmt5irGFrmol0MdWgWR/w8fOsTX2lYN4zxrmkt35gdOCiaSbSxVQzNo8b9YDYWZCba3Yax2z9AeJXGt3pla9pdhpN08VUs7HGwIUESNhgdpLCXTwFi14yTuuK9pI+WTW3p4upZmjaG3wDPWNT/9fnIesS9PsYfPS/sOYe9H+iZggIhcY3w85ZSG6O2WkKtmcB7JoNXZ+HKo3NTqNpf3P7YioidURkroh8LSKjzc5TplljIPU0YedjzU6Sv/TzsGAkhEdBp2fMTqNpV3BpMbUVwEQRic1zfy8R2SsicQ4UyCbAAqXUA0Cky8JqxpqpfwjhiavMTpK/JeMg5ZTuEUpzS65eM/0G6GV/h4hYgClAb4ziOFREIkWkhYjMz3MLB7YAd4jIMmC5i/N6N/9y0PQWqp5eDZfOmJ3mSofXwKavof2/ICLa7DSa9g+iXHzVi4jUA+Yrpay26Y7AOKVUT9v0iwBKqfEFzD8K2KCU+kNEZiqlBhbwvIeBhwGqVasWPW3atCLlTElJISQkpEjzOJs7ZAhOiSd60whO1ujBvqZPmJrl8uvhk5NJm01PIyqbjW0/IdcSaEoOs+kc5mfo3r37ZqVUm/weM2N00lrAUbvpBKD9VZ7/GzBORO4EDhX0JKXUF8AXAG3atFHdunUrUqgVK1ZQ1HmczR0yQDeOnFpBnaOzqdlrBNS9zrQkf78eS9+AtOMw7Beub9TDvBylICsri4SEBNLT0//xWFhYGIGBpftFkh93yOHqDIGBgURERODn5/juJDOKaX6XqhS4eqyUigXyXRvVXONQvaHUufAXzHsaHl0FvgHmhTkZC6s/glZDjQsLyriEhARCQ0OpV68ekueqrosXLxIaGmpSMvfK4coMSimSk5NJSEigfv36Ds9nxtH8BKC23XQEcNyEHFoBci2B0OdDSNpndG1nFpVj9AgVWAF6vm1ejlKUnp5O5cqV/1FItdIjIlSuXDnfrYOrMaOYbgQai0h9EfEH7gDmmpBDu5rGNxmD7v3xPiTFmRIhImE+HP8Ler/jVT1C6UJqvuK8B64+NepHYC3QVEQSRORBpVQ28ASwCNgNTFdK7XRlDq2Yek0wroqa/0zpd8939hD143+Axj2N8181zc25tJgqpYYqpWoopfyUUhFKqf/Y7l+olGqilGqolHrLlRm0EgitBjeNg0N/wtb/ll672Rkw61GUiLG7Qa+peYzZs2eza1fRR22YPXs2LVu2pFmzZrRo0YLZs2c7PK/ZZxZc5vZXQGkmu/Y+qN0Bfn8ZUpNd355SxoGvI2vZ1+RxCItwfZtakeXk5H/JcXGK6bZt2xg1ahRz5sxhz549zJ07l1GjRrF9+/Z/PDc7O7tYeUuDLqba1fn4QN+PIOMi/P6S69v78wPY9iN0e5HEate7vj3tH6ZOnUq7du1o3bo1jzzyyN+FMyQkhFdeeYXu3buzdu1aRo8eTWRkJC1btmTUqFGsWbOGuXPn8txzz9G6dWsOHDjAgQMH6NWrF9HR0XTp0oU9e/b8o73333+fMWPG/H3kvH79+rz44ou89957AHTr1o0xY8bQtWtXJk2aRHx8PD169KBt27aMHTv27+VcPoVt4MCBNGvWjLvuugtXn0dvz4xTozRPE94cOj0NrOA/pwAADa9JREFUf74Pre5w3aB1O2fBsjegxSDo+gKsXOmadjzEa/N2suv4hb+nc3JysFgsJVpmZM3yvNo3qsDHd+/ezU8//cTq1avx8/Pjscce44cffuCee+4hNTUVq9XKc889R1ZWFg8++CB79uxBRDh37hwVKlSgX79+9OnTh4EDjbMZe/TowWeffUbjxo1Zv349jz32GMuWLbuizZ07dzJq1Kgr7mvTpg1Tpkz5e/rcuXOstP0/9OvXjwcffJBHHnnkiucAbNmyhZ07d1KzZk06derE6tWr6dy5c4leM0fpNVPNMdePgkoNYP6zkFW0U0YckrAZZj0KtdtDv8l6P6lJli5dyubNm2nbti2tW7dm6dKlHDx4EACLxUJMjHEwsHz58gQGBvLQQw/xyy+/UK7cP4fXTklJYc2aNQwaNOjvtdwTJ07843lKqX8cPc9735AhQ/7+ffXq1QwaNAiAu++++4r52rVrR0REBD4+PrRu3ZpDhw4V74UoBr1mqjnGLwj6TITvbjPWUG942XnLPncUfrwDQsJhyA/gZ/5VPu4g7xpkaZwsr5Ti3nvvZfz4f17dHRgY+Peasa+vLxs2bGDp0qVMmzaNyZMn/2ONMzc3lwoVKrB169XHFouKimLTpk20bNny7/v++usvIiP/169RcHDwFfMUdOpSQMD/LjCxWCyluo9Vr5lqjmvQDVreAas+gsR/7vsqloyLRiHNToc7ZxijpWqm6dGjBzNnziQxMRGAM//f3t0HRXHfcRx//6QKltE0CTpqNBUMDyOjgg81oGBwBKuCbYw2as3UaJuSTkFrtNVxhiQOicWoUdTU8SHSGY2KVyeYJlpDPKNGbCNOnBCfKj6BtokPCeIzml//2Ot5wKEct/cA933N7Azs7d5+WH/3dff27rtXrnD27Nl6y127do2qqipGjhzJkiVL7AWzXbt2VFdXA8bRa3h4OFu2bAGMQn348OF6zzVz5kzmz59vP4o8c+YMb775Jq+88orTjIMGDcJisQCwYcMG9/5gE0kxFa4Z/obRSPqDae7fL+r7e2CZCt8chXEF0DHGlIii6Xr27Elubi5paWn07t2b1NRUp6fm1dXVpKen07t3b4YMGcLbb78NwPjx43nrrbeIj4+nvLycDRs2sHbtWvr06UNsbCxFRUX1nisuLo68vDwyMjKIiYkhIyODBQsWEBcX5zTj0qVLWb16NQMGDKCqqqpRf1dOTg7btnn4u0Fa6xY39evXT7vKarW6vI7Z/CGD1o3IcWi91q+21/rgOvc2tH228Tz/Wt20HF7izRxHjhxp8LGrV696LceD+EMOb2Rw9m8BHNQN1B05MhWui5sI3ZPg4xzj5nZN8flaOPCO0Z90wK/NzSeED0gxFa5TyrgYVXMT/jHH9fVPfgIfzTK+KjpcvgAnWgYppqJpwiIhaaZxN9N/Fzd+vW+OwZbJ0CEGxq6FVu59blIIfyHFVDTd4OkQFgUf/gHu3Hj48tcvwXu/MJqnTNxsXMgSooWQYiqa7gfBkL4EvjsHn/75wcvW3IJNE40b4k3YBD/q9uDlhWhmpJgK93QfBPEvwP7l8N8vnS+jNWz7PVT8E55dKTfEEy2SFFPhvtR5RvPmD6Ybnx2t69MF8OUW41tTsc96P5/wmZycHIqLjffUlyxZwo0bjXg7yIHWmtzcXCIjI4mKiiIlJYWvvmpc++Pdu3eTnp7ucuamkmIq3PfDx2D4fDh/0Lgds6MvLbD7TeMeTkkzna8vWqx58+YxbNgwoGnFdMWKFezfv5/Dhw9z4sQJ5syZw+jRo53eUqShtoDeIsVUmKPXWOgxFIpfh6u2W3pVfA7v/w6eTISMpdK8pBlYuXIlcXFxxMXFER4eTkpKCoWFhcyYMQMwvn30/+/Ql5eX2zsylZaWMmTIEPr168fw4cPt35qaPHkyFouF/Px8Lly4QEpKCikpKQDs3LmThIQE+vbty7hx47h27Vq9PHl5eSxbtszeSCUtLY3ExEQKCwuB+20BBw4cSElJCTt27CAmJobBgwezdetW+/O89tprTJkyhWeeeYaIiAjy8/NN33fS6ESYQykYtQjeSYDtf4S0N2DTBGjfBZ5f79s7nDZX22fXeh+67b27EOTmS7ZTLxjR8MXCzMxMMjMzqampYejQocyYMYMBAwbYe4vu3buXRx99lPPnz7Nv3z6SkpKoqakhKyuLoqIiOnTowObNm5k7dy7vvnv/LCU7O5vFixdjtVoJCwvj0qVL5ObmUlxcTGhoKHl5eSxevJicnBz7OlevXuX69ev06NGjVsb+/ftz9OhRAHtbwHnz5nHr1i0iIyPZtWsXTz31VK1OUwDHjh3DarVSXV1NdHQ0L7/8sku3cn4YKabCPI9FGH1IP3ndOCq9dwcmFkLo475OJlw0bdo0hg4dSkZGBmA0NqmurqaiooJx48axZ88e9u7dy5gxYzh+/DhlZWWkpqYCxul2586dH/j8Bw4c4MiRIwwaNAiAO3fukJCQ0Khs2qE9n2NbwGPHjhEeHk5kZCQAkyZNYtWqVfb1Ro0aRXBwMMHBwXTs2JGvv/6arl3Nu5ODFFNhrsQs433SS8dh0t+gQ5SvEzVfdY4gb3rpfvUFBQWcPXuW5cuX2+clJCSwbt06oqOjSUxMZPPmzZSUlLBo0SLOnTtHbGwsJSUljd6G1prU1FQ2btzY4DLt27cnNDSUU6dOERERYZ9/6NAhBg4cCNRuCwgPvquop9vzyXumwlxBreGFrTD1Y8915BceU1paysKFC1m/fj2tWt0vD8nJySxcuJDk5GT69OmD1WolODiYRx55hOjoaC5evGgvpjU1NU6vuDu253v66af57LPPOHnSuI34jRs3OHHiRL11Zs2aRXZ2Njdv3gSguLiYffv22ZtDO4qJieH06dOUl5cDPLBQe4IcmQrztetkTKLZWb58OVeuXLFfJOrfvz9r1qwhKSmJiooKkpOTCQoKolu3bsTEGC0T27Rpg8ViITs7m6qqKu7evcv06dOJja3d3Pqll15ixIgRdO7cGavVSkFBARMmTOD27dsA5ObmEhVV+0wmKyuLb7/9ll69ehEUFESnTp0oKiqibdu29bKHhISwatUqRo0aRVhYGIMHD6asrOyhf/PIkSNZs2YNXbp0adI+s2uonVRznqQFn3skR23Sgq82f8ghLfiEEKKFkmIqhBAmkGIqhJ8xziaFLzXl30CKqRB+JCQkhMuXL0tB9SGtNZcvXyYkxLW75MrVfCH8SNeuXamsrOTixYv1Hrt165bLL3BP8Iccns4QEhLi8gf6pZgK4Udat25NeHi408d2795NfHy8lxP5Zw5/yFCXnOYLIYQJpJgKIYQJpJgKIYQJVEu8aqiUugicdXG1MOCSB+I0twwgOeqSHLX5Qw5fZfix1rqDswdaZDFtCqXUQa11/0DPIDkkR3PI4Q8Z6pLTfCGEMIEUUyGEMIEU0/tWPXwRj/OHDCA56pIctflDDn/IUIu8ZyqEECaQI1MhhDCBFFMhhDCBFFMhhDCBFFMhhDCBFFMnlFKhSqm/KqVWK6V+6es8AEqpn9vyFCml0nycJVQpVaqUSvdhhlZKqTeUUsuUUr/yYY4nlVLblFLvKqVme3nbEUqptUopi8M8r4/dBnJ4fbw6y2Gb75XxGjDF1DbYv1FKldWZ/1Ol1HGl1EmHF8MYwKK1/g0w2svbdkpr/b4tz2TgeV/lsPkTUNiUDCbm+BnwBFADVPowRxTwodZ6CtDTm9vWWp/SWk+tM9ulseupHK6OVw/uD3BzvDZaQ3faa2kTkAz0Bcoc5gUB5UAE0AY4jPGCmAPE2ZZ5z8vb7gX8vc7U0WG9RUBfX+UAhgHjMV4k6T7MMRv4rW1diw9zPA5YgV3Aiz4aExaHn10au57K4ep49eD+cHu8NnYKmObQWus9SqnudWb/BDiptT4FoJTahHHEUwl0Bb7AhKN3V7attZ4P1DsdUUop4M/Adq31IR/mSAFCMQb1TaXUR1rr732QoxK4Y/v1nivbNznHTOBV23NZgHXe2nYDXBq7nsrh6nj14P5we7w2VsCc5jfgCaDC4fdK27ytwHNKqb8AH3h52w3JwvhfdqxSKtNXObTWc7XW04H3gNUmDkxX98dWYLhSahmwx6QMTcmxA8hWSq0Eznhz20qpx23bjVdKzbHNNmPsmpHDjPHqdg4Pjtd6AubItAHKyTyttb4OvOiLbTe0sNY6H8j3dQ77AloX+DKH1voG4Oz9MW/nKAPG+mjbl4HMOvPMGLtm5DBjvLqdw+GxAjezPFSgH5lWAt0cfu8KXAiAbUsO/8wh+8A/czRKoBfTz4FIpVS4UqoNxhvV2wJg25LDP3PIPvDPHI3jyatb/jQBG4H/cP+jNFNt80cCJzCuGs5taduWHP6ZQ/aBf+ZwZ5KuUUIIYYJAP80XQghTSDEVQggTSDEVQggTSDEVQggTSDEVQggTSDEVQggTSDEVAUkp1b1uuzch3CHFVAghTCDFVASyIGV0g/9KKbVTKdXW14FE8yXFVASySGCF1joW+A54zsd5RDMmxVQEstNa6y9sP5cC3X2YRTRzUkxFILvt8PM9pL+vcIMUUyGEMIEUUyGEMIG04BNCCBPIkakQQphAiqkQQphAiqkQQphAiqkQQphAiqkQQphAiqkQQphAiqkQQpjgf9Hr4VvLhgoaAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 360x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "x = np.linspace(-5,5,10)\n",
    "h = []\n",
    "num_abl1 = np.ndarray((16, 10))\n",
    "num_abl2 = np.ndarray((16, 10))\n",
    "\n",
    "for i in range(16):\n",
    "    h.append(10**-i)\n",
    "    num_abl1[i] = ableitung_o1(f, x, h[i])\n",
    "    num_abl2[i] = ableitung_o2(f, x, h[i])\n",
    "    \n",
    "ana_abl = df(x)\n",
    "diff1 = np.abs(num_abl1 - ana_abl)\n",
    "diff2 = np.abs(num_abl2 - ana_abl)\n",
    "\n",
    "diff1 = np.mean(diff1, axis=1)\n",
    "diff2 = np.mean(diff2, axis=1)\n",
    "\n",
    "plt.figure(figsize=(5,5))\n",
    "\n",
    "ax = plt.gca()\n",
    "ax.invert_xaxis()\n",
    "\n",
    "plt.loglog(h, diff1, label='erste Ordn.')\n",
    "plt.loglog(h, diff2, label='zweite Ordn.')\n",
    "\n",
    "plt.title('Genauigkeit genäherter Ableitungen')\n",
    "plt.xlabel('h')\n",
    "plt.ylabel('Genauigkeit')\n",
    "plt.grid()\n",
    "plt.legend()\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Aufgabenteil D\n",
    "\n",
    "1. Was zeigt der Plot? \n",
    "1. Wieso wird die Genauigkeit ab einem bestimmten Punkt für kleinere h wieder schlechter? "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Lösungshinweis\n",
    "\n",
    "zu 2.:<br> Erinnern Sie sich daran, dass Gleitkommazahlen in der Regel gerundet sind und selbst nur eine begrenzte Genauigkeit aufweisen. Diese Genauigkeit liegt bei c.a. $10^{-15}$. Wieso macht sie sich schon bei vergleichsweise großen $h$ bemerkbar? Schauen Sie sich die Formeln nochmal genau an."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "loesung",
     "hide-cell"
    ]
   },
   "source": [
    "### Lösungsvorschlag"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "loesung",
     "hide-cell"
    ]
   },
   "source": [
    "1. Der Plot zeigt, wie für kleinere $h$ sich die Genauigkeit der genäherten Ableitungen bis zu einem gewissen Punkt verbessert und dann wieder verschlechtert. Die Ableitung 2. Ordnung ist, für gleiche $h$, immer Genauer als die 1. Ordnung und erreicht maximale Genauigkeit ungefähr für $h=10^{-5}$ während die Ableitung 2. Ordnung maximale Genauigkeit bei ungefähr $h=10^{-8}$ erreicht.\n",
    "1. Das die Genauigkeit mit kleineren $h$ Werten steigt war erwartet, da die benutzten Formeln $\\mathcal{O}(h)$ und $\\mathcal{O}(h^2)$ haben. Das Problem liegt darin, dass zwei fast gleich große Zahlen voneinander subtrahiert werden. Dies resultiert in einem Verlust an Genauigkeit, der ab einem gewissen Punkt den Ganauigkeitsgewinn durch das Verkleinern von $h$ überwiegt."
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Tags",
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}