From 860eb53cf227d6da96ce035461d8b8dc9266acd2 Mon Sep 17 00:00:00 2001
From: penrose <penrose@mi.fu-berlin.de>
Date: Mon, 17 May 2021 01:46:26 +0200
Subject: [PATCH] nine-point-stencil added, error in boundary
---
nine_point_stencil.ipynb | 482 +++++++++++++++++++++++++++++++++++++++
1 file changed, 482 insertions(+)
create mode 100644 nine_point_stencil.ipynb
diff --git a/nine_point_stencil.ipynb b/nine_point_stencil.ipynb
new file mode 100644
index 0000000..160a682
--- /dev/null
+++ b/nine_point_stencil.ipynb
@@ -0,0 +1,482 @@
+{
+ "cells": [
+ {
+ "cell_type": "code",
+ "execution_count": 1,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import numpy as np\n",
+ "import scipy.sparse as sp\n",
+ "from scipy.sparse.linalg import spsolve\n",
+ "from matplotlib import pyplot as plt"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "**The following poisson problem with Dirichlet boundary condition is given:**\n",
+ "$$-\\Delta u = f \\quad in \\; \\Omega = (0,1)^2$$\n",
+ "$$ u = g \\quad on \\quad \\partial \\Omega \\quad$$"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 2,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "#given exact solution\n",
+ "def u(x,y):\n",
+ " return pow(x,4)*pow(y,5)-17*np.sin(x*y)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 3,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "#right-hand-side of the poisson equation\n",
+ "def fun(x,y):\n",
+ " return -(17*(x**2 + y**2)*np.sin(x*y) + 4*x**2*y**3*(5*x**2 + 3*y**2))"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "**First, the components of the following equation will be assembled:**\n",
+ "\n",
+ "$$A \\underline{u} = \\underline{f} + B\\underline{g} \\Longleftrightarrow A \\underline{u} - B\\underline{g} = \\underline{f} = -\\Delta u $$"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 4,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# matrix A only involves inner nodes. so for m inner nodes, A is an (m x m)-matrix\n",
+ "def sparse_A(h):\n",
+ " N = int(1/h)\n",
+ " m = pow(N-1,2)\n",
+ " A = np.zeros((m,m))+4*np.eye(m) - 2*np.eye(m,k=1) - 2*np.eye(m,k=-1) - 2*np.eye(m,k=N-1) - 2*np.eye(m,k=-(N-1))\n",
+ " A = A + np.eye(m,k=-N) + np.eye(m,k=-(N-2)) + np.eye(m,k=N-2) + np.eye(m,k=N)\n",
+ " A = pow(h,-4)*A\n",
+ " # insert zeros for neighbours of inner nodes that are next to the boundary, because of lexicographical order\n",
+ " for i in range(N-2):\n",
+ " #right boundary, successors:\n",
+ " A[(i+1)*(N-1)-1][(i+1)*(N-1)] = 0\n",
+ " #left boundary, predecessors:\n",
+ " A[(i+1)*(N-1)][(i+1)*(N-1)-1] = 0\n",
+ " #left boundary, diagonally below:\n",
+ " A[(i+1)*(N-1)][(i+1)*(N-1)-(N-2)] = 0\n",
+ " \n",
+ " for i in range(N-1):\n",
+ " #right boundary:\n",
+ " A[(i+1)*(N-1)-1][(i+1)*(N-1)-1-2] = 0\n",
+ " #left boundary\n",
+ " A[i*(N-1)][i*(N-1)+(2)] = 0\n",
+ " #right boundary, diagonally below:\n",
+ " A[i*(N-1)+N-2][i*(N-1)] = 0\n",
+ " #left boundary, diagonally above:\n",
+ " A[i*(N-1)][i*(N-1)+N-2] = 0\n",
+ " for i in range(N-3):\n",
+ " #right boundary, diagonally above:\n",
+ " A[i*(N-1)+N-2][i*(N-1)+N-2+N] = 0\n",
+ " \n",
+ " A = sp.csr_matrix(A)\n",
+ " return A"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 5,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "n = 2\n",
+ "N = pow(2,n)\n",
+ "h = 1/N"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 6,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "matrix([[1024., -512., 0., -512., 256., 0., 0., 0., 0.],\n",
+ " [-512., 1024., -512., 256., -512., 256., 0., 0., 0.],\n",
+ " [ 0., -512., 1024., 0., 256., -512., 0., 0., 0.],\n",
+ " [-512., 0., 0., 1024., -512., 0., -512., 256., 0.],\n",
+ " [ 256., -512., 256., -512., 1024., -512., 256., -512., 256.],\n",
+ " [ 0., 256., -512., 0., -512., 1024., 0., 256., -512.],\n",
+ " [ 0., 0., 256., -512., 0., 0., 1024., -512., 0.],\n",
+ " [ 0., 0., 0., 256., -512., 256., -512., 1024., -512.],\n",
+ " [ 0., 0., 0., 0., 256., -512., 0., -512., 1024.]])"
+ ]
+ },
+ "execution_count": 6,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "A = sparse_A(h).todense()\n",
+ "A"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 7,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "def vector_f(h):\n",
+ " N = int(1/h)\n",
+ " l = pow(N-1,2)\n",
+ " f = np.zeros(l)\n",
+ " for i in range(N-1):\n",
+ " for k in range(N-1):\n",
+ " f[k+i*(N-1)]=fun((k+1)/(N),(i+1)/(N)) \n",
+ " return f"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 8,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# for m inner nodes and l boundary nodes, B is a (m x l)-matrix\n",
+ "def sparse_B(h):\n",
+ " N = int(1/h)\n",
+ " m = pow(N-1,2)\n",
+ " l = 4*N\n",
+ " B = np.zeros((m,l))\n",
+ " \n",
+ " # since the lower and left boundary values are zero, only entries for the right\n",
+ " # and upper boundary nodes are needed\n",
+ " for i in range(N-1):\n",
+ " #right boundary, successor:\n",
+ " B[(i+1)*(N-1)-1][N + 2 * (i + 1)] = -2*pow(h,-4)\n",
+ " #upper boundary:\n",
+ " B[-(N-1)+i][-N+i] = -2*pow(h,-4)\n",
+ " #right above\n",
+ " B[-(N-1)+i][-N+i+1] = pow(h,-4)\n",
+ " #left above\n",
+ " B[-(N-2)+i][-N+i] = pow(h,-4)\n",
+ " for i in range(N-2):\n",
+ " #right boundary, diagonally above:\n",
+ " B[(i+1)*(N-1)-1][N + 2 * (i + 1)+2] = pow(h,-4)\n",
+ " #right boundary, diagonally below:\n",
+ " B[(i+2)*(N-1)-1][N + 2 * (i + 1)] = pow(h,-4)\n",
+ " \n",
+ " B = sp.csr_matrix(B)\n",
+ " return B"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 9,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "matrix([[ 0., 0., 0., 0., 0., 0., 0., 0.,\n",
+ " 0., 0., 0., 0., 0., 0., 0., 0.,\n",
+ " 0., 0., 0., 0., 0., 0., 0., 0.,\n",
+ " 0., 0., 0., 0., 0., 0., 0., 0.],\n",
+ " [ 0., 0., 0., 0., 0., 0., 0., 0.,\n",
+ " 0., 0., 0., 0., 0., 0., 0., 0.,\n",
+ " 0., 0., 0., 0., 0., 0., 0., 0.,\n",
+ " 0., 0., 0., 0., 0., 0., 0., 0.],\n",
+ " [ 0., 0., 0., 0., 0., 0., 0., 0.,\n",
+ " 0., 0., 0., 0., 0., 0., 0., 0.,\n",
+ " 0., 0., 4096., 0., -8192., 0., 4096., 0.,\n",
+ " 0., 0., 0., 0., 0., 0., 0., 0.],\n",
+ " [ 0., 0., 0., 0., 0., 0., 0., 0.,\n",
+ " 0., 0., 0., 0., 0., 0., 0., 0.,\n",
+ " 0., 0., 0., 0., 0., 0., 0., 0.,\n",
+ " -8192., 4096., 0., 0., 0., 0., 0., 0.],\n",
+ " [ 0., 0., 0., 0., 0., 0., 0., 0.,\n",
+ " 0., 0., 0., 0., 0., 0., 0., 0.,\n",
+ " 0., 0., 0., 0., 0., 0., 0., 0.,\n",
+ " 4096., -8192., 4096., 0., 0., 0., 0., 0.],\n",
+ " [ 0., 0., 0., 0., 0., 0., 0., 0.,\n",
+ " 0., 0., 0., 0., 0., 0., 0., 0.,\n",
+ " 0., 0., 0., 0., 0., 0., 0., 0.,\n",
+ " 0., 4096., -8192., 4096., 0., 0., 0., 0.],\n",
+ " [ 0., 0., 0., 0., 0., 0., 0., 0.,\n",
+ " 0., 0., 0., 0., 0., 0., 0., 0.,\n",
+ " 0., 0., 0., 0., 0., 0., 0., 0.,\n",
+ " 0., 0., 4096., -8192., 4096., 0., 0., 0.],\n",
+ " [ 0., 0., 0., 0., 0., 0., 0., 0.,\n",
+ " 0., 0., 0., 0., 0., 0., 0., 0.,\n",
+ " 0., 0., 0., 0., 0., 0., 0., 0.,\n",
+ " 0., 0., 0., 4096., -8192., 4096., 0., 0.],\n",
+ " [ 0., 0., 0., 0., 0., 0., 0., 0.,\n",
+ " 0., 0., 0., 0., 0., 0., 0., 0.,\n",
+ " 0., 0., 0., 0., 0., 0., 0., 0.,\n",
+ " 0., 0., 0., 0., 4096., -8192., 4096., 0.]])"
+ ]
+ },
+ "execution_count": 9,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "# example for B\n",
+ "n = 3\n",
+ "h = pow(2,-n)\n",
+ "B = sparse_B(h)\n",
+ "B.todense()[-10:-1]"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 10,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# vector g of length 4*N contains the boundary values which satisfy u(x,y)\n",
+ "def vector_g(h):\n",
+ " N = int(1/h)\n",
+ " l = 4*N\n",
+ " g = np.zeros(l)\n",
+ " \n",
+ " # upper boundary, where y=1\n",
+ " for i in range(N):\n",
+ " g[-N+i] = u((i+1)*h,1)\n",
+ " # right boundary, where x=1 \n",
+ " for i in range(N-1): \n",
+ " g[N+2+2*i] = u(1,(i+1)*h)\n",
+ " return g "
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 11,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "array([ 0. , 0. , 0. , 0. ,\n",
+ " 0. , 0. , -4.20489074, 0. ,\n",
+ " -8.11898416, 0. , -11.35055423, 0. ,\n",
+ " -4.20196106, -8.08773416, -11.27145267, -13.30500674])"
+ ]
+ },
+ "execution_count": 11,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "# example for g\n",
+ "n = 2\n",
+ "h = pow(2,-n)\n",
+ "g = vector_g(h)\n",
+ "g"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 12,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "n = 3\n",
+ "h = pow(2,-n)\n",
+ "N = pow(2,n)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 13,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "A = sparse_A(h)\n",
+ "f = vector_f(h)\n",
+ "B = sparse_B(h)\n",
+ "g = vector_g(h)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "**Now, the linear system Au = RHS will be solved.**"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 14,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "RHS = f+B.dot(g)\n",
+ "appr_u = sp.linalg.spsolve(A,RHS)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "#### **Error and Consistency**"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 15,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "def exact_solution(h):\n",
+ " N = int(1/h)\n",
+ " l = pow(N-1,2)\n",
+ " v = np.zeros(l)\n",
+ " for i in range(N-1):\n",
+ " for k in range(N-1):\n",
+ " v[k+i*(N-1)]=u((k+1)/(N),(i+1)/(N)) \n",
+ " return v"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 16,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "grid = np.zeros(7)\n",
+ "error_eukl = np.zeros(7)\n",
+ "error_inf = np.zeros(7)\n",
+ "\n",
+ "for i in range(5):\n",
+ " grid[i]=pow(2,-(i+2))\n",
+ " h = grid[i]\n",
+ " A = sparse_A(h)\n",
+ " f = vector_f(h)\n",
+ " B = sparse_B(h)\n",
+ " g = vector_g(h)\n",
+ " RHS = f+B.dot(g)\n",
+ " appr_u = sp.linalg.spsolve(A,RHS)\n",
+ " v = exact_solution(h)\n",
+ " x = (appr_u - v)\n",
+ " #print(appr_u[-1],v[-1])\n",
+ " error_eukl[i] = np.linalg.norm(x)\n",
+ " error_inf[i] = np.linalg.norm(x,ord = np.inf)\n",
+ " "
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 17,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "Text(0.5, 0, 'step size h')"
+ ]
+ },
+ "execution_count": 17,
+ "metadata": {},
+ "output_type": "execute_result"
+ },
+ {
+ "data": {
+ "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY8AAAEPCAYAAAC6Kkg/AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAAvJ0lEQVR4nO3deXyU9bn38c+VHRIIJEFlC0kQrRZFZKcKtj0ebbXa02qLWq0VBVzOaWuf56in7aM9XfR0sbVHj4iK+1K1tcWtHndsBWSpCygqBpCAihB2JITkev64Bxgmk3APM5OZJN/36zWvZO6Ze+4rreSb3/3bzN0RERFJRE6mCxARkY5H4SEiIglTeIiISMIUHiIikjCFh4iIJEzhISIiCcvLdAHtpaKiwquqqjJdhohIh7Jw4cJ17t4n9niXCY+qqioWLFiQ6TJERDoUM1sZ77huW4mISMIUHiIikjCFh4iIJEzhISIiCVN4iIhIwhQeIiKSsC4zVDcMd2fbzia27tjFlh2NbI58PfSgEgb07p7p8kREskaXCY+1Wxq47qmlbNnRyNaGXWyJBEPwddee481xtjcpLsjl1vNGMv7QivYvXEQkC1lX2QyqsO8QH3TB7+lRlEePojxKivLoUZi/5/ueRfl7Xyvc+31BXg4/mfUWy9dt43eTjuHLR/XN9I8iItJuzGyhu49scbyrhMeIESN94cIDm2G+aXsjk++az8IPNvCfpw/l3LGDUlydiEh2ai08ukyHudmBn1vaPZ97Jo/hC4cfxI//vJgbnn2PrhK6IiLxdJnwSFa3glymnzuCrx87gN8++y5Xz1pCU7wOEhGRLiB0h7mZnQycCQwEimJednefmMrCslF+bg6/PvNoKkoKuGV2Leu37eT6bwyjMC8306WJiLSrUOFhZv8OXAd8AiwDdqazqGxmZlz15SMoKy7g2qeWsml7I9PPHUFJYZcZuCYiErrlcRlwC3CZuzelsZ4OY+rEwZSXFHLFH9/g7Fvncsf5oygvKcx0WSIi7SJsn0dP4GEFx77OGDGAGeeO4J2PtnDm9DnUbdie6ZJERNpF2PB4GhibzkI6qi8ecTD3XjiGdVsb+PrNr/Dux1syXZKISNqFDY/LgC+b2VVmNsLMamIf6Swy242qKuOhaeMAOHP6HBaurM9wRSIi6RU2PBzYAvwceBV4L86jS/vMIT15ZNp4yooLOOe2eTy/9ONMlyQikjZhO8zvBMYDvwWW0oVHW7VlYFl3Hp42ju/cMZ+L7l7IL79+NF8fMSDTZYmIpFzY8DiBYKTVnekrJTwzOxK4BlgPPOfuj2S2or0qSgp5YMpYpt6zgB88/Dr123Zy0YQufVdPRDqhsLet1gFpvQ9jZjPNbK2ZLY45frKZvWNmy8zsysjhLwH/7e4XA+els64DUVKYx8zzR3HKUX35+ZNvc+1Tb2s5ExHpVMKGx++BS8wsncuZ3AmcHH3AzHKBmwjC4kjgrEir4x5gkpn9CihPY00HrDAvl9+fNZxvja3klpdq+fdH3mBXU3OmyxIRSYmwt616A0OBt8zsGWBDzOvu7lcnU4i7zzazqpjDo4Fl7l4LYGYPAqe7+7XApZFw+VMy102n3Bzjp6cPpaKkkN89+x4btjdy49nDKcrXciYi0rGFDY8fRn1/WJzXHUgqPFrRH1gV9bwOGBMJmf8AioFftXaymU0BpgBUVlamobz9MzO+90+HUV5cwP+btYRzb5/Hbd8eRWm3/IzUIyKSCqFuQ7l7zn4e6fpTOt5C6u7uK9x9iruf4+5/a6PuGe4+0t1H9unTJ00lhnPuuCr++6zhvLZqI9+8ZQ4fb96R0XpERJKx3/AwswIze9TMJrRHQTHqCFbx3W0AsCYDdaTEqUf3447zR/NB/Xa+fvMrLF+3LdMliYgckP2Gh7vvBP4pzHvTYD4wxMyqzawAmATMykAdKXPckAoenDKW7TubOOPmV1i8elOmSxIRSVjYQPg7aV7bysweAOYAh5tZnZlNdvddBEujPA28DTzk7kvSWUd7OHpALx6ZNo6i/FwmzZjLK8vWZbokEZGEhNrD3Mw+C/wZuCHy9UOCTvI93D2rx6GOHDnSFyw4sD3M0+WjTTs4b+Y8Vqzbzu8mHcOXj+qb6ZJERPaR7B7mbwKDCcJjJcHyJI1RDy1XcgAOKS3ioanjOGpAKZfev4h7567MdEkiIqGEHar7n8S0NCQ1enUv4N7JY7j0/kX86M+LWb91J//2xUMxizfQTEQkO4QKD3e/Js11dGndCnK55dwRXPHHN/jts++yflsD13zls+TkKEBEJDslvPG2mZUQzDivd3eNNU2R/Nwcfn3GMCpKCpkxu5b6bTu5/hvHUJCXiUFuIiJtC/2bycxOMrMFwEZgBbDJzF41sxPTVFuXk5Nj/MeXj+CqL32Gx9/4kMl3zWdrw65MlyUi0kKo8DCzk4AngBLgp8AlwM+AHsCTCpDUmjpxML8842heeX8959w6l/VbGzJdkojIPsIO1Z1DsBjiqdFDciOr7D4O9HL38WmrMgWycaju/jz71sdcev8i+vfuxt0XjGZA7+6ZLklEuphkh+oOA26KncsRef4/wDFJVygt/NORB3PvhWP4ZEsDZ9w8h3c/3pLpkkREgPDh0QD0bOW1HpHXJQ1GVZXx0NRxNLtz5vQ5LFxZn+mSRERCh8eLwE/NrDr6oJlVEmwH+0Jqy5JoR/TtyR8vHk9ZcQHn3DaPF5auzXRJItLFhQ2PK4BS4B0zm21mfzCzl4D3gF6R1yWNBpZ15+Fp4zj0oBIuvHsBf1pUl+mSRKQLC7ufx7vA0QTb0RYCxwJFBMuVHOPu76WtQtmjoqSQBy4ay5jqMi5/6HVue7k20yWJSBcVepKgu38I/J801iIh9CjK547vjOL7f3iNnz3xNuu27uSKkw/XciYi0q4SnmEumVeYl8t/n3UsvbsvZvpL71O/rYFf/MtR5OVqNrqItI/Q4WFm3wbOAioJbllFc3cfnMrCpG25OcbPvjqUipJCbnjuPeq3NXLj2cMpyk/XjsAiInuFCg8z+zHwE2Ax8BoampsVzIzvn3gY5SUFXD1rCefd/iq3fnskpd3yM12aiHRyYVsek4Eb3P376SxGDsx546ro3b2Ayx96jW/eMoe7LxjNQT1jG4ciIqkT9iZ5OfBYOguR5HxlWD9mnj+KD+q387WbX2H5Oi14LCLpEzY8XiJYokSy2PFD+vDARWPZvrOJM6e/wuLVmzJdkoh0UmHD43vAd8zsPDOrMLOc2Ecaa2zBzI43s+lmdpuZvdKe1852wwb24uFp4yjMy2XSjLm8smxdpksSkU4o7C/9d4GhwB3Ax+y7f3lK9jA3s5lmttbMFsccP9nM3jGzZWZ2JYC7v+zu0whW9L0r2Wt3NoP7lPDHi8fTr1cR598xn6fe/DDTJYlIJ5NNe5jfCdwI3L37gJnlAjcBJwJ1wHwzm+Xub0XecjZwYZrr6pAOKS3ioanjmHzXAi65fxE/++pQzhkzKNNliUgnkTV7mLv7bDOrijk8Gljm7rUAZvYgcDrwVmRRxk3uvjndtXVUvboXcO/kMVx6/yJ++Ohi1m/dyb9+4VDNRheRpB1QX4WZTTCz4lQXE0d/YFXU87rIMQiGD9/R1slmNsXMFpjZgk8++SRNJWa3bgW53HLuCL42vD/XP/Mu18xaQnNzuhuRItLZJbw8SeRW0gvAKGBRyiuKuVycYw7g7lfv72R3nwHMgGAnwdSW1nHk5+bw6zOHUV5SwK0vL6d+eyO/OXMYBXlazkREDsyBrm3VXvc96oCBUc8HAGva6dqdSk6O8cNTjqSipJBrn1rKxu07mf6tERQXankzEUncgf7p2V5/xc8HhphZtZkVAJOAWe107U5p6sTB/PKMo/n7snWcfetc6rclPVBORLqgAw2PlLc8zOwBYA5wuJnVmdlkd98FXAY8DbwNPOTuS1J97a7mGyMHcsu5I1n60RbOmP4Kqzd+mumSRKSDMfeu0RUwcuRIX7BgQabLyCqvLq9n8l3zKS7I4+7Jozns4B6ZLklEsoyZLXT3kbHH1WPahY2uLuOhqeNocufM6XNYuHJDpksSkQ4iVHhEliCZZmbPmdm7ZvZBzGNluguV9Diib0/+dPF4enfP55zb5vLCO2szXZKIdABhWx6/BP4H6EXQif1czOP5dBQn7WNgWXcenjaewX1KuOiuBTz6j7pMlyQiWS7sOM1vAT8NM7dCOqY+PQp5cMpYpty9kO//4XXWb93JhcfXZLosEclSYVseecDsdBYimdejKJ87vjOKLw09hJ898Tb/9deldJUBFSKSmLDh8QhwUjoLkexQlJ/LjWcfy9ljKrn5xfe54o9vsKupOdNliUiWCXvb6nLgPjObQTDnosWwHHdXv0cnkZtj/PyrQ6koKeT3z71H/bZGbjx7OEX5uZkuTUSyRNjw6AvUEKxoG70EuhNMGHRAv1k6ETPj8hMPo7y4gGseW8J5t7/Krd8eSWm3/EyXJiJZIGx43AFUAN8FlpKCzZ+kY/j2+Cp6Fxfwg4de45u3zOHuC0ZzUM+iTJclIhkWNjxGAue5+yPpLEay02nD+tGrWz7T7l3I16e/wj0XjKGqoj1W5BeRbBW2w/wD1Nro0iYc1of7LxrL1h27OGP6KyxevSnTJYlIBoUNj58BV5hZSTqLkex2zMBePDxtPIV5uUyaMZdX3l+X6ZJEJEPChsdJBHtprDCzx8zs7pjHXWmsUbLIoQeV8MjF4+hbWsT5M+fz18UfZrokEcmAsOFxHNAMbAGGAsfHeUgX0be0Gw9PG8fQ/j255L5F3D/vg0yXJCLtLFSHubtXp7sQ6Vh6dS/g3gvHcMl9i/iPR99k/dYGLvvCoZi11yaTIpJJWpJdDlj3gjxuPW8k/zK8P7955l1+8thbNDdrORORriBUy8PMKvf3HnfXvYsuKD83h9+cOYyy4gJu/9ty6rft5NdnDqMgT3+XiHRmYed5rGD/+5ZrhnkXlZNj/OiUI6goKeS//rqU1Rs/5fITD2P84HLdxhLppMKGxwW0DI9y4BSCZUt+msqipOMxMy4+YTCHlBby8yeWcs5t8ziqfylTJ9Zw8mcPIS9XLRGRziTpPczN7B5gpbv/KDUlhbrmCQSBtQR40N1f3N852sO8/exobOLP/1jNjNm11K7bRmVZdy46vpozRgykW4EaqCIdSTr3ML+XoGWSFDObaWZrzWxxzPGTzewdM1tmZldGDjuwFSgCtO1dlinKz2XS6EqeuXwi0781grLiAn78lyV87r+e5/fPvceGbVqsQKSjS0XL41zgBncvS/JzJhAEwt3uPjRyLBd4FziRICTmA2cBS9292cwOBq5393P29/lqeWSOuzN/xQZueel9nlu6lm75uXxz1EAmH1fNwLLumS5PRNrQWssj7GirCXEOFxBMGLwKeDm58sDdZ5tZVczh0cAyd6+N1PEgcLq7vxV5fQNQ2EbdU4ApAJWV+x0wJmliZoyuLmN0dRnvfLSFGbNruXfuSu6Zu5JTj+7L1AmDObJfz0yXKSIJCNXyMLNmWnaY7x5G8xJwjruvSbqYIDwej2p5nAGc7O4XRp6fC4wBnidYMqUXcLP6PDqeNRs/5Y6/L+f+eR+wbWcTEw7rw7QJNYzTCC2RrJJUywP4fJxjOwg6yj9KqrK2xfst4u7+J+BPabyupFm/Xt344SlHctkXhnDv3JXc8fcVnB01QutLQ/uSm6MQEclWYZcneSndhbSiDhgY9XwAkHQLR7JHabd8Lv38oUw+rppHIyO0Lrv/H1SWvcNFE2o4c8QAbX8rkoUOaLSVmeXEPlJdWMR8YIiZVZtZATAJmJWma0kGFeXnctboSp6NHqH158V87jqN0BLJRqF+6ZtZNzO7zszeN7MGoDHmkfS/bDN7AJgDHG5mdWY22d13AZcBTwNvAw+5+5JkryXZKzfHOHnoITx6yXj+MGUswwb24vpn3mX8dc/zk8eWULdhe6ZLFBHCd5jfAZwDPEYre5i7+09SXl0KqcO843rnoy3cMvt9Zr22Bge+cnRfpmiElki7aK3DPGx4rAd+4u6/T0dx7UHh0fGt2fgpM/+2nAdejRqhNbGGcTUaoSWSLsmGxxrg2+7+TDqKaw8Kj85j0/ZG7p23kjv+vpx1W3dy9IBSpk4YzMlDD9EILZEUSzY8fgEc7O6T01Fce1B4dD47Gpv406LV3PpyLcvXbWNQeXcuPF4jtERSKdnwyAVuBqoIOq83xL7H3WcmX2b6KDw6r6Zm55m3PuLml2p5fdVGyosLOH98FeeOG0Sv7gWZLk+kQ0s2PEYTDJE9qJW3uLtn9Z96Co/Oz92Zt7yeW156nxfe+YTuBcEaWhceX0P/Xt0yXZ5Ih5RseCwiWEPqSlofbbUyBXWmjcKja1n60WZmzK7dM0LrtGH9mDKhhiP6aoSWSCKSDY/twBnu/mQ6imsPCo+uaXXUCK3tO5uYeFgfpmqElkhoye7n8Q5QnNqSRNKvf69u/PjUI5lz5Rf5vycdzpI1mzj71nl89aa/8+SbH9LUnNyWBCJdVdiWx0nAL4HTsv32VGvU8hAIRmj9cVEdt86uZcX67Qwq785Fx9dwhkZoicSV7G2rl4FDgTKCzZliR1u5u09MRaHpovCQaE3Nzv8u+YjpL73P63WbNEJLpBXJhseLtNzPYx/uHm/Z9qyh8JB43J25tfXcMvt9XoyM0Jo0qpLJx1drhJYISYZHZ6DwkP15+8PN3Dq7llmvB6v+nzasH1Mm1vCZQzRCS7ouhYfCQ0JavfFTbn95OQ/OD0ZonXB4H6ZOGMzYmjKN0JIuJ+HwiOxbvsjdt7ayh/k+3H128mWmj8JDErVx+849uxyu37aTYQNKmTpxMCd9VmtoSddxIOHRDIx191db2cN8z1vRDHPpxHY0NvHIwjpufbmWleu3U1XenYsm1PD1YzVCSzq/AwmPicDCSMtjvyOpMrhVbSgKD0lWU7PzdGSE1ht1m6goiYzQGltFaff8TJcnkhbq81B4SIrsHqE1/aX3eendYITWWaMrueA4jdCSzifZobrD3P31Nl4/090fTrLGtFJ4SDq8/WFkDa3X12BohJZ0PsmGxw7gCne/IeZ4d+BGgo2isvrmr8JD0qluw3Zm/m3FnhFanz+8D1MnDmZMtUZoSceWis2g/p1gL4/z3f0TMzsWuB/oC1zi7veluOa26jkC+C5QATzn7jfv7xyFh7SHjdt3cs+cldz5SmSE1sBeTJtQwz9rhJZ0UEn3eZjZ54F7gFzgAeAS4DXgbHevTUGBM4FTgbXuPjTq+MnADZHr3ubu10W9lgPcGmaHQ4WHtKfYEVrVFcVcdHwNXzu2v0ZoSYeSkg5zMzsSWAgUAAuA8e7elKICJwBbgbt3h0dkB8N3gROBOmA+cJa7v2VmpxHsL3Kju9+/v89XeEgmNDU7f10cjNB6c3UwQus7n6vmW2MGaYSWdAipaHmcBNwJ7AKeAC4CHgcucPf1KSqyCng8KjzGAde4+0mR51cBuPu1Uec84e6ntPJ5U4ApAJWVlSNWruyQCwJLJ+DuzKldz/SXapkdNUJr8nHV9NMILclirYVHXsiTfwN8j2Ar2snuXm9mjwB3AW+Y2bnu/nwqC47oD6yKel4HjDGzE4CvEexu2OoGVe4+A5gBQcsjDfWJhGJmjB9cwfjBFby1ZjMzZr/Pna+s4K5XVnDS0EM47tAKRleXUVNRrA526RBChQdwMXBZdMe0uz9rZkcDdwD/m8BnJSLevyJ39xeBF9NwPZG0O7JfT343aTj/56TDue3l5Tz+xhqeeONDACpKChlTXcbo6jLG1JRx2EE9yFFHu2ShsL/wR7n7ktiDkdtVp5nZpakta486YGDU8wHAmjRdS6RdDejdnWtO+yxXf+VIatdt49Xl9cyrXc+85fU88WYQJqXd8hlVVcbYmiBQjuzbk7zcsBuAiqRPqPCIFxwxr9+UmnJamA8MMbNqYDUwCTg7TdcSyQgzY3CfEgb3KeGs0ZW4O3UbPmXe8npeXb6eV5fX8+zbHwNQUpjHiEG9g5ZJdRlHDSilME+jt6T9JXSrycx6A0OAotjXkl1V18weAE4AKsysDrja3W83s8sI5pfkAjP3F2QiHZ2ZMbCsOwPLunPGiAEAfLx5x54wmVdbz6+efgeAwrwcjq3svec21/CBvelWoDCR9As7SbAImAl8g/j9EGiGuUj7Wb+1gfkrNgS3upav560PN+MO+bnG0QN67ek3GTGoNz2KNCRYDlyyM8x/DnwH+L8EEwUvBXYA5xPMMP+uuz+VyoJTTeEhndmmTxtZtHIDcyO3ud6s28SuZifH4LP9SveEyaiqMnoXa492CS/Z8FgK/A64FWgERrr7oshrDwNr3P27Ka04xRQe0pVs37mLRSs38ury9cxdXs9rqzayc1czAJ85pAejI2EyurqMg3q0uAstskdS8zyASmCJuzeZWSNQHPXaTILhulkdHiJdSfeCPI4bUsFxQyqAYLmUN+o2Ma92Pa+uqOeRhXXcPSeYNFtTUbynz2R0dbmWlZdQwobHeqAk8v0qYBjwcuR5BaD/2kSyWFF+7p6WBkBjUzNL1mwOwiQyNPjB+cF83P69ujGmpixyq6ucqvLumrgoLYQNj7nAcOAp4I/AT82sB8FSJT8A/pae8kQkHfJzczhmYC+OGdiLqRMH09TsLP1oc2SuST0vvvMJf1q0GoCDehTuGRo8pqacQ/uUaOKihO7zGAlUuvufIqFxJ3AawfDZucAkd/8gnYUmS30eIuG5O+9/spW5tfV7RnR9vLkBgN7dg4mLY2rKGVNdxhF9e2q5+U4s5dvQmlkhUOjum5Mtrj0oPEQOnLvzQf32yFyTIExW1X8KQI/CPEZW9WZ0dTmjq8s4qn8pBXmaBd9ZJNth3oK7NwANSVUlIh2CmTGovJhB5cV8Y2SwYtCajZ8yf0V9pHWynhfe+QSAbvm5HDuoF6OrgjAZXtlLe5h0QoksyV4GnEKw1lTs2D5396tTXFtKqeUhkl6fbGlg/ordLZN6ln4UTFwsyM1h2MDSSL9JOccO6k1JYTrWUZV0SHaexz8TdJQXt/IW1wxzEYm2aXtjECYrggUfF6/ZTFOzk5tjDO3XkzE15YyuCiYuamOs7JVseCwG6glmli9198bUl5heCg+RzNrasItFKzcwLzIL/vVVm9jZ1IwZfOaQnntmwY+uLqOipDDT5UpEsuGxFfgXd38mHcW1B4WHSHbZ0djEa6s2Mq+2nldXrGfhyg3saAxmwQ/uU8zo6vI9S9H3LdVUskxJtsP8H0C/1JYkIl1ZUX4uY2vKGVtTDgxh565m3ly9ac9orsdeX8MDrwYzAAaWdWNMZDTXMQN7UVZcQM+ifI3qyqCwLY9RBHM7LnT3OekuKh3U8hDpWJqanbc/3MzcyCz4V1fUs3H7vnfMi/Jz6FmUT89u+fQsyot8zadnt7yo47HP975P4bN/ybY8FgLPAX8zs23AxpjX3d0HJVeiiMheuTnG0P6lDO1fyoXH19Dc7Cz7ZCtvrdnM5h2NbP60kc07dkW+NrL5013Ub9vJ8nXb9rzW1Nz2H8e7w6e0W+IB1KOLh0/Y8Pg1cBnB7aulwM60VSQiEkdOjnHYwT047OAeod7v7mzf2bQnWPYGTuR59Pc7gu/Xbd1JbQLh0y0/t9VWTZjWT34H3lI4bHicD/w02+dyiIjsZmYUF+ZRXJhH39LEz2+P8OlekNvmbbW2AqhHUV5GwydseDiQ1DazIiIdSSbC55MtDSxbu3XPe/eTPaHD50tD+6Z8Lk3Y8HgY+BJBv0fGmVkN8EOg1N3PyHQ9IiKxUhE+23Y2tRI4Lft7Nu9oZO2WHSxbu6tF+IytKc9YeDwF/NbMSoG/Ahti3+DuzydTiJnNBE4F1rr70KjjJwM3EKzge5u7X+futcBkM3skmWuKiGQrM6OkMI+Swjz6HcCWSdHh06dH6iddhg2PRyNfJ0ceuzlgka/JLk9yJ3AjcPfuA2aWC9wEnAjUAfPNbJa7v5XktUREOrXo8EmHsJ/6+bRcPYq7zzazqpjDo4FlkZYGZvYgcDqg8BARyaBQ4eHuL6W7kFb0J9j2drc6YIyZlQM/B4ab2VXufm28k81sCjAFoLKyMt21ioh0Gdm+LnK87cnc3dcD0/Z3srvPAGZAMMM8xbWJiHRZ2T5DpY5g/5DdBgBrMlSLiIhEZHt4zAeGmFm1mRUAk4BZGa5JRKTLy5rwMLMHgDnA4WZWZ2aT3X0XwbIoTwNvAw+5+5JM1ikiIlnU5+HuZ7Vy/EngyXYuR0RE2pA1LQ8REek4Qrc8zGwicBZQCRTFvOzu/sVUFiYiItkrVHiY2VTgZmA98B7QEPuWFNclIiJZLGzL4wfA/cAF7q69PEREuriwfR79gTsUHCIiAuHDYyFQk85CRESk4wgbHv8GfM/MJqSzGBER6RjC9nk8BvQEXjCz7bTcz8PdfVBKKxMRkawVNjyeI9izQ0REJPSS7OenuQ4REelANMNcREQS1mrLw8zOA55w9/WR79vk7nfv7z0iItI5tHXb6k5gLMGs8jv38zlO1N7jIiLSubUVHtXAh1Hfi4iIAG2Eh7uvjPe9iIiIOsxFRCRhCg8REUmYwkNERBKm8BARkYTtNzzMLNfMhplZn/YoKAwzqzGz283skUzXIiLSFYVpeTiwABieigua2UwzW2tmi2OOn2xm75jZMjO7ss2C3GvdfXIq6hERkcTtd20rd282s1VAcYqueSdwI1GTCs0sF7gJOBGoA+ab2SwgF7g25vwL3H1timoREZEDEHZV3VsI9vN4ItndBN19tplVxRweDSxz91oAM3sQON3drwVOTeZ6IiKSemHDowcwGKg1s78SzDyPXqLd3f3qJOroD6yKel4HjGntzWZWDvwcGG5mV0VCJt77pgBTACorK5MoT0REooUNj/+I+v6COK87kEx4WCufGZe7rwem7e9D3X0GMANg5MiR2o9ERLqOHZthUx1sXg3VEyCvMKUfH3Y/j3QP6a0DBkY9HwCsSfM1RUQ6psZPYfMa2LQKNq0OAmJ3UOx+3rB57/svnQ99DktpCWFbHuk2HxhiZtXAamAScHZmSxIRyYCmRtjyYRACm+pgc13LgNi+vuV5xX2gZ38oHxy0NEoHQGl/6DkAeg1s+f4kJRQeZnYqMBEoI1iq/SV3fyLBz3gAOAGoMLM64Gp3v93MLgOeJhhhNdPdlyTyuSIiWa+5GbatjYRBXRAGe76PBMSWj2hx176oNAiB0v7Qf0TwtXRgEBal/aFHP8gvatcfJVR4mFkP4HHgeGAXQXCUAz8ws5eBU919a5jPcvezWjn+JPBkmM8QEck67vDphqjbR3Uxt5LqYPOH0Ny473n53feGwOAvRoJhQORY5GthSWZ+pjaEbXn8AjgWOBd40N2bInMzJgE3R17/t/SUKCKSBRq2tGwlxAbErk/3PScnH3r2C0Jg4NjIbaRIKOwOhm69weKNGcpuYcPj68CP3P2+3QfcvQm4z8wqgH9H4SEiHVXjjiAAolsJm2I6oRs27XuO5UDJIUEgHDwUDjt5bwuiZyQcivtATudcQjBseJQDb7Xy2luR10VEsk/TrqADOt6IpN0th+3rWp7XvSIIgt7VUHXcvreRSvtDj76Qm9/+P0+WCBseywlmej8T57UvR14XEWlfzc2w7ZP4I5L2dEB/CN6873mFPfcGQb9j9nZG7wmIfpDfLSM/UkeRyPIkvzGzEuA+ghnmhxD0eVwIXJ6e8kSky3KHHRvjj0ja0wG9BppiVkzKK9obAjUnxNxKigREUc9M/ESdSthJgr+NLMn+feD8yGEDGoDr3P2G9JQnIp3Wzm2RfoVVLfsadj9v3LbvOTl5wbDU0v4wYFTLW0k9B0D3sg7ZAd3RhB2qWwr8J/ArYCzBPI96YK67b0hfeSLSIe1qCFoFcYesRo7t2BhzkkHJwUEI9PkMHPpPLYeslhwEObmZ+Ikkxn7Dw8zyCOZ1/Iu7PwY8lfaqRCR7NTcFE9laBENUQGyLs2tCt7K9k9sqx0ZCYeDeW0k9+kJeQfv/PHJAwuznscvMPgaa2qEeEckkd9i2Lub20ap9byVt+RA85tdBQcne1sEhR8V0QA8MOqALumfmZ5K0CNthfi9Bx7hmgIt0VO6wY1NM/0LsKKU10NSw73m5hXsnulUf30oHdKn6GbqYsOGxAjjbzOYDf6Hlfh64+8zUliYiCdm5veXto306o1fDzphVhCw3uF1U2h/6HwtHfKVlB3RxhYJBWggbHjdFvvYHRsR53QGFh0i67NoJW9bEmcsQNYz10zhjV4oPCkKgYggM/kLMXIb+0OMQdUDLAQkbHtVprUKkK2tugq0ft71u0ta1tFxptdfeEBg4qmUHdM9+Kd8ASGS3MKOtCoDvAfe7+/y0VyTSmbgHey/Erq4aHRBbPoTmXfuel1+8NwQOPjJ+B3QWrrQqXUeY0VY7zWwq8Gg71CPSsURv9dlaQOzase85Ofl7+xMGjY/fAd1BV1qVriPsbat/AEcBs9NYi0h22bPVZytzGWK3+oSolVYHQN+j4fAvxemA7rwrrUrXETY8fgA8YGYrgSfc3fd3gkhWi97qs7XJbvG2+ty90uqerT5jO6D7Qm627O4skj5h/yt/GCglGKa7y8xie+/c3QelujiRA9Jiq884k922fhxnpdXSvWGwe6vP6FtJPfu3+1afItkqbHg8R4uhHiIZEG+rz33WTFoVf6vPvG57Q2DwF1reSirtD4U9MvMziXRAYVfVPT/NdSTEzI4AvgtUAM+5+80ZLklSpWFr/CUxolsQjdv3PScnH3r2DUJg4JiW23yWDlAHtEiKtfvNWTObSbCx1Fp3Hxp1/GTgBiAXuM3dr2vtM9z9bWCameUAt6a5ZEmVXQ1RM5/r4u/PsCNmq08smMi2e8jqkH+OBENUi6H4IHVAi7Sz0OFhZsOBHwMTgF7AaHdfZGa/AGa7+19DftSdwI3A3VGfnUswi/1EoA6Yb2azCILk2pjzL3D3tWZ2GnBl5LMk05p2wdaPWtn/OfJ12yctz+teHgRD70HBsNXdq65qq0+RrBZ2P4/jgGeBWuB+4LKol5uBaUCo8HD32WZWFXN4NLDM3Wsj13sQON3dryVopcT7nFnALDN7IlJTvLqnAFMAKisrw5Qn8TQ3B3s8t5jHENVy2PJRy5VWC3vuDYG+w/a9jaStPkU6tLAtj+uAp4GvErQGosNjEXBeknX0B1ZFPa8DxrT2ZjM7AfgaUEgbK/26+wxgBsDIkSPV4R/Pnq0+21gzqc2tPvtD9cSWQ1ZLIyutikinFDY8jgW+5u5uZrG/hNcBfZKsI15PZqu/7N39ReDFJK/ZNeze6rO1NZPibfVpuXsDoP9IODJmyGrpgOB2kzqgRbqssOGxA2htJ5e+QGwvZ6LqgIFRzwcAa5L8zM5vn60+WwmIuFt9HhQEwO6tPmOHrJYcrJVWRaRNYcPjb8D3zOwvUcd2twwmA88nWcd8YIiZVQOrgUnA2Ul+ZscWu9VnvL6GuFt99o6EQPRWn1G3knr001afIpK0sOHxY+DvwOvAIwTB8W0zu55gf49RYS9oZg8AJwAVZlYHXO3ut5vZZQT9KrnATHdfEvqn6GjibvUZ09fQ2lafu0OgxVafkQ7oguLM/Ewi0qVY2GWqzOxY4FcEQ3VzCUZZvQxc7u7/SFuFKTJy5EhfsGBB+i8Ud6vP2L6G/Wz1uWdEkrb6FJHMMrOF7j4y9njoeR7uvgj4opkVAWXARnffvp/TOp94W33GBkRbW332Gw5HnBoJhaiA0FafItKBJDzD3N130NE7s92D5bYbNkPDlmBPhoZNUd9viXptU9Qtpv1s9Vl+KNScsPfW0u7JbiUHa6VVEelUus5vtE+Wwu+OigTDlpY7t8WTXwxFPffOgh4wKqYDeoC2+hSRLqnrhEduAVSOC2Y9F/YIQqGwR7AM957vo14r6KHWgohIK7rOb8eyGvjajExXISLSKWgpUhERSZjCQ0REEqbwEBGRhCk8REQkYQoPERFJmMJDREQSpvAQEZGEhV4YsaMzs0+AlZmuI4RSkt8fpb1lsub2uHaqr5GKzzvQzziQ8xI9p4JgkzgJJ9v/zQ9y9xYb/nWZ8OgozGyGu0/JdB2JyGTN7XHtVF8jFZ93oJ9xIOcleo6ZLYi3CqvE1xH/zYNuW2WjxzJdwAHIZM3tce1UXyMVn3egn3Eg53XE/yY7kg75v69aHiKSUmp5dA1qeYhIqmkRuS5ALQ8REUmYWh4iIpIwhYeIiCRM4SEiIglTeIhIuzGzGjO73cweyXQtkhyFh4iEYmYzzWytmS2OOX6ymb1jZsvM7Mq2PsPda919cnorlfbQdbahFZFk3QncCNy9+4CZ5QI3AScCdcB8M5sF5ALXxpx/gbuvbZ9SJd0UHiISirvPNrOqmMOjgWXuXgtgZg8Cp7v7tcCp7VyitCPdthKRZPQHVkU9r4sci8vMys1sOjDczK5Kd3GSPmp5iEgyLM6xVmceu/t6YFr6ypH2opaHiCSjDhgY9XwAsCZDtUg7UniISDLmA0PMrNrMCoBJwKwM1yTtQOEhIqGY2QPAHOBwM6szs8nuvgu4DHgaeBt4yN2XZLJOaR9aGFFERBKmloeIiCRM4SEiIglTeIiISMIUHiIikjCFh4iIJEzhISIiCVN4SJdlZseY2TVmVpbpWuKJ1NZuY+nN7E4zq2uv60nHpvCQruwY4GogK8MDuA0Yl+kiROLRwogiWcrd6wjWjhLJOmp5SKdlZoeZ2aOR3e92mNkHZvawmeWZ2fnAHZG3vmdmHnlURc7NM7OrzGypmTWY2Roz+42ZFUV9flXknEvM7PrIdbab2eNx9r2IV99JZvaKmW0ys62R3fj+X9Tr+9y2MrMXo+qMfVRFvW+imT1nZlvMbJuZPW1mQxP43224mb0c+VneMzOtgistKDykM3ucYG+Ji4GTgCuBBoL/7p8AfhZ535kEt4fGAR9Gjt0L/Ai4HziFYFe8ycB9ca5zFTAE+A5wKTAC+F8zy2+tMDOrIVhAcDnwTeA04HqguI2f55KoOscBxwHvAh8D9ZHPPQV4DtgKfAs4G+gBvGxmA+N8ZqyekZ/5XuB0goUPbzazz4c4V7oSd9dDj073ACoI9pU4rY33nB95z6Exx4+PHD8v5vg5kePHRJ5XRZ6/BeREve9zkeOT27j2GZH39GzjPdcE/0Rbff1G4FNgTNSxZcBzMe/rCawDfref/83ujNT0+ahjhZFzZ2T6/1M9suuhlod0VuuBWuA6M7vIzIYkcO7JwE7gj5HbV3lmlgf8b+T1CTHvf8Tdm3c/cfe/E/RVtNXZ/RrQCDxoZmeY2UEJ1IeZXUrQEjnP3edFjg0BBgP3xdS9nWA13Ni649nu7i9E/SwNwHtAZSL1Seen8JBOyd0dOBFYQHDL6V0zqzWzi0OcfhBQQHDrpzHqsTbyennM+z+O8xkf08Z2rO6+jOBWWg5wD/CRmc0zs4n7K87M/hm4AfiRuz8cUzfA7TF1NxLsJx5bdzwb4hxrAIriHJcuTKOtpNNy91rgPDMzYBjBvhP/Y2Yr3P2pNk5dD+wguH0VT+xOeQfHec/BBK2Ltup7AXjBzAoJbnX9J/CEmVW5+7p455jZEcBDwL3u/os4dUPQB/NsnNN3tlWPSCIUHtLpRVohr5nZ5QSd3kOBpwj+ogboFnPKX4ErgFJ3fy7EJc4ws2t237oys88RbMc6J2R9DcDzZlYC/AWoJuhn2IeZlRMMAngdmBLno94BVgCfdffrwlxb5EApPKRTMrOjCW7t/IGgEzmXoIN8F/B85G1vRb5eamZ3EdzeecPdX4zsmveImV0PvAo0E3SQfxm4wt3fjbpcD+DPZnYL0IfgNtl7wN1t1DeNoA/iSWAVQQf/VQStmsWtnHZf5H3/ChwbNKj2+Ie7N0T6Qv4S2RL2IYIQOhgYD3zg7te3VpNIIhQe0ll9BHwAXE7QCtgBvAmc6u4LAdz9dTO7huCv+IsI+h+qCf56/xbBL+kLgB8StFJWEGy3GtvHcS1wKMFopWLgBeAyd29so77XgS9Fzj2IYKjt34Bz3P3TVs75DMHIqSfivFYNrHD3J81sQqTm2whaVR8BcwmCVCQltA2tyAGKTMxbDlzk7rdluByRdqXRViIikjCFh4iIJEy3rUREJGFqeYiISMIUHiIikjCFh4iIJEzhISIiCVN4iIhIwhQeIiKSsP8PQ22RqOLrcUgAAAAASUVORK5CYII=\n",
+ "text/plain": [
+ "<Figure size 432x288 with 1 Axes>"
+ ]
+ },
+ "metadata": {
+ "needs_background": "light"
+ },
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "plt.loglog(grid,error_inf)\n",
+ "plt.loglog(grid, grid**2, label='2nd oder')\n",
+ "#plt.loglog(grid,grid, label='1st oder')\n",
+ "plt.ylabel('error in maximum-norm', fontsize=16)\n",
+ "plt.xlabel('step size h', fontsize=16)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 18,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "Text(0.5, 0, 'step size h')"
+ ]
+ },
+ "execution_count": 18,
+ "metadata": {},
+ "output_type": "execute_result"
+ },
+ {
+ "data": {
+ "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY8AAAEPCAYAAAC6Kkg/AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAAyFElEQVR4nO3deXxU9b3/8dcn+0JISMIaCARFK8WNTUVFbeutXq3eqq1irVVR1Nbuva12o7fL1dv+auutVkXFpVXUulRcrrauqFhlcaliUYwsAQUEBAIkhPD5/XEmMEwmyRkyk5kk7+fjMQ8y58yc8wHNfOb7+W7m7oiIiCQiK90BiIhI96PkISIiCVPyEBGRhCl5iIhIwpQ8REQkYUoeIiKSsJx0B9BVKisrfcSIEekOQ0SkW1mwYMFH7t4/9nivSR4jRoxg/vz56Q5DRKRbMbNl8Y6rbCUiIglT8hARkYQpeYiISMKUPEREJGHdInmY2Ugzu8XM7mvvmIiIdI20JQ8zm2lma8zszZjjJ5jZYjNbYmaXA7h7rbtPjX5dvGMiItI10jlU9zbgWuCOlgNmlg1cBxwP1AHzzGy2uy/qioCadzr1jTuCR8MO6hub2Nywg30H9GFov6KuCEFEpFtIW/Jw9zlmNiLm8ERgibvXApjZ3cCpwF4lDzObBkwDKK+qYfpDb7J5V2LYnSQ2N+5gS+MOtm5vjnud4rxs/njOOI7Zr9U8GRGRXinTJglWASuintcBh5lZBfAr4FAzu8Ldr4x3LPZi7j4DmAGQP3iU//W1VfTJz6GkIIc++TmUF+dRXV5ESUEOxXk59IkcD87n0qcgh7zsLH7xyCIuuG0eV512IF8YPyz1/woiIhku05KHxTnm7r4OuCTmYKtj7TmwqpT50/9tr4K65+LD+eqdC/nP+97gg40NfP1T+2IWL1QRkd4h00Zb1QHRX+2HAqvSFMsuJQW53PKVCZw2toqr//4OP3zwn+xo3pnusERE0ibTWh7zgFFmVgOsBM4Czk5vSIG8nCx++4WDGVxawHXPvMfqTY1ce/ahFOVl2j+hiEjqpXOo7izgJWB/M6szs6nuvgO4DHgCeBu4193fSleMscyM//zsJ/jlf4zh2cVrmDLjH3xU35jusEREupy5e/gXB4X+wUBB7LmWEVKZavz48Z7MVXX/vmg1X5+1kIF9C7j9/ImMqCxO2rVFRDKFmS1w9/Gxx0O1PMysIjJstoFgNNS7cR69yvGjB3LXRYezaVsTp18/l9dWfJzukEREukzYgv0twHEEk/r+BWxPWUTdyNjqftx/6SS+cusrnDXjJa47eyyfPmBgusMSEUm5UGUrM9sIfNPdb0t5RCmS7LJVtLWbG5l6+zzeXLmRX/7HgZx9WHVK7iMi0tU6VbYC1gOrkxtSz9G/JJ9ZFx3OMfv154cP/pPf/m0xifQliYh0N2GTxx+AS0wz49pUnJ/DTeeO58zxw/jD00v43l/eoElzQUSkhwrV5+HuV5vZEGCRmT0JbGj9Ep+e9Oi6mZzsLK46/UAGlxXw+yffZW19I3/80lj65GsuiIj0LGH7PP4duB/Ib+Ml7u7ZyQws2VLZ5xHPPfOW88MH3+SAwSXMPG8CA0pajW4WEcl4ne3zuJpg9vfBQL67Z8U8MjpxpMOZE6q5+dzxvLdmC6f9cS7vra1Pd0giIkkTNnlUA79093+6e1MqA+pJjvvEAO65+HAampo5/fq5zF+6Pt0hiYgkRdjk8SowJJWB9FQHDS3jgUuPpF9RHl+6+WUef/PDdIckItJpYZPHN4DvmdmRqQymp6quKOL+SycxekhfLr1zAbfPXZrukEREOiVs8vgrwVLpc8xsk5ktj3ksS12IPUN5cR53XXg4nzlgINNnv8WV//c2O3dqLoiIdE9hx5A+BeiTrpMK87K54ZxxTJ/9Jjc+V8uHGxv4zRkHk5eTaduqiIi0L2zy+CbQ6O4NqQymN8jOMn5x6hgGlxbymycWs3ZzIzd8eRx9C3LTHZqISGgdfuU1sxxgHXB86sPpHcyMrx23L1d/8WBeeX89X7zhJT7cqLwsIt1Hh8kjskHTaqA59eH0LqeNHcqt50+gbsM2Pv/HF3ln9eZ0hyQiEkrYYvufgQtTGUhvdfSo/txz8eE073ROv34u/6hdl+6QREQ6FDZ5LAUmmNk8M/uxmU01swuiHymMsRUzG21m95rZ9WZ2RlfeOxU+OaSUB746iYF9Czj3lld4+PVV6Q5JRKRdYde26mh52E6vbWVmM4GTgTXuPibq+AnANUA2cLO7X2Vm3wVecffnzWy2u5/S0fW7em2rvfHx1u1Mu2MBryxdz49POoALjx6Z7pBEpJdra22rsKOtapIcTzy3EexUeEfLATPLBq4j6KyvA+aZ2WzgT8B0MzsFqOiC2LpEWVEed0ydyHfufY1fPvo2qz5u4McnHUBWllbCF5HMEnZJ9pRPAnT3OWY2IubwRGCJu9cCRPZRP9XdrwS+FkkuD7R1TTObBkwDqK7uHrv7FeRm84cpYxlQsoiZL77P6k0N/PaLB1OQq7UnRSRzJLTRhJmNAY4BygmG785x9zdTEVhEFbAi6nkdcFgkyfwQKAZ+09ab3X0GMAOCslXqwkyu7Cxj+udGU1VWyK8ee5u19Y3c9OXxlBZpLoiIZIZQySMy1+M2YAoQXUNxM7sLOM/dUzGUN169xt19KZEWRU9lZlw0eSQDSwv43r2vc/oNc7n9golUlRWmOzQRkdCjraYDXwR+StD/URj586fAmZE/U6GOYE2tFkOBXjUU6ZSDh3D7BRNZvamBz1/3IotWbUp3SCIioZPHOcAv3P1X7r7M3Rsjf/4K+CVwborimweMMrMaM8sDzgJmp+heGeuIfSq475JJZGcZX7zxJV5496N0hyQivVzY5DEEeKmNc3NJwl4fZjYrco/9zazOzKZGZrdfBjwBvA3c6+5vdfZe3dH+g0p44KuTqCor5LxbX+HBV+vSHZKI9GJhO8xXAUcCT8Y5N4kklJLcfUobxx8DHuvs9XuCwaWF3HvJEVz8p/l8+57X+WBjA5cesw9mGsorIl0rbMvjTuBHZvYTMxtpZoWRUtIVwI8I5l1IFygtzOX2CybyuYOH8OvHF/PTh96iWfuCiEgXC9vy+BkwEvivyM8tDJgVOS5dJD8nm2vOPIQhpQXcOKeW1Zsa+N8ph2ouiIh0mVDLk+x6sdkngckE8zzWA8+5+6IUxZZU3WF5kr1x24vv81+PLOLQYWXc/JUJlBfnpTskEelBOrs8CQCRzupe2WGdqc47soaBfQv45j2vccb1c7nt/IlUVxSlOywR6eESnWE+CKgGCmLPufucZAUliTnxwMH0L8ln6u3zOe36F7n1vIkcOLQ03WGJSA8WqsPczKrM7BlgJcFw2mcjj2ei/pQ0Gj+inPsvPYL8nGzOnPESzy5ek+6QRKQHCzva6npgDPB94ETguMjjU1F/SprtO6CEB786iREVxUy9fT73zl/R8ZtERPZC2LLV0cA33F1DcjPcgL4F3HPx4Xz1zoV8/743+ODjBr7x6X01F0REkipsy2MboDpIN1FSkMvM8yZw2tgqfvfkO/zwwX+yo7mj/bxERMIL2/K4CfgywTIh0g3kZmfx2y8czJDSQq59ZgmrNzVy7dmHUpSX0BgJEZG4wn6SrAS+bGZPEywVsj72Be4+M5mBSeeZGd/77P4MLivgJ399k7Nm/IOZ502gsk9+ukMTkW4uY/YwT7WeOkkwrL8vWs3XZy1kQEkBt18wkZrK4nSHJCLdQFuTBMP2edR08BiZpDglRY4fPZC7Ljqc+sYdnH79XF5dviHdIYlINxYqeUT27mj3kepApfPGVvfj/ksn0Sc/hyk3/YMnF61Od0gi0k2FbXlID1FTWcz9l05iv4ElTPvTfO58WXlfRBKn5NEL9S/J5+5ph3PMfv350YNv8v+eWEwiC2SKiCh59FJFeTncdO54zpowjGufWcL3/vIGTZoLIiIhdctB/2Z2NPAlgvhHu/ukNIfULeVkZ3HlaQcyuLSQ3z35Dms2N3D9OePok98t/7cQkS6UMS0PM5tpZmvM7M2Y4yeY2WIzW2JmlwO4+/PufgnwCHB7OuLtKcyMb35mFL8+/SDmvreOM298iTWbGtIdlohkuIxJHsBtwAnRB8wsG7iOYDHG0cAUMxsd9ZKzCXYylE764oRh3PyV8bz/0RY+/8e5LFlTn+6QRCSDhU4eZvYVM3vczBaZWW3M473OBhLZDyR25vpEYIm717r7duBu4NRIPNXARnff1E7M08xsvpnNX7t2bWdD7PGO238Ad087nMYdzZxxw1zmL221kICICBB+P4+fALcCQ4DXgOdiHqnaCKoKiF5XvC5yDGBqJKY2ufsMdx/v7uP79++fohB7loOGlvHApUfSryiPs29+mcff/CDdIYlIBgrbMzoVuMbdv53KYOKIt464A7j79C6Opdeoriji/ksnceHt87j0zoVMP3k05x1Zk+6wRCSDhC1bVQAPpzKQNtQBw6KeDwVWpSGOXqe8OI87LzyczxwwkJ89vIgrH3ubnTs1F0REAmGTx3PAwakMpA3zgFFmVmNmecBZwOw0xNErFeZlc8M54zjn8GpunFPLt+55jcYdzekOS0QyQNiy1beAB8xsHW0vyd6pGWZmNgs4Fqg0szpgurvfYmaXEewjkg3MdPe3OnMfSUx2lvGLU8cwpKyQXz++mLWbG7nx3HH0LchNd2gikkaJLsne1ovd3TN6ZllvX5I9GR5YWMf373uDfQf04dbzJzC4tDDdIYlIirW1JHvYD/yf03bikF7itLFDGVBSwCV/XsBpf5zLbedPZP9BJekOS0TSIFTLoydQyyN53lq1kfNvnce2pmZmfHk8R+xTke6QRCRFOrsZlMgunxxSyoNfO5KBfQv4ysxXmP26BsCJ9Dah+ykio51OBPYHCmJOu7v/IpmBSWarKivkvkuOYNodC/jGrFdZvbGBC4+uwSze1BwR6WlCJQ8zGwK8AIwg6Pto+YSIrnkpefQyZUV53DF1It+59zV+9djbrNq4jR+fNJrsLCUQkZ4ubNnqN8BaoJogcRxGsG/5r4AlaA/zXqsgN5trp4zl/CNHcOuLS/n6rIU0NGkuiEhPF7ZsdTTwPXbP7t7p7kuBn0ZWvv1fIgsWSu+TlWVM/9wnqSor5JePvs2ydXP5xqdHcfwBA8lSK0SkR0pkeZJVkYmAW4B+UeeeJpjcJ73chUeP5IZzxrKpoYmL/7SAz1z9HLNeWa6WiEgPFDZ51AGVkZ/fA/4t6txEQLsHCQAnjBnMM989lj9MOZTi/ByueOCfHPU/z3DdM0vYuLUp3eGJSJKEnWF+A7DN3b9tZpcSbND0JNAEfBa40d2/ltJIO0nzPLqeu/PSe+u4YU4tc95ZS1FeNlMmVnPBUTVUlWl2ukh30NY8j7DJoxIod/d3Is+/DpwJFAGPAz9394xufSh5pNeiVZu46flaHo7MCfncwUOYNnkkBwzum+bIRKQ9nUoePYGSR2ZY+fE2Zr7wPrNeWc7W7c1M3q8/l0weyRH7VGiOiEgGSkryMLMsgr3EK4D57r4leSGmlpJHZtm4tYk/v7yMW19cykf1jRxYVcq0ySM5ccwgcrK18IFIpuh08jCzrwHTCRIHwAR3X2hmfwWedvf/TVawqaDkkZkampp58NWV3DSnltqPtjCsvJCLjh7JF8YNozAvO93hifR6nVrbyswuAq4B/krQ1xFdX3geOD0JMUovVJAbdKL//TvHcMM546jsk89PH3qLSVc9xe/+/g7r6hvTHaKIxBG2w/xtYLa7/yAyKbAJGB9peZwE3OLug1Ica6eo5dE9uDvzl23gxufe48m311CQm8UXxw/jwqNGUl1RlO7wRHqdzu7nUUOwm188W4CyvYxLZA9mxoQR5UwYUc6SNZuZMaeWWa8s58//WMaJBw7m4skjOWhoWbrDFOn1wvZMfkSwKGI8+wMrkxJNSGZ2rJk9b2Y3mNmxXXlv6Tr7Dijh12cczAs/+BTTJu/DnMVrOeXaF5ky4x88u3gNvWWkoEgmCps8HiZYxyp6AUSPzP/4NkFfSKeY2UwzW2Nmb8YcP8HMFpvZEjO7vOXeQD3B0vB1nb23ZLaBfQu4/MRPMPeKT/Gjfz+A9z/awnm3zuPEa57ngYV1NDXv7PgiIpJUYfs8KoC5wDDgZWBy5PkngDXAJHff2KlAzCYTJIQ73H1M5Fg28A5wPEGSmAdMAf7l7jvNbCBwtbt/qaPrq8+j59i+YyezX1/FjDnv8c7qeoaUFnDBUTWcNbGaPvmht6gRkRA6NdrK3dcB44ErgVyC9a1ygGuBIzqbOCL3mAOsjzk8EVji7rXuvh24Gzg1skAjwAYgv7P3lu4lLyeLM8YN5fFvTmbmeeMZVl7ELx99m0lXPsWvH/8XazZn9GIHIj1CRs0wN7MRwCNRLY8zgBPc/cLI8y8T7CXyNMGaWmXA9e7+bBvXmwZMA6iurh63bNmyFP8NJF1eXb6BGXNqefytD8nNyuK0sVVcNHkk+/Tvk+7QRLq1zo62Spd461W4uz8APNDRm919BjADgrJVkmOTDHJodT+uP2cc73+0hZufr+UvC+q4Z/4Kjj9gIBcfM5Jxw8vTHaJIj9Jm8jCzpxO4jrv7p5MQT6w6gn6WFkPZvSGVSCs1lcX86vMH8u3j9+OOuUu5/aVl/G3RasYP78fFx+zDpz8xQBtUiSRBe30eWQTf/FsenyDY9GkEUBj581iCobqp+m2cB4wysxozywPOAman6F7Sg1T2yec7/7Y/cy//FNM/N5oPNjZw0R3zOf53z3HPvOU07tAGVSKd0WbycPdj3f04dz+OYGmSJuBwdx/p7ke4+0jgiMjxazobiJnNAl4C9jezOjOb6u47gMsIJii+Ddzr7m919l7SexTn53D+kTU895/Hcs1Zh5Cfk80P7g82qLr+2ffYuE0bVInsjbBDdf8J/Mbd74hz7jzgu+5+YPLDSx4N1RUIlj95YclHzJhTy/PvfkSf/BymTBzGBUfVMLhUG1SJxOpsh/koYG0b59YA++5tYCJdycw4elR/jh7VnzdXbmTGnFpmvriUW19cyimHDOHiyfuw/6CSdIcpkvESWRhxsbv/R5xzDwH7ufsByQ8vedTykLasWL+VW154n3vmrWBbUzPH7d+faZP34fCR5dqgSnq9zm5DexZwJ0G/w33AamAgcAZBR/qX3P2epEacZEoe0pENW7bz538s47a5S1m3ZTsHDy3l4mP24bOfHES2RmhJL5WMzaA+A/wXwUzzXIKO8nnAdHd/KomxpoSSh4TV0NTMfQvquPn5Wpau28rwiiIuPHokXxg3lIJcbVAlvUvS9jCPbEVbCXwUtUxIxlPykEQ173T+9taH3DCnltdXfExFcR7nHjGCc48YTr/ivHSHJ9IlkpY8uislD9lb7s4r76/nxjm1PP2vNRTmZnPmhGFMPaqGYeXaoEp6toRHW5nZT4Gb3X1V5Of2uLv/orNBimQiM+OwkRUcNrKCxR8GG1Td+fIy7nhpKScdNISLJ49kTFVpusMU6VJttjzMbCfBpMBXIj+3x909o4vBanlIMn2wcRu3vriUu15eTn3jDo7ct4KLJ+/D0aMqNUJLehSVrZQ8JAU2NTRx18vLmfnC+6zZ3MgBg/ty8eSRnHTQYHKzw+61JpK5lDyUPCSFGnc089Crq5jxfC1L1tRTVVbI1KNqOHPCMIq1QZV0Y0oeSh7SBXbudJ7+1xpunPMe85ZuoLQwly8fPpyvTBpB/xLtWybdT8LJI9LPETazuLtn9NcrJQ/paguWbWDGnPf426LV5GYHux9edPRIaiqL0x2aSGh7s7bVzwmfPEQkxrjh/bjxy+OpXVvPTc+/z30L6pj1ynI+O3oQFx8zkkOr+6U7RJG9prKVSBdZs7mB2+cu5U8vLWNTww4mjijn4mNGctz+2qBKMldn17Y6zt2faef899z9/3UyxpRS8pBMUd+4g3vmreCW52tZtbGBYeWFHFZTwdjqfowb3o9RA/oomUjG6Gzy+Bg4xt1fj3Puu8D/qM9DJDFNzTt55I1VPPrGh7y6fAPrtmwHoKQgh0Or+zGuuh9jh5dxyLAySgpy0xyt9Fad3c/jXuBxM5vk7u9HXfRbwK+BryclSpFeJDc7i88fOpTPHzoUd2fZuq0sWLaBBcs3sHDZBn7/1Du4gxnsP7CEccODlsnY6n4MryjSZERJq7AtjyzgAeCTwCR3X2tm3wB+D3zT3f+Q0ihbx3MA8E2CBRqfcvfrO3qPWh7S3WxqaOL1FR8HCWXZBl5b/jGbG3cAUFGcx9hIMhk3vB8HVpVqxV9JiWQsyV4A/B0oAu4G/gf4jrv/PkkBzgROBta4+5io4ycQ7JGeTbDW1lVR57KAm9x9akfXV/KQ7q55p7NkTf2uZLJw+Qbe/2gLALnZxughpYyL9JuMHV6mbXUlKZIySdDMyoDngdHA9939t0kMcDJQD9zRkjzMLBt4BzgeqCPYP2SKuy8ys1OAy4Fr3f2ujq6v5CE90br6Rl5d/jELlgcJ5fUVH9O4I1iKbkhpwR6tkwMG99WSKZKwvVlV9442Tq0h2EXw4KjXuLt/pTMBuvscMxsRc3gisMTdayMx3Q2cCixy99nAbDN7FOgweYj0RBV98vnM6IF8ZvRAIOiEf/uDTbtbJ8s28MgbHwBQkJvFQUPLgmRS3Y+xw/tRrn1JZC+112E+mbYnCW4Bjo56nqrJIlXAiqjndcBhZnYscBqQDzzW1pvNbBowDaC6ujpFIYpkjtzsIEEcNLSM84+sAYIVgBcu+3hXZ/xNc2q5fmfwK1tTWbxriLCGCUsi2kwe7j6iC+NoS7z/i93dnwWe7ejN7j4DmAFB2SqpkYl0E4NLCznpoEJOOmgwEGyz+0bdRhZGSl3PLl7D/QvrgGCY8CHDynYlEw0TlrZk9NwMgpbGsKjnQ4FVaYpFpEcoyM1mYk05E2vKAXYNE25JJguWbeCap97dY5jw2Eipa9xwDROWQKjkYWYd1nzcfXnnw2llHjDKzGqAlcBZwNkpuI9Ir2VmjKgsZkRlMaeNHQrA5oYmXl+xcVep6+HXVnHXy8GveMsw4ZZy10FDNUy4Nwrb8lhKx/0anfq/x8xmAccClWZWB0x391vM7DLgicj1Z7r7W525j4h0rKQgl6NGVXLUqEogWGr+3TX1u1onC5dt4O+LVgOQk2V8sqqUsdW7y10aJtzzhZ0keB6tk0cFcBIwEviFu89MenRJpKG6Ism1fst2FkbmmyxYtoHX6z6moWnPYcItrZPRQzRMuLtK2WZQZvYnYJm7/7hTF0oxJQ+R1IoeJrxw+ccsXLaBlR9vA3YPE25JJmOry6joo82xuoNUJo/PAre6+5BOXSjFlDxEul70MOGFyzfw1qqNNDXvOUx47PCyyDDhErI1TDjjdHZhxPYMAAqScB0R6WHiDRP+58qNu/pNnnsnaphwfg6HRPpNxlb345DqMvpqmHDGCjvaanKcw3nAGOAKgiVLRETaVZCbzYQR5UwYsXuY8PL1W6PW6/qY/33qXXbGDBNuKXeN0DDhjBG2wzzefuYt/wWfA77k7hk9/0JlK5HuIXaY8KvLN7C5IVhNuLw4j4OHljKispjh5UVUVxRRXV7MsPJC8nM0XDgVOlu2Oi7OsQaCjvIPOxWZiEiUeMOEl6zdvZrwmys38vL769m6vXnXe8xgcN8ChpUXMbyiiOEVxcHPkedlRVrDK9m0h7mIdDvuzkf121m+fivL129h2bqtwc/rtrJs/VbWbm7c4/V9C3KorihieHlxpLVStKvlMri0UB317dibVXUPAt5x94YOLtwPOM7dH+h8mCIiHTMz+pfk078kn3HD+7U6v3X7Dlas38aydVsiCWYry9ZtZdEHm/jbog93jfiCYC+Uof0iCSWSWIKfg3JYUV6mr+KUHu39q7wKHAG8Ars2XvoYODpmL/P9gL/QyRnmIiLJUpSXw/6DSth/UEmrc807nQ82btvVStndYtnCwqj+lRb9S/Kj+ldaEkwx1eVFVPbJ67Ud+O0lj9h/EQP6oCQhIt1YdlbQ0hjar4hJMefcnY3bmlgWSSwr1m9l2bqgLPbSe+t48NWVRFf6i/Oyd/WzVJcXUV0R6cgvL6KqX2GPnlWv9piISISZUVaUR1lRHgcPK2t1vqGpmboN21i+fsvulsu6rdSu3cKzi9fu2sURgiQ1pKyA4eXFuzvyy4t2/dzdl7pX8hARCakgN5t9B/Rh3wF9Wp3budNZs7mxVT/L8vVbeeKtD1m/Zfsery8vzovqX9mzr2VASX7Gb8ql5CEikgRZWcag0gIGlRZw2MiKVuc3NwTlsBXrgxZLy8+vrtjAo//8gOadu+th+TlZu5JJdUXLkOOgBZMpc1o6Sh6fM7MxkZ+zCCYKnmJmh0S9ZmQqAhMR6UlKCnIZU1XKmKrSVueamneycsO2oLWyfivLI62XZeu28lLtug7ntES3YLpqTkub8zwis8rDcndPfypsh+Z5iEh3lIw5LSeOGbTXSWVvZpjX7NWdREQkaZIxp+XIfSqT3iJpM3m4+7Kk3klERJIuzJyWVOzs2C0HIZvZSDO7xczuS3csIiKZqmVOSyqWX8mY5GFmM81sjZm9GXP8BDNbbGZLzOxyAHevdfep6YlUREQyJnkAtwEnRB8ws2zgOuBEYDQwxcxGd31oIiISLWOSh7vPAdbHHJ4ILIm0NLYDdwOndnlwIiKyh4xJHm2oAlZEPa8DqsyswsxuAA41syvaerOZTTOz+WY2f+3atamOVUSk18j0Gebxennc3dcBl3T0ZnefAcyAYJ5HkmMTEem1QicPMzsGmAJUAwUxp93dP53MwCLqgGFRz4cCGb3drYhIbxCqbGVmFwPPAKcDZQQtguhHqspf84BRZlZjZnnAWcDsFN1LRERCCtvy+C5wF3BBpOM66cxsFnAsUGlmdcB0d7/FzC4DniDYR2Smu7+VivuLiEh4YZNHFXBrqhIHgLtPaeP4Y8BjqbqviIgkLmy5aQFaPVdERCLCJo9vAN8ys8mpDEZERLqHsGWrh4G+wDNmthXYEHPe3X14UiMTEZGMFTZ5PEWwEZSIiEi45OHu56U4DhER6UYyfXkSERHJQG22PMzsXOBRd18X+bld7n5HUiMTEZGM1V7Z6jbgcGBd5Of2OKDkISLSS3S0h/kHUT+LiIgAIfcw137mIiISTR3mIiKSMCUPERFJmJKHiIgkTMlDREQS1uEMczPLBsYAq9xdG4GLiGSYpuYmtjRtob6pni1NW1r9fGLNiRTnFif1nmGWJ3FgPnAS8Lek3l1EpJfa6TvZ2rR114f8Hh/82+vZumMr9dvjnGuq3/N92+vZvrP9rZbGDRxHTWlyZ1x0mDzcfaeZrQCSm7ZERLoZd6ehuWHPb/fb9/ymH/vhHv266HNbd2wNdc/87HyKc4spzi2mT24finOLGVA0gJrcml3Pi3OL6ZPXh6KcIvrk9dnjtcW5xVQWVib93yLsqro3Euzn8WgqdxMMy8xGAj8CSt39jHTHIyKZrWln064P7za/wcd8849X/tnStIVmb+7wftmWvesDvCi3iD65fSgtKKWqpKpVIoh9XUsiKM4JzuVm53bBv1DiwiaPEmAfoNbMHieYeR69RLu7+/QwFzKzmcDJwBp3HxN1/ATgGoK9ym9296vauoa71wJTzey+kPGLSDcTW9bZ49t9TFmno2/+jc2Noe5ZnFscfGjn7f5AryisaPfDvuWbfvS5/Ox8zCzF/0LpFTZ5/DDq5wvinHcgVPIgWCfrWqLWwop0yl8HHA/UAfPMbDZBIrky5v0XuPuakPcSkS6UzrJO9If7wKKBjMwbuUciiPvBH3W8KLeILNMA1LDC7ueRtH9Rd59jZiNiDk8ElkRaFJjZ3cCp7n4lQStFRFKoo9E6Lc/b+8afaFlnV2km8uHeUtZpVcJpo67f8rrcrMws6/R0YVseqVYFrIh6Xgcc1taLzawC+BVwqJldEUky8V43DZgGUF1dnbxoRTJA885mtu7Y2vpb/o49a/a7PvjjfOMPO1qnRbx6fWVh5R4f9rvOt/ONvzeUdXq6hJKHmZ0MHAOUEyzV/py7P5qEOOL9X9Tmtrfuvg64pKOLuvsMYAbA+PHjtY2upF10Wad+ez1bdmyJ25Eb++Ee73VhyzoF2QWtOmMHFQ3ao67f6oM/zjf+wpxClXW6k+1bYNMHsHkVDDsMcvKTevlQycPMSoBHgKOBHQSJowL4rpk9D5zs7vWdiKMOGBb1fCiwqhPXE0mqpuamNj/gO/rGH/vzTt/Z4f1yLKfVh3tZQdmusk6rETp5QUdv7DBNlXV6IHfYuj5ICpsij80fwKaVQbLYtCo417Bx93suWwCV+yY1jLAtj/8GxgJfBu529+ZIJ/dZwPWR89/oRBzzgFFmVgOsjFz37E5cT2SPsk573/JbSjnxvvG31PfDlHUMoyi3qNW39v6F/fccjRNVr4/7jT+vD3lZeSrr9EbNO6D+w0gSWLlnUoj+udXoMYM+A6HvYKjYB0YcBX2HBI+SwcHxJAubPE4Hfuzud7YccPdm4E4zqwS+T8jkYWazgGOBSjOrA6a7+y1mdhnwBMEIq5nu/lb4v4b0FO7Oth3b2h2GGXYW7rYd20Lds82yTlmc8fdxPuxV1pFQWspIu5JCnFZD/WpaVeyz84MP/75VUDUeDoj8XBL5s+/gIHF08XyQsMmjAljUxrlFkfOhuPuUNo4/BjwW9jqSWVrKOh0NwwwzQSuRsk70aJ1+Bf0YWjJ0jw/16M5blXUkJeKVkVpKR22VkVoUlEJJpIUw8JNRSaGl1TAEisohA1uhYZPH+wRDZv8e59y/R85LNxO2rBNmglaio3WiP8D7F/YPPVqn5WeN1pEu0dwUtAba61tot4w0pHUZqSUp9B0Med131adElif5rZn1Ae4kmGE+iKBv4kLgO6kJT2JlQlmnOLd4V1kn3iSseCN2VNaRjNOqjBSnbyFUGSm6byF9ZaSuFnaS4O/MrD/wbeC8yGEDGoGr3P2a1ITXc6RrtM7elnWiR/OorCPdSksZqVXfQktJKXKsMZEyUiQpZHAZqauFHapbCvwc+A1wOME8j/XAP9x9Q+rCS68wk7DSNVqnrWUXVNaRHi22jBSvbyFMGanm6NZ9C928jNTVwmwGlUMwr+Pz7v4w8H8pjyoFNm/fzMPvPdzmypmdLevEfrhHj9aJV87RaB2RGNu3dNC3sArq1xC6jBSdFHpBGamrhdnPY4eZrQY6XrAmgy3fvJwfvrB7fcd4k7D6FfRjWMmwVmWb6LH5sR/8KuuIdKCzZaSW0lG8MlLfKijspzJSGoTtMP8zQcd4tx1Ku0/pPjz6+UdV1hFJpqSXkaL6FlRGymhhk8dS4Gwzmwc8ROv9PHD3mckNLbkKcgqo7qvFEUVC61QZKVI2iltGGhIpI2XKuqyyN8L+17su8mcVMC7OeQcyOnmISETSykhjWg9RVRmp1wibPJK7c7qIpIbKSNJFwoy2ygO+Bdzl7vNSHpGIxKcykmSQMKOttpvZxcCDXRCPSO+jMpJ0Q2G/arwKHAjMSWEsIj1Pm2WkVbs36lEZSbqhsMnju8AsM1sGPOru2pVPpLE+Tkthb8pILUlBk9qk+wibPP4ClBIM091hZrG/Ee7uw5MdnEhauMPWde2voqoykvRyYZPHU7Szp7hIt9HcBJs/jEkKLX0NLT9/2LqMZFlBa6BkcPwyUt8qKBmkMpL0GmFX1T0vxXGIdN6uMtLKqP6EVXsmiHhlpJyC3YvkDZsYp29Bo5FEYnXL3wYzOwD4JlAJPOXu16c5JEml6DJS3A7nyM9xy0hlu/sXBh3Yeohq3yEqI4nshdDJw8wOBX4CTAbKgInuvtDM/huY4+6Ph7zOTIJdCde4+5io4ycA1xDsYX6zu1/V1jXc/W3gEjPLAm4K+3eQDNRSRtojKayi1R7PzTFL2u9RRtoXaia3TgolgyGvKD1/L5EeLux+HkcBTwK1wF3AZVGndwKXAKGSB3AbcC1wR9T1swmWQDkeqAPmmdlsgkRyZcz7L3D3NWZ2CnB55FqSiRrrYzqZY+cxdFRGqgrKSNGjkFo6o1VGEkmrsL99VwFPAP9B8IEenTwWAueGvaG7zzGzETGHJwJL3L0WwMzuBk519ysJWinxrjMbmG1mjxIkNOkqu8pIbXU4t4xG2tT6vQVlu/sTBh0Y0+E8WGUkkW4ibPIYC5zm7m5msaOuPgL6dzKOKmBF1PM64LC2XmxmxwKnAfm0s0y8mU0DpgFUV2tF3VCam3YngugyUvSw1c0ftl1G6jskUkY6pnVSUBlJpMcImzwagLZ+6wcDcXoqExLva2abQ4Pd/Vng2Y4u6u4zgBkA48eP11DjPcpIcfoWNq2CLWuJW0ZqKR0NO3zPUUgtj+IBKiOJ9CJhf9tfAL5lZg9FHWv5hJkKPN3JOOqAYVHPhwKrOnnN3qPNMlJMgohXRirst7s/YfBBrTucVUYSkTjCJo+fAC8CrwP3ESSOr5jZ1QT7e0zoZBzzgFFmVgOsBM4Czu7kNXuGHduh/sN2kkJ7ZaRBQVKoHKUykogkVdhJgq+b2WTgN8CPCMpMlwHPA8e4++KwNzSzWcCxQKWZ1QHT3f0WM7uMoFM+G5jp7m8l9Dfpjho3t9O30F4ZqXB36ailjBSdFFRGEpEUC/3p4u4LgU+bWQFQDnzs7lsTvaG7T2nj+GN04z3S97BzZ1BGaq9vYfMHHZSRhrQuI7W0FlRGEpE0S/irqbs30Jv7I5JRRuq/H4w8tnXfgspIItJN9M66xo7tsL0++ObfuDkYhdS4OXi+veXnluOb9uyMbq+M1LcqThkpMo9BZSQR6UF6z6fZmkXw632CpNBq45025JVAfh8oLI+UkQ5u3begMpKI9EK9J3nkFsPoUyC/JJIUoh99IL9v5Fyf3X9mZaU7ahGRjNR7kke/4XDy79IdhYhIj6Cv1iIikjAlDxERSZiSh4iIJEzJQ0REEqbkISIiCVPyEBGRhCl5iIhIwsy9d+yRZGZrgWXpjiOEUjq/uVZXS2fMXXHvZN8jGdfb22vszfsSfU8lwQ6jEk6m/84Pd/dWu8X2muTRXZjZDHeflu44EpHOmLvi3sm+RzKut7fX2Jv3JfoeM5vv7uMTja236o6/86CyVSZ6ON0B7IV0xtwV9072PZJxvb29xt68rzv+P9mddMt/X7U8RCSp1PLoHdTyEJFkm5HuACT11PIQEZGEqeUhIiIJU/IQEZGEKXmIiEjClDxEpMuY2Ugzu8XM7kt3LNI5Sh4iEoqZzTSzNWb2ZszxE8xssZktMbPL27uGu9e6+9TURipdofdsQysinXUbcC1wR8sBM8sGrgOOB+qAeWY2G8gGrox5/wXuvqZrQpVUU/IQkVDcfY6ZjYg5PBFY4u61AGZ2N3Cqu18JnNzFIUoXUtlKRDqjClgR9bwuciwuM6swsxuAQ83silQHJ6mjloeIdIbFOdbmzGN3XwdckrpwpKuo5SEinVEHDIt6PhRYlaZYpAspeYhIZ8wDRplZjZnlAWcBs9Mck3QBJQ8RCcXMZgEvAfubWZ2ZTXX3HcBlwBPA28C97v5WOuOUrqGFEUVEJGFqeYiISMKUPEREJGFKHiIikjAlDxERSZiSh4iIJEzJQ0REEqbkIb2WmR1iZj8zs/J0xxJPJLYuG0tvZreZWV1X3U+6NyUP6c0OAaYDGZk8gJuBI9IdhEg8WhhRJEO5ex3B2lEiGUctD+mxzGw/M3swsvtdg5ktN7O/mFmOmZ0H3Bp56btm5pHHiMh7c8zsCjP7l5k1mtkqM/utmRVEXX9E5D1fNbOrI/fZamaPxNn3Il58nzWzuWa20czqI7vx/TTq/B5lKzN7NirO2MeIqNcdY2ZPmdlmM9tiZk+Y2ZgE/t0ONbPnI3+Xd81Mq+BKK0oe0pM9QrC3xKXAZ4HLgUaC/+8fBX4Zed0XCMpDRwAfRI79GfgxcBdwEsGueFOBO+Pc5wpgFHA+8DVgHPA3M8ttKzAzG0mwgOD7wJnAKcDVQHE7f5+vRsV5BHAU8A6wGlgfue5JwFNAPXAOcDZQAjxvZsPiXDNW38jf+c/AqQQLH15vZseFeK/0Ju6uhx497gFUEuwrcUo7rzkv8pp9Y44fHTl+bszxL0WOHxJ5PiLyfBGQFfW6IyPHp7Zz7zMir+nbzmt+FvyKtnn+WmAbcFjUsSXAUzGv6wt8BPy+g3+z2yIxHRd1LD/y3hnp/m+qR2Y91PKQnmodUAtcZWYXmdmoBN57ArAduD9Svsoxsxzgb5Hzk2Nef5+772x54u4vEvRVtNfZ/RrQBNxtZmeY2YAE4sPMvkbQEjnX3V+OHBsF7APcGRP3VoLVcGPjjmeruz8T9XdpBN4FqhOJT3o+JQ/pkdzdgeOB+QQlp3fMrNbMLg3x9gFAHkHppynqsSZyviLm9avjXGM17WzH6u5LCEppWcCfgA/N7GUzO6aj4Mzs34BrgB+7+19i4ga4JSbuJoL9xGPjjmdDnGONQEGc49KLabSV9FjuXguca2YGHEyw78QfzWypu/9fO29dBzQQlK/iid0pb2Cc1wwkaF20F98zwDNmlk9Q6vo58KiZjXD3j+K9x8wOAO4F/uzu/x0nbgj6YJ6M8/bt7cUjkgglD+nxIq2Q18zsOwSd3mOA/yP4Rg1QGPOWx4EfAKXu/lSIW5xhZj9rKV2Z2ZEE27G+FDK+RuBpM+sDPATUEPQz7MHMKggGAbwOTItzqcXAUuCT7n5VmHuL7C0lD+mRzOwggtLOPQSdyNkEHeQ7gKcjL1sU+fNrZnY7QXnnDXd/NrJr3n1mdjXwCrCToIP834EfuPs7UbcrAf5qZjcC/QnKZO8Cd7QT3yUEfRCPASsIOvivIGjVvNnG2+6MvO7rwNigQbXLq+7eGOkLeSiyJey9BEloIDAJWO7uV7cVk0gilDykp/oQWA58h6AV0AD8EzjZ3RcAuPvrZvYzgm/xFxH0P9QQfHs/h+BD+gLgRwStlKUE263G9nFcCexLMFqpGHgGuMzdm9qJ73XgxMh7BxAMtX0B+JK7b2vjPZ8gGDn1aJxzNcBSd3/MzCZHYr6ZoFX1IfAPgkQqkhTahlZkL0Um5r0PXOTuN6c5HJEupdFWIiKSMCUPERFJmMpWIiKSMLU8REQkYUoeIiKSMCUPERFJmJKHiIgkTMlDREQSpuQhIiIJ+//7/fp7P7xKLAAAAABJRU5ErkJggg==\n",
+ "text/plain": [
+ "<Figure size 432x288 with 1 Axes>"
+ ]
+ },
+ "metadata": {
+ "needs_background": "light"
+ },
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "plt.loglog(grid,error_eukl)\n",
+ "plt.loglog(grid, grid**2, label='2nd oder')\n",
+ "plt.loglog(grid,grid, label='1st oder')\n",
+ "plt.ylabel('error in Euklidean norm', fontsize=16)\n",
+ "plt.xlabel('step size h', fontsize=16)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": []
+ }
+ ],
+ "metadata": {
+ "kernelspec": {
+ "display_name": "Python 3",
+ "language": "python",
+ "name": "python3"
+ },
+ "language_info": {
+ "codemirror_mode": {
+ "name": "ipython",
+ "version": 3
+ },
+ "file_extension": ".py",
+ "mimetype": "text/x-python",
+ "name": "python",
+ "nbconvert_exporter": "python",
+ "pygments_lexer": "ipython3",
+ "version": "3.7.8"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}
--
GitLab