{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# EART97012 \n",
" \n",
"# Geophysical Inversion \n",
" \n",
"## Lecture 2 "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Today's learning objectives \n",
" \n",
" \n",
"1. To introduce concepts underlying the computational tasks of optimisation and inversion\n",
"\n",
"\n",
"2. To review some more key linear algebra results, especially on the solvability of square problems, rank and range etc\n",
"\n",
"\n",
"3. To introduce some simple examples of over- and under-determined problems\n",
"\n",
" \n",
" \n",
"\n",
" \n",
"\n",
"Note that the linear algebra thory and associated simple code examples unpin huge parts of computer/computational science, data science, machine learning and AI, etc. So these are very valuable and universal concepts both within and outside of the relatively narrow focus of this module on inversion methods.\n",
" \n",
" \n",
" \n",
"I'll skip over some of the longer text sections to focus on explaining the linear algebra concepts using the whiteboard. Please read these through in your own time."
]
},
{
"cell_type": "markdown",
"metadata": {
"toc": true
},
"source": [
"Table of Contents \n",
""
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"scrolled": false
},
"outputs": [],
"source": [
"%matplotlib inline\n",
"%precision 6\n",
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"import scipy.linalg as sl\n",
"from pprint import pprint"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Introductory comments \n",
"\n",
"In this lecture we will consider some more important linear algebra theory, in particular what it tells us about solutions to *square* problems.\n",
"\n",
"We will introduce some simple examples of non-square systems - the case of over- and under-determined problems - with solution procedures coming in future weeks."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Examples (continued)\n",
"\n",
"At the end of the last lecture we ended with the curve-fitting example, where we had more equations/constraints that unknowns (an example of an *over-determined* problem) and introduced the concept of the least squares solution. \n",
"\n",
"Let's consider some more examples.\n",
"\n",
" \n",
"\n",
"## Inversion example (more complex)\n",
"\n",
"In this example from some current research we have a model for tides in the Bristol Channel and Severn Estuary, there are multiple numerical and physical \"parameters\" that go into the model, including incoming tidal boundary conditions, bathymetry etc, some of which are known to varying levels of uncertainty. \n",
"\n",
"In this example we consider bed roughness (or bottom friction) as the parameters we wish to invert for, given the data of time series of tidal heights at tide gauges indicated by the red dots in the following image which also shows the discretised domain and the computational mesh\n",
"\n",
" \n",
"\n",
"Using *a priori* information on the approximate distribution of sediment grain sizes on the seabed, we partition the domain into three zones: rock, gravel and sand, and assign one parameter value to each. In this case we therefore have more data than parameters.\n",
"\n",
"In this case an iterative approach is taken and the following figure shows how the inversion progresses\n",
"\n",
" \n",
"\n",
"$J$ is a misfit function between modelled tides and observations and we solve the inversion problem by seeking to minimse $J$ - we see it drop as the model fit to data improves with iteration. How the three parameter values vary with iteration is also shown.\n",
"\n",
"The quality of the prediction can then be compared at the green squares where data was not used in the inversion. We find that the error in prediction at these independent locations is indeed reduced within the same numerical model. We have also found that the prediction is improved when these physical values are used in a completely independent tidal model. This latter fact gives us some indication that our inversion is telling us something useful physically and not simply correcting for biases in the first model.\n",
"\n",
" \n",
"\n",
"For publications on this topic see \n",
"\n",
"- \n",
"\n",
"- "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Seismic inversion\n",
"\n",
"The more practical half of this module will focus on seismic inversion and similar applications, which you will already be somewhat familiar with from other modules."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Optimisation\n",
"\n",
"Optimisation, in a single real variable, seeks to find a value for $x$ such that a scalar function \n",
"\n",
"$$f(x)$$\n",
"\n",
"has a minimum value (or equivalently that $-f(x)$ has a maximum value). \n",
"\n",
"For example, $f(x)$ might be the energy of a particle at position $x$. The value of $x$ for which this is a minimum will provide a rest position for the particle. \n",
"\n",
"Or we could have a problem where we are seeking a design of something that minimises cost or maximises profit.\n",
"\n",
"There are many ways to solve such problems, and again there is no general solution. In practical real problems, there may be no solution, a single unique solution, several solutions or infinitely many solutions. \n",
"\n",
"### Optimisation as an inversion problem\n",
"\n",
"One way to address the optimisation problem is to try to solve \n",
"\n",
"$$f'(x)=0$$\n",
"\n",
"for $x$. With this approach note that we have converted an optimisation problem into an inversion problem (we are looking for the $x$ values that satisfy this equation - this will give us the [stationary points](https://en.wikipedia.org/wiki/Stationary_point) of the problem).\n",
"\n",
"\n",
"### Inversion as an optimisation problem\n",
"\n",
"In practice, it is often the case that in order to solve an inversion problem, it is instead converted into an optimisation problem. \n",
"\n",
"For example, we can define a scalar misfit function measuring in some sense the difference between model predictions and observations and seek to iteratively minimise this misfit. \n",
"\n",
"The least squares error is an example of this. We showed previously (*) that we could write down the general solution to this problem. But it's only feasible to use that approach for small cases. For larger problems we will need iterative solution methods such as in the tidal example above. \n",
"\n",
" \n",
"\n",
"(*) From last lecture an over-determined problem $V\\boldsymbol{a} = \\boldsymbol{y}$ has least-squares approximate solution $\\boldsymbol{a} = (V^TV)^{-1}V^T\\boldsymbol{y}$."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Optimisation - simple example\n",
"\n",
"Consider the problem: find $\\boldsymbol{x}\\equiv (x,y)$ which minimises the function\n",
"\n",
"$$ f(\\boldsymbol{x}) = 1+2x + 4y + x^2+2xy+3y^2$$\n",
"\n",
"\n",
"The following image shows a contour plot of the function, and the red star indicates the $(x,y)$ location of the minima\n",
"\n",
" \n",
"\n",
"A homework exercise asks you to compute the minima and to generate this image."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Optimisation (more complex)\n",
"\n",
"The following movie shows iterations from an optimisation problem where our task is to maximise the power generated by an array of 256 individual tidal turbines. At every iteration a shallow water solver computes the flow field and the power of the array. The design parameters are then the $(x,y)$ location of each of the turbines.\n",
"\n",
" \n",
"\n",
"For a publication on this see\n",
"\n",
"- \n",
"\n",
" \n",
" \n",
"\n",
"\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# More on forward and inverse problems \n",
"\n",
"## Abstract problem setting\n",
"\n",
"Building on the above motivation and various examples, let's think a bit more about the problem at hand.\n",
"\n",
"Our starting point is a model for a system that can be written as \n",
"\n",
"$$G(\\boldsymbol{m}) = \\boldsymbol{d}$$\n",
"\n",
"where $G$ is a function or operator that produces the output, or data, vector $\\boldsymbol{d}$ for given input parameters $\\boldsymbol{m}$. So $G$ can be thought of as a numerical model describing the underlying \"physics\", as well as things like the geometry and boundary conditions etc.\n",
"\n",
"Generation of the forward model often involves the discretisation of (partial differential equations) and actually solving the **forward model**: given $\\boldsymbol{m}$ compute $\\boldsymbol{d}$, may require linear or nonlinear solves so may be very costly, but in general is straightforward once a stable and robust $G$ has been constructed.\n",
"\n",
" \n",
"\n",
"The **inverse problem** asks the inverse model, for a given (e.g. set of observational data $\\boldsymbol{d}$) what can we conclude about the model parameters $\\boldsymbol{m}$, i.e. conceptually can we do thie following:\n",
"\n",
"$$\\boldsymbol{m} = G^{-1}(\\boldsymbol{d}) $$\n",
"\n",
"where in general we won't be able to explicitly form the inverse operator $G^{-1}$. Remember at the stage of thinking $G$ isn't necessarily just a matrix, it could be a complex computer code, and so writing $G^{-1}$ hides the complexities involved in actually performing this inversion in practice.\n",
"\n",
" \n",
"\n",
"Additionally, it will often be the case that the inverse equation has no exact solution at all (**non-existence**), or that there will be infinitely many solutions that all fit the data equally accurately (**non-uniqueness**), or that the solution for the model parameter $\\boldsymbol{m}$ varies significantly when only very small changes are made to the input data $\\boldsymbol{d}$ (**instability**). \n",
"\n",
"If you think back to the lecture I taught on PDEs - what we're saying here is that we will not in general, for perhaps a variety of reasons, have a well-posed problem to solve.\n",
"\n",
" \n",
"\n",
"We also need to contend with the fact that in the real world, observations always contain errors, as well as the fact that our representation of the physics $G$ will be imperfect, or our mesh resolution will be finite and some processes will be missing. We therefore need to instead think about the relation\n",
"\n",
"$$\\boldsymbol{d} = G(\\boldsymbol{m}) + \\boldsymbol{e}$$\n",
"\n",
"where $\\boldsymbol{e}$ represents the errors in the observations and/or model from all contributions which we may be able to estimate but do not know explicitly.\n",
"\n",
"We need to consider how such errors might impact upon the inversion process and its accuracy."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Real world problems\n",
"\n",
"Numerical inversion then is the theory and practice of solving equations that can be represented by $\\boldsymbol{d} = G(\\boldsymbol{m})$ or $\\boldsymbol{d} = G(\\boldsymbol{m}) + \\boldsymbol{e}$, and that produces useful solutions for $\\boldsymbol{m}$ efficiently, despite the problems of the non-existence, non-uniqueness, instability and data errors that will typically characterise these problems. \n",
"\n",
"Some key elements of real-world numerical inverse problems are: \n",
"\n",
"\n",
"- Inverse problems are typically characterised by a large number of parameters that may be conceptually infinite, although in practice we will only ever be able to use a finite number of model parameters in order to be able to solve the forward problem computationally. \n",
"\n",
"\n",
"- The number of observed data will always be finite. \n",
"\n",
"\n",
"- Observed data always contain measurement errors and noise.\n",
"\n",
"\n",
"- It will usually be impossible to estimate perfectly all the model parameters from data that are likely to be inaccurate, insufficient or inconsistent, but it is possible nonetheless to extract useful and sometimes powerful constraints about the values of these parameters and/or the relationships between them.\n",
"\n",
"\n",
"- We will almost always have additional prior information about the plausibility of model parameters and their inter-relationships. \n",
"\n",
"\n",
"- Most important problems will be non-linear.\n",
"\n",
"\n",
"- The computational cost of the forward problem is likely to be significant so that trial-and-error methods are likely to be impractical.\n",
"\n",
"\n",
"- Understanding and quantifying the uncertainty in the model parameters may be as important as estimating their values.\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Incomplete, inconsistent and inadequate data \n",
"\n",
"In the solution of any inverse problem, three important questions arise: \n",
"\n",
"\n",
"1. Does the solution exist? \n",
"\n",
"\n",
"2. Is the solution unique? \n",
"\n",
"\n",
"3. Is the solution stable? \n",
"\n",
"\n",
"From a physical point of view, we will often know that there must be a solution since we are likely to be studying real physical systems. But from the mathematical/computational point of view, given the errors in the observations as well as in $G$, there may be no model parameters at all that exactly fits the observed data. \n",
"\n",
"If there is a model parameter set that fits the data, say $G(\\boldsymbol{m}_1) = \\boldsymbol{d}$, then that solution will not be unique if there is also another model parameter set $\\boldsymbol{m}_2$ such that $G(\\boldsymbol{m}_2) = \\boldsymbol{d}$. In this case, it will be impossible to distinguish these two model parameter sets from the given data, though we may have good a priori reasons to prefer one over the other. \n",
"\n",
"The final question of stability is critical in inversion theory. Suppose that two different model parameter sets $\\boldsymbol{m}_1$ and $\\boldsymbol{m}_2$ generate two different datasets $\\boldsymbol{d}_1$ and $\\boldsymbol{d}_2$. Also assume that these two sets of parameters are very different, while the difference between the two datasets is small and lies within the noise level $\\boldsymbol{e}$. In this case, if we attempt to fit the observed data exactly, then the resulting parameter sets may vary dramatically but solely as a result of changes in the noise. \n",
"\n",
"In the early part of the twentieth century, it was widely considered that a mathematical problem was formulated correctly if all the three questions posed above had a positive answer. A mathematical problem was considered well-posed if its solution did exist, was unique, and was stable. In contrast, a problem was ill-posed, and was not considered physically or mathematically meaningful, if its solution did not exist, or was not unique, or was not a continuous function of the data (i.e. if for a small perturbation of the data there corresponded an arbitrarily large perturbation of the solution). \n",
"\n",
"It turns out that the majority of practical problems in many areas of mathematical physics and many other disciplines are ill-posed. Fortunately, such ill-posed problems can be both physically and mathematically meaningful, and they can be solved to provide useful, if necessarily imprecise, solutions. \n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Linear and nonlinear problems\n",
"\n",
"In linear problems, there is a linear relationship between the model parameters and the data. \n",
"\n",
"Linear problems by definition obey the relations\n",
"\n",
"$$\n",
"\\begin{align*}\n",
"G(\\boldsymbol{m}_1 + \\boldsymbol{m}_2) &= \\boldsymbol{d}_1 + \\boldsymbol{d}_2\\\\\n",
"G(\\alpha \\boldsymbol{m}) &= \\alpha G(\\boldsymbol{m})\n",
"\\end{align*}\n",
"$$\n",
"\n",
"or equivalently\n",
"\n",
"$$\n",
"\\begin{align*}\\boldsymbol{d}_{\\boldsymbol{m}_1 + \\boldsymbol{m}_2} &= \\boldsymbol{d}_{\\boldsymbol{m}_1} + \\boldsymbol{d}_{\\boldsymbol{m}_2}\\\\\n",
"\\boldsymbol{d}_{\\alpha\\boldsymbol{m}} &= \\alpha \\boldsymbol{d}_{\\boldsymbol{m}}\n",
"\\end{align*}\n",
"$$\n",
"\n",
"\n",
"Thus, in a linear problem, if the magnitude of the model parameters is changed then the magnitude of the observed data is changed by an equivalent factor (scalability), and if two models are summed then the new data produced will be simply the sum of the two original datasets (superposition). \n",
"\n",
" \n",
"\n",
"For all linear problems, our forward model can be written as \n",
"\n",
"$$G\\, \\boldsymbol{m} = \\boldsymbol{d}$$\n",
"\n",
"where now $G$ is a matrix. \n",
"\n",
"\n",
" \n",
"\n",
"\n",
"However, we have a crucial new aspect now - $G$ will NOT typically be a square matrix - that is, the number of equations and unknowns will typically not be the same. Even when $G$ is square, it will often also be singular so that the inverse problem will not often have an immediately straightforward solution. Remember that for a singular square matrix we may have zero or infinitely many solutions depending on the exact nature of the data/RHS vector - more on this next lecture.\n",
"\n",
" \n",
"\n",
"\n",
"In contrast, non-linear problems do not obey superposition and scalability, so that if the magnitude of the model parameters is varied, then the data will not simply be changed by an equivalent scaling in magnitude, and if two models are combined by addition then the resultant data will not be a simple addition of the two original datasets. \n",
"\n",
"Many important physical problems are non-linear, and in general the solution of non-linear problems is more difficult than the solution of linear problems. \n",
"\n",
"Fortunately it is often possible to solve non-linear problems by solving a sequence of linear problems that approximate the non-linear problem. "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Examples\n",
"\n",
"Two important examples of inversion problems (that are often cast as optimisation problems) are data assimilation and imaging.\n",
"\n",
"- [Data assimilation](https://en.wikipedia.org/wiki/Data_assimilation) (e.g. of the atmosphere in numerical weather prediction) seeks to combine a model with observational data for lots of different goals, e.g. to provide the best estimate for the initial condition of a weather forecast.\n",
"\n",
"\n",
"- [Imaging](https://en.wikipedia.org/wiki/Inverse_problem) (e.g. medical or geophysical) seeks to reveal hidden internal details (of bodies or the Earth) given observational data of travel times of some sort of wave (e.g. ultrasound scans) that is transmitted through the body being imaged. The underlying model $G$ is then often a wave equation solver.\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# An example - straight-ray tomography\n",
"\n",
"[\"Tomography is imaging by sections or sectioning, through the use of any kind of penetrating wave\"](https://en.wikipedia.org/wiki/Tomography).\n",
"\n",
"Most acoustic and seismic inverse problems are non-linear (since changing the model parameters - the structure of the body being imaged will change the ray path (e.g. due to refraction or reflections) the wave travels along, and hence change the arrival times (the data) in a nonlinear manner). \n",
"\n",
"However, if we assume ray that paths do not vary with the \"model\" (which we could justify if we assume that we don't change the $\\boldsymbol{m}$ values very much), and if the data are travel times and the model is composed of [slowness](https://en.wikipedia.org/wiki/Slowness_(seismology)) (reciprocal of velocity) values, then the problem is linear. \n",
"\n",
"Consider the following example\n",
"\n",
" \n",
"\n",
"\n",
"Suppose that a wall is assembled from an array of bricks as above, and that each brick is composed of a different type of clay. If the acoustic velocities of the different clays differ, then one might attempt to distinguish the different kinds of bricks by measuring the travel time of sound across the various rows and columns of bricks, in the wall. The data in this problem are $n = 8$ measurements of travel time (based on there being 4 rows and 4 columns) \n",
"\n",
"$$ \\boldsymbol{d} = [T_1 , T_2 , T_3 , \\ldots, T_8]^T$$\n",
"\n",
"We will assume that each brick is composed of a uniform material and that the travel time of sound across each brick is proportional to the width and height of the brick $h$. The proportionality factor is the brick’s slowness $s_i$, giving 16 model parameters \n",
"\n",
"$$ \\boldsymbol{s} = [s_1 , s_2 , s_3 , \\ldots, s_{16}]^T$$\n",
"\n",
"where the ordering of the model parameters is according to the numbering scheme in the schematic above. The travel time of acoustic waves through the rows and columns of a square array of bricks is measured with the acoustic source and receiver placed on the edges of the square. \n",
"\n",
"The set of linear simultaneous equations that we must now solve are: \n",
"\n",
"$$\n",
"\\begin{align*}\n",
"hs_1+hs_2 + hs_3 + hs_4 &= T_1\\\\\n",
"hs_5+hs_6 + hs_7 + hs_8 &= T_2\\\\\n",
"&\\vdots \\\\\n",
"hs_4+hs_8 + hs_{12} + hs_{16} &= T_8\n",
"\\end{align*}\n",
"$$\n",
"\n",
"where we order our firing of waves from left to right, starting at the top moving down, and then from top to bottom starting at the left and moving right.\n",
"\n",
"These algebraic equations can be arranged as a matrix equation $G\\, \\boldsymbol{m} = \\boldsymbol{d}$ of the form\n",
"\n",
"$$\n",
"\\begin{pmatrix}\n",
"h & h & h & h & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\\\\n",
"0 & 0 & 0 & 0 & h & h & h & h & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\\\\n",
"\\vdots & \\vdots & \\vdots & \\vdots & \\vdots & \\vdots & \\vdots & \\vdots & \\vdots & \\vdots & \\vdots & \\vdots & \\vdots & \\vdots & \\vdots & \\vdots \\\\\n",
"0 & 0 & 0 & h & 0 & 0 & 0 & h & 0 & 0 & 0 & h & 0 & 0 & 0 &h \\\\\n",
"\\end{pmatrix}\n",
"\\begin{pmatrix}\n",
"s_1\\\\\n",
"s_2\\\\\n",
"\\vdots \\\\\n",
"s_{16}\n",
"\\end{pmatrix}\n",
"=\n",
"\\begin{pmatrix}\n",
"T_1\\\\\n",
"T_2\\\\\n",
"\\vdots \\\\\n",
"T_8\n",
"\\end{pmatrix}.\n",
"$$\n",
"\n",
"This tomography problem is linear because the assumptions we have made were designed to make the travel times depend linearly upon the model parameters. In addition, slowness rather than velocity is used for the model parameters, and slowness has a linear relationship to travel time while velocity does not. \n",
"\n",
"X-ray tomography used in medical imaging, and radar tomography used in atmospheric physics, typically also assume ray paths that are model independent and employ model parameterisations that are related linearly to the data, and so are able to use linear inversion methods. In seismic tomography in the interior of the Earth, the independence of the ray path and the velocity model is generally a poor assumption - in reality ray paths change significantly when seismic models change, so that linear inversion is not normally adequate for seismic tomography. \n",
"\n",
"We will return to a version of this example in a later lecture as well as more on this topic in the practical parts of this module."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"How is this problem as formulated above different to the curve-fitting problem from earlier which could be written as\n",
"\n",
"$$\n",
"\\begin{pmatrix}\n",
"1 & x_0 & x_0^2 \\\\\n",
"1 & x_1 & x_1^2 \\\\\n",
"\\vdots & \\vdots & \\vdots \\\\\n",
"1 & x_5 & x_5^2\n",
"\\end{pmatrix}\n",
"\\begin{pmatrix}\n",
"a_0\\\\\n",
"a_1\\\\\n",
"a_2\n",
"\\end{pmatrix}\n",
"=\n",
"\\begin{pmatrix}\n",
"y_0\\\\\n",
"y_1\\\\\n",
"y_2\\\\\n",
"y_3\\\\\n",
"y_4\\\\\n",
"y_5\n",
"\\end{pmatrix},\n",
"$$\n",
"\n",
"which led to us developing the least squares solution?\n",
"\n",
".\n",
".\n",
".\n",
"\n",
"The answer is over- vs under-determined problems - the curve-fitting example is one where we have more equations than unknowns (over-determined), whereas this tomography example is one where we have more unknowns than equations (under-determined). We'll consider these different types of problems in more detail in the next few lectures."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Linear (or matrix) systems (a reminder on square systems)\n",
"\n",
"Consider the problem written in the form\n",
"\n",
"$$A\\boldsymbol{x}=\\boldsymbol{b},$$\n",
"\n",
"where for a given $m\\times n$ matrix $A$ and a given $m\\times 1$ right hand side vector $\\boldsymbol{b}\\;$ we want to find an $n\\times 1$ vector $\\boldsymbol{x}\\;$.\n",
"\n",
"Let's start by considering the square problem, i.e. where $m=n$. \n",
"\n",
"This was the case we considered in the *Numerical Methods module* for example where we considered topics such as the use of Gaussian elimination to solve problems of this type.\n",
"\n",
"In this module we are particularly interested in the case where the matrix may not be square, i.e. where the vectors $\\boldsymbol{x}$ and $\\boldsymbol{b}$ are of different lengths.\n",
"\n",
"In all situations it's important that we understand **whether our problem has a solution** and **whether that solution is unique**."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## A simple square example\n",
"\n",
"Recall that we can re-write a system of simultaneous (linear) equations in matrix form. \n",
"\n",
"For example, consider the following example:\n",
"\n",
"\\begin{eqnarray*}\n",
" 2x + 3y &=& 7 \\\\[5pt]\n",
" x - 4y &=& 3,\n",
"\\end{eqnarray*} \n",
"\n",
"this is completely equivalent to the following problem written in matrix form:\n",
"\n",
"$$\n",
"\\begin{pmatrix}\n",
" 2 & 3 \\\\\n",
" 1 & -4 \n",
"\\end{pmatrix} \n",
"\\begin{pmatrix}\n",
" x \\\\\n",
" y \n",
"\\end{pmatrix} =\n",
"\\begin{pmatrix}\n",
" 7 \\\\\n",
" 3 \n",
"\\end{pmatrix}.\n",
"$$\n",
"\n",
"Indeed when you see or are given a linear system, i.e. an $A$ and a $\\boldsymbol{b}$, it's often helpful to think about the set of linear equations it's representing. \n",
"[This, I think, may help us appreciate what row operations do (and are allowed) and why we need to treat the RHS vector consistently.]\n",
"\n",
"A collection of multiple linear equations for multiple unknowns is the general situation that will lead us to need to solve a matrix system; we may easily have a scenario where we have billions of unknowns."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## A simple example\n",
"\n",
"Let's go back to our simple example\n",
"\n",
"\\begin{eqnarray*}\n",
" 2x + 3y &=& 7 \\\\[5pt]\n",
" x - 4y &=& 3\n",
"\\end{eqnarray*} \n",
"\n",
"We can rearrange the first of these to\n",
"\n",
"$$y = -\\frac{2}{3} x + \\frac{7}{3}$$\n",
"\n",
"which we recognise as the equation of a straight line in 2D space: $y=mx +c$ where $m$ is the slope and $c$ is the intercept.\n",
"\n",
"Similarly the second equation can be rearranged into the standard form for a straignt line in 2D.\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"## Expanding the geometrical thinking\n",
"\n",
"Let's think through another interpretation which further utilises geometrical thinking. \n",
"\n",
"Consider all of $(x,y)$ space where $x$ and $y$ are independent variables that are free to vary arbitrarily - they will map out, or cover/span, all of 2D space. \n",
"\n",
"By stating that $ 2x + 3y = 7 $ we are now saying that $x$ and $y$ are not both independent - they cannot vary arbitrarily, one is effectively dependent on the other. \n",
"\n",
"For a linear system this one equation *restricts*, or *constrains*, us to a line in 2D (with a linear equation containing $x$, $y$ and $z$ variables constraining us to a plane in 3D etc). \n",
"\n",
"Indeed we can think about the equation as a constraint, we might even simply call it a constraint rather than an equation.\n",
"\n",
"For nonlinear problems we can still think of the model as a mapping, and still consider this in the form of a series of constraints over what potential output values our input data gets mapped to.\n",
"\n",
"Above we wrote this as $y\\equiv y(x)$ ($y$ is a function of $x$) to emphasise a point, but of course we could have just as easily rearranged to $x\\equiv x(y)$. If instead of $2x + 3y = 7$ we had the equation $2x = 7$ then of course we couldn't rearrange into the form $y\\equiv y(x)$, but trivially can into the form $x\\equiv x(y)$ (a constant in this case) - this just reflects the fact that this particular equation example represents a vertical line in 2D space ($m=\\infty$)!\n",
"\n",
"Of course our second equation can also be interpreted as a line in 2D space.\n",
"\n",
"Let's go back to our original example and plot the two lines that are defined by our two equations/constraints."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"scrolled": false
},
"outputs": [
{
"data": {
"text/plain": [
""
]
},
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAVcAAAFBCAYAAADDvuyeAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAA45UlEQVR4nO3dfXzN9f/H8cdrY1PmMhqRVESukiGj2FxfFQkp5FtfSVGkIl3oCn0rlEqR+LEV+6YLkavEFr6pXOQyX5FS5CIk5ovN9v798R6Gbc62c87nc85e99vt3Oyc89nnvBzz3Pu8P+8LMcaglFLKu0KcLkAppYKRhqtSSvmAhqtSSvmAhqtSSvmAhqtSSvmAhqtSSvlAIacL8JcyZcqYypUre3z8sWPHKFq0qO8K8qJAqVXr9K5AqRMCp9a81LlmzZoDxpiyFzxhjCkQt6ioKJMbiYmJuTreSYFSq9bpXYFSpzGBU2te6gRWmywyR7sFlFLKBzRclVLKBzRclVLKBzRclVLKBzRclVLKBwrMUCyl8io9PZ1du3Zx7Ngxp0s5R4kSJdiyZYvTZXgkUGo9v87ChQtz+eWXU7x48VyfS8NVqYs4cOAAIkK1atUICXHPh72jR49SrFgxp8vwSKDUmrlOYwzHjx9n9+7dALkOWPf8pAAiUkREvheR9SKyWUReyOIYEZE3RWS7iGwQkXpO1KoKjsOHDxMZGemqYFW+JyJceumlVKhQgf379+f6+93Wcj0JNDfGJItIYWCFiCwwxnyb6Zh2QNWM203Auxl/KuUTaWlpFC5c2OkylEMuueQSUlNTc/19rvpVnDHhITnjbuGM2/lbJXQC4jKO/RYoKSLlvVnH7t0wd65XT6kCnIg4XYJySF7/7d3WckVEQoE1QBVggjHmu/MOqQD8nun+rozH9mRxrn5AP4DIyEiSkpI8qmHSpGtISKjG4cO/0Lv3zlz/HfwtOTnZ47+bkwK1zhIlSnD06FHnCspGWlqaK+vKSqDUml2dJ06cyP3PblZzYt1wA0oCiUCt8x6fB9yc6f4SIOpi58vN2gKpqca0arXHgDHPPmtMerrH3+qIYJ637YTz6/zxxx+dKeQijhw54nQJHguUWrOrM6efAbJZW8B1LdfTjDGHRSQJaAtsyvTULuDKTPcrAn9487ULFYJhw/5LpUrleOklOHkS/vUv0E+GSilPuarPVUTKikjJjK8vAVoC/z3vsDnAPRmjBhoBfxtjLugSyK/QUHjvPXjwQXj1VXj0UdCNclUgefnll2nQoAHFixenbNmy3HrrrWzatOni3+hlEyZMIDo6muLFi1O8eHGio6OZN2+e3+uoXLkyInLBrUOHDj55PVeFK1AeSBSRDcAqYLEx5gsR6S8i/TOOmQ/sALYDk4GHfFVMSAhMmACDBsH48TBwIKSn++rVlPKupKQkHnroIb755huWLl1KoUKFaNmyJYcOHcr3uf/xj3/w/PPPe3RsxYoVeeGFF1i7di2rV6+mefPmdO7cmQ0bNuS7jtxYtWoVe/bsOXNbu3YtIkL37t1984JZ9RUE4y0/67mmpxszdKgxYEzfvsakpeXqVD4XqH2ZbhWsfa5Hjx41ISEhZs6cOcYYYz766CMTFhZmfv311zPHPPLII+aaa64xe/fuzfFcffr0Mc8991yeay1VqpSZOHHimfsVKlQwY8eOPeeYDRs2mPDwcLN582aPXyc3Ro4caUqUKGGOHTuWbZ2n5aXP1W0tV1cSsX2uzz4L778P994LaWlOV6VU7hw9epT09HRKlSoFQNeuXalduzYjR44EYMyYMcycOZOFCxcSGRnpkxrS0tJISEggOTmZxo0bn3k8OjqaVatWnXPs4MGD6du3LzVq1Djn8dGjRxMREZHjbfny5TnWYYxhypQp9OrVi0svvdR7f8FMXHtBy21E4MUXoXBhGDECUlIgPt5e/FIFy+DBsG6df1+zbl144438nWPQoEHUrVuX6OhowI7fHD16NB06dODaa69l1KhRLF26lKpVq+a73vNt3ryZli1bcuLECSIiIvjss8+oXbv2meejo6N55513ztyfPXs2P/zwAx999NEF5+rfv/9FP8pXqFAhx+cXL17ML7/8Qt++fXP5N/GcRkMuPfsshIfDsGGQmgozZkBYmNNVKZWzIUOGsGLFClasWEFoaOiZx1u3bk2DBg145plnmDt3Lg0aNMjy+0ePHs3o0aPP3D958iQiwpgxY848tmDBAm655ZYsv79q1aqsW7eOw4cP88knn9CnTx+SkpKoVasWAI0aNeKxxx7j0KFDFC1alMcff5wRI0Zw2WWXXXCu0qVLU7p06Ty9D6dNnjyZBg0aULdu3XydJycarnkwdKgN1Ecfha5dYdYsG7iqYMhvC9LfHn30URISEkhMTOSaa64557mlS5eyfv16jDE5dgWc31ocNmwYFSpU4JFHHjnzWE6txbCwMKpUqQJA/fr1WbVqFa+//jpTpkwBICoqirCwMFavXs0PP/xAoUKFGDBgQJbnOj/os5JT0O/fv5/PP/+cCRMm5HiO/NJwzaPBg23ADhgAnTvDp5/CJZc4XZVS5xo0aBAJCQkkJSVRvXr1c55bv349Xbp04a233mLevHkMHz6cRYsWZXme81uLxYoVo3Tp0mcCM7fS09M5efLkmfvh4eHceOONzJ07l+nTpzNjxoxs13PIb7fAtGnTCA8Pp0ePHnmq3VMarvnw0EM2YPv1g1tvhc8/hwDYPVgVEAMGDCA+Pp7Zs2dTqlQp9u7dC0BERAQHDx6kffv2DBkyhPvuu4+GDRtSp04dkpKSiImJ8WodTz75JLGxsVSrVo2jR48yY8YMkpKSLhjrGh0dzfjx42nVqhUdO3bM9nz56RYwxvD+++/To0cPny+BqOGaT3372oC9915o3x6++AICYNlKVQCcvkDUokWLcx5/+OGHWbx4MR07dmTEiBEA1KpVi27dujF8+HBWrlzp1Tr27t3L/fffz759+yhRogR16tRhwYIFtGnT5pzj6tatS0hICOPGjfPq62eWlJTEtm3b+OCDD3z2GqdpuHrBPffYgO3VC9q2hfnzoUQJp6tSBZ3J5ZTCf//73x4fO23atFwd68li2R9++CEPPPAANWvW9PjcuRUbG5vr9yWvNFy9pEcPO0yrRw9o1QoWLYKM4YRKqWykp6fz559/Mm3aNDZu3JirgHc7nUTgRXfcYS9srV8PLVrAgQNOV6SUuy1btozy5cszbdo0PvnkkzMTHIKBtly97PSFrc6doXlz+OoruPxyp6tSyp1iYmJID9IFO7Tl6gNt28K8ebB9O8TEwB6vr9mllHI7DVcfadECFi6E336DZs1g1y6nK1JK+ZOGqw81bQpffgl799qA3en+HWOUUl6i4epjjRvbftdDh2zY/vyz0xUppfxBw9UPGjaEJUvg2DHbgv3pJ6crUkr5moarn9SrB4mJdqnCZs3gxx+drkgp5Usarn5Uuzac3p03JgY2bnSyGqWUL2m4+lmNGvD113a6bEwMrF3rdEVKKV/QcHXAddfZgI2IsEO2vv/e6YqUUt6m4eqQa6+FZcvs+gMtW8I33zhdkVLKm1wVriJypYgkisgWEdksIoOyOCZGRP4WkXUZtxFO1OoNV11lA7ZcOWjd2rZmlQpUo0ePRkQYOHCgI68/YcIE6tSpQ/HixSlevDjR0dEXrBnrT64KV+AU8Jgx5nqgETBARGpkcdxyY0zdjNuL/i3RuypWtKFaqRK0a2fHxCoVaL799lsmT55MnTp1HKuhYsWKvPLKK6xdu5bVq1fTvHlzOnfuzIYNGxypx1XhaozZY4xZm/H1UWALkPM2jkGgfHk7iqBKFejYERYscLoiFQxmzZpFeHg4OzNNDRw0aBDXXnst+/bt89rr/P333/Ts2ZMpU6ZkuapVxYoVefvtt895bOPGjRQpUoQfvTgmsVOnTrRr144qVapw3XXXMWrUKIoVK+b1xb895apwzUxEKgM3At9l8XS0iKwXkQUi4ruVdf3o8svtONgaNeyKWnPmOF2RCnRdu3aldu3ajBw5EoAxY8Ywc+ZMFi5cmONmhLnVr18/unbtSvPmzbN8Pjo6mjVr1pzz2ODBg+nbty81apz7wXT06NFERETkeFu+fPlFa0pLSyMhIYHk5GQaN26c979cPrhyyUERiQA+AQYbY46c9/Ra4CpjTLKItAdmA1lutC4i/YB+AJGRkSSdHmTqgeTk5Fwd7y0vvFCIoUPr0KVLBM8+u4Vmzf686Pc4VWtuBWqdJUqU4OjRo2fuhw8bRoifBymn167NyVdeOeextLS0c+rKyjPPPEO3bt2oWLEiY8aMYe7cuZQrV+6i3+epadOm8dNPP/HOO+9w9OhR0tLSSElJOef89erVY/LkyWce++KLL1i7di1Tpky5oI6ePXvSvn37HF/ziiuuyLb+zZs307JlS06cOEFERAQffvghlStX9vjvm917euLEiVz/7Iq/tjzwlIgUBr4AFhljLrqZjoj8CtQ3xuS4NHX9+vXN6tWrPa7DFxu1eervv+1+XN99B3FxcPfdOR/vZK25Eah1btmyheuvv/7sAYMHw7p1/i2qbt0L9vT2ZOsUgMaNG/P9998zd+5c2rVrl+UxzzzzDKNGjcrxPImJiee8L1u3buXmm29m+fLlZ3aWjYmJoVatWud0A3zzzTc0adKEgwcPUrRoUWrWrMnAgQMZPHjwRWvPrZSUFH777TcOHz7MJ598wuTJk0lKSqJWrVoefX927+kFPwOZiMgaY0z98x93VctVRASYAmzJLlhFpBywzxhjRKQhtmvjoB/L9LkSJexyhR07Qu/ekJoKffo4XZU647yQc7OlS5eyfv16jDE5dgUMHjyYXr165XiuSpUqnXN/5cqVHDhw4JzgSktLY9myZUycOJFjx44RHh5OVFQUYWFhrF69mh9++IFChQoxYMCALF9j9OjRjB49Osc6FixYwC233JLlc2FhYWe2+65fvz6rVq3i9ddfZ8qUKTme0xdcFa5AE6A3sFFE1mU89hRQCcAYMxHoCjwoIqeA40AP47bmtxcUK2Y3OuzUye4sm5pqd5pVylPr16+nS5cuvPXWW8ybN4/hw4ezaNGiLI8tU6YMZcqUydX5O3fuTP365zbY7r33XqpWrcpTTz1FWFgYAOHh4dSpU4e5c+cyffp0ZsyYQeHChbM8Z//+/enevXuOr1uhgufXuNPT0zl58qTHx3uTq8LVGLMCkIsc8zbwdk7HBIuiRWHuXOjSBe6/3y768tBDTlelAsHOnTtp3749Q4YM4b777qNhw4bUqVPHq10zJUuWpGTJkuc8VrRoUUqXLn3Bx/CGDRsyYcIEWrVqRceOHbM9Z+nSpSldunSe6nnyySfp0KEDV155JUePHmXGjBkkJSU5NtbVVeGqLnTJJTB7NnTrBgMG2ID1QVeVCiKHDh2ibdu2dOzYkREj7BybWrVq0a1bN4YPH+7I0KTatWsTEhLCuHEXvYySZ3v37qVXr17s3buXEiVKUKdOHRYsWECbNm189po50XANAOHh8PHH9sLWo4/CyZMwbJjTVSm3Kl26NFu2bLngcX9sW53dFfWPPvqIBx54gJo1fTdyctq0aT47d15ouAaIsDBISIB77oEnn7QBOyJgJ/6qgiA9PZ0///yTadOm8eOPP/LJJ584XZJfabgGkEKFID7eBu1zz9kugpdecroqpbK2bNkymjdvTrVq1fjggw+ynL0VzDRcA0xoKEydagN21Cjbgr3ImGulHBETE0N6ejqA1yYtBBIN1wAUEgITJ9qAHTMGduyoQkwMSI7jLJRS/uTatQVUzkJC4K237AWuTz+tyEMPQUYjQSnlAtpyDWAiMHYs7N37GxMnViIlBd57z3YdKKWcpeEa4ETg/vt3ULVqJV580c7kmjrVXvxS3mOMQbTfpUDK6wRQ/S8YBETghRdsH+wzz9hRBPHxkM0MQ5VLoaGhpKamnpnOqQqW48ePZztdNycarkHk6afthIMnnrABm5BgA1flT8mSJdm3bx8VKlQgJEQvUxQUxhiOHz/O7t2787T+rYZrkHn8cRuogwbBHXfArFlQpIjTVQW2MmXKsGvXLrZu3ep0Kec4ceIERQLkHzdQaj2/zsKFCxMZGUnx4sVzfS4N1yD0yCM2YB980K6qNXu2XaNA5U1ISMgFy+25QVJSEjfeeKPTZXgkUGr1Zp36GSdI9e8PU6bA4sV2Xdhjx5yuSKmCRcM1iN13n93JICnJ7ixbACfJKOUYDdcg16sXzJgB33wDbdrYLWSUUr6n4VoA3HmnvbC1ejW0bAmHDjldkVLBT8O1gLj9dvj0U9iwAVq0gAM5bueolMovDdcCpGNHmDMH/vtfiI2Fffucrkip4KXhWsC0aQPz5sGOHRATA3/84XRFSgUnDdcCqHlzu3X3rl3QrBn8/rvTFSkVfDRcC6hbboEvv4T9+6FpU/j1V6crUiq4uCpcReRKEUkUkS0isllEBmVxjIjImyKyXUQ2iEg9J2oNBtHRsGQJHD5sA3b7dqcrUip4uCpcgVPAY8aY64FGwAARqXHeMe2Aqhm3fsC7/i0xuNSvD4mJ8L//2S4Cl02fVypguSpcjTF7jDFrM74+CmwBKpx3WCcgzljfAiVFpLyfSw0qdevaWVynTtmA3bzZ6YqUCnyuCtfMRKQycCPw3XlPVQAyX4LZxYUBrHKpVi34+mu7fUxMDKxf73RFSgU2yesq274kIhHA18AoY8yn5z03D3jZGLMi4/4SYKgxZk0W5+mH7TogMjIyKiEhweMakpOTiYiIyPtfwo+8WeuuXZcwZMgNnDwZyquvrqdatWSvnBcC5z3VOr0vUGrNS52xsbFrjDH1L3jCGOOqG1AYWAQMyeb5ScBdme5vBcpf7LxRUVEmNxITE3N1vJO8XevPPxtz1VXGlChhzLffeu+8gfKeap3eFyi15qVOYLXJInNc1S0gdpOiKcAWY8y4bA6bA9yTMWqgEfC3MWaP34osAK65xnYRXHYZtGoF//mP0xUpFXhcFa5AE6A30FxE1mXc2otIfxHpn3HMfGAHsB2YDDzkUK1B7aqrYNkyKF/ezupKSnK6IqUCi6t2IjC2HzXHLTYzmuED/FNRwVahgm3BtmgB7dvbdQlatnS6KqUCg9tarsplypWzrdaqVe3CL/PnO12RUoFBw1VdVNmysHQp1KwJnTvD5587XZFS7qfhqjxy2WV2quyNN0LXrnbxbaVU9jRclcdKlrQbHt50E/ToYbePUUplTcNV5Urx4na5wqZN7f5c06c7XZFS7qThqnItIsIuuN2yJfzjHzB5stMVKeU+Gq4qTy691A7NatcO+vWDCROcrkgpd9FwVXlWpAh89hl06gQDB8K47ObUKVUAabiqfAkPtyMHunaFxx6Dl192uiKl3MFVM7RUYCpcGGbOhLAweOopSEmBESNAcpxrp1Rw03BVXlGoEMTF2YB9/nk4eRJGjdKAVQWXhqvymtBQmDLFdhW8/LJtwb72mgasKpg0XJVXhYTAu+/aFuzYsTZgx4/XgFUFj4ar8joRG6iZA/add5yuSin/0nBVPiFiuwTCws52EfTs6XRVSvmPhqvyGRF7USs83F7k+v3362nWzF78UirY6ThX5VMi8NxzMHo0fPVVJHffDampTlellO9pG0L5xfDh8Pvv23n33SqkpkJCgm3RKhWstOWq/KZ79128+SbMng1dusCJE05XpJTvaLgqv3r4YZg0yW4Xc9tt8L//OV2RUr6h4ar8rl8/mDoVvvoKOnSA5GSnK1LK+zRclSPuvRfi4+323e3awZEjTleklHe5LlxFZKqI7BeRTdk8HyMif4vIuozbCH/XqLyjZ097YWvlSmjdGg4fdroipbzHdeEKTAPaXuSY5caYuhm3F/1Qk/KRbt3g449h7Vq7s8GhQ05XpJR3uC5cjTHLAP0vVoB07mwX3d60CZo3hz//dLoipfLPdeHqoWgRWS8iC0SkptPFqPzr0MFuG7N1K8TGwt69TlekVP6IMcbpGi4gIpWBL4wxtbJ4rjiQboxJFpH2wHhjTNVsztMP6AcQGRkZlZCQ4HENycnJRERE5KV8vwuUWj2pc926kgwfXpuyZU8yduw6ypZN8VN1ZwXT++kWgVJrXuqMjY1dY4ypf8ETxhjX3YDKwCYPj/0VKHOx46KiokxuJCYm5up4JwVKrZ7WuXy5MRERxlx7rTE7d/q2pqwE2/vpBoFSa17qBFabLDIn4LoFRKSciF0dVEQaYrs2DjpblfKmm2+GxYvhwAFo1gx++cXpipTKPdeFq4jMBFYC1URkl4j8U0T6i0j/jEO6AptEZD3wJtAj47eHCiKNGsGSJfD339C0KWzf7nRFSuWO6xZuMcbcdZHn3wbe9lM5ykFRUZCYaIdoNW0KS5dC9epOV6WUZ1zXclUqsxtugKQkSE+3XQSbspxaopT7aLgq16tZ0wZsaKgdprVundMVKXVxGq4qIFSvbtchuOQSO9Fg9WqnK1IqZxquKmBUqQJffw0lSkCLFvDtt05XpFT2NFxVQLn6ahuwZctCq1awfLnTFSmVNQ1XFXAqVbJdBBUqQNu2dkSBUm6j4aoC0hVX2Bbs1VdD+/bw5ZdOV6TUuTRcVcCKjLSt1mrV4NZbYd48pytS6iwNVxXQypa1kwtq14bbb7dLFyrlBhquKuCVLm3344qKsotvf/SR0xUppeGqgkTJkrbfNToa7roLPvzQ6YpUQafhqoJGsWKwcKGdJtu7N/zf/zldkSrINFxVUClaFL74wo6Bve8+mDTJ6YpUQaXhqoLOpZfC55/brWP694e33nK6IlUQabiqoFSkCHz6qd388JFHYOxYpytSBY2GqwpaYWF25ED37vD44zB6tNMVqYLEdYtlK+VNhQvbkQNhYfD003DyJDz/PNiNgpTyHQ1XFfQKFYJp02zQvvgipKTYVqwGrPIlDVdVIISGwvvv2xbsv/5lW7Bjx2rAKt/RcFUFRkgIvPuuDdjXX7ct2DfftI8r5W0arqpAEYHx4yE8HMaMsQE7caIGrPI+DVdV4IjAq6/aFuzo0TZgp0yxXQdKectFf1+LiF83MxaRqSKyX0Sy3OdTrDdFZLuIbBCRev6sTwUHERg50o4cmD4d7rkHTp1yuioVTDz5MPSDiIwXkVI+r8aaBrTN4fl2QNWMWz/gXT/UpIKQCDz3HLz8MsyYYRd8SU11uioVLDwJ14ZATWCbiDwsIj798GSMWQYcyuGQTkCcsb4FSopIeV/WpILbk0/CuHHw8cfQtSukpOgQApV/F+1zNcZsBFqKSGfgNeBBEXnMGLPA18VlowLwe6b7uzIe2+NMOSoYPPqo7YMdOBD27atF06Z2Cq1SeSXGGM8PFgkDHgWeAv4DDDHG/NfrRYlUBr4wxtTK4rl5wMvGmBUZ95cAQ40xa7I4th+264DIyMiohIQEj2tITk4mIiIib38BPwuUWgOhzi++KM+4cddRr95fjBy5iSJF0p0uKVuB8H6eFii15qXO2NjYNcaY+hc8YYzx+AaUBFoCbwJpQErG1yVycx4PXqcysCmb5yYBd2W6vxUof7FzRkVFmdxITEzM1fFOCpRaA6XOYcN+NCEhxsTEGHP0qNPVZC9Q3k9jAqfWvNQJrDZZZI4nowUGi8iHIvITcBCYCzQAxgN9gWrAjyJyU67iPu/mAPdkjBpoBPxtjNEuAeU1bdvuIz4eli+3W3cfOeJ0RSoQeTLO9TFgJfaq/LfAGmNMSqbn40RkGDAVe+ErX0RkJhADlBGRXcBzQGEAY8xEYD7QHtgO/A+4N7+vqdT57r7b9sHedZddeHvhQijlr/EyKih4ckHrSg/O83+AVxZ0M8bcdZHnDTDAG6+lVE66drUB27UrtGxp9+i67DKnq1KBwluT/v4EmnvpXEq5xm232V0NNm+G5s1h/36nK1KBwivhmtGv+7U3zqWU27RrZ/fl2rYNYmNhj/bwKw/ochVKeaBlS1iwAHbuhJgY2L3b6YqU22m4KuWhZs1g0SLbcm3a1AatUtnRcFUqF5o0ga++goMHbdju2OF0RcqtNFyVyqWGDWHpUjh61Abstm1OV6TcSMNVqTyoVw8SE+12MU2bwpYtTlek3EbDVak8qlMHkpLAGNuC3bjR6YqUm2i4KpUPNWrA11/bnWVjY+GHH5yuSLmFhqtS+VStGixbBkWL2okGq1Y5XZFyAw1Xpbzg2mttC7ZUKTsmduVKpytSTtNwVcpLKle2AXv55dC6tW3NqoJLw1UpL7ryShuqFSvaabNLljhdkXKKhqtSXla+vB1FcM010LGjXa5QFTwarkr5QGSkHQdbvTp06gRz5zpdkfI3DVelfKRMGdstUKcOdOkCn37qdEXKnzRclfKh0qXtWgQNGkD37vDvfztdkfIXDVelfKxECbuaVpMmdvuY+HinK1L+oOGqlB8UKwbz59u1YPv0gSlTnK5I+ZqGq1J+UrSo3dGgdWvo2xfefdfpipQvabgq5UeXXAKzZ9shWg89BOPHO12R8hUNV6X8rEgR+OQTO4Jg8GB47TWnK1K+4LpwFZG2IrJVRLaLyJNZPB8jIn+LyLqM2wgn6lQqP8LCICEBevSAoUNh5EinK1LeVsjpAjITkVBgAtAK2AWsEpE5xpgfzzt0uTGmo98LVMqLCheGDz6wfz77LKSkwAsvgIjTlSlvcFW4Ag2B7caYHQAikgB0As4PV6WCQmgo/N//2ZbsSy/ZnQ3+9S8NWEcYQ+jx4147ndvCtQLwe6b7u4CbsjguWkTWA38AjxtjNvujOKV8ITQU3nvPBuyrr9oW7LhxGrB+s327HXwcH881tWvbFXe8wG3hmtWPkznv/lrgKmNMsoi0B2YDVbM8mUg/oB9AZGQkSUlJHheSnJycq+OdFCi1ap0569YN9u+vwhtvVOSXX3bzyCPbCMnhqkigvJ/gvloLHTnC5YmJRC5eTInNmzEi/BUVxZ4aNdjmrTqNMa65AdHAokz3hwPDL/I9vwJlLnbuqKgokxuJiYm5Ot5JgVKr1nlx6enGPPGEMWBM377GpKVlf2ygvJ/GuKTWkyeN+ewzY26/3ZjChe2bXKuWMa+8YsyuXcaYvNUJrDZZZI7bWq6rgKoicjWwG+gB3J35ABEpB+wzxhgRaYgd8XDQ75Uq5QMi8MorEB5uRxCkpMDUqbbrQOWBMfDddxAXZxd2OHTILlk2cCD07g116/qs/8VV4WqMOSUiA4FFQCgw1RizWUT6Zzw/EegKPCgip4DjQI+M3x5KBQURe3ErLAxGjIDUVJsNhVz1v9XlfvnFDsWIj4dt2+zg4ttvt4HaqpVf3kzX/XMZY+YD8897bGKmr98G3vZ3XUr527PP2oB98knbgp0xw95X2fj7b5g1ywbq6T12YmLsG9i1KxQv7tdyXBeuSqmzhg2zXQSPPmrzYdYse19lSE21S47Fx8Pnn9uxbNWrw6hR0LMnXHWVY6VpuCrlcoMH2xbrgAHQubNddPuSS5yuykHGwNq1NlBnzIA//7Qrk/frZz/216/vinFsGq5KBYCHHrIB268f3HqrbaQVOL//Dh9+aDugt2yxb8htt9lAbdvWdX0mGq5KBYi+fW1+3HsvtG8Pw4YVgCEER4/apnpcnN2UzBi76vikSXZgcKlSTleYLQ1XpQLIPffYgO3VC4YOrUOTJnang6By6pTdfCw+3gbr8eN2K93nnrN/8WuvdbpCj2i4KhVgevSwAdu9ezFatbLXc1zcgPPc+vU2UD/8EPbuhZIl7bYNvXtDdLQr+lFzw3VLDiqlLq5LF3jxxc2sXw8tWsCBA05XlEd//AFjxsANN9gB/W++CTfdZBe83bvXbtfQuHHABStouCoVsBo3Psjnn9trO7GxsH+/0xV56Ngx2zpt0wauvBKeeMIOf3j7bRu2s2fb3x4BPuZMuwWUCmBt29p9uW67zY6XX7IEypd3uqospKfD0qX2Y//HH0Nysh2DOny47Ui+7jqnK/Q6DVelAlyLFrBggR1B0KyZDdgrr3S6qgw//gjx8TSaOtU2rYsXh+7dbaDecgs5LvsV4DRclQoCTZvCl1/apUibNbONxMqVHSpm/36YOdO2UtesgdBQjjVoQJE337RN7AIyAyJ4f20oVcA0bgxffQV//WUD9uef/fjix4/bVac6doQrrrDTyoyB11+H3bvZ+PLLcOedBSZYQcNVqaDSoIFttR47ZluzW7f68MXS0+0CKfffD+XK2TFi69bBY4/Bxo221Tp4sF3irwDSbgGlgsyNN9rJTC1bnr3IVaOGF19g27Yz26Lw669QtCjccYcdjxobq4vPZtCWq1JBqHZtOL1bSUwMbNiQzxMePAjvvGMH8193nV116rrrbMDu2wfTp9s012A9Q8NVqSB1/fX2U3tYmG1Qrl2byxOcPAmffWYXmS5f3i7LdewYvPaaXURl0SI7HbVoUZ/UH+i0W0CpIFa1qg3Y2Fg7ZGvRImjYMIdvMAa+/fbstih//WX7TB9++Oy2KMojGq5KBblrrrEB27y5/eS+cKEdWXCOHTvObouyfbu9qt+5s1+3RQk2+o4pVQBcdRV8/bVtvbZuDfPmQbMbDp/dFmX5cntgTAw89ZS9QOXnbVGCjYarUgVExYqQtDiVl5os5ECLeNJC5hCamrEtyujRdluUSpWcLjNoaLgqFeyMsWNO4+Ion5DAO3/+yaHQMkxM70fdN+6hySNRAbnqlNtpuCoVrH77za4+FR9vl84KD7d7xNxzD6ZBW6a0L8zmoTDrajsrVXmX64ZiiUhbEdkqIttF5MksnhcReTPj+Q0iUs+JOpVypaNHYdo0e/WqcmXbf3rZZfDee3Z91Fmz4NZbuaxcYZYssRf/77jDLp+qvMtV4SoiocAEoB1QA7hLRM6fW9IOqJpx6we869cilXKbU6fsEICePe2wqXvvta3W556zCwwsX26nqJYsec63lSoFixfboVl33mnXWlHe47ZugYbAdmPMDgARSQA6AT9mOqYTEGeMMcC3IlJSRMobY/b4v1ylnFN0+3aYO9duL713r03LPn3scn6NGnnUj1q8uB372rGjnQ+QkmJPofLPbeFaAfg90/1dwE0eHFMB0HBVwe+PP870ozbYuBEKF4YOHex41A4d8rR6f0QEzJ8PnTrZRm9qqt1pVuWP28I1q1+1Jg/H2ANF+mG7DoiMjCTp9GRrDyQnJ+fqeCcFSq1aZ96EHD9O2RUriPzyS0qtXYukp3Pk+uvZ2b8/R9q2JfX09q8rV+brdR5/PIQjR2py//2XsWnTT3Tu/IcXqrfc9p5mx6t1GmNccwOigUWZ7g8Hhp93zCTgrkz3twLlL3buqKgokxuJiYm5Ot5JgVKr1pkLp04Z89VXxvTpY0zRosaAMZUrG/PMM8Zs3WqM8U2dJ04Yc9tt9uXGjfPeeV3xnnogL3UCq00WmeO2lusqoKqIXA3sBnoAd593zBxgYEZ/7E3A30b7W1Ww2Lz57PbSu3bZTtG77rIf+2++2efbooSH2wEFd98NQ4bYPthhw3z6kkHLVeFqjDklIgOBRUAoMNUYs1lE+mc8PxGYD7QHtgP/A+51ql6lvOL0tihxcXbpqtBQu1/L2LF2XKqfV+8PC4OEBHtd7Mkn7eJYI0b4tYSg4KpwBTDGzMcGaObHJmb62gAD/F2XUl51/DjMmWMDddEiSEuDqCh44w3bUr38ckfLK1TINqDDwuyIrpQUeOklnciVG64LV6WCVno6rFhhA3XWLDhyxE74f/xx+7G/Zk2nKzxHaChMnWoHJIwaZVuwr76qAespDVelfO2nn85ui7Jzp11cumtX+7m7WTNXr94fEgKTJtkW7JgxtgX7xhsasJ7QcFXKFw4csItNx8XB99/blGrZ0jYBO3cOqNX7Q0Lg7bdtwL7xhh0H+/bbPr+2FvA0XJXylpMn7UKpcXF2VH5qqt3MaswYe/m9fHmnK8wzERg3zo4meOUV24KdNMnVjW7HabgqlR/G2MH78fFnt0UpVw4eecT2o95wg9MVeo0IvPyybcG+9JIN2KlTdZOC7OjbolReZLUtyu23237UFi2CNnFE4MUXbcA++6wN2Ph4e9FLnSs4fwKU8oW//jq7LcqKFTZpYmLg6aftun3Fijldod8884ztIhg61PZ+zJxpA1edpeGqVE5SU+1yfnFxdlxqSords3r0aLuM1JVXOl2hY554wgbsoEF28MOsWXlaNyZoabgqdT5jYPVqG6gzZ9or/2XKwAMP2PX46tXTsUgZHnnEtlgffNCuqvXZZ36fUOZaGq5Knfbbb/DBBzSYNMl+HR5u9z/p3RvattWOxWz0728Dtm9fuy7snDkBNdLMZzRcVcF25Ijd4yQuDjKWmkutU8d2KnbrdsHq/Spr991nf/f84x92WYR58wpUF3SWNFxVwXPqlN3fJD4eZs+28/yrVrWXwXv2ZN1vvxETE+N0lQGnd2/bgu3ZE1q3hgULCvbvJg1XVTAYA+vXn13Ob98+KF3aLr3fuzfcdNPZftTffnO21gB25522BdujB7RqZdekKV3a6aqcoeGqglumbVE4vS1Kx442UNu318vbPtClC3z6qR2d1qKF/ZBQEOnsYBV8jh2zA/xbt7ZDpYYOtVdYJkyAPXvs//zbb9dg9aHTF7b++187FPjQoYJ3MVBbrio4pKVBYqJtoX7yiQ3YypXtAP/evW2fqvKrNm3sha1bb4VHH61Lo0ZwxRVOV+U/Gq4qsG3aZAP1gw9sF0CJEnaRlN69oUkTXbrJYc2b2zkYbdqE06wZLF1acOZd6E+eCjz79tm17+rVs6tOjRtnv/7oI9i7F957D265RYPVJW65BcaM2cD+/Xb52l9/dboi/9CfPhUYjh+3Gzt16AAVKsCjj9rwHD8edu+GuXPtuNQiRZyuVGWhRo0jLFkChw9D06Z2rZtgp90Cyr3S02H5cjvA/+OPz26LMnSo/dh//fVOV6hyoX592y3QqhVnugiqVXO6Kt/RcFXus3Xr2X7UnTshIsKO6+nTx/6v1I/7AatuXXvdsUUL+0+5ZInrtg7zGv0pVe5w4IDdO+Smm6B6dbsqc/XqZwf8T5sGsbEarEGgVi34+mv7TxkTY+d2BCPX/KSKSGkRWSwi2zL+LJXNcb+KyEYRWSciq/1dp/KikyftsKnOne0WKA8/bB8bMwZ27bKXme++Gy691OlKlZdVrw7LltkVtGJjYc0apyvyPteEK/AksMQYUxVYknE/O7HGmLrGmPr+KU15jTHwn//YpZTKl7cLgX7/vV0UdN06e3vssYDeb0p5pkoV24ItXtx2E3z3ndMVeZebwrUTMD3j6+lAZ+dKUV73889cNX26Hcx/8822T7V9e9s6/f1321oNov2mlGeuvtq2YMuUsRe6VqxwuiLvcVO4Rhpj9gBk/Hl5NscZ4EsRWSMi/fxWncq9v/6yW4TefDNUqULl6dPtrKlp0+x41A8+sNN4dAvRAq1SJduCLV/eLpubmOh0Rd4hxhj/vZjIV0C5LJ56GphujCmZ6di/jDEX9LuKyBXGmD9E5HJgMfCwMWZZNq/XD+gHEBkZGZWQkOBxrcnJyURERHh8vJPcVKukplL6++8p9+WXXLZyJSGpqRy76ir2tmnDL40aUejqq50u8aLc9H7mJFDqBM9qPXQojCFDbmDv3iKMHLmJ+vX/8lN1Z+XlPY2NjV2TZRelMcYVN2ArUD7j6/LAVg++53ngcU/OHxUVZXIjMTExV8c7yfFa09ON+e47YwYONOayy4wBY8qWNWbQIGNWr7bPGxfU6SGt0/s8rXX/fmPq1DEmPNyYefN8W1NW8vKeAqtNFpnjpm6BOUCfjK/7AJ+ff4CIFBWRYqe/BloDm/xWoTrXzp0wapQdzH/TTTB5MrRsCV98YWdNvfEGREXpflPKY2XL2skFNWvaQSSfX5ACgcNN4fovoJWIbANaZdxHRK4QkfkZx0QCK0RkPfA9MM8Ys9CRaguqI0dg6lQ7fqZyZbsdyuWX2/n8e/eenaKq+02pPLrsMju5oF69s7vKBiLXzNAyxhwEWmTx+B9A+4yvdwB6SdnfTm+LEhdnt0U5ccJe9X/pJbunRwD0o6rAUrIkfPml/T3do4fd4fzuu52uKndcE67KZYyxY07j42HGjLPbovzzn3Zef8OG+nFf+VTx4nYfrltvhV69ICXFboAYKDRc1bl27z67LcqmTedui9Khg92BTik/iYiwC2537my3O0tJgX4BMgBTw1VBcjJ89pn92L9kiW21RkfDu+9C9+4Fd4c55QqXXmq3jLnjDnjgARuwAwc6XdXFabgWVGlp9rJsfLzdU+rYMdt3+uyz9jOYbouiXKRIEftjeueddgmKlBQYMsTpqnKm4VrQbNpkW6gffnh2W5S77rLL+TVpov2oyrXCw+3IgZ497fITJ0/C8OFOV5U9DdeCYO9emDnThuq6dVCokJ1n+MYb9mqBrt6vAkThwvb6auHC8NRTtgU7YoQ72wQarsHqf/+zI7Dj4+2YlrQ0uxT8+PG2pVq2rNMVKpUnhQrZdkJYGDz/vG3BjhrlvoDVcA0m6el2iaHT26IcPWq32hw6FO65xy6iqVQQCA2FKVNswL78sm3BvvaauwJWwzUIXPrbb/D003aVqd9+s+NXunWzw6d0WxQVpEJCYOJEG7Bjx9qAHT/ePQGr4RqoDhywU03j4mi4apX9SWvd2v4a79xZV+9XBYIIvPmmvdg1dqztInj3XXe0JzRcA8nJk3ZRlLg4mD/fTkutW5ftDz5IlWef1dX7VYEkYrsEwsNh9Gjbgn3/feeXCdZwdbvT26LEx8NHH9mN38uXh8GD7cf+OnXYlZREFQ1WVYCJwMiRNmCfe86uRTBtmr345RQNV7favt32ocbHw44d9mN+ly42UFu0cP7XslIuI2KHZWUepvXhh84t0Kbh6iaHDtnWaVwcrFxpf1qaN7e/im+/HYoVc7pCpVxv+HDbgn3sMduCTUiw9/1Nw9VpKSl26Z+4ONufmpICNWrAv/5lp6JUrOh0hUoFnCFD7CiChx+2axJ8/LH/58pouDrBGLuddHy8/bV68KAd1P/QQ3Y8at267hlPolSAGjjQBuwDD0CnTnZtIn8OotFw9aedO20/alwc/PST/VXaqZMN1FatdPV+pbysXz/73+qf/7QrZ86dC0WL+ue1NVx97e+/7WeS+Hi7fzBA06YwbJj9vFKihLP1KRXk7r3XBmyfPnZJjfnz/XP5QsPVF06dsvP54+PPboty3XV2W5RevezeU0opv+nVy3YR3H23nWuzYIHdSsaXNFy9xRj44Yez26Ls328Xmb7vPvsrs0ED7UdVykHdu9uA7d7dblL85Ze+XQdewzW/du06uy3K5s32X+/0tijt2+u2KEq5SOfO9sLWHXfYUY6LF/tugTgN17xITrbLosfF2dX8jYHGjXVbFKUCQIcOdtuYTp3sDvFffQXlynn/dVywvIElIt1EZLOIpItI/RyOaysiW0Vku4g86bcC09Ls54jevSEy0n7U37HDTgnZts1OUe3fX4NVqQDQurXd+PCXXyAmxu7L6W1uarluAroAk7I7QERCgQlAK2AXsEpE5hhjfvRZVRs32o/8p7dFKVnS9o737q3boigVwJo3h4ULbe9ds2b2Q6g3uablaozZYozZepHDGgLbjTE7jDEpQALQyevF7N1LxVmz4MYboU4deP11u4r/rFmwZw9MmgQ336zBqlSAu+UW2+964IAN2D17vDeNS4wxXjuZN4hIEvC4MWZ1Fs91BdoaY/pm3O8N3GSMyXKjXRHpB/QDiIyMjEpISPCohqsnT+aqGTM4Ur06+1q3Zn9sLKm+HreRD8nJyURERDhdxkVpnd4VKHWC+2vdujWCJ564gWrVDvHaa1ty9b2xsbFrjDEXdGX6tVtARL4Csuo6ftoY87knp8jisWx/Oxhj3gPeA6hfv76JiYnxpEy47jq+b9OGhvfcQ3HA7ZtMJyUl4fHfzUFap3cFSp3g/lpjYuw16Z9/3u61Ov0arsaYlvk8xS7gykz3KwJ/5POcF7riCv5XqZLXT6uUcq/ateHgwVSvnc81fa4eWgVUFZGrRSQM6AHMcbgmpZS6gGvCVURuF5FdQDQwT0QWZTx+hYjMBzDGnAIGAouALcBHxpjNTtWslFLZcc1QLGPMZ8BnWTz+B9A+0/35wHw/lqaUUrnmmparUkoFEw1XpZTyAQ1XpZTyAQ1XpZTyAQ1XpZTyAQ1XpZTyAQ1XpZTyAdct3OIrIvInsDMX31IGOOCjcrwtUGrVOr0rUOqEwKk1L3VeZYy5YD+DAhOuuSUiq7Na6caNAqVWrdO7AqVOCJxavVmndgsopZQPaLgqpZQPaLhm7z2nC8iFQKlV6/SuQKkTAqdWr9Wpfa5KKeUD2nJVSikf0HDNgafbfTvFsW3Gc0lEporIfhHZ5HQtORGRK0UkUUS2ZPy7D3K6pqyISBER+V5E1mfU+YLTNeVEREJF5AcR+cLpWnIiIr+KyEYRWSciF+zhl1sarjk7vd33MqcLOV+mbcbbATWAu0SkhrNVZWsa0NbpIjxwCnjMGHM90AgY4NL39CTQ3BhzA1AXaCsijZwtKUeDsIvbB4JYY0xdbwzH0nDNgYfbfTvFP9uMe4ExZhlwyOk6LsYYs8cYszbj66PYQKjgbFUXMlZyxt3CGTdXXjwRkYpAB+B9p2vxNw3XwFUB+D3T/V24MAgClYhUBm4EvnO4lCxlfNReB+wHFhtjXFkn8AYwFEh3uA5PGOBLEVkjIv3yezLXbPPiFC9s9+2UXG0zrjwnIhHAJ8BgY8wRp+vJijEmDagrIiWBz0SkljHGVX3aItIR2G+MWSMiMQ6X44kmxpg/RORyYLGI/DfjU1eeFPhw9cJ2307xzzbjBYyIFMYG64fGmE+drudijDGHRSQJ26ftqnAFmgC3iUh7oAhQXEQ+MMb0criuLGXs14cxZr+IfIbtestzuGq3QODSbca9TEQEmAJsMcaMc7qe7IhI2YwWKyJyCdAS+K+jRWXBGDPcGFPRGFMZ+/O51K3BKiJFRaTY6a+B1uTzl5WGaw6y2+7bDQJpm3ERmQmsBKqJyC4R+afTNWWjCdAbaJ4xHGddRqvLbcoDiSKyAftLdrExxtXDnAJAJLBCRNYD3wPzjDEL83NCnaGllFI+oC1XpZTyAQ1XpZTyAQ1XpZTyAQ1XpZTyAQ1XpZTyAQ1XpZTyAQ1XpZTyAQ1XpZTyAQ1XpTizMPpJEbkq02PjReRnEYl0sjYVmHSGllKcWVdgFfCDMeZ+EXkcu1ReE2PMNmerU4GowK+KpRTYBahF5CnsGhI/A09jV/vXYFV5oi1XpTIRkW+wS83daoxZ4HQ9KnBpn6tSGUSkOXADdiHyfQ6XowKctlyVAkTkBuBrYAh2z6cIY0wbZ6tSgUzDVRV4GSMEvgEmGWNeFJFawAZsn2uSo8WpgKXhqgo0ESkN/AdYZox5INPj/wYqGWOiHStOBTQNV6WU8gG9oKWUUj6g4aqUUj6g4aqUUj6g4aqUUj6g4aqUUj6g4aqUUj6g4aqUUj6g4aqUUj6g4aqUUj7w/03/1Lem16HwAAAAAElFTkSuQmCC\n",
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"x = np.linspace(-1,5,100)\n",
"y1 = -(2./3.)*x + (7./3.)\n",
"y2 = (1./4.)*x - (3./4.)\n",
"\n",
"fig = plt.figure(figsize=(5, 5))\n",
"\n",
"ax1 = fig.add_subplot(111)\n",
"\n",
"ax1.set_xlabel(\"$x$\", fontsize=14)\n",
"ax1.set_ylabel(\"$y$\", fontsize=14)\n",
"ax1.grid(True)\n",
"\n",
"ax1.plot(x,y1,'b', label='$2x+3y=7$')\n",
"ax1.plot(x,y2,'r', label='$x-4y=3$')\n",
"\n",
"ax1.legend(loc='best', fontsize=14)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The blue line maps out all combinations of $x$ and $y$ values that satisfy the first equation/constraint, and the red all the $x$ and $y$ values that satisfy the second.\n",
"\n",
"So in solving the system of simultaneous equations, or equivalently solving the corresponding matrix system, we are asking the question what $x$ and $y$ values satisfy BOTH equations (or constraints)?\n",
"\n",
"\n",
"From the image above we see that **there is a solution** (at the intersection) and it is clearly **unique** (there is only one intersection), for these particular equations. This is the only combination of $x$ and $y$ values that satisfy BOTH equations/constraints.\n",
"\n",
"\n",
"In 2D for this example, one constraint restricts our 2D space of $x,y$ values to a 1D subset, the other to another 1D subset, combined they restrict the 2D space down to 0D, i.e. a single point - our unique solution.\n",
"\n",
"\n",
"But of course in 2D there are two other possible scenarios - what are these? Can you construct suitable pairs of equations that demonstrate these possibilities, and the corresponding plots to help explain what is going in visually? .... see the exercises."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Solving a general linear system\n",
"\n",
"In this simple $2\\times 2$ case we can find the solution via the above geometric method, or we can substitute one equation into the other to solve the system:\n",
"\n",
"\\begin{eqnarray*}\n",
" 2x + 3y &=& 7 \\\\[5pt]\n",
" x - 4y &=& 3\n",
"\\end{eqnarray*} \n",
"\n",
"Lots of ways we could do this, but as $x$ is already on its own in the second equation perhaps easiest (as it delays the need to deal with fractions) to rearrange and substitute $x=3+4y$ into the first:\n",
"\n",
"$$ 2x + 3y = 7 \\quad\\implies\\quad 2(3+4y)+3y=7 \\quad\\implies\\quad \n",
"$$\n",
"$$(2\\times 4 + 3)y=7 - 2\\times 3 \\quad\\implies\\quad y = \\frac{1}{11} $$\n",
"\n",
"Then substitute this back into one of the equations to find $x$, again easiest option is the second:\n",
"\n",
"$$x=3+4y = 3+4\\times \\frac{1}{11} = \\frac{37}{11}$$\n",
"\n",
"But we can't really use either of these approaches (plotting or substitution) if our system is much larger in dimension - what if we had billions of unknowns (which on today's computers are problem sizes we would certainly be wanting to solve.)\n",
"\n",
"We can check our solution to the $2\\times 2$ problem by adding this point we just derived to our plot and comfirm if it agrees with the intersection."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"scrolled": false
},
"outputs": [
{
"data": {
"text/plain": [
""
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAVcAAAFBCAYAAADDvuyeAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAA+K0lEQVR4nO3dfXzN9f/H8cdr2FzMRStGxFIKITXUKDbXV0XCV6HSV0soUt9KSir0rahUKvryw6J9kxS5KmzRtysUufqK5DLLF8Xmcrb374/3zHA2uzjnfD5ne91vt3Ozc85nn89rx/Y87/P+vD/vtxhjUEop5V1BTheglFKFkYarUkr5gIarUkr5gIarUkr5gIarUkr5gIarUkr5QHGnC/CXyy67zEREROR6+6NHj1KmTBnfFeRFgVKr1uldgVInBE6t+alzzZo1B4wxFS94whhTJG6RkZEmLxISEvK0vZMCpVat07sCpU5jAqfW/NQJrDYeMke7BZRSygc0XJVSygc0XJVSygc0XJVSygc0XJVSygeKzFAspXJy5MgR9u/fT2pqKuXLl2fz5s1Ol3RRgVInBE6t59dZokQJKlWqRLly5fK8Lw1XVeQdOXKEP/74g6pVq1KqVClSUlIoW7as02VdVHJyckDUCYFTa9Y6jTEcP36cvXv3AuQ5YF3VLSAiJUXkBxFZJyIbReR5D9uIiLwpIttE5GcRudGJWlXhsX//fqpWrUrp0qUREafLUS4hIpQuXZqqVauyf//+PH+/21quJ4GWxpgUESkBfC0ii4wx32XZpgNQK+N2E/Buxr9K5UtqaiqlSpVyugzlUqVKlSI1NTXP3+eqlmvGBQ8pGXdLZNzOXyqhCzAjY9vvgAoiUsWbdezdC/Pne3WXyuW0xaqyk9/fDbe1XBGRYsAa4GpgojHm+/M2qQrsznJ/T8Zj+zzsKxaIBQgPDycxMTFXNUyaVJP4+Gv566/f6Nt3Z55/Bn9LSUnJ9c/mJLfWWb58eZKTkzPvp6WlnXPfrQKlTgicWrOr88SJE3n+3XVduBpj0oCGIlIBmCsi9YwxG7Js4ultxONCYMaYycBkgEaNGpno6Ohc1XDLLXDwYBJTp15J1apX8vzz4OaGTWJiIrn92Zzk1jo3b958zsmWQDz54naBUmt2dZYsWZIbbrghT/tyXbieYYz5S0QSgfZA1nDdA1yR5X414HdvHrt4cXjyyf9SvXplXnwRTp6Ef/7T3QGrlHIXV/W5ikjFjBYrIlIKaA3897zN5gH3ZIwauBk4bIy5oEugoIoVg8mT4aGH4JVX4NFHQRfKVW7z0ksv0bhxY8qVK0fFihW57bbb2LBhw8W/0csmTpxIgwYNKFeuHOXKlSMqKooFCxb4vY6IiAhE5IJbp06d/F6L21quVYDpGf2uQcBHxpjPRWQAgDHmPWAh0BHYBhwD+vmqmKAgmDgRgoNhwgRITYW33rKPK+UGiYmJDBw4kMaNG2OMYeTIkbRu3ZpNmzYRFhZWoH3fd999REREMGrUqItuW61aNV5++WVq1apFeno606dPp2vXrqxZs4YGDRoUqI68WLVqFWlpaZn39+3bR2RkJD179vRbDZk8zUNYGG8Fmc81Pd2YJ54wBozp39+YtLQ87crnCvNcmf6wadOmc+4fOXLEoUryxlOdycnJJigoyMybN88YY8xHH31kgoODzY4dOzK3eeSRR0zNmjVNUlJSjvu/9957zXPPPZfv+i655BLz3nvvZdZatWpVM378+HO2+fnnn01ISIjZuHFjvo+Tk9GjR5vy5cubo0eP5mr77P7vz/8dyQqdzzX/RGyf67PPwr/+Bf36QZY3R6VcIzk5mfT0dC655BIAunfvTv369Rk9ejQA48aN48MPP2Tx4sWEh4f7pIa0tDTi4+NJSUmhadOmmY9HRUWxatWqc7YdOnQo/fv3p27duuc8PnbsWEJDQ3O8rVy5Msc6jDFMmTKFPn36ULp0ae/9gLnktm4B1xKBF16AEiVg5Eg4dQri4uzJL1X4DB0Ka9f695gNG8IbbxRsH0OGDKFhw4ZERUUBdozm2LFj6dSpE1dddRVjxoxh+fLl1KpVq8D1nm/9+vVERUVx4sQJQkNDmTt3LvXr1898PioqinfeeSfz/qeffspPP/3ERx99dMG+BgwYcNGP8lWrVs3x+S+//JLffvuN/v375/En8Q6Nhjx69lkICYEnn7R9sLNm2T5ZpZw2bNgwvv76a77++muKFSuW+Xjbtm1p3LgxzzzzDPPnz6dx48Yev3/s2LGMHTs28/7JkycREcaNG5f52KJFi7j11ls9fv+1117L2rVr+euvv5gzZw733nsviYmJ1KtXD4Cbb76Zxx57jEOHDlGmTBkef/xxRo4cyaWXXnrBvsLCwgrcZ/z+++/TuHFjGjZsWKD95JeGaz488YQN1Ecfhe7dYfZsG7iq8ChoC9LfHn30UeLj40lISKBmzZrnPLd8+XLWrVuHMSbHroDzW4tPPvkkVatW5ZFHHsl8LKfWYnBwMFdffTUAjRo1YtWqVbz++utMmTIFgMjISIKDg1m9ejU//fQTxYsXZ9CgQR73dX7Qe5JT0O/fv5/PPvuMiRMn5rgPX9JwzaehQ23ADhoEXbvCJ5+AXp6unDBkyBDi4+NJTEykdu3a5zy3bt06unXrxltvvcWCBQsYPnw4S5Ys8bif81uLZcuWJSwsLDMw8yo9PZ2TJ09m3g8JCeGGG25g/vz5TJ8+nVmzZlGiRAmP31vQboFp06YREhJCr1698lW7N2i4FsDAgTZgY2Phttvgs88gAFYPVoXIoEGDiIuL49NPP+WSSy4hKSkJgNDQUA4ePEjHjh0ZNmwY999/P02aNKFBgwY+uVLuqaeeolOnTlxxxRUkJycza9YsEhMTLxjrGhUVxYQJE2jTpg2dO3fOdn8F6RYwxvCvf/2LXr16OXpVmIZrAfXvbwO2Xz/o2BE+/xwC4Co/VUicOUHUqlWrcx5/+OGH+fLLL+ncuTMjR44EoF69evTo0YPhw4fz7bfferWOpKQk+vTpQ1JSEuXLl6dBgwYsWrSIdu3anbNdw4YNCQoK4rXXXvPq8bNKTExk69atfPDBBz47Rm5ouHrBPffYgO3TB9q3h4ULoXx5p6tSRYHJ42WD//73v3O97bRp07y+7cyZM3nwwQe57rrrcr3vvIqJicnz6+ILGq5e0quXHabVqxe0aQNLlkDGUEOlirT09HT279/PO++8w/r16/MU8IFMLyLwojvvtCe21q2DVq3gwAGnK1LKeStWrKBWrVpMmzaNOXPmZF7gUNhpy9XLzpzY6toVWraEpUuhUiWnq1LKOdHR0Rw+fDggphz0Jm25+kD79rBgAWzbBtHRsM/rc3YppdxOw9VHWrWCxYth1y5o0QL27HG6IqWUP2m4+lDz5vDFF5CUZAN2p/tXjFFKeYmGq481bWr7XQ8dsmH7669OV6SU8gcNVz9o0gSWLYOjR20L9pdfnK5IKeVrGq5+cuONkJBgpyps0QI2bXK6IqWUL2m4+lH9+nBmdd7oaFi/3slqlFK+pOHqZ3Xrwldf2ctlo6Phxx+drkgp5Qsarg645hobsKGhdsjWDz84XZFSyts0XB1y1VWwYoWdf6B1a/jmG6crUkp5k6vCVUSuEJEEEdksIhtFZIiHbaJF5LCIrM24jXSiVm+oUcMGbOXK0Latbc0q5WZjx45FRBg8eLAjx584cSINGjSgXLlylCtXjqioqAvmjHULV4UrcBp4zBhTB7gZGCQidT1st9IY0zDj9oJ/S/SuatVsqFavDh062DGxSrnRd999x/vvv0+DBg0cq6FatWq8/PLL/Pjjj6xevZqWLVvStWtXfv75Z8dqyo6rwtUYs88Y82PG18nAZiDnJR4LgSpV7CiCq6+Gzp1h0SKnK1KBYvbs2YSEhLAzy+V/Q4YM4aqrruKPP/7w2nEOHz5M7969mTJlisdZrapVq3bBBNjr16+nZMmSbPLiuMMuXbrQoUMHrr76aq655hrGjBlD2bJlvT75tze4KlyzEpEI4Abgew9PR4nIOhFZJCK+m3XXjypVsuNg69a1M2rNm+d0RSoQdO/enfr16zN69GgAxo0bx4cffsjixYtzXIwwr2JjY+nevTstW7b0+HxUVBSrVq0657GhQ4fSv39/6tY998Pn2LFjCQ0NzfG2cuXKi9aUlpZGfHw8KSkpNG3aNP8/nI+4cspBEQkF5gBDjTFHznv6R6CGMSZFRDoCnwIeF2EXkVggFiA8PJzEM4NMcyElJSVP23vL888X54knGtCtWyjPPruZFi3+d9HvcarWvHJrneXLlyc5OTnzflpaGqcGDiTIzwOR0+vX5+TLL+d6+7S0NFJSUnjmmWfo0aMH1apVY9y4ccyfP5/KlSuf8zMVxLRp0/jll1945513SE5Otq/PqVPn7P/GG2/k/fffz3zs888/58cff2TKlCmZ33Pmud69e9OxY8ccj3n55ZdnW//GjRtp3bo1J06cIDQ0lJkzZxIREeGVnzdrnVmdOHEiz7+74oblELISkRLA58ASY8xFF9oRkR1AI2NMjlNTN2rUyKxevTrXdfhiEbfcOnzYrsf1/fcwYwbcfXfO2ztZa164tc7NmzdTp06dzPvJycmUffZZWLvWv4U0bJinNb2Tk5Mz50ht2rQpP/zwA/Pnz6dDhw4et3/mmWcYM2ZMjvtMSEg45/9oy5Yt3HLLLaxcuTJzZdno6Gjq1avH22+/nbndN998Q7NmzTh48CBlypThuuuuY/DgwQwdOvSCWgvq1KlT7Nq1i7/++os5c+bw/vvvk5iYSL169Qq87+zqPP93JCsRWWOMaXT+465quYqIAFOAzdkFq4hUBv4wxhgRaYLt2jjoxzJ9rnx5O11h587Qty+kpsK99zpdVRGTh5Bz2vLly1m3bh3GmBy7AoYOHUqfPn1y3Ff16tXPuf/tt99y4MCBc4IrLS2NFStW8N5773H06FFCQkKIjIwkODiY1atX89NPP1G8eHEGDRrk8Rhjx45l7NixOdaxaNEibr31Vo/PBQcHZy733ahRI1atWsXrr7/OlClTctynv7kqXIFmQF9gvYiszXjsaaA6gDHmPaA78JCInAaOA72M25rfXlC2rF3osEsXu7JsaqpdaVaprNatW0e3bt146623WLBgAcOHD2fJkiUet73sssu47LLL8rT/rl270qjRuY2yfv36UatWLZ5++mmCg4MBCAkJ4YYbbmD+/PlMnz6dWbNmUaJECY/7HDBgAD179szxuFWr5v48dnp6OidPnsz19v7iqnA1xnwNyEW2eRt4O6dtCosyZWD+fOjWDR54wE76MnCg01Upt9i1axcdO3Zk2LBh3H///TRp0oQGDRp4tfulQoUKVKhQ4ZzHypQpQ1hY2AUfw6OiopgwYQJt2rShc+fO2e4zLCyMsLCwfNXz1FNP0alTJ6644gqSk5OZNWsWiYmJrhzr6qpwVRcqVQo+/RR69IBBg2zAZnRjqSLs0KFDdOvWjc6dOzNypL2Opl69evTo0YPhw4c7MjSpYcOGBAUFXTAky5uSkpLo06cPSUlJlC9fngYNGrBo0SLatWvns2Pml4ZrAAgJgY8/tie2Hn0UTp6EJ590uirlpLCwMFavXn3ByRd/LFud3VnzmTNn8uCDD3Lddb4bHTlt2jSf7dvbNFwDRHAwxMfDPffAU0/ZgB0ZsBf+qsIiPT2d//3vf0ybNo3169f7JdwDhYZrACleHOLibNA+95ztInjxRaerUkXZihUraNmyJddeey1z5szxePVWUaXhGmCKFYOpU23AjhljW7AXGY+tlM9ER0eTnp7udBmupOEagIKC4L33bMCOGwfbt19NdDRIjuMslFL+5Nq5BVTOgoLgrbfsCa5PPqnGwIGgDQil3ENbrgFMBMaPh6SkXbz3XnVOnYLJk23XgVLKWRquAU4EHnhgO7VqVeeFF+yVXFOn2pNfSinn6J9gISACzz9v+2CfecaOIoiLg2yuPlRK+YGGayEyYoS94OAf/7ABGx9vA1cp5X96QquQefxxmDAB5s6FO++EEyecrkipoknDtRB65BF49134/HM7q9bx405XpFTRo+FaSA0YAFOmwJdf2nlhjx51uiKlrFGjRnllYmsAEeHjjz/2yr68TcO1ELv/fruSQWKiXVnWS6t+KBfZu3cvsbGxVKtWjeDgYKpWrcoDDzzAnj17nC7Nq+677z6P0xju27eP2267zYGKLk7DtZDr0wdmzYJvvoF27ewSMqpw2LFjB40aNWLDhg1Mnz6dbdu28cEHH7Bx40YaN27Mjh07CrT/U6dOeadQH6pcuTIhISFOl+GRhmsR8Le/wezZsHo1tG4Nhw45XVHhdGahvKCgICIiIpg5c6ZPj/fYY48RFBTE0qVLadWqFdWrVycmJoalS5cSFBR0zjIr0dHRDB48+JzvP781GB0dzUMPPcTjjz9OxYoVadasmcfj7t69my5duhAWFkbp0qWpXbs28fHxmc+vX7+e1q1bU6pUKcLCwrjvvvs4nMO7uqdWadaug1GjRjF9+nQWLFiAiCAimdMent8tcLFjnznWhAkTqFq1Kpdccgn9+vXj2LFj2daXXxquRcQdd8Ann8DPP0OrVnAgx+UcVV7NnDmT2NhYdu7ciTGGnTt3Ehsb67OAPXToEEuXLmXQoEGULl36nOdKly7NwIEDWbRoEX/++Wee9vvBBx9gjGHlypXMmDHD4zYDBw7k2LFjJCQksHHjRt54443M1QqOHTtG+/btCQ0N5YcffmDu3Ll888032a6nlRuPP/44PXv2pHXr1uzbt499+/Z5XEo7u2Pff//952y3cuVKNmzYwNKlS/n3v//N3LlzmTBhQr7ry46Ocy1COneGefOga1eIiYGlS8GLS9sXaSNGjLig9XPs2DFGjBhB7969vX68rVu3YozJdkXSunXrYoxh69atNGnSJNf7vfLKKxk/fnyO2+zcuZM777yT66+/PvN7zpg5cyYpKSnExcVlTuQ9efJkYmJi2LZtW+bCgnkRGhpKqVKlCAkJoXLlytlul9tjlytXjnfffZfixYtTp04devTowbJlyxg+fHiea8uJtlyLmHbtYMEC2L4doqPh99+drqhw2LVrV54e9xbJZiq0M2t2Zvd8diIjIy+6zZAhQxg9ejRRUVE888wzrFmzJvO5zZs306BBg3NWSGjatClBQUFs2rQpT7XkVW6PXbduXYpnuT788ssvZ//+/V6vR8O1CGrZ0i7dvWcPtGgBu3c7XVHgO39J6os9XlC1atVCRNi4caPH5zdv3oyIcNVVVwEQFBTE+Yskp6amXvB9ZcqUueix//73v/Pbb7/Rr18/fvnlF5o2bcqoUaMAG+rZBXp2j+e2tovJ7bHPX5VWRHwyJ62GaxF1663wxRewfz80bw4FPLFc5I0ZM8Zj3+eYMWN8crywsDBatWrFO++847E7YuLEiXTo0CFzldWKFSuyb9++c7Zbt25dvo9frVo1YmNj+eijj3jhhReYPHkyYFuF69atIznLuL9vvvmG9PT0bLswPNW2du3ac+4HBweTlpaWY035ObYvuSpcReQKEUkQkc0islFEhnjYRkTkTRHZJiI/i8iNTtRaGERFwbJl8NdfNmC3bXO6osDVu3dvJk+eTI0aNRARatSoweTJk33S33rGuHHjOH36NK1bt2b58uXs3r2bxMRE2rRpgzGGt98+uwJ9y5YtWbRoEfPmzWPLli0MGzaM3fn8yDJkyBAWL17M9u3bWbt2LYsXL6Zu3bqAfR3KlCnDPffcw/r161mxYgUPPvggt99+e7b9rS1btuSnn35i6tSpbNu2jVdeeYX//Oc/52wTERHBhg0b2LJlCwcOHPDYss3u2N26dctXX29BuSpcgdPAY8aYOsDNwCARqXveNh2AWhm3WOBd/5ZYuDRqBAkJcOyY7SLYssXpigJX79692bFjB+np6ezYscOnwQpQs2ZNVq9ezXXXXUffvn2pWbMmd999N3Xq1GHVqlXnnGi6//77M2/NmjUjNDSUO+64I1/HTU9P5+GHH6Zu3bq0adOG8PBwpk+fDtjW+pIlSzhy5AhNmjShS5cuREVFMXHixGz3165dO5577jlGjBhBZGQkO3bsYODAgeds88ADD1CnTh0aNWpExYoVLwjfnI49derUfP2cBWaMce0N+Axoc95jk4C7stzfAlS52L4iIyNNXiQkJORpeyd5o9b1642pVMmY8HBjNmwoeE2euPU13bRp0zn3jxw54lAleRModRoTOLVmV+f5vyNZAauNh8xxW8s1k4hEADcA35/3VFUg6+eZPRmPqQKoVw+++souHxMdDQXojlNK4dJxriISCswBhhpjjpz/tIdvMR4eQ0RisV0HhIeHZ17VkRspKSl52t5J3qz1lVdKMWzY9TRvXoxXXlnHtdemeGW/4N7XtHz58uecBElLSzvnvlsFSp0QOLVmV+eJEyfy/rvrqTnr5A0oASwBhmXzvHYLnMfbtf76qzE1ahhTvrwx333nvf269TXVbgHfC5RaC223gNjBaFOAzcaY17LZbB5wT8aogZuBw8aYfdlsq/KhZk3bRXDppdCmDXg4d6CUughXhSvQDOgLtBSRtRm3jiIyQEQGZGyzENgObAPeBwZmsy9VADVqwIoVUKWKvarLhZ/mvcoYjz1LSuX7d8NVfa7GmK/x3KeadRsD5H8WCJVrVavaFmyrVtCxo52XoHVrp6vyvhIlSnD8+PELLgJQCuD48eMXXNWVG25ruSqXqVzZtlpr1bITvyxc6HRF3lepUiX27t3LsWPHtAWrMhljOHbsGHv37qVSpUp5/n5XtVyVO1WsCMuXQ9u2dkat2bPt2lyFRbly5QD4/fffSU1N5cSJE5QsWdLhqi4uUOqEwKn1/DpLlChBeHh45u9IXmi4qly59FJ7qWy7dtC9u13doEcPp6vynnLlymX+ASUmJnLDDTc4XNHFBUqdEDi1erNO7RZQuVahgl3w8KaboFcvG7BKKc80XFWelCtnpyts3tyuz5VxSblS6jwarirPQkPthNutW8N998H77ztdkVLuo+Gq8qV0aTs0q0MHiI2FHCY9UqpI0nBV+VayJMyda0cODB4Mr2V3TZ1SRZCGqyqQkBA7NKt7d3jsMXjpJacrUsoddCiWKrASJeDDDyE4GJ5+Gk6dgpEjIY9r4ylVqGi4Kq8oXhxmzLABO2oUnDwJY8ZowKqiS8NVeU2xYjBliu0qeOkl24J99VUNWFU0abgqrwoKgnfftS3Y8eNtwE6YoAGrih4NV+V1IjZQswbsO+84XZVS/qXhqnxCxHYJBAef7SLw8WKoSrmKhqvyGRF7UiskxJ7k2r27Di1a2JNfShV2Os5V+ZQIPPccjB0LS5eGc/fdkJrqdFVK+Z62IZRfDB8Ou3dv4913ryY1FeLjbYtWqcJKW67Kb3r23MObb8Knn0K3bnDihNMVKeU7Gq7Krx5+GCZNssvF3H47HDvmdEVK+YaGq/K72FiYOhWWLoVOnSAlxemKlPI+DVfliH79IC7OLt/doQMcOeJ0RUp5l+vCVUSmish+EdmQzfPRInJYRNZm3Eb6u0blHb172xNb335rFz/86y+nK1LKe1wXrsA0oP1FtllpjGmYcXvBDzUpH+nRAz7+GH780a5scOiQ0xUp5R2uC1djzApA/8SKkK5d7aTbGzZAy5bwv/85XZFSBee6cM2lKBFZJyKLROQ6p4tRBdepk102ZssWiImBpCSnK1KqYMQY43QNFxCRCOBzY0w9D8+VA9KNMSki0hGYYIyplc1+YoFYgPDw8Mj4+Phc15CSkkJoaGh+yve7QKk1N3WuXVuB4cPrU7HiScaPX0vFiqf8VN1Zhen1dItAqTU/dcbExKwxxjS64AljjOtuQASwIZfb7gAuu9h2kZGRJi8SEhLytL2TAqXW3Na5cqUxoaHGXHWVMTt3+rYmTwrb6+kGgVJrfuoEVhsPmRNw3QIiUlnEzg4qIk2wXRsHna1KedMtt8CXX8KBA9CiBfz2m9MVKZV3rgtXEfkQ+Ba4VkT2iMjfRWSAiAzI2KQ7sEFE1gFvAr0y3j1UIXLzzbBsGRw+DM2bw7ZtTlekVN64buIWY8xdF3n+beBtP5WjHBQZCQkJdohW8+awfDnUru10VUrljutarkpldf31kJgI6em2i2CDx0tLlHIfDVfletddZwO2WDE7TGvtWqcrUuriNFxVQKhd285DUKqUvdBg9WqnK1IqZxquKmBcfTV89RWULw+tWsF33zldkVLZ03BVAeXKK23AVqwIbdrAypVOV6SUZxquKuBUr267CKpWhfbt7YgCpdxGw1UFpMsvty3YK6+Ejh3hiy+crkipc2m4qoAVHm5brddeC7fdBgsWOF2RUmdpuKqAVrGivbigfn244w47daFSbqDhqgJeWJhdjysy0k6+/dFHTleklIarKiQqVLD9rlFRcNddMHOm0xWpok7DVRUaZcvC4sX2Mtm+feH//s/pilRRpuGqCpUyZeDzz+0Y2Pvvh0mTnK5IFVUarqrQKV0aPvvMLh0zYAC89ZbTFamiSMNVFUolS8Inn9jFDx95BMaPd7oiVdRouKpCKzjYjhzo2RMefxzGjnW6IlWUuG6ybKW8qUQJO3IgOBhGjICTJ2HUKLALBSnlOxquqtArXhymTbNB+8ILcOqUbcVqwCpf0nBVRUKxYvCvf9kW7D//aVuw48drwCrf0XBVRUZQELz7rg3Y11+3Ldg337SPK+VtGq6qSBGBCRMgJATGjbMB+957GrDK+zRcVZEjAq+8YluwY8fagJ0yxXYdKOUtF32/FhG/LmYsIlNFZL+IeFznU6w3RWSbiPwsIjf6sz5VOIjA6NF25MD06XDPPXD6tNNVqcIkNx+GfhKRCSJyic+rsaYB7XN4vgNQK+MWC7zrh5pUISQCzz0HL70Es2bZCV9SU52uShUWuQnXJsB1wFYReVhEfPrhyRizAjiUwyZdgBnG+g6oICJVfFmTKtyeegpeew0+/hi6d4dTp3QIgSq4i/a5GmPWA61FpCvwKvCQiDxmjFnk6+KyURXYneX+nozH9jlTjioMHn3U9sEOHgx//FGP5s3tJbRK5ZcYY3K/sUgw8CjwNPAfYJgx5r9eL0okAvjcGFPPw3MLgJeMMV9n3F8GPGGMWeNh21hs1wHh4eGR8fHxua4hJSWF0NDQ/P0AfhYotQZCnZ9/XoXXXruGG2/8k9GjN1CyZLrTJWUrEF7PMwKl1vzUGRMTs8YY0+iCJ4wxub4BFYDWwJtAGnAq4+vyedlPLo4TAWzI5rlJwF1Z7m8Bqlxsn5GRkSYvEhIS8rS9kwKl1kCp88knN5mgIGOio41JTna6muwFyutpTODUmp86gdXGQ+bkZrTAUBGZKSK/AAeB+UBjYALQH7gW2CQiN+Up7vNvHnBPxqiBm4HDxhjtElBe0779H8TFwcqVdunuI0ecrkgFotyMc30M+BZ7Vv47YI0x5lSW52eIyJPAVOyJrwIRkQ+BaOAyEdkDPAeUADDGvAcsBDoC24BjQL+CHlOp8919t+2DvesuO/H24sVwib/Gy6hCITcntK7IxX7+D/DKhG7GmLsu8rwBBnnjWErlpHt3G7Ddu0Pr1naNrksvdboqFSi8ddHf/4CWXtqXUq5x++12VYONG6FlS9i/3+mKVKDwSrhm9Ot+5Y19KeU2HTrYdbm2boWYGNinPfwqF3S6CqVyoXVrWLQIdu6E6GjYu9fpipTbabgqlUstWsCSJbbl2ry5DVqlsqPhqlQeNGsGS5fCwYM2bLdvd7oi5VYarkrlUZMmsHw5JCfbgN261emKlBtpuCqVDzfeCAkJdrmY5s1h82anK1Juo+GqVD41aACJiWCMbcGuX+90RcpNNFyVKoC6deGrr+zKsjEx8NNPTlek3ELDVakCuvZaWLECypSxFxqsWuV0RcoNNFyV8oKrrrIt2EsusWNiv/3W6YqU0zRclfKSiAgbsJUqQdu2tjWrii4NV6W86IorbKhWq2Yvm122zOmKlFM0XJXysipV7CiCmjWhc2c7XaEqejRclfKB8HA7DrZ2bejSBebPd7oi5W8arkr5yGWX2W6BBg2gWzf45BOnK1L+pOGqlA+Fhdm5CBo3hp494d//droi5S8arkr5WPnydjatZs3s8jFxcU5XpPxBw1UpPyhbFhYutHPB3nsvTJnidEXK1zRclfKTMmXsigZt20L//vDuu05XpHxJw1UpPypVCj791A7RGjgQJkxwuiLlKxquSvlZyZIwZ44dQTB0KLz6qtMVKV9wXbiKSHsR2SIi20TkKQ/PR4vIYRFZm3Eb6USdShVEcDDEx0OvXvDEEzB6tNMVKW8r7nQBWYlIMWAi0AbYA6wSkXnGmE3nbbrSGNPZ7wUq5UUlSsAHH9h/n30WTp2C558HEacrU97gqnAFmgDbjDHbAUQkHugCnB+uShUKxYrB//2fbcm++KJd2eCf/9SAdYQxFDt+3Gu7c1u4VgV2Z7m/B7jJw3ZRIrIO+B143Biz0R/FKeULxYrB5Mk2YF95xbZgX3tNA9Zvtm2zg4/j4qhZv76dcccL3Baunn6dzHn3fwRqGGNSRKQj8ClQy+PORGKBWIDw8HASExNzXUhKSkqetndSoNSqdeasRw/Yv/9q3nijGr/9tpdHHtlKUA5nRQLl9QT31Vr8yBEqJSQQ/uWXlN+4ESPCn5GR7Ktbl63eqtMY45obEAUsyXJ/ODD8It+zA7jsYvuOjIw0eZGQkJCn7Z0UKLVqnReXnm7MP/5hDBjTv78xaWnZbxsor6cxLqn15Elj5s415o47jClRwr7I9eoZ8/LLxuzZY4zJX53AauMhc9zWcl0F1BKRK4G9QC/g7qwbiEhl4A9jjBGRJtgRDwf9XqlSPiACL78MISF2BMGpUzB1qu06UPlgDHz/PcyYYSd2OHTITlk2eDD07QsNG/qs/8VV4WqMOS0ig4ElQDFgqjFmo4gMyHj+PaA78JCInAaOA70y3j2UKhRE7Mmt4GAYORJSU202FHfVX6vL/fabHYoRFwdbt9rBxXfcYQO1TRu/vJiu++8yxiwEFp732HtZvn4beNvfdSnlb88+awP2qadsC3bWLHtfZePwYZg92wbqmTV2oqPtC9i9O5Qr59dyXBeuSqmznnzSdhE8+qjNh9mz7X2VITXVTjkWFweffWbHstWuDWPGQO/eUKOGY6VpuCrlckOH2hbroEHQtauddLtUKaercpAx8OOPNlBnzYL//c/OTB4baz/2N2rkinFsGq5KBYCBA23AxsbCbbfZRlqRs3s3zJxpO6A3b7YvyO2320Bt3951fSYarkoFiP79bX706wcdO8KTTxaBIQTJybapPmOGXZTMGDvr+KRJdmDwJZc4XWG2NFyVCiD33GMDtk8feOKJBjRrZlc6KFROn7aLj8XF2WA9ftwupfvcc/YHv+oqpyvMFQ1XpQJMr142YHv2LEubNvZ8josbcLm3bp0N1JkzISkJKlSwyzb07QtRUa7oR80L1005qJS6uG7d4IUXNrJuHbRqBQcOOF1RPv3+O4wbB9dfbwf0v/km3HSTnfA2Kcku19C0acAFK2i4KhWwmjY9yGef2XM7MTGwf7/TFeXS0aO2ddquHVxxBfzjH3b4w9tv27D99FP77hHgY860W0CpANa+vV2X6/bb7Xj5ZcugShWnq/IgPR2WL7cf+z/+GFJS7BjU4cNtR/I11zhdoddpuCoV4Fq1gkWL7AiCFi1swF5xhdNVZdi0CeLiuHnqVNu0LlcOeva0gXrrreQ47VeA03BVqhBo3hy++MJORdqihW0kRkQ4VMz+/fDhh7aVumYNFCvG0caNKfnmm7aJXUSugCi8bxtKFTFNm8LSpfDnnzZgf/3Vjwc/ftzOOtW5M1x+ub2szBh4/XXYu5f1L70Ef/tbkQlW0HBVqlBp3Ni2Wo8eta3ZLVt8eLD0dDtBygMPQOXKdozY2rXw2GOwfr1ttQ4daqf4K4K0W0CpQuaGG+zFTK1bnz3JVbeuFw+wdWvmsijs2AFlysCdd9rxqDExOvlsBm25KlUI1a8PZ1YriY6Gn38u4A4PHoR33rGD+a+5xs46dc01NmD/+AOmT7dprsGaScNVqUKqTh37qT042DYof/wxjzs4eRLmzrWTTFepYqflOnoUXn3VTqKyZIm9HLVMGZ/UH+i0W0CpQqxWLRuwMTF2yNaSJdCkSQ7fYAx8993ZZVH+/NP2mT788NllUVSuaMtVqUKuZk0bsGFh9pP7qFEziYiIICgoiIiICGbOnAnbt8MLL9iP+k2b2o/57dvDwoWwZw+MH6/BmkfaclWqCKhRA776Cho3nsnzz8cCxwDYuXMnsffcA+np9AbbQfv00/YElZ+XRSlsNFyVKiKqVYNixZ7mTLCecSw9nREVKtB73TqoXt2Z4gohDVelCjtj7JjTGTP4fe8uj5vsOnxYg9XLNFyVKqx27bKzT8XF2amzQkKoXro0O48du2DT6hqsXue6E1oi0l5EtojINhF5ysPzIiJvZjz/s4jc6ESdSrlScjJMmwYtW9rJBZ5+Gi69FCZPhqQkxkyeTOnSpc/7ptJ06TLGgWILN1e1XEWkGDARaAPsAVaJyDxjzKYsm3UAamXcbgLezfhXqaLp9Gk7qUBcnB2Xevy4XQrluefs8KmaNTM37d27NwAjRoxg165dVKtWnZIlxzBxYm9uvhnuusupH6LwcVW4Ak2AbcaY7QAiEg90AbKGaxdghjHGAN+JSAURqWKM2ef/cpVyTplt22D+fLu8dFKSXevl3nvtdH4335zt7P29e/fODFmwU6t27myvBzh1yu5CFZzbwrUqsDvL/T1c2Cr1tE1VQMNVFX6//57Zj9p4/XooUQI6dbIt1E6d8jV7f2ioHc7apYtdWTY11a40qwrGbeHq6a3W5GMbu6FILBALEB4eTuKZi61zISUlJU/bOylQatU68yfo+HEqfv014V98wSU//oikp3OkTh12DhjAkfbtST2z/Ou33xboOI8/HsSRI9fxwAOXsmHDL3Tt+rsXqrfc9ppmx6t1GmNccwOigCVZ7g8Hhp+3zSTgriz3twBVLrbvyMhIkxcJCQl52t5JgVKr1pkHp08bs3SpMffea0yZMsaAMRERxjzzjDFbthhjfFPniRPG3H67Pdxrr3lvv654TXMhP3UCq42HzHFby3UVUEtErgT2Ar2Au8/bZh4wOKM/9ibgsNH+VlVYbNx4dnnpPXvsVVJ33WU/9t9yi8+XRQkJgdmz4e67Ydgw2wf75JM+PWSh5apwNcacFpHBwBKgGDDVGLNRRAZkPP8esBDoCGzDXmrSz6l6lfKKM8uizJhhp64qVsyu1zJ+PNx2m99n7w8Ohvh4e17sqafs5FgjR/q1hELBVeEKYIxZiA3QrI+9l+VrAwzyd11KedXx4zBvng3UJUsgLQ0iI+GNN2xLtVIlR8srXtw2oIOD7YiuU6fgxRezHYCgPHBduCpVaKWnw9df20CdPRuOHLEX/D/+uP3Yf911Tld4jmLFYOpUOyBhzBjbgn3lFQ3Y3NJwVcrXfvnl7LIoO3fayaW7d7efu1u0cPXs/UFBMGmSbcGOG2dbsG+8oQGbGxquSvnCgQN2sukZM+CHH2xKtW5tm4BduwbU7P1BQfD22zZg33jDjoN9+22fn1sLeBquSnnLyZOwYIEN1IULbQrVr2+bfHffbZdKCVAi8NprdjTByy/bFuykSa5udDtOw1WpgjDGDt6Pizu7LErlyvDII7Yf9frrna7Qa0TgpZdsC/bFF23ATp1qT36pC+nLolR+bN8OH3xgQ3XbNjtc6o47bD9qq1aFNnFE7GowwcHw7LM2YOPi7Ekvda7C+RuglC/8+ac9yx8XZ8/6i9hlUUaMsMuilC3rdIV+88wztovgiSds78eHH9rAVWdpuCqVk9RUWLzY9qPOm2ebanXqwNixdhqpK65wukLH/OMfNmCHDLGDH2bPzte8MYWWhqtS5zMGVq+2gfrhh/bM/2WXwYMP2vn4brxRxyJleOQR22J96CE7q9bcuX6/oMy1NFyVOmPXLvjgAxpPmmS/DgmB22+3J6bat9eOxWwMGGADtn9/Oy/svHkBNdLMZzRcVdF25AjMmWNbqRlTzaU2aGA7FXv0gAoVHC0vUNx/v33vue8+Oy3CggVFqgvaIw1XVfScPg1ffmlPTH36qb3Ov1Ytexq8d2/W7tpFdHS001UGnL59bQu2d29o2xYWLSra700arqpoMAbWrTs7nd8ff0BYmJ16v29fuOmms/2ouzwvP60u7m9/sy3YXr2gTRs7J01YmNNVOUPDVRVuWZZF4cyyKJ0720Dt2FFPb/tAt27wySd2dFqrVvZDQlGkVwerwufoUTvAv21bO1TqiSfsGZaJE2HfPvuXf8cdGqw+dObE1n//a4cCHzpU9E4GastVFQ5paZCQYFuoc+bYgI2IsAP8+/a1farKr9q1sye2brsNHn20ITffDJdf7nRV/qPhqgLbhg02UD/4wHYBlC9vJ0np2xeaNdOpmxzWsqW9BqNduxBatIDly4vOdRf6m6cCzx9/2LnvbrzRzjr12mv2648+gqQkmDwZbr1Vg9Ulbr0Vxo37mf377fS1O3Y4XZF/6G+fCgzHj9uFnTp1gqpV4dFHbXhOmAB798L8+XZcasmSTleqPKhb9wjLlsFff0Hz5naum8JOuwWUe6Wnw8qVdoD/xx+fXRbliSfsx/46dZyuUOVBo0a2W6BNGzK7CK691umqfEfDVbnPli1n+1F37oTQUDuu59577V+lftwPWA0b2vOOrVrZ/8ply1y3dJjX6G+pcocDB+zaITfdBLVr21mZa9c+O+B/2jSIidFgLQTq1YOvvrL/ldHR9tqOwsg1v6kiEiYiX4rI1ox/L8lmux0isl5E1orIan/Xqbzo5Ek7bKprV7sEysMP28fGjYM9e+xp5rvvhtKlna5UeVnt2rBihZ1BKyYG1qxxuiLvc024Ak8By4wxtYBlGfezE2OMaWiMaeSf0pTXGAP/+Y+dSqlKFTsR6A8/2ElB1661t8ceC+j1plTuXH21bcGWK2e7Cb7/3umKvMtN4doFmJ7x9XSgq3OlKK/79VdqTJ9uB/PfcovtU+3Y0bZOd++2rdVCtN6Uyp0rr7Qt2Msusye6vv7a6Yq8x03hGm6M2QeQ8W+lbLYzwBciskZEYv1Wncq7P/+0S4TecgtcfTUR06fbq6amTbPjUT/4wF7Go0uIFmnVq9sWbJUqdtrchASnK/IOMcb472AiS4HKHp4aAUw3xlTIsu2fxpgL+l1F5HJjzO8iUgn4EnjYGLMim+PFArEA4eHhkfHx8bmuNSUlhdDQ0Fxv7yQ31SqpqYT98AOVv/iCS7/9lqDUVI7WqEFSu3b8dvPNFL/ySqdLvCg3vZ45CZQ6IXe1HjoUzLBh15OUVJLRozfQqNGffqrurPy8pjExMWs8dlEaY1xxA7YAVTK+rgJsycX3jAIez83+IyMjTV4kJCTkaXsnOV5rerox339vzODBxlx6qTFgTMWKxgwZYszq1fZ544I6c0nr9L7c1rp/vzENGhgTEmLMggW+rcmT/LymwGrjIXPc1C0wD7g34+t7gc/O30BEyohI2TNfA22BDX6rUJ1r504YM8YO5r/pJnj/fWjdGj7/3F419cYbEBmp602pXKtY0V5ccN11dhDJZxekQOBwU7j+E2gjIluBNhn3EZHLRWRhxjbhwNcisg74AVhgjFnsSLVF1ZEjMHWqHT8TEWGXQ6lUyV7Pn5R09hJVXW9K5dOll9qLC2688eyqsoHINVdoGWMOAq08PP470DHj6+2AnlL2tzPLosyYYZdFOXHCnvV/8UW7pkcA9KOqwFKhAnzxhX2f7tXLrnB+991OV5U3rglX5TLG2DGncXEwa9bZZVH+/nd7XX+TJvpxX/lUuXJ2Ha7bboM+feDUKbsAYqDQcFXn2rv37LIoGzacuyxKp052BTql/CQ01E643bWrXe7s1CmIDZABmBquClJSYO5c+7F/2TLbao2KgnffhZ49i+4Kc8oVSpe2S8bceSc8+KAN2MGDna7q4jRci6q0NHtaNi7Oril19KjtO332WfsZTJdFUS5SsqT9Nf3b3+wUFKdOwbBhTleVMw3XombDBttCnTnz7LIod91lp/Nr1kz7UZVrhYTYkQO9e9vpJ06ehOHDna4qexquRUFSEnz4oQ3VtWuheHF7neEbb9izBTp7vwoQJUrY86slSsDTT9sW7MiR7mwTaLgWVseO2RHYcXF2TEtamp0KfsIE21KtWNHpCpXKl+LFbTshOBhGjbIt2DFj3BewGq6FSXq6nWLozLIoycl2qc0nnoB77rGTaCpVCBQrBlOm2IB96SXbgn31VXcFrIZrIVB61y4YMcLOMrVrlx2/0qOHHT6ly6KoQiooCN57zwbs+PE2YCdMcE/AargGqgMH7KWmM2bQZNUq+5vWtq19G+/aVWfvV0WCCLz5pj3ZNX687SJ49113tCc0XAPJyZN2UpQZM2DhQntZasOGbHvoIa5+9lmdvV8VSSK2SyAkBMaOtS3Yf/3L+WmCNVzd7syyKHFx8NFHduH3KlVg6FD7sb9BA/YkJnK1BqsqwkRg9GgbsM89Z+cimDbNnvxyioarW23bZvtQ4+Jg+3b7Mb9bNxuorVo5/7aslMuI2GFZWYdpzZzp3ARtGq5ucuiQbZ3OmAHffmt/W1q2tG/Fd9wBZcs6XaFSrjd8uG3BPvaYbcHGx9v7/qbh6rRTp+zUPzNm2P7UU6egbl345z/tpSjVqjldoVIBZ9gwO4rg4YftnAQff+z/a2U0XJ1gjF1OOi7Ovq0ePGgH9Q8caMejNmzonvEkSgWowYNtwD74IHTpYucm8ucgGg1Xf9q50/ajzpgBv/xi30q7dLGB2qaNzt6vlJfFxto/q7//3c6cOX8+lCnjn2NruPra4cP2M0lcnF0/GKB5c3jySft5pXx5Z+tTqpDr188G7L332ik1Fi70z+kLDVdfOH3aXs8fF3d2WZRrrrHLovTpY9eeUkr5TZ8+tovg7rvttTaLFtmlZHxJw9VbjIGffjq7LMr+/XaS6fvvt2+ZjRtrP6pSDurZ0wZsz552keIvvvDtPPAargW1Z8/ZZVE2brT/e2eWRenYUZdFUcpFuna1J7buvNOOcvzyS99NEKfhmh8pKXZa9Bkz7Gz+xkDTprosilIBoFMnu2xMly52hfilS6FyZe8fxwXTG1gi0kNENopIuog0ymG79iKyRUS2ichTfiswLc1+jujbF8LD7Uf97dvtJSFbt9pLVAcM0GBVKgC0bWsXPvztN4iOtutyepubWq4bgG7ApOw2EJFiwESgDbAHWCUi84wxm3xW1fr19iP/mWVRKlSwveN9++qyKEoFsJYtYfFi23vXooX9EOpNrmm5GmM2G2O2XGSzJsA2Y8x2Y8wpIB7o4vVikpKoNns23HADNGgAr79uZ/GfPRv27YNJk+CWWzRYlQpwt95q+10PHLABu2+f9y7jEmOM13bmDSKSCDxujFnt4bnuQHtjTP+M+32Bm4wxHhfaFZFYIBYgPDw8Mj4+Plc1XPn++9SYNYsjtWvzR9u27I+JIdXX4zYKICUlhdDQUKfLuCit07sCpU5wf61btoTyj39cz7XXHuLVVzfn6XtjYmLWGGMu6Mr0a7eAiCwFPHUdjzDGfJabXXh4LNt3B2PMZGAyQKNGjUx0dHRuyoRrruGHdu1ocs89lAPcvsh0YmIiuf7ZHKR1eleg1AnurzU62p6T/vXXbV6r06/haoxpXcBd7AGuyHK/GvB7Afd5ocsv51j16l7frVLKverXh4MHU722P9f0uebSKqCWiFwpIsFAL2CewzUppdQFXBOuInKHiOwBooAFIrIk4/HLRWQhgDHmNDAYWAJsBj4yxmx0qmallMqOa4ZiGWPmAnM9PP470DHL/YXAQj+WppRSeeaalqtSShUmGq5KKeUDGq5KKeUDGq5KKeUDGq5KKeUDGq5KKeUDGq5KKeUDrpu4xVdE5H/Azjx8y2XAAR+V422BUqvW6V2BUicETq35qbOGMeaC9QyKTLjmlYis9jTTjRsFSq1ap3cFSp0QOLV6s07tFlBKKR/QcFVKKR/QcM3eZKcLyINAqVXr9K5AqRMCp1av1al9rkop5QPaclVKKR/QcM1Bbpf7dopjy4znkYhMFZH9IrLB6VpyIiJXiEiCiGzO+H8f4nRNnohISRH5QUTWZdT5vNM15UREionITyLyudO15EREdojIehFZKyIXrOGXVxquOTuz3PcKpws5X5ZlxjsAdYG7RKSus1VlaxrQ3ukicuE08Jgxpg5wMzDIpa/pSaClMeZ6oCHQXkRudrakHA3BTm4fCGKMMQ29MRxLwzUHuVzu2yn+WWbcC4wxK4BDTtdxMcaYfcaYHzO+TsYGQlVnq7qQsVIy7pbIuLny5ImIVAM6Af9yuhZ/03ANXFWB3Vnu78GFQRCoRCQCuAH43uFSPMr4qL0W2A98aYxxZZ3AG8ATQLrDdeSGAb4QkTUiElvQnblmmReneGG5b6fkaZlxlXsiEgrMAYYaY444XY8nxpg0oKGIVADmikg9Y4yr+rRFpDOw3xizRkSiHS4nN5oZY34XkUrAlyLy34xPXflS5MPVC8t9O8U/y4wXMSJSAhusM40xnzhdz8UYY/4SkURsn7arwhVoBtwuIh2BkkA5EfnAGNPH4bo8ylivD2PMfhGZi+16y3e4ardA4NJlxr1MRASYAmw2xrzmdD3ZEZGKGS1WRKQU0Br4r6NFeWCMGW6MqWaMicD+fi53a7CKSBkRKXvma6AtBXyz0nDNQXbLfbtBIC0zLiIfAt8C14rIHhH5u9M1ZaMZ0BdomTEcZ21Gq8ttqgAJIvIz9k32S2OMq4c5BYBw4GsRWQf8ACwwxiwuyA71Ci2llPIBbbkqpZQPaLgqpZQPaLgqpZQPaLgqpZQPaLgqpZQPaLgqpZQPaLgqpZQPaLgqpZQPaLgqRebE6CdFpEaWxyaIyK8iEu5kbSow6RVaSpE5r8Aq4CdjzAMi8jh2qrxmxpitzlanAlGRnxVLKbATUIvI09g5JH4FRmBn+9dgVfmiLVelshCRb7BTzd1mjFnkdD0qcGmfq1IZRKQlcD12IvI/HC5HBThtuSoFiMj1wFfAMOyaT6HGmHbOVqUCmYarKvIyRgh8A0wyxrwgIvWAn7F9romOFqcCloarKtJEJAz4D7DCGPNglsf/DVQ3xkQ5VpwKaBquSinlA3pCSymlfEDDVSmlfEDDVSmlfEDDVSmlfEDDVSmlfEDDVSmlfEDDVSmlfEDDVSmlfEDDVSmlfOD/AYdfKQe4kGYAAAAAAElFTkSuQmCC\n",
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"x = np.linspace(-1,5,100)\n",
"y1 = -(2./3.)*x + (7./3.)\n",
"y2 = (1./4.)*x - (3./4.)\n",
"\n",
"fig = plt.figure(figsize=(5, 5))\n",
"\n",
"ax1 = fig.add_subplot(111)\n",
"\n",
"ax1.set_xlabel(\"$x$\", fontsize=14)\n",
"ax1.set_ylabel(\"$y$\", fontsize=14)\n",
"ax1.grid(True)\n",
"\n",
"ax1.plot(x,y1,'b', label='$2x+3y=7$')\n",
"ax1.plot(x,y2,'r', label='$x-4y=3$')\n",
"ax1.plot(37./11,1./11, 'ko', label='Our solution')\n",
"\n",
"ax1.legend(loc='best', fontsize=14)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Success!!"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Matrix inverse\n",
"\n",
"For a square $n\\times n$ matrix $A$, its inverse is defined as the matrix $B$ with the same dimensions as $A$ that satisfies\n",
"\n",
"$$AB = BA = I$$\n",
"\n",
"where $I$ is the identity matrix of size $n\\times n$ (sometimes to explicitly indicate its dimension we would write $I_n$) - see a couple of cells below for the formal definition.\n",
"\n",
"That is if we pre- or post-multiply $A$ by this matrix we get the identity matrix. "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Such a matrix $B$, if its exists, is called the inverse of $A$ and denoted $A^{-1}$. Note that $(A^{-1})^{-1} = A$.\n",
"\n",
"The action of taking the inverse can be written $\\text{inv}(A)$, so for example $\\text{inv}(\\text{inv}(A)) = A$.\n",
"\n",
"Since $AB \\ne BA$ in general, then for non-square matrices it is certainly possible to find a matrix $B$ such that $AB=I$ but $BA\\ne I$. In this case we say that $B$ is a *right inverse*. Similar definition for *left inverse*. \n",
"\n",
"For square matrices a left inverse (it it exists) is also the right inverse and vice-versa."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Properties of the matrix inverse\n",
"\n",
"If a matrix has an inverse, then that inverse is unique.\n",
"\n",
"How could you prove a result like this:\n",
"\n",
"consider a matrix $A$ and suppose $B$ and $C$ are both inverses, i.e.\n",
"\n",
"$$AB = BA = I\\quad\\text{and}\\quad AC = CA = I$$\n",
"\n",
"then it must follow that\n",
"\n",
"$$B = BI = B(AC) = (BA)C = IC = C$$\n",
"\n",
"Hence the inverse is unique."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The following properties of inverses hold\n",
"\n",
"For invertible matrices $A$ and $B$ (of the same square size) and scalar $\\alpha\\ne 0$\n",
"\n",
"\n",
"1. $\\alpha A$ is invertible and $(\\alpha A)^{-1} = \\alpha^{-1} A^{-1}$\n",
"\n",
"\n",
"2. $AB$ is invertible and $(AB)^{-1} = B^{-1}A^{-1}$\n",
"\n",
"\n",
"3. $A^{-1}$ is invertible and $(A^{-1})^{-1}=A$\n",
"\n",
"\n",
"4. $A^{T}$ is invertible and $(A^{T})^{-1}=(A^{-1})^T$\n",
"\n",
"\n",
"As an exercise you can check these through some example using NumPy."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Using the inverse matrix to solve a linear system\n",
"\n",
"**IF** we have the inverse matrix, then to solve the matrix equation \n",
"\n",
"$$ A\\boldsymbol{x}=\\boldsymbol{b} $$\n",
"\n",
"we can simply multiply both sides by the inverse of the matrix $A$:\n",
"\n",
"\\begin{align}\n",
"A\\boldsymbol{x} & = \\boldsymbol{b}\\\\\n",
"\\implies A^{-1}A\\boldsymbol{x} & = A^{-1}\\boldsymbol{b}\\\\\n",
"\\implies I\\boldsymbol{x} & = A^{-1}\\boldsymbol{b}\\\\\n",
"\\implies \\boldsymbol{x} & = A^{-1}\\boldsymbol{b}\n",
"\\end{align}\n",
"\n",
"so we can find the solution $\\boldsymbol{x}$ by multiplying the inverse of $A$ with the RHS vector $\\boldsymbol{b}$."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### The inverse for our simple example\n",
"\n",
"Note that for our $2\\times 2$ case, when we went through the steps of substituting one equation into another to solve the system, this actually pretty much gave us enough information to write down the inverse matrix.\n",
"\n",
"Let's check that we do indeed have an inverse:"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"scrolled": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[[1. 0.]\n",
" [0. 1.]]\n",
"[[1. 0.]\n",
" [0. 1.]]\n"
]
}
],
"source": [
"# our matrix\n",
"A = np.array([[2, 3], [1, -4]])\n",
"\n",
"# I claim this is the inverse to that matrix\n",
"B = (-1./11)*np.array([[-4, -3], [-1, 2]])\n",
"\n",
"# let's check\n",
"print(A@B)\n",
"print(B@A)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Row operations and a reminder of Gaussian elimination\n",
"\n",
"What did we do when we used substitution to solve our problem\n",
"\n",
"$$\n",
"\\begin{align*}\n",
" 2x + 3y &= 7 \\\\[5pt]\n",
" x - 4y &= 3 \n",
"\\end{align*}\n",
" \\quad \\iff \\quad\n",
" \\begin{pmatrix}\n",
" 2 & 3 \\\\\n",
" 1 & -4 \n",
" \\end{pmatrix}\n",
" \\begin{pmatrix}\n",
" x \\\\\n",
" y \n",
" \\end{pmatrix}=\n",
" \\begin{pmatrix}\n",
" 7 \\\\\n",
" 3 \n",
" \\end{pmatrix} \n",
"$$\n",
"\n",
"?\n",
" \n",
"\n",
"We rearranged the second equation to yield $x=3+4y$ and substituted it into the first to obtain an equation that contained $y$'s only. \n",
"\n",
"We can achieve exactly the same thing if we subtract twice the second equation from the first: the $x$'s cancel, we get $3-2\\times (-4)$ (i.e. 11) $y$'s on the LHS, and on the RHS we get $7-2\\times 3$, i.e. 1. \n",
"\n",
"Hence, we have $11y=1$ and so $y=1/11$.\n",
"\n",
" \n",
"\n",
"We effectively performed what are called **row operations** on the matrix and the RHS vector. \n",
"\n",
"Recalling that each row of the matrix encodes an algebraic equation, we can:\n",
"\n",
"1. multiply each row by a non-zero scalar, \n",
"\n",
"2. add multiples of one row to another, \n",
"\n",
"3. swap rows, \n",
"\n",
"without changing the overall information that the rows (or the corresponding equations) are telling us (as long as we remember that we must apply the same operation to the RHS value, otherwise the equation would change!).\n",
"\n",
" \n",
"\n",
"For the example above, multiplying a row by a scalar or swapping rows would not change the two lines in the plot and so not change the solution. Adding a multiple of one row to another and replacing one row with this updated row **will** change one of the lines, **but** the intersection point won't change - let's convince ourselves of this. \n",
" \n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Row operations do not affect the corresponding solution\n",
"\n",
"Our simple example was the following\n",
"\n",
"$$\n",
"\\begin{align*}\n",
" 2x + 3y &= 7 \\\\[5pt]\n",
" x - 4y &= 3 \n",
"\\end{align*}\n",
" \\quad \\iff \\quad\n",
" \\begin{pmatrix}\n",
" 2 & 3 \\\\\n",
" 1 & -4 \n",
" \\end{pmatrix}\n",
" \\begin{pmatrix}\n",
" x \\\\\n",
" y \n",
" \\end{pmatrix}=\n",
" \\begin{pmatrix}\n",
" 7 \\\\\n",
" 3 \n",
" \\end{pmatrix} \n",
"$$\n",
"\n",
"We rearranged each to the following in order to easily draw two straight lines:\n",
"\n",
"\\begin{align*}\n",
" y &= -\\frac{2}{3}x + \\frac{7}{3} \\\\[5pt]\n",
" y &= \\frac{1}{4}x - \\frac{3}{4}\n",
"\\end{align*}\n",
"\n",
"This gives the first image below which is just a repeat of what we plotted above."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"scrolled": false
},
"outputs": [
{
"data": {
"text/plain": [
""
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAVcAAAFBCAYAAADDvuyeAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAA+K0lEQVR4nO3dfXzN9f/H8cdr2FzMRStGxFIKITXUKDbXV0XCV6HSV0soUt9KSir0rahUKvryw6J9kxS5KmzRtysUufqK5DLLF8Xmcrb374/3zHA2uzjnfD5ne91vt3Ozc85nn89rx/Y87/P+vD/vtxhjUEop5V1BTheglFKFkYarUkr5gIarUkr5gIarUkr5gIarUkr5gIarUkr5QHGnC/CXyy67zEREROR6+6NHj1KmTBnfFeRFgVKr1uldgVInBE6t+alzzZo1B4wxFS94whhTJG6RkZEmLxISEvK0vZMCpVat07sCpU5jAqfW/NQJrDYeMke7BZRSygc0XJVSygc0XJVSygc0XJVSygc0XJVSygeKzFAspXJy5MgR9u/fT2pqKuXLl2fz5s1Ol3RRgVInBE6t59dZokQJKlWqRLly5fK8Lw1XVeQdOXKEP/74g6pVq1KqVClSUlIoW7as02VdVHJyckDUCYFTa9Y6jTEcP36cvXv3AuQ5YF3VLSAiJUXkBxFZJyIbReR5D9uIiLwpIttE5GcRudGJWlXhsX//fqpWrUrp0qUREafLUS4hIpQuXZqqVauyf//+PH+/21quJ4GWxpgUESkBfC0ii4wx32XZpgNQK+N2E/Buxr9K5UtqaiqlSpVyugzlUqVKlSI1NTXP3+eqlmvGBQ8pGXdLZNzOXyqhCzAjY9vvgAoiUsWbdezdC/Pne3WXyuW0xaqyk9/fDbe1XBGRYsAa4GpgojHm+/M2qQrsznJ/T8Zj+zzsKxaIBQgPDycxMTFXNUyaVJP4+Gv566/f6Nt3Z55/Bn9LSUnJ9c/mJLfWWb58eZKTkzPvp6WlnXPfrQKlTgicWrOr88SJE3n+3XVduBpj0oCGIlIBmCsi9YwxG7Js4ultxONCYMaYycBkgEaNGpno6Ohc1XDLLXDwYBJTp15J1apX8vzz4OaGTWJiIrn92Zzk1jo3b958zsmWQDz54naBUmt2dZYsWZIbbrghT/tyXbieYYz5S0QSgfZA1nDdA1yR5X414HdvHrt4cXjyyf9SvXplXnwRTp6Ef/7T3QGrlHIXV/W5ikjFjBYrIlIKaA3897zN5gH3ZIwauBk4bIy5oEugoIoVg8mT4aGH4JVX4NFHQRfKVW7z0ksv0bhxY8qVK0fFihW57bbb2LBhw8W/0csmTpxIgwYNKFeuHOXKlSMqKooFCxb4vY6IiAhE5IJbp06d/F6L21quVYDpGf2uQcBHxpjPRWQAgDHmPWAh0BHYBhwD+vmqmKAgmDgRgoNhwgRITYW33rKPK+UGiYmJDBw4kMaNG2OMYeTIkbRu3ZpNmzYRFhZWoH3fd999REREMGrUqItuW61aNV5++WVq1apFeno606dPp2vXrqxZs4YGDRoUqI68WLVqFWlpaZn39+3bR2RkJD179vRbDZk8zUNYGG8Fmc81Pd2YJ54wBozp39+YtLQ87crnCvNcmf6wadOmc+4fOXLEoUryxlOdycnJJigoyMybN88YY8xHH31kgoODzY4dOzK3eeSRR0zNmjVNUlJSjvu/9957zXPPPZfv+i655BLz3nvvZdZatWpVM378+HO2+fnnn01ISIjZuHFjvo+Tk9GjR5vy5cubo0eP5mr77P7vz/8dyQqdzzX/RGyf67PPwr/+Bf36QZY3R6VcIzk5mfT0dC655BIAunfvTv369Rk9ejQA48aN48MPP2Tx4sWEh4f7pIa0tDTi4+NJSUmhadOmmY9HRUWxatWqc7YdOnQo/fv3p27duuc8PnbsWEJDQ3O8rVy5Msc6jDFMmTKFPn36ULp0ae/9gLnktm4B1xKBF16AEiVg5Eg4dQri4uzJL1X4DB0Ka9f695gNG8IbbxRsH0OGDKFhw4ZERUUBdozm2LFj6dSpE1dddRVjxoxh+fLl1KpVq8D1nm/9+vVERUVx4sQJQkNDmTt3LvXr1898PioqinfeeSfz/qeffspPP/3ERx99dMG+BgwYcNGP8lWrVs3x+S+//JLffvuN/v375/En8Q6Nhjx69lkICYEnn7R9sLNm2T5ZpZw2bNgwvv76a77++muKFSuW+Xjbtm1p3LgxzzzzDPPnz6dx48Yev3/s2LGMHTs28/7JkycREcaNG5f52KJFi7j11ls9fv+1117L2rVr+euvv5gzZw733nsviYmJ1KtXD4Cbb76Zxx57jEOHDlGmTBkef/xxRo4cyaWXXnrBvsLCwgrcZ/z+++/TuHFjGjZsWKD95JeGaz488YQN1Ecfhe7dYfZsG7iq8ChoC9LfHn30UeLj40lISKBmzZrnPLd8+XLWrVuHMSbHroDzW4tPPvkkVatW5ZFHHsl8LKfWYnBwMFdffTUAjRo1YtWqVbz++utMmTIFgMjISIKDg1m9ejU//fQTxYsXZ9CgQR73dX7Qe5JT0O/fv5/PPvuMiRMn5rgPX9JwzaehQ23ADhoEXbvCJ5+AXp6unDBkyBDi4+NJTEykdu3a5zy3bt06unXrxltvvcWCBQsYPnw4S5Ys8bif81uLZcuWJSwsLDMw8yo9PZ2TJ09m3g8JCeGGG25g/vz5TJ8+nVmzZlGiRAmP31vQboFp06YREhJCr1698lW7N2i4FsDAgTZgY2Phttvgs88gAFYPVoXIoEGDiIuL49NPP+WSSy4hKSkJgNDQUA4ePEjHjh0ZNmwY999/P02aNKFBgwY+uVLuqaeeolOnTlxxxRUkJycza9YsEhMTLxjrGhUVxYQJE2jTpg2dO3fOdn8F6RYwxvCvf/2LXr16OXpVmIZrAfXvbwO2Xz/o2BE+/xwC4Co/VUicOUHUqlWrcx5/+OGH+fLLL+ncuTMjR44EoF69evTo0YPhw4fz7bfferWOpKQk+vTpQ1JSEuXLl6dBgwYsWrSIdu3anbNdw4YNCQoK4rXXXvPq8bNKTExk69atfPDBBz47Rm5ouHrBPffYgO3TB9q3h4ULoXx5p6tSRYHJ42WD//73v3O97bRp07y+7cyZM3nwwQe57rrrcr3vvIqJicnz6+ILGq5e0quXHabVqxe0aQNLlkDGUEOlirT09HT279/PO++8w/r16/MU8IFMLyLwojvvtCe21q2DVq3gwAGnK1LKeStWrKBWrVpMmzaNOXPmZF7gUNhpy9XLzpzY6toVWraEpUuhUiWnq1LKOdHR0Rw+fDggphz0Jm25+kD79rBgAWzbBtHRsM/rc3YppdxOw9VHWrWCxYth1y5o0QL27HG6IqWUP2m4+lDz5vDFF5CUZAN2p/tXjFFKeYmGq481bWr7XQ8dsmH7669OV6SU8gcNVz9o0gSWLYOjR20L9pdfnK5IKeVrGq5+cuONkJBgpyps0QI2bXK6IqWUL2m4+lH9+nBmdd7oaFi/3slqlFK+pOHqZ3Xrwldf2ctlo6Phxx+drkgp5Qsarg645hobsKGhdsjWDz84XZFSyts0XB1y1VWwYoWdf6B1a/jmG6crUkp5k6vCVUSuEJEEEdksIhtFZIiHbaJF5LCIrM24jXSiVm+oUcMGbOXK0Latbc0q5WZjx45FRBg8eLAjx584cSINGjSgXLlylCtXjqioqAvmjHULV4UrcBp4zBhTB7gZGCQidT1st9IY0zDj9oJ/S/SuatVsqFavDh062DGxSrnRd999x/vvv0+DBg0cq6FatWq8/PLL/Pjjj6xevZqWLVvStWtXfv75Z8dqyo6rwtUYs88Y82PG18nAZiDnJR4LgSpV7CiCq6+Gzp1h0SKnK1KBYvbs2YSEhLAzy+V/Q4YM4aqrruKPP/7w2nEOHz5M7969mTJlisdZrapVq3bBBNjr16+nZMmSbPLiuMMuXbrQoUMHrr76aq655hrGjBlD2bJlvT75tze4KlyzEpEI4Abgew9PR4nIOhFZJCK+m3XXjypVsuNg69a1M2rNm+d0RSoQdO/enfr16zN69GgAxo0bx4cffsjixYtzXIwwr2JjY+nevTstW7b0+HxUVBSrVq0657GhQ4fSv39/6tY998Pn2LFjCQ0NzfG2cuXKi9aUlpZGfHw8KSkpNG3aNP8/nI+4cspBEQkF5gBDjTFHznv6R6CGMSZFRDoCnwIeF2EXkVggFiA8PJzEM4NMcyElJSVP23vL888X54knGtCtWyjPPruZFi3+d9HvcarWvHJrneXLlyc5OTnzflpaGqcGDiTIzwOR0+vX5+TLL+d6+7S0NFJSUnjmmWfo0aMH1apVY9y4ccyfP5/KlSuf8zMVxLRp0/jll1945513SE5Otq/PqVPn7P/GG2/k/fffz3zs888/58cff2TKlCmZ33Pmud69e9OxY8ccj3n55ZdnW//GjRtp3bo1J06cIDQ0lJkzZxIREeGVnzdrnVmdOHEiz7+74oblELISkRLA58ASY8xFF9oRkR1AI2NMjlNTN2rUyKxevTrXdfhiEbfcOnzYrsf1/fcwYwbcfXfO2ztZa164tc7NmzdTp06dzPvJycmUffZZWLvWv4U0bJinNb2Tk5Mz50ht2rQpP/zwA/Pnz6dDhw4et3/mmWcYM2ZMjvtMSEg45/9oy5Yt3HLLLaxcuTJzZdno6Gjq1avH22+/nbndN998Q7NmzTh48CBlypThuuuuY/DgwQwdOvSCWgvq1KlT7Nq1i7/++os5c+bw/vvvk5iYSL169Qq87+zqPP93JCsRWWOMaXT+465quYqIAFOAzdkFq4hUBv4wxhgRaYLt2jjoxzJ9rnx5O11h587Qty+kpsK99zpdVRGTh5Bz2vLly1m3bh3GmBy7AoYOHUqfPn1y3Ff16tXPuf/tt99y4MCBc4IrLS2NFStW8N5773H06FFCQkKIjIwkODiY1atX89NPP1G8eHEGDRrk8Rhjx45l7NixOdaxaNEibr31Vo/PBQcHZy733ahRI1atWsXrr7/OlClTctynv7kqXIFmQF9gvYiszXjsaaA6gDHmPaA78JCInAaOA72M25rfXlC2rF3osEsXu7JsaqpdaVaprNatW0e3bt146623WLBgAcOHD2fJkiUet73sssu47LLL8rT/rl270qjRuY2yfv36UatWLZ5++mmCg4MBCAkJ4YYbbmD+/PlMnz6dWbNmUaJECY/7HDBgAD179szxuFWr5v48dnp6OidPnsz19v7iqnA1xnwNyEW2eRt4O6dtCosyZWD+fOjWDR54wE76MnCg01Upt9i1axcdO3Zk2LBh3H///TRp0oQGDRp4tfulQoUKVKhQ4ZzHypQpQ1hY2AUfw6OiopgwYQJt2rShc+fO2e4zLCyMsLCwfNXz1FNP0alTJ6644gqSk5OZNWsWiYmJrhzr6qpwVRcqVQo+/RR69IBBg2zAZnRjqSLs0KFDdOvWjc6dOzNypL2Opl69evTo0YPhw4c7MjSpYcOGBAUFXTAky5uSkpLo06cPSUlJlC9fngYNGrBo0SLatWvns2Pml4ZrAAgJgY8/tie2Hn0UTp6EJ590uirlpLCwMFavXn3ByRd/LFud3VnzmTNn8uCDD3Lddb4bHTlt2jSf7dvbNFwDRHAwxMfDPffAU0/ZgB0ZsBf+qsIiPT2d//3vf0ybNo3169f7JdwDhYZrACleHOLibNA+95ztInjxRaerUkXZihUraNmyJddeey1z5szxePVWUaXhGmCKFYOpU23AjhljW7AXGY+tlM9ER0eTnp7udBmupOEagIKC4L33bMCOGwfbt19NdDRIjuMslFL+5Nq5BVTOgoLgrbfsCa5PPqnGwIGgDQil3ENbrgFMBMaPh6SkXbz3XnVOnYLJk23XgVLKWRquAU4EHnhgO7VqVeeFF+yVXFOn2pNfSinn6J9gISACzz9v+2CfecaOIoiLg2yuPlRK+YGGayEyYoS94OAf/7ABGx9vA1cp5X96QquQefxxmDAB5s6FO++EEyecrkipoknDtRB65BF49134/HM7q9bx405XpFTRo+FaSA0YAFOmwJdf2nlhjx51uiKlrFGjRnllYmsAEeHjjz/2yr68TcO1ELv/fruSQWKiXVnWS6t+KBfZu3cvsbGxVKtWjeDgYKpWrcoDDzzAnj17nC7Nq+677z6P0xju27eP2267zYGKLk7DtZDr0wdmzYJvvoF27ewSMqpw2LFjB40aNWLDhg1Mnz6dbdu28cEHH7Bx40YaN27Mjh07CrT/U6dOeadQH6pcuTIhISFOl+GRhmsR8Le/wezZsHo1tG4Nhw45XVHhdGahvKCgICIiIpg5c6ZPj/fYY48RFBTE0qVLadWqFdWrVycmJoalS5cSFBR0zjIr0dHRDB48+JzvP781GB0dzUMPPcTjjz9OxYoVadasmcfj7t69my5duhAWFkbp0qWpXbs28fHxmc+vX7+e1q1bU6pUKcLCwrjvvvs4nMO7uqdWadaug1GjRjF9+nQWLFiAiCAimdMent8tcLFjnznWhAkTqFq1Kpdccgn9+vXj2LFj2daXXxquRcQdd8Ann8DPP0OrVnAgx+UcVV7NnDmT2NhYdu7ciTGGnTt3Ehsb67OAPXToEEuXLmXQoEGULl36nOdKly7NwIEDWbRoEX/++Wee9vvBBx9gjGHlypXMmDHD4zYDBw7k2LFjJCQksHHjRt54443M1QqOHTtG+/btCQ0N5YcffmDu3Ll888032a6nlRuPP/44PXv2pHXr1uzbt499+/Z5XEo7u2Pff//952y3cuVKNmzYwNKlS/n3v//N3LlzmTBhQr7ry46Ocy1COneGefOga1eIiYGlS8GLS9sXaSNGjLig9XPs2DFGjBhB7969vX68rVu3YozJdkXSunXrYoxh69atNGnSJNf7vfLKKxk/fnyO2+zcuZM777yT66+/PvN7zpg5cyYpKSnExcVlTuQ9efJkYmJi2LZtW+bCgnkRGhpKqVKlCAkJoXLlytlul9tjlytXjnfffZfixYtTp04devTowbJlyxg+fHiea8uJtlyLmHbtYMEC2L4doqPh99+drqhw2LVrV54e9xbJZiq0M2t2Zvd8diIjIy+6zZAhQxg9ejRRUVE888wzrFmzJvO5zZs306BBg3NWSGjatClBQUFs2rQpT7XkVW6PXbduXYpnuT788ssvZ//+/V6vR8O1CGrZ0i7dvWcPtGgBu3c7XVHgO39J6os9XlC1atVCRNi4caPH5zdv3oyIcNVVVwEQFBTE+Yskp6amXvB9ZcqUueix//73v/Pbb7/Rr18/fvnlF5o2bcqoUaMAG+rZBXp2j+e2tovJ7bHPX5VWRHwyJ62GaxF1663wxRewfz80bw4FPLFc5I0ZM8Zj3+eYMWN8crywsDBatWrFO++847E7YuLEiXTo0CFzldWKFSuyb9++c7Zbt25dvo9frVo1YmNj+eijj3jhhReYPHkyYFuF69atIznLuL9vvvmG9PT0bLswPNW2du3ac+4HBweTlpaWY035ObYvuSpcReQKEUkQkc0islFEhnjYRkTkTRHZJiI/i8iNTtRaGERFwbJl8NdfNmC3bXO6osDVu3dvJk+eTI0aNRARatSoweTJk33S33rGuHHjOH36NK1bt2b58uXs3r2bxMRE2rRpgzGGt98+uwJ9y5YtWbRoEfPmzWPLli0MGzaM3fn8yDJkyBAWL17M9u3bWbt2LYsXL6Zu3bqAfR3KlCnDPffcw/r161mxYgUPPvggt99+e7b9rS1btuSnn35i6tSpbNu2jVdeeYX//Oc/52wTERHBhg0b2LJlCwcOHPDYss3u2N26dctXX29BuSpcgdPAY8aYOsDNwCARqXveNh2AWhm3WOBd/5ZYuDRqBAkJcOyY7SLYssXpigJX79692bFjB+np6ezYscOnwQpQs2ZNVq9ezXXXXUffvn2pWbMmd999N3Xq1GHVqlXnnGi6//77M2/NmjUjNDSUO+64I1/HTU9P5+GHH6Zu3bq0adOG8PBwpk+fDtjW+pIlSzhy5AhNmjShS5cuREVFMXHixGz3165dO5577jlGjBhBZGQkO3bsYODAgeds88ADD1CnTh0aNWpExYoVLwjfnI49derUfP2cBWaMce0N+Axoc95jk4C7stzfAlS52L4iIyNNXiQkJORpeyd5o9b1642pVMmY8HBjNmwoeE2euPU13bRp0zn3jxw54lAleRModRoTOLVmV+f5vyNZAauNh8xxW8s1k4hEADcA35/3VFUg6+eZPRmPqQKoVw+++souHxMdDQXojlNK4dJxriISCswBhhpjjpz/tIdvMR4eQ0RisV0HhIeHZ17VkRspKSl52t5J3qz1lVdKMWzY9TRvXoxXXlnHtdemeGW/4N7XtHz58uecBElLSzvnvlsFSp0QOLVmV+eJEyfy/rvrqTnr5A0oASwBhmXzvHYLnMfbtf76qzE1ahhTvrwx333nvf269TXVbgHfC5RaC223gNjBaFOAzcaY17LZbB5wT8aogZuBw8aYfdlsq/KhZk3bRXDppdCmDXg4d6CUughXhSvQDOgLtBSRtRm3jiIyQEQGZGyzENgObAPeBwZmsy9VADVqwIoVUKWKvarLhZ/mvcoYjz1LSuX7d8NVfa7GmK/x3KeadRsD5H8WCJVrVavaFmyrVtCxo52XoHVrp6vyvhIlSnD8+PELLgJQCuD48eMXXNWVG25ruSqXqVzZtlpr1bITvyxc6HRF3lepUiX27t3LsWPHtAWrMhljOHbsGHv37qVSpUp5/n5XtVyVO1WsCMuXQ9u2dkat2bPt2lyFRbly5QD4/fffSU1N5cSJE5QsWdLhqi4uUOqEwKn1/DpLlChBeHh45u9IXmi4qly59FJ7qWy7dtC9u13doEcPp6vynnLlymX+ASUmJnLDDTc4XNHFBUqdEDi1erNO7RZQuVahgl3w8KaboFcvG7BKKc80XFWelCtnpyts3tyuz5VxSblS6jwarirPQkPthNutW8N998H77ztdkVLuo+Gq8qV0aTs0q0MHiI2FHCY9UqpI0nBV+VayJMyda0cODB4Mr2V3TZ1SRZCGqyqQkBA7NKt7d3jsMXjpJacrUsoddCiWKrASJeDDDyE4GJ5+Gk6dgpEjIY9r4ylVqGi4Kq8oXhxmzLABO2oUnDwJY8ZowKqiS8NVeU2xYjBliu0qeOkl24J99VUNWFU0abgqrwoKgnfftS3Y8eNtwE6YoAGrih4NV+V1IjZQswbsO+84XZVS/qXhqnxCxHYJBAef7SLw8WKoSrmKhqvyGRF7UiskxJ7k2r27Di1a2JNfShV2Os5V+ZQIPPccjB0LS5eGc/fdkJrqdFVK+Z62IZRfDB8Ou3dv4913ryY1FeLjbYtWqcJKW67Kb3r23MObb8Knn0K3bnDihNMVKeU7Gq7Krx5+GCZNssvF3H47HDvmdEVK+YaGq/K72FiYOhWWLoVOnSAlxemKlPI+DVfliH79IC7OLt/doQMcOeJ0RUp5l+vCVUSmish+EdmQzfPRInJYRNZm3Eb6u0blHb172xNb335rFz/86y+nK1LKe1wXrsA0oP1FtllpjGmYcXvBDzUpH+nRAz7+GH780a5scOiQ0xUp5R2uC1djzApA/8SKkK5d7aTbGzZAy5bwv/85XZFSBee6cM2lKBFZJyKLROQ6p4tRBdepk102ZssWiImBpCSnK1KqYMQY43QNFxCRCOBzY0w9D8+VA9KNMSki0hGYYIyplc1+YoFYgPDw8Mj4+Phc15CSkkJoaGh+yve7QKk1N3WuXVuB4cPrU7HiScaPX0vFiqf8VN1Zhen1dItAqTU/dcbExKwxxjS64AljjOtuQASwIZfb7gAuu9h2kZGRJi8SEhLytL2TAqXW3Na5cqUxoaHGXHWVMTt3+rYmTwrb6+kGgVJrfuoEVhsPmRNw3QIiUlnEzg4qIk2wXRsHna1KedMtt8CXX8KBA9CiBfz2m9MVKZV3rgtXEfkQ+Ba4VkT2iMjfRWSAiAzI2KQ7sEFE1gFvAr0y3j1UIXLzzbBsGRw+DM2bw7ZtTlekVN64buIWY8xdF3n+beBtP5WjHBQZCQkJdohW8+awfDnUru10VUrljutarkpldf31kJgI6em2i2CDx0tLlHIfDVfletddZwO2WDE7TGvtWqcrUuriNFxVQKhd285DUKqUvdBg9WqnK1IqZxquKmBcfTV89RWULw+tWsF33zldkVLZ03BVAeXKK23AVqwIbdrAypVOV6SUZxquKuBUr267CKpWhfbt7YgCpdxGw1UFpMsvty3YK6+Ejh3hiy+crkipc2m4qoAVHm5brddeC7fdBgsWOF2RUmdpuKqAVrGivbigfn244w47daFSbqDhqgJeWJhdjysy0k6+/dFHTleklIarKiQqVLD9rlFRcNddMHOm0xWpok7DVRUaZcvC4sX2Mtm+feH//s/pilRRpuGqCpUyZeDzz+0Y2Pvvh0mTnK5IFVUarqrQKV0aPvvMLh0zYAC89ZbTFamiSMNVFUolS8Inn9jFDx95BMaPd7oiVdRouKpCKzjYjhzo2RMefxzGjnW6IlWUuG6ybKW8qUQJO3IgOBhGjICTJ2HUKLALBSnlOxquqtArXhymTbNB+8ILcOqUbcVqwCpf0nBVRUKxYvCvf9kW7D//aVuw48drwCrf0XBVRUZQELz7rg3Y11+3Ldg337SPK+VtGq6qSBGBCRMgJATGjbMB+957GrDK+zRcVZEjAq+8YluwY8fagJ0yxXYdKOUtF32/FhG/LmYsIlNFZL+IeFznU6w3RWSbiPwsIjf6sz5VOIjA6NF25MD06XDPPXD6tNNVqcIkNx+GfhKRCSJyic+rsaYB7XN4vgNQK+MWC7zrh5pUISQCzz0HL70Es2bZCV9SU52uShUWuQnXJsB1wFYReVhEfPrhyRizAjiUwyZdgBnG+g6oICJVfFmTKtyeegpeew0+/hi6d4dTp3QIgSq4i/a5GmPWA61FpCvwKvCQiDxmjFnk6+KyURXYneX+nozH9jlTjioMHn3U9sEOHgx//FGP5s3tJbRK5ZcYY3K/sUgw8CjwNPAfYJgx5r9eL0okAvjcGFPPw3MLgJeMMV9n3F8GPGGMWeNh21hs1wHh4eGR8fHxua4hJSWF0NDQ/P0AfhYotQZCnZ9/XoXXXruGG2/8k9GjN1CyZLrTJWUrEF7PMwKl1vzUGRMTs8YY0+iCJ4wxub4BFYDWwJtAGnAq4+vyedlPLo4TAWzI5rlJwF1Z7m8Bqlxsn5GRkSYvEhIS8rS9kwKl1kCp88knN5mgIGOio41JTna6muwFyutpTODUmp86gdXGQ+bkZrTAUBGZKSK/AAeB+UBjYALQH7gW2CQiN+Up7vNvHnBPxqiBm4HDxhjtElBe0779H8TFwcqVdunuI0ecrkgFotyMc30M+BZ7Vv47YI0x5lSW52eIyJPAVOyJrwIRkQ+BaOAyEdkDPAeUADDGvAcsBDoC24BjQL+CHlOp8919t+2DvesuO/H24sVwib/Gy6hCITcntK7IxX7+D/DKhG7GmLsu8rwBBnnjWErlpHt3G7Ddu0Pr1naNrksvdboqFSi8ddHf/4CWXtqXUq5x++12VYONG6FlS9i/3+mKVKDwSrhm9Ot+5Y19KeU2HTrYdbm2boWYGNinPfwqF3S6CqVyoXVrWLQIdu6E6GjYu9fpipTbabgqlUstWsCSJbbl2ry5DVqlsqPhqlQeNGsGS5fCwYM2bLdvd7oi5VYarkrlUZMmsHw5JCfbgN261emKlBtpuCqVDzfeCAkJdrmY5s1h82anK1Juo+GqVD41aACJiWCMbcGuX+90RcpNNFyVKoC6deGrr+zKsjEx8NNPTlek3ELDVakCuvZaWLECypSxFxqsWuV0RcoNNFyV8oKrrrIt2EsusWNiv/3W6YqU0zRclfKSiAgbsJUqQdu2tjWrii4NV6W86IorbKhWq2Yvm122zOmKlFM0XJXysipV7CiCmjWhc2c7XaEqejRclfKB8HA7DrZ2bejSBebPd7oi5W8arkr5yGWX2W6BBg2gWzf45BOnK1L+pOGqlA+Fhdm5CBo3hp494d//droi5S8arkr5WPnydjatZs3s8jFxcU5XpPxBw1UpPyhbFhYutHPB3nsvTJnidEXK1zRclfKTMmXsigZt20L//vDuu05XpHxJw1UpPypVCj791A7RGjgQJkxwuiLlKxquSvlZyZIwZ44dQTB0KLz6qtMVKV9wXbiKSHsR2SIi20TkKQ/PR4vIYRFZm3Eb6USdShVEcDDEx0OvXvDEEzB6tNMVKW8r7nQBWYlIMWAi0AbYA6wSkXnGmE3nbbrSGNPZ7wUq5UUlSsAHH9h/n30WTp2C558HEacrU97gqnAFmgDbjDHbAUQkHugCnB+uShUKxYrB//2fbcm++KJd2eCf/9SAdYQxFDt+3Gu7c1u4VgV2Z7m/B7jJw3ZRIrIO+B143Biz0R/FKeULxYrB5Mk2YF95xbZgX3tNA9Zvtm2zg4/j4qhZv76dcccL3Baunn6dzHn3fwRqGGNSRKQj8ClQy+PORGKBWIDw8HASExNzXUhKSkqetndSoNSqdeasRw/Yv/9q3nijGr/9tpdHHtlKUA5nRQLl9QT31Vr8yBEqJSQQ/uWXlN+4ESPCn5GR7Ktbl63eqtMY45obEAUsyXJ/ODD8It+zA7jsYvuOjIw0eZGQkJCn7Z0UKLVqnReXnm7MP/5hDBjTv78xaWnZbxsor6cxLqn15Elj5s415o47jClRwr7I9eoZ8/LLxuzZY4zJX53AauMhc9zWcl0F1BKRK4G9QC/g7qwbiEhl4A9jjBGRJtgRDwf9XqlSPiACL78MISF2BMGpUzB1qu06UPlgDHz/PcyYYSd2OHTITlk2eDD07QsNG/qs/8VV4WqMOS0ig4ElQDFgqjFmo4gMyHj+PaA78JCInAaOA70y3j2UKhRE7Mmt4GAYORJSU202FHfVX6vL/fabHYoRFwdbt9rBxXfcYQO1TRu/vJiu++8yxiwEFp732HtZvn4beNvfdSnlb88+awP2qadsC3bWLHtfZePwYZg92wbqmTV2oqPtC9i9O5Qr59dyXBeuSqmznnzSdhE8+qjNh9mz7X2VITXVTjkWFweffWbHstWuDWPGQO/eUKOGY6VpuCrlckOH2hbroEHQtauddLtUKaercpAx8OOPNlBnzYL//c/OTB4baz/2N2rkinFsGq5KBYCBA23AxsbCbbfZRlqRs3s3zJxpO6A3b7YvyO2320Bt3951fSYarkoFiP79bX706wcdO8KTTxaBIQTJybapPmOGXZTMGDvr+KRJdmDwJZc4XWG2NFyVCiD33GMDtk8feOKJBjRrZlc6KFROn7aLj8XF2WA9ftwupfvcc/YHv+oqpyvMFQ1XpQJMr142YHv2LEubNvZ8josbcLm3bp0N1JkzISkJKlSwyzb07QtRUa7oR80L1005qJS6uG7d4IUXNrJuHbRqBQcOOF1RPv3+O4wbB9dfbwf0v/km3HSTnfA2Kcku19C0acAFK2i4KhWwmjY9yGef2XM7MTGwf7/TFeXS0aO2ddquHVxxBfzjH3b4w9tv27D99FP77hHgY860W0CpANa+vV2X6/bb7Xj5ZcugShWnq/IgPR2WL7cf+z/+GFJS7BjU4cNtR/I11zhdoddpuCoV4Fq1gkWL7AiCFi1swF5xhdNVZdi0CeLiuHnqVNu0LlcOeva0gXrrreQ47VeA03BVqhBo3hy++MJORdqihW0kRkQ4VMz+/fDhh7aVumYNFCvG0caNKfnmm7aJXUSugCi8bxtKFTFNm8LSpfDnnzZgf/3Vjwc/ftzOOtW5M1x+ub2szBh4/XXYu5f1L70Ef/tbkQlW0HBVqlBp3Ni2Wo8eta3ZLVt8eLD0dDtBygMPQOXKdozY2rXw2GOwfr1ttQ4daqf4K4K0W0CpQuaGG+zFTK1bnz3JVbeuFw+wdWvmsijs2AFlysCdd9rxqDExOvlsBm25KlUI1a8PZ1YriY6Gn38u4A4PHoR33rGD+a+5xs46dc01NmD/+AOmT7dprsGaScNVqUKqTh37qT042DYof/wxjzs4eRLmzrWTTFepYqflOnoUXn3VTqKyZIm9HLVMGZ/UH+i0W0CpQqxWLRuwMTF2yNaSJdCkSQ7fYAx8993ZZVH+/NP2mT788NllUVSuaMtVqUKuZk0bsGFh9pP7qFEziYiIICgoiIiICGbOnAnbt8MLL9iP+k2b2o/57dvDwoWwZw+MH6/BmkfaclWqCKhRA776Cho3nsnzz8cCxwDYuXMnsffcA+np9AbbQfv00/YElZ+XRSlsNFyVKiKqVYNixZ7mTLCecSw9nREVKtB73TqoXt2Z4gohDVelCjtj7JjTGTP4fe8uj5vsOnxYg9XLNFyVKqx27bKzT8XF2amzQkKoXro0O48du2DT6hqsXue6E1oi0l5EtojINhF5ysPzIiJvZjz/s4jc6ESdSrlScjJMmwYtW9rJBZ5+Gi69FCZPhqQkxkyeTOnSpc/7ptJ06TLGgWILN1e1XEWkGDARaAPsAVaJyDxjzKYsm3UAamXcbgLezfhXqaLp9Gk7qUBcnB2Xevy4XQrluefs8KmaNTM37d27NwAjRoxg165dVKtWnZIlxzBxYm9uvhnuusupH6LwcVW4Ak2AbcaY7QAiEg90AbKGaxdghjHGAN+JSAURqWKM2ef/cpVyTplt22D+fLu8dFKSXevl3nvtdH4335zt7P29e/fODFmwU6t27myvBzh1yu5CFZzbwrUqsDvL/T1c2Cr1tE1VQMNVFX6//57Zj9p4/XooUQI6dbIt1E6d8jV7f2ioHc7apYtdWTY11a40qwrGbeHq6a3W5GMbu6FILBALEB4eTuKZi61zISUlJU/bOylQatU68yfo+HEqfv014V98wSU//oikp3OkTh12DhjAkfbtST2z/Ou33xboOI8/HsSRI9fxwAOXsmHDL3Tt+rsXqrfc9ppmx6t1GmNccwOigCVZ7g8Hhp+3zSTgriz3twBVLrbvyMhIkxcJCQl52t5JgVKr1pkHp08bs3SpMffea0yZMsaAMRERxjzzjDFbthhjfFPniRPG3H67Pdxrr3lvv654TXMhP3UCq42HzHFby3UVUEtErgT2Ar2Au8/bZh4wOKM/9ibgsNH+VlVYbNx4dnnpPXvsVVJ33WU/9t9yi8+XRQkJgdmz4e67Ydgw2wf75JM+PWSh5apwNcacFpHBwBKgGDDVGLNRRAZkPP8esBDoCGzDXmrSz6l6lfKKM8uizJhhp64qVsyu1zJ+PNx2m99n7w8Ohvh4e17sqafs5FgjR/q1hELBVeEKYIxZiA3QrI+9l+VrAwzyd11KedXx4zBvng3UJUsgLQ0iI+GNN2xLtVIlR8srXtw2oIOD7YiuU6fgxRezHYCgPHBduCpVaKWnw9df20CdPRuOHLEX/D/+uP3Yf911Tld4jmLFYOpUOyBhzBjbgn3lFQ3Y3NJwVcrXfvnl7LIoO3fayaW7d7efu1u0cPXs/UFBMGmSbcGOG2dbsG+8oQGbGxquSvnCgQN2sukZM+CHH2xKtW5tm4BduwbU7P1BQfD22zZg33jDjoN9+22fn1sLeBquSnnLyZOwYIEN1IULbQrVr2+bfHffbZdKCVAi8NprdjTByy/bFuykSa5udDtOw1WpgjDGDt6Pizu7LErlyvDII7Yf9frrna7Qa0TgpZdsC/bFF23ATp1qT36pC+nLolR+bN8OH3xgQ3XbNjtc6o47bD9qq1aFNnFE7GowwcHw7LM2YOPi7Ekvda7C+RuglC/8+ac9yx8XZ8/6i9hlUUaMsMuilC3rdIV+88wztovgiSds78eHH9rAVWdpuCqVk9RUWLzY9qPOm2ebanXqwNixdhqpK65wukLH/OMfNmCHDLGDH2bPzte8MYWWhqtS5zMGVq+2gfrhh/bM/2WXwYMP2vn4brxRxyJleOQR22J96CE7q9bcuX6/oMy1NFyVOmPXLvjgAxpPmmS/DgmB22+3J6bat9eOxWwMGGADtn9/Oy/svHkBNdLMZzRcVdF25AjMmWNbqRlTzaU2aGA7FXv0gAoVHC0vUNx/v33vue8+Oy3CggVFqgvaIw1XVfScPg1ffmlPTH36qb3Ov1Ytexq8d2/W7tpFdHS001UGnL59bQu2d29o2xYWLSra700arqpoMAbWrTs7nd8ff0BYmJ16v29fuOmms/2ouzwvP60u7m9/sy3YXr2gTRs7J01YmNNVOUPDVRVuWZZF4cyyKJ0720Dt2FFPb/tAt27wySd2dFqrVvZDQlGkVwerwufoUTvAv21bO1TqiSfsGZaJE2HfPvuXf8cdGqw+dObE1n//a4cCHzpU9E4GastVFQ5paZCQYFuoc+bYgI2IsAP8+/a1farKr9q1sye2brsNHn20ITffDJdf7nRV/qPhqgLbhg02UD/4wHYBlC9vJ0np2xeaNdOpmxzWsqW9BqNduxBatIDly4vOdRf6m6cCzx9/2LnvbrzRzjr12mv2648+gqQkmDwZbr1Vg9Ulbr0Vxo37mf377fS1O3Y4XZF/6G+fCgzHj9uFnTp1gqpV4dFHbXhOmAB798L8+XZcasmSTleqPKhb9wjLlsFff0Hz5naum8JOuwWUe6Wnw8qVdoD/xx+fXRbliSfsx/46dZyuUOVBo0a2W6BNGzK7CK691umqfEfDVbnPli1n+1F37oTQUDuu59577V+lftwPWA0b2vOOrVrZ/8ply1y3dJjX6G+pcocDB+zaITfdBLVr21mZa9c+O+B/2jSIidFgLQTq1YOvvrL/ldHR9tqOwsg1v6kiEiYiX4rI1ox/L8lmux0isl5E1orIan/Xqbzo5Ek7bKprV7sEysMP28fGjYM9e+xp5rvvhtKlna5UeVnt2rBihZ1BKyYG1qxxuiLvc024Ak8By4wxtYBlGfezE2OMaWiMaeSf0pTXGAP/+Y+dSqlKFTsR6A8/2ElB1661t8ceC+j1plTuXH21bcGWK2e7Cb7/3umKvMtN4doFmJ7x9XSgq3OlKK/79VdqTJ9uB/PfcovtU+3Y0bZOd++2rdVCtN6Uyp0rr7Qt2Msusye6vv7a6Yq8x03hGm6M2QeQ8W+lbLYzwBciskZEYv1Wncq7P/+0S4TecgtcfTUR06fbq6amTbPjUT/4wF7Go0uIFmnVq9sWbJUqdtrchASnK/IOMcb472AiS4HKHp4aAUw3xlTIsu2fxpgL+l1F5HJjzO8iUgn4EnjYGLMim+PFArEA4eHhkfHx8bmuNSUlhdDQ0Fxv7yQ31SqpqYT98AOVv/iCS7/9lqDUVI7WqEFSu3b8dvPNFL/ySqdLvCg3vZ45CZQ6IXe1HjoUzLBh15OUVJLRozfQqNGffqrurPy8pjExMWs8dlEaY1xxA7YAVTK+rgJsycX3jAIez83+IyMjTV4kJCTkaXsnOV5rerox339vzODBxlx6qTFgTMWKxgwZYszq1fZ544I6c0nr9L7c1rp/vzENGhgTEmLMggW+rcmT/LymwGrjIXPc1C0wD7g34+t7gc/O30BEyohI2TNfA22BDX6rUJ1r504YM8YO5r/pJnj/fWjdGj7/3F419cYbEBmp602pXKtY0V5ccN11dhDJZxekQOBwU7j+E2gjIluBNhn3EZHLRWRhxjbhwNcisg74AVhgjFnsSLVF1ZEjMHWqHT8TEWGXQ6lUyV7Pn5R09hJVXW9K5dOll9qLC2688eyqsoHINVdoGWMOAq08PP470DHj6+2AnlL2tzPLosyYYZdFOXHCnvV/8UW7pkcA9KOqwFKhAnzxhX2f7tXLrnB+991OV5U3rglX5TLG2DGncXEwa9bZZVH+/nd7XX+TJvpxX/lUuXJ2Ha7bboM+feDUKbsAYqDQcFXn2rv37LIoGzacuyxKp052BTql/CQ01E643bWrXe7s1CmIDZABmBquClJSYO5c+7F/2TLbao2KgnffhZ49i+4Kc8oVSpe2S8bceSc8+KAN2MGDna7q4jRci6q0NHtaNi7Oril19KjtO332WfsZTJdFUS5SsqT9Nf3b3+wUFKdOwbBhTleVMw3XombDBttCnTnz7LIod91lp/Nr1kz7UZVrhYTYkQO9e9vpJ06ehOHDna4qexquRUFSEnz4oQ3VtWuheHF7neEbb9izBTp7vwoQJUrY86slSsDTT9sW7MiR7mwTaLgWVseO2RHYcXF2TEtamp0KfsIE21KtWNHpCpXKl+LFbTshOBhGjbIt2DFj3BewGq6FSXq6nWLozLIoycl2qc0nnoB77rGTaCpVCBQrBlOm2IB96SXbgn31VXcFrIZrIVB61y4YMcLOMrVrlx2/0qOHHT6ly6KoQiooCN57zwbs+PE2YCdMcE/AargGqgMH7KWmM2bQZNUq+5vWtq19G+/aVWfvV0WCCLz5pj3ZNX687SJ49113tCc0XAPJyZN2UpQZM2DhQntZasOGbHvoIa5+9lmdvV8VSSK2SyAkBMaOtS3Yf/3L+WmCNVzd7syyKHFx8NFHduH3KlVg6FD7sb9BA/YkJnK1BqsqwkRg9GgbsM89Z+cimDbNnvxyioarW23bZvtQ4+Jg+3b7Mb9bNxuorVo5/7aslMuI2GFZWYdpzZzp3ARtGq5ucuiQbZ3OmAHffmt/W1q2tG/Fd9wBZcs6XaFSrjd8uG3BPvaYbcHGx9v7/qbh6rRTp+zUPzNm2P7UU6egbl345z/tpSjVqjldoVIBZ9gwO4rg4YftnAQff+z/a2U0XJ1gjF1OOi7Ovq0ePGgH9Q8caMejNmzonvEkSgWowYNtwD74IHTpYucm8ucgGg1Xf9q50/ajzpgBv/xi30q7dLGB2qaNzt6vlJfFxto/q7//3c6cOX8+lCnjn2NruPra4cP2M0lcnF0/GKB5c3jySft5pXx5Z+tTqpDr188G7L332ik1Fi70z+kLDVdfOH3aXs8fF3d2WZRrrrHLovTpY9eeUkr5TZ8+tovg7rvttTaLFtmlZHxJw9VbjIGffjq7LMr+/XaS6fvvt2+ZjRtrP6pSDurZ0wZsz552keIvvvDtPPAargW1Z8/ZZVE2brT/e2eWRenYUZdFUcpFuna1J7buvNOOcvzyS99NEKfhmh8pKXZa9Bkz7Gz+xkDTprosilIBoFMnu2xMly52hfilS6FyZe8fxwXTG1gi0kNENopIuog0ymG79iKyRUS2ichTfiswLc1+jujbF8LD7Uf97dvtJSFbt9pLVAcM0GBVKgC0bWsXPvztN4iOtutyepubWq4bgG7ApOw2EJFiwESgDbAHWCUi84wxm3xW1fr19iP/mWVRKlSwveN9++qyKEoFsJYtYfFi23vXooX9EOpNrmm5GmM2G2O2XGSzJsA2Y8x2Y8wpIB7o4vVikpKoNns23HADNGgAr79uZ/GfPRv27YNJk+CWWzRYlQpwt95q+10PHLABu2+f9y7jEmOM13bmDSKSCDxujFnt4bnuQHtjTP+M+32Bm4wxHhfaFZFYIBYgPDw8Mj4+Plc1XPn++9SYNYsjtWvzR9u27I+JIdXX4zYKICUlhdDQUKfLuCit07sCpU5wf61btoTyj39cz7XXHuLVVzfn6XtjYmLWGGMu6Mr0a7eAiCwFPHUdjzDGfJabXXh4LNt3B2PMZGAyQKNGjUx0dHRuyoRrruGHdu1ocs89lAPcvsh0YmIiuf7ZHKR1eleg1AnurzU62p6T/vXXbV6r06/haoxpXcBd7AGuyHK/GvB7Afd5ocsv51j16l7frVLKverXh4MHU722P9f0uebSKqCWiFwpIsFAL2CewzUppdQFXBOuInKHiOwBooAFIrIk4/HLRWQhgDHmNDAYWAJsBj4yxmx0qmallMqOa4ZiGWPmAnM9PP470DHL/YXAQj+WppRSeeaalqtSShUmGq5KKeUDGq5KKeUDGq5KKeUDGq5KKeUDGq5KKeUDGq5KKeUDrpu4xVdE5H/Azjx8y2XAAR+V422BUqvW6V2BUicETq35qbOGMeaC9QyKTLjmlYis9jTTjRsFSq1ap3cFSp0QOLV6s07tFlBKKR/QcFVKKR/QcM3eZKcLyINAqVXr9K5AqRMCp1av1al9rkop5QPaclVKKR/QcM1Bbpf7dopjy4znkYhMFZH9IrLB6VpyIiJXiEiCiGzO+H8f4nRNnohISRH5QUTWZdT5vNM15UREionITyLyudO15EREdojIehFZKyIXrOGXVxquOTuz3PcKpws5X5ZlxjsAdYG7RKSus1VlaxrQ3ukicuE08Jgxpg5wMzDIpa/pSaClMeZ6oCHQXkRudrakHA3BTm4fCGKMMQ29MRxLwzUHuVzu2yn+WWbcC4wxK4BDTtdxMcaYfcaYHzO+TsYGQlVnq7qQsVIy7pbIuLny5ImIVAM6Af9yuhZ/03ANXFWB3Vnu78GFQRCoRCQCuAH43uFSPMr4qL0W2A98aYxxZZ3AG8ATQLrDdeSGAb4QkTUiElvQnblmmReneGG5b6fkaZlxlXsiEgrMAYYaY444XY8nxpg0oKGIVADmikg9Y4yr+rRFpDOw3xizRkSiHS4nN5oZY34XkUrAlyLy34xPXflS5MPVC8t9O8U/y4wXMSJSAhusM40xnzhdz8UYY/4SkURsn7arwhVoBtwuIh2BkkA5EfnAGNPH4bo8ylivD2PMfhGZi+16y3e4ardA4NJlxr1MRASYAmw2xrzmdD3ZEZGKGS1WRKQU0Br4r6NFeWCMGW6MqWaMicD+fi53a7CKSBkRKXvma6AtBXyz0nDNQXbLfbtBIC0zLiIfAt8C14rIHhH5u9M1ZaMZ0BdomTEcZ21Gq8ttqgAJIvIz9k32S2OMq4c5BYBw4GsRWQf8ACwwxiwuyA71Ci2llPIBbbkqpZQPaLgqpZQPaLgqpZQPaLgqpZQPaLgqpZQPaLgqpZQPaLgqpZQPaLgqpZQPaLgqRebE6CdFpEaWxyaIyK8iEu5kbSow6RVaSpE5r8Aq4CdjzAMi8jh2qrxmxpitzlanAlGRnxVLKbATUIvI09g5JH4FRmBn+9dgVfmiLVelshCRb7BTzd1mjFnkdD0qcGmfq1IZRKQlcD12IvI/HC5HBThtuSoFiMj1wFfAMOyaT6HGmHbOVqUCmYarKvIyRgh8A0wyxrwgIvWAn7F9romOFqcCloarKtJEJAz4D7DCGPNglsf/DVQ3xkQ5VpwKaBquSinlA3pCSymlfEDDVSmlfEDDVSmlfEDDVSmlfEDDVSmlfEDDVSmlfEDDVSmlfEDDVSmlfEDDVSmlfOD/AYdfKQe4kGYAAAAAAElFTkSuQmCC\n",
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"x = np.linspace(-1,5,100)\n",
"\n",
"y1 = -(2./3.)*x + (7./3.)\n",
"y2 = (1./4.)*x - (3./4.)\n",
"\n",
"fig = plt.figure(figsize=(5, 5))\n",
"\n",
"ax1 = fig.add_subplot(111)\n",
"\n",
"ax1.set_xlabel(\"$x$\", fontsize=14)\n",
"ax1.set_ylabel(\"$y$\", fontsize=14)\n",
"ax1.grid(True)\n",
"\n",
"ax1.plot(x,y1,'b', label='$2x+3y=7$')\n",
"ax1.plot(x,y2,'r', label='$x-4y=3$')\n",
"ax1.plot(37./11,1./11, 'ko', label='Our solution')\n",
"\n",
"ax1.legend(loc='best', fontsize=14)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"But now as an example let's add 3 times the second equation to the first, multiply the second by -2, and swap:\n",
"\n",
"$$\n",
"\\begin{align*}\n",
" -2x + 8y &= -6 \\\\[5pt] \n",
"5x - 9y &= 16\n",
"\\end{align*}\n",
" \\quad \\iff \\quad\n",
" \\begin{pmatrix}\n",
" -2 & 8 \\\\\n",
" 5 & -9 \n",
" \\end{pmatrix}\n",
" \\begin{pmatrix}\n",
" x \\\\\n",
" y \n",
" \\end{pmatrix}=\n",
" \\begin{pmatrix}\n",
" -6 \\\\\n",
" 16 \n",
" \\end{pmatrix} \n",
"$$\n",
"\n",
"rearrange each individually into the simple equation of a straight line:\n",
"\n",
"\\begin{align*}\n",
" y &= \\frac{2}{8}x - \\frac{6}{8} \\\\[5pt]\n",
" y &= \\frac{5}{9}x - \\frac{16}{9}\n",
"\\end{align*}\n",
"\n",
"This will be our second plot below."
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"scrolled": false
},
"outputs": [
{
"data": {
"text/plain": [
""
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAVcAAAFBCAYAAADDvuyeAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAA8RklEQVR4nO3deZyNdf/H8deX7GONRraUZElSg0zd3dmzRsKtkKUaZYlb7iJC3ZkWKgqJElkjKbIvM9EtWbKTTPatZIkx1pnP74/v5Dc4M87MnHOu68x8no/HeTTnnGuu621MH9/zvb6LERGUUkr5VhanAyilVEakxVUppfxAi6tSSvmBFlellPIDLa5KKeUHWlyVUsoPbnI6QKAULlxYSpcu7fXxZ8+eJU+ePP4L5EPBklVz+law5ITgyZqWnOvXr/9TRIpc94aIZIpHWFiYpEZUVFSqjndSsGTVnL4VLDlFgidrWnIC68RDzdFuAaWU8gMtrkop5QdaXJVSyg+0uCqllB9ocVVKKT/INEOxUpKQkMDBgwc5e/bsldfy58/Pjh07HEzlvWDJGmw5s2XLxi233EK+fPmcjqSCkBZX4M8//8QYQ7ly5ciSxTbmz5w5Q968eR1O5p1gyRpMOUNCQjh37hyHDh0C0AKrUs113QLGmPHGmD+MMVuTed8YYz40xsQYYzYbY+5P7zVPnTpFaGjolcKqlDGG3LlzU7x4cf744w+n46gg5MZqMgFokML7DYGyiY8I4OP0XjA+Pp5s2bKl9zQqA8qVKxeXLl1yOoYKQq4rriKyAjiRwiHNgC8SJ0esBgoYY25N73WNMek9hcqA9PciE/n6a/LExPjsdMHY51ocOJDk+cHE145ce6AxJgLbuiU0NJTo6GiPJ8yfPz9nzpy56rX4+PjrXnOrYMkarDnPnz+f7O+Ok2JjY12ZyxO3Z71l6VIqvPUWJatXJ/rOO31yzmAsrp6aEh43AhORscBYgKpVq0rNmjU9nnDHjh3X3WgJlpsvEDxZgzVnzpw5ue+++xxM5Fl0dDTJ/U67jauzTpwIkZHwyCPsevlln+V0XbeAFw4CJZM8LwEcdiiLUiqYjRsHnTpB3bowbx7xuXL57NTBWFznAE8njhqoAfwlItd1CWQGb731FtWqVaN48eIUKVKEpk2bsnWrx0EWfjVq1CgqV65Mvnz5yJcvH+Hh4cybNy/gOUqXLo0x5rpH48aNA55FBYFRoyAiAho0gDlzIHdun57edcXVGDMN+BEoZ4w5aIx5xhjzvDHm+cRD5gO7gRhgHNDVoaiOi46OpmvXrixZsoTly5dz0003UbduXU6cSOl+oHc6duzI4MGDvTq2RIkSvPPOO/z888+sW7eO2rVr07x5czZv3pzuHKmxdu1ajhw5cuXx888/Y4yhdevWAc2hgsD770P37tCsGcyeDTlz+v4antYhzIiPlNZz3b59+3WvnT59Otnj3ebvrGfOnJEsWbLInDlzRERkxowZkj17dtm7d++VY1988UW544475OjRoymes0OHDjJo0KA0ZypYsKCMGTPmyvPixYtLZGTkVcds3rxZcuTIIdu2bUvzdVLy5ptvSv78+eXs2bOp+r5r/+49/X64QbCskSrisqxvvSUCIi1bily8eNVbup6r8ujMmTMkJCRQsGBBAFq2bMk999zDm2++CcCwYcOYNm0aCxcuJDQ01C8Z4uPjmT59OrGxsTz44INXXg8PD2f9+vVXHdurVy+effZZKlaseNXrkZGRhISEpPhYuXJlijlEhM8++4x27dqR28cf91SQEoE33oB+/eCpp2DaNPDj+PZgHC3gd716wfr1uciaNXDXrFIFhg9P3zl69uxJlSpVCA8PB+wYzcjISBo3bkyZMmUYMmQIy5cvp2zZsunOe60tW7YQHh7O+fPnCQkJYfbs2dxzzz1X3g8PD2fkyJFXnn/zzTds2LCBGTNmXHeu559//oYf5YsXL57i+0uWLGHPnj08++yzqfyTqAxJBAYMsKMCOnaETz/F3/+Da3HNIHr37s0PP/zADz/8QNYkvzT169enWrVqDBgwgLlz51KtWjWP3x8ZGUlkZOSV5xcuXMAYw7Bhw668tmDBAh5++GGP31+uXDk2btzIqVOnmDVrFh06dCA6OppKlSoBUKNGDV566SVOnDhBnjx56NOnDwMHDuTmm2++7lyFChWiUKFCafo5/G3cuHFUq1aNKlWqpOs8KgMQgf/8B957D557DsaMgUBMdffUV5ARHxm5z7Vr165StGhR2bFjx3XvLVu2THLnzi1ZsmSR9evXJ3uO48ePy65du648WrRoIT169Ljqtbi4OK8z1alTRzp37nzl+fnz5yV79uyyaNEiefvtt6VcuXJy8Zr+rr8NGTJE8uTJk+JjxYoVyV77999/l2zZssnYsWO9zpuU9rn6nmNZExJEevSwfazdu9vnKfBln6u2XINcz549mTlzJt9//z3ly5e/6r1NmzbRokULPvroI+bNm0e/fv1YtGiRx/Nc21rMmzcvhQoV4s40zlZJSEjgwoULV57nyJGDypUrM3fuXCZOnMjUqVOTXc8hvd0CEyZMIEeOHLRp0yZN2VUGkZAAXbvCJ5/ASy/B0KEQwOnMWlyDWLdu3Zg0aRJTp06lYMGCHD16FICQkBCOHz9Oo0aN6N27N507d6Z69epUrlzZLzNl+vbtS+PGjSlZsiRnzpxh6tSpREdHXzfWtXr16owaNYp69erRpEmTZM+Xnm4BEeHTTz+lTZs2QTEbTPlJfLztAvj8c+jb1/a1BnidCC2uQWz06NEANG3a9KrXe/TowZIlS2jSpAkDBw4EoFKlSrRq1Yp+/frx448/+jTH0aNHadeuHUePHiV//vxUrlyZBQsW8Oijj1513D333EOWLFl4//33fXr9pKKjo9m1axeTJ0/22zWUy12+bG9aTZkCgwbZhwML8GhxDWK2u8f7Oftffvml1+eeMGGCz4+dMWMGXbp04e677/b63KlVq1atKz8XlQldugTt2sGMGTBkCLz6qmNRtLgqv0pISODYsWNMmDCB7du3M2vWLKcjqYzq4kVo08bOuBo6FPr0cTSOFlflVytWrKB27dqUK1eOyZMnX5ngoJRPnT8PLVvCvHnw4YfQo4fTibS4Kv+qWbMmCQkJAEGxlqsKQnFx8PjjsHixHRkQEeF0IkCLq1IqmJ09C02bQnQ0jB9vlw90CS2uSqngdPo0NG4Mq1bBpEnQtq3Tia6ixVUpFXxOnbLrsK5fD9OnQ6tWTie6jhZXpVRwOXEC6teHzZth5kxo3tzpRB5pcVVKBY9jx6BePfjlFzvkysW7TGhxVUoFh99/hzp14Lff7LYs9es7nShFWlyVUu53+LAtrPv327GstWs7neiGtLgqpdztwAFbTI8ehYULIZk1hd1Gi6tSyr327LGF9cQJWLIEatRwOpHXtLgqpdwpJsYW1thYWLYMqlZ1OlGq6AaFymuRkZEYY+jevbsj1x81ahSVK1cmX7585MuXj/Dw8OvWjFUZxC+/wCOPwLlzsHx50BVW0OKqvLR69WrGjRtH5cqVHctQokQJ3nnnHX7++WfWrVtH7dq1ad68OZs3b3Ysk/KDrVuhZk27LmtUlN29MwhpcQ1iM2fOJEeOHOzfv//Kaz179qRMmTL8/vvvPrvOX3/9Rdu2bfnss888rmpVokSJ6xbA3rJlCzlz5mT79u0+y9GsWTMaNmzInXfeyV133cWQIUPImzevzxf/Vg7atAlq1bIbCH7/PSRucBmMtLgGsZYtW3LPPfcwdOhQAIYNG8a0adNYuHAhoaGhPrtOREQELVu2pHYyw1/Cw8NZu3btVa/16tWLZ599looVK171emRkJCEhISk+Vq5cecNM8fHxTJ8+ndjYWB588MG0/+GUe6xbZwtrrlywYgVcsydcsNEbWp706kWu9ev9vq/5VapUgeHDU/UtxhgiIyNp3Lgx5cuXZ8iQISxfvpyyZcv6LNa4ceOIiYlh0qRJyR4THh5+ZcsZgG+++YYNGzYwY8aM645N7+aDW7ZsITw8nPPnzxMSEsLs2bO55557vPiTKFf78Ue7VkChQraP9fbbnU6Ublpcg1z9+vW5//77GTBgAHPnzqVatWoejxswYABDhgxJ8VxRUVFXbV64c+dOXn31VVauXEn27NmT/b4aNWrw0ksvceLECfLkyUOfPn0YOHAgN99883XHpmfzQYBy5cqxceNGTp06xaxZs+jQoQPR0dFUCuKPj5neypXQqBGEhtrCWqqU04l8QourJ8OHc87Lfamctnz5crZu3YqIpNgV0KtXL9q1a5fiuUpd80v9448/8ueff15VuOLj41mxYgVjxozh7Nmz5MiRg7CwMLJnz866devYsGEDN910E926dfN4jcjISCIjI1PMsWDBAh5OZqB49uzZr2z3XbVqVdauXcsHH3zAZ599luI5lUtFRUGTJlCypB1ulcKnlmDjuuJqjGkAjACyAp+KyNvXvF8T+BbYk/jS1yLyRiAzusWmTZto0aIFQ4cOZdmyZfTr149FixZ5PLZw4cIULlw4Vedv3rw5Va8ZAtOpUyfKli3Lq6++eqU1myNHDu677z7mzp3LxIkTmTp1KtmyZfN4zvR2C1wrISGBCxcueH28cpElS+Cxx6BMGVtYfXifwA1cVVyNMVmBUUA94CCw1hgzR0SuveW8UkSS3/g+E9i3bx+NGjWid+/etG/fnn/+859UrlyZ6Ojoqz7ap0eBAgUoUKDAVa/lyZOHQoUKXfcxPDw8nBEjRlCvXj2aNEn+ryY93QJ9+/alcePGlCxZkjNnzjB16lSio6N1rGswmj8fWrSwN62WLIEiRZxO5HNuGy1QHYgRkd0ichGYDjRzOJPrnDhxggYNGtCkSRMGDhwIQKVKlWjVqhX9+vVzJFOVKlXIkiXLdUOyfOno0aO0a9eOcuXKUadOHdauXcuCBQto2LCh366p/OCbb+warJUq2T5WFxRWEfjf/2D7dt91BRo37fFujGkJNBCRZxOftwceEJHuSY6pCczCtmwPA31EZFsy54sAIgBCQ0PDpk+f7vG6+fPnv9KP97f4+HiyBnK0QDq4IWuzZs0oU6ZMisXVDTm9cW3OmJgY/vrrLwcTeRYbG0tISIjTMbzyd9Yi0dFUePNNYu+6i83vvstlh/NfumSIji7CrFkl2LkzH9Wq/c677+5I1Tlq1aq1XkSun0ImIq55AK2w/ax/P28PfHTNMfmAkMSvGwG7vDl3WFiYJGf79u3XvXb69Olkj3cbp7LGx8fL0aNH5e2335aiRYvKiRMnUjw+WH6m1+b09PvhBlFRUU5H8FpUVJTI5MkiWbKI/OMfIg7/Lhw7JjJkiEixYiIgUq6cyMcfi8yf/32qzwWsEw81x23dAgeBkkmel8C2Tq8QkdMiEpv49XwgmzEmdXdqlE+sWLGCW2+9lQkTJjBr1iyPs7eUAii6YAG0bw///KddNtChkTjbttmdt0uWhP79bc/E/PmwfTs8/zzkypXgs2u56oYWsBYoa4y5HTgEtAGeSnqAMaYo8LuIiDGmOrbf+HjAkypq1qxJQoLvfhlVBvXJJ5R/9127Pcs330Du3AG9fEICLFpk5+gsXgw5c8LTT8OLL8Ldd/vvuq4qriJy2RjTHViEHYo1XkS2GWOeT3x/DNASeMEYcxk4B7RJbJorpdzmo4/gxRc5XqMGN8+ZYytbgJw9C198ASNGwM6dUKwYDBliW66pHJWYJq4qrnDlo/78a14bk+TrkcDIQOdSSqXSe+9Bnz7QvDlbu3blkQAV1gMHYNQoGDsWTp60qxVOmQItW0IKEw19zm19rkqpjOCtt2xhbd0aZsxAkplU4kurV0ObNnZZgqFD7TrbP/wAa9bAU08FtrCCC1uuSqkgJgKvv24f7drB55/DTf4rM5cuwddfwwcfwE8/Qb580KsX9OgBt93mt8t6RYurUso3RODVV+Htt6FTJxg3zm8ry504YU8/ciQcPAh33mm7dzt0cGwgwnW0uCql0k/EdgO8/z506QKjR9sFr31s5057g2riRIiLsx/9R4+Gxo39crl00eKqlEqfhATo2dM2I3v0sNXPGJ+dXsQuPzB8OCxYADlyQNu29pIO7jp0Q1pclVJpl5BgR9+PGwcvvWTvJPmosJ47B5Mn26K6fbtdNOv11+3lbrnFJ5fwK5c1pFVGN3jwYJ8tbG2M4auvvvLJuVQaxMfDM8/Ywtq/v88K66FD9nQlS9oxqdmzw4QJsG8fDBwYHIUVtLgGtUOHDhEREUH58uXJnj07xYsX57nnnuPgwYNOR/Opjh07elzG8MiRIzRt2tSBRIrLl+00pwkTbHPyzTfTXVjXrbMf90uXtiO5Hn4YoqPh55/tjaocOXwRPHC0uAapPXv2ULVqVbZu3cqYMWOIiYlh8uTJbNu2jWrVqrF37950nf/ixYu+CepHRYsWJUew/R+XEVy6BE8+CVOn2pEBictepsXlyzBrFvzjH1CtGsydC926QUwMzJ4Njzzi0+7bgNLi6iNTpkyhdOnSZMmShdKlSzNlyhS/Xq9bt25kyZKFpUuXUrNmTUqVKkWtWrVYunQpWbJkuWqblZo1a9K9e/ervv/a1mDNmjV54YUX6NOnD0WKFOGhhx7yeN0DBw7QrFkzChUqRO7cuSlfvjxJl3LcsmULdevWJVeuXBQqVIiOHTumuFyfp1Zp0q6DwYMHM3HiRObNm4cxBmMM0dHRwPXdAje69t/XGjFiBMWLF6dgwYJ06tSJuLi4ZPOpa1y4AK1awVdf2ZEBr7ySptOcOmUncN15p505dfiwHat68KDtY73jDp+mdoQWVx+YMmUKERER7Nu3DxFh3759RERE+K3AnjhxgoULF9KtWzdyX7MIRu7cuenatSsLFizg5MmTqTrv5MmTERFWrlzJF1984fGYrl27EhcXR1RUFNu2bWP48OFXdiuIi4ujQYMGhISEsGbNGmbPns2qVavo3Llzmv6cAH369KF169bUrVuXI0eOcOTIEY9baXt77ZUrV7J161aWLl3Kl19+yezZsxkxYkSa82Uq58/b3QO+/daODPj3v1N9il277ICCEiXsyK3bbrMt1F277OD/fPl8H9spOlrAB/r3739d6ycuLo7+/fvTtm1bn19v165diAgVKlTw+H7FihUREXbt2kX16tW9Pu/tt9/Oe++9l+Ix+/bt44knnuDee++98j1/mzJlCrGxsUyaNOnK5o5jx46lVq1axMTEpLiBYnJCQkLIlSsXOXLkoGjRosked6Nr/70Yer58+fj444+56aabqFChAq1atbqy/5hKQVyc3T1g6VI7af+557z+VhH4+ecCvP8+fPcdZMtmexV69oT77vNfZKdpy9UH9u/fn6rXfcUk0xn19yJhyb2fnLCwsBse07NnT958803Cw8MZMGAA69evv/Lejh07qFy58lW75j744INkyZKF7duv3QbNt7y9dsWKFbkpyXTMYsWK8ccff/g1W9CLjbWj9JcuhfHjvS6s58/bw++9F156qQqrV8Nrr9m7/hMmZOzCClpcfeLaLalv9Hp6lS1bFmMM27Z53N2GHTt2YIyhTJkyAGTJkoVrV2W8dOnSdd+XJ0+eG177mWeeYc+ePXTq1Ilff/2VBx98kMGDBwO2qCdX0JN73dtsN+Ltta/dldYYo2vSpuT0aWjQAFautINOO3a84bccPQqDBkGpUnakFsB//vML+/fbgQUpfADJULS4+sCQIUM89n0OGTLEL9crVKgQjz76KKNHj/bYHTFq1CgaNmx4ZZfVIkWKcOTIkauO27RpU5qvX6JECSIiIpgxYwZvvPEGY8eOBWyrcNOmTZw5c+bKsatWrSIhISHZLgxP2TZu3HjV8+zZsxMfH59iprRcW93AqVNQv75dEWXaNLu0VAo2bLBDpkqVgjfegAcesI3dTZugUaOjgVzK1RW0uPpA27ZtGTt2LLfddhvGGG677TbGjh3rl/7Wv40cOZLLly9Tt25dvv/+ew4cOEB0dDT16tVDRBg58v+XvK1duzYLFixgzpw57Ny5k969e3PgwIE0Xbdnz54sXLiQ3bt3s3HjRhYuXEjFihUB+3PIkycPTz/9NFu2bGHFihV06dKFFi1aXLcBZNJsGzZsYPz48cTExPDuu+/yv//976pjSpcuzdatW9m5cyd//vmnx5ZtWq6tUnD8ONSpYweZzpxpRwh4EB9vNxeoWRPuv98Oq+rSBX791Q6rqlMneIdSpZcWVx9p27Yte/fuJSEhgb179/q1sAKUKVOGdevWcffddxMREcEdd9zBU089RYUKFVi7du1VN5o6d+585fHQQw8REhLC448/nqbrJiQk0KNHDypWrEi9evUIDQ1l4sSJgG2tL1q0iNOnT1O9enWaNWtGeHg448ePT/Z8jz76KIMGDaJ///6EhYWxd+9eunbtetUxzz33HBUqVKBq1aoUKVLkuuKb1murZBw7ZldE2bbt/7fBvsbp03YJgbvugscfh7177QStgwft6lRlywY6tAt52rUwIz5091fnBWvOTLX765EjIhUriuTKJbJ48XVv//abSK9eInnz2l1TH3pI5KuvRC5dciCrH6QlJ8ns/qpDsZRS1qFDtsV66JDdErVmTcAOpVq50g7y//Zbu0Rr69Z2XGq1ao4mdjUtrkop2L/fFtY//rBbpT70EBcuwJdf2hlTGzZAoULQrx907QrFizsd2P20uCqV2e3ebQvrqVOwZAl/3P4An/zXLkJ99ChUrGjnDbRtG/BdsYOaFlelMrNdu2xhjYtj1yfLeXvs/UyZYpcQaNjQfvSvVy/z3vFPDx0tkEiuGciuFGTM34ukiwzdVr48E/88RcSdy7mrzf1Mm2a3v9qxw3a71q+vhTWttOUKZM2alUuXLpE90HvvKtc7d+7cdbO6gtnfiwz9Pflkvwidzl8m/69beeute4mIsH2rKv205QoUKFCA33//XadBqitEhLi4OA4dOsQtwbL0vRdeeeX6RYaE8+TL15++fbWw+pK2XIHChQtz8OBBdu7ceeW18+fPkzNI5usFS9Zgy5ktWzZCQ0PJF+Tr4InAqlX2rv+hQ54XEzpwwL+LDGVGWlyxi4dcu8hKdHQ09wXJsj3BklVzBtbFi3bm6vDhdguVeiE/UhLwNPHZX4sMZWbaLaBUBnP8OERGwu23Q7t2dqrq7H+vYBH1eeuWIuTOleuq4/25yFBm5rriaoxpYIzZaYyJMcb09fC+McZ8mPj+ZmPM/U7kVMpttm+3i6aUKGF3T61YEebNgx2jltP8k4aYEiVou3EjY8eNC+giQ5mVq7oFjDFZgVFAPeAgsNYYM0dEkq603BAom/h4APg48b9KZToJCbB4sZ2aungx5MwJ7dvDiy9CpUrY2VZNm9vNqpYuhdBQ2rZtq8U0AFxVXIHqQIyI7AYwxkwHmgFJi2sz4IvEBRNWG2MKGGNuFZEj159OqYwpLg7mzCnGCy/AL7/Arbfa3a27dIHChRMP+u47eOIJ24RdsiTJGyoQjJsGSRtjWgINROTZxOftgQdEpHuSY74D3haRHxKfLwNeEZF1Hs4XAUQAhIaGhiXdpfRGYmNjCQkJSc8fJ2CCJavmTL9jx3LwzTfFmDu3GGfOZOOuu87QsuUBatY8RrZs////cuGVK6n4xhvElinD5nff5bLDIx7c/DNNKi05a9WqtV5Eql73hqelspx6AK2AT5M8bw98dM0x84B/JHm+DAi70blTWnLQk2BZIk0keLJqzrRbvVqkTRuRrFlFsmQRadlS5MMPf5aEBA8HT59uD6xRQ+TUqYBn9cSNP1NPfLnkoNtuaB0ESiZ5XgI4nIZjlAp6ly/DjBnw4INQo4adjtqrF/z2mx1idc89f10/NXXSJLsdy4MP2k7Y/PmdiK5w32iBtUBZY8ztxpjsQBtgzjXHzAGeThw1UAP4S7S/VWUgJ0/Cu+/CHXfAv/5lNwb48EO7yv+wYVC6dDLfOH683cSqZk1YsACS7ISrAs9VN7RE5LIxpjuwCMgKjBeRbcaY5xPfHwPMBxoBMUAc0MmpvEr50s6dtohOmGBvWNWqBSNHQpMmkOVGzaCPP7YLrT76KMyeDdeMZVWB56riCiAi87EFNOlrY5J8LUC3QOdSyh9E7Aip4cPtx/7s2e2n+l694N57vTzJiBH2G5o2tf0FOXL4L7DymuuKq1KZwblzMGWKLarbtsEtt8Drr8Pzz9uvvTZ0KLz8MrRoYbe/1pXdXEOLq1IBdPiwXeF/zBg7TbVKFdsN0KZN6huct02aZPtZ27SBL76ADLQ0YkagxVWpAFi/3s6i+vJLiI+HZs3sJ/l//jMNi1GLwKBB3D5+vJ2O9fnndtdA5SpaXJXyk8uX7W6pw4fDDz/Ym/fdukGPHlCmTBpPKgJ9+8K773KkUSNu1cLqWlpclfKxv/6CTz+Fjz6Cffvs6lQffGC3T0nXsFMR6N3bVusXXmBny5bcqoXVtdw2zlWpoLVrl10wpUQJ6NMHbrsNvv7avt6rVzoLa0KCbfYOHw49e8KoUV6Mz1JO0parUukgAtHRtmX63Xdw003w5JO2/t3vq8UwExIgIgI++8yODHj7bd01MAhocVUqDc6ftyOfRoyATZvsglP9+9tx/Lfe6sMLxcdD5852NMBrr9nxWlpYg4IWV6VS4ehRO4zq44/hjz/gnntsg/Kpp+xaqj516RI8/TRMnw7//S8MGODjCyh/0uKqlBc2brTdndOm2b2pGjeGf/8batf2U0Py4kVbsWfNgnfesd0BKqhocVUqGfHxth/1gw/g++8hd2547jl70+quu/x44QsXoFUrmDvXXrxXLz9eTPmLFlelrnH2bFZGjLCLqOzeDaVK2VmmzzwDBQv6+eLnztmprAsX2qlcL7zg5wsqf9HiqlSiPXtsQR07Npy4OLsk6jvvQPPmdhSA38XF2alby5bZgbLPPBOAiyp/0eKqMjURWLnSfvr+9ls72emRR44TGRlK9eoBDBIba9cWXLnSLjbw9NMBvLjyBy2uKlO6cMHO8x8+HDZsgEKFoF8/O5Rq164dVK8eGrgwf/0FjRrBTz/B5Ml2oKwKelpcVaZy7JgdSjV6tB1WVaGCfd6+vb1hBXZGVcCcPGkXuN6wwQ65atkygBdX/qTFVWUKW7bYAf+TJ9tWa4MGdihVvXoOjsk/ftwG2LbNDrl67DGHgih/0OKqMqyEBLuV1PDhdrX/XLmgY0c7NbVCBYfD/fEH1K0Lv/5qO3sbNHA4kPI1La4qw4mNhYkTbUt11y4oXtxOx3/uOdu36rgjR6BOHdi71w6krVvX6UTKD7S4qgxj/367od+4cXDqFFSvbmdUPfGEixbpP3jQTus6fNg2qx95xOlEyk+0uKqgJgI//mg/+n/9tX2tZUs7qalGDSeTebBvny2sx47B4sV2IK3KsLS4qqB06RJ89ZUdn7p2LRQoAC+9ZJc8LVXK6XQe7N5t98o+fdp2AAd0EK1yghZXFVSOH4exY+3H/8OH7Rz/0aPtmPs8eZxOl4xff7Ut1nPn7Owrny30qtxMi6sKCjt22I/+X3xh11KtV8/2rTZo4PIF+bdvtzev4uMhKgoqV3Y6kQoQLa7KtRISbNfk8OGwaJHderp9ezuUqlIlp9N5YcsWW1izZrXbFVSs6HQiFUBaXJXrxMXBpEm2qP7yi13Z/7//hS5doEgRp9N56eefbfM6Vy5YvtzPaxQqN9Liqlzj4EG7797YsXDihO2anDQJWreG7NmdTpcKa9bYKa358tnCmuZ9tFUw0+KqHLdmjW2lzpxpuwKaN7dTUx96KAi3i1q1ynYEFyliC+tttzmdSDnENcXVGFMI+BIoDewFWovISQ/H7QXOAPHAZRGpGriUylcuX7bjUocPt+NU8+Wzfandu0Pp0k6nS6Pvv7f7vxQrZgtriRJOJ1IOctN91r7AMhEpCyxLfJ6cWiJSRQtr8Dl5Et59F+64A/71LzvF/sMPbZfAsGFBXFiXLoWGDW1L9fvvtbAq97RcgWZAzcSvJwLRwCtOhVG+tXMnfPBBWZYutTesata0Y1UbN7Y304PawoW2L+Ouu2yRveUWpxMpF3BTyzVURI4AJP43ud9QARYbY9YbYyIClk6lmoitNY0bQ/nysGDBrbRubXdSjYqyK+wFfWGdO9duzVKxov1DaWFViYyIBO5ixiwFinp4qz8wUUQKJDn2pIhctx2cMaaYiBw2xtwCLAF6iMiKZK4XAUQAhIaGhk2fPt3rrLGxsYSEhHh9vJPclvXChSwsXRrKV1+VYO/ePBQseJFmzQ5Ru/YuSpZ0/21/b3+ehVesoOIbbxBbtiyb332Xy3nzBiDd/3Pb33tKgiVrWnLWqlVrvccuShFxxQPYCdya+PWtwE4vvmcw0Meb84eFhUlqREVFpep4J7kl66FDIgMGiBQuLAIi994r8vnnIufP2/fdkvNGvMo5bZpI1qwi4eEip075PZMnwfLzFAmerGnJCawTDzXHTd0Cc4AOiV93AL699gBjTB5jTN6/vwbqA1sDllB5tH69nTlVujQMGWIXe4qKsjuXdOxoZ1ZlKJMmQdu2dqzYokWQP7/TiZQLuemG1tvADGPMM8B+oBXYbgDgUxFpBIQCs40d/HgTMFVEFjqUN1OLj7cL6H/wAfzwA4SEwAsvwIsvZvAx8599Zlfdrl3b/gBcu1qMcppriquIHAfqeHj9MNAo8evdwL0BjqaS+OsvGD/eDp/au9e2Vt9/Hzp3zgQNuNGj7ZqGDRrYQbq5cjmdSLmYa4qrcreYGPjoI1tYY2Phn/+E996zN8qD/o6/N0aMsCtwN21qp5JluL4O5WtaXFWyROxiTsOH2xFHN90EbdrYmVRhYU6nC6B33oG+fe1+MVOnBtlCB8opWlzVdc6fh+nTbVHdtAkKF4YBA2yf6q23Op0uwP77Xxg40P6rMmmS/RdGKS/c8DfFGFNeRH4JRBjlrN9/h48/to8//rBrpn76KTz1VCbsXhSB116zwx+eftr2h2SK/g/lK978M7zBGDMWGCweFlJRwW/jRtulOHUqXLxoZ1T16mXXeQ66Val8QQReftkudvDss/DJJy7f7kC5kTe/MdWBu4Fdxpgexhj95zsDiI+HOXPsnnn33QczZtgRRjt3wnffQd26mbew3jlypC2sXbtqYVVpdsPfGhHZIiJ1gWeBF4EtxpiGfk+m/OLMGTuMqlw5e6f/t9/s/ZqDB+1CKpl6wfyEBOjalRJff20XlB05UgurSjOvf3NE5BtsC3YiMN0YM98YU95fwZRv7dkDvXvblfB69rTri8yYYXd8fvllKHjdKg6ZTHy8bbqPGcO+p56y48wyZdNd+Upqb33mBtZjC2w3YLMxZgzwmoj85etwKn1E7Oyp4cPhm29sI6xVK9ufWr26w+Hc5PJl6NQJJk+GgQPZU7Mmt2lhVel0w5arMaaXMWaKMeZX4DgwF6gGjMB2FZQDthtjHvBrUuW1ixftqKGqVe1g/+hoeOUV23qdOlUL61UuXbLrBEyeDG++Ca+/ri1W5RPetFxfAn4EPgZWA+tF5GKS978wxrwCjMd2GyiHHDtm77+MGgVHj0KFCvZ5u3aQO7fT6Vzo4kU7fnX2bBg6FPr0cTqRykBuWFxFpKQX5/kciEx/HJUWu3fnYfJk2/i6cMHuNtKrl93ZWRthyTh/3vaRfPedHYf24otOJ1IZjK+mmxwDavvoXMoLCQmwYIFdlWrZsmrkymWX9+vZ07ZYVQrOnYPHH7fLBX78MTz/vNOJVAbkk+KauGDs9744l0pZbCx88YVtbP36KxQvDs89t5u33rqDm292Ol0QOHvW7i8TFWWnnz3zjNOJVAalg/iCxP79dshUyZJ21bsCBezNqT174Kmn9mth9caZM7bPJDra/gulhVX5ka5C4WIisHq1HUo1a5Z9rUULO749PNzRaMHnr79sYV2zxv6r9K9/OZ1IZXBaXF3o0iX46itbVNessa3U3r2he3coVcrpdEHo5EmoX98u8TVjhv0XSik/0+LqIsePw7hxdtbloUN2KuqoUXZRpiDYONOd/vzTDpvYvt3uHtCkidOJVCahxdUFduywN6i++MLeyK5bF8aOtbuJ6NT2dPj9d/vDjImx+101aOB0IpWJaHF1iAgsXmyHUi1aZHcNad/eDqWqVMnpdBnA4cN2zcT9+2HePLuhoFIBpMU1wOLi7NTUESNsi7VoUbvYfZcuUKSI0+kyiAMHbDE9ehQWLoSHH3Y6kcqEtLgGyKFDti/1k0/s/ZX777fdAK1b6153PrV3ry2sx4/bjwY6rEI5RIurn61ZY+/6z5xpZ1U1b26npv7jHzo11ed++82u/n3mDCxdCtWqOZ1IZWJaXP3g8mW7Fsjw4bBqFeTNCz162MfttzudLoPaudO2WC9cgOXL7fYKSjlIi6sPnTxpZ1R+9JHt9itTxhbYTp0gXz6n02Vg27bZm1d/7wWudwSVC2hx9YGdO+3WKRMm2BtWtWrZAtukiW4Y6nebNtnhVtmy2RZred0cQ7mDFtc0EoFly2zLdN48yJ7dbkHdqxfce6/T6TKJ9evtBIE8eWxhLVvW6URKXaHFNZXOnbNT04cPh61b7V5UgwbBCy9AaKjT6TKRn36CRx+1c4OjorQzW7mOa+b/GGNaGWO2GWMSjDFVUziugTFmpzEmxhjTN1D5jhyB116zc/uffdZ+3P/8cztGffBgLawB9cMPtsVauDCsWKGFVbmSm1quW4EWwCfJHWCMyQqMAuoBB4G1xpg5IrLdX6F+/tm2UqdPt6MAmja1q1I98ogOpXJEdLTtzC5RwvbLFC/udCKlPHJNcRWRHQAm5YpVHYgRkd2Jx04HmgE+La7x8bByZWEGDbINo5AQ+7G/Rw+4805fXkmlypIl0KyZbakuW2antynlUq4prl4qDhxI8vwg4PNdZ/v2hWHDKlG6NLz/PnTuDPnz+/oqKlXmz7dLBZYrZycI6Fxh5XLG7tASoIsZsxTw1NzoLyLfJh4TDfQRkXUevr8V8KiIPJv4vD1QXUR6JHO9CCACIDQ0NGz69Ole5dy/Pxe//JKVOnVig2IoVWxsLCFBsCZhWnPe/L//cffgwZy9/XY2DR3KZT//S5fRf55OCJasaclZq1at9SJy/X0iEXHVA4gGqibzXjiwKMnzfkA/b84bFhYmqREVFZWq450ULFnTlHPmTJGbbhKpXl3k5ElfR/IoQ/88HRIsWdOSE1gnHmqOa0YLeGktUNYYc7sxJjvQBpjjcCblL1OnQps28MADtr+1QAGnEynlNdcUV2PM48aYg9jW6TxjzKLE14sZY+YDiMhloDuwCNgBzBCRbU5lVn40caJd4Pbhh+2ygTp/WAUZ19zQEpHZwGwPrx8GGiV5Ph+YH8BoKtDGjbML3NapY3cQyJ3b6URKpZprWq5KAXbTsIgIuyXL3LlaWFXQ0uKq3OODD+wWt489ZtdszJnT6URKpZkWV+UOb79t9w9v2dLuK67bM6ggp8VVOUsE3ngD+vWzy4pNm2aXD1QqyLnmhpbKhERgwACIjIQOHeCzz3QBXJVhaHFVzhCBl1+GYcPguedgzBjIoh+kVMahxVUFngj07Gm3a+jWzW7joIVVZTD6G60CKyEBnn/eFtbeve1/tbCqDEh/q1XgxMfblcbHjrU3sIYN00VxVYal3QIqMC5fpsJbb9l1WAcNsg8trCoD0+Kq/O/SJWjbltBly+zIgH79nE6klN9pt4DyrwsXoHVrmDmTmBde0MKqMg1tuSr/OX/ezriaNw8++oiDlSqhu+SozEJbrso/4uLsflfz58Mnn9g1A5TKRLTlqnzv7Fm7TW50NIwfDx07Op1IqYDT4qp86/RpaNwYVq2CSZOgbVunEynlCC2uyndOnYKGDWHtWrsAS+vWTidSyjFaXJVvnDgB9evD5s12ycDmzZ1OpJSjtLiq9Dt2DOrVg19+sYtcN27sdCKlHKfFVaXP0aNQty789hvMmWNbr0opLa4qHQ4dspsIHjhgh1zVquV0IqVcQ4urSpv9+6F2bfj9d7v19cMPO51IKVfR4qpSb88eW1hPnoQlS6BGDacTKeU6WlxV6sTE2MIaG2tXuAoLczqRUq6kxVV575dfbGG9dAmWL4cqVZxOpJRraXFV3tm61d68MsZOa737bqcTKeVqunCLurGNG6FmTbszqxZWpbyixVWlbN062xWQOzesWAHlyzudSKmg4JriaoxpZYzZZoxJMMZUTeG4vcaYLcaYjcaYdYHMmOmsXm0nCOTPbwvrnboaq1LeclOf61agBfCJF8fWEpE//Zwnc1u5Eho1gqJF7c2rkiWdTqRUUHFNcRWRHQBGN61zXlQUNGkCpUrZ4VbFijmdSKmg45pugVQQYLExZr0xJsLpMBnO4sW2xXr77fbmlRZWpdLEiEjgLmbMUqCoh7f6i8i3icdEA31ExGN/qjGmmIgcNsbcAiwBeojIimSOjQAiAEJDQ8OmT5/uddbY2FhCQkK8Pt5Jvspa6McfqTRoEHGlSrFp2DAuFSiQ/nBJBMvPVHP6XrBkTUvOWrVqrReR6+8TiYirHkA0UNXLYwdjC/ENjw0LC5PUiIqKStXxTvJJ1tmzRbJlEwkLEzl+PP3n8yBYfqaa0/eCJWtacgLrxEPNCapuAWNMHmNM3r+/Bupjb4Sp9Jg5E1q1gvvvh6VLoVAhpxMpFfRcU1yNMY8bYw4C4cA8Y8yixNeLGWPmJx4WCvxgjNkErAHmichCZxJnEFOmQJs2dvGVJUvAx10BSmVWbhotMBuY7eH1w0CjxK93A/cGOFrG9fnn8Mwz8Mgj8N13kCeP04mUyjBc03JVATZ2LHTubCcJzJunhVUpH9PimhmNHAlduti9rubMsVNblVI+pcU1s3nvPejRw+7O+vXXkDOn04mUypC0uGYmkZHQpw+0bg0zZkD27E4nUirD0uKaGYjA4MHQvz+0a2dHCGTL5nQqpTI014wWUH4iYovqW29Bp04wbpxdl1Up5VdaXDMyEfjPf2w/a5cuMHo0ZNEPK0oFghbXjCohAXr2tCMDevSAESPsFi1KqYDQ4poRJSTA88/bLoCXXoKhQ7WwKhVg+hkxo4mPt7Ouxo2DV1/VwqqUQ7TlmpFcvgwdOsDUqfD66/Daa1pYlXKIFteM4tIleOop+OorOzKgb1+nEymVqWlxzQDMxYvQsqWdyvr++/DvfzsdSalMT4trsDt3jkoDB8JPP8FHH0H37k4nUkqhxTW4xcVBs2YUWrMGPvkEInRLMaXcQkcLBKvYWLuq1bJl7Hz5ZS2sSrmMtlyD0enTdofW1ath8mSOFitGeaczKaWuoi3XYHPqFNSvb/tYp02zIwSUUq6jLddgcvy4LaxbtthNBZs3dzqRUioZWlyDxbFjdkuWnTvhm29st4BSyrW0uAaDo0ehTh3YswfmzoV69ZxOpJS6AS2ubnfoENSubf87fz7UrOl0IqWUF7S4utm+fbawHjsGixbBQw85nUgp5SUtrm61e7ctrKdOwZIl8MADTidSSqWCFlc32rXLFta4OFi2DMLCnE6klEolLa5us2OHvXl16RJERUHlyk4nUkqlgU4icJMtW+CRR+xOAtHRWliVCmJaXN1iwwaoVctuef3993D33U4nUkqlg2uKqzFmqDHmF2PMZmPMbGNMgWSOa2CM2WmMiTHGZIwVodeutX2sefLAihVQrpzTiZRS6eSa4gosASqJSGXgV6DftQcYY7ICo4CGQEXgSWNMxYCm9LVVq+zMq4IFbYu1TBmnEymlfMA1xVVEFovI5cSnq4ESHg6rDsSIyG4RuQhMB5oFKqPPrVgBjz4KoaH269KlnU6klPIR1xTXa3QGFnh4vThwIMnzg4mvBZ9ly6BBAyhRwrZYS3j6t0QpFayMiATuYsYsBYp6eKu/iHybeEx/oCrQQq4JZ4xpBTwqIs8mPm8PVBeRHslcLwKIAAgNDQ2bPn2611ljY2MJCQnx+vjUKLhmDZVee41zxYuzadgwLhUqlK7z+TOrL2lO3wqWnBA8WdOSs1atWutFpOp1b4iIax5AB+BHIHcy74cDi5I87wf08+bcYWFhkhpRUVGpOt5rc+eKZM8uUqWKyLFjPjml37L6mOb0rWDJKRI8WdOSE1gnHmqOa7oFjDENgFeAx0QkLpnD1gJljTG3G2OyA22AOYHKmG6zZ0OLFnb86rJlULiw04mUUn7imuIKjATyAkuMMRuNMWMAjDHFjDHzAcTe8OoOLAJ2ADNEZJtTgVPlyy+hVSuoWhWWLoV0dgUopdzNNdNfReTOZF4/DDRK8nw+MD9QuXxi0iTo2NGuajVvHuTN63QipZSfuanlmjGNHw8dOth1WBcs0MKqVCahxdWfxoyBZ56x+159952dgaWUyhS0uPrLiBHwwgvQpInd8ypXLqcTKaUCSIurPwwdCr162ZEBs2ZBzpxOJ1JKBZgWV18bMgRefhn+9S+YPh2yZ3c6kVLKAVpcfUUEBg6EAQOgfXuYPNkuH6iUypRcMxQrqIlAv37wzjvQuTOMHQtZszqdSinlIC2u6SUCvXvD8OHw/PMwahRk0Q8ESmV2WgXSIyEBune3hbVnTxg9WgurUgrQ4pp2CQnQpYstqP/5D3zwARjjdCqllEtocU2L+Hjbt/rpp/YG1jvvaGFVSl1F+1xT6/JlePppmDYN3ngDXnvN6URKKRfS4poaly7Bk0/aiQHvvGPHsyqllAdaXL114QK0bg1z5tj+1V69nE6klHIxLa7eOHcOnnjCrmo1ahR07ep0IqWUy2lxvZG4OGjWzO4c8OmndpUrpZS6AS2uKYmNtatarVwJEybYG1lKKeUFLa7JyBobC48+Cj/9ZNcJePJJpyMppYKIFldPTp7k3v/8B2Ji7N5XTzzhdCKlVJDR4urJf/9LyG+/2SFXjz3mdBqlVBDS4urJkCFsvOMO7tfCqpRKI53+6kmuXJyuVMnpFEqpIKbFVSml/ECLq1JK+YEWV6WU8gMtrkop5QdaXJVSyg+0uCqllB+4ZpyrMWYo0BS4CPwGdBKRUx6O2wucAeKByyJSNYAxlVLKK25quS4BKolIZeBXoF8Kx9YSkSpaWJVSbuWa4ioii0XkcuLT1UAJJ/MopVR6uKa4XqMzsCCZ9wRYbIxZb4yJCGAmpZTyWkD7XI0xS4GiHt7qLyLfJh7TH7gMTEnmNA+JyGFjzC3AEmPMLyKyIpnrRQB/F+BYY8zOVMQtDPyZiuOdFCxZNadvBUtOCJ6sacl5m6cXjYikP46PGGM6AM8DdUQkzovjBwOxIjLMD1nWBUufbrBk1Zy+FSw5IXiy+jKna7oFjDENgFeAx5IrrMaYPMaYvH9/DdQHtgYupVJKecc1xRUYCeTFftTfaIwZA2CMKWaMmZ94TCjwgzFmE7AGmCciC52Jq5RSyXPNOFcRuTOZ1w8DjRK/3g3cG6BIYwN0HV8Ilqya07eCJScET1af5XRVn6tSSmUUbuoWUEqpDEOLawqMMa2MMduMMQnGGNfd6TTGNDDG7DTGxBhj+jqdJznGmPHGmD+MMa6++WiMKWmMiTLG7Ej8e+/pdCZPjDE5jTFrjDGbEnO+7nSmlBhjshpjNhhjvnM6S0qMMXuNMVsS7/msS+/5tLimbCvQAvA4jtZJxpiswCigIVAReNIYU9HZVMmaADRwOoQXLgMviUgFoAbQzaU/0wtAbRG5F6gCNDDG1HA2Uop6AjucDuEln02t1+KaAhHZISKpmXgQSNWBGBHZLSIXgelAM4czeZQ4yeOE0zluRESOiMjPiV+fwRaE4s6mup5YsYlPsyU+XHnzxBhTAmgMfOp0lkDT4hq8igMHkjw/iAsLQbAyxpQG7gN+cjiKR4kftTcCfwBLRMSVOYHhwMtAgsM5vOHTqfWuGYrlFG+m5LqU8fCaK1svwcYYEwLMAnqJyGmn83giIvFAFWNMAWC2MaaSiLiqT9sY0wT4Q0TWG2NqOhzHG15PrfdGpi+uIlLX6QxpdBAomeR5CeCwQ1kyDGNMNmxhnSIiXzud50ZE5JQxJhrbp+2q4go8BDxmjGkE5ATyGWMmi0g7h3N5lDimHhH5wxgzG9v1lubiqt0CwWstUNYYc7sxJjvQBpjjcKagZowxwGfADhF53+k8yTHGFElssWKMyQXUBX5xNJQHItJPREqISGns7+dytxZWf0yt1+KaAmPM48aYg0A4MM8Ys8jpTH9LXPu2O7AIe+NlhohsczaVZ8aYacCPQDljzEFjzDNOZ0rGQ0B7oHbicJyNia0ut7kViDLGbMb+I7tERFw9zCkI+Hxqvc7QUkopP9CWq1JK+YEWV6WU8gMtrkop5QdaXJVSyg+0uCqllB9ocVVKKT/Q4qqUUn6gxVUppfxAi6tSXFkY/YIx5rYkr40wxvxmjAl1MpsKTjpDSymurCuwFtggIs8ZY/pgl8p7SER2OZtOBaNMvyqWUmAXoDbGvIpdQ+I3oD92tX8trCpNtOWqVBLGmFXYpeaaisgCp/Oo4KV9rkolMsbUBu7FLkT+u8NxVJDTlqtSgDHmXuB7oDd2z6cQEXnU2VQqmGlxVZle4giBVcAnIvKGMaYSsBnb5xrtaDgVtLS4qkzNGFMI+B+wQkS6JHn9S6CUiIQ7Fk4FNS2uSinlB3pDSyml/ECLq1JK+YEWV6WU8gMtrkop5QdaXJVSyg+0uCqllB9ocVVKKT/Q4qqUUn6gxVUppfzg/wBSvF4Vvf5wJQAAAABJRU5ErkJggg==\n",
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"x = np.linspace(-1,5,100)\n",
"\n",
"y1 = (2./8.)*x - (6./8.)\n",
"y2 = (5./9.)*x - (16./9.)\n",
"\n",
"fig = plt.figure(figsize=(5, 5))\n",
"\n",
"ax1 = fig.add_subplot(111)\n",
"\n",
"ax1.set_xlabel(\"$x$\", fontsize=14)\n",
"ax1.set_ylabel(\"$y$\", fontsize=14)\n",
"ax1.grid(True)\n",
"\n",
"ax1.plot(x,y1,'b', label='$2x+3y=7$')\n",
"ax1.plot(x,y2,'r', label='$x-4y=3$')\n",
"ax1.plot(37./11,1./11, 'ko', label='Our solution')\n",
"\n",
"ax1.legend(loc='best', fontsize=14)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"So the nature of each individual line is completely different, BUT their intersection, i.e. the point that is a member of both lines, is unchanged.\n",
"\n",
"So if we just want to find the solution to this problem it doesn't matter if we solve the original system\n",
"\n",
"$$\n",
" \\begin{pmatrix}\n",
" 2 & 3 \\\\\n",
" 1 & -4 \n",
" \\end{pmatrix}\n",
" \\begin{pmatrix}\n",
" x \\\\\n",
" y \n",
" \\end{pmatrix}=\n",
" \\begin{pmatrix}\n",
" 7 \\\\\n",
" 3 \n",
" \\end{pmatrix} \n",
"$$\n",
"\n",
"or this updated one\n",
"\n",
"$$\n",
" \\begin{pmatrix}\n",
" -2 & 8 \\\\\n",
" 5 & -9 \n",
" \\end{pmatrix}\n",
" \\begin{pmatrix}\n",
" x \\\\\n",
" y \n",
" \\end{pmatrix}=\n",
" \\begin{pmatrix}\n",
" -6 \\\\\n",
" 16 \n",
" \\end{pmatrix} \n",
"$$\n",
"\n",
"As we have seen in Numerical Methods, the value here comes when we can turn the matrix on the LHS into one that is trivial to solve the corresponding system for, e.g. the identity matrix (a diagonal matrix being almost as easy), or if we know how to do back substitution a lower triangular matrix (cf. Gaussian elimination)."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Problems with a zero RHS\n",
"\n",
"Question: what happens to the system undergoing row operations in the case when $\\boldsymbol{b} = \\boldsymbol{0}$?\n",
"\n",
"[It turns out (later) that this is an important question we can ask of our $A$ matrix irrespective of what the actual $\\boldsymbol{b}$ vector is.]\n",
"\n",
"\n",
" \n",
"\n",
"Suppose you do row operations to a system with a zero RHS ($A\\boldsymbol{x}=\\boldsymbol{0}$) and are able to end up with this system\n",
"\n",
"$$\n",
" \\begin{pmatrix}\n",
" 1 & 0 \\\\\n",
" 0 & 1 \n",
" \\end{pmatrix}\n",
" \\begin{pmatrix}\n",
" x \\\\\n",
" y \n",
" \\end{pmatrix}=\n",
" \\begin{pmatrix}\n",
" 0 \\\\\n",
" 0 \n",
" \\end{pmatrix} \n",
"$$\n",
"\n",
"what does this tell you about possible solutions to the problem $A\\boldsymbol{x}=\\boldsymbol{0}$?\n",
"\n",
" \n",
" \n",
"\n",
"\n",
"OK, now suppose instead of the final form above you end up with the following\n",
"\n",
"$$\n",
" \\begin{pmatrix}\n",
" 1 & 0 \\\\\n",
" 0 & 0 \n",
" \\end{pmatrix}\n",
" \\begin{pmatrix}\n",
" x \\\\\n",
" y \n",
" \\end{pmatrix}=\n",
" \\begin{pmatrix}\n",
" 0 \\\\\n",
" 0 \n",
" \\end{pmatrix} \n",
"$$\n",
"\n",
"what does this tell you about possible solutions to the problem $A\\boldsymbol{x}=\\boldsymbol{0}$?\n",
"\n",
" \n",
" \n",
"\n",
"Note that for a matrix $A$ the space $\\{\\boldsymbol{x} | A\\boldsymbol{x} = \\boldsymbol{0}\\}$ (which means the space of all possible $\\boldsymbol{x}$ vectors such that $A\\boldsymbol{x} = \\boldsymbol{0}$) is called the **null space** of $A$. \n",
"\n",
" \n",
"\n",
"So note that as well as solving linear systems, row operations give us a method for calculating the null space of a matrix (for this case you do not actually need to include the RHS in the augmented matrix, as we know it can never change from the zero vector under row operations), \n",
"\n",
"and the dimension of the null space helps us define the rank of a matrix (see later).\n",
"\n",
" \n",
"\n",
"Using row operations to find the null space and to calculate the rank are mentioned again at the end, and a worked example is given in the homework, which you should read through.\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### (Reduced) row echelon form of a matrix\n",
"\n",
"During the steps of Gaussian elimination where we perform row operation on a matrix to transform it one for which\n",
"\n",
"\n",
"- all non-zero rows are above any rows of all zeros (NB. if we do have rows of all zeros then the matrix is singular), and\n",
"\n",
"\n",
"- the leading coefficient in the row (first non-zero entry from the left) is always strictly to the right of the leading coefficient of the row above\n",
"\n",
"\n",
"If a matrix is in this form we say that it is in [**row echelon form**](https://en.wikipedia.org/wiki/Row_echelon_form).\n",
"\n",
"\n",
" \n",
"\n",
"- If in addition the leading coefficient in any non-zero row is unity and every column containing a leading '1' has zeros everywhere else within that column, we say the matrix is in **reduced row echelon form**.\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Solving our simple example using the inverse matrix\n",
"\n",
"So we have enough information now to solve our problem\n",
"\n",
"$$\n",
"\\begin{align*}\n",
" 2x + 3y &= 7 \\\\[5pt]\n",
" x - 4y &= 3 \n",
"\\end{align*}\n",
" \\quad \\iff \\quad\n",
" \\begin{pmatrix}\n",
" 2 & 3 \\\\\n",
" 1 & -4 \n",
" \\end{pmatrix}\n",
" \\begin{pmatrix}\n",
" x \\\\\n",
" y \n",
" \\end{pmatrix}=\n",
" \\begin{pmatrix}\n",
" 7 \\\\\n",
" 3 \n",
" \\end{pmatrix} \n",
"$$\n",
"\n",
"\n",
"Recalling the solution we calculated through substitution: $x=37/11$ and $y=1/11$, we can use to check our numerical solution.\n",
"\n",
"Let's calculate \n",
"\n",
"$$\\boldsymbol{x} = A^{-1}\\boldsymbol{b}$$"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {
"scrolled": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"x = [3.363636 0.090909]\n",
"\n",
"Our computed solution solves the linear system (Ax == b): True\n",
"\n",
"and agrees with the exact solution we computed by hand: True\n"
]
}
],
"source": [
"# form the LHS matrix\n",
"A = np.array([[2., 3.], [1., -4.]])\n",
"\n",
"# form the RHS vector\n",
"b = np.array([7., 3.])\n",
"\n",
"# the inverse matrix from above:\n",
"invA = (-1./11)*np.array([[-4, -3], [-1, 2]])\n",
"\n",
"# compute and print our solution\n",
"x = invA@b\n",
"print('x = ',x)\n",
"\n",
"# Check our two numerical solutions using the numpy allclose function\n",
"print('\\nOur computed solution solves the linear system (Ax == b): ', \n",
" np.allclose(A@x, b))\n",
"\n",
"# and check further that this agrees with the solution we computed by hand\n",
"print('\\nand agrees with the exact solution we computed by hand: ', \n",
" np.allclose(x, [37./11., 1./11.]))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## The determinant \n",
"\n",
"So for *square* linear systems (and also when finding eigenvalues/vectors), we need a method to determine whether a given matrix is invertible or not - the determinant of the matrix tells us this.\n",
"\n",
"The [*determinant*](https://en.wikipedia.org/wiki/Determinant) is a quantity of a square matrix which can be computed directly from its entries.\n",
"\n",
"It can be written as either $\\det(A)$ or $|A|$.\n",
"\n",
"Geometrically it gives the scaling of the volume (area in 2D) of a shape under the linear transformation represented by the matrix.\n",
"\n",
"If the mapped volume has zero volume, how does this fact alone tell you that the matrix can't have an inverse?\n",
"\n",
"[since zero volume in the mapping implies that (at least) two distinct points are mapped to the same point - we therefore do not have a well-defined inverse]."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Consider the matrix\n",
"\n",
"$$\n",
"\\left(\n",
" \\begin{array}{rr}\n",
" 2 & 3 \\\\\n",
" 4 & 6 \\\\\n",
" \\end{array}\n",
"\\right),\n",
"$$\n",
"\n",
"for a $2 \\times 2$ matrix we can compute the determinant in our heads (or with a calculator!) - we multiply the two main diagonal entries and subtract the multiple of the off-diagonal entries:\n",
"\n",
"$$\\det(A) = (2\\times 6) - (3\\times 4) = 0$$\n",
"\n",
"Let's check this using a SciPy function for computing determinants:"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {
"scrolled": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"-6.661338147750939e-16\n",
"True\n"
]
},
{
"data": {
"text/plain": [
"False"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"A = np.array( [[2, 3], [4, 6] ] )\n",
"print(sl.det(A))\n",
"print(np.allclose(0,sl.det(A)))\n",
"sl.det(A)==0"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {
"scrolled": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"-11.0\n"
]
}
],
"source": [
"# and for our linear system example above:\n",
"A = np.array([[2, 3], [1, -4]])\n",
"print(sl.det(A))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The former determinant is zero (to round off) and this means that the corresponding matrix does not have an inverse, \n",
"\n",
"the latter is safely non-zero and hence is why we were able to find the inverse without problem."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Properties of the determinant\n",
"\n",
"1. For square, equal size matrices $|AB| = |A|\\,|B|$\n",
"\n",
"\n",
"2. $|A^T| = |A|$\n",
"\n",
"\n",
"3. $|\\alpha A| = \\alpha^n |A|$, where $n$ is the dimension of the square matrix $A$\n",
"\n",
"\n",
"4. If any row (or column) is made up of zero entries only, then the determinant of that matrix is zero\n",
"\n",
"\n",
"5. The determinant is equal to the product of the eigenvalues of the matrix (where we need to multiply multiple times if an eigenvalue has an algebraic multiplicity greater than one) - so a determinant of zero is equivalent to having at least one zero eigenvalue!\n",
"\n",
"\n",
"6. For an upper or lower triangular matrix, the determinant is equal to the product of the main diagonal entries (as the eigenvalues in this case are the diagonal entries!)\n",
"\n",
"\n",
" \n",
"\n",
"Exercise: you can check that all of these seem to be true through considering some concrete examples."
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {
"scrolled": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"det(A) = -6.661338147750939e-16\n",
"\n",
"Marix A is singular: True\n"
]
}
],
"source": [
"# Let's check our maths!\n",
"A = np.array([[2., 3.], [4., 6.]])\n",
"print('det(A) = ', sl.det(A))\n",
"print('\\nMarix A is singular: ', np.isclose(sl.det(A), 0.))"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {
"scrolled": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"A = \n",
"[[0 1 0]\n",
" [1 2 3]\n",
" [4 5 6]]\n"
]
}
],
"source": [
"# can you do this one in your head?\n",
"A = np.array([[0, 1, 0], [1, 2, 3], [4, 5, 6]])\n",
"print('A = ')\n",
"print(A)"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {
"scrolled": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"det(A) = 6.0\n"
]
}
],
"source": [
"print('det(A) = ', sl.det(A))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### The inverse matrix\n",
"\n",
"The existence of an inverse to a square matrix is equivalent to whether or not that matrix's determinant is zero or not.\n",
"\n",
"If the determinant is non-zero then it does have an inverse.\n",
"\n",
"For a $2\\times 2$ matrix an expression for the inverse of a matrix\n",
"\n",
"$$A = \n",
"\\begin{pmatrix}\n",
" a & b\\\\\n",
" c & d\n",
"\\end{pmatrix}$$\n",
"\n",
"is\n",
"\n",
"$$A^{-1} = \\frac{1}{|A|}\n",
"\\begin{pmatrix}\n",
" d & -b\\\\\n",
" -c & a\n",
"\\end{pmatrix}$$\n",
"\n",
"How is the inverse existing/not-existing depending on whether the determinant is zero or not consistent with this expression?\n",
"\n",
"Similar expressions exist for larger square matrices, but they are harder to evaluate."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Using SciPy to compute the inverse matrix\n",
"\n",
"A quick review of some useful SciPy/NumPy functions.\n",
"\n",
"Let's try a $3\\times 3$ example."
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {
"scrolled": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[[ 0.142857 -0.075188 0.022556]\n",
" [-0.285714 0.518797 -0.255639]\n",
" [ 0.142857 -0.285714 0.285714]]\n",
"\n",
" [[ 0.142857 -0.075188 0.022556]\n",
" [-0.285714 0.518797 -0.255639]\n",
" [ 0.142857 -0.285714 0.285714]]\n",
"133.00000000000003\n",
"\n",
" 133.0\n",
"array([[ 1.0000e+00, 0.0000e+00, 0.0000e+00],\n",
" [-2.2204e-16, 1.0000e+00, 2.2204e-16],\n",
" [-8.3267e-17, -1.1102e-16, 1.0000e+00]])\n",
"\n",
"\n",
"array([[ 1.0000e+00, 0.0000e+00, 0.0000e+00],\n",
" [-2.2204e-16, 1.0000e+00, 2.2204e-16],\n",
" [-8.3267e-17, -1.1102e-16, 1.0000e+00]])\n",
"\n",
" True\n",
"array([[ 1.0000e+00, 0.0000e+00, 0.0000e+00],\n",
" [-2.2204e-16, 1.0000e+00, 2.2204e-16],\n",
" [-8.3267e-17, -1.1102e-16, 1.0000e+00]])\n",
"\n",
"\n",
"array([[ 1.0000e+00, 0.0000e+00, 0.0000e+00],\n",
" [-2.2204e-16, 1.0000e+00, 2.2204e-16],\n",
" [-8.3267e-17, -1.1102e-16, 1.0000e+00]])\n",
"\n",
"\n",
"array([[ 1.0000e+00, 0.0000e+00, 0.0000e+00],\n",
" [-2.2204e-16, 1.0000e+00, 2.2204e-16],\n",
" [-8.3267e-17, -1.1102e-16, 1.0000e+00]])\n",
"\n",
"\n",
"[[ 1.4286 -0.1504 0.0226]\n",
" [-1.7143 2.594 -1.0226]\n",
" [ 0.1429 -1.1429 2. ]]\n"
]
}
],
"source": [
"A = np.array([[10., 2., 1.],\n",
" [6., 5., 4.],\n",
" [1., 4., 7.]])\n",
"\n",
"# the inverse of the matrix A - computed using a scipy algorithm\n",
"print(sl.inv(A))\n",
"\n",
"# or via numpy\n",
"print('\\n', np.linalg.inv(A))\n",
"\n",
"# see https://docs.scipy.org/doc/scipy/reference/tutorial/linalg.html \n",
"# for comments on scipy.linalg vs numpy.linalg\n",
"#sl.inv?\n",
"\n",
"# the determinant of the matrix A\n",
"print(sl.det(A))\n",
"\n",
"# OR\n",
"print('\\n', np.linalg.det(A))\n",
"\n",
"# it sometimes makes our life a bit easier if we don't always print out all \n",
"# significant figures. We can do this in notebooks with \n",
"%precision 4\n",
"\n",
"# Multiply A with its inverse using numpy's dot. Note that due to\n",
"# roundoff errors the off diagonal values are not exactly zero.\n",
"pprint(np.dot(A,sl.inv(A)))\n",
"\n",
"# print a blank line just to make output easier to read\n",
"print('\\n')\n",
"\n",
"# check the result with the @ operator as well\n",
"pprint(A @ sl.inv(A))\n",
"\n",
"# use allclose to confirm this is indeed the identity (to round-off)\n",
"print('\\n', np.allclose(A @ sl.inv(A), np.eye(np.shape(A)[0])))\n",
"\n",
"\n",
"# three ways to achieve the same thing\n",
"pprint(np.dot(A,sl.inv(A)))\n",
"print('\\n')\n",
"pprint(A.dot(sl.inv(A)))\n",
"print('\\n')\n",
"pprint(A @ sl.inv(A))\n",
"print('\\n')\n",
"\n",
"# and how not to do it.\n",
"print(A*sl.inv(A))\n",
"# note that the * operator simply does operations \"element-wise\"\n",
"# which is not the same thing as matrix-matrix multiplication\n",
"# (so * in numpy is the same as \".*\" in matlab)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can ignore the very small numbers which are due to floating point round-off errors - `np.allclose` can be used to confirm that these are indeed the identity matrix."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Singular matrices\n",
"\n",
"If $|A|=0$ then $A$ is termed a [*singular matrix*](http://mathworld.wolfram.com/SingularMatrix.html) and it does not have an inverse.\n"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {
"scrolled": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Marix A is singular: True\n"
]
}
],
"source": [
"# example from above\n",
"A = np.array([[2., 3.], [4., 6.]])\n",
"print('Marix A is singular: ', np.isclose(sl.det(A), 0.))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"What happens if we try to use SciPy to compute the inverse of a matrix with a zero (or very small) determinant?"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {
"scrolled": true
},
"outputs": [
{
"ename": "LinAlgError",
"evalue": "singular matrix",
"output_type": "error",
"traceback": [
"\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
"\u001b[1;31mLinAlgError\u001b[0m Traceback (most recent call last)",
"\u001b[1;32m\u001b[0m in \u001b[0;36m\u001b[1;34m\u001b[0m\n\u001b[1;32m----> 1\u001b[1;33m \u001b[0msl\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0minv\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mA\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m",
"\u001b[1;32m~\\anaconda3\\lib\\site-packages\\scipy\\linalg\\basic.py\u001b[0m in \u001b[0;36minv\u001b[1;34m(a, overwrite_a, check_finite)\u001b[0m\n\u001b[0;32m 966\u001b[0m \u001b[0minv_a\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0minfo\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mgetri\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mlu\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mpiv\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mlwork\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mlwork\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0moverwrite_lu\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 967\u001b[0m \u001b[1;32mif\u001b[0m \u001b[0minfo\u001b[0m \u001b[1;33m>\u001b[0m \u001b[1;36m0\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 968\u001b[1;33m \u001b[1;32mraise\u001b[0m \u001b[0mLinAlgError\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34m\"singular matrix\"\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 969\u001b[0m \u001b[1;32mif\u001b[0m \u001b[0minfo\u001b[0m \u001b[1;33m<\u001b[0m \u001b[1;36m0\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 970\u001b[0m raise ValueError('illegal value in %d-th argument of internal '\n",
"\u001b[1;31mLinAlgError\u001b[0m: singular matrix"
]
}
],
"source": [
"sl.inv(A)"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {
"scrolled": false
},
"outputs": [
{
"ename": "LinAlgError",
"evalue": "Singular matrix",
"output_type": "error",
"traceback": [
"\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
"\u001b[1;31mLinAlgError\u001b[0m Traceback (most recent call last)",
"\u001b[1;32m\u001b[0m in \u001b[0;36m\u001b[1;34m\u001b[0m\n\u001b[1;32m----> 1\u001b[1;33m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mlinalg\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0minv\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mA\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m",
"\u001b[1;32m<__array_function__ internals>\u001b[0m in \u001b[0;36minv\u001b[1;34m(*args, **kwargs)\u001b[0m\n",
"\u001b[1;32m~\\anaconda3\\lib\\site-packages\\numpy\\linalg\\linalg.py\u001b[0m in \u001b[0;36minv\u001b[1;34m(a)\u001b[0m\n\u001b[0;32m 543\u001b[0m \u001b[0msignature\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;34m'D->D'\u001b[0m \u001b[1;32mif\u001b[0m \u001b[0misComplexType\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mt\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;32melse\u001b[0m \u001b[1;34m'd->d'\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 544\u001b[0m \u001b[0mextobj\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mget_linalg_error_extobj\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0m_raise_linalgerror_singular\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 545\u001b[1;33m \u001b[0mainv\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0m_umath_linalg\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0minv\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0ma\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0msignature\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0msignature\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mextobj\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mextobj\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 546\u001b[0m \u001b[1;32mreturn\u001b[0m \u001b[0mwrap\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mainv\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mastype\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mresult_t\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mcopy\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;32mFalse\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 547\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n",
"\u001b[1;32m~\\anaconda3\\lib\\site-packages\\numpy\\linalg\\linalg.py\u001b[0m in \u001b[0;36m_raise_linalgerror_singular\u001b[1;34m(err, flag)\u001b[0m\n\u001b[0;32m 86\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 87\u001b[0m \u001b[1;32mdef\u001b[0m \u001b[0m_raise_linalgerror_singular\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0merr\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mflag\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 88\u001b[1;33m \u001b[1;32mraise\u001b[0m \u001b[0mLinAlgError\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34m\"Singular matrix\"\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 89\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 90\u001b[0m \u001b[1;32mdef\u001b[0m \u001b[0m_raise_linalgerror_nonposdef\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0merr\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mflag\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
"\u001b[1;31mLinAlgError\u001b[0m: Singular matrix"
]
}
],
"source": [
"np.linalg.inv(A)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Remember - read the error message!"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Summary\n",
"\n",
"\n",
"- We have discussed more examples of inversion and closely related optimisation problems \n",
"\n",
"\n",
"- We introduced a simple tomography example as a demonstration for where an under-determined problem might arise\n",
"\n",
"\n",
"- We have reviewed some more linear algebra theory, mainly starting to think about the solvability of square systems - the aim here was to generate some insight that will be useful when considering non-square cases next.\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"scrolled": false
},
"outputs": [],
"source": []
}
],
"metadata": {
"celltoolbar": "Slideshow",
"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.13"
},
"latex_envs": {
"LaTeX_envs_menu_present": true,
"autoclose": true,
"autocomplete": true,
"bibliofile": "biblio.bib",
"cite_by": "apalike",
"current_citInitial": 1,
"eqLabelWithNumbers": true,
"eqNumInitial": 1,
"hotkeys": {
"equation": "Ctrl-E",
"itemize": "Ctrl-I"
},
"labels_anchors": false,
"latex_user_defs": false,
"report_style_numbering": false,
"user_envs_cfg": true
},
"toc": {
"base_numbering": 1,
"nav_menu": {},
"number_sections": true,
"sideBar": true,
"skip_h1_title": false,
"title_cell": "Table of Contents",
"title_sidebar": "Contents",
"toc_cell": true,
"toc_position": {},
"toc_section_display": true,
"toc_window_display": false
}
},
"nbformat": 4,
"nbformat_minor": 1
}