diff --git a/rb-intro/solutions/Sheet2_Solutions.ipynb b/rb-intro/solutions/Sheet2_Solutions.ipynb
new file mode 100644
index 0000000000000000000000000000000000000000..0eb222f526544b9cb2e504692389172e73638f15
--- /dev/null
+++ b/rb-intro/solutions/Sheet2_Solutions.ipynb
@@ -0,0 +1,907 @@
+{
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "# pyMOR RB-Intro\n",
+ "\n",
+ "# Exercise Sheet 2 - Solutions"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Problem 1"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Problem 1 a)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Choose an arbitrary model `m` from Sheet 1, Problem 2. Add a `CubicParameterSpace` as the `parameter_space` of the `Model` using:\n",
+ "\n",
+ "```python\n",
+ "m = m.with_(parameter_space=CubicParameterSpace(...))\n",
+ "``` \t\n",
+ "\n",
+ "While the `ParameterType` only defines the components of a `Parameter` and their shape, the parameter space defines a subset $\\mathcal{P} \\subset \\mathbb{R}^P$ of valid values for each component."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Solution"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "First we define again the models from Sheet 1, Problem 2:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 1,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stderr",
+ "output_type": "stream",
+ "text": [
+ "00:01 |WARNING|BitmapFunction: Image R.png not in grayscale mode. Converting to grayscale.\n",
+ "00:01 |WARNING|BitmapFunction: Image B.png not in grayscale mode. Converting to grayscale.\n"
+ ]
+ }
+ ],
+ "source": [
+ "from pymor.basic import *\n",
+ "set_log_levels({'pymor': 'WARN'})\n",
+ "\n",
+ "# Problem 2 a)\n",
+ "domain = RectDomain(bottom='neumann')\n",
+ "neumann_data = ExpressionFunction('-cos(pi*x[...,0])**2*neum', 2, (), parameter_type= {'neum': ()})\n",
+ "diffusion = ExpressionFunction(\n",
+ " '1. - (sqrt( (np.mod(x[...,0],1./K)-0.5/K)**2 + (np.mod(x[...,1],1./K)-0.5/K)**2) <= 0.3/K) * 0.999',\n",
+ " 2, (),\n",
+ " values={'K': 10}\n",
+ ")\n",
+ "problem = StationaryProblem(\n",
+ " domain=domain,\n",
+ " diffusion=diffusion,\n",
+ " neumann_data=neumann_data\n",
+ ")\n",
+ "m_a, _ = discretize_stationary_cg(problem, diameter=1/100)\n",
+ "\n",
+ "# Aufgabe 2 b)\n",
+ "diffusion = ExpressionFunction(\n",
+ " '1. - (sqrt( (np.mod(x[...,0],1./K)-0.5/K)**2 + (np.mod(x[...,1],1./K)-0.5/K)**2) <= 0.3/K) * (1 - diffu)',\n",
+ " 2, (),\n",
+ " values={'K': 10},\n",
+ " parameter_type= {'diffu': ()}\n",
+ ")\n",
+ "problem = StationaryProblem(\n",
+ " domain=domain,\n",
+ " diffusion=diffusion,\n",
+ " neumann_data=neumann_data\n",
+ ")\n",
+ "m_b, _ = discretize_stationary_cg(problem, diameter=1/100)\n",
+ "\n",
+ "# Aufgabe 2 c)\n",
+ "diffusion_R = BitmapFunction('R.png', range=[1, 0])\n",
+ "diffusion_B = BitmapFunction('B.png', range=[1, 0])\n",
+ "diffusion = LincombFunction(\n",
+ " [ConstantFunction(1., 2), diffusion_R, diffusion_B],\n",
+ " [1., ExpressionParameterFunctional('-(1 - R)', {'R': ()}), ExpressionParameterFunctional('-(1 - B)', {'B': ()})]\n",
+ ")\n",
+ "problem = StationaryProblem(\n",
+ " domain=domain,\n",
+ " diffusion=diffusion,\n",
+ " neumann_data=ConstantFunction(-1, 2)\n",
+ ")\n",
+ "m_c, _ = discretize_stationary_cg(problem, diameter=1/100)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Currently, the only parameter space implemented in pyMOR is `CubicParameterSpace`, which represents spaces $\\mathcal{P}$ which are the Cartesian product of closed intervals. When instantiating `CubicParameterSpace`, the corresponding `ParameterType` along with (the same) lower and upper bounds for parameter components have to be provided. Alternatively, dictionary of intervals for the individual parameter components can be provided via the `ranges` argument:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 2,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "m_a = m_a.with_(parameter_space=CubicParameterSpace(m_a.parameter_type, -10, 10))\n",
+ "m_b = m_b.with_(parameter_space=CubicParameterSpace(m_b.parameter_type, \n",
+ " ranges={'diffu': [0.0001, 1.], 'neum': [-10, 10]}))\n",
+ "m_c = m_c.with_(parameter_space=CubicParameterSpace(m_c.parameter_type, 0.0001, 1.))"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Here we ant to set the `parameter_space` attribute of the models `m_a`, `m_b`, `m_c` to an appropriate space. However, all `Models` (and `Operators`) in pyMOR are immutable, which means that you are not allowed to assign to attributes of these objects. Thus, we have to create *new* `Models` in which the `parameter_space` attribute has been given the desired (new) value. To easily achieve this, each immutable object in pyMOR has a `with_` method that can be used in the form\n",
+ "\n",
+ "```\n",
+ "obj.with_(attrib1=value1, attrib2=value2)\n",
+ "```\n",
+ "\n",
+ "which will return a copy of `obj` where the attributes `attrib1`, `attrib2` have been assigned the values `value1` and `value2`. Thus, in the statement\n",
+ "\n",
+ "```\n",
+ "m_a = m_a.with_(...)\n",
+ "```\n",
+ "\n",
+ "a new, modified copy of `m_a` is created and assigned the same name `m_a`, discarding the original `Model` object.\n",
+ "\n",
+ "The immutability of `Models` and `Operators` in pyMOR ensures that those objects cannot unexpectedly change when passed to a function or assigned to another name. (Note that in pyMOR all function calls pass arguments by reference and assignment never creates copies.) The underscore in the method name `with_` is required as `with` is a keyword in Python."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "In the following we will fix one of the three `Models` we have defined:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 3,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# m = m_a\n",
+ "# m = m_b\n",
+ "m = m_c"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Problem 1 b)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "The main use of parameter spaces is that they allow to sample `Parameters` from the set $\\mathcal{P}$ they define, e.g. for training or validating a reduced order model.\n",
+ "\n",
+ "Solve the model for 30 different `Parameters` generated by the `sample_uniformly` and `sample_randomly` methods of `CubicParameterSpace`.\n",
+ "Store all solutions in a single `VectorArray` `U`. Then, visualize `U`. Moreover, compute the $H^1$ norms and semi norms of the vectors in `U`.\n",
+ "\n",
+ "**Hints:** Use the `solution_space` of the model to create an empty `VectorArray` of the right type using its `empty` method. Then use its `append` methods to add the solution vectors. The norms can be computed using the `h1_0_norm` and `h1_0_semi_norm` methods of the model."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Solution"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "We use `sample_randomly`:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 4,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "U = m.solution_space.empty()\n",
+ "for mu in m.parameter_space.sample_randomly(30):\n",
+ " U.append(m.solve(mu))"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Each `VectorArray` returned by `m.solve` contains a single vector (since `m` is stationary). After appending all these vectors to `U` we and up with an array of 30 vectors:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 5,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "30"
+ ]
+ },
+ "execution_count": 5,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "len(U)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Note that `appending` one `VectorArray` `V` to another array `U` will append (shallow) copies of the vectors in `V` to `U`. So modifying `U` will not affect `V`.\n",
+ "\n",
+ "A `VectorArray` containing multiple vectors can be visualized as a time series:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 6,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "application/vnd.jupyter.widget-view+json": {
+ "model_id": "5f379428ad59414eb4a1fa5225ad9a02",
+ "version_major": 2,
+ "version_minor": 0
+ },
+ "text/plain": [
+ "ThreeJSPlot(children=(HBox(children=(Renderer(children=(Renderer(camera=PerspectiveCamera(position=(0.0, 0.0, …"
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "m.visualize(U)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Finally, we compute the norms:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 7,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "H^1-Norm: [0.55094866 0.54427404 0.62515333 0.60929032 0.54537415 0.63415985\n",
+ " 0.56591976 0.61230459 0.56571096 0.57038588 0.58929855 0.57455573\n",
+ " 0.54938834 0.57832506 0.65540032 0.5807036 0.60386921 0.53527094\n",
+ " 0.62674714 0.55087898 0.59563011 0.62360886 0.56590276 0.56524171\n",
+ " 0.57989481 0.53575488 0.53441772 0.54179062 0.63666902 0.64178785]\n",
+ "H^1_0-Norm: [0.53802876 0.53149403 0.61061681 0.59586326 0.53255317 0.62094288\n",
+ " 0.55265167 0.59789037 0.5522074 0.55671071 0.57560529 0.56077105\n",
+ " 0.5364431 0.56461142 0.64176098 0.56705885 0.59051618 0.52294491\n",
+ " 0.61245129 0.5378326 0.58175874 0.61027344 0.55249642 0.55174501\n",
+ " 0.56618618 0.52340308 0.52213937 0.52916328 0.62212136 0.62757799]\n"
+ ]
+ }
+ ],
+ "source": [
+ "print('H^1-Norm:', m.h1_0_norm(U))\n",
+ "print('H^1_0-Norm:', m.h1_0_semi_norm(U))"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "**Warning: `m` also has norms `m.h1_norm` und `m.h1_semi_norm` corresponding to the inner product `Operators` `m.h1_product` and `m.h1_semi_product`, which only differ in the boundary treatment for the underlying system matrices. While for computing norms of $H^1_0$ functions both versions will yield the same results, the computation of Riesz representatives requires choosing the right product for the space ($H^1$ or $H^1_0$) for which the Riesz representatives are to be computed.**"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Aufgabe 2"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Problem 2 a)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Let `U` be an arbitrary `VectorArray` and `V` a length 1 `VectorArray` of the same dimension.\n",
+ "Derive a linear equation system for the coefficients of the orthogonal projection (in Euclidean norm) of the vector contained in `V` onto the linear space spanned by the basis vectors contained in `U`.\n",
+ " \n",
+ " "
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Solution"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Let $v_{proj}$ be the orthogonal projection of $v$ onto the linear space spanned by the vectors in $U$. Thus, $v - v_{proj}$ is orthogonal to this space and in particular to each vector in $U$. With $u_i$, $i = 1, \\ldots, N$ being the vectors in $U$, we have:\n",
+ "\n",
+ "$$(v - v_{proj}, u_i) = 0, \\qquad i = 1, \\ldots, N.$$\n",
+ "\n",
+ "Let $\\lambda_j$, $j = 1, \\ldots, N$ be the coefficients of $v_{proj}$ w.r.t. the $u_j$ basis, i.e. $\\sum_{j=1}^N \\lambda_j u_j = v_{proj}$. Then:\n",
+ "\n",
+ "$$ \\sum_{j=1}^N \\lambda_j (v_i, v_j) = (v, u_i), \\qquad i = 1, \\ldots, N.$$\n",
+ "\n",
+ "With $M_{i,j} := (v_i, v_j)$, $R := [(v, u_1), \\ldots, (v, u_N)]^T$ and $\\Lambda := [\\lambda_1, \\ldots, \\lambda_N]^T$, we obtain the linear equation system\n",
+ "\n",
+ "$$ M \\cdot \\Lambda = R,$$\n",
+ "\n",
+ "which determine $\\lambda_i$ and, thus, $v_{proj}$."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "collapsed": true
+ },
+ "source": [
+ "## Problem 2 b)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Assemble and solve the equation system for the `VectorArray` `U` from Problem 1 (b) and a `VectorArray` `V` containing an additional solution of `m`. Form the linear combination of the obtained coefficients with the basis vectors in `U`. Visualize the result together with `V`. Also visualize the projection error.\n",
+ "\n",
+ "**Hints:** Use the `gramian`, `inner` and `lincomb` methods of `VectorArrayInterface`. `m.solve` can be passed a tuple of `VectorArrays` to simultaneously plot multiply functions."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Solution"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Every `VectorArray` has a `gramian` method, which computes precisely the Gramian Matrix $M$ defined above w.r.t the Euclidean inner product. `gramian` is a special case of the more general `inner` method, which in case of `U.inner(V)` returns the matrix $(u_i, v_j)$ of all Euclidean inner products between the vectors in `U` und `V`. \n",
+ "The linear combination $\\sum_{j=1}^N \\lambda_j u_j$ can efficiently formed using the `lincomb` method of `U`, which is passed a 1d oder 2d NumPy array of coefficients $\\lambda_j$.\n",
+ "\n",
+ "Thus, we can compute $v_{proj}$ as follows:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 8,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "V = m.solve(m.parameter_space.sample_randomly(1)[0])\n",
+ "M = U.gramian()\n",
+ "rhs = U.inner(V)\n",
+ "v = np.linalg.solve(M, rhs)\n",
+ "V_proj = U.lincomb(v.T) # a row vector of coefficients is expected"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "When we pass a tuple of `VectorArrays` to `visualize`, we will get multiple plots for each individual array. The `separate_colorbars` parameter defines if all plots should share the same color bar or not:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 9,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "application/vnd.jupyter.widget-view+json": {
+ "model_id": "1ae4d6115d6b44f3a6cedb7d8bf307a8",
+ "version_major": 2,
+ "version_minor": 0
+ },
+ "text/plain": [
+ "ThreeJSPlot(children=(HBox(children=(Renderer(children=(Renderer(camera=PerspectiveCamera(position=(0.0, 0.0, …"
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "m.visualize((V, V_proj, V - V_proj), separate_colorbars=True)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Problem 2 c)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Repeat parts (a) and (b) for the $H^1_0$-orthogonal projection onto `U`. Compute Euclidean, the $H^1_0$ and $L^2$ norms of the projection error for both projections.\n",
+ "\n",
+ "**Hint:** Use the `Operators` in `m.products` as the `product` argument to `gramian` and `inner`."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Solution"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 10,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "application/vnd.jupyter.widget-view+json": {
+ "model_id": "072d29e1b9f844e997de2661e95e6871",
+ "version_major": 2,
+ "version_minor": 0
+ },
+ "text/plain": [
+ "ThreeJSPlot(children=(HBox(children=(Renderer(children=(Renderer(camera=PerspectiveCamera(position=(0.0, 0.0, …"
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "M = U.gramian(product=m.products['h1_0_semi'])\n",
+ "rhs = U.inner(V, product=m.products['h1_0_semi'])\n",
+ "v = np.linalg.solve(M, rhs)\n",
+ "V_h1_proj = U.lincomb(v.T)\n",
+ "\n",
+ "m.visualize((V, V_h1_proj, V - V_h1_proj), separate_colorbars=True)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 11,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "eukld.-Norm: eukld.: 4.16e-07 H^1_0: 2.43e-08 \n",
+ "H^1_0-Norm: eukld.: 1.25e-05 H^1_0: 1.72e-07 \n",
+ "L^2-Norm: eukld.: 3.97e-07 H^1_0: 2.44e-08 \n"
+ ]
+ }
+ ],
+ "source": [
+ "for norm_name, norm in [('eukld.', None), ('H^1_0', m.h1_0_semi_norm), ('L^2', m.l2_0_norm)]:\n",
+ " print('{:<16}'.format(norm_name + '-Norm:'), end='')\n",
+ " for proj_name, VV in [('eukld.', V_proj), ('H^1_0', V_h1_proj)]:\n",
+ " if norm is None:\n",
+ " err = ((V - VV).l2_norm() / V.l2_norm())[0]\n",
+ " else:\n",
+ " err = (norm(V - VV) / norm(V))[0]\n",
+ " print('{}: {:.2e} '.format(proj_name, err), end='')\n",
+ " print()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Problem 2 d)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Compute 100 random solutions of `m`.\n",
+ "Compute the maximum projection error for these solution vectors onto the first `N` vectors in `U`.\n",
+ "Plot the projection error in dependence of `N`."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Solution"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "First, we compute the 100 test snapshots:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 12,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "...................................................................................................."
+ ]
+ }
+ ],
+ "source": [
+ "V = m.solution_space.empty()\n",
+ "for mu in m.parameter_space.sample_randomly(100, seed=99):\n",
+ " V.append(m.solve(mu))\n",
+ " print('.', end='', flush=True)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "To optimize the computation of the projection matrix and the right-hand side for varying basis sizes, we first compute those for the full basis and then extract appropriate sub-matrices. Similarly to NumPy arrays, we can slice or index a `VectorArrays` to obtain a view onto a subset of vectors in the original array. (No copy is created when slicing or indexing.)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 13,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "M = U.gramian(product=m.h1_0_semi_product)\n",
+ "rhs = U.inner(V, product=m.h1_0_semi_product)\n",
+ "errors = []\n",
+ "for N in range(len(U)):\n",
+ " try:\n",
+ " v = np.linalg.solve(M[:N, :N], rhs[:N, :])\n",
+ " except np.linalg.LinAlgError:\n",
+ " v = np.zeros((N, len(V)))\n",
+ " V_proj = U[:N].lincomb(v.T)\n",
+ " errors.append(np.max(m.h1_0_semi_norm(V - V_proj)))"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Depending on the problem, it can happen that the vectors in `U` can become (numerically) linear dependent. In this case, `np.linalg.solve` will raise a `np.linalg.LinalgError`. We catch this error in the `try`-`except` block and set `V_proj` to 0 in this case."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 14,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "<Figure size 640x480 with 1 Axes>"
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "from matplotlib import pyplot as plt\n",
+ "plt.figure()\n",
+ "plt.semilogy(errors)\n",
+ "plt.show()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Problem 2 e)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Compute the condition number of the projection matrices in dependence of `N` and plot the result."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Solution"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "We compute the condition numbers of the system matrices with `numpy.linalg.cond`, where we need to take care not to pass a $0 \\times 0$-Matrix:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 15,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAD9CAYAAAC7iRw+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xt0lfWd7/H3N3cI4Z4ECIQAcle5NKLWG7W1QqtFrceK7WmnZcrYqTOd03XOVKczx+msmU5npjN11tQZh1a0TmdEj0M7qAiVLhUcbQ2iYACREC5JyJUQyIVcdvb3/JGtTZFASHby7MvntVZW9vPbl+f75Ek+PPz2b/9+5u6IiEjySAm6ABERGV4KfhGRJKPgFxFJMgp+EZEko+AXEUkyCn4RkSSj4BcRSTIKfhGRJJMW7Rc0s5nAt4Ex7n5npO064POR/S1w949Ge78iItI//briN7P1ZlZnZqVnta8wswNmVmZm9wO4e7m7r+n9OHff4e73As8BP4lW8SIicvH6e8X/OPBD4In3G8wsFXgYuAmoBErMbJO77zvP69wDrDnP/QBMnDjRi4qK+lmaiIgAvPnmmw3unnuhx/Ur+N19u5kVndW8DChz93IAM9sArALOGfxmVgiccvfmC+2vqKiInTt39qc0ERGJMLOj/XncYN7cLQAqem1XAgVmNsHMHgGWmNkDve5fAzzW14uZ2Voz22lmO+vr6wdRloiInE/U39x19xPAvedof/ACz1sHrAMoLi7WlKEiIkNkMFf8VcC0XttTI20iIhLDBhP8JcBsM5thZhnA3cCm6JQlIiJDpb/DOZ8EXgfmmlmlma1x9xBwH7AV2A887e57h65UERGJhv6O6lndR/tmYHNUKxIRkSGlKRtERJJM1Ef1iIjIh7k7HaEwp9u7aGkP0dIRoqU9xOkPbnfR3B4if0wWdxVPu/ALDoKCX0TkIoS6wzS2dnKitZOmti5Oneni1JlOTp3p6rX9m6+mtq4Pwj4UvvBI9WsvmajgFxEZDm2dIY41tlHf3EF9cwcNLe9/7/yt7ca2TryP/E5NMcaMSP/ga9zIDGZMzGZ0Vjo5WWmMykojJzONnKx0RmVGtrPSyMnsuT87M42MtKHvgVfwi0hSOdHSQVldC4fqWymra6GsvoVDdS1UNZ350GMz01LIzckkNyeTaeNHsqRwXM/2qAwmjMpk7Ih0Ro9IZ+zInqAflZmGmQVwVBdHwS8icamrO0xbZzdnOrtp6wxxpuv9292/dbu5vYvDDa2RsG/hZFvXB68xIj2VWXnZFBeN4+7caczIzSYvJ4vcnEwmjsqImyC/WAp+EYkrja2d3PzQduqbO/r9nAnZGczKHcWKSydzSd4oZuVmc0neKKaMGUFKSuIF+4Uo+EUkrry4r4b65g5+7/qZ5OZkMjIjjZEZqWSlpzIyo+drREYqI9JTGZmRRnZmKjlZ6UGXHVMU/CISV7aU1jB13AjuXzkvIbthhoM+wCUicaO5vYv/LjvBioWTFPqDoOAXkbjx0oF6OrvDrLh0UtClxDUFv4jEja2lNUwclcnSwnFBlxLXFPwiEhfau7p56UAdn1yYn5QjcaJJwS8iceHVgw20dXazYqG6eQYr6sFvZjPN7FEze6ZX23Iz22Fmj5jZ8mjvU0QS35a9NeRkpXHVzAlBlxL3+rsQy3ozqzOz0rPaV5jZATMrM7P7Ady93N3XnPUSDrQAWfQsyi4i0m+h7jDb9tfyifn5wzKXTaLr70/wcWBF7wYzSwUeBlYCC4DVZragj+fvcPeVwLeA7wysVBFJVm8cbqSprYub1c0TFf0KfnffDjSe1bwMKItc4XcCG4BVfTw/HLl5EsgcYK0ikqS27K0hKz2FG+bkBl1KQhjM/5kKgIpe25VAgZlNMLNHgCVm9gCAmd1hZv8K/Bvww3O9mJmtNbOdZrazvr5+EGWJSCIJh52te2u4YU4uIzJSgy4nIUR9ygZ3PwHce1bbRmDjBZ63DlgHUFxcfOHVCkQkKeyubKL2dIc+tBVFg7nirwJ6LxMzNdImIhI1W/bWkJZi3Dg3P+hSEsZggr8EmG1mM8wsA7gb2BSdskREetap3Vpaw9WzJjBmpGbYjJb+Dud8EngdmGtmlWa2xt1DwH3AVmA/8LS77x26UkUk2bxX28KRE23q5omyfvXxu/vqPto3A5ujWpGISMSW0hrM4KYF6uaJJn0SQkRi1pa9NXykcBx5OVlBl5JQFPwiEpOOnWhjf/VpdfMMAQW/iMSkrXtrAPRp3SGg4BeRmLRlbw0LJo9m2viRQZeScBT8IhJz6k63s+vYSXXzDBEFv4jEnF/sq8UdBf8QUfCLSMzZureGGROzmZ03KuhSEpKCX0Riyqm2Ll4/dIKbF07CTEssDgUFv4jElF++W0so7OrmGUIKfhGJKVtKa5g0OovLC8YEXUrCUvCLSMxo6wyx/WA9Ny/MJyVF3TxDRcEvIjFj+3v1tHeFuVndPENKwS8iMWNLaQ3jRqazrGh80KUkNAW/iMSEzlCYX75bxyfm55OWqmgaSlH/6ZrZTDN71Mye6dU238weMbNnzOxr0d6niMS/18tP0Nwe0mieYdDfhVjWm1mdmZWe1b7CzA6YWZmZ3Q/g7uXuvqb349x9v7vfC9wFXBOt4kUkcWwprSE7I5VrLpkYdCkJr79X/I8DK3o3mFkq8DCwElgArDazBX29gJl9BngeLdwiIr20d3Xz018d5bk9x1k+L4+s9NSgS0p4/V2Ba7uZFZ3VvAwoc/dyADPbAKwC9vXxGpuATWb2PPAfAy1YRBLDydZO/u1XR/nJa0c40drJomlj+aOPzw66rKTQr+DvQwFQ0Wu7ErjSzCYAfwUsMbMH3P2vzWw5cAeQSR9X/Ga2FlgLUFhYOIiyRCSWVTS28eirh3mqpIIzXd3cOC+P37t+JstmjNcUDcNkMMF/Tu5+Arj3rLaXgZcv8Lx1wDqA4uJij3ZdIhKs0qpT/Ov2cja/U02KwarFBay9fiZz8nOCLi3pDCb4q4BpvbanRtpERABwd3YcbGDd9nJeLWtgVGYav3vtDL58zQwmjdE6ukEZTPCXALPNbAY9gX83cE9UqhKRYePuHGts4+iJNjpDYbq6w3R2h+l4/3Yo/Jv2UJjObqcj1E17V5iOrm7aQ910dIVpj7S1d3XTEer53toR4mRbF3k5mdy/ch73XFnI6Kz0oA856fUr+M3sSWA5MNHMKoEH3f1RM7sP2AqkAuvdfe+QVSoiUdHaEWJ3ZRNvHWvirWMneetYEydaO/v9/Iy0FDJTU8hMTyUrPYXMtBSy0lMjXymMGZFOVnoKWWmpZKansKRwHKsWTyEzTaN1YkV/R/Ws7qN9MxqeKRKz3J3DDa28dayJXcdOsutYEwdqThOOvIs2Mzeb5XPzWDp9LLPzcshKTyE9NYWMtBQyzvqenppCeqrpDdgEEPU3d0UkGGc6u3mvtpkDNc28W9PMgdrT7Dt+mpNtXQCMykxj8bSx3PexS1gyfRyLp45lXHZGwFVLEBT8InGmO+wcOdH6m4CvOc2BmmaONrbhkSv5rPQU5uTn8MkFk1hSOJYlheO4JG8UqZrqWFDwi8SsUHeYY41tHKxroayuhYO1zR/c7giFAUgxKJqQzfzJo7ltSQHzJuUwd9JoCsePVMhLnxT8IjHgcEMr+6tPc7C2hYN1zZTVtVBe30pnd/iDxxSMHcEleaO4auYE5k7KYd6kHGbn5TAiQ2+aysVR8IsE6GBtM3+z5V227a/7oG3a+BHMzsvhhjm5zM7PYXbeKGbljWJUpv5cJTr0myQSgLrT7fxg23s8VVJBdkYa//uTc1g+N4+ZudmMzNCfpQwt/YaJDKOWjhDrtpfzo+3ldHWH+eLVRfzBjZcwYVRm0KVJElHwiwyDru4wT5VU8NC2gzS0dPDpyyfzfz45l6KJ2UGXJklIwS8yhNydX+yr5W+2vEt5fSvLisbzoy9+hCWF44IuTZKYgl9kCLg7O4+e5G+3vEvJkZPMys3mR18s5hPz8/TJVwmcgl8kSjpDYX5VfoJt+2vZtq+W46fayc3J5Lu3X8ZdxVO1gLjEDAW/yCA0tXXy0oE6tu2r45X36mnpCJGVnsJ1s3P5xidmc8vlU8jWMEyJMfqNFLlIRxpa2ba/lhf31bLz6Em6w05uTia3LprMJ+bnc80lE7VurMQ0Bb/IeYTDzsG6Ft48epJdx07y5tGTHG5oBWDepBx+f/ksPjE/n8sKxpCiKRIkTij4RXo5daaLtyua2BUJ+rePNdHcEQJgfHYGSwvH8qWrp/Px+flMGz8y4GpFBibqwW9mM4FvA2Pc/c6+2kRixWtlDWzafZxdx05ysK4F957Jz+bk5/CZxVNYWjiOj0wfx/QJIzUiRxJCf1fgWg/cAtS5+6W92lcA/0jPClw/dvfvuXs5sMbMnnn/cedqE4kF//V2Ff/rqbcZlZnGksJx3HJ5T9AvmjaGHC0RKAmqv1f8jwM/BJ54v8HMUoGHgZuASqDEzDa5+75oFykyFN4P/SuKxvPYl6/QHDmSNPo1sNjdtwONZzUvA8rcvdzdO4ENwKoo1ycyJBT6kswG84mSAqCi13YlUGBmE8zsEWCJmT0AcK62s5nZWjPbaWY76+vrB1GWyPkp9CXZRf033t1PAPdeqO0cz1sHrAMoLi72aNclAgp9ERjcFX8VMK3X9tRIm0hMUuiL9BhM8JcAs81shpllAHcDm6JTlkh0KfRFfqNfwW9mTwKvA3PNrNLM1rh7CLgP2ArsB552971DV6rIwCj0RX5bv/4C3H11H+2bgc1RrUgkihT6Ih+meWIlYSn0Rc5NwS8J6fk91Qp9kT4o+CXh1Jxq5/7/3MPiaWMV+iLnoOCXhOLu/OnPS+kKh/nB5xYr9EXOQcEvCeWF0hq27a/lmzfNYfqE7KDLEYlJCn5JGE1tnfzf/9rLpQWj+co1M4IuRyRm6f/BkjC+u3k/J9s6+clXrtDC5iLnob8OSQj/XdbA0zsrWXv9TBZOGRN0OSIxTcEvce9MZzcPbHyHogkj+cbHZwddjkjMU1ePxL2Htr3HscY2nvzqVWSlpwZdjkjM0xW/xLXSqlP8aEc5q5dN4+pZE4IuRyQuKPglbnV1h/njZ/YwYVQm96+cH3Q5InFDXT0St3684zD7qk/zyBeWMmaEFkYX6S9d8UtcOtzQykPb3mPFwkmsuHRy0OWIxJVhCX4zW2BmT5vZv5jZncOxT0lc7s4DG/eQkZbCd1YtDLockbgz4OA3s/VmVmdmpWe1rzCzA2ZWZmb3R5pXAv/k7l8DvjiIekV4qqSCX5U38iefmk/+6KygyxGJO4O54n8cWNG7wcxSgYfpCfoFwGozWwD8G3C3mf0doKEXMmB1p9v5q837uWrmeO6+YtqFnyAiHzLg4Hf37UDjWc3LgDJ3L3f3TmADsMrd69z968D9QMOAq5Wk9+CmvXSEwvz1HZdjZkGXIxKXoj2qpwCo6LVdCVxpZkXAnwDZwN+d64lmthZYC1BYWBjlsiTeHW5o5Zk3K3ihtIY/XjGXGRM186bIQA3LcE53P0Ik1M/zmHXAOoDi4mIfhrIkhrk7B+taeOGdGl4orebdmmYAbpyXx1evmxlwdSLxLdrBXwX07nidGmkTuSB3Z+/x02wprWFzaTXl9a2YQfH0cfzZLQtYcekkCsaOCLpMkbgX7eAvAWab2Qx6Av9u4J4o70MSRKg7TEtHiPKGVraW1vBCaQ3HGttIMbhq5gS+/NEibl44iTyN3BGJqgEHv5k9CSwHJppZJfCguz9qZvcBW4FUYL27741KpRI3mtu7eKqkgvqWDprbQzS3hzh9povm9q7fbLd30dbZ/cFz0lONj86ayNc/NoubFkxifHZGgEcgktgGHPzuvrqP9s3A5gFXJHGtpSPE7zxWwptHT5KRmkJOVhqjR6STk5VGTlYa+aOzIrfTGZ3V0543OpPrZudq2gWRYaK5eiRq2jpDfOWxEt6uaOKfP7+UT12mqRREYpHm6pGoONPZzVceL2Hn0UYe+txihb5IDFPwy6C1d3Xzu0+U8MbhRn7wucXcumhK0CWJyHmoq0cGpb2rm68+sZPXDp3g+3cuYtXigqBLEpEL0BW/DFhHqJt7f/omOw428DefvZzPfmRq0CWJSD8o+GVAOkNhfv+nu3j5QD1/fcdl3FWsCdNE4oWCXy5aV3eYr//HLn75bh1/edulrF6muZVE4omCXy5KV3eYP3zyLV7cV8t3PrOQL1w1PeiSROQiKfil30LdYf7oqbd5obSGP7tlAV/6aFHQJYnIACj4pV9Onenim0/v5vk91Xz7U/NZc+2MoEsSkQHScE45r6MnWnnsv4/w9M4K2jq7+daKeXz1ek2LLBLPFPzyIe7OG4cb+fGrh9m2v5a0FOPWRVNYc+0MFk4ZE3R5IjJICn75QGcozPPvHOfRVw9TWnWacSPTue9jl/A/r5quqZFFEoiCX2hq6+Tff32MJ14/Qu3pDmblZvPd2y/j9iUFjMhIDbo8EYmyYQl+M7sO+Hxkfwvc/aPDsV85t3DYKW9oZXdFE78+fIJnd1dzpqub62ZP5HufvZwbZueSkqKFzEUS1WAWYlkP3ALUufulvdpXAP9Iz0IsP3b377n7DmCHmd1GzypdMozqmtvZXXGK3RVNvF3RxO7KJprbQwBkZ6Ryy+WTWXPdDOZNGh1wpSIyHAZzxf848EPgifcbzCwVeBi4CagESsxsk7vvizzkHmDNIPYp/VBadYrXDjX0hHzFKaqazgCQmmLMm5TDrYumsHjaWBZPG8us3FGk6upeJKkMZgWu7WZWdFbzMqDM3csBzGwDsArYZ2aFwCl3bx7oPqVv7s6rZQ3880uHeL38BADTxo9g6fRxfPmaIhZPG8vCKWPUZy8iUe/jLwAqem1XAldGbq8BHuvriWa2FlgLUFiouV/6Kxx2tu6t4Z9fPsQ7VafIH53Jn356PrctKWDiqMygyxORGDRso3rc/cEL3L8OWAdQXFzsw1JUHOsMhfn521U88sohyutbKZowku/dcRm3Ly0gM01X9SLSt2gHfxXQe37eqZE2iZK2zhAb3qjgxzvKOX6qnQWTR/PDe5aw8tLJ6qsXkX6JdvCXALPNbAY9gX83PW/oyiCEusPUnG7nZ7uqeOy1IzS2drJsxni+e8dl3DAnFzMFvoj032CGcz4JLAcmmlkl8KC7P2pm9wFb6RnOud7d90al0gR1prObmtPtVJ86Q+3pdqpPtVN7KvI9st3Q0kE40vn18Xl5fG35LIqLxgdbuIjErcGM6lndR/tmYPOAK0oi33vhXR555dCH2nOy0pg0OotJY7KYk5/D5DFZ5I/Jonj6eOZOygmgUhFJJJqyISDtXd38+6+PcuWM8XzuimlMGt0T7pNGZ5GdqdMiIkNHCROQbftraW4P8Qc3zuba2RODLkdEkogWYgnIxl1VTBqdxdWzJgRdiogkGQV/AOqbO3jlvXpuW1KgIZgiMuwU/AHYtPs43WHnjqUFQZciIklIwR+AjbsquaxgDHPyNUJHRIafgn+YHahpZu/x07raF5HAKPiH2ca3Kj9Yw1ZEJAgK/mHUHXZ+/lYVy+fmauZMEQmMgn8YvXaogdrTHdyxdGrQpYhIElPwD6ONu6oYnZXGjfPygi5FRJKYgn+YtHSE2FJawy2LppCVrvnyRSQ4Cv5hsqW0hjNd3XxWo3lEJGAK/mGycVcl0yeMZGnhuKBLEZEkNyzBb2bLzWyHmT1iZsuHY5+xpKrpDK+Xn+COJVO1aIqIBG7AwW9m682szsxKz2pfYWYHzKzMzO6PNDvQAmTRswB7Uvn5W1W4w+1L1M0jIsEbzBX/48CK3g1mlgo8DKwEFgCrzWwBsMPdVwLfAr4ziH3GHXdn465KlhWNp3DCyKDLEREZePC7+3ag8azmZUCZu5e7eyewAVjl7uHI/SeBpPrk0p7KUxyqb9UUDSISM6K9EEsBUNFruxK40szuAG4GxgI/PNcTzWwtsBagsLAwymUFZ+OuSjLSUvjU5ZODLkVEBBimFbjcfSOw8QKPWQesAyguLvbhqGuodYbCPLunmpsW5DM6Kz3ockREgOiP6qkCpvXanhppS0qvvFdPY2unxu6LSEyJdvCXALPNbIaZZQB3A5uivI+4sXFXJRNHZXDd7NygSxER+cBghnM+CbwOzDWzSjNb4+4h4D5gK7AfeNrd90an1PjS1NbJL/fX8ZlFBaSn6nNyIhI7BtzH7+6r+2jfDGwecEUJ4rk91XR2hzWaR0Riji5Fh8jGXZXMzc9h4ZTRQZciIvJbFPxD4HBDK7uONXHH0gJN0SAiMUfBPwR+tquSFIPbNEWDiMQgBX+UhcPOxrequOaSieSPzgq6HBGRD1HwR1nJkUYqT57hs1peUURilII/ikLdYb7/iwPkZKXxyYX5QZcjInJOwzJlQ7J4+KVDlBw5yUOfW8zIDP1oRSQ26Yo/SnYeaeQff/kety8p0Ju6IhLTFPxRcOpMF9/Y8DZTx43kL1YtDLocEZHzUn/EILk73/7ZO9ScbueZe68mR7NwikiM0xX/ID3zZiXP7anmmzfNYYkWUheROKDgH4TDDa08uGkvV80cz703zAq6HBGRflHwD1BnKMw3NrxFemoKP/jcYlJTNDWDiMQH9fEP0N+/eIA9lad45AsfYfKYEUGXIyLSb7riH4BXDzbwr6+Uc8+Vhay4dFLQ5YiIXJRhCX4zm29mj5jZM2b2teHY51BpbO3km0+/zazcbP7s0wuCLkdE5KINZgWu9WZWZ2alZ7WvMLMDZlZmZvcDuPt+d78XuAu4ZnAlB8fd+eNndtPU1sU/rV7KiIzUoEsSEblog7nifxxY0bvBzFKBh4GVwAJgtZktiNz3GeB54nh1rp/+6ijb9tdx/8p5LNACKyISpwYc/O6+HWg8q3kZUObu5e7eCWwAVkUev8ndVwKfH+g+g3Sgppm/fH4/y+fm8uVrioIuR0RkwKI9qqcAqOi1XQlcaWbLgTuATPq44jeztcBagMLCwiiXNTj1zR384ZNvkZOVzt/duUiraolIXBuW4Zzu/jLw8gUesw5YB1BcXOxDX1WfdVB58gxvHG6k5EgjbxxppLy+FYDHv3wFuTmZQZUmIhIV0Q7+KmBar+2pkbaYFQ47ZfUt/PpwIyWRsK8+1Q7A6Kw0rigaz13F07h+dq769UUkIUQ7+EuA2WY2g57Avxu4J8r7iIpTZ7r49s/e4dWyBpraugDIy8lk2YzxH3zNycshRZ/IFZEEM+DgN7MngeXARDOrBB5090fN7D5gK5AKrHf3vVGpNMr+384KnttTzZ0fmcqVkaAvHD9S/fcikvAGHPzuvrqP9s3EwZDNZ/dUc2nBaL7/PxYFXYqIyLBKyikbKhrb2F3RxC2XTwm6FBGRYZeUwf/cnmoAPn3Z5IArEREZfkkZ/M/uPs6SwrFMGz8y6FJERIZd0gX/ofoW9lWfVjePiCStpAv+53ZXY6ZuHhFJXkkV/O7Os3uOc0XReCaNyQq6HBGRQCRV8B+obaasroVbF6mbR0SSV1IF/3O7q0kxWKlVs0QkiSVN8Ls7z+05zkdnTWTiKE20JiLJK2mCv7TqNEdOtHHrIr2pKyLJLWmC/7k9x0lLMW5eqG4eEUluSRH8Pd081Vw/J5exIzOCLkdEJFBJEfy7jjVR1XSGWy5XN4+ISFIE/7O7j5ORlsJNC/KDLkVEJHDDEvxmNtPMHjWzZ4Zjf711h53N71Tzsbm55GSlD/fuRURizoCD38zWm1mdmZWe1b7CzA6YWZmZ3Q/g7uXuvmawxQ7EG4cbqWvu0Ie2REQiBnPF/ziwoneDmaUCDwMrgQXAajNbMIh9DNpze44zIj2VG+flBVmGiEjMGHDwu/t2oPGs5mVAWeQKvxPYAKwaRH2DEuoO80JpDR+fn8fIjGgvLywiEp+i3cdfAFT02q4ECsxsgpk9AiwxswfO9UQzW2tmO81sZ319fVSKee3QCRpbO9XNIyLSy7BcBrv7CeDeCzxmHbAOoLi42KOx3+f2HCcnM40b5uRG4+VERBJCtK/4q4BpvbanRtqGXWcozJbSGm5amE9WemoQJYiIxKRoB38JMNvMZphZBnA3sCnK++iXHQfrOd0e4lattCUi8lsGM5zzSeB1YK6ZVZrZGncPAfcBW4H9wNPuvjc6pV6cZ3cfZ+zIdK65ZGIQuxcRiVkD7uN399V9tG8GNg+4oiho7+rmxX213LpoChlpSfHhZBGRfkvIVHz5QB2tnd1aUF1E5BwSMvif3V3NxFEZXDVzfNCliIjEnIQL/taOEL98t5aVl04mLTXhDk9EZNASLhm37a+lvSusD22JiPQh4YL/uT3VTBqdRfH0cUGXIiISkxIq+E+d6eKVA/V8+vLJpKRY0OWIiMSkhAr+5vYublqQz2fUzSMi0qeEmrJy6riRPPz5pUGXISIS0xLqil9ERC5MwS8ikmQU/CIiSUbBLyKSZBT8IiJJRsEvIpJkFPwiIklGwS8ikmTMPSrrmkeVmdUDR89qngg0BFDOUNNxxZ9EPbZEPS5I3GM7+7imu3vuhZ4Uk8F/Lma2092Lg64j2nRc8SdRjy1RjwsS99gGelzq6hERSTIKfhGRJBNPwb8u6AKGiI4r/iTqsSXqcUHiHtuAjitu+vhFRCQ64umKX0REoiDmg9/MVpjZATMrM7P7g64nmszsiJm9Y2Zvm9nOoOsZKDNbb2Z1Zlbaq228mb1oZgcj3+NyLcw+ju3Pzawqct7eNrNPBVnjQJjZNDN7ycz2mdleM/tGpD2uz9t5jiuuz5mZZZnZG2a2O3Jc34m0zzCzX0fy8Skzy+jX68VyV4+ZpQLvATcBlUAJsNrd9wVaWJSY2RGg2N3jenyxmV0PtABPuPulkba/BRrd/XuRf7DHufu3gqxzIPo4tj8HWtz9+0HWNhhmNhmY7O67zCwHeBO4Dfgd4vi8nee47iKOz5mZGZDt7i1mlg68CnwD+Caw0d03mNkjwG53/5cLvV6sX/EvA8rcvdzdO4ENwKqAa5KzuPt2oPGs5lXATyK3f0LPH1/c6ePY4p67V7v7rsjtZmAmWNEZAAACMklEQVQ/UECcn7fzHFdc8x4tkc30yJcDNwLPRNr7fb5iPfgLgIpe25UkwEnsxYFfmNmbZrY26GKiLN/dqyO3a4D8IIsZAveZ2Z5IV1BcdYeczcyKgCXAr0mg83bWcUGcnzMzSzWzt4E64EXgENDk7qHIQ/qdj7Ee/InuWndfCqwEvh7pVkg43tOfGLt9ihfvX4BZwGKgGvj7YMsZODMbBfwn8Efufrr3ffF83s5xXHF/zty9290XA1Pp6Q2ZN9DXivXgrwKm9dqeGmlLCO5eFfleB/yMnpOZKGoj/a3v97vWBVxP1Lh7beSPMAz8iDg9b5G+4v8E/t3dN0aa4/68neu4EuWcAbh7E/AScDUw1szSInf1Ox9jPfhLgNmRd64zgLuBTQHXFBVmlh158wkzywY+CZSe/1lxZRPwpcjtLwH/FWAtUfV+MEbcThyet8ibhY8C+939H3rdFdfnra/jivdzZma5ZjY2cnsEPQNe9tPzD8CdkYf1+3zF9KgegMiwq4eAVGC9u/9VwCVFhZnNpOcqHyAN+I94PTYzexJYTs9MgbXAg8DPgaeBQnpmWr3L3ePuTdI+jm05PV0GDhwBfq9Xv3hcMLNrgR3AO0A40vwn9PSHx+15O89xrSaOz5mZXU7Pm7ep9FywP+3ufxHJkQ3AeOAt4Avu3nHB14v14BcRkeiK9a4eERGJMgW/iEiSUfCLiCQZBb+ISJJR8IuIJBkFv4hIklHwi4gkGQW/iEiS+f9F+MTYKqOUfQAAAABJRU5ErkJggg==\n",
+ "text/plain": [
+ "<Figure size 432x288 with 1 Axes>"
+ ]
+ },
+ "metadata": {
+ "needs_background": "light"
+ },
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "conds = []\n",
+ "for N in range(1, len(U)):\n",
+ " conds.append(np.linalg.cond(M[:N, :N]))\n",
+ "plt.figure()\n",
+ "plt.semilogy(range(1, len(U)), conds)\n",
+ "plt.show()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Problem 3"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "In this excercise we will consider the following greedy algorithm:"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### Strong Greedy Algorithm\n",
+ "**Input:** Set of vectors $\\mathcal{M} \\subset V$, target basis size $N$ \n",
+ "\n",
+ "**Output:** Basis $b_1,\\dots,b_N$\n",
+ "\n",
+ "**for** $n := 1$ **to** $N$:\n",
+ "\n",
+ "$\\qquad b_n := \\operatorname{argmax}_{m\\in \\mathcal{M}} \\|m - P_{\\operatorname{span}(\\{b_1,\\dots,b_{n-1})\\}}(m)\\|$\n",
+ " "
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Problem 3 a)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Implement the strong greedy algorithms for a `VectorArray` `V` containing the vectors in $\\mathcal{M}$, a given `Operator` representing the inner product on $V$ and a target basis size `N`.\n",
+ " \n",
+ " "
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Solution"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "To compute the norm of vectors in a `VectorArray` with respect to an arbitrary inner product, we can use the `norm` method of the array:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 16,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "def greedy(M, product, N):\n",
+ " basis = M.space.empty()\n",
+ " max_errors = []\n",
+ " for n in range(N):\n",
+ " lhs = basis.gramian(product)\n",
+ " rhs = basis.inner(M, product)\n",
+ " try:\n",
+ " m = np.linalg.solve(lhs, rhs)\n",
+ " except np.linalg.LinAlgError:\n",
+ " break\n",
+ " M_proj = basis.lincomb(m.T)\n",
+ " M_err = M - M_proj\n",
+ " errors = M_err.norm(product)\n",
+ " max_errors.append(np.max(errors))\n",
+ " basis.append(M[np.argmax(errors)])\n",
+ " return basis"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Problem 3 b)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Apply your algorithm to the `VectorArray` `U` with `N=30`. For the resulting basis compute the projection error and condition numbers as in Problem 2 (d) and (e). Compare your results with the results from Problem 2."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Solution"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "First, we compute a reduced basis using the just implemented greedy method:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 17,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "greedy_basis = greedy(U, m.h1_0_product, 30)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Now we compute the approximation error as usual:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 18,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAD+CAYAAAAzmNK6AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xl4VOX5//H3M5N9IRtJWBMgAbKwBAhLjCyCLC4otVgFq4h1oe7tt/6qbS3aarXVtu5alUWtFQXrClZUQJA9UfY1QIBgIAkhZE8myfP740xiQAKTZJIzM7lf1zVXJicz59yHucgn5zyb0lojhBCi47GYXYAQQghzSAAIIUQHJQEghBAdlASAEEJ0UBIAQgjRQUkACCFEByUBIIQQHZQEgBBCdFBe7XUgpVQg8BJQDazSWr/dXscWQgjxY626AlBKzVdK5Smldpy1fYpSaq9SKksp9aB98zXAEq31bcBVrTmuEEKI1mvtLaCFwJTGG5RSVuBF4DIgCZihlEoCegBH7S+rbeVxhRBCtFKrbgFprVcrpXqdtXkEkKW1PgiglFoEXA3kYITAFs4TPEqp24HbAQIDA4clJCS0pkQhhOhQMjMzC7TWkY68ti3aALrzw1/6YPziHwk8B7yglLoC+KSpN2utXwVeBUhNTdUZGRltUKIQQngmpdRhR1/bbo3AWusyYHZ7HU8IIcT5tUU30GNAz0bf97BvE0II4ULaIgA2A32VUr2VUj7A9cDHzdmBUmqqUurV06dPt0F5QgghoJW3gJRS7wDjgM5KqRxgrtZ6nlLqbuBzwArM11rvbM5+tdafAJ+kpqbe1pr6hBAXZrPZyMnJobKy0uxSRDP4+fnRo0cPvL29W7yP1vYCmtHE9mXAstbsWwjRPnJycggODqZXr14opcwuRzhAa83JkyfJycmhd+/eLd6PTAUhRAdXWVlJRESE/PJ3I0opIiIiWn3VJgEghJBf/m7IGZ+ZSwZAaxqB62rryFj8N7b+b0EbVCaEEJ7DJQNAa/2J1vr2kJCQ5r9XKUL3LMJr44uUV9e0QXVCiPZ2+eWXU1RURFFRES+99FLD9lWrVnHllVc65RirVq1i3bp1TtlXS2VnZzNgwIB2O55LBkBrWC0K38E/JVnv5+3/fWN2OUIIJ1i2bBmhoaE/CgBnam0AaK2pq6tzYkVtz+MCAKDnxUbnpJObF3O0sNzkaoQQ5/PUU0/x3HPPAfCrX/2K8ePHA7BixQpuuOEGAHr16kVBQQEPPvggBw4cICUlhQceeACA0tJSpk+fTkJCAjfccANaawC++uorhgwZwsCBA7nllluoqqo6Y18AGRkZjBs3juzsbF555RX++c9/kpKSwpo1a86oMT8/n4kTJ5KcnMytt95KbGwsBQUFZGdn079/f2666SYGDBjA0aNHWb58OWlpaQwdOpRrr72W0tJSADIzMxk7dizDhg1j8uTJ5ObmNmwfPHgwgwcP5sUXX2w45pgxY9iyZUvD9xdffDFbt2516r99u00F0a7C+2CLGsRlJzbwl2W7efnnw8yuSAi38OgnO9n1fbFT95nUrRNzpyY3+fPRo0fz97//nXvvvZeMjAyqqqqw2WysWbOGMWPGnPHaJ598kh07djT8Yly1ahXfffcdO3fupFu3bqSnp7N27VpSU1O5+eab+eqrr+jXrx833XQTL7/8Mvfff/85a+jVqxdz5swhKCiI3/zmNz/6+aOPPsr48eN56KGH+N///se8efMafrZ//37eeOMNRo0aRUFBAY899hhffvklgYGB/PWvf+Uf//gHDz30EPfccw8fffQRkZGRvPvuu/z+979n/vz5zJ49mxdeeIExY8Y0hBrAL37xCxYuXMgzzzzDvn37qKysZPDgwc36t78Ql7wCcMZIYO+BP2GwymLbju2sO1DgxOqEEM40bNgwMjMzKS4uxtfXl7S0NDIyMlizZg2jR4++4PtHjBhBjx49sFgspKSkkJ2dzd69e+nduzf9+vUDYNasWaxevbrFNX7zzTdcf/31AEyZMoWwsLCGn8XGxjJq1CgANmzYwK5du0hPTyclJYU33niDw4cPs3fvXnbs2MHEiRNJSUnhscceIycnp6Fdoz7obrzxxob9XnvttXz66afYbDbmz5/PzTff3OL6m+KSVwBOGQmcPA2+epQZQd/x6Md9WHrvxXhZXTLvhHAZ5/tLva14e3vTu3dvFi5cyEUXXcSgQYNYuXIlWVlZJCYmXvD9vr6+Dc+tVis1Nefv/OHl5dVwr94Zo58DAwMbnmutmThxIu+8884Zr9m+fTvJycmsX7/+jO1FRUVN7jcgIICJEyfy0Ucf8d5775GZmdnqWs/mub8Rw/tAl0H8vNN37D1Rwn82HTG7IiFEE0aPHs3TTz/NmDFjGD16NK+88gpDhgz5UV/34OBgSkpKLri//v37k52dTVZWFgBvvfUWY8eOBYzbPfW/TN9//32H9p2ens57770HwPLlyzl16tQ5Xzdq1CjWrl3bcNyysjL27dtH//79yc/PbwgAm83Gzp07CQ0NJTQ0lG++MTqsvP32mSvl3nrrrdx7770MHz78jKsOZ/HcAABInkZo4Vamxtby9+X7OFVWbXZFQohzGD16NLm5uaSlpREdHY2fn985b/9ERESQnp7OgAEDzrhffjY/Pz8WLFjAtddey8CBA7FYLMyZMweAuXPnct9995GamorVam14z9SpU/nggw/O2Qg8d+5cli9fzoABA1i8eDFdunQhODj4R8eNjIxk4cKFzJgxg0GDBpGWlsaePXvw8fFhyZIl/Pa3v2Xw4MGkpKQ09DhasGABd911FykpKQ0N2PWGDRtGp06dmD27bWbSV2cf0JW0ekGYkwfg+aHkpf2RtK8TmTkihj9Pa78+tkK4g927dzt0q6Ujq6qqwmq14uXlxfr16/nlL395Rg+dtvL9998zbtw49uzZg8Xy47/Xz/XZKaUytdapjuzfs68AIuKgy0Cijn7Gz0fG8PbGw+zOdW4PByGE5zty5AjDhw9n8ODB3Hvvvbz22mttfsw333yTkSNH8vjjj5/zl78zuOQVgFJqKjA1Pj7+tv3797duZ6ufhhV/5vScLYx7dT/9uwTzzm2jZO4TIezkCsB9eeQVQGumgviR5J8AEHJoGb+e1J8NBwv5bMfx1u9XCCHcnEsGgFNFxEH0QNj5ITNHxJDQJZjHl+6m0lZrdmVCCGEqzw8AMMYE5GzCWnKMR65K5lhRBf/6+qDZVQkhhKk6SAAYt4HY9TGj+kRwxcCuvPx1FseKKsytSwghTNQxAqD+NtCuDwF46PIEtIYnlu02uTAhhCPaYzpoZ1i4cCF333232WU4rGMEAEDy1XB0I5zOoUdYAHPGxvHptlw2HjxpdmVCiAtoj+mgLzSFhCdyyQBwxmRwP5L0w20ggDlj4+gW4sfcj3eSebiQ6hr3msdbCE/RVtNBN7Z582YGDRrU8L76RVcWLlzIVVddxfjx45kwYUJDPcOHD2fQoEHMnTu3YR///ve/GTFiBCkpKdxxxx3U1hodSRYsWEC/fv0YMWIEa9euBaCkpITevXtjs9kAKC4uPuN7V+G5k8GdrXM8RA8wbgOl3Ym/j5VHrkrml29/y09fXo+/t5VhsWGM7B3OyD4RDO4Zgq+X9cL7FcKTfPYgHN/u3H12GQiXPdnkj9tiOuiLL774jPfNnj2b1157jbS0NB588MEzfvbtt9+ybds2wsPDWb58Ofv372fTpk1orbnqqqtYvXp1wxTOa9euxdvbmzvvvJO3336biRMnMnfuXDIzMwkJCeGSSy5hyJAhBAcHM27cOJYuXcq0adNYtGgR11xzDd7e3k76R3UOlwyANpM0DVY+BqePQUh3JiV3YdPvJrA5u5ANBwvZcPAkf/9iHwC+XhaGxIQysncEI/uEMzQmDD9vCQQhnO3s6aCHDh3aMB10/ZXB+dRPBw00TAfdOACKioooKSkhLS0NgJkzZ/Lpp582/HzixImEh4cDxkRvy5cvZ8iQIYBxdbF//362bdtGZmYmw4cPB6CiooKoqCg2btzIuHHjiIyMBOC6665j3z7jd8itt97K3/72N6ZNm8aCBQvaZfRwc3WsAEi2B8Duj2HULwGICPJlyoCuTBnQFYCi8mo2HTICYeOhkzy3Yj/6K/CxWnj4ykRuTOtl4gkI0cbO85d6W2nv6aDPdvZ0zg899BB33HHHGa95/vnnmTVrFk888cQZ2z/88MMm95uenk52djarVq2itra2Xdf6dZRLtgG0mc59ISoZdjb9oYUG+DApuQt/nJrE0ntHs+WPk5g3K5W+0UG8uubgOe8vCiFax9nTQTcWGhpKcHAwGzduBGDRokVNvnby5MnMnz+/YRnHY8eOkZeXx4QJE1iyZAl5eXkAFBYWcvjwYUaOHMnXX3/NyZMnsdlsLF68+Iz93XTTTcycObPNZvNsrY4VAGCMCTi6AYq/d+jlIf7eTEiMZsaIGI4WVnAgv7SNCxSi43H2dNBnmzdvHrfddhspKSmUlZXR1DQzkyZNYubMmaSlpTFw4ECmT59OSUkJSUlJPPbYY0yaNIlBgwYxceJEcnNz6dq1K4888ghpaWmkp6f/6Irlhhtu4NSpU8yYMaN5/yDtxCUng6vX6umgz6VgP7yQClP+CqPmOPy2Y0UVpD+5gocuS+COsXHOrUkIE3WEyeBKS0sJCgoCjIbk3Nxcnn322TY/7pIlS/joo49466232mT/rZ0MrmO1AcAPt4F2fdisAOge6k9Cl2BW7MmTABDCzSxdupQnnniCmpoaYmNjWbhwYZsf85577uGzzz5j2bJlbX6slnLJAGg0HXTbHCB5Gqz8CxTnQqeuDr9tfEIU/1p9kNMVNkL8Xas7lxCiaddddx3XXXddux7z+eefb9fjtYRLtgE4dTroc0maBmijN1AzjE+IorZOs3pfftvUJYRJXPlWsDg3Z3xmLhkAbS6yH0Qlnbc30LkMiQkjNMCblXvy2qgwIdqfn58fJ0+elBBwI1prTp48iZ+fX6v245K3gNpF0jRY9QSUHIfgLg69xWpRjOsXyap9+dTWaawWWVVMuL8ePXqQk5NDfr5c2boTPz+/hgFwLdVxAyB5Gqz6izE30MjbHX7bJQlRfLjle7YcLWJYbFgbFihE+6gfiCU6no55Cwggsj9EJsLOD5r1trH9IrEo5DaQEMLtddwAAGNQ2JH1xtxADgoN8CE1NpwVEgBCCDfXsQNg8HWgFGxu3iRNlyREsSu3mOOnK9uoMCGEaHsdOwDCekHCFZCxAKrLHH7b+IQoALkKEEK4tY4dAABpd0NlEWz5j8Nv6RcdRPdQfwkAIYRbkwDoORK6D4MNL0OdY6uCKaUYnxDF2qwCKm21bVygEEK0DZcMgDZZErLpg8GoO6HwAOz/3OG3jU+IosJWywZZU1gI4aZcMgDafCqIsyVdDZ16wPoXHX5LWlwEft4W6Q4qhHBbLhkA7c7qbQwGy14DuVsdeouft5X0uM6s2JsnQ+iFEG5JAqDe0FngHQjrX3L4LZckRHG0sIKsPFkkRgjhfiQA6vmHwpCfw473jWmiHXCJdAcVQrgxCYDGRs2BuhqHB4Y1XiRGCCHcjQRAY+F97APD5kN1uUNvGZ8QRcbhU5yusLVxcUII4VwSAGdLuwsqTsHWdxx6uSwSI4RwVxIAZ4tJg25DYMNLDg0Mk0VihBDuSgLgbErBqLvgZBZkfXHBl5+9SIwQQrgLCYBzSZ4Gwd1g/QsOvfyShCgKy6rZcrSojQsTQgjnkQA4F6s3jLwDDq2G49sv+PKx/SKxWpTcBhJCuBUJgKYMmwXeAQ4NDAsN8GFYTJh0BxVCuBUJgKb4hxkDw7YvNhaOv4D6RWJyT1e0Q3FCCNF6EgDnM7J+YNjrF3zphERjVPDKPdIdVAjhHlwyANp1OujziYiD/pfD5nlgO/9f9n2jZJEYIYR7cckAaPfpoM8n7U6oKISti877MlkkRgjhblwyAFxKbDp0HezQwLDxibJIjBDCfUgAXIhSxrrBBfvgwFfnfWlaH1kkRgjhPiQAHJE0DYK6QObC876sfpGYr/bIIjFCCNcnAeAILx/oPwUOfg2155/1c1JyNDmnKkh/cgWPfLyTdQcKqKl1bLF5IYRoT15mF+A24iYYVwBHN0Gv9CZfdu2wnnhbLXy24zjvbDrCwnXZhAV4c2liNJOTu3Bx3874eVvbr24hhGiCBICjeo8BZTXaAc4TABaL4pqhPbhmaA/Kq2v4em8+n+88zv92HmdxZg4BPlbG9Y9kcnIXLkmIopOfdzuehBBC/EC58r3q1NRUnZGRYXYZP5g3CWqq4I6vm/3W6po61h88yec7j7N85wkKSqvwtiou6R/FzJExjOkbicWi2qBoIURHopTK1FqnOvJauQJojrgJsOoJKCuAwM7NequPl4Wx/SIZ2y+SP189gO+OnOKzHcf58LtjLN91gh5h/swYEcO1qT2ICvZroxMQQogfyBVAc+RkwOsT4KfzYOB0p+yyqqaW5TtP8J+NR1h/8CReFsXEpGhmjowhPa6zXBUIIZpFrgDaSrchxiRxWV85LQB8vaxMHdyNqYO7cTC/lHc2HWFxZg6f7ThObEQAM0bEMH1YDzoH+TrleEIIUU+uAJpr8c1weD383x5jkFgbqLTV8r8dx/nPxiNsyi7E26q4bEBX7p3Ql/iooDY5phDCM8gVQFuKmwA7P4ATO6HLgDY5hJ+3lWlDujNtSHf2nyjhP5uOsDgjh6Xbc7l+eE/uv7QfkcFyRSCEaB0ZCNZcceONrxeYFsJZ+kYHM3dqMqseGMcNI2N4d/NRxj21kme/3E95dU271CCE8EwSAM0V0h0iE412gHbUOciXP109gOW/GsPovpH888t9jHtqFYs2HZGRxkKIFpEAaIn4CXBkPVSXtfuh+0QG8cqNw1gyJ40eYf48+N/tXP7cGlbsOSHzDwkhmkUCoCXixkNtNWSvNa2E1F7hvP/Li3jl50Ox1WpuWZjBzNc2sj3H5EV0hBBuQwKgJWIvAi+/dmsHaIpSiikDurL8V2P409XJ7D1RwtQXvuHRT3ZSWydXA0KI85MAaAlvfyMEDqwwuxIAvK0WbkrrxdcPjGNWWiwL1mZzzzvfyspkQojzkgBoqbgJxiIxRUfNrqRBsJ83j149gD9ckciy7ceZNX8TxZXnn75aCNFxSQC0VPwE46vJt4HO5dbRfXj2+hS+PXKKn72ynhPFlWaXJIRwQRIALRWZAMHd2r07qKOuTunO/JuHc7SwnGteWseB/FKzSxJCuJh2CwClVB+l1Dyl1JL2OmabUgrix9tXCXPNAVmj+0ay6PY0qmpqmf7yOr47csrskoQQLsShAFBKzVdK5Smldpy1fYpSaq9SKksp9eD59qG1Pqi1/kVrinU5cROg6jQcyzS7kiYN7BHC+7+8iE7+3sx8baMsWC+EaODoFcBCYErjDUopK/AicBmQBMxQSiUppQYqpT496xHl1KpdRZ9xoCwu2Q7QWGxEIEvmXERcVCC3vpnB4gzXabgWQpjHoQDQWq8GCs/aPALIsv9lXw0sAq7WWm/XWl951sPhPzuVUrcrpTKUUhn5+fkOn4gpAsKh21CXbQdoLDLYl0W3p5HWJ4IHlmzjxZVZMnJYiA6uNW0A3YHGf0rm2Ledk1IqQin1CjBEKfVQU6/TWr+qtU7VWqdGRka2orx2Ej8Bvv8Wys/OR9cT5OvF/JuHc9Xgbjz1+V4e/WSXhIAQHVi7NQJrrU9qredoreO01k+013HbXNwE0HVwcJXZlTjEx8vCM9elMDu9FwvXZbN0e67ZJQkhTNKaADgG9Gz0fQ/7to6l+zDwDXGZUcGOsFgUf7giicSunXhi2R4ZMSxEB9WaANgM9FVK9VZK+QDXAx87oyil1FSl1KunT7vBxGZWL+gz1ggAN7qdYrUoHr4ykWNFFcz75pDZ5QghTOBoN9B3gPVAf6VUjlLqF1rrGuBu4HNgN/Ce1nqnM4rSWn+itb49JCTEGbtre/EToPgY5O81u5JmuSiuM5OTo3lxZRZ5MlpYiA7H0V5AM7TWXbXW3lrrHlrrefbty7TW/ez39R9v21JdWDuvEuZMv7s8EVttHU997l7hJYRoPZkKwhlCYyCir1t0Bz1bbEQgt6T3Zsm3ObKWgBAdjASAs8RPgMNrwVZhdiXNdtf4eMIDfPjTpzulW6gQHYhLBoBbNQLXi5sANZVweJ3ZlTRbJz9v/m9SfzZnn2LZ9uNmlyOEaCcuGQBu1wgM0CsdrD5u1R20seuG9yShSzB/WbZbuoUK0UG4ZAC4JZ9AiElzy3YAMLqF/vHKJOkWKkQHIgHgTPETIH83nHbP8XAXxXdmYlI0L0m3UCE6BAkAZ4qzrxJ2cKW5dbTC7y9PpLq2jqeXS7dQITydSwaAWzYCA0QnQ1AXt70NBNCrcyCz03uzODOHHcfc7N9fCNEsLhkAbtkIDMYqYXHjjSuAOvdtSL27vluozBYqhEdzyQBwa/0mQ8Up2PWR2ZW0WCc/b349qR+bsgv5bId0CxXCU0kAOFviVIgeAF/MBZv7NqRelyrdQoXwdBIAzmaxwuS/wOkjsOFFs6tpMS+rhYevTCLnVAXz10q3UCE8kQRAW+gzFvpfAWv+ASUnzK6mxdLjO3NpYjQvrsgir8R9r2aEEOfmkgHgtr2AGpv0Z6ipghV/NruSVvn9FUa30L9/vs/sUoQQTuaSAeC2vYAai4iDkXfAd/+G3G1mV9NivTsHclNaL97LPMru3GKzyxFCOJFLBoDHGPMA+IfB579zq9XCznbP+HiCfb144rM9ZpcihHAiCYC25B8Kl/wOstfAnqVmV9NioQE+3DuhL6v35fP1vnyzyxFCOIkEQFsbNhsiE2D5H4w2ATd1Y1osMeEB/GXpbmrr3PdqRgjxAwmAtmb1gsmPw6lDsOlVs6tpMV8vK7+dksDeEyUsyTxqdjlCCCeQAGgP8ZdC/ET4+ikoKzC7mha7fGAXhsaE8vfl+yirqjG7HCFEK7lkAHhEN9CzTX4cqkth5V/MrqTFlFL8/ook8kqqeHX1QbPLEUK0kksGgEd0Az1bZH8Y/gvIXAAndpldTYsNiw3jioFdeXX1QU7ImgFCuDWXDACPNe4h8A2G5b93626h/29Kf2rq6vjHchkcJoQ7kwBoTwHhMPZBY93g/V+YXU2LxUbI4DAhPIEEQHsbfiuExxmDw2ptZlfTYveMj6eTnzd/Wbbb7FKEEC0kAdDevHyMBuGT+yFjvtnVtFhogA/3jI9nzf4CGRwmhJuSADBDvynQe6zRI6i80OxqWkwGhwnh3iQAzKCUsWZAVTG8MwMOrnLLRmFfLysPXmYMDlucIYPDhHA3LhkAHjkO4GxdBsCVz0DhQXjzavjXGNi22O3aBS4bYB8c9oUMDhPC3bhkAHjkOIBzGTYL7t8OVz0PNZXw31vh2RRY9wJUukfvmvrBYfkyOEwIt+OSAdChePvB0Jvgzo0w8z0I62WME/hnMix/GIq/N7vCCxoWG8YVg2RwmBDuRgLAVVgs0G8yzF4Kt60w5g9a/wI8MxA+mAMndppd4Xn9dnICNXV1/PWzPeQVV1JSaZOGYSFcnNIu3PiYmpqqMzIyzC7DPKeyYcMr8O2bUFNhtBkMm2V2VU167NNdvP7NmQvI+3lbCPTxIsDXanz1sRLoa3y9NDGa6cN6oJQyqWIhPI9SKlNrnerQayUA3EB5Ifz3dsj6AiY9BhfdY3ZF52SrrePLXSc4VW6jvLqGsqpa42t1DeVVtcbX6lrKqmooKK3mSGE5o/t25smfDqJ7qL/Z5QvhESQAPFFNNXxwO+z8AEb/Bsb/wehO6qbq6jRvbzrCE8t2o4DfXZHIzBExcjUgRCs1JwCkDcBdePnAT+cZDcZrnoZlD0BdndlVtZjForhxVCyf3z+GlJhQfv/BDn4+byNHC8vNLk2IDkMCwJ1YrDD1OUi7Gza/Bh/OgVr37nvfMzyAf/9iJH/5yUC2Hj3N5GdW89b6bOqkAVmINicB4G6UMtoBxv8Btr0L790ENvfueqmUYubIGD7/1RiGxYbx8Ec7mfn6Bg6fLDO7NCE8mgSAO1IKxjwAlz0Fe5fCf66FqhKzq2q17qH+vHnLCP7200HsPFbMlGfWsHDtIbkaEKKNuGQAdIipIJxh5O0w7RXIXgtvTnPrieXqKaX42fCeLP/1GEb1CeeRT3Yx8/UNVNpqzS5NCI/jkgHQYaaCcIaUGfCzN+H4Nlh4BZQcN7sip+ga4s/8m4fz2LQBbDhYyKJNR8wuSQiP45IBIJop8Uq4YTGcOgzzp8C+5cbI4dJ8t+4ppJTihpExjOgVzstfH5CrACGczMvsAoST9BkHN30Eb0832gTqKSsEdoagKAiMgqBoCIo0vkb0hX6TzKrYIUop7ru0Lze8vpH3Mo5yU1ovs0sSwmNIAHiSnsPh3u8gbzeU5UFp/eMElOUbX/P3Gl/r7NNOz/4fxKaZW/cFXBQXQWpsGC+vOsB1w3vi62U1uyQhPIIEgKcJCIde6ed/jdZGMLww3BhP4OIBoJTi/kv78fN5G3kvI4cbR8WaXZIQHkHaADoipSA4GlJmwq6PoeSE2RVdUHp8BMNiw3h5ZRZVNdIWIIQzSAB0ZMNvNW4Fffum2ZVckFKK+yb05fvTlSzJzDG7HCE8ggRAR9Y5HvpcApkL3GJKidF9OzMkJpSXVh6gusZ9ezcJ4SokADq6EbdB8THYu8zsSi6o/irgWFEF738rVwFCtJYEQEfXdzJ06gGbXze7EoeM7RfJ4J6hvLgyS64ChGglCYCOzuoFqbPh0NeQv8/sai5IKcX9E/qSc6qC/8pVgBCtIgEgYOgssPq4zVXAuP6RDOoRwgsrs7DVylWAEC0lASCMkcFJ02DrO1BVanY1F2SMCzCuAj749pjZ5QjhtiQAhGH4rVBVDNvfM7sSh1zSP0quAoRoJZcMAJkO2gQ9R0CXgbDpdWOksItTSnHv+L4cKSznw+/kKkCIlnDJAJDpoE2gFAy/DfJ2wpH1ZlfjkAmJUQzo3okXVmZRI1cBQjSbSwaAMMnAa8E3xG0ag+uvAg7YA5KHAAATq0lEQVSfLOejLd+bXY4QbkcCQPzAJwCG3OA28wMBTEyKJqmrXAUI0RISAOJMDfMDvWF2JQ5RSnHvhL4cKijjk21yFSBEc0gAiDNFxEHceMhwj/mBACYlRZPQJZjnV2RRKwvIC+EwWQ9A/Njw22DRDGN+oKSrzK7mgiwWY46gX779Lb96dwvdQv3x8bLga3/4eFnwsdq/2p+H+HszJCYMHy/5G0h0XBIA4sf6TYaQnsZiMW4QAACTk7swKSmaFXvyqK6po9qB9oAgXy/G9o9kYmI04/pHEhrg0w6VCuE6JADEj1msxvxAX/3JWEIysr/ZFV2QxaJ49abUhu+11lTX1hlhYA+EKltdw7bc05Ws2HOCL3fnsXRbLlaLYnivMC5NjGZiUjSxEYEmno0Q7UNpFx70k5qaqjMyMswuo2MqK4B/JMKw2XD538yups3U1Wm2HTvNl7tO8OXuE+w5XgJA36ggLk2K5tLEaIb0DMViUSZXKoRjlFKZWuvUC79SAkCcz39vh72fwa93g2+Q2dW0i6OF5Xy52wiDjQcLqanT3D6mD7+7PNHs0oRwSHMCQFrARNPq5wfa9q7ZlbSbnuEBzE7vzdu3jiLz4YlcMbArb60/TFF5tdmlCeF0EgCiaT2GQ5dBxshgF75SbCsh/t7cPT6eClstizYfNbscIZxOAkA0TSljyci8XbBnKVQUdbggSOzaibQ+Eby5LltGGguPI72AxPkNmA5f/BHevcH43uIF/uEQ2BkCIiAg3P7V/giMhOAu0KkbBHcDbz9z63eC2em9uP2tTD7feYIrBnU1uxwhnEYCQJyfTwDM/h/kboHyk2c+yk5C3h7jeUUh6HP8hewfbg+DrsbX+kdwN+iRCv6h7X9OzTQhMZqY8AAWrD0kASA8igSAuLCoBONxPnV1UFkEZflQ/D2U5ELxMeN5sf157hbj5/V8O8HIO2DUncaVhIuyWhSzLurFnz/dxdajRQzu6fqhJYQjpBuoaF81VVByHE5lG43Luz8Gn2AYeTuk3e2yQVBSaSPtiRVcmhjFM9cPMbscIZok3UCF6/LyhbBY6DMWrnsLfrkO4ifAmn/AMwPhi7nGIDQXE+znzbWpPVi6PZcTxZVmlyOEU0gACHNFJ8PP3oA710PfSbD2WXhmECx/GErzL/z+dnTzRb2oqdP8e8Nhs0sRwikkAIRriEqEaxfAnRsg4XJY/wI8Owg+/z2U5pldHQCxEYFMSIjm7Y1HqLTVml2OEK0mASBcS1QC/PR1uGsTJE6FDS/BP5Ph39Mhc6HpYXBLei8Ky6r5WJagFB5AGoGFayvIgswFsPsTKDoMKOg5EhKvhIQrIbx3u5ajteayZ9cA8Nl9o1FKJokTrsUlG4GVUtOUUq8ppd5VSk1qr+MKN9c5HiY/DvdthTlrYdyDUF0Gy/8Az6XAy+mw8gk4vr1dRikrpZid3os9x0tYf/Bkmx9PiLbk0BWAUmo+cCWQp7Ue0Gj7FOBZwAq8rrV+0oF9hQFPa61/caHXyhWAaFLhIWN6ij1L4ch6QENoLAz5uTF9hX9Ymx260lbLRU+uYGhMGK/PcugPLSHaTVtcASwEppx1ECvwInAZkATMUEolKaUGKqU+PesR1eitf7C/T4iWC+8NF90Nt3wGv9kPU58z1jNe+Tj8c4AxfUUbtRf4eVuZOSKGr/ac4PDJsjY5hhDtwaEA0FqvBgrP2jwCyNJaH9RaVwOLgKu11tu11lee9chThr8Cn2mtv23qWEqp25VSGUqpjPx81+oGKFxUUCQMmwU3fmDcJuo3GdY9b4wrWPYAFDl/Js8b02KxKsXCddlO37cQ7aU1bQDdgcb/s3Ls25pyD3ApMF0pNaepF2mtX9Vap2qtUyMjI1tRnuiQugyA6fPh7gwYOB0y5httBR/eZTQoO0l0Jz+uHNSVxRk5lFTanLZfIdpTuzUCa62f01oP01rP0Vq/0l7HFR1URBxc/SLcuwVSfwE7lsCLw2HxbDi+wymHmJ3em9KqGhZn5Dhlf0K0t9ZMBncM6Nno+x72bUK4jtCexprGY34D61+EzfNg53+h12hj2mrvAOPhEwDe/uAdaH9e//A35icK6mJMdW394b/M4J6hDIsN44312cy6qBdWWTdYuJnWBMBmoK9SqjfGL/7rgZnOKEopNRWYGh8f74zdCQFBUTDxUbj4ftj0Guz6GE7ngK0cqsvBVnbu6azPoIx1EIK6QHA0BHXhyZAA/p1TxY4vjjA4oS/4hYJfiPHwCTQW1RHCRTnaDfQdYBzQGTgBzNVaz1NKXQ48g9ENdL7W+nFnFifdQEW70RpqqxsFQvkPz8tPQulxKDkBpfZHyXEoPYEuzUPpJqaFUFbw6/RDIDR+dBkMCVdAyPmazYRovuZ0A5WRwEK0Rl0dC7/M5N2Vm/nXNTHEBNig8rT9UdzoeaNHRaERIgDdhxkjmhOnQue+5p6L8AgSAEK0o6LyakY98RVXDe7G36YPduxN+XuN6S32fArff2dsi0z4IQy6DpbbR6JFmhMAsiKYEK0UGuDDNUN7sCQjh37RwdyU1gsfrwt0sIvsbzzG/MYYp7BnqREG3/wD1jwNITHGLaKEKyBmFFi92+dkRIfiklcAjRqBb9u/f7/Z5QhxQfklVfxm8Va+3pdPn8hAHr4yiUv6R134jWcrK4C9nxlhcGAl1FYZS2f2GQvxlxqPkB7OPwHhMeQWkBAm0Fqzcm8ef/50N4cKyhifEMUfrkikT2RQy3ZYVWKEQNaXxqPY3ss6MuGHMIi9yFhlTQg7CQAhTFRdU8fCdYd47qssqmpqmZ3em7vHx9PJrxW3cbQ22g2yvjDC4PA6o9eSd4AxpiH+UojoA17+4O1nbPfyM8Yx1H+1+ki7QgcgASCEC8grqeTpz/eyODOHiEAf/t/kBKYP64HFGQPGqssg+xvY/4URCqeyHXiT+iEQrN5GIFi8zvPcB+LGQ+otRqgItyABIIQL2ZZTxKOf7CLz8CkGdg/hkauSGBYb7tyDFB4yZj+tqQCb/VFT2fTXWpvxqLM1/byyGPJ3Q0hPuOR3MOg6sFidW7dwOrcPAGkEFp5Ga83HW7/niWV7OF5cSbcQPwJ8vQjwsRLgYyXQx4sAXy8CfawE+Ni3+1oJC/BhYPcQ+ncJxttqwgquB1bCl49A7haISoIJf4R+U+RWkgtz+wCoJ1cAwtOUVdXwxvpsDuaXUV5dQ3l1LeVVtZTZn5dV1VBRbXxf1+i/pq+XhYHdQ0jpGUpKTCgpPUPpHurfPktS1tXBrg9hxZ+h8CD0HGVMqxEzqu2PLZpNAkAIN6e1pqqmjrziKrbkFLH1aBFbjhax49hpqmqMOYs6B/mS0tMeCj3DiI0IICTAm2Bfr7YJhlobfPsmfP1XYyRzv8uMK4LoJOcfS7SYBIAQHspWW8ee3BK2HD3Fd/ZQOJh/5qpkXhZFaIA3If7ehAX4EBrgTWiAD6H+3oQF2r/392l4Tf3PA32sjgVHdRlseBnWPmt0VU2ZCeMeMmZeFaaTABCiAzldbmPbsSKOn66kqNzGqfJqiipsFJVX27+3cbq8mlPlNipsTUxcx5nBUR8YPcMDSOraiaRunYiPCsLPu1EjcHkhrPm7MbuqssAVTxtrMgtTSQAIIc6p0lZLUbmN0/UBUWHjdLmNogojLBp/X1hmI7ugrCE0rBZFXGQgiV07kdS1E4n2R2RtHnx0JxxaDSk3wOVPG2sqCFO4fQBILyAhXENtnebwyTJ255awO7eYXbnF7M4tJvd0ZcNrOgf5EtfZj+vK3mZayTsc847luYiHyfOJQSlQgFKq4WtcZCCj4iIY3iucIF+ZjszZ3D4A6skVgBCu6VRZdaNAKOFoYTl1WjOoMoP7ip/CBxsvBt3Dat+xaA0ajdZQU6s5WFCKrVZjtSgG9wghLS6CtD6dGRYbhr+PjDNoLQkAIYR5Th+DJbfA0Q3GKOLJT5wxkriiupaMw4WsP3CS9QdPsi3nNLV1Gh+rhZSeoYyKiyCtTwRDYkLPbHMQDpEAEEKYq9YGX/0J1j1nrG1w7RsQ3vucLy2tqmFzdiEb7IGw49hp6rQx9mFE73BG9+3MxfGRJHYNbp9xD25OAkAI4Rr2LIMP54AGpr1oLHZzAacrbGw+VMjaAwV8s7+A/XmlAHQO8uHi+M5c3DeS0X07E91J5ic6FwkAIYTrOHUYFs8yVj4bdReMmgNWX/vEc/WTz3mD5dxTXRw/Xck3WQWs2Z/P2qwCCkqrAegbFcTFfTszum9nLorrLLeL7CQAhBCupaYKlj8Mm/7V9GuU1QgDq48RDBavH805pIGa2jqqa+qoqq3DVlOHBsrxRwWE0SksirCIKCwBYeB/jken7tCpa5ueqtncfknIRt1AzS5FCOEMXr5w+d8g6WooPPDDrKO11cajruaH57X253W2H+1GAd72RyBQU6cpKK6gpKCQ8uJ8SssOUfH9Tjpby/GrLT13LbEXG6OXk64G3xYu1uMh5ApACOERKm21rNqbz0dbjvHVnjxqa2wkhWl+kuDPxN6+9PSrguPbYOs7cDILvAONEEiZCbHpTd6CcjdyC0gI0aEVV9r4fMdxPtryPesOFFCnIblbJy4b0IWE6GCS6/YQffC/WHZ9AFXFEBprBMHg6yGsV5vUVF1Th7dVtXlPJgkAIYSwyyuu5NNtuXy05Rhbc043bPexWugfYeUn/lsYX/klsac3odDUxKTjNeQGiE4Gqw/a6k21tlJWY6G0xkJpjaLUZqHEpii1aUqraiiuqKG40kZxhY3iyhr7V2PKjfqfVdfUERsRwKWJ0UxMiiY1NgyvNljjQQJACCHO4XS5jQMFpRzIK+VAfhkH8ks5kF/K4ZPlRNfl8xPrN0y3fk1vywmH9lerFTV4UYuFOizUKisaC1pZ0Bar0bDd8NWLXB3Gd2WR7KvrynHvGLrHDyJ14ADGJkQ5bVoMCQAhhGiG6po6jhSWG4GQV4LtSAbBNYUEWusI9KrD31pLgFXjb63Fz1KHn6UWX1WDr6rFV9XiY9V4UYfSdUaDtq6Fulqo/76u1mjUPp2DLtiHqippOHaZ9uUQ3SgJ7IV/1wRi+qcQHpNsrMDWgiU43b4XkBBCtCcfLwvxUUHERwVBchegb5sdS2ltrN9csI/a/L0UH9qB//e7iSjeSXTWKiwHjD/KC+/LJjwsrM3qAAkAIYRoX0pBcDQER2PtPZquI4zNWmsO5haQ+V0GJ4/uZ05oaJuX4pIBIOMAhBAdjVKKuG6RxHW7DLisXY7pkh1ftdafaK1vDwkJMbsUIYTwWC4ZAEIIIdqeBIAQQnRQEgBCCNFBSQAIIUQHJQEghBAdlASAEEJ0UBIAQgjRQbn0XEBKqXzgcAvf3hkocGI5ZvO08wHPOydPOx/wvHPytPOBH59TrNY60pE3unQAtIZSKsPRCZHcgaedD3jeOXna+YDnnZOnnQ+07pzkFpAQQnRQEgBCCNFBeXIAvGp2AU7maecDnndOnnY+4Hnn5GnnA604J49tAxBCCHF+nnwFIIQQ4jw8LgCUUlOUUnuVUllKqQfNrscZlFLZSqntSqktSim3XCNTKTVfKZWnlNrRaFu4UuoLpdR++9e2Xf7IiZo4n0eUUsfsn9MWpdTlZtbYHEqpnkqplUqpXUqpnUqp++zb3fkzauqc3PJzUkr5KaU2KaW22s/nUfv23kqpjfbfee8qpXwc3qcn3QJSSlmBfcBEIAfYDMzQWu8ytbBWUkplA6laa7ftv6yUGgOUAm9qrQfYt/0NKNRaP2kP6zCt9W/NrNNRTZzPI0Cp1vppM2trCaVUV6Cr1vpbpVQwkAlMA27GfT+jps7pZ7jh56SUUkCg1rpUKeUNfAPcB/wa+K/WepFS6hVgq9b6ZUf26WlXACOALK31Qa11NbAIuNrkmgSgtV4NFJ61+WrgDfvzNzD+c7qFJs7HbWmtc7XW39qflwC7ge6492fU1Dm5JW0otX/rbX9oYDywxL69WZ+RpwVAd+Boo+9zcOMPvBENLFdKZSqlbje7GCeK1lrn2p8fB6LNLMZJ7lZKbbPfInKb2yWNKaV6AUOAjXjIZ3TWOYGbfk5KKatSaguQB3wBHACKtNY19pc063eepwWAp7pYaz0UY6HQu+y3HzyKNu5Fuvv9yJeBOCAFyAX+bm45zaeUCgLeB+7XWhc3/pm7fkbnOCe3/Zy01rVa6xSgB8Ydj4TW7M/TAuAY0LPR9z3s29ya1vqY/Wse8AHGB+8JTtjv09bfr80zuZ5W0VqfsP8HrQNew80+J/t95feBt7XW/7VvduvP6Fzn5O6fE4DWughYCaQBoUopL/uPmvU7z9MCYDPQ194q7gNcD3xsck2topQKtDdgoZQKBCYBO87/LrfxMTDL/nwW8JGJtbRa/S9Ku5/gRp+TvYFxHrBba/2PRj9y28+oqXNy189JKRWplAq1P/fH6OyyGyMIpttf1qzPyKN6AQHYu3Q9A1iB+Vrrx00uqVWUUn0w/uoH8AL+447npJR6BxiHMXPhCWAu8CHwHhCDMevrz7TWbtGw2sT5jMO4raCBbOCORvfPXZpS6mJgDbAdqLNv/h3GPXN3/YyaOqcZuOHnpJQahNHIa8X44/09rfWf7L8jFgHhwHfAz7XWVQ7t09MCQAghhGM87RaQEEIIB0kACCFEByUBIIQQHZQEgBBCdFASAEII0UFJAAghRAclASCEEB2UBIAQQnRQ/x8HsUpwsXdeVAAAAABJRU5ErkJggg==\n",
+ "text/plain": [
+ "<Figure size 432x288 with 1 Axes>"
+ ]
+ },
+ "metadata": {
+ "needs_background": "light"
+ },
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "M_greedy = greedy_basis.gramian(m.h1_0_semi_product)\n",
+ "rhs = greedy_basis.inner(V, m.h1_0_semi_product)\n",
+ "greedy_proj_errors = []\n",
+ "for N in range(len(greedy_basis)):\n",
+ " try:\n",
+ " v = np.linalg.solve(M_greedy[:N, :N], rhs[:N, :])\n",
+ " except np.linalg.LinAlgError:\n",
+ " v = np.zeros((N, len(V)))\n",
+ " V_proj = greedy_basis[:N].lincomb(v.T)\n",
+ " greedy_proj_errors.append(np.max(m.h1_0_semi_norm(V - V_proj)))\n",
+ "\n",
+ "plt.figure()\n",
+ "plt.semilogy(errors, label='without greedy')\n",
+ "plt.semilogy(greedy_proj_errors, label='with greedy')\n",
+ "plt.legend()\n",
+ "plt.show()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Schließlich noch die Konditionen der Projektions-Matrizen:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 19,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAD9CAYAAAC7iRw+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xd8VtX9wPHPNxsyWAkrJCQgeyUQlggiylARkKICToq1WK3a1rZY2x/apa2jWmdRhpMhoiIgIigyxRBm2BsC2RCyyHzO74/7gAEJxOSZyff9euWV557n3nu+lyf55nDuueeIMQallFJ1h4+7A1BKKeVamviVUqqO0cSvlFJ1jCZ+pZSqYzTxK6VUHaOJXyml6hhN/EopVcdo4ldKqTrGz9EnFJE2wJNAA2PMOHvZQOBOe32djTFXO7pepZRSVVOlFr+IzBSRDBFJvqh8hIjsFZEDIjIVwBhzyBgzueJ+xpg1xpgpwGLgHUcFr5RS6qeraot/NvAq8O65AhHxBV4DhgIpQKKILDLG7LrMeSYCky/zPgDh4eEmJiamiqEppZQCSEpKyjLGRFxpvyolfmPMahGJuai4D3DAGHMIQETmAqOBSyZ+EYkGzhhj8q5UX0xMDJs2bapKaEoppexE5GhV9qvJzd1I4HiF7RQgUkSaiMibQLyIPFHh/cnArMpOJiIPiMgmEdmUmZlZg7CUUkpdjsNv7hpjsoEplyifdoXjpgPTARISEnTKUKWUcpKatPhPAFEVtlvZy5RSSnmwmrT4E4F2IhKLlfDHY928dYrS0lJSUlIoKipyVhXKSYKCgmjVqhX+/v7uDkUpRRUTv4jMAQYD4SKSAkwzxswQkYeBLwFfYKYxZqezAk1JSSE0NJSYmBhExFnVKAczxpCdnU1KSgqxsbHuDkcpRdVH9UyopHwpsNShEVWiqKhIk74XEhGaNGmC3rBXynN41ZQNmvS9k35uSnkWh4/qUUop9WPGGIrLbOQWlZJfVEZ+cRn5RWXknn9dSl5RGc0aBHF7QtSVT1gDmvgd6KabbuLDDz8E4MMPP+RXv/oVAKtWreL5559n8eLFNa5j1apVBAQEcPXV7pvu6MiRI4wcOZLk5OQr76xULVNWbuNUQQnZBSXkFJZy5mwpZ86WcOZsaYXtH75yCkvPJ/sy25VHql9zVbgmfm+ydKl1u+PIkSO8/vrr5xO/I61atYqQkJBqJ35jDMYYfHy8qpdPKacrLCnj2KlCMvOKycwrJiv/3PeSC7ZPFZZgKsnfvj5Cg3r+578a1Q8gNjyYsCB/QoP8CAnyIzTQj9Agf0IC7dtBfoQGWu8HB/oR4Of8301N/FX03HPPERgYyCOPPMJvfvMbtm3bxtdff83XX3/NjBkz+OCDD85PNTF16lQOHjxIXFwcQ4cO5eabbyY/P59x48aRnJxMr169eP/99xERVq5cyeOPP05ZWRm9e/fmjTfeIDAw8Py5wsPD2bRpE48//jizZ8/mzTffxNfXl/fff59XXnmFgQMHno8xMzOTiRMncvLkSfr3789XX31FUlIS+fn5DB8+nL59+5KUlMTSpUvZu3cv06ZNo7i4mLZt2zJr1ixCQkJISkrit7/9Lfn5+YSHhzN79mxatGhBUlISP//5zwEYNmzY+ToHDRrEf//7X+Li4gC45ppreO211+jRo4drPyClqig7v5gDGfkczCzgQEY+BzLzOZiRz4mcsz/aN9DPh4jQQCJCA4lqXJ/46EbWdkgATUICaVjPn7B6/jSsbyX6kEA/r7in5ZWJ/+nPd7LrZK5Dz9m5ZRjTbulS6fsDBw7khRde4JFHHmHTpk0UFxdTWlrKmjVrGDRo0AX7PvvssyQnJ7N161bAaqVv2bKFnTt30rJlSwYMGMC6detISEjgvvvuY+XKlbRv35577rmHN954g8cee+ySMcTExDBlyhRCQkJ4/PHHf/T+008/zZAhQ3jiiSdYtmwZM2bMOP/e/v37eeedd+jXrx9ZWVn8/e9/Z8WKFQQHB/Ovf/2LF198kSeeeIJf//rXfPbZZ0RERDBv3jyefPJJZs6cyaRJk3j11VcZNGgQv//978+fd/LkycyePZuXXnqJffv2UVRUpElfuURpuY3CknLOlpRTWFLG2dJzr8sveJ1XVMrhrAJ7ss/ndGHp+XPU8/elbdNgEmIaMT4iitiIYJqGBhERGkh4SIDXJPKfyisTvzv06tWLpKQkcnNzCQwMpGfPnmzatIk1a9bw3//+94rH9+nTh1atWgEQFxfHkSNHCA0NJTY2lvbt2wNw77338tprr1Wa+K9k7dq1fPLJJwCMGDGCRo0anX+vdevW9OvXD4DvvvuOXbt2MWDAAABKSkro378/e/fuJTk5maFDhwJQXl5OixYtyMnJIScn5/wfuLvvvpsvvvgCgNtuu42//e1vPPfcc8ycOZP77ruvWrErVVWnCkoY/tJqMvOKq3xMk+AA2kaEMKJrC65qGkLbiGCuahpCywb18PGpfYn9Srwy8V+uZe4s/v7+xMbGMnv2bK6++mq6d+/ON998w4EDB+jUqdMVjw8MDDz/2tfXl7Kyssvu7+fnh81mA3DI08rBwcHnXxtjGDp0KHPmzLlgnx07dtClSxc2bNhwQXlOTk6l561fvz5Dhw7ls88+Y/78+SQlJdU4VqUu56tdaWTmFfPLQW2ICA2kfoAf9QN8CfL3pX6A9VUvwJd6/r7UD/AjONCX0CB9arwivcP3EwwcOJDnn3+eQYMGMXDgQN58803i4+N/9F/B0NBQ8vKuOPs0HTp04MiRIxw4cACA9957j2uvvRawunXOJdGPP/64SuceMGAA8+fPB2D58uWcPn36kvv169ePdevWna+3oKCAffv20aFDBzIzM88n/tLSUnbu3EnDhg1p2LAha9euBeCDDz644Hz3338/jzzyCL17977gfxlKOcOy5DRaNarH1Bs7cv/ANkzsG82Y+EhGdG3OoPYRJMQ0pkvLBrSJCKF5gyBN+pegif8nGDhwIKmpqfTv359mzZoRFBR0wc3Vc5o0acKAAQPo2rXrBf3hFwsKCmLWrFncdtttdOvWDR8fH6ZMsSY2nTZtGo8++igJCQn4+vqeP+aWW27hk08+IS4ujjVr1lxwvmnTprF8+XK6du3KRx99RPPmzQkNDf1RvREREcyePZsJEybQvXt3+vfvz549ewgICGDBggX88Y9/pEePHsTFxbF+/XoAZs2axUMPPURcXBzmoiENvXr1IiwsjEmTJlX9H1OpasgrKmXdgWxGdGleK/veWfE0rHqWSocNOYhc/EvsCRISEszFC7Hs3r27Sl0qdVlxcTG+vr74+fmxYcMGHnzwwfM3mJ3p5MmTDB48mD179lQ6TFQ/P+UIi7ad5JE5W1gwpT8JMY3dHY5jHV4N74yCXvfBLS9V6xQikmSMSbjSfl7Zx68u7dixY9x+++3YbDYCAgJ46623nF7nu+++y5NPPsmLL76ozwYop/syOY3wkEB6RteyLsWzp+GTKdCkLQz/h9Or08Rfi7Rr144tW7a4tM577rmHe+65x6V1qrqpqLScb/ZmMCY+svaNxFnyOOSnw+SvICD4yvvXkDbRlFJeYe3+LApLyhnRpbm7Q3Gs7fMheQEMngqRPV1SpcMTv4i0EZEZIrKgQtlgEVkjIm+KyGBH16mUqv2W7UwjNMiPfm2auDsUx8k5Bkt+B1H94JrfuqzaKiV+EZkpIhkiknxR+QgR2SsiB0RkKoAx5pAxZvJFpzBAPhCEtSi7UkpVWVm5jRW707mhUzOXzGXjErZyq1/fGBj7P/DxvfIxDlLVf8HZwIiKBSLiC7wG3Ah0BiaISOdKjl9jjLkR+CPwdPVCVUrVVd8fPkVOYSnDa1M3z/r/wtF1cNO/oVGMS6uuUuI3xqwGTl1U3Ac4YG/hlwBzgdGVHG+zvzwNBF5qn9rgpptuOj+9weuvv36+fNWqVYwcOdKNkV1o9uzZPPzww+4OQ6kqW7YzjSB/H65tH+HuUBzj5Fb4+h/QeTT0uOQCh05Vk/8zRQLHK2ynAJEi0kRE3gTiReQJABEZKyL/A94DXr3UyUTkARHZJCKbvHWZvqVLl9KwYcMfJX5HutJUD0rVNjab4cudaVzbPoJ6Aa7rDnGakkJY+AsIDoeRL4EbHkRzeGeZMSbbGDPFGNPWGPOMvWyhMeaXxpg7jDGrKjluujEmwRiTEBHheX/Vn3vuufOTsf3mN79hyJAhAHz99dfceeedgDXNQlZW1gXTMp97cvfctMwdO3bkzjvv/NHTrwCJiYl07979/HFdu3YFrBb6qFGjGDJkCNdff/35eHr37k337t2ZNm3a+XO8//779OnTh7i4OH75y19SXl4OWE/etm/fnj59+rBu3ToA8vLyiI2NpbTUmq0wNzf3gm2lPMG2lBzSc4sZ0bWWdPN89X+QtQ/GvAH13fMQWk3G8Z8AKi4T08pe5nxfTIW0HY49Z/NucOOzlb7tjGmZr7nmmguOmzRpEm+99Rb9+/dn6tSpF7y3efNmtm/fTuPGjVm+fDn79+/n+++/xxjDqFGjWL169fmplNetW4e/vz+/+tWv+OCDDxg6dCjTpk0jKSmJBg0acN111xEfH09oaCiDBw9myZIljBkzhrlz5zJ27Fj8/XVuE+U5lu1Mw89HGNKhmbtDqbl9yyHxLej3ELS9zm1h1KTFnwi0E5FYEQkAxgOLHBOW57l4Wub+/fufn5b5UvP1XOzctMw+Pj7np2WuKCcnh7y8PPr37w/AxIkTL3h/6NChNG5stQ6WL1/O8uXLiY+Pp2fPnuzZs4f9+/ezcuVKkpKS6N27N3FxcaxcuZJDhw6xceNGBg8eTEREBAEBAdxxxx3nz3v//fcza9YswPpfgc63ozyJMYYvk9Po37YJDep7eYOkIAs+ewiadobr/8+toVSpxS8ic4DBQLiIpADTjDEzRORh4EvAF5hpjNnptEgrukzL3FlcPS3zxS6eVvmJJ57gl7/85QX7vPLKK9x7770888wzF5R/+umnlZ53wIABHDlyhFWrVlFeXn6+e0kpT7AvPZ8j2YX8YlAbd4dSM8bAokegKAfu/gT8g9waTlVH9UwwxrQwxvgbY1oZY2bYy5caY9rb+/OdP8GEmzl6WuaKGjZsSGhoKBs3bgRg7ty5le47fPhwZs6cSX5+PgAnTpwgIyOD66+/ngULFpCRkQHAqVOnOHr0KH379uXbb78lOzub0tJSPvroowvOd8899zBx4kRt7SuPsyw5DREY2tnLu3k2vwN7l8D106C5+xtXteRJCNdw9LTMF5sxYwa/+MUviIuLo6CggAYNGlxyv2HDhjFx4kT69+9Pt27dGDduHHl5eXTu3Jm///3vDBs2jO7duzN06FBSU1Np0aIFTz31FP3792fAgAE/+h/KnXfeyenTp5kwwfXDypS6nGU70+gV3Yimoe5tIVeLzQb5GXB4DSx7AmKvhX6/cndUgE7L7FHy8/MJCQkBrBvEqampvPzyy06vd8GCBXz22We89957TqujLnx+yrGOZRcy6Llv+PPNnbh/oId19dhskL4DclMhPw3y0iEv1ZpoLS/V2s5PB2ONqiOoITy4HhpEOjUsnZbZCy1ZsoRnnnmGsrIyWrduzezZs51e569//Wu++OILli5d6vS6lPopvtyZBuB5T+uWFsH8u2H/8gvL64dDaHPrq2mXH16HNofIXhDW0j3xXoImfg9yxx13XDDixhVeeeUVl9anVFUt25lG5xZhRDWu7+5QflB6FuZOhINfww1PQcwgCG0GwU3BL8Dd0VWZVyV+Y0ztXG6tlvPE7kTl2TJyi9h87DS/uaG9u0P5QUkhzBlvrZQ16lXoebe7I6o2r7m5GxQURHZ2tiYRL2OMITs7m6AgL7w5p9xm+a50jMFzntYtKYAPb7eS/pjXvTrpgxe1+Fu1akVKSgreOo9PXRYUFESrVq3cHYbyIl/uTCM2PJh2TUPcHQoU51tJ/9gGuPV/0MO13bHO4DWJ/9wDVEqp2u1MYSkbDmZz/8A27u/aLc6D98dBSiKMfQu6jXNvPA7iNYlfKVU3rNyTTpnNuL+bpygX3v8ZnEiCcTOgy63ujceBNPErpTzKsuQ0mocF0T3y0g8wusTZHCvpp26F22ZD51Hui8UJvObmrlKq9issKWP1/kyGd2mGj4+bunnOnob3xkDqNrj93VqX9EFb/EopD7J6XyZFpTaGu6ubp/CUlfQzdsMd70OHEVc+xgtpi18p5TGWJafRqL4/fWLcsEDJqcPw7ijI2AN3fFBrkz5oi18p5SFKymys3JPBiC7N8fN1YZs0Yw+sfRF2LADfAJjwIVx1g+vqdwOHJ34RaQM8CTQwxoyzl3UCHgXCgZXGmDccXa9SyrttOJRNXlGZ60bznNwCa16A3YvBvx70exD6PwxhLVxTvxtV6c+qiMwUkQwRSb6ofISI7BWRAyIyFcAYc8gYM7nifsaY3caYKcDtwABHBa+Uqj2WJacRHODLgKvCnVvR0fXw3liYPth6EnfQ7+GxZBj+jzqR9KHqLf7ZwKvAu+cKRMQXeA0YCqQAiSKyyBiz61InEJFRwIOA8+b+VUp5naLSchYkpbB4+0kGd2xKkL+v4ysxBg6uhNUvwLH11kya10+D3vdDUJjj6/NwVUr8xpjVIhJzUXEf4IAx5hCAiMwFRgOXTPzGmEXAIhFZAnxY3YCVUrXD6YIS3vvuKO+sP0J2QQk9ohry2PXtHFuJrRz2LLG6dFK3Qlgk3PhviL8bAjxo1k8Xq0kffyRwvMJ2CtBXRJoA/wDiReQJY8wzIjIYGAsEApec+F1EHgAeAIiOjq5BWEopT3b8VCEz1h5mXuJxzpaWM6RjU345qA19Yhs7boqGkkLY+gF89zqcOgSN28CoV6D7eK+aPtlZHH5z1xiTDUy5qGwVsOoKx00HpoO1Apej41JKuVfyiTP8b/Uhlu5IxUdgdFwkDwxqQ/tmoY6rJC8dvp8Om2ZYD2JF9oJxs6DTKPDVQYzn1ORf4gQQVWG7lb1MKaUAa1ruNfuzmL76EGsPZBES6Mf918QyaUAszRs4cKrujN2w4VXYPh/KS6HjzdYIneh+4O6J3jxQTRJ/ItBORGKxEv54YKJDolJKuYwxhmOnCjmaXUhJmY3Schsl5TaKz70us/1QXmajpNxQXFZOUamN4tJyisrKKS61UWQvKyotp7jM+l5QXMbpwlKahgYy9caOTOwbTViQv6MCh0OrrIR/YAX41bP67vs/BE3aOqaOWqpKiV9E5gCDgXARSQGmGWNmiMjDwJeALzDTGLPTaZEqpRyioLiMbSk5bDmWw5Zjp9lyLIfsgpIqHx/g50Ogrw+B/r4E+fsQ6OdDkL+v/cuHBvX8CfL3IcjPl0B/H+KjGzE6riWBfg4arZOXDge+gu/egPRka9nD6/4MvSdDfTc88euFqjqqZ0Il5Uup5GatUsr9jDEczipgy7EcNh87zeZjOexNy8Vmv4vWJiKYwR2a0rN1Q9o1DSXI3wd/Xx8C/HwIuOi7v68P/r7i+jnyC7LhyBrr6/AayNprlUd0spZA7HYb+OsKbz+F3u1QqpY4W1LOvvQ89qblsSctj73puew6mcvpwlIAQgL9iItqyMPXXUV860bEtWpIo2APHOFy9jQcWfdDos+wdyT4B0Pr/hA3EWIHQsue2n9fTZr4lfIy5TbDkeyCHxJ8Wi570/I4eqqQc0tSB/n70L5ZKMM6Nyc+uiHx0Y24qmkIvu6a6vhK0pJhx0dw8GtI2wEYq88+ui90/QvEDoKW8eDroPsDdZwmfqU8VFm5jWOnCtmfkc+BjHz2p+edf11cZgPARyCmSTCdWoQxJj6Sjs1D6dA8jOjG9T03yZ+Tl2Yl+23zIH0H+PhBVD8YPNVK9JG9wC/Q3VHWSpr4lfIAh7MK2J2ay/70fPZn5HEgI59DmQWUlNvO7xPZsB5XNQ2hX5smdGgeSsfmobRrGkq9ACdMceAsJQXWk7Tb5sKhb8DYrAR/43PQdSwEO3meHgVo4lfKrfan5/GvZXtYsTvjfFlU43q0axrKte0jaNcslHZNQ2jbNISQQC/9dbXZrP76bXNh9yIoyYcG0XDNb6HHeAh38DQN6oq89CdJKe+WkVvEf1bsY17icYID/Hh8WHsGd2hKm4hg6gfUkl9LWzmsexkS34bcExAYZi1Y3mM8RF8NProOlLvUkp8wpbxDfnEZ01cf4q3Vhygtt3FP/xh+PeQqmoTUsr7s/Az4eLI17XHb62HY36DDTda898rtNPEr5QKl5TbmJR7npRX7ycov5ubuLfj9sA7EhAe7OzTHO7oePpoERTkw+nWIv9PdEamLaOJXyomMMSzflc6/lu3hUGYBfWIa89Y9vYiPbuTu0BzPGFj/Cqx4ChrFwF0fQ/Ou7o5KXYImfqWcwBjDpqOn+feyPSQeOU3biGDeuieBGzo1df2Tr65wNgc+ewj2LIbOo60nauvgAifeQhO/Ug5SUmbju0PZrNidzopd6Zw8U0REaCD/vLUbtye0cu0C4q6Uug3m3wNnUmDEs9B3ij5R6+E08StVAzmFJXyzN4MVuzL4dl8m+cVlBPn7MLBdBI/e0I6R3VsS7K3DMK/EGNj8Diz9gzX+ftIXENXH3VGpKqilP5FKOc+RrAJW7E7nq13pbDp6mnKbISI0kFt6tOCGTs0YcFW4c9aN9SQlhbDkt7BtDrQdAmPf0oevvIgmfqUuw2Yz7M/IJ+noaTYfO03S0dMczioAoGPzUH41uC03dGpGt8gG+Hj6FAlVYQzYyqD0LJQVQ5n9e8Xt4jz4+u/W4ieDn4BBvwefWv6HrpbRxK9UBWfOlrL1eA6b7Yl+67Ec8orLAGgcHEDP6Ibc278113dqRlTjWrJYd34mLJsK+5ZBaaE1jcKV1G8Cdy+0WvvK6zg88YtIG+BJoIExZlxlZUp5ivUHsli07SSbj51mf0Y+xliTn7VvFsqouJb0jG5Er9aNaN2kfu0akWMMJH8MX/zBasXHTYT64eAXZM1v71fhq+K2fz1rhat6tXBIah1R1RW4ZgIjgQxjTNcK5SOAl7FW4HrbGPOsMeYQMFlEFpzb71JlSnmCz7ae4DfzthIS6Ed8dCNGdrcSfY+oBoQ6aolAT5SXbvXR71lszWs/5nVo2sndUSkXqWqLfzbwKvDuuQIR8QVeA4YCKUCiiCwyxuxydJBKOcO5pN87pjGzJvWuPXPkXI4x1mRpy6Za/fZD/wr9HgLfOnDt6ryqLr24WkRiLiruAxywt+YRkbnAaEATv/J4dTLpnzkBix+D/cshqi+Mfk1nxqyjavJESSRwvMJ2ChApIk1E5E0gXkSeALhU2cVE5AER2SQimzIzM2sQllKXV+eSvjGw+V14vZ+1lOGIZ60x95r06yyH/8QbY7KBKVcqu8Rx04HpAAkJCcbRcSkFdTDp5xyDRY9Yi57EDIRR/4XGbdwdlXKzmvzUnwCiKmy3spcp5ZHqTNIvL4NjG6wbt1vet8pufgF6/VznwFdAzRJ/ItBORGKxEv54YKJDolLKwWp90i8tgsPfWitc7f0CCrOtoZftR1hz4TeMdneEyoNUdTjnHGAwEC4iKcA0Y8wMEXkY+BJrOOdMY8xOp0WqVDXV2qRfnGfdqN292Ppekm+tctV+OHQcCVfdAIEh7o5SeaCqjuqZUEn5UmCpQyNSyoFqXdIvOgO7P7e+Dn4D5cXWQ1ddfwadRkHsQPCrZat5KYfz8t8CpSpXa5K+zQZH11r99bs+g7Iia7Hy3pOh0y3W0EydK0f9BF76m6DU5S3Znur9Sf9MCmz90Er4OUetbpy4iRB3F0T21DnvVbV54W+DUpeXdqaIqR9vJy6qofcl/bJi2LPESvYHvwYMxA6CIX+2+u0DasnEcMqtvOg3QqkrM8bw50+TKbXZ+M8dcd6T9NN2wOb3YMd8OHsawlrBtX+wWviNYtwdnaplvOS3Qqmq+SI5jRW70/nTTR1p3STY3eFcWUEWfPkkbJ8LvgFWn338XRB7rfbbK6fRxK9qjZzCEv7vs510jQzj5wNi3R3O5dlssOU9+Or/oKQABj4O/R+C+o3dHZmqAzTxq1rjn0t3c7qwhHd+3tuzFzbP2A2fPwbHv4PWA2DkfyCig7ujUnWIJn5VK6w7kMX8TSk8OLgtXVo2cHc4l1ZSCKv/DetfsUbojH7d6sPX0TnKxTTxK693tqScJxbuIKZJfR693kNnnNz/FSz5nTUsM+4uax784CbujkrVUZr4ldd7acU+jp0qZM4v+hHk72E3RHNTrUVPdn0K4e3hviUQc427o1J1nCZ+5dWST5zhrTWHmNAniv5tPagFbbPBphmw8q/W2Pzr/gwDHtHpFJRH0MSvvFZpuY0/LNhOk5BApt7oQevFZh+Ezx6GY+uhzXXWlMhN2ro7KqXO08SvvNbbaw6zKzWXN+/qSYN6HrAwus0G3/8PVjxtjcnXm7fKQ2niV17pcFYBL63Yx4guzRnRtYW7w7mwld9uGNzyMoS1dHdUSl2SSxK/iHQGngKygZXGmAWuqFfVTsYYnli4nQA/H54e3cW9wWgrX3mhaj/lIiIzRSRDRJIvKh8hIntF5ICITLUX3wi8Yox5ELinBvEqxbzE43x36BR/uqkTzcKC3BdI9kGYfbM1aid2IDz0HcTfqUlfebyaPN44GxhRsUBEfIHXsBJ9Z2CCvbX/HjBeRJ4DPGjohfI2GblF/GPpbvq1acz43lFXPsAZbDb47g14YwCk77Ra+RPna9eO8hrV7uoxxqwWkZiLivsAB4wxhwBEZC4w2hjzDPCQ/Q/DwurWqdS0RTspLrPxzNjuiDta1ln7YdEj2pevvJqj+/gjgeMVtlOAvvY/EH8CgoHnLnWgiDwAPAAQHa0LQ6sLHc4qYEHScb5ITuMPIzoQG+6CmTdt5ZCxC45vhGMbrbl1co5BYAPty1dezSU3d40xR7An9cvsMx2YDpCQkGBcEJbyYMYY9mfk88WONL5ITmVPWh4AQzo25RcD2zin0uI8SNkEx7+3knzKJijOtd4LaWYtcdh3CnQZC2EeMJJIqWpydOI/AVTseG1lL1PqiowpSQYYAAAZj0lEQVQx7DyZy7LkNJYmp3IoswARSGjdiL+M7MyIrs2JbFjPsZXaymH9fyF5IaQng7EBAk07Q7dxENUPovtCw9baule1hqMTfyLQTkRisRL+eGCig+tQtURZuY384jIOZRXwZXIaXySncexUIT4C/do0YdLVMQzv0pymzhq5c/Y0fHw/HFgB0f2tOfGj+kKrBKjX0Dl1KuUBqp34RWQOMBgIF5EUYJoxZoaIPAx8CfgCM40xOx0SqfIaeUWlzEs8TmZ+MXlFZeQVlZF7tpS8otIftotKKSwpP3+Mv69wddtwHrquLUM7N6dxcIBzg0xLhnl3wpkT1nz4vSZpi17VGTUZ1TOhkvKlwNJqR6S8Wn5xGffNSiTp6GkCfH0IDfIjrJ4/oUF+hAb50SwsyP7an7Agq7xpWCAD20W4btqF5I+tp2wDw2DSUojq45p6lfIQOmWDcpjCkjJ+PiuRrcdzeP3OntzUzcNugJaXwcqnrIVQovrC7e9CaHN3R6WUy2niVw5xtqScn89OZNPRU7w8Pt7zkn5BNiyYBIe/hd73w/BnwM/J3UlKeShN/KrGikrLuf/dRL4/fIr/3BHHLT087IGmk1th3t2Qnw6jX4P4u9wdkVJupYlf1UhRaTm/eHcT6w9m8/y4HoyOi3R3SBfaNhc+fxTqh8PPl0FkT3dHpJTbaeJX1VZcVs6U95NYsz+Lf4/rzs96tXJ3SD8oL4Xlf4aNb0LMQBg3C0Ii3B2VUh5BE7+qlpIyG796fzOr9mbyzNhu3J7gpgnTLlZ0BrbNg8S3IGsf9HvIWtjcV3/UlTpHfxvUT1ZabuOhDzezck8Gfx/TlQl9PGBupfSd8P1bsH0+lBZAZC+4433odIu7I1PK42jiVz9JabmNR+Zs4atd6Tw9qgt39WvtvmDKSmD3Ikh8G45tAL8g6DoOek/WvnylLkMTv6qysnIbj83byhfJafxlZGfuvTrGPYGcSYGk2ZD0DhRkQKNYGPZ3iLsT6jd2T0xKeRFN/KpKzpwt5S+fJrNkeypP3tSJydfEujaA8jI4+DVsfgf2LgVjoP0Ia0x+2yHgU5M1hZSqWzTxq8s6ml3ArHVHmL/pOIUl5fxxREd+MchJ0yJfzBg4uQW2z7OmWSjIhPpNYMCj1tw6jdzYzaSUF9PEr37EGMP3h0/x9trDrNidjp+PcEuPlky+JpYuLRs4P4DTR2HHfGt0TvZ+axHz9iOgx3i4aqg+catUDWniV+eVlNlYsuMkM9YeJvlELo3q+/PwdVdxd7/Wzpsa+Zyzp2Hnp9aonGPrrbLWA+Dqh6HzGJ0mWSkH0sSvyCks4YONx3h3wxHSc4tpGxHMP2/txq3xkdQL8HVu5UfWWg9Z7fsSyksgvD0M+Qt0u027cpRyEpckfhEZCNxpr6+zMeZqV9SrLs1mMxzKKmDb8Rw2Hs7m822pnC0tZ2C7cJ79WXeubReBj48T56Y3xrpRu/p5q3VfPxwSJkOPO6BFnM6Lr5ST1WQhlpnASCDDGNO1QvkI4GWshVjeNsY8a4xZA6wRkTFYq3QpF8rIK2Lb8TNsO57D1uM5bEvJIa+oDIDgAF9Gdm/B5IGxdGwe5txAjIF9y2D1c3AiCcIi4cbnoOfd4O/gJRWVUpWqSYt/NvAq8O65AhHxBV4DhgIpQKKILDLG7LLvMhGYXIM6VRUknzjD+oNZVpI/foYTOWcB8PUROjYP5ZYeLYmLakhcVEPaRoTg68zWPYDNZj1otfp5SN9hrV97y8vQYwL4BTq3bqXUj9RkBa7VIhJzUXEf4IAx5hCAiMwFRgO7RCQaOGOMyatunapyxhjWHsji9W8OsuFQNgBRjevRs3UjJg2IIS6qIV1aNnB+n31F5WWwc6GV8LP2QpN2MOZNaxFzXxettqWU+hFH9/FHAscrbKcAfe2vJwOzKjtQRB4AHgCIjvaAuV+8hM1m+HJnGq+vOsiOE2doFhbIn2/uxJj4SMJD3NSaLiuxxt6vfRFOHYKmnWHcTGt0jo8L//AopS7JZaN6jDHTrvD+dGA6QEJCgnFJUF6spMzGp1tP8Oa3BzmUWUBMk/o8O7Ybt/aMJNDPTcm1tAi2vAfrXoYzx60btXd8AB1u0idrlfIgjk78J4CK8/O2spcpByksKWPu98d5e80hTp4ponOLMF6dGM+NXVs4v6++MiUFsGkWrP+vtcpVVF8Y+R+46gYdoaOUB3J04k8E2olILFbCH491Q1fVQFm5jbTcIj7ZfIJZ649wqqCEPrGN+efYblzbPgJxV3ItOgPfT4cNr8PZUxB7LfxsBsRcowlfKQ9Wk+Gcc4DBQLiIpADTjDEzRORh4Eus4ZwzjTE7HRJpLXW2pJy03CJSz5wlPbeI1DNFpJ+xf7dvZ+UXY7N3fl3fsSkPDm5LQowbZ6EsyIaNb8DG6VB8BtoNh0GPQ1Qf98WklKqymozqmVBJ+VJgabUjqkOe/WIPb3578EfloUF+NA8LonmDINo3C6VFgyCaNQgioXVjOjQPdUOkdnnpsOEVSJxpLXbSaZSV8Fv0cF9MSqmfTKdscJOi0nI+2HiUvrGNuaN3FM3DrOTePCyI4EAP/FiOJ8K7o6CsyFrsZOBvoWknd0ellKoGD8wwdcOK3enkFZXx6yHtuKZduLvDuby8NJh3F4Q0hbsWQpO27o5IKVUDmvjdZOHmEzQPC6J/2ybuDuXyyoph3t1QnAt3a9JXqjbQwdVukJlXzLf7MhkTH+m+IZhVtfT3kPI9jHkDmnVxdzRKKQfQxO8Gi7adpNxmGNsz0t2hXN6mmdZShwN/B13GuDsapZSDaOJ3g4WbU+gW2YD2zdw4QudKjm6ApX+AdsPguifdHY1SyoE08bvY3rQ8dp7M9ezW/pkTMP8eaBgNY9/S+XWUqmX05q6LLdyScn4NW49UWmSN4Ck9C/ct1iUPlaqFNPG7ULnN8OmWEwzuEOG+mTMvxxhY/Bs4uRnGfwgRHdwdkVLKCbSrx4XWH8wiPbeYsT1buTuUS/t+Omz7EK6dCh1vdnc0Sikn0cTvQgs3nyAsyI8hHZu6O5QfO7wGlj1hTaF87R/dHY1Syok08btIfnEZy5LTGNmjJUH+HnazNOcYfHSv9XDWrf/TufOVquX0N9xFliWncba0nJ952miekkKYeyeUl8L4ORDk5AXXlVJup4nfRRZuTqF1k/r0jG7k7lB+kJcOn06BtB3ws7ch/Cp3R6SUcgGXjOoRkcHA34CdwFxjzCpX1OspTuScZcOhbB67vr37Fk05x2aDw6usFbP2LgVbGdzwNLQf7t64lFIuU5OFWGYCI4EMY0zXCuUjgJexFmJ52xjzLGCAfCAIawH2OuXTLScwBm6Nd2M3T146bH0fkt6BnKNQrzH0nQK97oPwdu6LSynlcjVp8c8GXgXePVcgIr7Aa8BQrASfKCKLgDXGmG9FpBnwInBnDer1KsYYFm5OoU9MY6Kb1Hdt5TYbHPoakmbD3i+s1n3MQLj+/6DTLeDngc8SKKWcriYrcK0WkZiLivsAB4wxhwBEZC4w2hizy/7+aaBOZZvtKWc4mFnALwa2cV2leemw5T1rgrWcY1C/CfR7EHreq617pZTD+/gjgeMVtlOAviIyFhgONMT6X8KPiMgDwAMA0dHRDg7LfRZuTiHAz4eburdwfmX5GbD2P5A4A8qL7a37adq6V0pdwCU3d40xC4GFV9hnOjAdICEhwbgiLmcrKbPx+fZUhnZuRliQv/MqKjwF6162nrwtK4IeE+Gax7R1r5S6JEcn/hNAVIXtVvayOunbfZmcKihx3tj9ojOw4XXY8BqU5EO3cdZ0CzosUyl1GY5O/IlAOxGJxUr444GJDq7DayzcnEJ4SAAD20U49sQlBbDxf1YrvyjH6soZ/Cdo1tmx9SilaqWaDOecAwwGwkUkBZhmjJkhIg8DX2IN55xpjNnpkEi9TE5hCSt3Z3BXv9b4+zroObnSImtVrLUvQkGmfZGUP0HLeMecXylVJ9RkVM+ESsqXAkurHVEtsXh7KiXlNscsuGKzWWPwv3kG8k5C7LUw5M8Q1afm51ZK1Tk6H7+TLNycQodmoXRpWcO5b7IPwuePwpE10KoPjP0fxA5yTJBKqTpJE78THM4qYPOxHJ64sWP1p2goL4PvXoNv/gm+AXDLy9Y4fHdP+aCU8nqa+J3gk80p+AiMqe4UDWk74LOHIXWrNT/+zS9AmIcu1aiU8jqa+B3MZjMs3HKCAVeF0yws6KcdXFYMq5+zHsIKagjjZkGXW7WVr5RyKE38DpZ45BQpp8/y+LCfuF7tsY2w6NeQtRe6j4cRz0D9xs4JUilVp2nid6CychvPL99LaJAfw7o0q9pBxfnw9d+scfkNWsGdH0O7G5wbqFKqTtPE70CvfXOQxCOneemOOOoHVOGf9sg6+GQKnDkGfR6wZs0MDHV+oEqpOk0Tv4NsOnKKl1fu49b4yKrd1N39OSz4OTSIgknLoHV/5weplFJo4neIM2dLeXTuVlo1qs9fR3e58gGb34PPH4HIXjBxvvblK6VcStfcrSFjDE9+soO03CJeHh9H6JVm4Vz/Cix6GNoMhns+06SvlHI5bfHX0IKkFBZvT+X3wzsQf7mF1I2xbuKueQE6j4Gx03WOfKWUW2jir4HDWQVMW7STfm0aM+XatpXvaCuHJb+DpFnWGrc3vwg+vi6LUymlKtLEX00lZTYenbsFf18f/nNHHL4+lTxkVVYCnzwAOz+Ba35rjdzRB7KUUm6kib+aXvhqL9tTzvDmXb1o0aDepXcqKYB5d8PBlTD0bzDgEdcGqZRSl6CJvxrW7s/if98eYmLfaEZ0bX7pnc6ehg9uhxObYNQr0PMe1waplFKVcEniF5FOwKNAOLDSGPOGK+p1hlMFJfx2/lbaRgTzl5srWfEqLw3eGwvZ++G2d6DzKNcGqZRSl1Ht4ZwiMlNEMkQk+aLyESKyV0QOiMhUAGPMbmPMFOB2YEDNQnYfYwx/WLCNnMJSXpnQk3oBl7hBe/oozBwOp49YY/Q16SulPExNxvHPBkZULBARX+A14EagMzBBRDrb3xsFLMGLV+d6/7ujrNidwdQbO9L5UgusGAOfPgiFp+HeRdD2OtcHqZRSV1DtxG+MWQ2cuqi4D3DAGHPIGFMCzAVG2/dfZIy5EbizunW60960PP6+ZDeDO0QwaUDMpXfa8REcXQfD/gqtElwan1JKVZWj+/gjgeMVtlOAviIyGBgLBFJJi19EHgAeAIiOjnZwWDWTmVfMI3O2EBrkz3Pjelx6Va2iXFj+Z2jZE+L1Rq5SynO55OauMWYVsOoK+0wHpgMkJCQY50dVaRyknD7L94dPkXjkFN8fOcWhzAIAZk/qTURoJU/brnoW8jNgwhzw0ZkwlFKey9GJ/wQQVWG7lb3MY9lshgOZ+Ww8fIpEe7JPPVMEQFiQH71jGnN7QhSD2kVcul8fIH0XbHwTet1rTbymlFIezNGJPxFoJyKxWAl/PDDRwXU4xJmzpTz5yQ7WHsgip7AUgKahgfSJbXz+q33TUHwqeyL3HGNg6eMQFAbXT3NB5EopVTPVTvwiMgcYDISLSAowzRgzQ0QeBr4EfIGZxpidDonUwT7adJzF21MZ16sVfe2JPrpx/Uv331/OjgXWDd2RL+lMm0opr1DtxG+MmVBJ+VK8YMjm59tT6RoZxvO39aj+Sc7f0I3XJ3OVUl6jTt6FPH6qkG3HcxjZvWXNTvTtvyA/HW5+QWfbVEp5jTqZ+BdvTwXg5m4tqn+S9F3w3RtWS19v6CqlvEidTPyfbztJfHRDohrXr94JjIGlv9cbukopr1TnEv/BzHx2pebWrJsn+WM4utaaWz+4ieOCU0opF6hziX/xtlREatDNU5QLXz4JLeKg572ODU4ppVygTs3Hb4zh8+0n6R3TmOYNgqp3knM3dMd/qDd0lVJeqU61+Pem53EgI59belSzmydjt/2G7t3QSm/oKqW8U51K/Iu3peIjcGNlq2ZdzrkbuoGhcP1TDo9NKaVcpc509RhjWLz9JFe3DSc8pJKJ1i4n+WM4sgZuflFv6CqlvFqdafEnn8jlSHYht/Soxk3dc0/otugBve5zeGxKKeVKdSbxL95+Ej8fYXiXn9jNk7IJ3hpiraN7kz6hq5TyfnUi8VvdPKkMah9Bw/oBVTuorBhWPA0zhkJpIdy9EKJ6OzdQpZRygTrRx7/5WA4ncs7yu2Htq3ZA6nb4ZApk7IS4u2DEPyGogXODVEopF6kTif/zbScJ8PNhaOdml9+xvBTW/scaq1+/CUyYBx1GXP4YpZTyMi5J/CLSBngSaGCMGeeKOs8ptxmW7kjlug4RhAb5V75jxm6rlZ+6FbqOg5ue0/n1lVK1UrX7+EVkpohkiEjyReUjRGSviBwQkakAxphDxpjJNQ22Or4/fIqMvOLKH9qylcPal+B/g+DMcbj9XRg3Q5O+UqrWqsnN3dnABf0gIuILvAbcCHQGJohI5xrUUWOLt5+knr8vQzo2/fGbWQdg5ghYMQ3aDYNfbYTOo10fpFJKuVBNVuBaLSIxFxX3AQ4YYw4BiMhcYDSwq7r11ERZuY0vktO4vlNT6gdcdKkHv4E5E8AvAMa+Bd1ug5+67KJSSnkhRw/njASOV9hOASJFpImIvAnEi8gTlzpQRB4QkU0isikzM9Mhwaw/mM2pgpIfd/OkboN5d0HjNlYrv/vtmvSVUnWGS27uGmOygSlX2Gc6MB0gISHBOKLexdtPEhrox7XtI34oPHUY3h8H9RrBXQsgrAarcCmllBdydOI/AURV2G5lL3O5kjIby5LTGNqlGUH+9qdt8zPh/bFgK4W7FkNYDdfcVUopL+Torp5EoJ2IxIpIADAeWOTgOqpkzf5McovKuOXcSlvF+fDh7ZCbChPnQ0QHd4SllFJuV5PhnHOADUAHEUkRkcnGmDLgYeBLYDcw3xiz0zGh/jSfbztJw/r+DLgq3How66N7rTH642ZCVB93hKSUUh6hJqN6JlRSvhRYWu2IHKCotJyvdqVzS4+WBPgKfPprOLACbnkZOt7kztCUUsrtauUkbav2ZlBQUm4tqL7yadg2Bwb/SadUVkopamni/3xbKuEhAfTPWmDNvdNrElz7B3eHpZRSHqHWJf6C4jJW7knn8Va78f1yKnQcCTe/oOP0lVLKrtbNzrlidzpx5cncfuzfENUXfva2Lp6ilFIV1LoW/5bEtbwd8CLSOBYmzAH/eu4OSSmlPEqtSvy5aYd4MOWP2PyDkbs+1hk2lVLqEmpV4s+3+ZMV3I7UkR9Aw6grH6CUUnVQrerjb9kyipZ/+MrdYSillEerVS1+pZRSV6aJXyml6hhN/EopVcdo4ldKqTpGE79SStUxmviVUqqO0cSvlFJ1jCZ+pZSqY8QYh6xr7lAikgkcvag4HMhyQzjOptflfWrrtdXW64Lae20XX1drY0zElQ7yyMR/KSKyyRiT4O44HE2vy/vU1murrdcFtffaqntd2tWjlFJ1jCZ+pZSqY7wp8U93dwBOotflfWrrtdXW64Lae23Vui6v6eNXSinlGN7U4ldKKeUAHp/4RWSEiOwVkQMiMtXd8TiSiBwRkR0islVENrk7nuoSkZkikiEiyRXKGovIVyKy3/69kTtjrK5Kru0pETlh/9y2ishN7oyxOkQkSkS+EZFdIrJTRB61l3v153aZ6/Lqz0xEgkTkexHZZr+up+3lsSKy0Z4f54lIQJXO58ldPSLiC+wDhgIpQCIwwRizy62BOYiIHAESjDFePb5YRAYB+cC7xpiu9rJ/A6eMMc/a/2A3Msb80Z1xVkcl1/YUkG+Med6dsdWEiLQAWhhjNotIKJAEjAHuw4s/t8tc1+148WcmIgIEG2PyRcQfWAs8CvwWWGiMmSsibwLbjDFvXOl8nt7i7wMcMMYcMsaUAHOB0W6OSV3EGLMaOHVR8WjgHfvrd7B++bxOJdfm9YwxqcaYzfbXecBuIBIv/9wuc11ezVjy7Zv+9i8DDAEW2Mur/Hl5euKPBI5X2E6hFnyIFRhguYgkicgD7g7GwZoZY1Ltr9OAZu4MxgkeFpHt9q4gr+oOuZiIxADxwEZq0ed20XWBl39mIuIrIluBDOAr4CCQY4wps+9S5fzo6Ym/trvGGNMTuBF4yN6tUOsYqz/Rc/sUf7o3gLZAHJAKvODecKpPREKAj4HHjDG5Fd/z5s/tEtfl9Z+ZMabcGBMHtMLqDelY3XN5euI/AURV2G5lL6sVjDEn7N8zgE+wPszaIt3e33qu3zXDzfE4jDEm3f5LaAPewks/N3tf8cfAB8aYhfZir//cLnVdteUzAzDG5ADfAP2BhiLiZ3+ryvnR0xN/ItDOfuc6ABgPLHJzTA4hIsH2m0+ISDAwDEi+/FFeZRFwr/31vcBnbozFoc4lRrtb8cLPzX6zcAaw2xjzYoW3vPpzq+y6vP0zE5EIEWlof10Pa8DLbqw/AOPsu1X58/LoUT0A9mFXLwG+wExjzD/cHJJDiEgbrFY+gB/wobdem4jMAQZjzRSYDkwDPgXmA9FYM63ebozxupuklVzbYKwuAwMcAX5ZoV/cK4jINcAaYAdgsxf/Cas/3Gs/t8tc1wS8+DMTke5YN299sRrs840xf7XnkblAY2ALcJcxpviK5/P0xK+UUsqxPL2rRymllINp4ldKqTpGE79SStUxmviVUqqO0cSvlFJ1jCZ+pZSqYzTxK6VUHaOJXyml6pj/B/LXniaSmSDpAAAAAElFTkSuQmCC\n",
+ "text/plain": [
+ "<Figure size 432x288 with 1 Axes>"
+ ]
+ },
+ "metadata": {
+ "needs_background": "light"
+ },
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "greedy_conds = []\n",
+ "for N in range(1, len(U)):\n",
+ " greedy_conds.append(np.linalg.cond(M_greedy[:N, :N]))\n",
+ "plt.figure()\n",
+ "plt.semilogy(range(1, len(U)), conds, label='without greedy')\n",
+ "plt.semilogy(range(1, len(U)), greedy_conds, label='with greedy')\n",
+ "plt.legend()\n",
+ "plt.show()"
+ ]
+ }
+ ],
+ "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.2"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}