{ "cells": [ { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "using LinearAlgebra, PyPlot" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Review: Solving ODEs via eigenvectors\n", "\n", "If we have a simple scalar ODE:\n", "\n", "$$\n", "\\frac{dx}{dt} = a x\n", "$$\n", "\n", "then the solution is\n", "\n", "$$\n", "x(t) = e^{at} x(0)\n", "$$\n", "\n", "where $x(0)$ is the initial condition.\n", "\n", "If we have an $m\\times m$ system of ODEs\n", "\n", "$$\n", "\\frac{d\\vec{x}}{dt} = A\\vec{x}\n", "$$\n", "\n", "we know that if $A = X \\Lambda X^{-1}$ is diagonalizable with eigensolutions $A\\vec{x}_k = \\lambda_k \\vec{x}_k$ ($k=1,2,\\ldots,m$), then we can write the solution as:\n", "\n", "$$\n", "\\vec{x}(t) = c_1 e^{\\lambda_1 t} \\vec{x}_1 + c_2 e^{\\lambda_2 t} \\vec{x}_2 + \\cdots\n", "$$\n", "\n", "where the $\\vec{c}$ coefficients are determined from the initial conditions\n", "\n", "$$\n", "\\vec{x}(0) = c_1 \\vec{x}_1 + c_2 \\vec{x}_2 + \\cdots\n", "$$\n", "\n", "i.e. $\\vec{c} = X^{-1} \\vec{x}(0)$ where $X$ is the matrix whose columns are the eigenvectors and $\\vec{c} = (c_1, c_2, \\ldots, c_m)$.\n", "\n", "## Matrix exponential, first guess:\n", "\n", "It sure would be nice to have a formula as simple as $e^{at} x(0)$ from the scalar case. Can we **define the exponential of a matrix** so that \n", "\n", "$$\n", "\\vec{x}(t) = \\underbrace{e^{At}}_\\mbox{???} \\vec{x}(0) \\, ?\n", "$$\n", "\n", "But what is the exponential of a matrix? \n", "\n", "We can guess at least one case. For **eigenvectors, the matrix A acts like a scalar λ**, so we should have $e^{At} \\vec{x}_k = e^{\\lambda_k t} \\vec{x}_k$! \n", "\n", "This turns out to be exactly correct, but let's take it a bit more slowly." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Writing ODE solution in matrix form\n", "\n", "Another way of saying this is that we'd like to write the solution $x(t)$ as $\\mbox{(some matrix)} \\times \\vec{x}(0)$. This will help us to understand the solution as a *linear operation on the initial condition* and manipulate it algebraically, in much the same way as writing the solution to $Ax=b$ as $x = A^{-1} b$ helps us work with matrix equations (even though we rarely compute matrix inverses explicitly in practice).\n", "\n", "To do so, let's break down\n", "\n", "$$\n", "\\vec{x}(t) = c_1 e^{\\lambda_1 t} \\vec{x}_1 + c_2 e^{\\lambda_2 t} \\vec{x}_2 + \\cdots\n", "$$\n", "\n", "into steps.\n", "\n", "1. Compute $\\vec{c} = X^{-1} \\vec{x}(0)$. That is, write the initial condition in the basis of eigenvectors. (In practice, we would solve $X \\vec{c} = \\vec{x}(0)$ by elimination, rather than computing $X^{-1}$ explicitly!)\n", "\n", "2. Multiply each component of $\\vec{c}$ by $e^{\\lambda t}$.\n", "\n", "3. Multiply by $X$: i.e. multiply each coefficient $c_k e^{\\lambda_k t}$ by $\\vec{x}_k$ and add them up.\n", "\n", "In matrix form, this becomes:\n", "\n", "$$\n", "\\vec{x}(t) = X \\underbrace{\\begin{pmatrix} e^{\\lambda_1 t} & & & \\\\\n", " & e^{\\lambda_2 t} & & \\\\\n", " & & \\ddots & \\\\\n", " & & & e^{\\lambda_m t} \\end{pmatrix}}_{e^{\\Lambda t}} \\underbrace{X^{-1} \\vec{x}(0)}_\\vec{c}\n", "= \\boxed{ e^{At} \\vec{x}(0) }\n", "$$\n", "\n", "where we have *defined* the \"matrix exponential\" of a diagonalizable matrix as:\n", "\n", "$$\n", "e^{At} = X e^{\\Lambda t} X^{-1}\n", "$$\n", "\n", "Note that we have defined the exponential $e^{\\Lambda t}$ of a *diagonal matrix* $\\Lambda$ to be the diagonal matrix of the $e^{\\lambda t}$ values.\n", "\n", "* Equivalently, $e^{At}$ is the matrix with the **same eigenvectors as A but with eigenvalues λ replaced by** $e^{\\lambda t}$.\n", "\n", "* Equivalently, **for eigenvectors, A acts like a number λ**, so $e^{At} \\vec{x}_k = e^{\\lambda_k t} \\vec{x}_k$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Example\n", "\n", "For example, the matrix\n", "\n", "$$\n", "A = \\begin{pmatrix}\n", "0 & 1 \\\\\n", "1 & 0\n", "\\end{pmatrix}\n", "$$\n", "\n", "has two eigenvalues $\\lambda_1 = +1$ and $\\lambda_2 = -1$ (corresponding to exponentially *growing* and *decaying* solutions to $d\\vec{x}/dt = A\\vec{x}$, respectively). The corresponding eigenvectors are:\n", "\n", "$$\n", "\\vec{x}_1 = \\begin{pmatrix} 1 \\\\ 1 \\end{pmatrix} , \\; \\vec{x}_2 = \\begin{pmatrix} 1 \\\\ -1 \\end{pmatrix} .\n", "$$\n", "\n", "Hence, the matrix exponential should be:\n", "\n", "$$\n", "e^{At} = \\underbrace{\\begin{pmatrix} 1 & 1 \\\\ 1 & -1 \\end{pmatrix}}_X\n", " \\underbrace{\\begin{pmatrix} e^t & \\\\ & e^{-t} \\end{pmatrix}}_{e^{\\Lambda t}}\n", " \\underbrace{\\begin{pmatrix} 1 & 1 \\\\ 1 & -1 \\end{pmatrix}^{-1}}_{X^{-1}}\n", " = \\begin{pmatrix} 1 & 1 \\\\ 1 & -1 \\end{pmatrix}\n", " \\begin{pmatrix} e^t & \\\\ & e^{-t} \\end{pmatrix}\n", " \\left[ \\frac{1}{2} \\begin{pmatrix} 1 & 1 \\\\ 1 & -1\\end{pmatrix} \\right]\n", " = \\frac{1}{2} \n", " \\begin{pmatrix} e^t & e^{-t} \\\\ e^t & -e^{-t} \\end{pmatrix}\n", " \\begin{pmatrix} 1 & 1 \\\\ 1 & -1 \\end{pmatrix}\n", " = \\frac{1}{2} \n", " \\begin{pmatrix} e^t + e^{-t} & e^t - e^{-t} \\\\ e^t - e^{-t} & e^t + e^{-t}\\end{pmatrix}\n", " = \\begin{pmatrix} \\cosh(t) & \\sinh(t) \\\\ \\sinh(t) & \\cosh(t) \\end{pmatrix}\n", "$$\n", "\n", "In this example, $e^{At}$ turns out to have a very nice form! In general, no one ever, ever, calculates matrix exponentials analytically like this except for toy $2\\times 2$ problems or *very* special matrices. (I will never ask you to go through this tedious algebra on an exam.)\n", "\n", "The computer is pretty good at computing matrix exponentials, however, and in Julia this is calculated by the `exp(A*t)` function. (There is a famous paper: [19 dubious ways to compute the exponential of a matrix](http://www.cs.cornell.edu/cv/researchpdf/19ways+.pdf) on techniques for this tricky problem.) Let's try it:" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "2×2 Matrix{Float64}:\n", " 1.54308 1.1752\n", " 1.1752 1.54308" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "t = 1\n", "[cosh(t) sinh(t)\n", " sinh(t) cosh(t)]" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "2×2 Matrix{Float64}:\n", " 1.54308 1.1752\n", " 1.1752 1.54308" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "exp([0 1; 1 0]*t)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Yup, it matches for $t=1$.\n", "\n", "What happens for larger $t$, say $t=20$?" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "2×2 Matrix{Float64}:\n", " 2.42583e8 2.42583e8\n", " 2.42583e8 2.42583e8" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "t = 20\n", "[cosh(t) sinh(t); sinh(t) cosh(t)]" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "2×2 Matrix{Float64}:\n", " 2.42583e8 2.42583e8\n", " 2.42583e8 2.42583e8" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "exp([0 1; 1 0]*20)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For large $t$, the $e^t$ exponentially growing term takes over, and $\\cosh(t) \\approx \\sinh(t) \\approx e^t/2$:\n", "\n", "$$\n", "e^{At} = \\begin{pmatrix} \\cosh(t) & \\sinh(t) \\\\ \\sinh(t) & \\cosh(t) \\end{pmatrix}\n", "\\approx \\frac{e^t}{2} \\begin{pmatrix} 1 & 1 \\\\ 1 & 1 \\end{pmatrix}\n", "$$" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "2×2 Matrix{Float64}:\n", " 2.42583e8 2.42583e8\n", " 2.42583e8 2.42583e8" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "exp(20)/2 * [1 1; 1 1]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "But we could have seen this from our eigenvector expansion too:\n", "\n", "$$\n", "\\vec{x}(t) = c_1 e^t \\begin{pmatrix} 1 \\\\ 1 \\end{pmatrix} + c_2 e^{-t} \\begin{pmatrix} 1 \\\\ -1 \\end{pmatrix} \\approx c_1 e^t \\begin{pmatrix} 1 \\\\ 1 \\end{pmatrix}\n", "$$\n", "\n", "where $c_1$ is the coefficient of the initial condition: (nearly) every initial condition should give $\\vec{x}(t)$ proportional to $(1,1)$ for large $t$, except in the very special case where $c_1 = 0$.\n", "\n", "In fact, since **these** two eigenvectors are an **orthogonal basis** (not by chance: we will see later that it happens because $A^T = A$), we can get $c_1$ just by a dot product:\n", "\n", "$$\n", "c_1 = \\frac{\\vec{x}_1 ^T \\vec{x}(0)}{\\vec{x}_1 ^T \\vec{x}_1} = \\frac{\\vec{x}_1 ^T \\vec{x}(0)}{2}\n", "$$\n", "\n", "and hence\n", "\n", "$$\n", "\\vec{x}(t) \\approx c_1 e^t \\vec{x}_1 = \\frac{e^t}{2} \\vec{x}_1 \\vec{x}_1^T \\vec{x}(0) = \\frac{e^t}{2} \\begin{pmatrix} 1 & 1 \\\\ 1 & 1 \\end{pmatrix} \\vec{x}(0) \n", "$$\n", "\n", "which is the same as our approximation for $e^{At}$ above." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Series definition of a matrix exponential\n", "\n", "Just plugging in $t=1$ above, we see that we have defined the [matrix exponential](https://en.wikipedia.org/wiki/Matrix_exponential) by\n", "\n", "$$\n", "e^{A} = X e^{\\Lambda} X^{-1}\n", "$$\n", "\n", "This works (for a diagonalizable matrix $A$, at least), but it is a bit odd. It doesn't *look* much like any definition of $e^x$ for scalar $x$, and it's not clear how you would extend it to non-diagonalizable (defective) matrices.\n", "\n", "Instead, we can **equivalently** define matrix exponentials by starting with the **Taylor series** of $e^x$:\n", "\n", "$$\n", "e^x = 1 + x + \\frac{x^2}{2!} + \\frac{x^3}{3!} + \\cdots + \\frac{x^n}{n!} + \\cdots\n", "$$\n", "\n", "It is quite natural to define $e^A$ (for **any square** matrix $A$) by the **same series**:\n", "\n", "$$\n", "e^A = I + A + \\frac{A^2}{2!} + \\frac{A^3}{3!} + \\cdots + \\frac{A^n}{n!} + \\cdots\n", "$$\n", "\n", "This involves only familiar matrix multiplication and addition, so it is completely unambiguous, and it converges because the $n!$ denominator grows faster than $A^n \\sim \\lambda^n$ for the biggest $|\\lambda|$.\n", "\n", "Let's try summing up 100 terms of this series for a random $A$ and comparing it to both Julia's `expm` and to our formula in terms of eigenvectors:" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "5×5 Matrix{Float64}:\n", " -0.818911 -0.837843 -0.833094 -0.465218 -0.111137\n", " 0.255629 1.16573 -1.20799 0.466481 -0.295185\n", " -1.69869 -0.807249 -1.03416 -1.17902 0.329974\n", " 1.62197 0.379368 0.15949 -0.0619288 -0.998135\n", " -0.654186 -0.666233 1.09477 -0.164912 0.640717" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "A = randn(5,5)" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "5×5 Matrix{Float64}:\n", " 0.181729 -1.11919 0.148658 -0.412992 0.246856\n", " 3.02992 4.29519 -2.92774 1.88224 -2.13708\n", " -2.20443 -1.05848 1.66506 -1.03794 1.36019\n", " 2.09946 0.722445 -1.26534 1.2118 -1.76688\n", " -2.77878 -2.14487 2.53331 -1.26725 3.12757" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "exp(A)" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "5×5 Matrix{Float64}:\n", " 0.181729 -1.11919 0.148658 -0.412992 0.246856\n", " 3.02992 4.29519 -2.92774 1.88224 -2.13708\n", " -2.20443 -1.05848 1.66506 -1.03794 1.36019\n", " 2.09946 0.722445 -1.26534 1.2118 -1.76688\n", " -2.77878 -2.14487 2.53331 -1.26725 3.12757" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "series = I + A # first two terms\n", "term = A\n", "for n = 2:100\n", " term = term*A / n # compute Aⁿ / n! from the previous term Aⁿ⁻¹/(n-1)!\n", " series = series + term\n", "end\n", "series" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "5×5 Matrix{Float64}:\n", " 0.181729 -1.11919 0.148658 -0.412992 0.246856\n", " 3.02992 4.29519 -2.92774 1.88224 -2.13708\n", " -2.20443 -1.05848 1.66506 -1.03794 1.36019\n", " 2.09946 0.722445 -1.26534 1.2118 -1.76688\n", " -2.77878 -2.14487 2.53331 -1.26725 3.12757" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "λ, X = eigen(A)\n", "X * Diagonal(exp.(λ)) / X" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "5×5 Matrix{Float64}:\n", " 0.181729 -1.11919 0.148658 -0.412992 0.246856\n", " 3.02992 4.29519 -2.92774 1.88224 -2.13708\n", " -2.20443 -1.05848 1.66506 -1.03794 1.36019\n", " 2.09946 0.722445 -1.26534 1.2118 -1.76688\n", " -2.77878 -2.14487 2.53331 -1.26725 3.12757" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "real(X * Diagonal(exp.(λ)) / X) # get rid of tiny imaginary parts" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Hurray, they all match, up to roundoff errors! (Though the eigenvector method doesn't realize that the result is real, and we see tiny imaginary parts due to roundoff errors.)\n", "\n", "But why does the eigenvector definition match the series definition? They look quite different, but they are not! We can see this simply by looking at what the series does to an eigenvector:\n", "\n", "## Series definition for eigenvectors\n", "\n", "Even simpler, the key fact is that the eigenvalues of $e^A$ are $e^\\lambda$. We can see this from the series definition:\n", "\n", "If $Ax = \\lambda x$, thend\n", "$$\n", "e^A x = \\left(I + A + \\frac{A^2}{2!} + \\cdots\\right)x = \\left(1 + \\lambda + \\frac{\\lambda^2}{2!} + \\cdots\\right) x = e^\\lambda x\n", "$$\n", "from the series definition of $e^\\lambda$.\n", "\n", "It follows that $e^A$ has the same eigenvectors as $A$ and the eigenvalues become $e^\\lambda$.\n", "\n", "If $A$ is diagonalizable, this means $e^A = X e^\\Lambda X^{-1}$: we get the same result as before!" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Matrix exponentials and derivatives\n", "\n", "In first-year calculus, we learn that $\\frac{d}{dt} e^{at} = a e^{at}$. The same thing works for matrices!\n", "\n", "$$\n", "\\boxed{\\frac{d}{dt} e^{At} = A e^{At}}\n", "$$\n", "\n", "You can derive this in various ways. For example, you can plug $e^{At}$ into the series definition and take the derivative term-by-term.\n", "\n", "This is why $\\vec{x}(t) = e^{At} \\vec{x}(0)$ solves our ODE:\n", "\n", "1. It satisfies $d\\vec{x}/dt = A\\vec{x}$, since $\\frac{d}{dt} e^{At} \\vec{x}(0) = A e^{At} \\vec{x}(0)$\n", "\n", "2. It satisfies the initial condition: $e^{A\\times0} \\vec{x}(0) = \\vec{x}(0)$, since from the series definition we can see that $e^{A\\times0}=I$.\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Products of matrix exponentials\n", "\n", "In high school, you learn that $e^x e^y = e^{x+y}$. (In fact, exponentials $a^x$ are essentially the *only* functions that have this property.)\n", "\n", "However, this is **not** in general true for matrices:\n", "\n", "$$\n", "\\boxed{e^A e^B \\ne e^{A + B} }\n", "$$\n", "\n", "unless $AB = BA$ (unless they **commute**).\n", "\n", "This can be seen from the series definition: if you multiply together the series for $e^A$ and $e^B$, you can only re-arrange this into the series for $e^{A + B}$ if you are allowed to re-order products of $A$ and $B$. For example, the $(A+B)^2=(A+B)(A+B)$ term gives $A^2 +AB+BA +B^2$ (not $A^2 +2AB +B^2$!), which requires both orders $BA$ and $AB$.\n", "\n", "Let's try it:" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "5×5 Matrix{Float64}:\n", " 0.323484 -1.42524 0.0588823 -0.336001 0.390326\n", " -0.881855 4.45288 -4.68822 -1.15012 -2.35137\n", " 0.308648 -1.45824 4.01309 1.19521 1.27798\n", " -0.481631 1.64395 -4.93559 -1.08891 -1.70655\n", " 1.03601 -2.15709 9.89571 3.02326 4.39306" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "B = randn(5,5)\n", "exp(A) * exp(B)" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "5×5 Matrix{Float64}:\n", " -0.572969 -2.44162 1.64648 -0.782038 1.11299\n", " 0.681504 3.88447 -1.7922 1.01024 -1.21317\n", " -0.830851 -2.25456 1.85245 -0.607258 0.658459\n", " 1.08542 4.28899 -2.02442 1.72467 -0.937318\n", " -2.03497 -4.54824 5.62928 -1.55535 3.48544" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "exp(A + B)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "They are not even close!" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "However, since $A$ and $2A$ commute ($A\\times2A=2A^2 = 2A \\times A$), we *do* have $e^{A}e^{2A}=e^{3A}$:" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "5×5 Matrix{Float64}:\n", " -47.5951 -36.0198 39.0739 -23.0627 36.4263\n", " 226.625 155.117 -183.007 105.17 -174.234\n", " -88.0121 -54.1272 70.0116 -39.1407 68.523\n", " 85.1844 51.5045 -67.7549 37.6019 -67.0832\n", " -174.238 -113.934 140.131 -79.2429 135.861" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "exp(A) * exp(2A)" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "5×5 Matrix{Float64}:\n", " -47.5951 -36.0198 39.0739 -23.0627 36.4263\n", " 226.625 155.117 -183.007 105.17 -174.234\n", " -88.0121 -54.1272 70.0116 -39.1407 68.523\n", " 85.1844 51.5045 -67.7549 37.6019 -67.0832\n", " -174.238 -113.934 140.131 -79.2429 135.861" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "exp(3A)" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "5×5 Matrix{Float64}:\n", " -47.5951 -36.0198 39.0739 -23.0627 36.4263\n", " 226.625 155.117 -183.007 105.17 -174.234\n", " -88.0121 -54.1272 70.0116 -39.1407 68.523\n", " 85.1844 51.5045 -67.7549 37.6019 -67.0832\n", " -174.238 -113.934 140.131 -79.2429 135.861" ] }, "execution_count": 16, "metadata": {}, "output_type": "execute_result" } ], "source": [ "exp(2A) * exp(A)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Inverses of matrix exponentials\n", "\n", "As a special case of the above, since $A$ and $-A$ commute, we have $e^A e^{-A} = e^{A-A} = I$, so:\n", "\n", "$$\n", "\\boxed{\\left(e^A\\right)^{-1} = e^{-A}}\n", "$$\n", "\n", "For example" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "5×5 Matrix{Float64}:\n", " 3.2008 1.52405 2.99 1.83183 0.523262\n", " 1.20707 0.966905 1.95747 0.701775 0.110564\n", " 3.10572 1.80607 4.97103 3.13859 0.60015\n", " -3.26614 -1.27767 -2.92433 -0.296963 0.488793\n", " -0.167361 0.0365872 -1.21242 -0.553745 0.572403" ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" } ], "source": [ "inv(exp(A))" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "5×5 Matrix{Float64}:\n", " 3.2008 1.52405 2.99 1.83183 0.523262\n", " 1.20707 0.966905 1.95747 0.701775 0.110564\n", " 3.10572 1.80607 4.97103 3.13859 0.60015\n", " -3.26614 -1.27767 -2.92433 -0.296963 0.488793\n", " -0.167361 0.0365872 -1.21242 -0.553745 0.572403" ] }, "execution_count": 18, "metadata": {}, "output_type": "execute_result" } ], "source": [ "exp(-A)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Matrix exponentials as propagators\n", "\n", "From above, we had $\\vec{x}(t) = e^{At} \\vec{x}(0)$ solving $d\\vec{x}/dt = A\\vec{x}$ given the initial condition at $t=0$.\n", "\n", "However, there is nothing that special about $t=0$. We could instead have given $\\vec{x}(t)$ and asked for $\\vec{x}(t+\\Delta t)$ and the result would have been similar:\n", "\n", "$$\n", "\\boxed{ \\vec{x}(t+\\Delta t) = e^{A\\Delta t} \\vec{x(t)} } = e^{A\\Delta t} e^{A t} \\vec{x}(0) = \n", "e^{A(t + \\Delta t)} \\vec{x}(0)\\, .\n", "$$\n", "\n", "Viewed in this way, the matrix $T = e^{A\\Delta t}$ can be thought of as a \"propagator\" matrix: it takes the solution at any time $t$ and \"propagates\" it forwards in time by $\\Delta t$.\n", "\n", "The *inverse* of this propagator matrix is simply $T^{-1} = e^{-A\\Delta t}$, which propagates *backwards* in time by $\\Delta t$. \n", "\n", "If we multiply by this propagator matrix repeatedly, we can get $\\vec{x}$ at a whole sequence of time points:\n", "\n", "$$\n", "\\vec{x}(0), \\vec{x}(\\Delta t), \\vec{x}(2\\Delta t), \\ldots =\n", "\\vec{x}(0), T \\vec{x}(0), T^2 \\vec{x}(0), \\ldots\n", "$$\n", "\n", "which is nice for plotting the solutions as a function of time! Let's try it for our [two masses and springs example](ODEs.ipynb):" ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjkAAAGxCAYAAACA4KdFAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8/fFQqAAAACXBIWXMAAA9hAAAPYQGoP6dpAACFHUlEQVR4nO2deXwT1fr/P5M0FAoUwbIUUmhZBNmRTVy5eFnclXtRXBECCCgKuFyXewXcisqifpUiEhHvTwWvirtCVUARwQJFkB2h0IEWLEsLlJY0md8fp6eTtkk7SWZNnvfrxevMpMnkcDJz5jPPeRZBkiQJBEEQBEEQUYbN6A4QBEEQBEFoAYkcgiAIgiCiEhI5BEEQBEFEJSRyCIIgCIKISkjkEARBEAQRlZDIIQiCIAgiKiGRQxAEQRBEVBJndAeMwufz4ciRI2jYsCEEQTC6OwRBEARBKECSJJw+fRotW7aEzVazrSZmRc6RI0eQkpJidDcIgiAIggiD3NxcOJ3OGt8TsyKnYcOGANggJSYmqnpsj8eDlStXYsiQIXA4HKoeO9qgsVIOjZVyaKyUQ2OlHBqr0NBqvIqKipCSklJxH6+JmBU5fIkqMTFRE5GTkJCAxMREuhBqgcZKOTRWyqGxUg6NlXJorEJD6/FS4mpCjscEQRAEQUQlJHIIgiAIgohKSOQQBEEQBBGVkMghCIIgCCIqsazIKSsrw7///W+kpaWhXr16aNu2LZ599ln4fD6ju0YQBEEQhAmwbHTVSy+9hAULFmDJkiXo0qULNm7ciNGjR6NRo0Z4+OGHje4eQRAEQRAGY1mR8+uvv+Lmm2/G9ddfDwBITU3Fhx9+iI0bNxrcM4IgCIIgzIBlRc4VV1yBBQsWYM+ePbjooovw+++/Y+3atXj11VcDvr+0tBSlpaUV+0VFRQBYHL/H41G1b/x4ah83GqGxUg6NlXJorJRDY6UcGqvQ0Gq8QjmeIEmSpOq364QkSXjqqafw0ksvwW63w+v14oUXXsCTTz4Z8P0zZszAzJkzq73+wQcfICEhQevuEgRBEAShAsXFxbjzzjtRWFhYazJfy4qcpUuX4rHHHsMrr7yCLl26YMuWLZgyZQrmzp2LUaNGVXt/IEtOSkoKCgoKNMl4nJmZicGDB1NWzFqgsVIOjZVyaKyUQ2OlHBqr0NBqvIqKipCUlKRI5Fh2ueqxxx7DE088gZEjRwIAunXrhoMHDyI9PT2gyImPj0d8fHy11x0Oh2Ynq5bHjjZorJRDY6UcGivl0Fgph8YqNNQer1COZVmRU1xcXK3Eut1upxByQntEEVi3jm1fdhlQSxVcgiAIwhgsK3JuvPFGvPDCC2jdujW6dOmC7OxszJ07F2PGjDG6a0Q043YD48YBfJVXEIC33wZcLmP7RRAEQVTDsiLn//7v//Cf//wHkyZNwrFjx9CyZUvcf//9eOaZZ4zuGhGtiGJlgQOw7fvvB4YOJYsOQRDhQdZhzbCsyGnYsCFeffXVoCHjBKE669ZVFjgcrxf49VdgxAj9+0QQhLUh67CmWLasA0HoitsNlDu5B2TkSPYegiAIpdRkHRZF4/oVRZDIIYjaEEVg/PjAVhyOz0cTE0EQofHzzzVbh4mIIZFDELWxdy8TMbXh9QL79mnfH4IgrI/bDdx5Z/C/k3VYFUjkEERtdOgAVElXAJut+mt2O9C+vX79IgjCmnDrcE2QdVgVSOQQRG04ncBtt8n7djuwcCH7Jwjy62+9RVERBEHUDlmHdcOy0VUEoRtuN7Bsmbyfni5HPqSlAddcAzgcwNVXG9M/giCsRYcO7AHJ3x+HW4b9xQ9ZhyOGLDkEUROBnI6ffFI2Ie/fz1qPB+jYkdbQCYKoHacTuPhieZ9bh996S37NZiPrsAqQyCGImghkVuYmZFFka+YcWkMnCEIJJSXAn3+y7TlzgJwcZh0eOxa49Vb2+k03sSSjRESQyDEjogisWgVkZbGWbprGwc3K/nATck0CiCAIIhj//jdQWsq2H3sMWLFC/tsFF7D2s8+ANm3IOhwhJHLMhtvNTuxBg4B+/VhLJ7pxOJ3ADTfI+3a7bEIOFHVFa+gEQdSEKAJz58r7/hZgUQSWLAn8NyIsSOSYCe7/UdU64POx1+lENwa7nbUPPCCblQEmdBYulP8OAK+/TmvoBEEEZ+/e6gkAuQWYrMOqQyLHTLz2WvCwQp+P/Z3Qny1bWPuPf1QXMC4XEz7cxNy3r44dIwjCcnToUP01bgEm67DqkMgxC6LIHNBqYt48suboTWEhEzEA0KNH4Pc4nUDv3mx72zZdukUQhEVp2LDyvv8SOLcOcz9AQaAIqwghkWMWApkwq0JmS/35/nvWJicDTZoEf1+3bqz95hsSogRBBOePP1jbogULLPFfAgfYNg8l79KFqpFHCIkcsxDIhFkVMlvqi9sNjBjBtvPyanb+PnWKtZ98Qo7iBEEEZ+tW1vbqBQwcGNhKc801rN2zh+XgIsKGRI5ZaNECqFdP3rfbgVGjKq/PTpyof79ilUBJAINFOYgi8N578j5FRBAEEQxeXTw1Nfh7UlOBxETg/Hk2t9BcEjYkcszCV18B586xE/v775kJ8913gYMHgcaN2XveeIOsBHoRSpQDRUQQBKEEtxv473/Z9oIFwedymw1o3pxtjx1L834EkMgxA243MHw42y4qYgLH34TJl0IAshLoRShRDhQRQRBEbVStPC5JNVuH/R+SaN4PGxI5RlPbskhNORUI7XA6meWM4x8BEei9CxfK+1RzhiCIqoRqHaZ5XxVI5BhNbSc+WQmM46qrWFu/PnDgQM1RDi6XXHPm8ccpIoIgiMqQddgQSOQYTW0nM1kJjGPvXtZ26gSkpNT+/u7dWXv8uHZ9IgjCmjidLJiEU5t1OD1d2XuJGiGRYzROZ+V12kAns7+VYNo0shLohb81TQlt27KWVxcmCILwhz8sXX999fw4VXn0UaBOHba9Zg3N+2FCIscMJCWx9rrrgp/4fOmEWxcI7eFjrdREzEXO/v3a9IcgCGvDH4CuvLJ2q4zNJs89Z85o268ohkSOGfj9d9YOGxb8xO/Th7W//EIe9nrBM5NeeKGy97drx9pDhyiBF0EQ1eEih88VtXHRRaylh9uwIZFjBjZtYm1ycvD3cCFUUEA5E/TA7QbWrWPbjzyibLxbtADi45kj+fr12vaPIAjrwUUOt/rWBl8q37NHm/7EACRyjOb114EjR9j27bcHvpmKIvDQQ/I+5UzQlqr5LJSO9zvvAKWlbHvgQBKiBEHInD4N/PUX21ZqyeEiZ/16mu/DhESOkYgiMHWqvB/sZkoZdfUlnPEOVxgRBBEbcF+9xEQmeJTALThZWWTBDxMSOUai9GZKORP0JZzxJiFKEERN8MriRUXKBIsoAnPnyvv04BQWJHKMpEMHQBAqvxboZspz5fjfeClngnY4ncDdd8v7SnJUkBAlCCIYoshqVXGUCBZ6cFIFEjlG4nQCvXrJ+zXdTF0u4Lvv2HZiIjBmjD59jFVat2btjTfWns8CICFKEERwwinTQA9OqkAix2i8Xta++GLtN9Mrr2QnfVERcPSoLt2LWQ4dYu1llykXKi4X8PHHbDslhZJ3EQTBCJRQtDbBQtnuVYFEjpFIkqzkhw+v/eStWxdIS2PbO3dq27dY5+BB1nKLjlJ692Ztfn51UzNBELGJ01l5LlFapsHlAu64g20/9BA9OIWBpUXO4cOHcffdd+PCCy9EQkICevbsiU0854wVOHoUOHuWKfTUVGWfufhi1n7+OTmgaQkXOW3ahPa55GTmZ+XxyOGiBEEQ586x9u23lS2Bc7p2Ze3Jk5p0K9qxrMg5efIkLr/8cjgcDnz77bfYsWMH5syZgwsuuMDorimHW3Fat2ZJ5JTAM+m+9hqFFGqF1ysLyFAtOQ4HSwoIkAglCIJRUiI/9Nx6a2hLTvxBiz94ESERZ3QHwuWll15CSkoKFi9eXPFaqlJriFkINcW3KAIrV8r73EN/6FBap1WTvDygrAyIiwNatgz98ykp7BiiKC9fEQQRuxw+zNq6dYEmTUL7LImciLCsyPniiy8wdOhQjBgxAmvWrEGrVq0wadIkjBs3LuD7S0tLUcqz0QIoKioCAHg8HnhUrjPEj1fbcW1ZWbAD8DZrBp+CPgg7dyIugId+2a5dkJo3D7e7hqJ0rPRE2L8fcQAkpxNlPl/IvjX2li1hA+DNyVH0uyrFjGNlVmislENjpZxwx0rIyZHnlLKy0L60ZUs4AEi5uSgrKWH+PBZBq3MrlONZVuTs378fGRkZmDZtGp566in89ttveOihhxAfH49777232vvT09Mxc+bMaq+vXLkSCQkJmvQxMzMz6N9aZ2ai55tvAgBsS5dia5MmODR4cI3Hq1tQgCGCAMFP6PhsNvxw8CBKvvlGnU4bRE1jpTepX3+NHgBO1qmDn8MY165lZWgHYP9PP2GHBtZFM42V2aGxUg6NlXJCHSvn6tXoDaCgXj2sC3VO8Xpxo90OW1kZfvx//w8lTZuG9nkToPa5VVxcrPi9giRVNQ1Ygzp16qBPnz5Yx4soAnjooYeQlZWFX3/9tdr7A1lyUlJSUFBQgMTERFX75vF4kJmZicGDB8PhcFR/gygirn17CH4WAsluR9nevbUuOwnvvAP7hAkQAEg2G7wZGZBGj1a1/3pS61jpjLB4MRtfSYIEwPvWWyGPr23uXNifeAK+O+6Ad8kS1fpmtrEyMzRWyqGxUk64Y2V7+WXY//1v+O6+G9533gn5e+M6doRw4AC8c+fCd8stlnFP0OrcKioqQlJSEgoLC2u9f1vWkpOcnIzOnTtXeu3iiy/GJ598EvD98fHxiA/g3OtwODS7sIMeOyen2hKI4PXCcfCgHCIejPvvB+bNA3bvhrBkCeL8M/NaGC1/B8WIIjBxYkXSLgFA3KRJwHXXheUoaNu+HbajR1WfkEwxVhaBxko5NFbKCXmsyosw29q0gS2cMa5TBwBgnzYN9kcfZflzLBROrva5FZLAVO1bdebyyy/H7t27K722Z88etAk15NcIIs1k2akTa8v9igiVUCuN+pYtrN26lSLgCIKQIy3DeeARRblQJ0A1rELEsiJn6tSpWL9+PV588UXs27cPH3zwARYuXIgHHnjA6K7VjtMJvPCCvK80MRSnbVvW8qq2hDqokUZdFIFXXpH3aUIiCILP1fXqhf7ZcEpCEBVYVuT07dsXy5cvx4cffoiuXbviueeew6uvvoq77rrL6K4pY8AA1rZqFVpiKEAOOech6IQ6qJFGnYrqEQThj9sN/PEH2x4zJnTLrtJCzkRALCtyAOCGG27Atm3bUFJSgp07dwYNHzclubms7dgxdBMmWXK0Y9QoeULJygp93ZuK6hEEwRFFYPx4eT8cy67TyUo6cEK1/Mc4lhY5loaLnJSU0D/LRc6ePfJxCHXIz2emYbsd6Nkz9M9TUT2CIDhqWXbHjGFtYmLolv8Yh0SOUUQiclatYm1JCat5RY6t6lEeBYHk5OoWGaW4XMAVV7DtefNoQiKIWEUty26rVqwtKgKSktTpW4xAIscoDh1ibagiRxQBf+dqcmxVF55+nU8q4cKTAPrlZiIIIsZwOllaCk64S01NmrCSEID8IEYogkSOUYRrySHHVm1RS+TwIp35+ZEdhwgdUYTwv/+h5dq1JP4J4+nShbWXXx7+UpMgyHMSn6MIRZDIMYpwRQ45tmqLWiInOZm1eXmRHYdQjigC06YBKSmIu+su9J09G3Ht2tFyLmEs/EGnW7fIfPP4Z0m4hwSJHCM4exY4eZJth1psjTu28gggQSDHVjXhpuBwqo/7Q5YcfXG7gdatmQ+UH4IkseiWrCyDOkbEPPxBhz/4hAtZcsKCRI4R+E/EXbuG/qTpcgGzZrHtgQPJsVVN1LbkkMjRHh6mG6wMn88HXHopWXQIY+BzAH/wCRc+J5ElJyRI5OiNKALTp8v74ToO8/Dmv/5SrWsE2Jo5UFErJmz4hEbLVdoTyE+tKuSgTxiFWpYcbq0nS05IkMjRG7Uch1u3Zu3Bg8GfYInQcLvlLNJ33hnZkz8XOadOsVB/QjsCZYQNBDnoE0agtiVn+3YS6yFAIkdv1ErRzUXO6dNAYaE6fYtl1MhM6s8FFwC86j0tWWnLihWVhb4gwHvHHZCqXmc2GznoE/ri88nXf6SWnM2bWbtzJxX+DQESOXrjdAKDBsn74eZNSEiQk0IdPKhe/2IVtUPzBYGcj/WgqjgFAEGA74UXsGXSJFSycUoSE0SxiCgCH33E/pEVQD9OnADKyth2s2bhH0cUZT9MgJZfQ4BEjhE0acLaBx6ILEV3mzas5YkFifDRIjSf/868OB+hPoHEqc8H4c8/caxXr8q/qSTF5o2BR57dfjv717o1WQH0gvvjJCVF5udH+dHChkSOEfAT/6qrIgv99vfLISLD6QQmTJD3Iy2C53YD2dlse/x4uqloRYcO1V+z2yG1a4cGeXkQYv3GkJUFjB1beTmPwur1g1txGzSITFxTfrSwIZFjBGp523NLDmV2VYfu3VkbSWZSoPoSSqxaEPSgaqSJnzg9k5wMKZZvDG430L9/4L9RWL0+LF3K2pycyPxoqPBv2JDI0RtJUk/k8MR1y5aRI5oaHD3K2s6dI5s8yLSsD243MGCAvP/oo5XEaUlSErwZGZWfgGPlxlBb7iCA/Dq0RhSBxYvl/UjH2+ViedEA4OWXKT+aQkjk6M3p00BxMduOROSIIvC//8n7NGFFzrFjrI3EQRAg07IeBLqJV8l2DADS6NHAhg1sRxCA227TqYMGoyR3EAB4vRB42gRCXfburS4yI33Y4YV/z58P/xgxBokcveFWnIYNgfr1wz+OFhdQrMMtOc2bR3YcblqORQuCXoRiLevTh90cJIn9DrHwIBAsd1Cg1zZu1L4/sUgQf7GIHnYoYjNkSOTojVpLVWQtUB+1LDkAMyUvW8a209LItKw2oZ7/XLg+9lhsLO1WzR1kswGLFgEvvVTtrfZ//xt1Cwp07FyM4HSypW9OpMEMAGVSDwMSOXqjZorvBQvkfTUuoFhHTZEDsLpkgFyMlVAPp5OFQ3NqOv9FEfjtN3k/2pd2A+UOAoChQ5lVqwqC14v6dNPUBoeDtS+9FFkwA4dq4oUMiRy9UUvkAMC4cbKyX76crAWRotZyFYeLpVOnaA1dC7xe1o4ZU/MNJNaWdoPkDsK+fQEtYJLdjrNqzEdEdfiD0+DB6jyA0nJVyJDI0Zs9e1jboIE6x+Nh5DyrJhEe58/LFhe1LDkXXADExbFtKqSqLpIErFnDtocNq/kGEmtLuzWVjgngL+YLZPUhIsfnk697teYUEjkhQyJHT9xuZlLn22r4BfCibVSZNjK4T4LdLmcqjhSbDWjalG3zJzpCHV56Sba8jRxZ87UUazlGnM7KxSCrLuW5XCyBaHlZGHtGBoaMGwfBP9yZiJxTp+SHTz4PRAq3uJ0+DZw9q84xoxwSOXqhVYI4LnKi1b9AL/gNs2nT6k/9kcCf4EjkqIcoAk89Je8r8bHxzzHy4ovRvbT7++/ysvhnnwVfyjt+vGJTkCTYJ02ieURN+DXfqFFkJR38adCA1S0EyJqjEBI5eqFVgjj+dEaWnMjYvp21F1yg7nFJ5KhPuD42vXuzlifRjEbcbqBXL3m/oCCwxSrAGArR7KdkBGovVQFU+DcMSOTohVZ+AbRcFTluN3DvvWx71y51w4v5BEc+OeoRbv4RHu3288/RabEIlCAxmIUriANy1PopGYHa0ZocLnK+/z46z2OVIZGjF04nkJEh76sV8k3LVZERyo0hHMgnR31atJBDcwHl1xK3UmRnR2eunFCsxdxPqdxBWRIEeOfPj14/JSPQSuSUlLB2xozoPI9VhkSOntx8s7y9b586fgH+y1U11akhAqN1nSlarlKfnTsBj4f5Jvzwg7L8I6IIpKfL+9GYKydUa7HLBcycCQAoSkmBNHiwxh2MMbQQOaLIRDonGs9jlSGRoyf8pE9KkmuQRAq35BQXAzt2qHPMWELr8GISOeqzciVru3cHBg1SZn2IhaKpTifw3HPyvhILV/mDUaNDhxDXvj1ZBdREC5ETazmfVIBEjp7wCB41T/oPPpC3u3enSSpUtA4v5r/1vn30tKUGbjcrzQCwwptKz/dYyZXTsSNrL7qodguXKFZYcgBAIKuAumghcmrKgUQEhESOnvCTXq2MulXD0mmSCg+XS3ZKXbxY3fDi9etZu3s3rZ9HSlX/qVDSMFTxQYEgRGeunK1bWXvFFbX/32LBumUkWogcpxOYPFnep3I+tUIiR0/UPulpklIPnu3Yv6BepIgiMGuWvE8iNDIiPd9dLjm/zo03RmeuHC5yunev/b2xYt0yCh7xGqjyeyTwSNALL1SnHlaUQyJHT9QWOTRJqYMkabd+TiJUPdQ433munGgtSLl5M2uV1KIqt25J5WMqAazoL1kFIsftBv78k23XlpE7VHgI+alTQMuW6h03SokKkZOeng5BEDBlyhSju1Izat9IY8UErzWnTrFoHUD99XMSoerhdFZeng3HVN+uHWv5DSiaePNN4NAhtn3HHcpurC4XynbuhDcuDgIAXHWVlj2MDbR2I2jWjM31Xm+lrNVEYCwvcrKysrBw4UJ0V2KeNRotrAUul7wk8re/kekyHPjvkpgI1K2r3nFjrWaSHpTXW8J114Vnqm/blrUnTjBxGy2IIvDQQ/J+KDfWtDSc4gkW3W5aTo0UrS24Dod8HUSrRVJFLC1yzpw5g7vuugtvv/02GjdubHR3aodHV6nleMzhTrMnTqh73FhBbYdwf1wu4OKL2fZ775EIjRTuc1Jb5fFgNGgg/87RZM2J8MZaFh/PNl5+mRzkI0UPCy6VdlBMnNEdiIQHHngA119/Pf7+97/j+eefr/G9paWlKC0trdgvKioCAHg8Hnj4UoVK8ONVPW7csWMQAJQ1aQJJze9s1gwOAJIookzl/4vWBBsrPREOH0YcAF/TpvBq0A97ixaw7dyJMq83ot/dDGNlNHGbN7NrqGnTGseyprGyt20L29GjKNu9G5IVLMBKSE1FnCBA8MuhItntKGvTRl6KDUJZTg6a/f67/ILPB+n++1GmNAdRDKHoGmzeHLYpU2CfOxcA+x288+dDat681t9CKfbmzWHbtg1loqjuvURltJqzQjmeZUXO0qVLsXnzZmRlZSl6f3p6Omb65YTgrFy5Egm8qqvKZGZmVtq/Pi8PcQBWbd+OYhWtLo6iIlwHQCgowHeffw6ff8p7i1B1rPQkddUq9ACQ7/Mh65tvVD9+n/Pn0QrAzp9/xv5GjSI+npFjZSRpX3+N7uVLKfY778SWSZNwqJYsvYHG6pL4eKQAEBcvxu5z51DCTf8W58p27dCk3HLjs9nw+4QJOLR1q2z9CkLStm24PECxzg3vv4/j3bpp1l8rU9s1mBwXh34AClu3xvpnnmHnmIpzSy+vF60B7F6zBvsuvFC142qF2nNWcXGx4vcKkmS9WgC5ubno06cPVq5ciR49egAABg4ciJ49e+LVV18N+JlAlpyUlBQUFBQgMTFR1f55PB5kZmZi8ODBcHDBcfYsHOVLap7jx4GGDdX7QklCXGIihNJSeHbvBtLS1Du2xgQcK52xPfss7M8/D+/48fC98Yb6x3/oIdgXLID36afhmz497OOYYawMQxQR1749S1hXjmS3o2zv3oDWhprGyj58OGxffcWOYbPBm5EBafRobfuvA3H9+0PIzoZ3xgz47r1XsRWmLCcHdTt2rG4FCjK2sYzSa9C2cCHsDz4I3003wfvxx6r3w/bkk7DPmQPvQw/BN3u26sdXC63mrKKiIiQlJaGwsLDW+7clLTmbNm3CsWPH0JuHgwLwer346aef8MYbb6C0tBR2u73SZ+Lj4xHP1539cDgcmt0wKh2b52FxOOA4exZo0kTdL2vVCti/H46jR1m2U4uh5e9QKwUFANiykl2LPpQX6bSfPKnK8Q0dK6PIyanmcyJ4vXAcPFijqK82VqIIfP21fAyfD3GTJjFHZivf0CWpwv/GPmIE7KE86KSmYsukSeg5fz4TOoIA4a234LDQw5Le1HoNls/3tmbNYNPiWi0v52M/dkybOUtl1J6zQjmWJR2Pr7nmGmzbtg1btmyp+NenTx/cdddd2LJlSzWBYwoWLWKtx8PqVqnt2OdfqJMIDa2qBXP4cki5mCLCQK109tFa++fYMeD0aTZGPIIsBA4NHgzfww+znVtuIQf5SOHXevkDjuqQ47FiLClyGjZsiK5du1b6V79+fVx44YXoyiONzIQeFZB5oU4SOaHDc4vEaWTY5GvmlNMifJxOOccNEH46+2jNXbR3L2tbtw47DYL0t7+xjY0bKYw8Uv76i7Va+XvxZI9//km/VS1YUuRYDj2eHrnI+fVXOulDwe1mkzoATJqkTegsWXIiR5LkFAzvvBN+OnunE5g/X96Plto/v/3G2kj+H3v2sDY3l8LII4Vf61qJnLVrWUu/Va1EjchZvXp1UKdjw9GjcmxuLms/+YROeqXoVeCURE7kHDzIlmMcDuDuuyO7md9/v7w0+cUX1l+acbuBRx9l2+vWhXXt1y0ogP1f/5JfoDprkaGlyBFFwD+AgX6rGokakWNqnE5g6FB5X+2nR1EEPvpI3qeTXhl61ZbyX66yXjCjOeBh0J07M6ETKamprD1/PvJjGUkkldn9aJCXVylyDUB0+CoZhZYih2rihQSJHL3gy0ljxqhfOTZanSm1Ri//DD7RlZQAIeR3IPz4+WfW+vvlRELr1qzl/lhWRaUb3pnk5IpCnRVEg6+SUWgpcqLVr0wjSOToBT/p+/ZVf/2fTvrwcDqBp56S97Xyz6hfH+DpC2jJKnTcboDnAlm+XJ2l2JQU1vJlXqui0rVfkpQEb0ZG5WNFg6+SERQXyw8zWogcqokXEiRy9EJLZR+tzpR6MGAAa9u3V9/CxhEEirAKl6p+U2Eux1SDixyrW3KcTsA/gWUE1740ejSwZYv8wvDhkfcvFuFzfZ066iZ99cflAgYNYtvp6db3K9MQEjl6wW9uWnnb33+/nGDwm2/opFcKz5HTrp22opBn5dy5U7vviEa08j/gy1VWt+QAAA/9rlcPOHAgsmu/Wzc5ueLmzZH3LRbxf6CtGnCiJtyvzMS1q8wAiRy94Ce+lnVG+MRdVqbdd0QbWlWG98ftBnbtYtv33EORb6Gg1VJstCxXASzyDGBCnf+/IqFPH9YuW0bBC+Ggdfg4h89ZfA4jAkIiRw98PoAX5NTyxG/ZkrVHjmj3HdGG1tmOtVpuiRWcTuDee+V9tZZiuRg4fJgtU1oZLnLatFHneNxy9vbblI4iHHbvZm2DBtp+D4kcRZDI0YNTp+SJQ0tLDomc0OEiRytLDoV7Rs4FF7D2H/9Qz2+qvEAnJIlZQKx8I+cih1tyI0EUmXM3h9JRhIbbDfDyGGHmLFIMfzDjcxgREBI5esDNlw0bMmc0raDSDqHDn4K0suRQ5Fvk8Bw511+vjt+UKAITJsj7Vr+Rc+dpNSw5JMrDp2rOIkDb84osOYogkaMHWjsdc8iSEzpaL1fxcE/ugCgIFPkWCpIkO8CqZW2Lthu5mstVJMrDR+/zikSOIkjk6IEeTseALHLIkqMcrZerALa8wkuOXHYZRb6Fwty5bLkXAG68UR3zf7TdyNUUOSTKw0fv84rPWSdOUIRVDZDI0QO9vO35chVZcpTh82lvyeFcdBFrz5zR9nuiCVEEHntM3ldrWSmabuRlZfJ4qFHuAmAi/N//ZttDh5IoV4reSfqaNGEiCiC/nBogkaMHei9XHTtGyl4JJ08yczIANG2q7Xfx356SASpHy3IlLhfwxBNs+9ZbrXsjnztXXiLp3189R9drrmHtpk3W9VUyApcL6NKFbb/7rrbnlc0mz1skcoJCIkcP9FquatqUKXtJArKztf2uaICvZTdurK1DOECVyMOhQ4fqr6lp/r/4YtYWFqpzPL0RReDJJ+V9NR2ot21j7V9/URh5qJw+zVpuvdUS8supFRI5esCjH9QyJwdj8WLZMjFgAE1MtbF9O2sbN9b+u7jApSKdynE6ZeskoH65En7svDx1jqc3Wjm6iqIcBg1YP/pMb7i1VuuHWkAWOatX0+8TBBI5WuN2Ax9+yLbnzNFOeFRNOkcTU8243cDtt7Pt/fu1F4QNGsjWIrLmKKO0VDbDL12qfm0xq0cjduhQvWyAGpauaIs+05PSUuDsWbath8gpKmLtSy+RxS0IJHK0RM9stzQxKUfvfBYAuxnRklVo7NzJHGsvuAC47Tb1HTiTk1l76hRw7py6x9YDpxP4+9/lfbUsXdEWfaYnPLO9zQY0aqTtd4kisGGDvE8PtgEhkaMhwr59+gkPmpiUY5QgpErkocGTAHbvrk2hw0aNWFFLwLpLVnyp9YEH1LN08Sgh//nEqtFnesOv7caNq8/HaqOlY34UQSJHQ6T27fUTHtEUFqs1RglCsuSExtq1rG3XTpvjC4L1l6x4Tqyrr1b3Wne5gB075PnkuuvUO3Y0wy05eixVabVcGWWQyNESpxNYsEDeV9txsiouFwspBYDLL7duWKzWVBWEWuez4JDIUY7bzQpEAiwUVytfA75kZVWRw/vNc2SpSceOQOfObHvTJvWPH41wS06TJtp/l9NZ2UFc6/uLRSGRozX//Ke8vXu39sKja1fW8icKIjAuF9CrF9tesEAfQUjLVcrQ05fNyhFWkiSLHP8oNDXp04e1//sf+XooQc/IKgAYNYq1jRur75gfJZDI0Rr+1N6ggXZmd3+oSKdyeLkA/rSqNWTJUYaePlNcHPz6q/Vu4sePs2geQDuRc/48a997j6J3lKC3yOEh5IWFslWSqASJHK3RK9sxh092hYVyKCMRGL1KOnD4xEcip2b09DXIzWXtsmXWu4nzB5mmTbVJZimKbFw4FL1TO3r65ADyfcXnIwtxEEjkaI1edas4iYlA/fps24omeL0oLpbrSGlZnNMfKu2gDKcT+Mc/5H2tfA1EEfj0U3nfajdxLnK08McBKC1FOOjpkwOwBLNcUFFph4CQyNEavUo6cPwjRmjJKjh8QoiPBxo21Oc7ablKOfHxrB0zRjtfA6uH4GotcigtRejovVwFUGmHWiCRozV6L1cB1g+L1QP/pSotcrAEgk98omgda4FR8Bw5t96qXbSI1W/iO3eyVqukc5SWInT0Xq4CSOTUAokcrdF7uQog52MlcJGj11IVAPzwA2sLCqzn/6EnpaUsRwugbXV4pxOYPVvet1IIrtsNvPoq2/7wQ+3OJZcLeO01tt22LTB0qDbfEy3k57OW1xDUAxI5NUIiR2v0Xq4CyJKjBD4h6OV0LIrA00/L+1bz/9CT9HT5JnHZZdqKwYceki0VGzZYIwS3alkSLUPsATmA4c8/SZzXhNsN7NnDtu++W79x4nMYiZyAkMjRGiOWq7glJzubbqLB0NuSQ06cyhBF4Nln5X2txaDdLt8k7HZtvkNt9DyXSJwrw8gCyXwOI8fjgJDI0RojLDm7drF29Wp68gqG3pYcq/t/6IURzsBWM/freS6ROFeGkeNktfNXZ0jkaI3elhxRlNPhA/TkFQy9c+RwJ06OXqUkrEaHDtVf01oM8psE96cwO3qWiyFxrgwjx4lETo2QyNEavR2P6clLGUY4HrtcQN++bHv+fGv4f+hNcnLlxHZ6OANb8SZx883y9r592p1LgSqSZ2SQOK+K01l5WU9PJ3Z+/h48SA+zASCRoyWSpH/eBHryUgafDPQKH+dQ6vWa2bWLlRKoVw/4/nt96vG0aMFaK4kcbnVq1gxITdX2u1wu9vDExedll2n7fVblqqtYm5qqbx2pVatYS1GbAbGsyElPT0ffvn3RsGFDNGvWDLfccgt2795tdLcqU1goR4noJXJoWaR23G5WLBUA7rlH30mBnwdUQDUwmZms7doVuOYafZ+ErShyuEDTmrZtgSuuYNuLFpHFIBD8gbZ1a/3mW1EEnnxS3if3hGpYVuSsWbMGDzzwANavX4/MzEyUlZVhyJAhOGumek18qapuXX2z3LpcQKdObHvJEloW8cfIKAiAKpHXhNsNTJvGtjdu1E98kshRBi8X8+qrZDEIhBHZjsk9oVbijO5AuHz33XeV9hcvXoxmzZph06ZNuIqbDf0oLS1FKa/YC6CoqAgA4PF44PF4VO1bxfHeeYe1JSWQ2rSBNyMD0ujRqn5XMOytW8O2axfKSkogqfz/UxM+Vmr/BsEQdu5EXIBJoWzXLkg6+OfYGjWCHYDvr7/gDfH/rPdY6YooIm78eAh+uV+k++9H2aBBYT0VhzJWwoUXIg6AlJ+PMouMre3wYXYeNWsW8nlUFUVjJYqI+/prVCzu+nwR/T5Wpaaxsv31F/tNGjeO+DdRTGoq4mw2CH5zmmS3o6xNG8AE57JWc1Yox7OsyKlKYWEhAKBJkMJo6enpmDlzZrXXV65ciYSEBNX7U7egAHFz5lTsCz4fbBMnItNuR4kOTsg9vV60AbB39Wrs0TJrrEpk8mUKjalbUIAhgiDfTAH4bDb8cPAgSr75RvPvb3PkCHoCOLprF34L8/v0Gis9Sdq2DZdXEZ+C14sN77+P4926hX1cJWOVeOAA/gagNDcXK3Q4B9Sgy/r1aA/gz7NnsUOlPtc0Vlr9PlYl0Fh13bQJ7QD8efKkar+JElpPnIieb74JAYAkCNgyYQIObd0ql0YxAWrPWcXFxYrfK0hS1aQUyvB4PMjPz0dxcTGaNm0aVFzogSRJuPnmm3Hy5En8/PPPAd8TyJKTkpKCgoICJCYmqtofj8eD7Llzcfl//lPtb2WZmZCuvlrV7wuEbfp02NPT4b3/fvj+7/80/75w8Xg8yMzMxODBg+FwOHT5Ttvjj8NenhJfstvhnT9fNwub8MkniLvjDvguuwze1atD+qwRY6Ubooi49u2rP5Hu3Ru2JUfxWOXnw9G6NSQAZXv2aO/IqwL2e+6BbdkyeF95Bb6HH47oWIrGSuXfx6rUNFb20aNhe/99eF98Eb5HH9W1X/Ybb4RtxQp4Z8yA76mndP3umtBqzioqKkJSUhIKCwtrvX+HZMk5c+YM3n//fXz44Yf47bffKokGp9OJIUOGYPz48ejLw2R14sEHH8TWrVuxdu3aoO+Jj49HPK9s7IfD4dDkhnEmORlSFYsB7HbEdeoE6HGDSklhX5mfD7sFboha/Q4B6dePtT16QPjqK8TpOUmXL4nZTp6ELcz/r65jpRdpacDIkcAHH7B9ux3CW2/BkZYW0WEVjVX50rcAwNGpE3PcN7sfW3kKBHurVqpd3zWOVVoaG5dx4yqSNQrp6RH/PlYl4FidOgUAsDdrpv+c27Yt+26v15TzvdpzVijHUux4PG/ePKSmpuLtt9/GoEGD8Omnn2LLli3YvXs3fv31V0yfPh1lZWUYPHgwhg0bhr1794bV+VCZPHkyvvjiC6xatQpOEz1RlCQlQbr2WvkFvYv/8dIOVL+qOjxHTqdO+j+FcosnOR4HZ9Qo/UJwRRGYMEHet0p0Sl4ea/V0PHa5gClT5P0nniDnY3+McDzmWNF5XicUW3LWrVuHVatWoVuQ9dd+/fphzJgxWLBgAdxuN9asWYMOgbKXqoQkSZg8eTKWL1+O1atXI82ETxRSmzZs4+67WdFBPW+ovEgnVSKvjt4lHfzxDyGXJP3z9JgZbokdPFi/a6Wm6BQTPTRVw4joKlGUK5IDsiAcOtTcY6UXXOQY4bpBIicoikXO//73P0Xvi4+Px6RJk8LukFIeeOABfPDBB/j888/RsGFD5Jdf9I0aNUK9evU0/34lCDwXyiWX6D8JcEtOfj6btK1SfFAP9C7p4A+fAMvKgNOnAZX9wSzLvHnAoUNs+957gZISfSw5PHmmv9Axe/LMkpKKpZGKPFx6YFVBqBdkyTElYeXJGTBgQEUItlFkZGSgsLAQAwcORHJycsW/ZcuWGdqvSpw8yVojlH2zZsxK4POxauSEjBElHTgJCSxvEkBLVhxRBB55RN7Xc8mIJ8/kFjVBMH/yzNdfl7d79NBvyYiyqQfH55PneyNEDn9gI5FTjbBEzoYNG1BSUlLt9aKiIjz22GMRd0oJkiQF/Hfffffp8v2K4JYcI0TOu+/K1Zz796e1c3+MXK4CKOtxVYyoPO6Py8WsRwAwebK5nY6NzHDLBaG/VXjkSO2/1wqcOiWfw0YuV/EHOKKCkETO8OHDMWvWLAiCgGMBBvPs2bOYO3euap2zOoJRlhyjs/qaHSMtOQA5H1fFiMrjVeE+fX4Ro6bE6Ay3LhdzCu/Yke2//z5lPwbka7lBg8oFZvWCz2Vnz7J/RAUhhZC3adMGX331FSRJQo8ePXDhhReiR48e6NGjB7p3746tW7cimQoQyvAndb3Nl7R2HhxJMt6Sw9Pj790LDBliTB/MRNVknHpHIgIAT9D511/6fWc4dOjAltSqpKbQfclozx55mxyQjZvrOQ0asKK2586x+a08pJwIUeTMmzcPAHMuXrt2LY4cOYLs7Gxs2bIFy5cvh8/nw8svv6xJRy2H1ys7B+ptybGiM6VenD3LJgLAGJHjdgPr17PtyZOZf46Zl0f04JNPWNumDVtmbd9e/5slzwpudpHjdALDh8tjZoQgrGl5MVZFjpFOxwATvs2bMysbiZxKhFXW4ezZs4iLYx+9+eabVe1QtOA4e1ZOBNi4sb5fztfOeeIuKzhT6gVfqkpIYE8/elJ1GVGS6AnY7WZjALDoqj//BAYO1L8fVhE5gHwDGzECmDtX/3Mn0EOUzRbbD1E8L5wGJYIU06yZLHKIChT75Bzi4Z1AhcCpicMxnp+lzpkzbKNhQ30yHFfF5QJeeYVtX3UVWQs4Ri5VGe1PYTa46PMrymmY7xgXOQUF+n93qHAh1quXMeK4akQawH67FSv074sZcLuBqVPZ9tq1xvkncb+cn38m/0s/FIucvn37Yty4cfjtt9+CvqewsBBvv/02unbtik8//VSVDlqVCpFjYE0v8MSNFMUjY2SOHArBrYyZRB/3yTl+XN/cM+HAhZiRhXeHDq0ucsaPB7KyjOuTEVQV6oBxQp27R8ydS87gfihertq5cydefPFFDBs2DA6HA3369EHLli1Rt25dnDx5Ejt27MD27dvRp08fvPLKK7jWv6RBDOI4fZptGClyKOtxdXbvZq0RSfj4E/D48fLNPZaXEc3iRAvIvhSSxB4KjBQQtcFFDhdmRhBIoPp8wKWXWqP2l1qYJchDFOWM4QA5g/uh2JLTpEkTzJ49G0eOHEFGRgYuuugiFBQUVNSouuuuu7Bp0yb88ssvMS9wAJNYcnjW4xMnWJbUWMftZvV2AOCHH4x50nG5mLABgJ49Y+dmEAins7KgMcKJluNwyL5zZvfL4f0zUogFskoC7OY6fnzsLJeYxTprdK4pExOy43HdunUxfPhwDB8+XIv+RA2msORccAGL3ikpYQX9TFjfSzeC+X8Y8aTDc4zEej6LN96QHTYFgdV3M1L0NW3Kstaa3S/HDJacQFZJjs/Halxxn0COKALr1rHttDTgzBkmEqxsaeDjMHYs27fZjBHqZrKKmoywMh4TtWMKS44g0JIVx0z+H5QMkN3wHnpI3pcklsnXSAuAFXLlnD8PFBaybaOX1Fwulg4hUJHZOXOAjAzmo/PRR6zSe+vWwO23s3/9+gGDBrHXqoohq+Fyyf6PbrcxQt3prFwaxUirqMnQTORs2rRJq0NbAocZRA4gL1kdOWJsP4zGLGZlQPb/OHmyuvCKFcxoXrdCGDkXxnY7s9QaTd++lW+uHEkCJk1iYub229kNt+rvzd/3+ONMBFl5iYvP99xKawRjxrC2fn0WSh7LS+F+aCZybr31Vq0ObQnqcZOy0dW/uSUn1kUONytzjDIrA7LwlSQ5IiLWCFQ+wWjzOhc5GzaY94bLBdiFFwb2iTGChx+OvC9vvWXtiCAj6xRy/Es7GG3lMxFhJQPk3HbbbQFflyQJJ2I4bFlYvBjJPKvtCy8AqanGqWpuyfn1V+Cf/4xt86XLBTz6KBMWK1YAf/+7Mf2oU4clIjxzhk2ORlv79MbtZokq/TGDeZ0/CLz7LvDee+aMEjKDP05VavLPCQXutGy1iKCyMnkJ0chruXFjIC6O9efYMSAlxbi+mIiIRM7333+P//73v2hQJXOsJEn46aefIuqYZRFF2CdORMUqtdFZbXNzWfvRR8DHH5tz4taLsjLZctK9u6FdQZMmTOQcPx5bzoGB8orYbEyE9+1rbL++/VbeN2sILrfkmEnkAGxO6d4d6N8/8LKUUoI5LZsZXogZ0D+7vT+CwPJ/HTnCkp6SyAEQ4XLVwIED0aBBA1x99dWV/g0cOBC9evVSq4/WYu9eCGZxcBVFJmw4sV6NnD8F22zG1Zjh8O+PNYtnsPwqRkeamdFHKBA8Go0XeTUTffsCb78d2BEZYK/ffz974Hr66eDvmzfPWnMUv4YbNWKWFCPhS1Y86SkRniXn9OnTaNiwYY1Zjb/77ruwO2VpOnSAZLNVFjpG+RpQIb3K8JIOSUnG+0pxkRNrEVaBMuIa7YsDWCME1+0GnnmGbX/3nXGRPDXhcjHr16+/sv3UVOYECwADBsjzzogRzNn4kUeY6PHH62WfHzFCr15Hhhn8cThc5FD9qgrCsuRceeWVyM/PV7sv0YHTCe+bb6JiqjTS18BMEUVmgD/d8InASGIxjDwrS07G6M+sWcaLbqeTWRc4ZvAR8sdMdb5qw+lkAmXECGbd4dtVx9LpZKHmgZyWR460jhOy0RXI/SGRU42wRE6fPn3Qv39/7Nq1q9Lr2dnZuO6661TpmJWR/vEP2Sdn1y7jnracTmD+fHnfbBO33hhZnLMqsbZc5XYH99fo00f//gSCO0PHxQEHDpjLSmKmPE9qwp2WqwodKy2tkyXH1IQlchYtWoQxY8bgiiuuwNq1a7Fnzx7cdttt6NOnD+Lj49Xuo/UoP+mlhATjrSb33y+HE375pbkmbr0xsjhnVaLNkiOKwKpVcvK3jz5ir4ki267qbMwxk2WRO/OWlTH/CjMRzVZZlwv48MPqr1tFxJlJ5PC5jUROBWF7SU2fPh116tTB4MGD4fV6MXToUGRlZeGSSy5Rs3+WRODe9mY46QG2Lv7XX4Fzk8QSZlquiiZLzuzZwL/+FTh8uKqfiz9G5ioKREIC+1dczK4XI4q4BsMs5QO04rLL2P/JDL6MoWLG5SpyPK4gLEtOXl4eHnroITz33HPo3LkzHA4HRo4cSQKHYyZlD7DU6YAcTh6rmGm5KhosOaLIlngeeyx4fpSaBM769eazLJo567HLJZ+7X39tvrGLBC7i/COu0tOtIeLMNN/TclU1whI5bdu2xc8//4z//e9/2LRpEz799FNMmjQJL730ktr9syZ8ucoMJz0g50uIdZFDlhz1mD2bnVeLFoX3+WnTjM2LEwwzixxJknOydOlibF+0wOVigpnzxBPWcD42o8jJzbWGP5MOhCVyFi9ejOzsbFx//fUAgKFDh2LVqlV47bXXMGnSJFU7aEUEftIbmRjKHxI5DDP65HC/FSvxyiuVb0ahYrOxUgBmhPvlmLESeVER4PGwbbMlA1QDUWTimWMV52MzLVf98ANrCwutXSZDRcISOSNHjqz22iWXXIJ169Zh9erVkfbJ+phJ2QMkcji8EnskGVnVYtUq1ubnW2syEkXmfxMudjtbljDrMoSZLTlceNWvD9SrZ2xftMCqEWRmme9FkZWt4VhFJGqMqhXeUlNT8csvv6h5SGtSblKWyJJjHhYtAvLy2PYttxgrKkRRTuoGWGsy+uWXmv1snnqKRVNNmCAnXOSvr1pl/urIZhY5Zi3poBZWjSDjFuKyMmP7YVWRqDGql7FtbJYbu4EIZlH2HO54fPgwO+ljDVFkIoJjtKiw6mTkdrMkbYGYMAE4eJAVpB0xAsjIYIJm1Sr59YEDzWvB4ZhZ5HBLTrRWmObOx/7ZyIOdb2bB7ZYfHocPN/bhyaoiUWNUFzkEzOd43Lw5O9m9XmDzZqN7oz9mExVWnIx4xt2q2GzAyy8zURMoo60VhI0/3EpiZpETrZYcgFn5cnLka+H99827nFv1mjD64YmLRE60pRkIExI5WsAjIMxi1Xr3XdmCc+ml5pwwtMRsosKKk1EgoQgAS5dG5oRsNriVxIyOx1x4Raslx58//5S3jRYPwTDbwxPARCK3fk2dau6lYZ0gkaMBplquMtvThhE4nXLKfsAc5S1cLqBHD7a9aJH5J6ONG6u/ZrezoovRhBWWq6LZkgNYpyK82R6eOB06sLa42Nh+mISwMx7/8MMP+OGHH3Ds2DH4qqjZd955J+KOWRozOR7X9LRhZsuB2nTuzNqrrwb+3/8zx/+9RQvg998rJ0AzI6Jo3sKaamNmkRMrlhwuHsye/djpBB5/nF0HgDkengAgOZm1R44Y2w+TEJYlZ+bMmRgyZAh++OEHFBQU4OTJk5X+xTSSJIcUmiFvglmfNvQmP5+13bsbPwlxrJIQMNhSlVkKa6oJFxBnzpjPchArlpyq2Y8FwRziIRCXX87aiy4yT+Rgy5as5dGkMU5YlpwFCxbg3XffxT333KN2f6zPmTMQeCihGZar+IQxfrx8ozLrhKElXOS0aGFsP/zhIsfspR0CCeJoFcoffyxvd+zIrh0z3LiA2LHkAGzMd+9miSdvvNE8v0FV+LWbmmqeOZUsOZUIy5Jz/vx5XHbZZWr3JWTmz5+PtLQ01K1bF71798bPP/9sdJcqTnpvXJx5bl4uF1uiAdiNyawThpaYUeRYpX7VH39U3jeLWV5tzJZqoCqxYsnhDB3K2g0bzPMbVMVMVnsOFzn5+cFrysUQYYmcsWPH4oMPPlC7LyGxbNkyTJkyBU8//TSys7Nx5ZVX4tprr8WhQ4cM7RfK/ZHsZWWIa9/ePJFMvE7Q4cPmyPirN2YUOVZYrnK7gfLyLQBYRlWzmOXVxozRMv7wootGJ53Ti507WXv0qHnDyM0UZMLhc1xZmTmjBHUmrOWqkpISLFy4EN9//z26d+8Oh8NR6e9z585VpXM1MXfuXLhcLowdOxYA8Oqrr2LFihXIyMhAenp6tfeXlpaitLS0Yr+oqAgA4PF44OH1YCJFFBH3/PPgbqSCzwfp/vtRNmiQ8U+9LVogThAgnDsHT16eaUzefOxV+w2CEJefDwGAJylJrv9jMEJiIuIA+AoK4FXQJ73GqgJRRNz48RD8RLE0bx7KJk0yzRgGI6yxSk1FnM0GwU/oSHY7ytq0Mfz/KyxahLjTp1mfrrkG3owMSKNHq3Js3c8rJYgi4h5+uGIuhUnm0qpjZfvrL9gBeBs1gs9E4xfXtCmEv/6C59AhQ1OZaHVuhXK8sETO1q1b0bNnTwDAH1VM2YIOkSLnz5/Hpk2b8ESViI8hQ4Zg3bp1AT+Tnp6OmTNnVnt95cqVSEhIUKVfSdu24fIqVhLB68WG99/H8W7dVPmOSBjSpAnqHT+OdR98gFM8zNAkZGZmandwrxc3Hj0KAcCP27ejxCQOec3278cAAEU5OVjzzTeKP6fpWPmRtG0bLq9i2TDT+ayEUMeq9cSJ6PnmmxAASIKALRMm4NDWrcDWrdp0UAF1Cwow5IEHKvYFnw+2iRORabejRMWlK73OKyWY/dzjY9V7+3Y4Aew4ehT7Q7iGtWZg/fpo9Ndf2PjllzjGa/YZiNrnVnEI4fGCJFlv7eLIkSNo1aoVfvnll0q+QS+++CKWLFmC3bt3V/tMIEtOSkoKCgoKkJiYqE7HRBFx7dtXfxLcu9d4Sw4A+8CBsK1bh7IPPoD0z38a3R0ATJFnZmZi8ODB1SyCqnHsGBxOJyRBQNnZs0Bc2JkTVEXIykLc5ZdDSklBmX/ysyDoMlb+iCLi2rWrbMkx0flcE5GMle2++2D/4AN4J0+Gb84cjXqoHGH1asQNGVLt9bLMTEhXXx3x8XU/r5Rg0rm06ljZr7sOtu+/R9k770C6+27D+lUV+403wrZiBbzTpsH34IOGjZlW51ZRURGSkpJQWFhY6/3bHLN9mFS1GkmSFNSSFB8fj/j4+GqvOxwO9QY/LQ1YuBDS/fdD8Hoh2e0Q3noLjrQ0dY4fKWlpwLp1iBNFwCyTWTmq/g5VKXfsFZo2hcNM1ZubNwfAkkeG8n/XdKz8SUtjyf64ddRs57MCwhqrtm0BAHafD3YzXCcXXxwwb0xcp06qXse6nVdKKJ9L/aNCzXTuVYxVuU9OXLNm5ppTz54FANjnzoX91VcNjxJU+9wK5VhhZzw+deoU5syZg7Fjx2LcuHGYO3cuCgsLwz1cSCQlJcFutyOfO5OWc+zYMTQvv3EYhsuFsr17sfa559hTh5kcNFNTWfvLL+aNVtACMzodA7Kz4tmzgJ+V0TRIEnDgANt++unodTiuCl8CMovTphkzduuBy1U5su/WW43rSzDMGF0limyO55gtSlBnwhI5GzduRLt27TBv3jycOHECBQUFmDdvHtq1a4fNOhSArFOnDnr37l1tnS8zM9MUoe1wOtm6sdkmIV4t9/PPzRutoAVmFTmNGsmJGs0YRj5rlpxQLD0dWLHC2P7oBb9hmUXkACyJJQBceWXsiE2AWbHKLWv4/Xdj+xIIM0ZXWaUshk6EJXKmTp2Km266CTk5Ofj000+xfPlyHDhwADfccAOmTJmichcDM23aNCxatAjvvPMOdu7cialTp+LQoUOYMGGCLt9vOURRzpUDxJa637WLtQ0bGtuPqths8uRotjByUWTWG04snS9ms+QAciLAzp3N9/CkNbzG2//+Z67zz+MByqN0TSVyOnSoXiomWpN3KiBsS86//vUvxPk5cMbFxeHxxx/HxkCF/DTg9ttvx6uvvopnn30WPXv2xE8//YRvvvkGbdq00eX7LYfZc4Bohdst15b59FPzWa/MmhAwlp8Gucgx028Sa4kA/fF6WZuRYS4LtH8JIzPUKeQ4nZVrzcXK8mYQwhI5iYmJAZPu5ebmoqGOT8uTJk1CTk4OSktLsWnTJlx11VW6fbfliMUaVrwCO79ZS5L5rBFmTQgYKMVAtJ8vHH9LjlmCT2OppIM/ogh8+aW8byaLIr9mL7iAXRtmgq9o2O3Mry5WljcDEJbIuf322+FyubBs2TLk5uZCFEUsXboUY8eOxR133KF2Hwk14DWsODZb9Kt7K1ivuCXHbM7gVR9WYulpkAvPkhIghHwcmhKrlhwzWxTN6I/D4f6HXi9gpohSAwgrhHz27NkQBAH33nsvyspTjDscDkycOBGz+NIAYT5cLuC//wXWrGFLONGu7rn1qkroramsEXyinDMHmDfP8FDPCrKyWOt0snOmffvYEDgA0KABUKcOcP48Exf16xvdo9gVOWa+hvlyppkiqzh16rBzpaCABQ/E2nnjR1iWnDp16uC1117DyZMnsWXLFmRnZ+PEiROYN29ewFw0hIno2pW1Zlse0QKzW69EEVi/Xt43kyl+5UrW9uoFDBxonjHTA0Ewn19OrC5XmfkaNrMlB6Bq5OWEnScHABISEtCtWzd0795dtdIIhMa0a8daBRl2o4K77pK3f//dHFYSjllN8W438MorbPurr8zj6KknZgojl6TYteQA7JrlOXIefdQ81/D+/aytW9fYfgSDixyTlLExCsXLVdOmTcNzzz2H+vXrY9q0aTW+V48CnUSYcDOv0TdSveCVm+vUAbp0MbYvVeGhnv5Cx2hTPHfW5nBn7aFDzfH0rBdmCiM/fZotnQGxKXIAoG9fYPlywAR1mABAWLwYeO45tvPFF+xBwCzii9OyJWtj3JKjWORkZ2dXVP7Mzs4O+j49CnQSEcAtObt2seSAKSnG9kdr/BMBmu3c5JlsuTneDM69NTlrx6LIMcNyFRdaCQnsXyxy8cWs3bHD2H6AFUy1T5xYPWrTbA8CZMkBEILIWbVqVcX2kiVL4HQ6YasSkixJEnJ5Vl3CnPz0E2vPnWNlHszi6KoVZs12zBkxgv0GqanAzz8bP0ma0bpkBGay5MSqP44/nTuzdtcuJsKrpsPQkQZ5eZUKhwIw54MAWXIAhOmTk5aWhoIAF/+JEyeQZpICakQARBF44AF530yOrlphdpHDnRZLS80xQTqdlWsEmcG6ZATcJ+f3342/PmLZH4fTti0rgHnuXGVnfQM4k5wMyQo5x8iSAyBMkSMFSZB15swZ1DWrExZhjbwxamN2keOfDNAsiefq1GHtuHGxVSfJH35NmKHOG1lygCVLWBkFALjiCkN/j5KkJHgzMuQXzBTx5Q8XOfv3Gy/UDSSkPDnc4VgQBDzzzDOVIqq8Xi82bNiAnj17qtpBQkXMnHNCK/hTDL/gzYa/Jae42Bw5WbZsYe2tt5pv4tYDUQSWLZP3ucXTKJ8LLrhiNambCZ3hpdGjmeOxKDKH6JtuMqQfNfLzz6w9epQJ9Wh3TQhCSJac7OxsZGdnQ5IkbNu2rWI/Ozsbu3btQo8ePfDuu+9q1FUiYnjOCe6AKwjmfAJRE7Nbcho0YGZ4wBy5i4qLgd272XazZsb2xSjMFNrvdgMvvsi2P/ssNsP5zWqBLixkLXeKNhOiCDz1lLwfC64JQQjJksOdj0ePHo3XXnsNiYmJmnSK0BCXiz0R3nUXkJYW/cre7CJHEJg15+hRFsljdLTbiy/KN/h+/WLz6c8sFs9gtdfMFsWjNWb5PfzxeFhoP2DOZIAUJVlBWD45ixcvJoFjZXgh04MH5fwb0Qp/cjEwGqNWuF+O0eHKoihbDYDYffpzOiuPg1HO12a1YOgNt0D7X8NGW6B5BXJBYAU6zUYsFmQOQli1q5599tka//7MM8+E1RlCJ1q1YgUYT59mEyYPz4w2Fi2Sk4fdeqt5rRJmqURe0zJNjD394YEHgCeeYNs7dgAXXaR/H8xowTAKl4tZFrt3Z/sjRxrbH/5AYsYK5IAsDMeNY9e0WZ2jdSAskbN8+fJK+x6PBwcOHEBcXBzatWtHIsfsCALQqRMrwrhzZ3SKHFFkVgiO0c6jNcHN3UZbcjp0qP5arN5U69cH4uOZQ7hR9fj4jWrsWLYfwzcqAEC3buxaOXGCCe8ePQzrisAtOWZcquK4XMz5eMkSYOJEcz7g6UBYNnx/h+Ps7Gz88ccfyMvLwzXXXIOpU6eq3UdCCzp1Yu3XX0fncoSVTP1mseQ4nUCjRvJ+rObIAdiDgBnqV7lcAHcN+OGHmL1RVcCF+N69xvbD7MU5OdwpmjtJxyCqOSokJibi2WefxX/+8x+1DkloydmzrF282Pg8IFpgpTVps1hyTp2SJ8Mvv4zdHDkcM2Q99niAoiK23bWrcf0wC3zZkESOMnggQwxXIlDVG/PUqVMojGHFaBl4bgdONDqYOp3AqFHyvpmtEmax5OzcydpWrYAbbjDnWOmJGepX8e+22YDGjY3rh1nglpw1awydrwR+rfJr16yQyAnPJ+f111+vtC9JEvLy8vDf//4Xw4YNU6VjhIbEioNp27asHTYMePtt8/7fzGLJ4cUPo9FHKxzMYMnh392kiTkdXPWGC5sVK4xNcGc1S44oGl7zyyjCEjnz5s2rtG+z2dC0aVOMGjUKTz75pCodIzQkVqI2eGRV377mFTiAeSw5Gzaw1uhcPWbBDD45VNJBRhRZxCTHyGACq4icVq2Yf9n58+xcat7c6B7pTlgi58CBA2r3g9CTWIna4CKnVStj+1EbfKI8dIhN5Eb8Dm43s3YBzE/rssti2x8HMJclJ5aLc3JMlODOMstVDgcraXPkCFuyikGRE3u2K4LhcgGjR7PtMWOi84bGTdtmFzlr1rA2N9cYJ/BgtYGiyUcrHMzgk0MiR8ZMwQRWCCHncMvsN9/E5DWt2JLDi3MqYe7cuWF1htCZK69kT+1mDKtWA27JMbOFShRZoT+OESZ4Ez0hmwozWHJouUrGRBZogQtfK4icsjLWTp8OzJxp3qSoGqFY5GRnZyt6n8CLPxLmp1cv1mZlMStCNPlinD8PHDvGts1syTGDwIgVH61QMYNPDllyKuNyscjQr78G/vMf427WXHx6vcZ8v1JEEdi8Wd43c1JUjVAscnhxTiKK4I6mZ88CqanRpfDz8lhbp465bxBmEBhOJ0v//tZb8vdHo49WqJAlx5x06cJEDl8y0pnWmZmylfiWW8w9b8ZKJG0NhO2Tc+rUKcyZMwdjx47FuHHjMG/ePMqRYyVEEZg0Sd6Ptlw5/v44ZrYuchM8xygTPLd2DR1KSQA5/j45VW8UekGWnOrw1BD79+v/3aKInvPno2JGMfu82aFD9fkvxqy0YYmcjRs3ol27dpg3bx5OnDiBgoICzJ07F+3atcNmf9MYYV6sVPYgHLZuZa0Vbg4ul1x4cNEiYwTGnj2svfrqmHnCqxW+XFVaKmcI1xuy5FQnLY21BkT5Cvv2QQhmGTEjTifgX0syBq20YYmcqVOn4qabbkJOTg4+/fRTLF++HAcOHMANN9yAKVOmqNxFQhPMFKmgNm43qyINMH8jK5Ss4JYUoywGu3eztmNHY77fjNSvz5Y7AWDbNmP6QJac6nBLzoEDul8vUvv2kKxmGfEPGvrjj5iz0oZtyfnXv/6FuDjZpScuLg6PP/44Nm7cqFrnCA3hyyT+WVT/7/+sr/B5OLT/5GdmczKHP6kb4f8hSbIlh9cGIoB33mEO7ABwxRX6i2VJkp3neYQMAbRuzdriYmDLFn2/2+nEvptvlvetYBlJTJSjwPj5HEOEJXISExNx6NChaq/n5uaiYcOGEXeK0AmXiz0N+ZvlzS4GasOqy3D8SZ0vT+jJsWNyYc66dfX/fjNSNXeQEb4X8+ezAp0AS85oBYukHvz3v/J2nz66j8tJ/iDQtat1/NcMXOIzmrBEzu233w6Xy4Vly5YhNzcXoihi6dKlGDt2LO644w61+0hoSUqKfAFMnWr9iuRWXYbjlhwjRI5/XquOHa39+6uF0WJZFIGHHpL3ze7gqhcmEJ/x/IGgfXtzW3D8iWGRE1ZZh9mzZ0MQBNx7770oKzejOhwOTJw4EbNmzVK1g4TGiCKwaZO8H2keBVEE1q1j22lpwJkzTHjoNRmYKGFYSBi1XCWKwCuvyPsxmEcjIEaH9pshf5IZMcG41CkqYhtWcgbnIicnx9BuGEFYlpw6dergtddew8mTJ7FlyxZkZ2fjxIkTmDdvHuLj49XuYyVycnLgcrmQlpaGevXqoV27dpg+fTrOx+BaoyrUlEchFEQReOwxtl5+++3sX79+wKBB+luH7rpL3s7OtoY52ajlKrV+/2iDi2XuZCoI+oplCv0NjAkstXVOn2YbVnIGj2FLTlgi59y5cyguLkZCQgK6deuGRo0aYeHChVi5cqXa/avGrl274PP58NZbb2H79u2YN28eFixYgKeeekrz745K1Jg0Zs9m4mb27MDRDj4fMzHrZVLOzWVtQgLQrZs+3xkpRi1XdehQ/TW6mTJcLuDhh9n23XfrK5adTuC+++R9Kzi46gEXn/5zls7jUmHJsaLI2bYt5pY8w1quuvnmmzF8+HBMmDABp06dQv/+/eFwOCry5UycOFHtflYwbNgwDBs2rGK/bdu22L17NzIyMjB79uygnystLUVpaWnFflH5ierxeODhzn0qwY+n9nE1oXlzCBkZsE+YAEGSIAHwvvACpObNZafHGhDmzoX9iSdQa7o9nw/eefPgq7KcqcVYCfv3Iw6AlJJSsZxqei64AA4AUkEByoKMhSbnVfPmiGvdGkJ5IIFkt8M7f77i39+sqDVWttatYQfgKy6GV+fxsHXqxL570CB4Fy1iN3IN+mCp+QoA7r0XQqNGiBsxAlKrVii7917dzlWPx1Phk1N2wQWQLDJmwm+/sZv9gQOQ2rSBNyMDEi/QrCFlOTlI2rYNZZ07s6z6KhHKuRqWyNm8eTPmzZsHAPj444/RvHlzZGdn45NPPsEzzzyjqcgJRGFhIZrUUigtPT0dM2fOrPb6ypUrkZCQoEm/MjMzNTmu6jRvjh5//ztSMzMhALA/+SS25Obi0ODBNX6s0Z49uFqJwCnHNm8efm7RAoUBwpTVHKvW33+PXgCOJSRg/TffqHZcLXGcOYPrAAinT+O7zz+Hz+EI+l5VzytJwnUFBXAA2Dx5Mv7q0QMlSUmsYnEUEOlYtcrNRR8Ax/fswTqdx6Tzhg3oAOBA/fr4Y+tWOcGlRlhmvgIQf/o0hgFAXh6+/eILSHFh3crC4ury5aqsnBwcs8B1UregAENmzKjYF3w+2CZORKbdzq51jb6z7Vdfof3nn+NySYL0zDPYMmlSrfcUpRQXFyt+ryBJoWdTSkhIwK5du9C6dWvcdttt6NKlC6ZPn47c3Fx07NgxpA5Eyp9//olLLrmkosREMAJZclJSUlBQUIDExERV++TxeJCZmYnBgwfDUcPNyjSIIuLat4fg59An2e0o27s3qBlYWLy4wvoTCpLNVukpQouxss2cCfsLL8A7dix88+erckzNkSTE1a8PoawMnv37A467JufViRNwtGjBjn/qFFviiwLUGishMxNx118PqVs3lPk76OuA/f77YVu8GN6ZM+F78knNvsdy8xXArpdGjSCUlMCzcyfQrp0uX+vxeIA2bZBQUICydesg9emjy/dGgrB6NeKGDKn2ellmJqSrr1b/++bMgf3JJ6s9/NZ2TwmFoqIiJCUlobCwsNb7d1jyt3379vjss89w6623YsWKFZg6dSoA4NixY2ELhhkzZgS0tPiTlZWFPn4n1ZEjRzBs2DCMGDGiRoEDAPHx8QGdoh0Oh2YXtpbHVpWcnGoRC4LXC8fBg/Jarj+iCEycGDzb6N13A5dfzrIOVz2uz4e4iROB666rdLKrOlbla872tDTYrTD+nKQkID8fjlOnAo97OaqOFc93lZwMR6NG6hzTREQ8VuUCUDh+XP9r+fhxAIC9RQtdzmPLzFectDRg5044cnOBTp10+1qh3NUhrkULwArjdfHFASMF4zp1Ur//r7wCBBHkNd5TQiSU8zQskfPMM8/gzjvvxNSpU3HNNddgwIABANjST69evcI5JB588EGMHDmyxvek+q3pHTlyBH/7298wYMAALPQvcEiETqBwWZuNpbUPxGuvVQ/jBFg0yEsvsSgrgF1A48dXf6/Px47hH7qsJgcPsrZNG22OrxXlIkfXMPI//2StTk/CloMnyiwoYKJez2KvVNKhZtq2BXbu1Ddi6OxZxPFIXqv8LtxZe9w4+RzWwllbFIF//Sv43w0KaAhL5Pzzn//EFVdcgby8PPTo0aPi9WuuuQa33nprWB1JSkpCksKT5vDhw/jb3/6G3r17Y/HixbBVjQ4iQoNfBPffz8KHASZELr2Uve4fVZKVxaKoqmKzAevXA337yq/xwpP9+1e3+sybxyJXmjdX///DqxNbLXuvERFWfKx4PSCiMnxOOn+e5XzSM6M7FeesGW4R0LMaebnwlOLjITRooN/3RorLxUTIjBnAsGHaRAquXRvUui/Z7RAMig4MWx20aNECvXr1qiQw+vXrh04amw2PHDmCgQMHIiUlBbNnz8Zff/2F/Px85Ofna/q9UY/LBfz6a+Un1aqh37Nns9w3gZg2rbLA4fTtCzzySPXXtcrF8vbb8hLMbbdZK3svv5mtW6dfmKeVqrUbQUKCLJbLl490gyw5NcOF+YYN+l0v/BxIStLXqqcGfO7mKTbUxO0GAlQ7kAAcGDqU+eIYlK/MciaQlStXYt++ffjxxx/hdDqRnJxc8Y+IkDNnqitxn4/539x9t7wMVRWbTc4nEoiHH66ei6em5bBwEUVgwgR532qp8PmT+5tv6pNA0e0Gli1j2/PmWUsQ6oUgyCJDz2XEsjLg5Em2TSInMHypdfVq3RKOCvwc4MuYVqJjR9YGyhodCVVLbXBsNnhnzcLWiRMNze9kOZFz3333QZKkgP+ICAmUZRUAvvoKeP/9wJ+x2diSVk0ncaCK5+XLYcLixZH12R+j6w1FgiiyyZqjtUCrOjFJkrUEoZ74++XoBbcYCIJcQZqQEUUgI0Pe1+uBpvxBRLKi8GzTBqhThxViDlBgO2yC+WguXQpp2jT1vidMLCdyCA1xOgMvLQVDEJgfjhIzZJDlMPukSair1s3Dyqnw9S6vYGVBqDdGWHL4dzVpUvnhgGAYdP4K3HpkxVQL/nPhRx+pIwhFEZgzJ/B3lQckGQ2JHKIygZaWgvHyy4H9cIIRYDlM8HpRPy8vhA7WgNMJXHutvG+lVPh6CzQT1ACyDFzk6OmTQ07HNWPE+et2w/bCCwAA4euvrbm8W6cOa//1L3WW+AI9nAHA1KmmmXdJ5BCVCVQbpio2GxM4jz4a2rEDTEySzYazavpT8YiHiRNZ/h8rFOcE2Lj7m3a1FmhOJzB5sn7fZ2WMWK4ip+Oa4fMUx2bT9vwtX97lyU8FKy7viiLw++/yvhpLfIHcG2rz0dQZEjlEdVwulmvm0Ucrn8SCwF47eDC4E3JNVK3sDACShGbZ2ZH3mcPN1cOGWe+GPWoUaxMT9RFo/Kn38sutJQj1xojlKrLk1I7Lxa5zAJg5U9vzNxqWd9VeEne7gUGDKr9mt9fuo6kzJHKIwDidLFnfoUNs/fajj9j2K69EdgIPHVpJ5AiShB4ZGeo8EUmStZPb8RvamTOAHtGCfKwuvdRUk5Lp4CJnxw79ntzJkqMMnrKkvGimZkTD8q6a/wceuOAvmmw25ndpsoclEjlEzTidwIgR7J8aN8IAT0Q2n0926IuE48flyc6Kye34sojPJ4cPawlPomZFQagn27ax9qefdAtVrsjiG6AUDeEHz4Kfk6Pt95RbofktXdJ6eUwLnE6WnoITyRJ1IMuWzwecPRtZHzWARA6hLwGeJnw2GyQ1brTc7NqqFVCvXuTH0xuHA7jgAratR9ZjLiytKAj1QhQrixo9QpXdboCnVnjzTWs6uOqFXiIHYBaK8geRsi+/NJ3FQhETJsj5chYtCv//0KFD9ddMatkikUPoS4CcOUfUCjX87TfWpqSoczwj4EtWWvt/SBJZcpSgty8G5S8KDS5yeL06LfH5gFOn2HbXrtp/n1bwItdHjoR/jE8/rbxvYssWiRxCf1wu9uTVpQsAwPnLL4hr3z6yJ1a3G5gyhW1v2GDdp1/ug6G1JSc/Hzh3jvlHUe234OjtixENDq56wovw/vWX9kslp05B4LX9rJjxmNO9O2szM8MTz6Ioz7X+DB0aUbe0gmY3wjh27KjYFCJZBqjqBGflp1+9inS+/jprJYndyK0qCrXG6QTeeEPe1zrUPhocXPXkggvkoqkbNmj7XeXWVU9CgpxvxorwvGThlsPYsydw+R+TCnESOYQxqBnOGE1Pv3osV4ki8NJL8r7VanzpzcSJsgPwTz9p64vhdDIRxaH8RTXjdgOnT7Ptv/9dW7Fefk2e17MSvdqIovyAA4R37Qd6r4mFOIkcwhgCPbHabOFdKNH09KvHcpXeJSSiAS4+HQ7tv2vkSHl7xw5rOrjqgd7+S+XXZGmjRtocXw8ifSB0u4H77qv8msmFOIkcwhh4SGaVxIBYsSK8Yz33nLxv8ouuRvRYrrJQZIRp0LO0A//t69UDLrpI+++zKnpbcLklJzFRm+PrQSQPhBbKjeMPiRzCOKokBozoSaxzZ9a2b2/t7L16LFc5nZXDxq0sCvVCz6zHlAhQGXpbcKNhuapq1nlBUH7tWyg3jj8kcgjj2LuXORz7E+6T2PbtrLV69l69HI/5xPTWW9YWhXqhZ/0qKumgjEB19rQU69FgyQHYtf7FF2w7MREYPVrZ5zZurP6aBSzAJHII4+jQgWUO9SfciyYri7VWFjiA/PSem6udb8GZM8DRo2z7ttusP2Z6QJYcc+JyAV9/zbaTkrQV69Hgk8MZOhRISGAZ4pcsqX2uEUXgiSeqvz5rlunnDxI5hHE4nfBmZFQWOmPGhH4ctxv4/HO2/dJL1g6H/vFH1hYUaFdCgJcMaNJEzrBM1IyePjlc5JAlRxn9+rG2oIDlftKKaFiu4jgcctLUMWNqn2sCLVUBcmJBE0MihzAUafRorFy4ED7+NPD226Hd3KMpQ6woAk89Je9rFdpN5RxCR09LDl+uIkuOMho3lnPlaJn5OFqWqwA2p+zZI+/XNtfExVV/zQJLVQCJHMIkCIcPyzuh3NyjKUeOXv8XKucQOnr65NByVWgIgj7lHfLzWVs1/YIVCSWNhNsNXH115dcsFKxAIocwnAZ5eRDCzdvSoUPlCC3AMk8Y1dArWuT331lLN1HlGGHJoeUq5fDyDloV6nS7KwRUv5degsALqFoVpXnKLBo27g+JHMJwziQnh++A7HRWfsqw0BNGNXi0CEeLonduN/Dee2x7/nxr+y/piRE+OSRClaNlNfIqS+KCJME+aZI1l8Q5VecaIHCeMouGjftDIocwnJKkJHgzMipbZNLTld/ci4pYO3269cOhXS5g2DC2PXOmuv+XaPJf0hv/5SqtlyvIkhM6XORs2KD++RzgRi9YdUncHyV5ynjUqj8Ws5STyCFMgTR6NDB1qvzCE08oszKcOwds3cq2x4yxpgWnKmlprD1/Xt3jRpP/kt5wq4rHI9dK0gqy5IQOd6ZftUr9qMQASzuSxW70AQnml/Prr2w7K8uyYeP+kMghzIEoAq++Ku8rdT5esQIoK2M3BB4SaXVatGAtd3RUi2iq8aU3CQlA3bpse9s27b6nrAw4cYJtkyVHGaJYuaip2lGJTicwebJ8eJsN3vnzLXWjD0ig+QBgtdPuuw/o3z+w1dICYeP+kMghTIGwb1/oVga3Gxg+nG0XFADvvKNdB/VEK5HjdFa2llnZf0lv3G6gpIRtX3WVdr5MXOAAQHGxNt8RbehhoezZEwDg69MHmQsXMsuz1QmUMRpgY7lkSWCBY8GHIhI5hCmQ2rcP7O1fv37gDwTy+o8W/5LmzVnLsxKrSadOrO3Xz/r+S3pR1ZdJq/xFALBggbzdti05hitBDwslvxY7dUJJNC0julzAhx8qe68WgRA6QCKHMAf8qcJul1/z+VgtqkATfTT7l2hlyQFk34W+fS03WRmGXueaKDJnc46WYiqa0CMqsVzkSM2aqXdMs3DZZdXTcFTFZgPWr7fkQxGJHMI8uFzM6c3/ggs20Uezf4m/yFE7kocSAYaOXudaNAt3rXG5gL/9jW2np6t/M+aWHG5ljSacTuCRR2p+z7Rp7MHIgpDIIczFmTPKMnE6ncCVV8r70eRfwifS8+dZAT01oZIOocMtBVx8C4I251o0JbY0gosvZq3a1wwQ3ZYcAHj44cBOyAB7/eGH9e2PipDIIcyF0kycbjewZo28r8XTm1HUrQvwSsdqL1lxkUOWnNBwuYApU9j2nXdqc645ncANN8j70STc9UDLhIDRbMkBArsLAGx/4UJLn4MkcghzUfWpGaieibOqIygAPPlkdPkuaOGXc/IkcOoU2yZLTujwMeNRVlrAbyb33kuO4aGiZf2qaLfkAOxcy8lhuYZ++421UXAOksghzEegTJzjx8vZN9eti37fBf7E+OOP6ok3bsVp3LhyqDKhDD1KOxw7xtp+/Sz99GwIWtWvKiuTEzRGqyWH43QCAwcy/5uBA6PiHLS0yCktLUXPnj0hCAK2bNlidHcItQhWL+XSS1mSqpEjq38m2nwXeI6U555TL4ProkWsPXlS/aywsYAeRTq5yIlmi4FWcEvOkSPqZgv/6y/2oGWzURZqC2JpkfP444+jZcuWRneDUJtgmTiDJamKNt8FUQQ2bZL31QglFsXKYbYUnhw6/vWrtIJETvg0bQrUq8fmh9xc9Y7L/XGSkqr7rBCmx7Ii59tvv8XKlSsxe/Zso7tCqE2wTJzB+PBDy68bVyJITRmBLzepeMyoWuLTGn9LjlZFOknkhI8gAPyhN1BhyXCJdqfjKCfO6A6Ew9GjRzFu3Dh89tlnSEhIUPSZ0tJSlJaWVuwXlVeu9ng88Hg8qvaPH0/t40YjQcfq3nuBiy9G3BVXQKjhhiLZ7Sjr04cVTowWUlMRJwiV/t+S3Q5PmzbAjh3hnVepqYgD4B+gLNntKGvTJrrGrhxNrsHERDgAoKwMnhMngMRE9Y4NAOfPw3HyJADA07ixbr9LtMxXwuLFsP/5JwQA0p13wltUpEr5BWHbNsQB8CUmRs1Y6YVW4xXK8SwnciRJwn333YcJEyagT58+yFHoZJaeno6Z/tlEy1m5cqVioRQqmZmZmhw3Ggk2Vp1vugkdPv886Of23XgjdmzdKlcijxL8/98+mw2/T5iAQzt2AAj/vLouIQGOcl+fimNG4dj5o/Y1eH18POJKS7H6449RzCPgVKLuiRMYCvbbfPPrr8otmSph5fmqbkEBhkyYUCHiBUmCbeJEZNrtEZVhaJ2ZiZ7z57Nj/vILdj/+ODB4sKXHygjUHq/iEOq6CZKkld01NGbMmBFQhPiTlZWFdevWYdmyZfjpp59gt9uRk5ODtLQ0ZGdno2d5EbVABLLkpKSkoKCgAIkqP5F5PB5kZmZi8ODBcDgcqh472qh1rEQRce3bQ6jqiAxAstlQtm9f9Pji+LNrFxzdu0OqVw9l27cDTmdk59WZM3A0aQIAKPv0U0g9e0bnuJWj1TUY164dhNxclP3yCyS1M8Bu2QJHv36QWrRA2aFD6h67BqJhvhJWr0bckCHVXi/LzIR09dXhHTTA3CPZ7Vj51lu48o47LDtWeqLVuVVUVISkpCQUFhbWev82jSXnwQcfxMhAUTN+pKam4vnnn8f69esRHx9f6W99+vTBXXfdhSVLlgT8bHx8fLXPAIDD4dDsZNXy2NFG0LFKS2P+Offfz3xIOHY7hLfegiMtTb9O6kl5pIhw7hwcSUmA39iEdV7xm+aFFyLu1ltV6qT5Uf0aTEoCcnMRd+pUpd9EFcqXqoRmzQyZNyw9X118MbN8+T8M2e2I69Qp/N8pJ6dalKfg9aJ+Xp61x8oA1B6vUI5lGpGTlJSEJAVmxddffx3PP/98xf6RI0cwdOhQLFu2DP3799eyi4RRuFwsd86+fawq+dmzLFw8ii0RaNiQ/Tt9moXEduwY2fH27mVthw6R9y2W0TJXDjkdhw8PVhg/XhYmkUZc8ijPKpacs8nJEXaW0BPTiByltG7dutJ+gwYNAADt2rWDM5pverGO0xndoiYQrVoBu3YBhw9HLnJ4FFU05RIyAi3DyEnkRIbLBfTqBfTuzSKt7r03suNx4TR2LNu32eCdPz8iHx9CfywbQk4QUU+rVqw9fDjyY5HIUQctEwKSyImcXr1Y7TdJUicHlMsF8MCU1atVidYi9MXyIic1NRWSJNXodEwQlkRNkfPHH6wtdz4mwoSLnK1b1U+keOAAawP4DhIKEQSAW/vVKO9w+rScfbxXr8iPR+iO5UUOQUQtaokctxvYsIFtT5lC5RwiYc8e1n71lbqlMdxu4KOP2PbLL9NvFAlqViM/coS1iYlAuWsEYS1I5BCEWeHZWyMROVUrtlM5h/ARRZZdm6PWWFb9jSSJfqNIUFPk8GuPygdZFhI5BGFWuCVn587wb3iBip1SOYfw0Ko0Bv1G6sJFzvr1kQtFbskhkWNZSOQQhFnhRTp37QLatIGweHHox+jQgfkp+BNtFdv1IlDhWDXGUqvjxirct2nlysiXFLnI4Q8chOUgkUMQZkQUgfR0ed/ng33SJNQNNarH6QSGD5f3o61iu544ncCcOfK+WmPpdAKvv67+cWMRUawsaiJdUiRLjuUhkUMQZiTAEgbPthoyF1zA2lGjmJ9CNFVs15vJk2XL2Pr16o3ltdeytk4dZomg3yg81F76I5FjeUjkEIQZCbCEEXa2VT7BDx5M1oFIsdvlPDZxKuZS5TlykpOBlBT1jhtrqL30RyLH8pDIIQgzwrOtciLJtkolHdSFixwuTNSAEgGqQ4DrJqKlP17zTU1BS+gKiRyCMCsuF/DPf7LtRx8NL9vq2bPy0yg5sqpD8+asPXpUvWOSyFEPl4tZLQFg5szwl/4WLQJyc9n2P/5BuYssCokcgjAzXbqwtrxCdcj8+SdrGzaUM7cSkUGWHPPTrRtrT5wI7/OiyByWOZRfyrKQyCEIM8Nzfhw8GN7nFyxg7enT6mbojWXIkmN+2rVjLRf5oRLEgVkI93iEYZDIIQgz06YNa8PJ3iqKssgB6GlULciSY3740my4UVVB8ktJXDwRloFEDkGYGX9LTtUny9rQKkNvrEOWHPPDxcj+/aFfNwBzVL79dnmfchdZFhI5BGFmnE4WIVJaGrrlIFA0FWXSjRyy5Jif1q3ZdVNSImcODxUeNn7bbZRfysKQyCEIM+NwVDw92pYvDy3jcatWQHy8vE9Po+pAlhzz8957sgXn0kvD80XjfnCXX07XjIWh4H+CMDvlQsX+8MMYIgjwer2Vq1YHQxSZBchuB777DujUiSZrNfC35EhSdd+NUPH5gL/+qnxsInyqVnXnvmhDh4Z2/nORw/3iCEtClhyCMDOiKCfzAyBIEuyTJilzHt6xg7UXXQT8/e8kcNSCCxGPRx7jSDhxQrY6lJREfrxYR63SDjwRIIkcS0MihyDMjJ/A4QhKJ2x+A+7cWeVOxTjvvy9vd+8eeVh+Roa83aEDhflHihqlHc6dk5cQSeRYGhI5BGFmAoSySkon7Kws1pIFRz2CLYWEG5YvisCMGeodj5BLO/gLnVB90bgVp0EDucAtYUlI5BCEmXE6gaefrtj1ldewqnXCdruBDz9k26+/TtYBtVC7yrXaxyMYLhfwyy9sOy4OGDUqtM/ziKyWLSP3uSIMhUQOQZidxx6r2Fw1b17tNayqWhskiawDaqF2lesgSecozF8F+vUD6tUDysqAAweUf87tBu6+m23v2UMPCBaHRA5BmJ3ERCA5GQDgKC2t/f1kHdAOvhTChYkgRBaW73QCQ4bI+xTmrx42G3O6B4Bdu5R9hj8g+CfRpAcES0MihyCsQKdOAIAGhw/X/l6yDmiLywX8619se/jwyJPEJSWxdvx4SjqnNuXXDXbvVvZ+ekCIOkjkEIQV6NgRgEKR43QCN9wg75N1QH34zbOoKPJj5eWx9sor6TdSG/47/fijMmuM2suRhOGQyCEIK1Aucprs2KFssuaWnAceIOuAFpQvH1YIlEjgx+DHJNTjyBHWfvstCwWvzb/G6WQPBBybjR4QLA6JHIKwAuVVyJN27EBc+/a1T9bZ2ay97TaaoLWARI75EcXK14nS8PybbpK3d++mBwSLQyKHIMyOKAL/938Vu0Jtk/Xx40BuLttu2lSHDsYgXJAcPw6cPx/+cc6dA06dqnxMQh3C9a/Zs4e1qam0TBUFkMghCLMT6mT94ovydteuFAKrBU2asOKpAJCfH/5x+Gfj4ynpnNqE61/DRQ6PzCIsDYkcgjA7oUzWogjMmyfvUwZdbbDZ5GrkkSxZcZ8RSjqnPuGG+5PIiSpI5BCE2SmfrKXyyVqqabLeu7dyjg+AQmC1Qg2/HPLH0RaXC3jkEbZ9663K/Gt+/521tNQbFZDIIQgr4HLBO2cOAEDq2TP4ZN2uXfXXKARWG0jkWIMrr2StEqHvdgPffce2Z8ygpd4ogEQOQVgEafBgAICwcydLVR+IqnlbKEeOdnBh8uuv4S8H8iR1DRuq0yeiOt27s3b7dmD//uDvo3IoUYllRc7XX3+N/v37o169ekhKSsLw4cON7hJBaEuHDvAkJEAoKQHefTfw5PvFF6zt0wdYtYpy5GgJt8L897/KcrBUxe0G5s9n20uWkNVAK77/nrVeL/NvCzbOlO04KrGkyPnkk09wzz33YPTo0fj999/xyy+/4M477zS6WwShLTYbSho3ZtvjxlW/sbrdcsXyTZuAP/8kC45WiCLw5ZfyfqgO3lVrJJHVQBtEkY0rp6bficqhRCWWEzllZWV4+OGH8corr2DChAm46KKL0LFjR/zzn/80umsEoS2iWLmsg/+ETaZ2fYnUwZusBvoQyjg7ncDQofI+LfVGBXFGdyBUNm/ejMOHD8Nms6FXr17Iz89Hz549MXv2bHTp0iXo50pLS1HqV8G5qNx3wePxwOPxqNpHfjy1jxuN0Fgpx7trFxzVXvSibNcuQJIQF2AyL9u1CxIPdY4hND+vUlMRZ7OxxIzlSHY7ytq0AZR8Z6SfV5GovgZDHGd7WRlsALxTp8I3eTITOH7vi+qx0gCtxiuU4wmSVPVxxNwsXboUd9xxB1q3bo25c+ciNTUVc+bMwcqVK7Fnzx40adIk4OdmzJiBmTNnVnv9gw8+QEJCgtbdJoiIqVtQgCHjxkHwu2R9goDMt98GAAwZOxb+xnafzYbMhQtRwqtcE6qS9uWX6F6+XOiz2fD7xIk4VO4crujz33yD7gsXhv15QhmtMzPRMyMDgs8HCcDW8eORc9111d8oSRh2772IP30aq2fPRiEtU5mW4uJi3HnnnSgsLERiYmKN7zWNyAkmQvzJysrCnj17cNddd+Gtt97C+HLzfGlpKZxOJ55//nnc77/+6kcgS05KSgoKCgpqHaRQ8Xg8yMzMxODBg+FwVHv2JvygsVKOx+PB7scfR88336wQM5IgwLtgAeDzwT5xovy6zQZvRgak0aON6q6h6HJeSRLiGjWCUFICz48/AldcEdrnd+yAo2dPSPXro2zbNsOWRWLiGszNRVz//hAKClC2ciWkgQOrv2fDBjiuvJJZek6dYlmoqxATY6UiWo1XUVERkpKSFIkc0yxXPfjggxg5cmSN70lNTcXp06cBAJ07d654PT4+Hm3btsWhQ4eCfjY+Ph7xAU5ah8Oh2cmq5bGjDRorZRzr1Ys5R5Y/mwiShLiJE6v5hwgA4q67Ti49EKNofl61bg3s2QOHIIQ+1uUlHYS0NDjS0jToXGhE9TXYti0waBDw0UeIW7oUuPjiyqLS7WbO/AAErxeOZctqjEqM6rHSALXHK5RjmUbkJCUlIUmBWb13796Ij4/H7t27cUX5k5PH40FOTg7atGmjdTcJwlAa5OVVWq4CUN2xkr+2bx85TWqN08nKAPCCqKHAncJTUtTtExEYXhrlnXdYCoaFC5mQqRrpBjCn/aFD6fqJAiwXXZWYmIgJEyZg+vTpWLlyJXbv3o2JEycCAEaMGGFw7whCW84kJ0OqWscqEBT6qg9coIQjcvhnSORojygCH30k7/tHJlKkW1RjGktOKLzyyiuIi4vDPffcg3PnzqF///748ccf0ZjnECGIKKUkKQnejAy2RBXIggNQ6KuecIESTqg+iRz9qEnI8AK4/n+nh4SowXKWHICtx82ePRtHjx5FUVERMjMzawwfJ4hoQho9Gvjww+Bv+PBDynKsF2TJsQZcyPhjszEh43QC/nnW6CEhqrCkyCGImOeyy6pnZwXYBD1ggP79iVX4jZBEjrlxOpkPjv81I0nAihXM6dh/KSs9nR4SoggSOQRhRZxO4O23K0/aNhs9gepNuMtVkkQiR2+GDq0ucsaPr5wpHACefJIyhUcRlvTJIQgC7Glz6FBWBRtgFhwSOPrCBcrx48zvo0MHZZ87dQo4e1azbhEBCOSXE8ivjfvq0LUUFZAlhyCsjNMJjBjB/tGkrD8ffyxvd+qkvJL4a6+F9zkifAL55QSCnI6jChI5BEEQ4RBKheuqn3v22dA/R0QG98upSeiQ03HUQSKHIAgiHMLNrxJpBXMifFwuikyMMUjkEARBhEOg5Q8lSx0dOlSPjKMlEv2gyMSYgkQOQRBEOFQNSxYEZUsdTifQtau8T0sk+kKRiTEFRVcRBEGEi8vFbpBjxgBduihf6vB4WDt7NnD77XRz1RuKTIwZSOQQBEFEQr9+rM3NZb42gZZC/PH5gAMH2Pbw4XRzNQoemUhENbRcRRAEEQlt27K2sBA4caL29x85ApSWAnFxlAiQIDSGRA5BEEQk1KsHtGzJtv/8s/b38/e0acOEDkEQmkEihyAIIlLatWPt/v21v5eLHP4ZgiA0g0QOQRBEpPAlq8zM2pP6ZWeztnlzbftEEASJHIIgiIg5dYq177zDlqGClWlwu4E33mDb/+//UTkHgtAYEjkEQRCRIIrAF1/I+8HKNIhi5YrXkkTlHAhCY0jkEARBRILSMg3hloEgCCJsSOQQBEFEgtLyDuGWgSAIImxI5BAEQUQCL+/ACVYiwOmsXLWcyjkQhOaQyCEIgogUlwu4/nq2/eSTwcs7NGvG2qFDgZwcqnhNEBpDIocgCEIN+vdn7aFDwd+zbRtrhwwhCw5B6ACJHIIgCDXo1o2169YFj5javJm1LVro0yeCiHFI5BAEQajBzp2s/fPPwLlyMjLYEhUA3HMP5cghCB0gkUMQBBEpogj8+9/yftVcOaIIPPhg8L8TBKEJJHIIgiAipbYcOJQjhyAMgUQOQRBEpNSWA6dDh+qfoRw5BKE5JHIIgiAihefK8Rc6CxbIEVTJyUC9evLfKEcOQegCiRyCIAg1cLmAXbsAh4PtOxyyz8133wHnzgH16wPff085cghCJ0jkEARBqEWHDkC7dmz7vvtYlNV99wE33sheO3uWCRyy4BCELpDIIQiCUAtRBHbvlvd9PmDJksoFPCmqiiB0g0QOQRCEWgSqSF4ViqoiCN0gkUMQBKEWgaKsqkJRVQShG5YUOXv27MHNN9+MpKQkJCYm4vLLL8eqVauM7hZBELGO0wnMmlXze2bNIp8cgtAJS4qc66+/HmVlZfjxxx+xadMm9OzZEzfccAPy8/ON7hpBELFOnz6R/Z0gCNWwnMgpKCjAvn378MQTT6B79+7o0KEDZs2aheLiYmzfvt3o7hEEEevUtGRFS1UEoStxRncgVC688EJcfPHFeO+993DJJZcgPj4eb731Fpo3b47evXsH/VxpaSlKS0sr9ouKigAAHo8HHo9H1T7y46l93GiExko5NFbKMXSsmjeHkJEB+6RJELxeSAAEAJLdDu/8+ZCaNwdM9BvSeaUcGqvQ0Gq8QjmeIEm1hQKYj8OHD+Pmm2/G5s2bYbPZ0Lx5c3z99dfo2bNn0M/MmDEDM2fOrPb6Bx98gISEBA17SxBELFK3oAD18/JQFh+PuNJSnE1ORklSktHdIgjLU1xcjDvvvBOFhYVITEys8b2mETnBRIg/WVlZ6N27N2655RZ4PB48/fTTqFevHhYtWoQvvvgCWVlZSE5ODvjZQJaclJQUFBQU1DpIoeLxeJCZmYnBgwfDwbOfEgGhsVIOjZVyaKyUQ2OlHBqr0NBqvIqKipCUlKRI5JhmuerBBx/EyJEja3xPamoqfvzxR3z11Vc4efJkxX9u/vz5yMzMxJIlS/DEE08E/Gx8fDzi4+Orve5wODQ7WbU8drRBY6UcGivl0Fgph8ZKOTRWoaH2eIVyLNOInKSkJCQpMOUWFxcDAGxVHPtsNht8Pp8mfSMIgiAIwnpYLrpqwIABaNy4MUaNGoXff/8de/bswWOPPYYDBw7g+uuvN7p7BEEQBEGYBMuJnKSkJHz33Xc4c+YMBg0ahD59+mDt2rX4/PPP0aNHD6O7RxAEQRCESTDNclUo9OnTBytWrDC6GwRBEARBmBjLWXIIgiAIgiCUQCKHIAiCIIiohEQOQRAEQRBRCYkcgiAIgiCiEhI5BEEQBEFEJSRyCIIgCIKISiwZQq4GvGQXr0auJh6PB8XFxSgqKqLU37VAY6UcGivl0Fgph8ZKOTRWoaHVePH7tpLSmzErck6fPg0ASElJMbgnBEEQBEGEyunTp9GoUaMa32OaKuR64/P5cOTIETRs2BCCIKh6bF7hPDc3V/UK59EGjZVyaKyUQ2OlHBor5dBYhYZW4yVJEk6fPo2WLVtWq2NZlZi15NhsNjidTk2/IzExkS4EhdBYKYfGSjk0VsqhsVIOjVVoaDFetVlwOOR4TBAEQRBEVEIihyAIgiCIqIREjgbEx8dj+vTpiI+PN7orpofGSjk0VsqhsVIOjZVyaKxCwwzjFbOOxwRBEARBRDdkySEIgiAIIiohkUMQBEEQRFRCIocgCIIgiKiERA5BEARBEFEJiRyVmT9/PtLS0lC3bl307t0bP//8s9FdMpwZM2ZAEIRK/1q0aFHxd0mSMGPGDLRs2RL16tXDwIEDsX37dgN7rB8//fQTbrzxRrRs2RKCIOCzzz6r9HclY1NaWorJkycjKSkJ9evXx0033QRRFHX8X+hDbWN13333VTvPLr300krviZWxSk9PR9++fdGwYUM0a9YMt9xyC3bv3l3pPXRuMZSMFZ1bjIyMDHTv3r0iud+AAQPw7bffVvzdjOcUiRwVWbZsGaZMmYKnn34a2dnZuPLKK3Httdfi0KFDRnfNcLp06YK8vLyKf9u2bav428svv4y5c+fijTfeQFZWFlq0aIHBgwdX1BeLZs6ePYsePXrgjTfeCPh3JWMzZcoULF++HEuXLsXatWtx5swZ3HDDDfB6vXr9N3ShtrECgGHDhlU6z7755ptKf4+VsVqzZg0eeOABrF+/HpmZmSgrK8OQIUNw9uzZivfQucVQMlYAnVsA4HQ6MWvWLGzcuBEbN27EoEGDcPPNN1cIGVOeUxKhGv369ZMmTJhQ6bVOnTpJTzzxhEE9MgfTp0+XevToEfBvPp9PatGihTRr1qyK10pKSqRGjRpJCxYs0KmH5gCAtHz58op9JWNz6tQpyeFwSEuXLq14z+HDhyWbzSZ99913uvVdb6qOlSRJ0qhRo6Sbb7456GdidawkSZKOHTsmAZDWrFkjSRKdWzVRdawkic6tmmjcuLG0aNEi055TZMlRifPnz2PTpk0YMmRIpdeHDBmCdevWGdQr87B37160bNkSaWlpGDlyJPbv3w8AOHDgAPLz8yuNW3x8PK6++uqYHzclY7Np0yZ4PJ5K72nZsiW6du0ak+O3evVqNGvWDBdddBHGjRuHY8eOVfwtlseqsLAQANCkSRMAdG7VRNWx4tC5VRmv14ulS5fi7NmzGDBggGnPKRI5KlFQUACv14vmzZtXer158+bIz883qFfmoH///njvvfewYsUKvP3228jPz8dll12G48ePV4wNjVt1lIxNfn4+6tSpg8aNGwd9T6xw7bXX4v3338ePP/6IOXPmICsrC4MGDUJpaSmA2B0rSZIwbdo0XHHFFejatSsAOreCEWisADq3/Nm2bRsaNGiA+Ph4TJgwAcuXL0fnzp1Ne07FbBVyrRAEodK+JEnVXos1rr322ortbt26YcCAAWjXrh2WLFlS4bxH4xaccMYmFsfv9ttvr9ju2rUr+vTpgzZt2uDrr7/G8OHDg34u2sfqwQcfxNatW7F27dpqf6NzqzLBxorOLZmOHTtiy5YtOHXqFD755BOMGjUKa9asqfi72c4psuSoRFJSEux2ezU1euzYsWrKNtapX78+unXrhr1791ZEWdG4VUfJ2LRo0QLnz5/HyZMng74nVklOTkabNm2wd+9eALE5VpMnT8YXX3yBVatWwel0VrxO51Z1go1VIGL53KpTpw7at2+PPn36ID09HT169MBrr71m2nOKRI5K1KlTB71790ZmZmal1zMzM3HZZZcZ1CtzUlpaip07dyI5ORlpaWlo0aJFpXE7f/481qxZE/PjpmRsevfuDYfDUek9eXl5+OOPP2J+/I4fP47c3FwkJycDiK2xkiQJDz74ID799FP8+OOPSEtLq/R3OrdkahurQMTyuVUVSZJQWlpq3nNKE3fmGGXp0qWSw+GQ3G63tGPHDmnKlClS/fr1pZycHKO7ZiiPPPKItHr1amn//v3S+vXrpRtuuEFq2LBhxbjMmjVLatSokfTpp59K27Ztk+644w4pOTlZKioqMrjn2nP69GkpOztbys7OlgBIc+fOlbKzs6WDBw9KkqRsbCZMmCA5nU7p+++/lzZv3iwNGjRI6tGjh1RWVmbUf0sTahqr06dPS4888oi0bt066cCBA9KqVaukAQMGSK1atYrJsZo4caLUqFEjafXq1VJeXl7Fv+Li4or30LnFqG2s6NySefLJJ6WffvpJOnDggLR161bpqaeekmw2m7Ry5UpJksx5TpHIUZk333xTatOmjVSnTh3pkksuqRSGGKvcfvvtUnJysuRwOKSWLVtKw4cPl7Zv317xd5/PJ02fPl1q0aKFFB8fL1111VXStm3bDOyxfqxatUoCUO3fqFGjJElSNjbnzp2THnzwQalJkyZSvXr1pBtuuEE6dOiQAf8bbalprIqLi6UhQ4ZITZs2lRwOh9S6dWtp1KhR1cYhVsYq0DgBkBYvXlzxHjq3GLWNFZ1bMmPGjKm4vzVt2lS65pprKgSOJJnznBIkSZK0sRERBEEQBEEYB/nkEARBEAQRlZDIIQiCIAgiKiGRQxAEQRBEVEIihyAIgiCIqIREDkEQBEEQUQmJHIIgCIIgohISOQRBEARBRCUkcgiCIAiCiEpI5BAEQRAEEZWQyCEIwrQMHDgQU6ZMMbobBEFYFCrrQBCEKRg4cCB69uyJV199teK1EydOwOFwoGHDhrr3Z8qUKcjJycFnn32m+3cTBKEOZMkhCMK0NGnSxBCBAwBZWVno16+fId9NEIQ6kMghCMJw7rvvPqxZswavvfYaBEGAIAjIycmptlw1cOBATJ48GVOmTEHjxo3RvHlzLFy4EGfPnsXo0aPRsGFDtGvXDt9++23FZyRJwssvv4y2bduiXr166NGjBz7++OOgffF4PKhTpw7WrVuHp59+GoIgoH///lr+9wmC0AgSOQRBGM5rr72GAQMGYNy4ccjLy0NeXh5SUlICvnfJkiVISkrCb7/9hsmTJ2PixIkYMWIELrvsMmzevBlDhw7FPffcg+LiYgDAv//9byxevBgZGRnYvn07pk6dirvvvhtr1qwJeHy73Y61a9cCALZs2YK8vDysWLFCm/84QRCaQj45BEGYgkA+OVVfGzhwILxeL37++WcAgNfrRaNGjTB8+HC89957AID8/HwkJyfj119/Rbdu3ZCUlIQff/wRAwYMqDju2LFjUVxcjA8++CBgXz777DOMHTsWBQUF2vxnCYLQhTijO0AQBBEK3bt3r9i22+248MIL0a1bt4rXmjdvDgA4duwYduzYgZKSEgwePLjSMc6fP49evXoF/Y7s7Gz06NFD5Z4TBKE3JHIIgrAUDoej0r4gCJVeEwQBAODz+eDz+QAAX3/9NVq1alXpc/Hx8UG/Y8uWLSRyCCIKIJFDEIQpqFOnDrxer6rH7Ny5M+Lj43Ho0CFcffXVij+3bds23Hrrrar2hSAI/SGRQxCEKUhNTcWGDRuQk5ODBg0aoEmTJhEfs2HDhnj00UcxdepU+Hw+XHHFFSgqKsK6devQoEEDjBo1KuDnfD4ftm7diiNHjqB+/fpo1KhRxH0hCEJ/KLqKIAhT8Oijj8Jut6Nz585o2rQpDh06pMpxn3vuOTzzzDNIT0/HxRdfjKFDh+LLL79EWlpa0M88//zzWLZsGVq1aoVnn31WlX4QBKE/FF1FEARBEERUQpYcgiAIgiCiEhI5BEEQBEFEJSRyCIIgCIKISkjkEARBEAQRlZDIIQiCIAgiKiGRQxAEQRBEVEIihyAIgiCIqIREDkEQBEEQUQmJHIIgCIIgohISOQRBEARBRCUkcgiCIAiCiEr+P6n24C/qa92wAAAAAElFTkSuQmCC", "text/plain": [ "Figure(PyObject
)" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "C = [ 0 0 1 0\n", " 0 0 0 1\n", " -0.02 0.01 0 0\n", " 0.01 -0.02 0 0 ]\n", "Δt = 1.0\n", "T = exp(C*Δt) # propagator matrix\n", "\n", "x₀ = [0.0,0,1,0] # initial condition\n", "\n", "# loop over 300 timesteps and keep track of x₁(t)\n", "x = x₀\n", "x₁ = [ x₀[1] ]\n", "for i = 1:300\n", " x = T*x # repeatedly multiply by T\n", " push!(x₁, x[1]) # & store current x₁(t) in the array x₁\n", "end\n", "\n", "plot((0:300)*Δt, x₁, \"r.-\")\n", "xlabel(\"time \\$t\\$\")\n", "ylabel(\"solution \\$x_1(t)\\$\")\n", "grid()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "(This is **not** an approximate solution. It is the *exact* solution, up to the computer's roundoff errors, at the times $t=0,\\Delta t, 2\\Delta t, \\ldots$. Don't confuse it with approximations like [Euler's method](https://en.wikipedia.org/wiki/Euler_method).)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Stability of solutions in eᴬ vs. Aⁿ\n", "\n", "It is important to compare and contrast the two cases we have studied:\n", "\n", "* Multiplying by $A^n$ (e.g. in **linear recurrence** equations $x_{n+1} = Ax_n$) corresponds to multiplying each eigenvector by $\\lambda^n$, which:\n", " - blows up if $|\\lambda| > 1$\n", " - decays if $|\\lambda| < 1$\n", " - oscillates if $|\\lambda|=1$\n", " - If $\\lambda = 1$ you have a steady-state vector (a stable \"attractor\" if the other eigenvalue have $|\\lambda| < 1$).\n", " - large-$n$ behavior dominated by biggest $|\\lambda|$\n", "\n", "versus\n", "\n", "* Multiplying by $e^{At}$ (e.g. **linear differential** equations $dx/dt=Ax$ ), corresponds to multiplying each eigenvector by $e^{\\lambda t}$, which\n", " - blows up if $\\operatorname{Re}(\\lambda) > 0$\n", " - decays if $\\operatorname{Re}(\\lambda) < 0$\n", " - oscillates if $\\operatorname{Re}(\\lambda) = 0$ (purely imaginary λ)\n", " - If $\\lambda = 0$ you have a steady-state solution (a stable \"attractor\" if the other eigenvalue have $\\operatorname{Re}(\\lambda) < 0$). \n", " - large-$t$ behavior dominated by biggest $\\operatorname{Re}(\\lambda)$. (Note: not biggest magnitude! Remember, $0 > -1 > -2$.)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Relating $e^{At}$ and $A^n$\n", "\n", "These two cases are **related by the propagator matrix** $T = e^{A\\Delta t}$! Solving the ODE for long time, or multiplying by $e^{At}$ for large $t$, corresponds to **repeatedly multiplying** by $T$!\n", "\n", "What are the eigenvalues of $T$ for a diagonalizable $A = X \\Lambda X^{-1}$? Well, since \n", "\n", "$$\n", "T = e^{A \\Delta t} = X e^{\\Lambda \\Delta t} X^{-1} = \n", "X \\begin{pmatrix} e^{\\lambda_1 \\Delta t} & & & \\\\\n", " & e^{\\lambda_2 \\Delta t} & & \\\\\n", " & & \\ddots & \\\\\n", " & & & e^{\\lambda_m \\Delta t} \\end{pmatrix} X^{-1}\n", "$$\n", "\n", "the eigenvalues of $T$ are just $e^{\\lambda \\Delta t}$ (the equation above is precisely the diagonalization of $T$).\n", "\n", "Equivalently, for an eigenvector $\\vec{x}_k$ of $A$, $T\\vec{x}_k = e^{\\lambda_k \\Delta t} \\vec{x}_k$, so $\\vec{x}_k$ is **also an eigenvector** of $T$ with eigenvalue $e^{\\lambda_k \\Delta t}$. Let's check:" ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "5-element Vector{Float64}:\n", " 0.23567369540695432\n", " 0.2899190400271146\n", " 0.9169902702871335\n", " 2.0483367693538845\n", " 6.99042092346493" ] }, "execution_count": 20, "metadata": {}, "output_type": "execute_result" } ], "source": [ "eigvals(exp(A*Δt))" ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "5-element Vector{Float64}:\n", " 0.235673695406957\n", " 0.28991904002711205\n", " 0.916990270287135\n", " 2.0483367693538863\n", " 6.990420923464918" ] }, "execution_count": 21, "metadata": {}, "output_type": "execute_result" } ], "source": [ "λ = eigvals(A)\n", "exp.(λ * Δt)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Yup, they match (although the order is different: Julia gives the eigenvalues in a somewhat \"random\" order)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "What does this mean for stability of the solutions?\n", "\n", "For example, if $A$ has an real eigenvalue with $\\lambda < 0$, a decaying solution, then $T$ has an eigenvalue $e^{\\lambda \\Delta t} < 1$, which is also decaying when you multiply by $T$ repeatedly!\n", "\n", "It is easy to verify that going from $\\lambda \\to e^\\lambda$ turns the **conditions for growing/decaying ODE (eᴬᵗ) solutions into the rules for growing/decaying Aⁿ solutions!**." ] } ], "metadata": { "@webio": { "lastCommId": null, "lastKernelId": null }, "kernelspec": { "display_name": "Julia 1.7.1", "language": "julia", "name": "julia-1.7" }, "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", "version": "1.7.1" } }, "nbformat": 4, "nbformat_minor": 2 }