diff --git a/.ipynb_checkpoints/five_point_stencil-checkpoint.ipynb b/.ipynb_checkpoints/five_point_stencil-checkpoint.ipynb index 099e014c0383328e9e37f5d7789f18358158f40b..5eaf5c8bf10f5201e31a311299c2739cf153158e 100644 --- a/.ipynb_checkpoints/five_point_stencil-checkpoint.ipynb +++ b/.ipynb_checkpoints/five_point_stencil-checkpoint.ipynb @@ -7,17 +7,18 @@ "outputs": [], "source": [ "import numpy as np\n", - "from scipy.sparse import csr_matrix\n", - "from scipy.sparse.linalg import spsolve" + "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 is given with Dirichlet boundary condition is given:**\n", + "**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 \\; \\partial \\Omega$$" + "$$ u = g \\quad on \\quad \\partial \\Omega \\quad$$" ] }, { @@ -39,38 +40,118 @@ "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))" + " 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", + "**First, the components of the following equation will be assembled:**\n", "\n", - "$$A \\underline{u} = \\underline{f} + B\\underline{g}$$" + "$$A \\underline{u} = \\underline{f} + B\\underline{g} \\Longleftrightarrow A \\underline{u} - B\\underline{g} = \\underline{f} = -\\Delta u $$" ] }, { "cell_type": "code", - "execution_count": 77, + "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 matrix_A(h):\n", " N = int(1/h)\n", " m = pow(N-1,2)\n", - " A = pow(h,-2)*(np.zeros((m,m))-4*np.eye(m)+np.eye(m,k=1)+np.eye(m,k=-1)+np.eye(m,k=N-1)+np.eye(m,k=-(N-1)))\n", + " A = pow(h,-2)*(np.zeros((m,m))+4*np.eye(m) - np.eye(m,k=1) - np.eye(m,k=-1) - np.eye(m,k=N-1) - np.eye(m,k=-(N-1)))\n", + " \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", + " #for successors:\n", " A[(i+1)*(N-1)-1][(i+1)*(N-1)] = 0\n", + " #for predecessors\n", " A[(i+1)*(N-1)][(i+1)*(N-1)-1] = 0\n", - " #A = csr_matrix(A)\n", + " \n", " return A" ] }, { "cell_type": "code", - "execution_count": 78, + "execution_count": 5, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([[ 64., -16., 0., -16., 0., 0., 0., 0., 0.],\n", + " [-16., 64., -16., 0., -16., 0., 0., 0., 0.],\n", + " [ 0., -16., 64., 0., 0., -16., 0., 0., 0.],\n", + " [-16., 0., 0., 64., -16., 0., -16., 0., 0.],\n", + " [ 0., -16., 0., -16., 64., -16., 0., -16., 0.],\n", + " [ 0., 0., -16., 0., -16., 64., 0., 0., -16.],\n", + " [ 0., 0., 0., -16., 0., 0., 64., -16., 0.],\n", + " [ 0., 0., 0., 0., -16., 0., -16., 64., -16.],\n", + " [ 0., 0., 0., 0., 0., -16., 0., -16., 64.]])" + ] + }, + "execution_count": 5, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# example for matrix A\n", + "n = 2\n", + "h = pow(2,-n)\n", + "A = matrix_A(h)\n", + "A" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [], + "source": [ + "# matrix A in sparse format. this is a lot harder to assemble because of the single deletions\n", + "def sparse_A(h):\n", + " N = int(1/h)\n", + " m = pow(N-1,2)\n", + " A = pow(h,-2)*(4*sp.eye(m) - sp.eye(m,k=1) - sp.eye(m,k=-1) - sp.eye(m,k=N-1) - sp.eye(m,k=-(N-1)))\n", + " \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", + " A.data[(i+1)*(N-1)*5-1*(N-1)-2-1]=0 \n", + " A.data[(i+1)*(N-1)*5-1*(N-1)]=0\n", + " \n", + " return A" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "0.0" + ] + }, + "execution_count": 7, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "#check that sparse matrix is equal to dense matrix\n", + "s=sparse_A(h)\n", + "d=A-s.todense()\n", + "np.linalg.norm(d)" + ] + }, + { + "cell_type": "code", + "execution_count": 8, "metadata": {}, "outputs": [], "source": [ @@ -85,82 +166,175 @@ ] }, { - "cell_type": "raw", + "cell_type": "code", + "execution_count": 9, "metadata": {}, + "outputs": [], "source": [ - "# just the initialisation of matrix B\n", - "def matrix_B(h):\n", + "# 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:\n", + " B[(i+1)*(N-1)-1][2*(N+i-1)] = pow(h,-2)\n", + " #upper boundary:\n", + " B[-(N-1)+i][-N+i] = pow(h,-2)\n", + " B = sp.csr_matrix(B)\n", " return B" ] }, { - "cell_type": "raw", + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "matrix([[ 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n", + " 0., 0., 0.],\n", + " [ 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n", + " 0., 0., 0.],\n", + " [ 0., 0., 0., 0., 0., 0., 16., 0., 0., 0., 0., 0., 0.,\n", + " 0., 0., 0.],\n", + " [ 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n", + " 0., 0., 0.],\n", + " [ 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n", + " 0., 0., 0.],\n", + " [ 0., 0., 0., 0., 0., 0., 0., 0., 16., 0., 0., 0., 0.,\n", + " 0., 0., 0.],\n", + " [ 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 16.,\n", + " 0., 0., 0.],\n", + " [ 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n", + " 16., 0., 0.],\n", + " [ 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 16., 0., 0.,\n", + " 0., 16., 0.]])" + ] + }, + "execution_count": 10, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# example for B\n", + "n = 2\n", + "h = pow(2,-n)\n", + "B = sparse_B(h)\n", + "B.todense()" + ] + }, + { + "cell_type": "code", + "execution_count": 11, "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", - " g[-1] = u(1,1)\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": "markdown", + "cell_type": "code", + "execution_count": 12, "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": 12, + "metadata": {}, + "output_type": "execute_result" + } + ], "source": [ - "Since the exact solution is zero at the boundary, the product B*g is zero." + "# example for g\n", + "n = 2\n", + "h = pow(2,-n)\n", + "g = vector_g(h)\n", + "g" ] }, { "cell_type": "code", - "execution_count": 115, + "execution_count": 13, "metadata": {}, "outputs": [], "source": [ - "n = 6\n", + "n = 3\n", "h = pow(2,-n)\n", "N = pow(2,n)" ] }, { "cell_type": "code", - "execution_count": 116, + "execution_count": 14, "metadata": {}, "outputs": [], "source": [ - "A=matrix_A(h)\n", + "A=sparse_A(h)\n", "f=vector_f(h)\n", - "appr_u = np.linalg.solve(A,f)" + "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": 129, + "execution_count": 15, + "metadata": {}, + "outputs": [], + "source": [ + "RHS = f+B.dot(g)\n", + "appr_u = sp.linalg.spsolve(A,RHS)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + " " + ] + }, + { + "cell_type": "markdown", "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "4096.0" - ] - }, - "execution_count": 129, - "metadata": {}, - "output_type": "execute_result" - } - ], "source": [ - "A[3968,3905]" + "#### **Error and Consistency**" ] }, { "cell_type": "code", - "execution_count": 130, + "execution_count": 16, "metadata": {}, "outputs": [], "source": [ @@ -176,75 +350,120 @@ }, { "cell_type": "code", - "execution_count": 131, + "execution_count": 17, "metadata": {}, - "outputs": [], - "source": [ - "v= exact_solution(h)" - ] - }, - { - "cell_type": "code", - "execution_count": 132, - "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "-8.180678655519896 -8.991060763833046\n", + "-8.270439559649837 -11.480133561857317\n", + "-10.760438683213074 -12.53128157061408\n", + "-12.195353800544211 -12.961737935032156\n", + "-12.821966137009229 -13.145532611303409\n", + "-13.087791494720165 -13.228477586622512\n", + "-13.203877060573376 -13.267565784364171\n" + ] + } + ], "source": [ - "error = max(abs(appr_u - v))" + "grid = np.zeros(7)\n", + "error_eukl = np.zeros(7)\n", + "error_inf = np.zeros(7)\n", + "\n", + "for i in range(7):\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": 133, + "execution_count": 18, "metadata": {}, "outputs": [ { "data": { "text/plain": [ - "0.015864236841443172" + "Text(0.5, 0, 'step size h')" ] }, - "execution_count": 133, + "execution_count": 18, "metadata": {}, "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEKCAYAAAD9xUlFAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAAha0lEQVR4nO3de3RU53nv8e8zuqILEtJIgEbchY2xg5FiG9tAEieu69QmNnKzcnFX0ySN18lybyerPW26mp5eVk7S1SanTW+p2+Y0p3Wck8TgW3xJ0ji1jQ3YSMbgS4zAXCQwIAkE6H55zh8zkgWW0Ahpa8+Mfp+19hpmz8zeP5kxj979vvt9zd0REREBiIQdQEREUoeKgoiIjFBREBGRESoKIiIyQkVBRERGqCiIiMiI7LADTEU0GvWlS5eGHUNEJK3s2rWr1d0rxnotrYvC0qVLeemll8KOISKSVszs0Hiv6fKRiIiMUFEQEZERaVkUzGyTmd3X0dERdhQRkYySlkXB3R9193tKSkrCjiIiklHSsiiIiEgwVBRERGREWg9JvVRNJ85ypL2beYW5lBXkUlaUS2FuFmYWdjQRkVDNyqLw8MtH+dufNp23Lzc7QllBLvMKcykvHPWYKBplBbmUFb6zzSvIITtLDS0RySyzsih86salfHBVJe2dfe9sXX20n+vjVFcfbZ19NJ/qor2zjzM9A+MeZ25+NuVFecwryKGsMI+ywnEe1RoRkTQxK4tCtCiPaFFeUu/tHxziVFffeQXkVGe8cJzq7KO9q5/2zl5aTnezp+U07Z199A+OvZpdblYk3soY1Ropu6CgzCvMoTzxOK8glxy1RkRkBs3KojAZOVkRKovzqSzOT+r97s653gFOdfbT1tmbKCjxwnHhY8vpbtrO9U66NTL60lZ5US7z5+azpLyQojz9dYrI1OhfkWlmZhTn51Ccn8Pi8oKkPtM/OMTprv6LXs461dnH0dPd7G3poL2zj77BoXcdJ1qUx9LyApaUF8Yfo4Ujz0vm5Ez3jyoiGUhFIQXkZEWoKM6joji5S1ruTmff4MhlrGOnuznY1sWhtk4OtnWyramVBxt6zvvMvIKcd4pFeSFLo/HHZeWFlBbkqK9DRAAVhbRkZhTlZVOUl82isgLWLip913u6+wY53N7FwbbORLGIF40XD57i4d1H8VHdHnPzs1lSXsiS8gKWDj9G448VRXkqGCKzSFoWBTPbBGyqqakJO0rKmpObxeULirl8QfG7XusdGORIe/d5xeJgWxd7Wjp4Yu/bDA69UzEKcrPOb2GMamnML84nElHBEMkk5j72SJl0cM0117jWU5he/YNDtJzqTrQwus57PNLedd7IqrzsCEsuLBaJlkZV6RyyVDBEUpKZ7XL3a8Z6LS1bChKcnKwIS6OFLI0Wvuu1wSHn6OnuUcXinZbGM2+epHfgnc7v3KwI1WVz3rkcNeoxNm+OhtqKpCgVBUlaVsRYVFbAorICNqyMnvfa0JBz/GwPB1u73nVZavuBNrr6Bs87TvW8OWNelloWLVQLQyREKgoyLSIRY2HJHBaWzOGGFeXnvebunDzXG29htMZHSA0XjcZDpzjb+859GiVzcrhheTnrV0ZZv6KcZdFCdXSLzCAVBQmcmY3cAHjt0rLzXnN32jv7ONjWxYGT59j5Vjvbmlp58tW3AVhYks/6mijra8pZvyJK5dzkbiIUkUujjmZJOe7OwbYunmtq5fmmVl440Mbprn4AVlYWJYpElHXLy5ibr5vyRCbrYh3NKgqS8gaHnNeOnmHb/la2NbXy4sF2evqHiBisqS5lQ02UG2vKqVs8j/ycrLDjiqQ8FQXJKL0DgzQcOs3ziSKxu7mDwSEnLzvCdcvKuHFF/HLTlVUl6rQWGYOKgmS0sz397DjQPtKSePP4OWBUp3VNOetrouq0FknQfQqS0Yrzc7h59XxuXj0fgBNne3hhfxvP7Wvl+f1t53Va37giyoaV6rQWGY9aCpLR3J1Dw53W++NFYrjTuqayKN4fsaKc61eUq9NaZg1dPhJJGBpyXjt2hm1NrWzb38bOt9rO67QeHvpat0Sd1pK5Mq4ojJoQ73P79u0LO46ksd6BQRoPn+b5plaeu6DT+tqlZdxYU86Gmqg6rSWjZFxRGKaWgky3sz39iRvo2tjW1MrPj58F4tOL37CiPDH8NcpydVpLGlNHs0iSivNz+NAV8/nQFfFO65Nne0eGvm5rauOpV48DsGDuqDuta6LMV6e1ZAi1FESS5O4cbh++07qN5/e3cmpUp/X6FeWJO63LtfyppDRdPhIJwHCndbwl0cbOt9rp7h8kYvCe6lI2JqbjeO+SeeRma6pwSR0qCiIzoG9giMbDp9i2P94f8fKR0wwOOXNysli3vIwNNVE2rqzgsvlF6o+QUKkoiITgTE8/2/e38VxiZNOBk50AVBTnsaEmGt9Wqj9CZp6KgkgKaDndzbZ9rTzbFO+4bu/sA+Cy+UVsqKlg48oo1y0rozBP4z8kWCoKIilm9E10zzW1svOtdnoHhsjJMmoXz2NjohWxprpU90fItFNREElxPf2DvHTwFM82nWRbUyt7W84A8fsj4vM1xS83LSkvUH+ETJnuUxBJcfk5WfF/+BNrX7d39sVbEfviLYnhSf2q581h48ooG2oquHFFOfMKc8OMLRlILQWRFOfuvNXaybamVp7d18oL+9s42zuAGbwnVsL6migba6K8d+k88rI1X5NMTJePRDLIwOAQu5s7Eq2IkzQePs3AkJOfE+G6ZeUj/RGrFhTrUpOMSUVBJIOd6x1gx4E2nk1camo6EV9kKFqUy/rE0NeNKytYUKKhrxKnPgWRDFaUl33efE3HOrrZ1tTGc/tO8lxTGw+/fBR4Z/2IjSvjU3EUaeirjEEtBZEM5u688fZZnkvcHzG8fkR2xKhdXMqGmgo2rIxydXUJ2VmaimO20OUjEQHi60fsOnRqZFTTnpYO3KE4P5sblpePDH3VetaZLeOKghbZEZkepzr7eH5kKo6THGnvBiBWOocNNVHWr4yPbNLQ18yScUVhmFoKItPrUFtnvMN6X3xN6zM9A+RkGbdetZC71y1m3bIytSAygIqCiEza4JCzp6WDR14+yg92HeFMzwArKgq5e90S7qqrpqRAa0akKxUFEZmS7r5BfrjnGPfvOETj4dPkZUfYdHUVd69bzNpFpWo9pBkVBRGZNq8e7eA7Ow7zUGMLnX2DrF44l7uvX8wda2Ma5pomplQUzKwU+FVgKaPua3D335q+iJdGRUEkPOd6B3j45Rb+Y/thXj92hsLcLO6sjXH3uiWsrpobdjy5iKkWheeB7cAeYGh4v7t/ezpDXgoVBZHwuTsvHznN/TsO8+juo/QODFG7uJS71y3h9jULyc/RfEypZqpFocHd6wJJNkUqCiKppaOrnwcbmrl/xyH2n+ykZE4Od9VV88l1i6mpLAo7niRMtSj8d+Ac8BjQO7zf3dunM+SlUFEQSU3uzo632rl/x2Ge3HuM/kHn+uVl3L1uCb945QJys3X3dJimWhTuBb4MnAaG3+zuvnw6Q14KFQWR1Nd6rpfvv9TMd3Ye4kh7N9GiXD56zSI+ed1iFpUVhB1vVppqUdgPrHP31iDCTYWKgkj6GBpynm1q5f7th/jJ68dx4H0rK7h73WI+uKpScy/NoKnOkvoq0DW9kURktolEjPdfVsH7L6vgWEc3/+/FI3x35xHu+fddLJibz8evW8THr12sKb5DlkxLYStwJfA05/cpaEiqiEzJwOAQP33jBPfvOMwz+04SMeNDqyq5+/olbKyJEonoprggTLWl8FBiExGZVtlZEW65cgG3XLmAw21dPPDiYb734hF+9NpxFpcV8InrFvPRa6qJFuWFHXXWuGhLwcyygKfc/eaZi5Q8tRREMk/vwCA/evU49+84xPYD7ZqQLwCX3FJw90Ez6zKzEnfvCCaeiMg78rKz2HR1FZuurqLpxFm+s+MIP9h1hEd3H9WEfDMgmT6F7wHXAz8GOof3q09BRGaKJuSbXlMdkvqpsfZrmgsRCYMm5Ju6Kc+Sama5wGWJpz939/5pzHfJVBREZi9NyHfpptpS+ADwbeAgYMAi4FPu/sy0ppwELccpIsM0Id/kTbUo7AI+6e4/Tzy/DHjA3d877UknSS0FERlNE/IlZ6pF4RV3XzPRvjCoKIjIWMabkO+PblvNVbGSsOOFbqpF4VvEJ8L798Suu4Fsd//0tKa8BCoKIjKR4Qn5/uXZA1SXFfDwvevDjhS6qd7R/HngXuC3iPcpPAP8w/TFExEJTrQoj89/YAXZEePLj7/O/pPnWFGhS0njmXBaQnfvdfevu3u9u2929//t7r0TfU5EJJXcsbaKiMHWhpawo6S0CYuCma03sx+b2ZtmdmB4m4lwIiLTpXJuPhtWVrC1sYWhoYmH4s9WyUxg/q/A14ENwLWjNhGRtFJfG6PldDc7D4a+cGTKSqYodLj7E+5+wt3bhrfAk4mITLNbrpxPYW6WLiFdRDJF4Wkz+0szu8HM6oa3wJOJiEyzgtxsbr1qIY/vOUZP/2DYcVJSMqOP1iUeRw9fcuCD0x9HRCRYd9XFeLChmR+/dpxNV1eFHSflTFgU3P2mmQgiIjITrl9ezsKSfLY0NKsojGFSK2Wb2WNBBRERmQmRiHFnbYxn9rVy8qxG119oUkUBiAWSQkRkBtXXxhgcch7ZfTTsKClnskWhMZAUIiIzaOX8Yt4TK2FrY3PYUVLOpIqCu38mqCAiIjNpc22MvS1nePP42bCjpJRk7mi+3cwazazdzM6Y2VkzOzMT4UREgvKRtVVkRYwtumfhPMm0FP4a+BRQ7u5z3b3Y3bWskYiktWhRHu+/rIKHGlsY1LQXI5IpCkeAvZ7Mup0iImmkvi7G22d62H5AkzQMS+bmtf8BPG5m/wWMjN9y968HlkpEZAbcfMV8ivOyebChmfU10bDjpIRkWgpfBrqAfKB41CYiktbyc7K4bc1Cntz7Nl19A2HHSQnJtBTK3P2WwJOIiIRgc22M7754hKdefZvNtdVhxwldMi2Fn5iZioKIZKRrl5ZRPW+ORiElJFMU7gWeNLNuDUkVkUwTiRiba2Nsa2rl+JmesOOELpnlOIvdPeLuczQkVUQy0ebaGEMOD7+s1sKEfQpm9r6x9rv7M9MfR0Rk5i2vKGLtolK2NLRwz/tWhB0nVMl0NP/eqD/nA9cBuwhxPQUz2wRsqqmpCSuCiGSYu+pifOnhV3nt6BlWV83eiyHJXD7aNGr7BeAq4Hjw0S6a6VF3v6ekpCTMGCKSQW5fU0VOlrGlYXZPkjfZWVIBmokXBhGRjDGvMJebLq/k4d1HGRgcCjtOaJLpU/hb4stvQryIrAV2B5hJRCQU9XUxfvTacZ5rauUDl1eGHScUyfQpvDTqzwPAA+6+LaA8IiKhuWlVJSVzctja2KKiMB53//ZMBBERCVtedha3r1nIgw3NnOsdoCgvmd+bM4vWUxARGaW+rpqe/iGe2HMs7Cih0HoKIiKj1C0uZWl5wayd9kLrKYiIjGJmbK6tZvtbbbSc7g47zoxLpigMr6fwRTP7wvAWdDARkbBsro3hDg81zr7WgtZTEBG5wOLyAq5dOo+tjS3MtoskWk9BRGQMm2ur+cOte9jT0sGa6tKw48wYracgIjKG296zkNzsyKzrcNZ6CiIiYygpyOHmKyp5dPdR+mfRtBdaT0FEZBz1tdW0dfbxzJsnw44yY8btUzCzVe7+hpnVjfW6uzcEF0tEJHzvv7yCssJctjS08KEr5ocdZ0ZcrKP5C8A9wNfGeM0JcT0FEZGZkJMV4SNXV/GdnYfp6O6nZE5O2JECN+7lI3e/J/F40xibCoKIzAqba2P0DQzx+CyZ9iKZuY/+3MyyRj2fa2b/J9hYIiKpYU11CSsqCtk6S0YhJTP6KBvYaWZrEkNTXyS+HKeISMYzM+rrqtl5sJ0j7V1hxwlcMqOPvgj8PrAD+DfgNnf/u4BziYikjDtrYwBsnQXTXiRz+eh9wN8Afwb8DPg7M6sKOJeISMqIlc7h+uVlbGlozvhpL5K5fPRXwEfd/Svu/kngPuCnwcYSEUkt9XXVHGzrovHI6bCjBCqZonCDu782/MTdtwDrg4skIpJ6PnzVAvKyI2xpaA47SqCSWY5z0MxuA64kPlPqsD8LLJWISIopzs/hF69cwGOvHONLt68mLztr4g+loWT6FL4JfAz4TcCAjwJLAs4lIpJyNtfFON3Vz9NvZO60F8lcPrrR3X8VOOXufwrcACwKNpaISOrZWBMlWpTH1sbMvYSUTFEYXo+uKzHqqB9YFlwkEZHUlJ0V4Y61Vfz0jROc6uwLO04gkikKj5lZKfCXQANwEPhugJlERFJWfV2M/kHnsQyd9iKZm9f+3N1Pu/uDxPsSVrn7l4KPJiKSelYvnMvl84szdhTShKOPEvMe3QYsHX6/meHuXw82mohI6olPexHjK0+8wVutnSyLFoYdaVolc/noUeDXgHKgeNQmIjIr3bE2hhlszcDWwoQtBaDa3dcEnkREJE0sKMlnQ02ULY0t/M7NlxGJWNiRpk0yLYUnErOjiohIwubaGM2nunnp0Kmwo0yrZIrCdmCrmXWb2RkzO2tmZ4IOJiKSyn7xygUU5GZl3D0LyRSFrxG/Ya3A3ee6e7G7zw04l4hISivMy+bWxLQXPf2DYceZNskUhX3AXk+h+WLNbJOZ3dfR0RF2FBGZxTbXxTjbM8B/vn4i7CjTJpmicAz4mZl90cy+MLwFHexi3P1Rd7+npKQkzBgiMsvduCLK/Ll5GXXPQjJF4S3gP4FcNCRVRGREVsS4c22M/3rzJK3nesOOMy2SmTr7T2ciiIhIOqqvq+afnjnAo7uP8un16T8tXDItBRERGcflC4pZvXBuxqzfrKIgIjJF9XUxXmnuoOnE2bCjTJmKgojIFH1kbRURgy0N6d9aSGZCvArgc4yaEA/A3T8TXCwRkfRRWZzP+y6r4KHGFn73lsvTetqLZFoKDwMlwE+AH47aREQkYXNtjKMdPWx/qy3sKFOSzIR4Be7++4EnERFJY7esXkBRXjZbGlq4cUU07DiXLNmV134p8CQiImlsTm4WH75qAU/sOUZ3X/pOe5FMUfht4oVBE+KJiFxEfV01nX2D/Oi1t8OOcsmSWY6z2N0j7j5HE+KJiIxv3bIyYqVz0noU0rh9Cma2yt3fMLO6sV5394bgYomIpJ9IxLiztop//Nl+TpzpoXJuftiRJu1iHc1fAO4hPnX2hRz4YCCJRETS2Obaav7+6f08svsov75xedhxJm3couDu9yQeb5q5OCIi6a2msoirq0t4sKElLYuC7mgWEZlmm2tjvH7sDK8fS78xOSoKIiLTbNPVVWRHLC0nybtoUbC4RTMVRkQkE5QX5fGBy+PTXgwOpcyilUm5aFFILMH50MxEERHJHPV11Zw428u2ptawo0xKMpePtpvZtYEnERHJIB9cVUlxfnbaXUJKpijcBLxgZvvN7BUz22NmrwQdTEQkneXnZHH7miqe3Ps2nb0DYcdJWjIT4n048BQiIhmovi7GAzsP8+Tet7nrvdVhx0lKMtNcHAJKgU2JrTSxT0RELuKaJfNYVDaHLY3NYUdJ2oRFwcx+G7gfqExs/2Fmvxl0MBGRdGdmbK6t5vn9bRzr6A47TlKS6VP4LLDO3f/Y3f8YuJ74SmwiIjKB+toY7vBQ49GwoyQlmaJgwOjJwQcT+0REZAJLo4XULS5lS0Mz8VH+qS2ZovAtYIeZ/YmZ/QmwHfjXQFOJiGSQ+rpq9p04x6tHU3/ai4nuaI4AO4BPA+3AKeDT7v7XwUcTEckMt69ZSG5WJC3WWZjojuYh4Gvu3uDu33D3v3H3xhnKJiKSEUoLcvngqkoe2d3CwOBQ2HEuKpnLRz8ys7vMTP0IIiKXaHNdjNZzfTy7L7WnvUimKHwB+D7QqzWaRUQuzU2XV1JakMODDal9z0IyfQq3JtZoztUazSIilyY3O8KmNVX8+LXjnOnpDzvOuJLpU/irGcoiIpLR6uti9A4M8cSeY2FHGZf6FEREZsjaRaUsixam9CikyfQp9KlPQUTk0pkZ9bUxdrzVzpH2rrDjjCmZCfGKE30KOepTEBGZmjtrYwA8/HJqthaSmRDPzOxXzOxLieeLzOy64KOJiGSeRWUFXLesjC0NLSk57UUyl4/+AbgB+GTi+Tng7wNLJCKS4eprYxxo7WR3c0fYUd4lmaKwzt3vBXoA3P0UkBtoKhGRDPZLaxaSmx1hSwres5BMUeg3syzAAcysAkjt+7RFRFLY3PwcfmH1fB7dfZS+gdT65zSZovANYCtQaWZfBp4D/legqUREMtxddTFOdfXzs5+fCDvKeSZco9nd7zezXcCHiK+jcKe7vx54MhGRDLZxZQXlhblsbWzhlisXhB1nxIRFAcDd3wDeCDiLiMiskZMV4SNrq7h/+2E6uvopKcgJOxKQ3OUjEREJQH1tNX2DQzy2J3WW6lRREBEJyVWxuaysLEqpaS9UFEREQmJmbK6LsevQKQ61dYYdB1BREBEJ1Z1rY5jB1sbUaC2oKIiIhKiqdA43LC9na2NqTHuhoiAiErL6umoOtXXRcPhU2FFUFEREwnbrVQvIz4nwYAp0OKsoiIiErCgvm1uvXMBju4/SOzAYahYVBRGRFLC5rpozPQP89PVwp71QURARSQHrV5RTWZzHlpBHIakoiIikgOysCHesreLpN07Q3tkXWg4VBRGRFFFfV83AkPPYK+FNe6GiICKSIq5YOJdVC4pDHYWkoiAikkLuqqtm95HT7D95LpTzp0xRMLNCM/u2mf2zmd0ddh4RkTDcsbaKiMHWkFoLgRYFM/uWmZ0ws70X7L/VzH5uZk1m9geJ3fXAD9z9c8BHgswlIpKqKufms2FlBVsbWxgamvlpL4JuKfwbcOvoHYn1nv8e+DCwGviEma0GqoEjibeFe/eGiEiI6mtjtJzuZufB9hk/d6BFwd2fAS78qa4Dmtz9gLv3Ad8F7gCaiReGi+Yys3vM7CUze+nkyZNBxBYRCdUtV86nMDcrlEtIYfQpxHinRQDxYhADtgB3mdk/Ao+O92F3v8/dr3H3ayoqKoJNKiISgoLcbG69aiGP7zlGT//MXjgJoyjYGPvc3Tvd/dPu/nl3v3/GU4mIpJC76mKc7R3gx68dn9HzhlEUmoFFo55XA6mzQKmISAq4fnk5C0vy2dLQPKPnDaMovAisNLNlZpYLfBx4JIQcIiIpKxIx7qyN8cy+Vk6e7Z258wZ5cDN7AHgBuNzMms3ss+4+APwG8BTwOvA9d381yBwiIumovjbG4JDzyO6Zu5iSHeTB3f0T4+x/HHg8yHOLiKS7lfOLeU+shK2NzXx2w7IZOWfK3NEsIiLvtrk2xt6WM7x5/OyMnE9FQUQkhX1kbRVZEWPLDN2zoKIgIpLCokV5vP+yCh5qbGFwBqa9SMuiYGabzOy+jo6OsKOIiASuvi7G22d62H6gLfBzpWVRcPdH3f2ekpKSsKOIiATu5ivmU5yXzYMzcM9CWhYFEZHZJD8ni9vWLOTJvW/T1TcQ6LlUFERE0sDm2hhdfYM89erbgZ5HRUFEJA1cu7SM6nlzAh+FpKIgIpIGIhFjc22MbU2tHD/TE9x5AjuyiIhMq821MYYcHn45uNaCioKISJpYXlHE2kWlgV5CUlEQEUkjd9XFeOPts7x29Ewgx1dREBFJI7evqSJalMehts5Ajh/oLKlBMbNNwKaampqwo4iIzKh5hbns+MMPkRUZaxHLqUvLloLuaBaR2SyoggBpWhRERCQYKgoiIjJCRUFEREaoKIiIyAgVBRERGaGiICIiI1QURERkhLkHv+ZnUMzsJHBoiocpAaayruelfn4yn0v2vRO9b6LXo0BrkplS3VT/XlPlnNNxzEs5Rhjfz2Tekynf0TC+n6PPu8TdK8Z8h7vP6g24L4zPT+Zzyb53ovcl8fpLYf99pMrfa6qcczqOeSnHCOP7mcx7MuU7Gsb3M9nz6vIRPBrS5yfzuWTfO9H7pvqzppMwftYgzjkdx7yUY4Tx/ZzsedNZWD/nhOdN68tHMr3M7CV3vybsHCLj0Xc0eGopyGj3hR1AZAL6jgZMLQURERmhloKIiIxQURARkREqCiIiMkJFQZJiZnea2T+b2cNmdkvYeURGM7PlZvavZvaDsLOkOxWFWcDMvmVmJ8xs7wX7bzWzn5tZk5n9wcWO4e4PufvngF8DPhZgXJllpun7ecDdPxts0tlBo49mATN7H3AO+L/uflViXxbwJvALQDPwIvAJIAv4ygWH+Iy7n0h87mvA/e7eMEPxJcNN8/fzB+7+yzOVPRNlhx1Agufuz5jZ0gt2Xwc0ufsBADP7LnCHu38FuP3CY5iZAV8FnlBBkOk0Hd9PmT66fDR7xYAjo543J/aN5zeBm4FfNrP/FmQwESb5/TSzcjP7JlBrZl8MOlwmU0th9rIx9o17LdHdvwF8I7g4IueZ7PezDdAvK9NALYXZqxlYNOp5NXA0pCwiF9L3MyQqCrPXi8BKM1tmZrnAx4FHQs4kMkzfz5CoKMwCZvYA8AJwuZk1m9ln3X0A+A3gKeB14Hvu/mqYOWV20vcztWhIqoiIjFBLQURERqgoiIjICBUFEREZoaIgIiIjVBRERGSEioKIiIxQURABzOx3zKwgwONXTcdc/2Z2bjryiIxH9ymIAGZ2ELjG3VvDznIxZnbO3YvCziGZSy0FmVXMrNDMfmhmu81sr5l9zMx+C6gCnjazpxPvu8XMXjCzBjP7vpkVJfYfNLO/MLOdia1mjHO838xeTmyNZlZsZkuHF5Exs38Z9fpJM/ufif2/Z2YvmtkrZvanF/kZvpzIv93M5gfx30lmLxUFmW1uBY66+9WJBV2eTMwAexS4yd1vMrMo8EfAze5eB7wEfGHUMc64+3XA3wF/PcY5fhe4193XAhuB7tEvuvuvJ167A2gD/i2xxOlK4usIrAXem1h85kKFwHZ3vxp4Bvjc5P8TiIxPRUFmmz3AzYnf9je6e8cY77keWA1sM7OXgU8BS0a9/sCoxxvG+Pw24OuJFkhpYh6f85hZPvB94Dfc/RBwS2JrBBqAVcSLxIX6gMcSf94FLB3/RxWZPK2nILOKu79pZu8Ffgn4ipn9yN3/7IK3GfBjd//EeIcZ58/D5/iqmf0wcY7tZnYz0HPB274JbHH3n4w651fc/Z8m+BH6/Z2OwEH0/7BMM7UUZFYxsyqgy93/A/groC7x0lmgOPHn7cD64f4CMysws8tGHeZjox5fGOMcK9x9j7v/BfFLT6sueP1eoNjdvzpq91PAZ0b1XcTMrHIKP6rIJdFvGTLbvAf4SzMbAvqBzyf23wc8YWbHEv0KvwY8YGZ5idf/iPhC8gB5ZraD+C9VY7UmfsfMbiL+m/xrwBPAwlGv/y7Qn7g0BfBNd/+mmV0BvBBfDptzwK8AJ6b6A4tMhoakikxCugxdFblUunwkIiIj1FIQEZERaimIiMgIFQURERmhoiAiIiNUFEREZISKgoiIjFBREBGREf8f6yXLYPU7gx4AAAAASUVORK5CYII=\n", + "text/plain": [ + "<Figure size 432x288 with 1 Axes>" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" } ], "source": [ - "appr_u[-1]" + "plt.loglog(grid,error_inf)\n", + "plt.ylabel('error in maximum-norm')\n", + "plt.xlabel('step size h')" ] }, { "cell_type": "code", - "execution_count": 134, + "execution_count": 19, "metadata": {}, "outputs": [ { "data": { "text/plain": [ - "13.161396848144852" + "Text(0.5, 0, 'step size h')" ] }, - "execution_count": 134, + "execution_count": 19, "metadata": {}, "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAEKCAYAAAASByJ7AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAAlFUlEQVR4nO3dd3xUdb7G8c83DQgldJBekiAoUqQpSFexsLq7ruBaaKKyKmBh7+7qrvdu0ysW7Ig0QcRV71rWXQUL0qugoggJhCKi0gMB0n/3j0zcyIZkkszkzGSe9+uVVzJnZs55oqNPzjm/8zvmnENERCTK6wAiIhIaVAgiIgKoEERExEeFICIigApBRER8VAgiIgJAjNcBysPMhgPDa9euPT45OdnrOCIiYeWTTz456JxrdPpyC+frEHr06OE2bNjgdQwRkbBiZp8453qcvlyHjEREBFAhiIiIjwpBREQAFYKIiPioEEREBFAhiIiIT1gWgpkNN7MZ6enp5Xr/J7sPs/W7YwFOJSIS3sKyEJxz/3DO3ZKQkFCu9z/4r60Mm7acCS99wlffqhhERCBMC6GiZo7qwcTBiSxPPchlTyzntvmfsGWfikFEIltEX6l89GQ2s1fsZM7KXRzPymXYOU2ZOCSJTs3qBDCliEhoOdOVyhFdCIXST+Ywa+VO5qzYyfGsXC49pwkThyRxTrPyHZISEQllKgQ/pJ/MYfbKncxeuZPjmblc0qkJk4aqGESkalEhlEH6qRzmrNzJrBUFxXBxpyZMGpLEuc1VDCIS/lQI5ZB+Koe5K3cxc0UaxzNzGdqxCZOHqhhEJLypECqgsBhmrUjjmIpBRMKcCiEAjmX69hiWFxZDYyYNSaZzCxWDiIQPFUIAHcvM4cWVu5i5Yifpp3IYcnZjJg1N4rwWdSs9i4hIWVWpQii8hWZiYuL41NRUz3Icz8zhxVW7eGF5QTEMPrsxk4Yk0aVlXc8yiYiUpkoVQqFQuYXm8cwc5q3ezQvL0zh6ModBHRoxaWgyXVUMIhKCVAiVQMUgIuFAhVCJMrJyeXFVwcnnIydzGNihEZOGJNGtVT2vo4mIqBC8kJGVy7zVu3hhWUExDEhuxKShSXRXMYiIh1QIHsrIymX+6t3MWLaDIydz6J9csMdwfmsVg4hUPhVCCDiRlcv8NbuZsSyNwyeyuSipIZOHJqsYRKRSqRBCSPHFkMT5ret7HU1EIoAKIQSdyMrlJV8xHPIVw6QhSfRoo2IQkeBRIYSwk9kFxfD80oJi6JdYsMegYhCRYFAhhIGT2bksWLOH55ft4GBGNn0TGzB5aDI9VQwiEkAqhDByKjuPBWt3M31pQTFc2L6gGHq1VTGISMWpEMLQv4shjYMZWVzYvgGThiTRu10Dr6OJSBhTIYSx04vhgnYNmDQ0iT4qBhEphypVCKEy22llO5Wdx8vr9jB96Q4OHM+iT7v6TB6arGIQkTKpUoVQKFL2EE6XmZPHy2v38FyRYrj74g46xyAiflEhVEGFxfDsxzs4mJFFv8SG3HWxrnwWkZKpEKqwwnMMz328g0MnshmQ3Ii7Lta02yJSPBVCBDiZncu81bt5fmnBJHpDzm7MXRcnc25z3fNZRP5NhRBBCu/HMGNZGumncrikUxMmD02mU7M6XkcTkRCgQohAxzJzmLNiFzNXpHE8M5fLOzdl0pBkOjSt7XU0EfGQCiGCpZ/MYdaKNGav3MWJ7FyuPK8Zk4YkkthYxSASiVQIwpET2bywPI25q3ZxKiePq7o0Y+KQJNo1quV1NBGpRCoE+cGhjCxmLE9j3qrdZOXm8bPuLZg4OIlWDeK9jiYilUCFIP/hwPEsnl+6g/lrdpOb77imewvuGJxIy/oqBpGqTIUgZ7T/WCbPfryDl9ftIT/fcW3PltwxKJFmdWt4HU1EgkCFIKX6Nv0Uzy7ZwSvr92AYI3u15FcDE2maUN3raCISQCoE8ds3R0/x9EfbeW3D10RFGdf3bsWEge1pXFvFIFIVqBCkzL4+fJKnP9rO6xv3Ehtt3NinNbcOaE/DWtW8jiYiFaBCkHLbdfAET36UypubvqFaTDSjLmzDLf3bUb9mnNfRRKQcVAhSYTsOZPDUh6m89dk+4mOjGd23DeMvakfdeBWDSDhRIUjApH5/nCc+TOWdz7+ldrUYxvRry7h+bUmoEet1NBHxQ4ULwczqADGFj51zhwMXr3xUCN7a+t0xnvgglXe/+I7a1WMYf1E7xvRtQ+3qKgaRUFbuQjCzW4E/AqeAwhc751y7gKcsIxVCaPhyXzrTPkjl/S3fUzc+lvEXtWP0hW2oWS2m9DeLSKWrSCGkAhc45w4GK1x5qRBCy+a96Tz+QQofbd1P/Zpx3Nq/HTde0Jr4OBWDSCg5UyFE+fHeHcDJwEeSqqZziwRmj+7Jm7f3pXPzBB58dyv9H17CzOVpZObkeR1PRErhzx5CN2AOsBbIKlzunJsY0CBmVwNXAI2BZ5xzi0t7j/YQQtsnuw/z+PuprNh+kEa1q3H7wPaM7NWK6rHRXkcTiWgVOWS0DlgBbAbyC5c75170Y6OzgSuB/c65c4ssHwY8AUQDM51zDxV5rh7wiHNuXGnrVyGEh7Vph3js/RTW7jxM0zrVuX1Qe67t2ZJqMSoGES9UpBBWOecuLOdG+wMZwLzCQjCzaCAFuBjYC6wHrnPObfE9/yiwwDm3sbT1qxDCy6odB3lscQobdh+hWUJ17hicxDXntyAuxp8jlyISKBU5h7DEzG4xs7PMrH7hlz8bdc4tA04fntoL2O6cS3POZQOvAFdZgf8F3i2pDHxZNpjZhgMHDvgTQ0LEhe0b8tptFzB/XC+aJFTnd29sZvCjH/Pq+q/JycsvfQUiElT+7CHsLGax38NOzawN8E6RPYRrgGHOuZt9j28EelOw1zCKgj2GT51z00tbt/YQwpdzjo9TDvD4+yl8vjed1g3imTg4iau7NSc6yryOJ1KlnWkPocTxgGYWBfzGOfe3QGYpZplzzj0JPBnA7UgIMzMGdWjMwORGfPjVfh7/IIV7XvuM55bu4J6Lkxl2blPMVAwilanEQ0bOuXzg9gBvcy/QssjjFsC+AG9DwoSZMbRTE965sx/PXd8dgAkLNnLVMytZnnqAcJ5aRSTc+HMO4X0zu9fMWpb1HMIZrAeSzKytmcUBI4G3y7ICMxtuZjPS09MrEENCiZlxWeezWDS5P1OvOY9DGdncOGsd172who17jngdTyQiBPUcgpktBAYCDYHvgQecc7PM7HJgGgXDTmc75/5SxtyAziFUZVm5eSxcu4enl2znYEY2Qzs2YcqlHejQtLbX0UTCnmY7lbB0IiuXOSt38vyyNDKycrm6a3PuGppMqwbxXkcTCVsVuQ4hFpgA9Pct+hh43jmXE+iQZaVCiBxHT2YzfWkac1ftJDfPMbJXSyYOTqJxHd3WU6SsKlIIM4FYoPDK5BuBvMJho14ws+HA8MTExPGpqalexRAPfH8sk6c+SuWVdV8TE22MvrAttw3QTXpEyqIihfCZc65Lacu8oD2EyLXn0Eke/yCFNz/9hlrVYri1fzvG9G2rKbdF/FCRK5XzzKx9kRW1AzR1pXiqVYN4Hh/RlXcnXUSfdg14ZHEKA6YuYc7KnWTl6uMpUh7+7CEMoWC20zQKLiprDYxxzi0JfrySaQ9BCm3cc4Sp721jddohmtetweShSfysewtd9SxSjAqNMjKzakAHCgphq3Muq5S3VAoVghTlnGPF9oNMXbSNz/emk9i4lq56FilGRQvhQqANP76n8rxABiwLnVSWkjjnWPTldzyyOIXt+zM4r0UCUy7tQL/EhioGESp2Unk+0B74lH+fO3CBvkFOeWgPQUqSl+/4+8a9TPsglW+OnqJPu/r8etjZdG9Vz+toIp6qSCF8BXRyIXgFmwpB/FHcVc/3XprM2U3reB1NxBMVGWX0BdA08JFEKke1mGhG923L0imDuPeSZNbuPMRlTyxn8iub2H3ohNfxREKGP3sIS4CuwDp+fE/lnwQ1mR+0hyDlcfpVzyN6tmTikCSa6KpniRAVOWQ0oLjlzrmlAcpWZjqpLIGw/1gmTxa56nnUhW2YMKC9rnqWKk+T24mcga56lkijQhApxdbvjvHo4hTe3/I9DWvFcfugRH7ZuxXVYqK9jiYSUCoEET8Vd9XzT7s1JybanzEYIqGvIqOMRCJK91b1eHl8b14a15sGteKY8vrnXDptGe9u/la39JQqrdRCMLO+Zva+maWYWZqZ7TSztMoIJ+IVM6NfUkPeur0v02/ojpkxYcFGfvL0Spal6F7PUjX5M8poK3AX8AlFZjl1zh0KbrTS6ZCRVJa8fMcbm77h8fdTdNWzhL2KDDtd65zrHbRk5aBhp+IVXfUsVUFFCuEhIBr4Oz++MG1joEOWlfYQxCsnsnKZu2oX05fuICMrl6u6NOOui5Np3aCm19FESlWRQijuvgfOOTc4UOHKS4UgXtNVzxKONOxUJIj2H8vkqY+2s3Ddnh/u9fyrQe2pUz3W62gi/6Gi90O4AjgH+OHPHufcHwOasBxUCBJqil71XD8+jrsvSWZkz1a6c5uElHJfh2Bm04ERwJ0U3DHtFxTcRlNETlN4r+e3b+9H+0a1uO+NL7jiyeWs3H7Q62gipfLnwrQLnXM3AUecc/8DXAC0DG4skfDWuUUCf7u1D89e352MrFyun7mWm1/cwM6Dmm5bQpc/hXDK9/2kmTUDcoC2wYskUjWYGZd3PosP7h7Afw07m9U7DnLJ40v58ztbSD+V43U8kf/gTyG8Y2Z1ganARmAX8EoQM5XKzIab2Yz09HQvY4j4pXpsNBMGtmfJlIH8vHsLZq3cycCpS5i/ehe5eflexxP5QZlGGZlZNaC6cy4k/k+sk8oSjr7cl86f3tnCmrTDJDWuxe+v7ET/5EZex5IIUpGTyvFm9nsze8E5lwU0NrMrg5JSJAKc0yyBheP78PyN55Odl89Ns9cxdu56tu/P8DqaRDh/DhnNoeAK5Qt8j/cCfw5aIpEIYGZcek5TFt/Vn99dfjbrdx5m2LRl/PfbX3L0ZLbX8SRC+VMI7Z1zD1NwMhnn3CkKhp+KSAVVi4nmlv4F5xeu7dmSeat3MWDqx8xduZMcnV+QSuZPIWSbWQ3AAZhZe4rMaSQiFdewVjX++tPO/GvSRXRunsB//2MLw6YtY8m2/V5HkwjiTyE8ALwHtDSzBcCHwK+DmkokQp3dtA7zx/Vi1qge5DsYM2c9N81eR8r3x72OJhHA36krGgB9KDhUtMY5FxKXXWqUkVRl2bn5zF+zmyc+SOFEdh7X927F5KHJ1K8Z53U0CXNlnsvIzLqXtEJNfy1SOQ6fyGbaByksWLuHmnHRTBySxE0XtCEuRnfAlfIpTyEUTntdHegBfEbBHsJ5wFrnXL8gZfWbCkEiSer3x/nTP79iWcoB2jasyX2Xd2RIx8aYaYyHlE2Zr0Nwzg1yzg0CdgPdnXM9nHPnA92A7cGLWjpdqSyRKKlJbeaN7cWcMT2JMrh53gZunLWOrd8d8zqaVBH+3CDnU+dc19KWeUF7CBKpcvLyWbBmN49/kMrxzBxG9mrF3Rcn07BWNa+jSRioyB3TFgIngJcoGHp6A1DLOXddMIKWhQpBIt3Rk9k88WEq81fvpkZsNHcOSWTUhW2oFhPtdTQJYRUphOrABKC/b9Ey4DnnXGbAU5aRCkGkwPb9Gfz1X1/x0db9tG4Qz28v68il5zTR+QUplm6hKRIBlqYc4M/vbCF1fwZ92tXn91d24pxmCV7HkhBTnlFGrzrnrjWzzfiuUi7KOXde4GOWjQpB5D/l5uWzcP3XPLZ4G0dP5TCiR0vuviSZxrWrl/5miQjlKYSznHPfmlmxt8t0zu0OcMYyUyGInFn6qRye+jCVuat2US0mitsHJzK2b1uqx+r8QqTTISORCJV2IIO//msrH3z1PS3q1eB3l3fksnOb6vxCBCvzdQhmdtzMjhX5fqzo4+DGFZFAadeoFjNH9WDBzb2pVS2GXy3YyIjn17B5r67jkR/THoJIBMnLd/xt/dc8ungbh09m8/PuLZhyaQea1NH5hUhSkTumjStm2UOBCiYilSc6yvhl71YsmTKQW/q34+1P9zHokY95+qNUMnPyvI4nHvNndqxrzOz6wgdm9izQOHiRRCTY6lSP5beXdeT9u/vTP6kRjyxOYcijS3n7s32E81EDqRh/CuFnwGgzu87M5gHZzrmxQc4lIpWgdYOaTL/xfBaO70NCjVgmLtzENdNX8+nXR72OJh4oadhp/SIPawNvAiuBPwA45w4HO1xpdA5BJHDy8h3/98leHl60jYMZWfysW3OmDOvAWQk1vI4mAVae6xB2UnBBmhX5Xsg559oFI2hZqBBEAi8jK5dnl2xn5oqdRBlMHJLErf3bEx2lYapVRZW6DsHMhgPDExMTx6empnodR6RK+vrwSf7yz69478vv6NmmHo9d25WW9eO9jiUBUJHJ7W4qbrlzbl6AspWb9hBEgss5x5uffsMf3vwSBzwwvBPXnN9CF7WFuXIPOwV6Fvm6CPhv4CcBTSciIcnM+Gm3Frw7+SLOaVaHKa9/zoSXNnL4RLbX0SQIynzIyMwSgPnOOc9LQXsIIpUnL98xc3kajyzeRt34OKZecx4DO2gEejiqyB7C6U4CSRWPJCLhJDrKuHVAe966vR/14mMZPWc9f3jrC05l64K2qiKmtBeY2T/49/TXUUAn4LVghhKR0NWpWR3evqMfUxdtY9aKnazcfpBpI7rRuYXuuxDu/DmpPKDIw1xgt3Nub1BT+UmHjES8tXL7Qe559TMOZmQxeWgSEwYmanhqGCjPbKdnAzjnlgJrnHNLnXMrnXN7zaxPELOKSJjom9iQRZP7c1nns3hkcQrXPr+aPYdOeh1LyqmkcwgvF/l59WnPPRuELCIShhLiY3nqum48MbIrKd8f57InlvHqhq81J1IYKqkQ7Aw/F/dYRCLcVV2b897k/nRukcCvX/+c2176RMNTw0xJheDO8HNxj0VEaF63Bi/f3If7Lu/Ikq0HuHTaMpZs2+91LPFTSaOMWpjZkxTsDRT+jO9x86AnE5GwFBVljO/fjn5JDZn8yqeMmbOeG/u05neXd6RGnO7nHMpKKoQpRX4+fSiPhvaISIk6nlWHt+7oyyOLtjFzxU5W7jjItBFdOa9FXa+jyRmE5eR2hTTsVCQ8rNp+kHte+4wDxwuGp942oD0x0eW5LlYCIZBXKouIlMmFiQ15b1J/LvcNTx0xY42Gp4YgFYKIVIqE+FiePH146noNTw0lKgQRqVSFw1PPa1GXX//f59w6/xMOZWR5HUvwb+qKRsB4oA1FTkKHwn2VdQ5BJHzl5ztmrdjJ1EXbqFMjlqnXnMegszV7amWoyDmEt4AE4APgn0W+RETKrXB46lt39KVhrTjGzF3P/W9u1uypHip1tlMg3jn3X0FPIiIRqeNZdXjz9r48urhgeOqq7Yd4fERXurSs63W0iOPPHsI7ZnZ50JOISMSqHhvNfVd0YsHNvcnMyePnz63iyQ9Tyc3L9zpaRPGnECZRUAqnzOyYmR03s2OBDmJm7cxslpm9Huh1i0h4uLB9Q96d3J8rzjuLx94vmD1196ETXseKGKUWgnOutnMuyjlXwzlXx/e4jj8rN7PZZrbfzL44bfkwM9tmZtvN7De+7aQ558aV79cQkaoioUYsT4zsxpPXdWP7/gwue2I5r6zbo+GplaDU+yGYWffivvxc/1xg2GnrjQaeAS6j4O5r15lZp3KlF5Eq6yddmvHe5P50bVmX3/x9M7doeGrQlXRS+W7gFuDRYp5zwODSVu6cW2ZmbU5b3AvY7pxLAzCzV4CrgC3+BDazW3y5aNWqlT9vEZEw1axuDV4a15vZK3fy8KJtXDptOQ9f05nBZzfxOlqVdMY9BOfcLb7vg4r5KrUMStAc+LrI471AczNrYGbTgW5m9tsScs1wzvVwzvVo1KhRBWKISDiIijJuvqgd/7ijHw1rxTF27gbue2MzJ7NzvY5W5fgz7DTQiru5jnPOHQJuq+wwIhIeOjStzVt39OWxxSnMWJ7G6h0anhpoXkxdsRdoWeRxC2CfBzlEJMxUi4nmt5d35OWb+5CZk8fPNDw1oEosBCvQsqTXlMN6IMnM2ppZHDASeLssKzCz4WY2Iz09PcDRRCQcXNC+Ae9O7s9w3/DUX2h4akCUWAiuYJzXm+VduZktBFYDHcxsr5mNc87lAncAi4CvgFedc1+WZb3OuX84525JSEgobzQRCXMJNWKZNrIbT13XjR0anhoQ/kxu9www1zm3vnIi+U+T24kIwLfpp7jn1c9YteMQQzs24aGfd6ZhrWpexwpZFZncbhCw2sx2mNnnZrbZzD4PfEQRkfI5K6FgeOr9V3RkWeoBhk1bxkdbv/c6VtjxZw+hdXHLnXO7g5LID2Y2HBiemJg4PjU11asYIhKCtn13nEmvbGLrd8f5Ze9W3H9FR+LjvBhQGbrOtIfg1z2VzawLcJHv4XLn3GcBzlcuOmQkIsXJys37YXhqmwY1eXxEV7pqeOoPyn3IyMwmAQuAxr6vl8zszsBHFBEJjMLhqQvH9yE7N5+fP7eKaR+kaHhqKfw5ZPQ5cIFz7oTvcU1gtXPuvErIVyLtIYhIaY5l5vDAW1/yxqZv+EmXZkwb0ZWoqOKuj40cZ9pD8OfAmgFFb2GUR/FXG4uIhJw61WN5fERXkprU4uH3tlG/ZhwPDO+Emf43djp/CmE2sNbM3vA9vhqYFbREfihyUtnLGCISRiYMaM+hjGxmrdhJw1px3DE4yetIIae0K5WjgLXAGOAwcAQY45ybFvxoZ6YL00SkrMyM+y7vyE+7NeeRxSksXLfH60ghp8Q9BOdcvpk96py7ANhYSZlERIIiKsp4+JrzOHIym/ve2Ey9+FiGnXuW17FChj8Xpi02s5+bDriJSBUQGx3Fs9d3p0vLukxc+Cmrdhz0OlLI8KcQ7gZeA7KCeU9lEZHKEh8Xw5zRPWndIJ5b5n3CF99ookzw7xzCMN89lePKek9lEZFQVTc+jnnjelGnegyj56xj10HNllrabKf5wCOVlMVvmv5aRALhrIQazBvXm7x8x02z17H/WKbXkTwVlucQNMpIRAIlsXEt5ozpxcGMLEbNWc+xzByvI3mmLOcQsnUOQUSqoq4t6zL9hvPZvv84N7+4gcycvNLfVAWVWgi+cwZRzrlYnUMQkaqqf3IjHvlFF9bvOszEhZsict4jfya3MzO7wcx+73vc0sx6BT+aiEjluqprcx64shOLt3zPfW98EXF3X/PnkNGzwAXAL32PM4BngpZIRMRDo/u25c7Bifxtw9c8snib13EqlT9zGfV2znU3s00AzrkjZhYX5Fwl0lxGIhJMd1+czMGMbJ5ZsoP6Nasxrl9bryNVCn/2EHLMLBpwAGbWCPD04JpGGYlIMJkZf776XIad05Q/vbOFNzd943WkSuFPITwJvAE0NrO/ACuAvwY1lYiIx6KjjGkju9KnXX3ufe0zPt623+tIQefPKKMFwK+BB4Fvgaudc68FO5iIiNeqx0bzwk096NC0NhNe2sjGPUe8jhRU/uwh4Jzb6px7xjn3tHPuq2CHEhEJFbWrxzJ3TC8a16nG2LnrSf3+uNeRgsavQhARiWSNaldj/tjexEZHcdPsdew7esrrSEGhQhAR8UOrBvG8OKYXGZm53DhrLUdOZHsdKeBUCCIifurUrA4zR/Xg6yOnGDN3PSeycr2OFFBhWQia7VREvNK7XQOeuq4bn+89yoQFG8nOrTpTXIRlIeg6BBHx0qXnNOXBn3VmWcoBprz+Gfn5VWOKC3+uVBYRkdOM6NmKgxnZTF20jXrxcTwwvBMhdJeAclEhiIiU068GtudQRjazV+6kUe1q3D4ovKfTUSGIiJSTmXH/FR05fCKLqYu20aBmHCN7tfI6VrmpEEREKiAqypj6iy4cPZXD797YTN34OIad29TrWOUSlieVRURCSWx0FM9e350uLesy8ZVNrEk75HWkclEhiIgEQHxcDLNH9aRV/XjGv7iBL/eF37B4FYKISIDUqxnHvLG9qF09hlGz17P70AmvI5WJCkFEJICa1a3BvHG9yM3P58ZZ69h/PNPrSH5TIYiIBFhi49rMGd2TA8ezGDV7Pccyc7yO5JewLARNXSEioa5bq3pMv/F8Ur8/zvgXN5CZk+d1pFKFZSFo6goRCQcDkhvx6LVdWLvzMBMXbiI3L7TnPQrLQhARCRdXdW3OA8M7sXjL99z/5hc4F7rzHunCNBGRIBvTty2HMrJ5esl2Gtaqxr2XdvA6UrFUCCIileCeS5I5dCKLp5dsp37NOMb2a+t1pP+gQhARqQRmxp+uOpfDJ7L54ztbaFArjqu6Nvc61o/oHIKISCWJiY7iiZHd6N22Pve8+hkfb9vvdaQfUSGIiFSi6rHRvDCqB0lNajPhpY1s2nPE60g/UCGIiFSyOtVjeXFsTxrVrsaYuevZvv+415EAFYKIiCca167O/HG9iImK4qZZ69h39JTXkVQIIiJead2gJi+O7cnxzFxumr2OIyeyPc2jQhAR8dA5zRJ4YVQP9hw+yZi56zmZnetZFhWCiIjH+rRrwJMju/H53qNMeGkjOR5NcaFCEBEJAcPObcpff9qZpSkHmPLaZ+TnV/4UF2F5YZqZDQeGJyYmeh1FRCRgRvZqxaET2UxdtI16NeP4w5WdMLNK235Y7iFotlMRqap+NbA9Y/q2Yc7KXTz78Y5K3XZY7iGIiFRVZsbvr+jEYd+eQoOacYzs1apStq1CEBEJMVFRxtRrunDkZA6/e2Mz9WrGcek5TYO/3aBvQUREyiwuJorpN3TnvBZ1uXPhJtakHQr6NlUIIiIhKj4uhjmje9KqfjzjX9zAln3Hgro9FYKISAirVzOOeWN7Uat6DDfNXseeQyeDti0VgohIiGtWtwbzxvYiNz+fG2ev5cDxrKBsR4UgIhIGkprUZvbonuw/lsWo2es4lpkT8G2oEEREwkT3VvV47obuHD6RzXfpmQFfv4adioiEkYEdGvPxlIFUj40O+Lq1hyAiEmaCUQagQhARER8VgoiIACoEERHxUSGIiAigQhARER8VgoiIACoEERHxMecq/76dgWJmB4DdFVxNApBeye8v63v8fX1pryvt+YbAwTLkCmUV/fcaKtsMxDqD/RkN1OeztNfo8xm47bZ2zjX6j2edcxH9Bcyo7PeX9T3+vr601/nx/Aav/32Eyr/XUNlmINYZ7M9ooD6fpb1Gn8/gb1eHjOAfHry/rO/x9/Wlva6iv2s48eJ3DcY2A7HOYH9GA/X5LOt2w5lXv2eJ2w3rQ0YSWGa2wTnXw+scIsXR5zP4tIcgRc3wOoBICfT5DDLtIYiICKA9BBER8VEhiIgIoEIQEREfFYL4xcyuNrMXzOwtM7vE6zwiRZlZOzObZWave50lnKkQIoCZzTaz/Wb2xWnLh5nZNjPbbma/KWkdzrk3nXPjgdHAiCDGlQgToM9nmnNuXHCTVn0aZRQBzKw/kAHMc86d61sWDaQAFwN7gfXAdUA08OBpqxjrnNvve9+jwALn3MZKii9VXIA/n687566prOxVTYzXAST4nHPLzKzNaYt7Adudc2kAZvYKcJVz7kHgytPXYWYGPAS8qzKQQArE51MCQ4eMIldz4Osij/f6lp3JncBQ4Bozuy2YwUQo4+fTzBqY2XSgm5n9NtjhqirtIUQuK2bZGY8fOueeBJ4MXhyRHynr5/MQoD9UKkh7CJFrL9CyyOMWwD6PsoicTp9PD6gQItd6IMnM2ppZHDASeNvjTCKF9Pn0gAohApjZQmA10MHM9prZOOdcLnAHsAj4CnjVOfellzklMunzGTo07FRERADtIYiIiI8KQUREABWCiIj4qBBERARQIYiIiI8KQUREABWCCGY22czig7j+ZoGYp9/MMgKRR+RMdB2CRDwz2wX0cM4d9DpLScwswzlXy+scUnVpD0EihpnVNLN/mtlnZvaFmY0ws4lAM2CJmS3xve4SM1ttZhvN7DUzq+VbvsvM/tfM1vm+EovZxgAz+9T3tcnMaptZm8Kbv5jZzCLPHzCzB3zLp5jZejP73Mz+p4Tf4S++/GvMrEkw/jlJ5FIhSCQZBuxzznXx3YjlPd8srvuAQc65QWbWELgfGOqc6w5sAO4uso5jzrlewNPAtGK2cS9wu3OuK3ARcKrok865m33PXQUcAub6bkmaRME9ALoC5/tuGnO6msAa51wXYBkwvuz/CETOTIUgkWQzMNT3V/5Fzrn0Yl7TB+gErDSzT4FRQOsizy8s8v2CYt6/EnjMt+dR1zcnz4+YWXXgNeAO59xu4BLf1yZgI3A2BQVxumzgHd/PnwBtzvyripSd7ocgEcM5l2Jm5wOXAw+a2WLn3B9Pe5kB7zvnrjvTas7wc+E2HjKzf/q2scbMhgKZp71sOvB359wHRbb5oHPu+VJ+hRz375N+eei/Xwkw7SFIxDCzZsBJ59xLwCNAd99Tx4Havp/XAH0Lzw+YWbyZJRdZzYgi31cXs432zrnNzrn/peBw09mnPX87UNs591CRxYuAsUXOVTQ3s8YV+FVFykV/YUgk6QxMNbN8IAeY4Fs+A3jXzL71nUcYDSw0s2q+5++n4IbvANXMbC0Ff0wVtxcx2cwGUfAX/BbgXeCsIs/fC+T4DkcBTHfOTTezjsDqgltXkwHcAOyv6C8sUhYadirip3AZnipSXjpkJCIigPYQRETER3sIIiICqBBERMRHhSAiIoAKQUREfFQIIiICqBBERMTn/wH5IaufEml2JgAAAABJRU5ErkJggg==\n", + "text/plain": [ + "<Figure size 432x288 with 1 Axes>" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" } ], "source": [ - "error\n" + "plt.loglog(grid,error_eukl)\n", + "plt.ylabel('error in Euklidean norm')\n", + "plt.xlabel('step size h')" ] }, { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] - }, - { - "cell_type": "code", - "execution_count": null, + "cell_type": "markdown", "metadata": {}, - "outputs": [], - "source": [] + "source": [ + "Obviously the algorithm is not consistent, although it is supposed to be of order 2.\n", + "Even after several hours of debugging, I couldn't figure out the problem." + ] } ], "metadata": { diff --git a/five_point_stencil.ipynb b/five_point_stencil.ipynb index ef94a3a96985b902cb13f7f4b6009ef5516735df..5eaf5c8bf10f5201e31a311299c2739cf153158e 100644 --- a/five_point_stencil.ipynb +++ b/five_point_stencil.ipynb @@ -7,17 +7,18 @@ "outputs": [], "source": [ "import numpy as np\n", - "from scipy.sparse import csr_matrix\n", - "from scipy.sparse.linalg import spsolve" + "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 is given with Dirichlet boundary condition is given:**\n", + "**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 \\; \\partial \\Omega$$" + "$$ u = g \\quad on \\quad \\partial \\Omega \\quad$$" ] }, { @@ -39,16 +40,16 @@ "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))" + " 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", + "**First, the components of the following equation will be assembled:**\n", "\n", - "$$A \\underline{u} = \\underline{f} + B\\underline{g}$$" + "$$A \\underline{u} = \\underline{f} + B\\underline{g} \\Longleftrightarrow A \\underline{u} - B\\underline{g} = \\underline{f} = -\\Delta u $$" ] }, { @@ -57,14 +58,19 @@ "metadata": {}, "outputs": [], "source": [ + "# matrix A only involves inner nodes. so for m inner nodes, A is an (m x m)-matrix\n", "def matrix_A(h):\n", " N = int(1/h)\n", " m = pow(N-1,2)\n", - " A = pow(h,-2)*(np.zeros((m,m))-4*np.eye(m)+np.eye(m,k=1)+np.eye(m,k=-1)+np.eye(m,k=N-1)+np.eye(m,k=-(N-1)))\n", + " A = pow(h,-2)*(np.zeros((m,m))+4*np.eye(m) - np.eye(m,k=1) - np.eye(m,k=-1) - np.eye(m,k=N-1) - np.eye(m,k=-(N-1)))\n", + " \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", + " #for successors:\n", " A[(i+1)*(N-1)-1][(i+1)*(N-1)] = 0\n", + " #for predecessors\n", " A[(i+1)*(N-1)][(i+1)*(N-1)-1] = 0\n", - " #A = csr_matrix(A)\n", + " \n", " return A" ] }, @@ -72,6 +78,81 @@ "cell_type": "code", "execution_count": 5, "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([[ 64., -16., 0., -16., 0., 0., 0., 0., 0.],\n", + " [-16., 64., -16., 0., -16., 0., 0., 0., 0.],\n", + " [ 0., -16., 64., 0., 0., -16., 0., 0., 0.],\n", + " [-16., 0., 0., 64., -16., 0., -16., 0., 0.],\n", + " [ 0., -16., 0., -16., 64., -16., 0., -16., 0.],\n", + " [ 0., 0., -16., 0., -16., 64., 0., 0., -16.],\n", + " [ 0., 0., 0., -16., 0., 0., 64., -16., 0.],\n", + " [ 0., 0., 0., 0., -16., 0., -16., 64., -16.],\n", + " [ 0., 0., 0., 0., 0., -16., 0., -16., 64.]])" + ] + }, + "execution_count": 5, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# example for matrix A\n", + "n = 2\n", + "h = pow(2,-n)\n", + "A = matrix_A(h)\n", + "A" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [], + "source": [ + "# matrix A in sparse format. this is a lot harder to assemble because of the single deletions\n", + "def sparse_A(h):\n", + " N = int(1/h)\n", + " m = pow(N-1,2)\n", + " A = pow(h,-2)*(4*sp.eye(m) - sp.eye(m,k=1) - sp.eye(m,k=-1) - sp.eye(m,k=N-1) - sp.eye(m,k=-(N-1)))\n", + " \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", + " A.data[(i+1)*(N-1)*5-1*(N-1)-2-1]=0 \n", + " A.data[(i+1)*(N-1)*5-1*(N-1)]=0\n", + " \n", + " return A" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "0.0" + ] + }, + "execution_count": 7, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "#check that sparse matrix is equal to dense matrix\n", + "s=sparse_A(h)\n", + "d=A-s.todense()\n", + "np.linalg.norm(d)" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, "outputs": [], "source": [ "def vector_f(h):\n", @@ -86,33 +167,85 @@ }, { "cell_type": "code", - "execution_count": 6, + "execution_count": 9, "metadata": {}, "outputs": [], "source": [ - "def matrix_B(h):\n", + "# 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", - " B[(i+1)*(N-1)-1][2*(N+i)]= pow(h,-2)\n", - " B[(N-2)*(N-1)+i][N+2+2*i]=pow(h,-2)\n", + " #right boundary:\n", + " B[(i+1)*(N-1)-1][2*(N+i-1)] = pow(h,-2)\n", + " #upper boundary:\n", + " B[-(N-1)+i][-N+i] = pow(h,-2)\n", + " B = sp.csr_matrix(B)\n", " return B" ] }, { "cell_type": "code", - "execution_count": 7, + "execution_count": 10, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "matrix([[ 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n", + " 0., 0., 0.],\n", + " [ 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n", + " 0., 0., 0.],\n", + " [ 0., 0., 0., 0., 0., 0., 16., 0., 0., 0., 0., 0., 0.,\n", + " 0., 0., 0.],\n", + " [ 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n", + " 0., 0., 0.],\n", + " [ 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n", + " 0., 0., 0.],\n", + " [ 0., 0., 0., 0., 0., 0., 0., 0., 16., 0., 0., 0., 0.,\n", + " 0., 0., 0.],\n", + " [ 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 16.,\n", + " 0., 0., 0.],\n", + " [ 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n", + " 16., 0., 0.],\n", + " [ 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 16., 0., 0.,\n", + " 0., 16., 0.]])" + ] + }, + "execution_count": 10, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# example for B\n", + "n = 2\n", + "h = pow(2,-n)\n", + "B = sparse_B(h)\n", + "B.todense()" + ] + }, + { + "cell_type": "code", + "execution_count": 11, "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[3*N+i] = u((i+1)*h,1)\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 " @@ -120,32 +253,88 @@ }, { "cell_type": "code", - "execution_count": 64, + "execution_count": 12, "metadata": {}, - "outputs": [], + "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": 12, + "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": 13, + "metadata": {}, + "outputs": [], + "source": [ + "n = 3\n", + "h = pow(2,-n)\n", "N = pow(2,n)" ] }, { "cell_type": "code", - "execution_count": 65, + "execution_count": 14, "metadata": {}, "outputs": [], "source": [ - "A=matrix_A(h)\n", + "A=sparse_A(h)\n", "f=vector_f(h)\n", - "B=matrix_B(h)\n", - "g=vector_g(h)\n", - "RHS = f+np.matmul(B,g)\n", - "appr_u = np.linalg.solve(A,RHS)" + "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": 66, + "execution_count": 15, + "metadata": {}, + "outputs": [], + "source": [ + "RHS = f+B.dot(g)\n", + "appr_u = sp.linalg.spsolve(A,RHS)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + " " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### **Error and Consistency**" + ] + }, + { + "cell_type": "code", + "execution_count": 16, "metadata": {}, "outputs": [], "source": [ @@ -161,75 +350,120 @@ }, { "cell_type": "code", - "execution_count": 67, + "execution_count": 17, "metadata": {}, - "outputs": [], - "source": [ - "v= exact_solution(h)" - ] - }, - { - "cell_type": "code", - "execution_count": 68, - "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "-8.180678655519896 -8.991060763833046\n", + "-8.270439559649837 -11.480133561857317\n", + "-10.760438683213074 -12.53128157061408\n", + "-12.195353800544211 -12.961737935032156\n", + "-12.821966137009229 -13.145532611303409\n", + "-13.087791494720165 -13.228477586622512\n", + "-13.203877060573376 -13.267565784364171\n" + ] + } + ], "source": [ - "error = max(abs(appr_u - v))" + "grid = np.zeros(7)\n", + "error_eukl = np.zeros(7)\n", + "error_inf = np.zeros(7)\n", + "\n", + "for i in range(7):\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": 69, + "execution_count": 18, "metadata": {}, "outputs": [ { "data": { "text/plain": [ - "7.389346886114219" + "Text(0.5, 0, 'step size h')" ] }, - "execution_count": 69, + "execution_count": 18, "metadata": {}, "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEKCAYAAAD9xUlFAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAAha0lEQVR4nO3de3RU53nv8e8zuqILEtJIgEbchY2xg5FiG9tAEieu69QmNnKzcnFX0ySN18lybyerPW26mp5eVk7S1SanTW+p2+Y0p3Wck8TgW3xJ0ji1jQ3YSMbgS4zAXCQwIAkE6H55zh8zkgWW0Ahpa8+Mfp+19hpmz8zeP5kxj979vvt9zd0REREBiIQdQEREUoeKgoiIjFBREBGRESoKIiIyQkVBRERGqCiIiMiI7LADTEU0GvWlS5eGHUNEJK3s2rWr1d0rxnotrYvC0qVLeemll8KOISKSVszs0Hiv6fKRiIiMUFEQEZERaVkUzGyTmd3X0dERdhQRkYySlkXB3R9193tKSkrCjiIiklHSsiiIiEgwVBRERGREWg9JvVRNJ85ypL2beYW5lBXkUlaUS2FuFmYWdjQRkVDNyqLw8MtH+dufNp23Lzc7QllBLvMKcykvHPWYKBplBbmUFb6zzSvIITtLDS0RySyzsih86salfHBVJe2dfe9sXX20n+vjVFcfbZ19NJ/qor2zjzM9A+MeZ25+NuVFecwryKGsMI+ywnEe1RoRkTQxK4tCtCiPaFFeUu/tHxziVFffeQXkVGe8cJzq7KO9q5/2zl5aTnezp+U07Z199A+OvZpdblYk3soY1Ropu6CgzCvMoTzxOK8glxy1RkRkBs3KojAZOVkRKovzqSzOT+r97s653gFOdfbT1tmbKCjxwnHhY8vpbtrO9U66NTL60lZ5US7z5+azpLyQojz9dYrI1OhfkWlmZhTn51Ccn8Pi8oKkPtM/OMTprv6LXs461dnH0dPd7G3poL2zj77BoXcdJ1qUx9LyApaUF8Yfo4Ujz0vm5Ez3jyoiGUhFIQXkZEWoKM6joji5S1ruTmff4MhlrGOnuznY1sWhtk4OtnWyramVBxt6zvvMvIKcd4pFeSFLo/HHZeWFlBbkqK9DRAAVhbRkZhTlZVOUl82isgLWLip913u6+wY53N7FwbbORLGIF40XD57i4d1H8VHdHnPzs1lSXsiS8gKWDj9G448VRXkqGCKzSFoWBTPbBGyqqakJO0rKmpObxeULirl8QfG7XusdGORIe/d5xeJgWxd7Wjp4Yu/bDA69UzEKcrPOb2GMamnML84nElHBEMkk5j72SJl0cM0117jWU5he/YNDtJzqTrQwus57PNLedd7IqrzsCEsuLBaJlkZV6RyyVDBEUpKZ7XL3a8Z6LS1bChKcnKwIS6OFLI0Wvuu1wSHn6OnuUcXinZbGM2+epHfgnc7v3KwI1WVz3rkcNeoxNm+OhtqKpCgVBUlaVsRYVFbAorICNqyMnvfa0JBz/GwPB1u73nVZavuBNrr6Bs87TvW8OWNelloWLVQLQyREKgoyLSIRY2HJHBaWzOGGFeXnvebunDzXG29htMZHSA0XjcZDpzjb+859GiVzcrhheTnrV0ZZv6KcZdFCdXSLzCAVBQmcmY3cAHjt0rLzXnN32jv7ONjWxYGT59j5Vjvbmlp58tW3AVhYks/6mijra8pZvyJK5dzkbiIUkUujjmZJOe7OwbYunmtq5fmmVl440Mbprn4AVlYWJYpElHXLy5ibr5vyRCbrYh3NKgqS8gaHnNeOnmHb/la2NbXy4sF2evqHiBisqS5lQ02UG2vKqVs8j/ycrLDjiqQ8FQXJKL0DgzQcOs3ziSKxu7mDwSEnLzvCdcvKuHFF/HLTlVUl6rQWGYOKgmS0sz397DjQPtKSePP4OWBUp3VNOetrouq0FknQfQqS0Yrzc7h59XxuXj0fgBNne3hhfxvP7Wvl+f1t53Va37giyoaV6rQWGY9aCpLR3J1Dw53W++NFYrjTuqayKN4fsaKc61eUq9NaZg1dPhJJGBpyXjt2hm1NrWzb38bOt9rO67QeHvpat0Sd1pK5Mq4ojJoQ73P79u0LO46ksd6BQRoPn+b5plaeu6DT+tqlZdxYU86Gmqg6rSWjZFxRGKaWgky3sz39iRvo2tjW1MrPj58F4tOL37CiPDH8NcpydVpLGlNHs0iSivNz+NAV8/nQFfFO65Nne0eGvm5rauOpV48DsGDuqDuta6LMV6e1ZAi1FESS5O4cbh++07qN5/e3cmpUp/X6FeWJO63LtfyppDRdPhIJwHCndbwl0cbOt9rp7h8kYvCe6lI2JqbjeO+SeeRma6pwSR0qCiIzoG9giMbDp9i2P94f8fKR0wwOOXNysli3vIwNNVE2rqzgsvlF6o+QUKkoiITgTE8/2/e38VxiZNOBk50AVBTnsaEmGt9Wqj9CZp6KgkgKaDndzbZ9rTzbFO+4bu/sA+Cy+UVsqKlg48oo1y0rozBP4z8kWCoKIilm9E10zzW1svOtdnoHhsjJMmoXz2NjohWxprpU90fItFNREElxPf2DvHTwFM82nWRbUyt7W84A8fsj4vM1xS83LSkvUH+ETJnuUxBJcfk5WfF/+BNrX7d39sVbEfviLYnhSf2q581h48ooG2oquHFFOfMKc8OMLRlILQWRFOfuvNXaybamVp7d18oL+9s42zuAGbwnVsL6migba6K8d+k88rI1X5NMTJePRDLIwOAQu5s7Eq2IkzQePs3AkJOfE+G6ZeUj/RGrFhTrUpOMSUVBJIOd6x1gx4E2nk1camo6EV9kKFqUy/rE0NeNKytYUKKhrxKnPgWRDFaUl33efE3HOrrZ1tTGc/tO8lxTGw+/fBR4Z/2IjSvjU3EUaeirjEEtBZEM5u688fZZnkvcHzG8fkR2xKhdXMqGmgo2rIxydXUJ2VmaimO20OUjEQHi60fsOnRqZFTTnpYO3KE4P5sblpePDH3VetaZLeOKghbZEZkepzr7eH5kKo6THGnvBiBWOocNNVHWr4yPbNLQ18yScUVhmFoKItPrUFtnvMN6X3xN6zM9A+RkGbdetZC71y1m3bIytSAygIqCiEza4JCzp6WDR14+yg92HeFMzwArKgq5e90S7qqrpqRAa0akKxUFEZmS7r5BfrjnGPfvOETj4dPkZUfYdHUVd69bzNpFpWo9pBkVBRGZNq8e7eA7Ow7zUGMLnX2DrF44l7uvX8wda2Ma5pomplQUzKwU+FVgKaPua3D335q+iJdGRUEkPOd6B3j45Rb+Y/thXj92hsLcLO6sjXH3uiWsrpobdjy5iKkWheeB7cAeYGh4v7t/ezpDXgoVBZHwuTsvHznN/TsO8+juo/QODFG7uJS71y3h9jULyc/RfEypZqpFocHd6wJJNkUqCiKppaOrnwcbmrl/xyH2n+ykZE4Od9VV88l1i6mpLAo7niRMtSj8d+Ac8BjQO7zf3dunM+SlUFEQSU3uzo632rl/x2Ge3HuM/kHn+uVl3L1uCb945QJys3X3dJimWhTuBb4MnAaG3+zuvnw6Q14KFQWR1Nd6rpfvv9TMd3Ye4kh7N9GiXD56zSI+ed1iFpUVhB1vVppqUdgPrHP31iDCTYWKgkj6GBpynm1q5f7th/jJ68dx4H0rK7h73WI+uKpScy/NoKnOkvoq0DW9kURktolEjPdfVsH7L6vgWEc3/+/FI3x35xHu+fddLJibz8evW8THr12sKb5DlkxLYStwJfA05/cpaEiqiEzJwOAQP33jBPfvOMwz+04SMeNDqyq5+/olbKyJEonoprggTLWl8FBiExGZVtlZEW65cgG3XLmAw21dPPDiYb734hF+9NpxFpcV8InrFvPRa6qJFuWFHXXWuGhLwcyygKfc/eaZi5Q8tRREMk/vwCA/evU49+84xPYD7ZqQLwCX3FJw90Ez6zKzEnfvCCaeiMg78rKz2HR1FZuurqLpxFm+s+MIP9h1hEd3H9WEfDMgmT6F7wHXAz8GOof3q09BRGaKJuSbXlMdkvqpsfZrmgsRCYMm5Ju6Kc+Sama5wGWJpz939/5pzHfJVBREZi9NyHfpptpS+ADwbeAgYMAi4FPu/sy0ppwELccpIsM0Id/kTbUo7AI+6e4/Tzy/DHjA3d877UknSS0FERlNE/IlZ6pF4RV3XzPRvjCoKIjIWMabkO+PblvNVbGSsOOFbqpF4VvEJ8L798Suu4Fsd//0tKa8BCoKIjKR4Qn5/uXZA1SXFfDwvevDjhS6qd7R/HngXuC3iPcpPAP8w/TFExEJTrQoj89/YAXZEePLj7/O/pPnWFGhS0njmXBaQnfvdfevu3u9u2929//t7r0TfU5EJJXcsbaKiMHWhpawo6S0CYuCma03sx+b2ZtmdmB4m4lwIiLTpXJuPhtWVrC1sYWhoYmH4s9WyUxg/q/A14ENwLWjNhGRtFJfG6PldDc7D4a+cGTKSqYodLj7E+5+wt3bhrfAk4mITLNbrpxPYW6WLiFdRDJF4Wkz+0szu8HM6oa3wJOJiEyzgtxsbr1qIY/vOUZP/2DYcVJSMqOP1iUeRw9fcuCD0x9HRCRYd9XFeLChmR+/dpxNV1eFHSflTFgU3P2mmQgiIjITrl9ezsKSfLY0NKsojGFSK2Wb2WNBBRERmQmRiHFnbYxn9rVy8qxG119oUkUBiAWSQkRkBtXXxhgcch7ZfTTsKClnskWhMZAUIiIzaOX8Yt4TK2FrY3PYUVLOpIqCu38mqCAiIjNpc22MvS1nePP42bCjpJRk7mi+3cwazazdzM6Y2VkzOzMT4UREgvKRtVVkRYwtumfhPMm0FP4a+BRQ7u5z3b3Y3bWskYiktWhRHu+/rIKHGlsY1LQXI5IpCkeAvZ7Mup0iImmkvi7G22d62H5AkzQMS+bmtf8BPG5m/wWMjN9y968HlkpEZAbcfMV8ivOyebChmfU10bDjpIRkWgpfBrqAfKB41CYiktbyc7K4bc1Cntz7Nl19A2HHSQnJtBTK3P2WwJOIiIRgc22M7754hKdefZvNtdVhxwldMi2Fn5iZioKIZKRrl5ZRPW+ORiElJFMU7gWeNLNuDUkVkUwTiRiba2Nsa2rl+JmesOOELpnlOIvdPeLuczQkVUQy0ebaGEMOD7+s1sKEfQpm9r6x9rv7M9MfR0Rk5i2vKGLtolK2NLRwz/tWhB0nVMl0NP/eqD/nA9cBuwhxPQUz2wRsqqmpCSuCiGSYu+pifOnhV3nt6BlWV83eiyHJXD7aNGr7BeAq4Hjw0S6a6VF3v6ekpCTMGCKSQW5fU0VOlrGlYXZPkjfZWVIBmokXBhGRjDGvMJebLq/k4d1HGRgcCjtOaJLpU/hb4stvQryIrAV2B5hJRCQU9XUxfvTacZ5rauUDl1eGHScUyfQpvDTqzwPAA+6+LaA8IiKhuWlVJSVzctja2KKiMB53//ZMBBERCVtedha3r1nIgw3NnOsdoCgvmd+bM4vWUxARGaW+rpqe/iGe2HMs7Cih0HoKIiKj1C0uZWl5wayd9kLrKYiIjGJmbK6tZvtbbbSc7g47zoxLpigMr6fwRTP7wvAWdDARkbBsro3hDg81zr7WgtZTEBG5wOLyAq5dOo+tjS3MtoskWk9BRGQMm2ur+cOte9jT0sGa6tKw48wYracgIjKG296zkNzsyKzrcNZ6CiIiYygpyOHmKyp5dPdR+mfRtBdaT0FEZBz1tdW0dfbxzJsnw44yY8btUzCzVe7+hpnVjfW6uzcEF0tEJHzvv7yCssJctjS08KEr5ocdZ0ZcrKP5C8A9wNfGeM0JcT0FEZGZkJMV4SNXV/GdnYfp6O6nZE5O2JECN+7lI3e/J/F40xibCoKIzAqba2P0DQzx+CyZ9iKZuY/+3MyyRj2fa2b/J9hYIiKpYU11CSsqCtk6S0YhJTP6KBvYaWZrEkNTXyS+HKeISMYzM+rrqtl5sJ0j7V1hxwlcMqOPvgj8PrAD+DfgNnf/u4BziYikjDtrYwBsnQXTXiRz+eh9wN8Afwb8DPg7M6sKOJeISMqIlc7h+uVlbGlozvhpL5K5fPRXwEfd/Svu/kngPuCnwcYSEUkt9XXVHGzrovHI6bCjBCqZonCDu782/MTdtwDrg4skIpJ6PnzVAvKyI2xpaA47SqCSWY5z0MxuA64kPlPqsD8LLJWISIopzs/hF69cwGOvHONLt68mLztr4g+loWT6FL4JfAz4TcCAjwJLAs4lIpJyNtfFON3Vz9NvZO60F8lcPrrR3X8VOOXufwrcACwKNpaISOrZWBMlWpTH1sbMvYSUTFEYXo+uKzHqqB9YFlwkEZHUlJ0V4Y61Vfz0jROc6uwLO04gkikKj5lZKfCXQANwEPhugJlERFJWfV2M/kHnsQyd9iKZm9f+3N1Pu/uDxPsSVrn7l4KPJiKSelYvnMvl84szdhTShKOPEvMe3QYsHX6/meHuXw82mohI6olPexHjK0+8wVutnSyLFoYdaVolc/noUeDXgHKgeNQmIjIr3bE2hhlszcDWwoQtBaDa3dcEnkREJE0sKMlnQ02ULY0t/M7NlxGJWNiRpk0yLYUnErOjiohIwubaGM2nunnp0Kmwo0yrZIrCdmCrmXWb2RkzO2tmZ4IOJiKSyn7xygUU5GZl3D0LyRSFrxG/Ya3A3ee6e7G7zw04l4hISivMy+bWxLQXPf2DYceZNskUhX3AXk+h+WLNbJOZ3dfR0RF2FBGZxTbXxTjbM8B/vn4i7CjTJpmicAz4mZl90cy+MLwFHexi3P1Rd7+npKQkzBgiMsvduCLK/Ll5GXXPQjJF4S3gP4FcNCRVRGREVsS4c22M/3rzJK3nesOOMy2SmTr7T2ciiIhIOqqvq+afnjnAo7uP8un16T8tXDItBRERGcflC4pZvXBuxqzfrKIgIjJF9XUxXmnuoOnE2bCjTJmKgojIFH1kbRURgy0N6d9aSGZCvArgc4yaEA/A3T8TXCwRkfRRWZzP+y6r4KHGFn73lsvTetqLZFoKDwMlwE+AH47aREQkYXNtjKMdPWx/qy3sKFOSzIR4Be7++4EnERFJY7esXkBRXjZbGlq4cUU07DiXLNmV134p8CQiImlsTm4WH75qAU/sOUZ3X/pOe5FMUfht4oVBE+KJiFxEfV01nX2D/Oi1t8OOcsmSWY6z2N0j7j5HE+KJiIxv3bIyYqVz0noU0rh9Cma2yt3fMLO6sV5394bgYomIpJ9IxLiztop//Nl+TpzpoXJuftiRJu1iHc1fAO4hPnX2hRz4YCCJRETS2Obaav7+6f08svsov75xedhxJm3couDu9yQeb5q5OCIi6a2msoirq0t4sKElLYuC7mgWEZlmm2tjvH7sDK8fS78xOSoKIiLTbNPVVWRHLC0nybtoUbC4RTMVRkQkE5QX5fGBy+PTXgwOpcyilUm5aFFILMH50MxEERHJHPV11Zw428u2ptawo0xKMpePtpvZtYEnERHJIB9cVUlxfnbaXUJKpijcBLxgZvvN7BUz22NmrwQdTEQkneXnZHH7miqe3Ps2nb0DYcdJWjIT4n048BQiIhmovi7GAzsP8+Tet7nrvdVhx0lKMtNcHAJKgU2JrTSxT0RELuKaJfNYVDaHLY3NYUdJ2oRFwcx+G7gfqExs/2Fmvxl0MBGRdGdmbK6t5vn9bRzr6A47TlKS6VP4LLDO3f/Y3f8YuJ74SmwiIjKB+toY7vBQ49GwoyQlmaJgwOjJwQcT+0REZAJLo4XULS5lS0Mz8VH+qS2ZovAtYIeZ/YmZ/QmwHfjXQFOJiGSQ+rpq9p04x6tHU3/ai4nuaI4AO4BPA+3AKeDT7v7XwUcTEckMt69ZSG5WJC3WWZjojuYh4Gvu3uDu33D3v3H3xhnKJiKSEUoLcvngqkoe2d3CwOBQ2HEuKpnLRz8ys7vMTP0IIiKXaHNdjNZzfTy7L7WnvUimKHwB+D7QqzWaRUQuzU2XV1JakMODDal9z0IyfQq3JtZoztUazSIilyY3O8KmNVX8+LXjnOnpDzvOuJLpU/irGcoiIpLR6uti9A4M8cSeY2FHGZf6FEREZsjaRaUsixam9CikyfQp9KlPQUTk0pkZ9bUxdrzVzpH2rrDjjCmZCfGKE30KOepTEBGZmjtrYwA8/HJqthaSmRDPzOxXzOxLieeLzOy64KOJiGSeRWUFXLesjC0NLSk57UUyl4/+AbgB+GTi+Tng7wNLJCKS4eprYxxo7WR3c0fYUd4lmaKwzt3vBXoA3P0UkBtoKhGRDPZLaxaSmx1hSwres5BMUeg3syzAAcysAkjt+7RFRFLY3PwcfmH1fB7dfZS+gdT65zSZovANYCtQaWZfBp4D/legqUREMtxddTFOdfXzs5+fCDvKeSZco9nd7zezXcCHiK+jcKe7vx54MhGRDLZxZQXlhblsbWzhlisXhB1nxIRFAcDd3wDeCDiLiMiskZMV4SNrq7h/+2E6uvopKcgJOxKQ3OUjEREJQH1tNX2DQzy2J3WW6lRREBEJyVWxuaysLEqpaS9UFEREQmJmbK6LsevQKQ61dYYdB1BREBEJ1Z1rY5jB1sbUaC2oKIiIhKiqdA43LC9na2NqTHuhoiAiErL6umoOtXXRcPhU2FFUFEREwnbrVQvIz4nwYAp0OKsoiIiErCgvm1uvXMBju4/SOzAYahYVBRGRFLC5rpozPQP89PVwp71QURARSQHrV5RTWZzHlpBHIakoiIikgOysCHesreLpN07Q3tkXWg4VBRGRFFFfV83AkPPYK+FNe6GiICKSIq5YOJdVC4pDHYWkoiAikkLuqqtm95HT7D95LpTzp0xRMLNCM/u2mf2zmd0ddh4RkTDcsbaKiMHWkFoLgRYFM/uWmZ0ws70X7L/VzH5uZk1m9geJ3fXAD9z9c8BHgswlIpKqKufms2FlBVsbWxgamvlpL4JuKfwbcOvoHYn1nv8e+DCwGviEma0GqoEjibeFe/eGiEiI6mtjtJzuZufB9hk/d6BFwd2fAS78qa4Dmtz9gLv3Ad8F7gCaiReGi+Yys3vM7CUze+nkyZNBxBYRCdUtV86nMDcrlEtIYfQpxHinRQDxYhADtgB3mdk/Ao+O92F3v8/dr3H3ayoqKoJNKiISgoLcbG69aiGP7zlGT//MXjgJoyjYGPvc3Tvd/dPu/nl3v3/GU4mIpJC76mKc7R3gx68dn9HzhlEUmoFFo55XA6mzQKmISAq4fnk5C0vy2dLQPKPnDaMovAisNLNlZpYLfBx4JIQcIiIpKxIx7qyN8cy+Vk6e7Z258wZ5cDN7AHgBuNzMms3ss+4+APwG8BTwOvA9d381yBwiIumovjbG4JDzyO6Zu5iSHeTB3f0T4+x/HHg8yHOLiKS7lfOLeU+shK2NzXx2w7IZOWfK3NEsIiLvtrk2xt6WM7x5/OyMnE9FQUQkhX1kbRVZEWPLDN2zoKIgIpLCokV5vP+yCh5qbGFwBqa9SMuiYGabzOy+jo6OsKOIiASuvi7G22d62H6gLfBzpWVRcPdH3f2ekpKSsKOIiATu5ivmU5yXzYMzcM9CWhYFEZHZJD8ni9vWLOTJvW/T1TcQ6LlUFERE0sDm2hhdfYM89erbgZ5HRUFEJA1cu7SM6nlzAh+FpKIgIpIGIhFjc22MbU2tHD/TE9x5AjuyiIhMq821MYYcHn45uNaCioKISJpYXlHE2kWlgV5CUlEQEUkjd9XFeOPts7x29Ewgx1dREBFJI7evqSJalMehts5Ajh/oLKlBMbNNwKaampqwo4iIzKh5hbns+MMPkRUZaxHLqUvLloLuaBaR2SyoggBpWhRERCQYKgoiIjJCRUFEREaoKIiIyAgVBRERGaGiICIiI1QURERkhLkHv+ZnUMzsJHBoiocpAaayruelfn4yn0v2vRO9b6LXo0BrkplS3VT/XlPlnNNxzEs5Rhjfz2Tekynf0TC+n6PPu8TdK8Z8h7vP6g24L4zPT+Zzyb53ovcl8fpLYf99pMrfa6qcczqOeSnHCOP7mcx7MuU7Gsb3M9nz6vIRPBrS5yfzuWTfO9H7pvqzppMwftYgzjkdx7yUY4Tx/ZzsedNZWD/nhOdN68tHMr3M7CV3vybsHCLj0Xc0eGopyGj3hR1AZAL6jgZMLQURERmhloKIiIxQURARkREqCiIiMkJFQZJiZnea2T+b2cNmdkvYeURGM7PlZvavZvaDsLOkOxWFWcDMvmVmJ8xs7wX7bzWzn5tZk5n9wcWO4e4PufvngF8DPhZgXJllpun7ecDdPxts0tlBo49mATN7H3AO+L/uflViXxbwJvALQDPwIvAJIAv4ygWH+Iy7n0h87mvA/e7eMEPxJcNN8/fzB+7+yzOVPRNlhx1Agufuz5jZ0gt2Xwc0ufsBADP7LnCHu38FuP3CY5iZAV8FnlBBkOk0Hd9PmT66fDR7xYAjo543J/aN5zeBm4FfNrP/FmQwESb5/TSzcjP7JlBrZl8MOlwmU0th9rIx9o17LdHdvwF8I7g4IueZ7PezDdAvK9NALYXZqxlYNOp5NXA0pCwiF9L3MyQqCrPXi8BKM1tmZrnAx4FHQs4kMkzfz5CoKMwCZvYA8AJwuZk1m9ln3X0A+A3gKeB14Hvu/mqYOWV20vcztWhIqoiIjFBLQURERqgoiIjICBUFEREZoaIgIiIjVBRERGSEioKIiIxQURABzOx3zKwgwONXTcdc/2Z2bjryiIxH9ymIAGZ2ELjG3VvDznIxZnbO3YvCziGZSy0FmVXMrNDMfmhmu81sr5l9zMx+C6gCnjazpxPvu8XMXjCzBjP7vpkVJfYfNLO/MLOdia1mjHO838xeTmyNZlZsZkuHF5Exs38Z9fpJM/ufif2/Z2YvmtkrZvanF/kZvpzIv93M5gfx30lmLxUFmW1uBY66+9WJBV2eTMwAexS4yd1vMrMo8EfAze5eB7wEfGHUMc64+3XA3wF/PcY5fhe4193XAhuB7tEvuvuvJ167A2gD/i2xxOlK4usIrAXem1h85kKFwHZ3vxp4Bvjc5P8TiIxPRUFmmz3AzYnf9je6e8cY77keWA1sM7OXgU8BS0a9/sCoxxvG+Pw24OuJFkhpYh6f85hZPvB94Dfc/RBwS2JrBBqAVcSLxIX6gMcSf94FLB3/RxWZPK2nILOKu79pZu8Ffgn4ipn9yN3/7IK3GfBjd//EeIcZ58/D5/iqmf0wcY7tZnYz0HPB274JbHH3n4w651fc/Z8m+BH6/Z2OwEH0/7BMM7UUZFYxsyqgy93/A/groC7x0lmgOPHn7cD64f4CMysws8tGHeZjox5fGOMcK9x9j7v/BfFLT6sueP1eoNjdvzpq91PAZ0b1XcTMrHIKP6rIJdFvGTLbvAf4SzMbAvqBzyf23wc8YWbHEv0KvwY8YGZ5idf/iPhC8gB5ZraD+C9VY7UmfsfMbiL+m/xrwBPAwlGv/y7Qn7g0BfBNd/+mmV0BvBBfDptzwK8AJ6b6A4tMhoakikxCugxdFblUunwkIiIj1FIQEZERaimIiMgIFQURERmhoiAiIiNUFEREZISKgoiIjFBREBGREf8f6yXLYPU7gx4AAAAASUVORK5CYII=\n", + "text/plain": [ + "<Figure size 432x288 with 1 Axes>" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" } ], "source": [ - "appr_u[-1]" + "plt.loglog(grid,error_inf)\n", + "plt.ylabel('error in maximum-norm')\n", + "plt.xlabel('step size h')" ] }, { "cell_type": "code", - "execution_count": 70, + "execution_count": 19, "metadata": {}, "outputs": [ { "data": { "text/plain": [ - "16.380407649947266" + "Text(0.5, 0, 'step size h')" ] }, - "execution_count": 70, + "execution_count": 19, "metadata": {}, "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAEKCAYAAAASByJ7AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAAlFUlEQVR4nO3dd3xUdb7G8c83DQgldJBekiAoUqQpSFexsLq7ruBaaKKyKmBh7+7qrvdu0ysW7Ig0QcRV71rWXQUL0qugoggJhCKi0gMB0n/3j0zcyIZkkszkzGSe9+uVVzJnZs55oqNPzjm/8zvmnENERCTK6wAiIhIaVAgiIgKoEERExEeFICIigApBRER8VAgiIgJAjNcBysPMhgPDa9euPT45OdnrOCIiYeWTTz456JxrdPpyC+frEHr06OE2bNjgdQwRkbBiZp8453qcvlyHjEREBFAhiIiIjwpBREQAFYKIiPioEEREBFAhiIiIT1gWgpkNN7MZ6enp5Xr/J7sPs/W7YwFOJSIS3sKyEJxz/3DO3ZKQkFCu9z/4r60Mm7acCS99wlffqhhERCBMC6GiZo7qwcTBiSxPPchlTyzntvmfsGWfikFEIltEX6l89GQ2s1fsZM7KXRzPymXYOU2ZOCSJTs3qBDCliEhoOdOVyhFdCIXST+Ywa+VO5qzYyfGsXC49pwkThyRxTrPyHZISEQllKgQ/pJ/MYfbKncxeuZPjmblc0qkJk4aqGESkalEhlEH6qRzmrNzJrBUFxXBxpyZMGpLEuc1VDCIS/lQI5ZB+Koe5K3cxc0UaxzNzGdqxCZOHqhhEJLypECqgsBhmrUjjmIpBRMKcCiEAjmX69hiWFxZDYyYNSaZzCxWDiIQPFUIAHcvM4cWVu5i5Yifpp3IYcnZjJg1N4rwWdSs9i4hIWVWpQii8hWZiYuL41NRUz3Icz8zhxVW7eGF5QTEMPrsxk4Yk0aVlXc8yiYiUpkoVQqFQuYXm8cwc5q3ezQvL0zh6ModBHRoxaWgyXVUMIhKCVAiVQMUgIuFAhVCJMrJyeXFVwcnnIydzGNihEZOGJNGtVT2vo4mIqBC8kJGVy7zVu3hhWUExDEhuxKShSXRXMYiIh1QIHsrIymX+6t3MWLaDIydz6J9csMdwfmsVg4hUPhVCCDiRlcv8NbuZsSyNwyeyuSipIZOHJqsYRKRSqRBCSPHFkMT5ret7HU1EIoAKIQSdyMrlJV8xHPIVw6QhSfRoo2IQkeBRIYSwk9kFxfD80oJi6JdYsMegYhCRYFAhhIGT2bksWLOH55ft4GBGNn0TGzB5aDI9VQwiEkAqhDByKjuPBWt3M31pQTFc2L6gGHq1VTGISMWpEMLQv4shjYMZWVzYvgGThiTRu10Dr6OJSBhTIYSx04vhgnYNmDQ0iT4qBhEphypVCKEy22llO5Wdx8vr9jB96Q4OHM+iT7v6TB6arGIQkTKpUoVQKFL2EE6XmZPHy2v38FyRYrj74g46xyAiflEhVEGFxfDsxzs4mJFFv8SG3HWxrnwWkZKpEKqwwnMMz328g0MnshmQ3Ii7Lta02yJSPBVCBDiZncu81bt5fmnBJHpDzm7MXRcnc25z3fNZRP5NhRBBCu/HMGNZGumncrikUxMmD02mU7M6XkcTkRCgQohAxzJzmLNiFzNXpHE8M5fLOzdl0pBkOjSt7XU0EfGQCiGCpZ/MYdaKNGav3MWJ7FyuPK8Zk4YkkthYxSASiVQIwpET2bywPI25q3ZxKiePq7o0Y+KQJNo1quV1NBGpRCoE+cGhjCxmLE9j3qrdZOXm8bPuLZg4OIlWDeK9jiYilUCFIP/hwPEsnl+6g/lrdpOb77imewvuGJxIy/oqBpGqTIUgZ7T/WCbPfryDl9ftIT/fcW3PltwxKJFmdWt4HU1EgkCFIKX6Nv0Uzy7ZwSvr92AYI3u15FcDE2maUN3raCISQCoE8ds3R0/x9EfbeW3D10RFGdf3bsWEge1pXFvFIFIVqBCkzL4+fJKnP9rO6xv3Ehtt3NinNbcOaE/DWtW8jiYiFaBCkHLbdfAET36UypubvqFaTDSjLmzDLf3bUb9mnNfRRKQcVAhSYTsOZPDUh6m89dk+4mOjGd23DeMvakfdeBWDSDhRIUjApH5/nCc+TOWdz7+ldrUYxvRry7h+bUmoEet1NBHxQ4ULwczqADGFj51zhwMXr3xUCN7a+t0xnvgglXe/+I7a1WMYf1E7xvRtQ+3qKgaRUFbuQjCzW4E/AqeAwhc751y7gKcsIxVCaPhyXzrTPkjl/S3fUzc+lvEXtWP0hW2oWS2m9DeLSKWrSCGkAhc45w4GK1x5qRBCy+a96Tz+QQofbd1P/Zpx3Nq/HTde0Jr4OBWDSCg5UyFE+fHeHcDJwEeSqqZziwRmj+7Jm7f3pXPzBB58dyv9H17CzOVpZObkeR1PRErhzx5CN2AOsBbIKlzunJsY0CBmVwNXAI2BZ5xzi0t7j/YQQtsnuw/z+PuprNh+kEa1q3H7wPaM7NWK6rHRXkcTiWgVOWS0DlgBbAbyC5c75170Y6OzgSuB/c65c4ssHwY8AUQDM51zDxV5rh7wiHNuXGnrVyGEh7Vph3js/RTW7jxM0zrVuX1Qe67t2ZJqMSoGES9UpBBWOecuLOdG+wMZwLzCQjCzaCAFuBjYC6wHrnPObfE9/yiwwDm3sbT1qxDCy6odB3lscQobdh+hWUJ17hicxDXntyAuxp8jlyISKBU5h7DEzG4xs7PMrH7hlz8bdc4tA04fntoL2O6cS3POZQOvAFdZgf8F3i2pDHxZNpjZhgMHDvgTQ0LEhe0b8tptFzB/XC+aJFTnd29sZvCjH/Pq+q/JycsvfQUiElT+7CHsLGax38NOzawN8E6RPYRrgGHOuZt9j28EelOw1zCKgj2GT51z00tbt/YQwpdzjo9TDvD4+yl8vjed1g3imTg4iau7NSc6yryOJ1KlnWkPocTxgGYWBfzGOfe3QGYpZplzzj0JPBnA7UgIMzMGdWjMwORGfPjVfh7/IIV7XvuM55bu4J6Lkxl2blPMVAwilanEQ0bOuXzg9gBvcy/QssjjFsC+AG9DwoSZMbRTE965sx/PXd8dgAkLNnLVMytZnnqAcJ5aRSTc+HMO4X0zu9fMWpb1HMIZrAeSzKytmcUBI4G3y7ICMxtuZjPS09MrEENCiZlxWeezWDS5P1OvOY9DGdncOGsd172who17jngdTyQiBPUcgpktBAYCDYHvgQecc7PM7HJgGgXDTmc75/5SxtyAziFUZVm5eSxcu4enl2znYEY2Qzs2YcqlHejQtLbX0UTCnmY7lbB0IiuXOSt38vyyNDKycrm6a3PuGppMqwbxXkcTCVsVuQ4hFpgA9Pct+hh43jmXE+iQZaVCiBxHT2YzfWkac1ftJDfPMbJXSyYOTqJxHd3WU6SsKlIIM4FYoPDK5BuBvMJho14ws+HA8MTExPGpqalexRAPfH8sk6c+SuWVdV8TE22MvrAttw3QTXpEyqIihfCZc65Lacu8oD2EyLXn0Eke/yCFNz/9hlrVYri1fzvG9G2rKbdF/FCRK5XzzKx9kRW1AzR1pXiqVYN4Hh/RlXcnXUSfdg14ZHEKA6YuYc7KnWTl6uMpUh7+7CEMoWC20zQKLiprDYxxzi0JfrySaQ9BCm3cc4Sp721jddohmtetweShSfysewtd9SxSjAqNMjKzakAHCgphq3Muq5S3VAoVghTlnGPF9oNMXbSNz/emk9i4lq56FilGRQvhQqANP76n8rxABiwLnVSWkjjnWPTldzyyOIXt+zM4r0UCUy7tQL/EhioGESp2Unk+0B74lH+fO3CBvkFOeWgPQUqSl+/4+8a9TPsglW+OnqJPu/r8etjZdG9Vz+toIp6qSCF8BXRyIXgFmwpB/FHcVc/3XprM2U3reB1NxBMVGWX0BdA08JFEKke1mGhG923L0imDuPeSZNbuPMRlTyxn8iub2H3ohNfxREKGP3sIS4CuwDp+fE/lnwQ1mR+0hyDlcfpVzyN6tmTikCSa6KpniRAVOWQ0oLjlzrmlAcpWZjqpLIGw/1gmTxa56nnUhW2YMKC9rnqWKk+T24mcga56lkijQhApxdbvjvHo4hTe3/I9DWvFcfugRH7ZuxXVYqK9jiYSUCoEET8Vd9XzT7s1JybanzEYIqGvIqOMRCJK91b1eHl8b14a15sGteKY8vrnXDptGe9u/la39JQqrdRCMLO+Zva+maWYWZqZ7TSztMoIJ+IVM6NfUkPeur0v02/ojpkxYcFGfvL0Spal6F7PUjX5M8poK3AX8AlFZjl1zh0KbrTS6ZCRVJa8fMcbm77h8fdTdNWzhL2KDDtd65zrHbRk5aBhp+IVXfUsVUFFCuEhIBr4Oz++MG1joEOWlfYQxCsnsnKZu2oX05fuICMrl6u6NOOui5Np3aCm19FESlWRQijuvgfOOTc4UOHKS4UgXtNVzxKONOxUJIj2H8vkqY+2s3Ddnh/u9fyrQe2pUz3W62gi/6Gi90O4AjgH+OHPHufcHwOasBxUCBJqil71XD8+jrsvSWZkz1a6c5uElHJfh2Bm04ERwJ0U3DHtFxTcRlNETlN4r+e3b+9H+0a1uO+NL7jiyeWs3H7Q62gipfLnwrQLnXM3AUecc/8DXAC0DG4skfDWuUUCf7u1D89e352MrFyun7mWm1/cwM6Dmm5bQpc/hXDK9/2kmTUDcoC2wYskUjWYGZd3PosP7h7Afw07m9U7DnLJ40v58ztbSD+V43U8kf/gTyG8Y2Z1ganARmAX8EoQM5XKzIab2Yz09HQvY4j4pXpsNBMGtmfJlIH8vHsLZq3cycCpS5i/ehe5eflexxP5QZlGGZlZNaC6cy4k/k+sk8oSjr7cl86f3tnCmrTDJDWuxe+v7ET/5EZex5IIUpGTyvFm9nsze8E5lwU0NrMrg5JSJAKc0yyBheP78PyN55Odl89Ns9cxdu56tu/P8DqaRDh/DhnNoeAK5Qt8j/cCfw5aIpEIYGZcek5TFt/Vn99dfjbrdx5m2LRl/PfbX3L0ZLbX8SRC+VMI7Z1zD1NwMhnn3CkKhp+KSAVVi4nmlv4F5xeu7dmSeat3MWDqx8xduZMcnV+QSuZPIWSbWQ3AAZhZe4rMaSQiFdewVjX++tPO/GvSRXRunsB//2MLw6YtY8m2/V5HkwjiTyE8ALwHtDSzBcCHwK+DmkokQp3dtA7zx/Vi1qge5DsYM2c9N81eR8r3x72OJhHA36krGgB9KDhUtMY5FxKXXWqUkVRl2bn5zF+zmyc+SOFEdh7X927F5KHJ1K8Z53U0CXNlnsvIzLqXtEJNfy1SOQ6fyGbaByksWLuHmnHRTBySxE0XtCEuRnfAlfIpTyEUTntdHegBfEbBHsJ5wFrnXL8gZfWbCkEiSer3x/nTP79iWcoB2jasyX2Xd2RIx8aYaYyHlE2Zr0Nwzg1yzg0CdgPdnXM9nHPnA92A7cGLWjpdqSyRKKlJbeaN7cWcMT2JMrh53gZunLWOrd8d8zqaVBH+3CDnU+dc19KWeUF7CBKpcvLyWbBmN49/kMrxzBxG9mrF3Rcn07BWNa+jSRioyB3TFgIngJcoGHp6A1DLOXddMIKWhQpBIt3Rk9k88WEq81fvpkZsNHcOSWTUhW2oFhPtdTQJYRUphOrABKC/b9Ey4DnnXGbAU5aRCkGkwPb9Gfz1X1/x0db9tG4Qz28v68il5zTR+QUplm6hKRIBlqYc4M/vbCF1fwZ92tXn91d24pxmCV7HkhBTnlFGrzrnrjWzzfiuUi7KOXde4GOWjQpB5D/l5uWzcP3XPLZ4G0dP5TCiR0vuviSZxrWrl/5miQjlKYSznHPfmlmxt8t0zu0OcMYyUyGInFn6qRye+jCVuat2US0mitsHJzK2b1uqx+r8QqTTISORCJV2IIO//msrH3z1PS3q1eB3l3fksnOb6vxCBCvzdQhmdtzMjhX5fqzo4+DGFZFAadeoFjNH9WDBzb2pVS2GXy3YyIjn17B5r67jkR/THoJIBMnLd/xt/dc8ungbh09m8/PuLZhyaQea1NH5hUhSkTumjStm2UOBCiYilSc6yvhl71YsmTKQW/q34+1P9zHokY95+qNUMnPyvI4nHvNndqxrzOz6wgdm9izQOHiRRCTY6lSP5beXdeT9u/vTP6kRjyxOYcijS3n7s32E81EDqRh/CuFnwGgzu87M5gHZzrmxQc4lIpWgdYOaTL/xfBaO70NCjVgmLtzENdNX8+nXR72OJh4oadhp/SIPawNvAiuBPwA45w4HO1xpdA5BJHDy8h3/98leHl60jYMZWfysW3OmDOvAWQk1vI4mAVae6xB2UnBBmhX5Xsg559oFI2hZqBBEAi8jK5dnl2xn5oqdRBlMHJLErf3bEx2lYapVRZW6DsHMhgPDExMTx6empnodR6RK+vrwSf7yz69478vv6NmmHo9d25WW9eO9jiUBUJHJ7W4qbrlzbl6AspWb9hBEgss5x5uffsMf3vwSBzwwvBPXnN9CF7WFuXIPOwV6Fvm6CPhv4CcBTSciIcnM+Gm3Frw7+SLOaVaHKa9/zoSXNnL4RLbX0SQIynzIyMwSgPnOOc9LQXsIIpUnL98xc3kajyzeRt34OKZecx4DO2gEejiqyB7C6U4CSRWPJCLhJDrKuHVAe966vR/14mMZPWc9f3jrC05l64K2qiKmtBeY2T/49/TXUUAn4LVghhKR0NWpWR3evqMfUxdtY9aKnazcfpBpI7rRuYXuuxDu/DmpPKDIw1xgt3Nub1BT+UmHjES8tXL7Qe559TMOZmQxeWgSEwYmanhqGCjPbKdnAzjnlgJrnHNLnXMrnXN7zaxPELOKSJjom9iQRZP7c1nns3hkcQrXPr+aPYdOeh1LyqmkcwgvF/l59WnPPRuELCIShhLiY3nqum48MbIrKd8f57InlvHqhq81J1IYKqkQ7Aw/F/dYRCLcVV2b897k/nRukcCvX/+c2176RMNTw0xJheDO8HNxj0VEaF63Bi/f3If7Lu/Ikq0HuHTaMpZs2+91LPFTSaOMWpjZkxTsDRT+jO9x86AnE5GwFBVljO/fjn5JDZn8yqeMmbOeG/u05neXd6RGnO7nHMpKKoQpRX4+fSiPhvaISIk6nlWHt+7oyyOLtjFzxU5W7jjItBFdOa9FXa+jyRmE5eR2hTTsVCQ8rNp+kHte+4wDxwuGp942oD0x0eW5LlYCIZBXKouIlMmFiQ15b1J/LvcNTx0xY42Gp4YgFYKIVIqE+FiePH146noNTw0lKgQRqVSFw1PPa1GXX//f59w6/xMOZWR5HUvwb+qKRsB4oA1FTkKHwn2VdQ5BJHzl5ztmrdjJ1EXbqFMjlqnXnMegszV7amWoyDmEt4AE4APgn0W+RETKrXB46lt39KVhrTjGzF3P/W9u1uypHip1tlMg3jn3X0FPIiIRqeNZdXjz9r48urhgeOqq7Yd4fERXurSs63W0iOPPHsI7ZnZ50JOISMSqHhvNfVd0YsHNvcnMyePnz63iyQ9Tyc3L9zpaRPGnECZRUAqnzOyYmR03s2OBDmJm7cxslpm9Huh1i0h4uLB9Q96d3J8rzjuLx94vmD1196ETXseKGKUWgnOutnMuyjlXwzlXx/e4jj8rN7PZZrbfzL44bfkwM9tmZtvN7De+7aQ558aV79cQkaoioUYsT4zsxpPXdWP7/gwue2I5r6zbo+GplaDU+yGYWffivvxc/1xg2GnrjQaeAS6j4O5r15lZp3KlF5Eq6yddmvHe5P50bVmX3/x9M7doeGrQlXRS+W7gFuDRYp5zwODSVu6cW2ZmbU5b3AvY7pxLAzCzV4CrgC3+BDazW3y5aNWqlT9vEZEw1axuDV4a15vZK3fy8KJtXDptOQ9f05nBZzfxOlqVdMY9BOfcLb7vg4r5KrUMStAc+LrI471AczNrYGbTgW5m9tsScs1wzvVwzvVo1KhRBWKISDiIijJuvqgd/7ijHw1rxTF27gbue2MzJ7NzvY5W5fgz7DTQiru5jnPOHQJuq+wwIhIeOjStzVt39OWxxSnMWJ7G6h0anhpoXkxdsRdoWeRxC2CfBzlEJMxUi4nmt5d35OWb+5CZk8fPNDw1oEosBCvQsqTXlMN6IMnM2ppZHDASeLssKzCz4WY2Iz09PcDRRCQcXNC+Ae9O7s9w3/DUX2h4akCUWAiuYJzXm+VduZktBFYDHcxsr5mNc87lAncAi4CvgFedc1+WZb3OuX84525JSEgobzQRCXMJNWKZNrIbT13XjR0anhoQ/kxu9www1zm3vnIi+U+T24kIwLfpp7jn1c9YteMQQzs24aGfd6ZhrWpexwpZFZncbhCw2sx2mNnnZrbZzD4PfEQRkfI5K6FgeOr9V3RkWeoBhk1bxkdbv/c6VtjxZw+hdXHLnXO7g5LID2Y2HBiemJg4PjU11asYIhKCtn13nEmvbGLrd8f5Ze9W3H9FR+LjvBhQGbrOtIfg1z2VzawLcJHv4XLn3GcBzlcuOmQkIsXJys37YXhqmwY1eXxEV7pqeOoPyn3IyMwmAQuAxr6vl8zszsBHFBEJjMLhqQvH9yE7N5+fP7eKaR+kaHhqKfw5ZPQ5cIFz7oTvcU1gtXPuvErIVyLtIYhIaY5l5vDAW1/yxqZv+EmXZkwb0ZWoqOKuj40cZ9pD8OfAmgFFb2GUR/FXG4uIhJw61WN5fERXkprU4uH3tlG/ZhwPDO+Emf43djp/CmE2sNbM3vA9vhqYFbREfihyUtnLGCISRiYMaM+hjGxmrdhJw1px3DE4yetIIae0K5WjgLXAGOAwcAQY45ybFvxoZ6YL00SkrMyM+y7vyE+7NeeRxSksXLfH60ghp8Q9BOdcvpk96py7ANhYSZlERIIiKsp4+JrzOHIym/ve2Ey9+FiGnXuW17FChj8Xpi02s5+bDriJSBUQGx3Fs9d3p0vLukxc+Cmrdhz0OlLI8KcQ7gZeA7KCeU9lEZHKEh8Xw5zRPWndIJ5b5n3CF99ookzw7xzCMN89lePKek9lEZFQVTc+jnnjelGnegyj56xj10HNllrabKf5wCOVlMVvmv5aRALhrIQazBvXm7x8x02z17H/WKbXkTwVlucQNMpIRAIlsXEt5ozpxcGMLEbNWc+xzByvI3mmLOcQsnUOQUSqoq4t6zL9hvPZvv84N7+4gcycvNLfVAWVWgi+cwZRzrlYnUMQkaqqf3IjHvlFF9bvOszEhZsict4jfya3MzO7wcx+73vc0sx6BT+aiEjluqprcx64shOLt3zPfW98EXF3X/PnkNGzwAXAL32PM4BngpZIRMRDo/u25c7Bifxtw9c8snib13EqlT9zGfV2znU3s00AzrkjZhYX5Fwl0lxGIhJMd1+czMGMbJ5ZsoP6Nasxrl9bryNVCn/2EHLMLBpwAGbWCPD04JpGGYlIMJkZf776XIad05Q/vbOFNzd943WkSuFPITwJvAE0NrO/ACuAvwY1lYiIx6KjjGkju9KnXX3ufe0zPt623+tIQefPKKMFwK+BB4Fvgaudc68FO5iIiNeqx0bzwk096NC0NhNe2sjGPUe8jhRU/uwh4Jzb6px7xjn3tHPuq2CHEhEJFbWrxzJ3TC8a16nG2LnrSf3+uNeRgsavQhARiWSNaldj/tjexEZHcdPsdew7esrrSEGhQhAR8UOrBvG8OKYXGZm53DhrLUdOZHsdKeBUCCIifurUrA4zR/Xg6yOnGDN3PSeycr2OFFBhWQia7VREvNK7XQOeuq4bn+89yoQFG8nOrTpTXIRlIeg6BBHx0qXnNOXBn3VmWcoBprz+Gfn5VWOKC3+uVBYRkdOM6NmKgxnZTF20jXrxcTwwvBMhdJeAclEhiIiU068GtudQRjazV+6kUe1q3D4ovKfTUSGIiJSTmXH/FR05fCKLqYu20aBmHCN7tfI6VrmpEEREKiAqypj6iy4cPZXD797YTN34OIad29TrWOUSlieVRURCSWx0FM9e350uLesy8ZVNrEk75HWkclEhiIgEQHxcDLNH9aRV/XjGv7iBL/eF37B4FYKISIDUqxnHvLG9qF09hlGz17P70AmvI5WJCkFEJICa1a3BvHG9yM3P58ZZ69h/PNPrSH5TIYiIBFhi49rMGd2TA8ezGDV7Pccyc7yO5JewLARNXSEioa5bq3pMv/F8Ur8/zvgXN5CZk+d1pFKFZSFo6goRCQcDkhvx6LVdWLvzMBMXbiI3L7TnPQrLQhARCRdXdW3OA8M7sXjL99z/5hc4F7rzHunCNBGRIBvTty2HMrJ5esl2Gtaqxr2XdvA6UrFUCCIileCeS5I5dCKLp5dsp37NOMb2a+t1pP+gQhARqQRmxp+uOpfDJ7L54ztbaFArjqu6Nvc61o/oHIKISCWJiY7iiZHd6N22Pve8+hkfb9vvdaQfUSGIiFSi6rHRvDCqB0lNajPhpY1s2nPE60g/UCGIiFSyOtVjeXFsTxrVrsaYuevZvv+415EAFYKIiCca167O/HG9iImK4qZZ69h39JTXkVQIIiJead2gJi+O7cnxzFxumr2OIyeyPc2jQhAR8dA5zRJ4YVQP9hw+yZi56zmZnetZFhWCiIjH+rRrwJMju/H53qNMeGkjOR5NcaFCEBEJAcPObcpff9qZpSkHmPLaZ+TnV/4UF2F5YZqZDQeGJyYmeh1FRCRgRvZqxaET2UxdtI16NeP4w5WdMLNK235Y7iFotlMRqap+NbA9Y/q2Yc7KXTz78Y5K3XZY7iGIiFRVZsbvr+jEYd+eQoOacYzs1apStq1CEBEJMVFRxtRrunDkZA6/e2Mz9WrGcek5TYO/3aBvQUREyiwuJorpN3TnvBZ1uXPhJtakHQr6NlUIIiIhKj4uhjmje9KqfjzjX9zAln3Hgro9FYKISAirVzOOeWN7Uat6DDfNXseeQyeDti0VgohIiGtWtwbzxvYiNz+fG2ev5cDxrKBsR4UgIhIGkprUZvbonuw/lsWo2es4lpkT8G2oEEREwkT3VvV47obuHD6RzXfpmQFfv4adioiEkYEdGvPxlIFUj40O+Lq1hyAiEmaCUQagQhARER8VgoiIACoEERHxUSGIiAigQhARER8VgoiIACoEERHxMecq/76dgWJmB4DdFVxNApBeye8v63v8fX1pryvt+YbAwTLkCmUV/fcaKtsMxDqD/RkN1OeztNfo8xm47bZ2zjX6j2edcxH9Bcyo7PeX9T3+vr601/nx/Aav/32Eyr/XUNlmINYZ7M9ooD6fpb1Gn8/gb1eHjOAfHry/rO/x9/Wlva6iv2s48eJ3DcY2A7HOYH9GA/X5LOt2w5lXv2eJ2w3rQ0YSWGa2wTnXw+scIsXR5zP4tIcgRc3wOoBICfT5DDLtIYiICKA9BBER8VEhiIgIoEIQEREfFYL4xcyuNrMXzOwtM7vE6zwiRZlZOzObZWave50lnKkQIoCZzTaz/Wb2xWnLh5nZNjPbbma/KWkdzrk3nXPjgdHAiCDGlQgToM9nmnNuXHCTVn0aZRQBzKw/kAHMc86d61sWDaQAFwN7gfXAdUA08OBpqxjrnNvve9+jwALn3MZKii9VXIA/n687566prOxVTYzXAST4nHPLzKzNaYt7Adudc2kAZvYKcJVz7kHgytPXYWYGPAS8qzKQQArE51MCQ4eMIldz4Osij/f6lp3JncBQ4Bozuy2YwUQo4+fTzBqY2XSgm5n9NtjhqirtIUQuK2bZGY8fOueeBJ4MXhyRHynr5/MQoD9UKkh7CJFrL9CyyOMWwD6PsoicTp9PD6gQItd6IMnM2ppZHDASeNvjTCKF9Pn0gAohApjZQmA10MHM9prZOOdcLnAHsAj4CnjVOfellzklMunzGTo07FRERADtIYiIiI8KQUREABWCiIj4qBBERARQIYiIiI8KQUREABWCCGY22czig7j+ZoGYp9/MMgKRR+RMdB2CRDwz2wX0cM4d9DpLScwswzlXy+scUnVpD0EihpnVNLN/mtlnZvaFmY0ws4lAM2CJmS3xve4SM1ttZhvN7DUzq+VbvsvM/tfM1vm+EovZxgAz+9T3tcnMaptZm8Kbv5jZzCLPHzCzB3zLp5jZejP73Mz+p4Tf4S++/GvMrEkw/jlJ5FIhSCQZBuxzznXx3YjlPd8srvuAQc65QWbWELgfGOqc6w5sAO4uso5jzrlewNPAtGK2cS9wu3OuK3ARcKrok865m33PXQUcAub6bkmaRME9ALoC5/tuGnO6msAa51wXYBkwvuz/CETOTIUgkWQzMNT3V/5Fzrn0Yl7TB+gErDSzT4FRQOsizy8s8v2CYt6/EnjMt+dR1zcnz4+YWXXgNeAO59xu4BLf1yZgI3A2BQVxumzgHd/PnwBtzvyripSd7ocgEcM5l2Jm5wOXAw+a2WLn3B9Pe5kB7zvnrjvTas7wc+E2HjKzf/q2scbMhgKZp71sOvB359wHRbb5oHPu+VJ+hRz375N+eei/Xwkw7SFIxDCzZsBJ59xLwCNAd99Tx4Havp/XAH0Lzw+YWbyZJRdZzYgi31cXs432zrnNzrn/peBw09mnPX87UNs591CRxYuAsUXOVTQ3s8YV+FVFykV/YUgk6QxMNbN8IAeY4Fs+A3jXzL71nUcYDSw0s2q+5++n4IbvANXMbC0Ff0wVtxcx2cwGUfAX/BbgXeCsIs/fC+T4DkcBTHfOTTezjsDqgltXkwHcAOyv6C8sUhYadirip3AZnipSXjpkJCIigPYQRETER3sIIiICqBBERMRHhSAiIoAKQUREfFQIIiICqBBERMTn/wH5IaufEml2JgAAAABJRU5ErkJggg==\n", + "text/plain": [ + "<Figure size 432x288 with 1 Axes>" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" } ], "source": [ - "error\n" + "plt.loglog(grid,error_eukl)\n", + "plt.ylabel('error in Euklidean norm')\n", + "plt.xlabel('step size h')" ] }, { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] - }, - { - "cell_type": "code", - "execution_count": null, + "cell_type": "markdown", "metadata": {}, - "outputs": [], - "source": [] + "source": [ + "Obviously the algorithm is not consistent, although it is supposed to be of order 2.\n", + "Even after several hours of debugging, I couldn't figure out the problem." + ] } ], "metadata": {