M368K Workbook - Fixed Point algorithms for linear systems.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Fixed point iterations\n",
    "\n",
    "by [Richard Tsai](https://www.researchgate.net/profile/Richard_Tsai3)\n",
    "\n",
    "$f:\\mathbb{R}^d\\mapsto\\mathbb{R}^d$\n",
    "\n",
    "* Need to come up with $g(x)$ such that \n",
    "$$ f(x^*)=0 \\iff x^* = g(x^*) $$\n",
    "\n",
    "* Introduce the fixed point iteration for $g$: $$ x_{k+1} := g(x_k).$$\n",
    "\n",
    "**What should $g(x)$ be?** \n",
    "So that $x_k\\longrightarrow x^*$\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Linear systems of equations\n",
    "Let $$f(x)=Ax-b,$$\n",
    "where $A$ is an $d\\times d$ invertible matrix, and $x,b$ are two vectors in $\\mathbb{R}^d$.\n",
    "We want to find the solution of $f(x)=0$.\n",
    "\n",
    "Typically, we write the above linear system of equations as \n",
    "$$ Ax =b $$\n",
    "\n",
    "We will consider fixed point iterations with $g(x)=Gx+c$.\n",
    "\n",
    "### General strategy for creating fixed-point iterations for systems by splitting the matrix\n",
    "\n",
    "A general strategy in coming up with fixed point iterations is to *split* $A$ into the sum of two matrices, separate the equations accordingly by the equal sign\n",
    "$$\n",
    "A=M-N \\implies Ax=b \\iff Mx=Nx+b\n",
    "$$\n",
    "and define the iterations\n",
    "$$\n",
    "Mx^{(k+1)} = Nx^{(k)} +b, \n",
    "$$\n",
    "\n",
    "leading to the general form\n",
    "\n",
    "$$x^{(k+1)}=M^{-1}Nx^{(k)}+M^{-1}b \\iff x^{(k+1)}=G x^{(k)}+ c $$\n",
    "\n",
    "**Of course, $M$ has to be very simple and easy to invert.**\n",
    "\n",
    "**Both iterative schemes can be used to solve *nonlinear* equations**\n",
    "\n",
    "** Learniing objective **:\n",
    "* Know how convergence of a fixed point iterations is related to the eigenvalues of the matrices involved\n",
    "\n",
    "\n",
    "## Three fixed point iterations to solve linear systems of equations\n",
    "\n",
    "We first discuss two particular fixed-point iterations (two different $g$'s):\n",
    "\n",
    "* Jacobi iterations (J)\n",
    "* Gauss-Seidel iterations (G-S)\n",
    "* Successive Over Relaxations (SOR): this method builds on top of (G-S)\n",
    "\n",
    "**Learning objectives**:\n",
    "\n",
    "* Know how to write Julia code for the two methods\n",
    "* Know the general convergence theory. \n",
    "    * For this class: Burden and Faires: Lemma 7.18, Theorem 7.19, Theorem 7.21\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Jacobi iterations\n",
    "\n",
    "In the Jacobi iterations, $M=D$ is the diagonal of $A$ (and $N$ is $A$ with its diagonal culled out): $$G=G_J:=D^{-1}(L+U)$$ \n",
    "\n",
    "$$x^{(k+1)}=G_{J} x^{(k)}+ c_{J}$$\n",
    "\n",
    "* **Very easy parallelization**\n",
    "\n",
    "\n",
    "$$ \n",
    "a_{11} \\mathbf{x_1^{(k+1)}} + a_{12} x_2^{(k)}+a_{13} x_3^{(k)}+...+a_{1n} x_n^{(k)}=b_1 \n",
    "$$\n",
    "\n",
    "$$ \n",
    "a_{21} x_1^{(k)} + a_{22} \\mathbf{x_2^{(k+1)}}+a_{13} x_3^{(k)}+...+a_{1n} x_n^{(k)}=b_2 \n",
    "$$\n",
    "\n",
    "$$\\vdots$$\n",
    "\n",
    "\n",
    "Solving nonlinear equations: $ f_1(x) = 0,f_2(x)=0,...,f_n(x)=0, x=(x_1,x_2,...,x_n) $\n",
    "\n",
    "$$ f_1(\\mathbf{x_1^{(k+1)}}, x_2^{(k)}, ..., x_n^{(k)})=0 $$\n",
    "$$ f_2({x_1^{(k)}}, \\mathbf{x_2^{(k+1)}}, ..., x_n^{(k)})=0 $$\n",
    "$$\\vdots$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[-0.161604, -0.293657, -0.049588]\n",
      "[0.173608, 0.129596, 0.235815]\n",
      "[0.00333625, -0.0739547, -0.0658961]\n",
      "[0.13456, 0.0921354, 0.0825413]\n",
      "[0.0568105, -0.00223245, -0.0357369]\n",
      "[0.111168, 0.065812, 0.0326414]\n",
      "[0.0770452, 0.0239654, -0.0160352]\n",
      "[0.100017, 0.0525354, 0.0141558]\n",
      "[0.0852462, 0.034333, -0.00633711]\n",
      "[0.0950345, 0.0464692, 0.00676908]\n",
      "[0.0886754, 0.0386165, -0.00194712]\n",
      "[0.0928608, 0.0437987, 0.0037021]\n",
      "[0.0901296, 0.0404229, -2.18793e-5]\n",
      "[0.091922, 0.0426409, 0.00240575]\n",
      "[0.0907501, 0.0411918, 0.000811544]\n",
      "[0.0915182, 0.042142, 0.0018534]\n",
      "[0.0910156, 0.0415205, 0.00117035]\n",
      "[0.0913448, 0.0419277, 0.00161722]\n",
      "[0.0911293, 0.0416612, 0.00132446]\n",
      "[0.0912704, 0.0418357, 0.00151609]\n",
      "[0.091178, 0.0417215, 0.00139058]\n",
      "[0.0912385, 0.0417963, 0.00147275]\n",
      "[0.0911989, 0.0417473, 0.00141894]\n",
      "[0.0912249, 0.0417794, 0.00145417]\n",
      "[0.0912079, 0.0417584, 0.0014311]\n",
      "\n",
      "Checking the residual: [-0.000111221, -0.00020635, -0.000286973]\n"
     ]
    }
   ],
   "source": [
    "#Example: Jacobi iterations\n",
    "\n",
    "A=[10. 2. 3.; 4. 15. 6.; 7. 8. 19.]\n",
    "\n",
    "#What happens if the diagonal elements are smaller, but A is still invertible?\n",
    "#A=[1. 2. 3.; 4. 5. 6.; 7. 8. 19.]\n",
    "\n",
    "\n",
    "M=zeros(3,3)\n",
    "\n",
    "M=diagm( diag(A),0)\n",
    "\n",
    "N=M-A\n",
    "\n",
    "Minv = inv(M)\n",
    "\n",
    "b=ones(3)\n",
    "x=rand(3)\n",
    "\n",
    "#the following for-loop is only for illustrating teh concept of Jacobi iteration\n",
    "#in practice, one should never allocate memory for M or Minv\n",
    "for k=1:25\n",
    "   x = Minv*(N*x+b) \n",
    "   println(x)\n",
    "end\n",
    "\n",
    "println(\"\\nChecking the residual: \", A*x-b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The eigenvalues of Minv*N is important. \n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "3-element Array{Float64,1}:\n",
       " -0.65478 \n",
       "  0.223815\n",
       "  0.430965"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "println(\"The eigenvalues of Minv*N is important. \")\n",
    "eigvals(Minv*N)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "JacobiIterations! (generic function with 1 method)"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#a more appropirate way to implementing the JAcobi iterations\n",
    "function JacobiIterations!(x, A, b, numitrs)\n",
    "   \n",
    "   xnext=copy(x)\n",
    "   sum=0.0\n",
    "    \n",
    "   for itr=1:numitrs \n",
    "   for i=1:size(A,1)\n",
    "        sum=0.0\n",
    "        for j in [1:i-1; i+1:size(A,2)]\n",
    "            sum+=A[i,j]*x[j]\n",
    "        end\n",
    "        xnext[i] = (b[i]-sum)/A[i,i]\n",
    "   end\n",
    "    \n",
    "   x.=xnext\n",
    "   #println(itr,\": \", x)\n",
    "   end \n",
    "end\n",
    "\n",
    "\n",
    "function JacobiIterations!(x, A, b, tol; max_itrs=100000)\n",
    "   \n",
    "   xnext=rand(size(x))\n",
    "   sum=0.0\n",
    "\n",
    "   itr=0\n",
    "    \n",
    "   while norm(xnext-x)>tol && itr<max_itrs\n",
    "        x.=xnext\n",
    "        \n",
    "    for i=1:size(A,1)\n",
    "        sum=0.0\n",
    "        for j in [1:i-1; i+1:size(A,2)]\n",
    "            sum+=A[i,j]*x[j]\n",
    "        end\n",
    "        xnext[i] = (b[i]-sum)/A[i,i]\n",
    "    end\n",
    "   itr+=1\n",
    "   end \n",
    "    \n",
    "   println(\"Jacobi finishes in \",itr,\" iterations.\")\n",
    "\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Jacobi finishes in 0 iterations.\n",
      "[0.0306215, 0.98599, 0.363916]\n",
      "[2.36994, 16.0958, 14.0167]\n"
     ]
    }
   ],
   "source": [
    "x=rand(3)\n",
    "\n",
    "JacobiIterations!(x, A, b, 25)\n",
    "\n",
    "println(x)\n",
    "println(A*x-b)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Gauss-Seidel iterations\n",
    "\n",
    "In the Gauss-Seidel iterations, $M$ is the lower triangular part of $A$ (including the diagonal)\n",
    "\n",
    "$$ A= D-L-U, M=D-L, N=U \\implies G_{GS}:=(D-L)^{-1}U$$\n",
    "\n",
    "$$x^{(k+1)}=G_{GS} x^{(k)}+ c_{GS}$$\n",
    "\n",
    "* This algorithm is of **sequential nature**:\n",
    "    * Solve for first variable, $x_1$ in the first equation, using the previous guess for all other variables    \n",
    "$$ \n",
    "a_{11} \\mathbf{x_1} + a_{12} x_2^{(k)}+a_{13} x_3^{(k)}+...+a_{1n} x_n^{(k)}=b_1 \n",
    "$$\n",
    "    Obtain the value of $\\mathbf{x_1^{(k+1)}}$.\n",
    "    \n",
    "    * Solve for the second variable, $x_2$, in the second equation, using the value $\\mathbf{x_1^{(k+1)}}$, and the previous guess for all other variables\n",
    "$$ \n",
    "a_{21} \\mathbf{x_1^{(k+1)}} + a_{22} \\mathbf{x_2}+a_{13} x_3^{(k)}+...+a_{1n} x_n^{(k)}=b_2 \n",
    "$$\n",
    "    Obtain the value of $\\mathbf{x_2^{(k+1)}}$.\n",
    "$$ \\vdots$$\n",
    "\n",
    "Solving nonlinear equations: $ f_1(x) = 0,f_2(x)=0,...,f_n(x)=0, x=(x_1,x_2,...,x_n) $\n",
    "\n",
    "$$ f_1(\\mathbf{x_1^{(k+1)}}, x_2^{(k)}, ..., x_n^{(k)})=0 $$\n",
    "$$ f_2(\\mathbf{x_1^{(k+1)}}, \\mathbf{x_2^{(k+1)}}, ..., x_n^{(k)})=0 $$\n",
    "$$\\vdots$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "GaussSeidelIterations! (generic function with 1 method)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function GaussSeidelIterations!(x, A, b, tol; max_itrs=100000)\n",
    "   \n",
    "   sum=0.0\n",
    "   itr=0\n",
    "   xold=rand(size(x))\n",
    "    \n",
    "   while norm(xold-x)>tol && itr<max_itrs #for itr=1:numitrs \n",
    "    xold.=x\n",
    "    for i=1:size(A,1)\n",
    "        sum=0.0\n",
    "            for j in [1:i-1; i+1:size(A,2)]\n",
    "            sum+=A[i,j]*x[j]\n",
    "        end\n",
    "        x[i] = (b[i]-sum)/A[i,i]\n",
    "    end\n",
    "    itr+=1\n",
    "   end\n",
    "   \n",
    "   println(\"Gauss-Seidel finishes in \",itr,\" iterations.\")\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Gauss-Seidel is faster than Jacobi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Gauss-Seidel finishes in 10 iterations.\n",
      "[0.0912146, 0.0417667, 0.00144021]\n",
      "[-1.38889e-7, 6.43824e-7, 0.0]\n",
      "\n",
      "Jacobi finishes in 36 iterations.\n",
      "[0.0912144, 0.0417664, 0.00143994]\n",
      "[-3.51563e-6, -6.52255e-6, -9.07119e-6]\n"
     ]
    }
   ],
   "source": [
    "x0=rand(3)\n",
    "\n",
    "x=copy(x0)\n",
    "GaussSeidelIterations!(x, A, b, 0.000001)\n",
    "\n",
    "println(x)\n",
    "println(A*x-b,\"\\n\")\n",
    "\n",
    "x=copy(x0)\n",
    "JacobiIterations!(x, A, b, 0.000001)\n",
    "\n",
    "println(x)\n",
    "println(A*x-b)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Successive over relaxation (SOR)\n",
    "\n",
    "SOR is obtained by *averaging* the Gauss-Seidel update, $x^{(k+1)}_{GS}$, with the approximation computed in the previous iteration, $x^{(k)}$:\n",
    "\n",
    "$$ x^{(k+1)}_{SOR}:= \\omega x^{(k+1)}_{GS}+(1-\\omega) x^{(k)}. $$\n",
    "\n",
    "\n",
    "Three important theorems about SOR:\n",
    "\n",
    "Theorem (Kahan): Let $a_{ii}\\neq 0, i=1,2,\\cdots,n$. **If SOR converges, then $0<\\omega<2$.**\n",
    "\n",
    "Theorem (Ostrowski-Reich): **SOR converges for all $0<\\omega<2$ iff $A$ is invertible, symmetric, positive definite.**\n",
    "\n",
    "**Choose $\\omega$ to make $\\rho(G_{SOR})$ as small as possible.**\n",
    "\n",
    "* This was pioneered by David Young, here at UT Ausitn in 1950'\n",
    "\n",
    "Theorem (Young): $\\omega_{opt}$ as a function of $\\rho(G_J)$.\n",
    "\n",
    "For more in-depth exposition about the choice of $\\omega$, check out GGK, Sec 11.3.4.\n",
    "\n",
    "## Symmetric Successive over relaxation (SSOR)\n",
    "\n",
    "Involves two SOR updates:\n",
    "\n",
    "* One SOR step in which the G-S update involves solving each row of the linear system in **increasing order** (row 1, than row 2,..., row n). \n",
    "    * In other words, the G-S iteration is defined by $A=D-L-U=M-N$, with $M=D-L$.\n",
    "\n",
    "* One SOR step in which the G-S update involves solving each row of the linear system in **decreasing order** (row n, than row (n-1),..., row 1).\n",
    "    * In other words, the G-S iteration is defined by $A=D-L-U=M-N$, with $M=D-U$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "SOR! (generic function with 1 method)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function SOR!(x, A, b, ω, tol; max_itrs=100000)\n",
    "   \n",
    "   sum=0.0\n",
    "   itr=0\n",
    "   xold=rand(size(x))\n",
    "    \n",
    "   diff=xold-x\n",
    "    \n",
    "   while norm(diff)>tol && itr<max_itrs \n",
    "    xold.=x\n",
    "    for i=1:size(A,1)\n",
    "        sum=0.0\n",
    "        for j=1:i-1 \n",
    "            sum+=A[i,j]*x[j]\n",
    "        end\n",
    "        for j=i+1:size(A,2) \n",
    "            sum+=A[i,j]*x[j]\n",
    "        end\n",
    "\n",
    "        x[i] = ω*(b[i]-sum)/A[i,i]+(1.0-ω)*xold[i]\n",
    "    end\n",
    "    diff.=xold.-x\n",
    "    itr+=1\n",
    "   end\n",
    "   \n",
    "   println(\"SOR finishes in \",itr,\" iterations.\")\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Gauss-Seidel finishes in 24 iterations.\n",
      "[0.171247, 0.0528541, -0.0443975]\n",
      "[-7.40741e-13, -2.06857e-12, 0.0]\n",
      "\n",
      "Jacobi finishes in 485 iterations.\n",
      "[0.171247, 0.0528541, -0.0443975]\n",
      "[3.15925e-12, 6.17173e-12, 8.66729e-12]\n",
      "\n",
      "SOR finishes in 50 iterations.\n",
      "  0.101188 seconds (36.02 k allocations: 1.548 MiB)\n",
      "[0.171247, 0.0528541, -0.0443975]\n",
      "[-9.62452e-13, -6.09579e-12, 2.88547e-12]\n",
      "\n"
     ]
    }
   ],
   "source": [
    "#Example: Jacobi iterations\n",
    "\n",
    "A=[6. 2. 3.; 4. 11. 6.; 7. 8. 14.]\n",
    "\n",
    "\n",
    "M=zeros(3,3)\n",
    "\n",
    "M=diagm( diag(A),0)\n",
    "\n",
    "N=M-A\n",
    "\n",
    "Minv = inv(M)\n",
    "\n",
    "tol=1e-12\n",
    "\n",
    "b=ones(3)\n",
    "x=rand(3)\n",
    "\n",
    "x0=rand(3)\n",
    "\n",
    "x=copy(x0)\n",
    "GaussSeidelIterations!(x, A, b, tol)\n",
    "\n",
    "println(x)\n",
    "println(A*x-b,\"\\n\")\n",
    "\n",
    "x=copy(x0)\n",
    "JacobiIterations!(x, A, b, tol)\n",
    "\n",
    "println(x)\n",
    "println(A*x-b,\"\\n\")\n",
    "\n",
    "#SOR with non-optimal weight\n",
    "x=copy(x0)\n",
    "@time SOR!(x, A, b, 0.75, tol)\n",
    "\n",
    "println(x)\n",
    "println(A*x-b,\"\\n\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example: solving a two-point boundary value problem for Poisson's equation\n",
    "\n",
    "$$ -u_{xx} =f(x),~~~0<x<1,~~~u(0)=u(1)=0 $$\n",
    "\n",
    "* Discretize the unit interval: $x_j = j\\Delta x, j=0,1,2,\\cdots,N,\\Delta x =1/N$\n",
    "\n",
    "* Approximate the differential equation on the grid nodes:\n",
    "    $$ -u_{xx} (x_j)=f(x_j) \\approx -(u_{j+1}-2u_j+u_{j-1})/\\Delta x^2 = f(x_j),~~~j=1,2,\\cdots,N-1 $$\n",
    "\n",
    "* Set the boundary conditions: $u_0=0=u_N$.\n",
    "\n",
    "This leads to a linear system to solve: $AU=F$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  0.682930 seconds (257.98 k allocations: 13.603 MiB, 0.87% gc time)\n",
      "\n",
      "\n",
      "Gauss-Seidel finishes in 4494 iterations.\n",
      "  0.210010 seconds (683.78 k allocations: 131.004 MiB, 6.99% gc time)\n",
      "\n",
      "\n",
      "Jacobi finishes in 11250 iterations.\n",
      "  0.405807 seconds (1.68 M allocations: 326.109 MiB, 3.64% gc time)\n",
      "\n",
      "\n",
      "SOR finishes in 1580 iterations.\n",
      "  0.125506 seconds (26.05 k allocations: 1.146 MiB)\n",
      "λ_J = 0.9980267284282716\n",
      "\n",
      "optimal ω for this case is 1.8818383898322273\n",
      "\n",
      "SOR finishes in 1123 iterations.\n",
      "  0.032533 seconds (61 allocations: 2.938 KiB)\n"
     ]
    }
   ],
   "source": [
    "N=50\n",
    "xgrid=linspace(0,1,N+1)\n",
    "Δx=1/N\n",
    "\n",
    "f=x-> 1.0 #x.^2+0.0*sin(20*x)\n",
    "\n",
    "#The vector U stores (u1,u2,u3,cdots, u_(N-1))\n",
    "U=zeros(N-1)\n",
    "\n",
    "#N-2 because we exclude x=0 and X=1.\n",
    "A=spdiagm( (-ones(N-2), 2.0*ones(N-1),-ones(N-2)), [-1,0,1])\n",
    "#A=-2.*diagm(ones(N-1))+diagm(ones(N-2),1)+diagm(ones(N-2),-1)\n",
    "F=collect(f.(xgrid[2:end-1]))\n",
    "\n",
    "b=F.*Δx^2\n",
    "\n",
    "@time U=-A\\b\n",
    "\n",
    "Ugs=copy(U)\n",
    "Uj=copy(U)\n",
    "Usor=copy(U)\n",
    "\n",
    "println(\"\\n\")\n",
    "@time GaussSeidelIterations!(Ugs, A, b, 1e-10)\n",
    "\n",
    "println(\"\\n\")\n",
    "@time JacobiIterations!(Uj, A, b, 1e-10)\n",
    "\n",
    "println(\"\\n\")\n",
    "@time SOR!(Usor, A, b, 1.5, 1e-10)\n",
    "\n",
    "M=diagm(diag(A),0)\n",
    "N=M-A\n",
    "λ=svdvals(M\\N)\n",
    "ωopt=2.0/(1.0+sqrt(1.0-λ[1]^2))\n",
    "\n",
    "println(\"λ_J = \", λ[1],\"\\n\")\n",
    "println(\"optimal ω for this case is \", ωopt,\"\\n\")\n",
    "Usor=copy(U)\n",
    "@time SOR!(Usor, A, b, ωopt-0.2, tol)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\u001b[1m\u001b[36mINFO: \u001b[39m\u001b[22m\u001b[36mRecompiling stale cache file /Users/ytsai/.julia/lib/v0.6/PyPlot.ji for module PyPlot.\n",
      "\u001b[39m"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjsAAAGgCAYAAABMn6ZGAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzs3Xt8U/X9P/BXkqbphTa9XymUa1ug5U5bQAWEgoCAOrlJVcTqNnFj0/0m00237zY2neI2nfdNBRQFRUWwioDKpRcKLW2hLddC7xfaJr0maXJ+f6QNIKVQaPrJ5fV8PPKHbU54BZC8es7nvD8ySZIkEBERETkouegARERERNbEskNEREQOjWWHiIiIHBrLDhERETk0lh0iIiJyaCw7RERE5NBYdoiIiMihsewQERGRQ2PZISIiIofGskNEREQOjWWHiIiIHJqL6AAimEwmlJeXw8vLCzKZTHQcIiIiug6SJKGxsRFhYWGQy6//fI1Tlp3y8nJERESIjkFEREQ3oKSkBP3797/u5ztl2fHy8gJg/s3y9vYWnIaIiIiuh1arRUREhOVz/Ho5ZdnpvHTl7e3NskNERGRneroEhQuUiYiIyKGx7BAREZFDY9khIiIih8ayQ0RERA6NZYeIiIgcGssOEREROTSWHSIiInJoLDtERETk0Fh2iIiIyKGx7BAREZFDY9khIiIih8ayQ0RERA6NZacX5Zdp8IfP8/FZdpnoKERERELkljbgl5uz8f2JGtFRLFh2etHB07V4P+0cPsw8LzoKERGREJ9ll+PznHJ8eqRUdBQLlp1eNC8uDACQWVyHKm2b4DRERER9y2iS8GVuOQDgzo7PRFvAstOLwn3cMWGgLyQJ2JFbIToOERFRnzpUXIfqRh283Vxwy/AA0XEsWHZ62fy4UADA9o5mS0RE5Cy2HzV/9s0eGQKVi0JwmotYdnrZ3LhQyGVA9vkGlNS1iI5DRETUJ9qNJnyVXwkAuHO07VzCAlh2el2QlxsSBvsDAHbk8VIWERE5h4OnL6CuWQ8/T1dMHuIvOs5lWHasoLPRdp7OIyIicnSdn3lzY0PgorCtemFbaRzEnJEhcJHLcKxcizM1TaLjEBERWZWu3Yivj5kvYc23obuwOrHsWIGvpyumDjOvQv+Sd2UREZGD23eiFtq2dgR7qzAx0k90nCuw7FhJ53yBL46WQ5IkwWmIiIisp/MO5HmxYVDIZYLTXIllx0pmjQyGq0KOU9VNKKpqFB2HiIjIKlr1Rnx7vAoAMH90qOA0XWPZsRJvNyWmRQUC4EJlIiJyXHuLqtGsNyLcxx1jI3xEx+kSy44Vdd6V9WVuBS9lERGRQ+r8gf7O0WGQyWzvEhbAsmNVt8cEwV2pwLkLLcgr04iOQ0RE1KuadO3YU1gNALjTRi9hASw7VuXh6oLbY4IA8FIWERE5nm+PV0HXbsLgAE+MCPUWHeeqWHasrPNS1o7cCphMvJRFRESOo/MH+fk2fAkLYNmxutuGB8JL5YJyTRuOnK8XHYeIiKhXaFoM+OFkDQDgzjjbvYQFsOxYnZtSgVkjgwHwUhYRETmOr49VwmCUEB3ihWHBXqLjdItlpw9YLmXlVcLIS1lEROQAOgcJ2toO511h2ekDU4cGwMdDidomHTLOXBAdh4iI6KbUNulw4FQtAGC+jV/CAlh2+oRSIccdo8x/GTqbMBERkb36Kr8SJgkY3V+Ngf6eouNcE8tOH+lcvPVVfiX07SbBaYiIiG6c5S4sG9zhvCssO30kfrA/Avqp0NBisJz6IyIisjeVmjYcKq4DAMyzg0tYAMtOn1HIZZbrmryURURE9mpHXgUkCZgY6YswH3fRca4Ly04f6iw73xyrQpvBKDgNERFRz9nbJSyAZadPjRvgizC1G5p07fiuqEZ0HCIioh4pqWtBTkkD5DLgjtgQ0XGuG8tOH5LLZZhv2Qmdl7KIiMi+dC7DSBzijyAvN8Fprp9Vy059fT2Sk5OhVquhVquRnJyMhoaGbo/R6XR4/PHHERAQAE9PTyxYsAClpaWW7x89ehTLli1DREQE3N3dERMTg3/+85/WfBu96s6O037fFlShWdcuOA0REdH12360AsDFzzJ7YdWys3z5cuTk5CA1NRWpqanIyclBcnJyt8esWbMG27Ztw+bNm7F//340NTVh/vz5MBrNa1wOHz6MwMBAbNy4EceOHcPTTz+NtWvX4pVXXrHmW+k1o8K9MSjAE20GE745Xik6DhER0XU5UdWIggotlAoZ5oyyn0tYACCTJMkq+xcUFBRgxIgRSE9PR3x8PAAgPT0diYmJKCwsRFRU1BXHaDQaBAYGYsOGDViyZAkAoLy8HBEREdi5cydmz57d5a/12GOPoaCgAHv27Ony+zqdDjqdzvLfWq0WERER0Gg08Pbu+y3pX/72BF7+9iRuGx6I9x6a1Oe/PhERUU89n1qI/3x3GrNGBOOt+ycIyaDVaqFWq3v8+W21MztpaWlQq9WWogMACQkJUKvVOHjwYJfHHD58GAaDAUlJSZavhYWFYdSoUVc9BjCXJD8/v6t+f926dZZLaWq1GhERETfwjnrPojHhAIB9J2tQ06i7xrOJiIjEMpkkfJ5jXq/T+RlmT6xWdiorKxEUFHTF14OCglBZ2fXlm8rKSri6usLX1/eyrwcHB1/1mLS0NHz88cd49NFHr5pl7dq10Gg0lkdJSUkP3knviwzwxJgIH5gk7oRORES271BxHcoaWuGlcsHtMVd+ttu6Hped5557DjKZrNtHVlYWAEAmk11xvCRJXX69O1c75tixY1i4cCH+8Ic/YNasWVc9XqVSwdvb+7KHaHeNNTfjz3PKBCchIiLq3mcdZ3XuiA2Bm1IhOE3PufT0gNWrV2Pp0qXdPicyMhK5ubmoqqq64ns1NTUIDg7u8riQkBDo9XrU19dfdnanuroakydPvuy5x48fx4wZM5CSkoJnnnmmp29DuHlxofjTl8dxtFSDMzVNGBzYT3QkIiKiK+jajdiRa7+XsIAbOLMTEBCA6Ojobh9ubm5ITEyERqNBZmam5diMjAxoNJorikun8ePHQ6lUYteuXZavVVRUID8//7Jjjh07hunTp+OBBx7AX/7yl56+BZsQ0E+FW4cFALjYmImIiGzNd0U10La1I8TbDfGD/UXHuSFWW7MTExODOXPmICUlBenp6UhPT0dKSgrmz59vuROrrKwM0dHRlkKkVquxatUqPPHEE9i9ezeys7OxYsUKxMbGYubMmQAuFp1Zs2bh17/+NSorK1FZWYmaGvubSLyo41LWZ9llsNJNcURERDfls2zzcosFY8KgkPdsGYqtsOqcnU2bNiE2NhZJSUlISkpCXFwcNmzYYPm+wWBAUVERWlpaLF9bv349Fi1ahMWLF2PKlCnw8PDA9u3boVCYrxFu2bIFNTU12LRpE0JDQy2PiRMnWvOtWMWsEcHwcFXgfF0LjpzvftgiERFRX9O0GrC7oBqA/V7CAqw4Z8eW3eh9+tbw649y8Gl2Ge5PHIg/LRwlNAsREdGlPjp0Hr/9JA9RwV5IXXNLj28w6m02N2eHrs/CjktZ24+Ww2A0CU5DRER00baOS1gLx4YJLzo3g2VHsClD/BHQT4X6FgN+OGF/646IiMgxlTe0Iv1MHQBgoR1fwgJYdoRzUcixoGMndN6VRUREtuKLjqG38YP8EO7jLjjNzWHZsQGLxprLzjfHKtHYZhCchoiI6OJdWJ13Dtszlh0bEBuuxuBAT+jaTfj62JWDGImIiPpSQYUWhZWNcFXIMXdUqOg4N41lxwbIZDLLLX3cPoKIiET7rOOzaHp0INQeSsFpbh7Ljo3oLDsHTtWiWtsmOA0RETkrk0nCFx1rSO9ygEtYAMuOzRjg74HxA31hki4uCiMiIuprGWfrUKFpg7ebC6ZF2d8O511h2bEhi8Z03pXFS1lERCRG58LkubGhdrnDeVdYdmzIvLgwuMhlyC/T4lR1o+g4RETkZNoMRuzMqwDgGHdhdWLZsSF+nq64bXggAOCzbF7KIiKivrW3sBqNunaEqd0wKdJPdJxew7JjYyw7oeeUwWRyum3LiIhIoG2WHc7DIbfTHc67wrJjY2bGBKOfygWl9a04fL5edBwiInISDS167C0y73DuKHdhdWLZsTHurgrMHhkC4OIiMSIiImvbmVcJg1FCdIgXokK8RMfpVSw7NqizUX+ZWwF9O3dCJyIi6+v8AdvRzuoALDs2KXGIP4K8VNC0GvBdxylFIiIiaymtb0FmcR1kMmBBxxgUR8KyY4MUcpllJ/RPj/BSFhERWde2js+ahEH+CFXb9w7nXWHZsVH3jO8PANhdWIW6Zr3gNERE5KgkScLWI6UAgJ90fPY4GpYdGxUT6o1R4d4wGCV8wYnKRERkJVnn6nHuQgs8XRW4IzZEdByrYNmxYT8ZZ27YnY2biIiot23NMn/GzI0NhYeri+A01sGyY8MWjAmHUmHePqKwUis6DhEROZhWvRE7OraHcNRLWADLjk3z83TFjGjzjrOfHObZHSIi6l1fH6tEk64dEX7umOhA20P8GMuOjfvJ+AgAwLbschiMnLlDRES9Z2vHD9L3jOvvUNtD/BjLjo2bFhUIf09X1Dbp8MOJGtFxiIjIQZQ1tOLA6VoA5rLjyFh2bJxSIbdsDrqVl7KIiKiXbDtSCkkCEgb7IcLPQ3Qcq2LZsQOdi8a+LahCPWfuEBHRTZIkyfIDdOdyCUfGsmMHYkK9MTKsY+bO0XLRcYiIyM4dPleP4gst8HBV4I5Rjjlb51IsO3ai8+zOJ5y5Q0REN6nzs2RubCg8VY45W+dSLDt2YsHoMLjIZcgt1aCoslF0HCIislOteiO+PGqerePoC5M7sezYCf9+qoszd3h2h4iIbtA3xyvRqGtHf193xA9y3Nk6l2LZsSOdl7I+PVKGds7cISKiG+Ass3UuxbJjR6ZHB12cuXOSM3eIiKhnyhtasf+Uc8zWuRTLjh1RKuRYOIYzd4iI6MZsyy6DJAHxg/wwwN+xZ+tcimXHzlhm7hyvRkMLZ+4QEdH1uXy2jvOc1QFYduzOiDBvjAj1ht5ownbO3CEiout05HwDztY2w8NVgbmxoaLj9CmWHTvU2ch5KYuIiK5X52fGHaOcY7bOpVh27NDCMeaZO0dLNThRxZk7RETUvTaDEV92XA24Z3y44DR9j2XHDvn3U2F658wdnt0hIqJr+PqYebZOuI87Egb5i47T51h27JRl5k42Z+4QEVH3LLN1xjvPbJ1LsezYqelRQfDzdEVNow77TtaKjkNERDaqQnPpbB3nu4QFsOzYLVcXORaOCQPAhcpERHR1nx4xz9aZNMgPA/09RccRgmXHjnVeytp1vAp1zZy5Q0REl5MkCVuySgA432ydS7Hs2LGRYWrEhquhN5rwKTcHJSKiH0k/U4fiCy3op3LBPCebrXMplh07t2RiBABg86ESSJIkOA0REdmSzYfOAwDuHB3mdLN1LsWyY+cWjgmDu1KBU9VNOHK+XnQcIiKyEQ0tenyVXwkAWDYpQnAasVh27JyXmxLz4synJj/MLBGchoiIbMW27DLo202ICfVGbLhadByhWHYcQGdj35FbAW2bQXAaIiISTZIkbO74AXjZpAjIZM43W+dSVi079fX1SE5OhlqthlqtRnJyMhoaGro9RqfT4fHHH0dAQAA8PT2xYMEClJZ2vfj2woUL6N+/P2Qy2TVf15GNG+CLoUH90Gow4oscbg5KROTsckoaUFTVCJWLHAtHO+dsnUtZtewsX74cOTk5SE1NRWpqKnJycpCcnNztMWvWrMG2bduwefNm7N+/H01NTZg/fz6MRuMVz121ahXi4uKsFd9uyGQyLO1YqPzRIV7KIiJydp1ndebFhkLtoRScRjyrlZ2CggKkpqbi7bffRmJiIhITE/HWW2/hyy+/RFFRUZfHaDQavPPOO3jxxRcxc+ZMjB07Fhs3bkReXh6+/fbby5772muvoaGhAU8++aS13oJduXtcf7gq5Mgr0yC/TCM6DhERCdKka8f2XPNZ/qWTBghOYxusVnbS0tKgVqsRHx9v+VpCQgLUajUOHjzY5TGHDx+GwWBAUlKS5WthYWEYNWrUZcccP34cf/rTn/D+++9DLr/2W9DpdNBqtZc9HI2fpyuSRgYDuHirIREROZ/tR8vRojdicKAnJkb6io5jE6xWdiorKxEUFHTF14OCglBZWXnVY1xdXeHre/kfTnBwsOUYnU6HZcuW4YUXXsCAAdfXWNetW2dZN6RWqxER4Zi34C3raPCfZ5ejVX/lZT8iInJ8mzPNP/AunciFyZ16XHaee+45yGSybh9ZWVkA0OVvsiRJPf7Nv/SYtWvXIiYmBitWrLju49euXQuNRmN5lJQ45rqWxMH+iPBzR6OuHTvyKkTHISKiPna8XIujpRooFTLcPc55t4f4sR6PU1y9ejWWLl3a7XMiIyORm5uLqqqqK75XU1OD4ODgLo8LCQmBXq9HfX39ZWd3qqurMXnyZADAnj17kJeXh61btwKAZWpwQEAAnn76afzxj3+84nVVKhVUKtX1vUE7JpfLsHTiALzwdRE+OnTeqfdBISJyRh91LGOYNSIYAf0c/3PvevW47AQEBCAgIOCaz0tMTIRGo0FmZiYmTZoEAMjIyIBGo7EUlx8bP348lEoldu3ahcWLFwMAKioqkJ+fj+effx4A8Mknn6C1tdVyzKFDh/DQQw9h3759GDJkSE/fjsP5yfj+eGnXCRwqrsep6kYMDfISHYmIiPpAm8GIbdllAIAlE7kw+VJWW7MTExODOXPmICUlBenp6UhPT0dKSgrmz5+PqKgoAEBZWRmio6ORmZkJAFCr1Vi1ahWeeOIJ7N69G9nZ2VixYgViY2Mxc+ZMAMCQIUMwatQoy2PQoEGWX6+rNULOJtjbDdOjzL8PvA2diMh5fJVfAW1bO8J93HHL0GuflHAmVp2zs2nTJsTGxiIpKQlJSUmIi4vDhg0bLN83GAwoKipCS0uL5Wvr16/HokWLsHjxYkyZMgUeHh7Yvn07FAqFNaM6lM6Jyp8cKYOunQuViYicQeeWQUsmRkAu58LkS8kkJ9wqW6vVQq1WQ6PRwNvbW3ScXtduNGHK3/egSqvDq8vHWfbOIiIix3SmpgkzXvwechlw4KkZCFW7i45kFTf6+c29sRyQi0KOxRPMZ3c4c4eIyPF1LluYFhXksEXnZrDsOKjOsrPvZC1K6lqu8WwiIrJX+nYTPjli3kNyyUTHnCN3s1h2HFSEnwduGWZeoPZxFhcqExE5qt0FVaht0iPQS4UZ0bxRpyssOw6ss+F/nFWCdqNJcBoiIrKGzR2XsH4yvj+UCn6sd4W/Kw5s1ohg+Hm6okqrw/cnakTHISKiXlZa34IfTpr/fV/KS1hXxbLjwFQuCtwzLhzAxVsSiYjIcXycVQpJAiYP8cdAf0/RcWwWy46D67yUtbeoGlXaNsFpiIiotxhNErZkXZytQ1fHsuPghgZ5YWKkL4wmCZt5doeIyGHsKaxGhaYNPh5KzB4ZIjqOTWPZcQIrEgYCAD7MPM+FykREDmJj+jkA5lEjbkruMtAdlh0nMGdUCPw9XVGpbcO3BdWi4xAR0U06d6HZcuPJffHc9PNaWHacgMpFYbme2/mTABER2a9NGebp+LcND+TC5OvAsuMklk0aAJkM2H+qFmdqmkTHISKiG9RmMFqGxXYuU6Dusew4iQg/D8yIMk/W7PyJgIiI7M+O3Ao0tBgQ7uPOicnXiWXHiaxINP8EsCWrBK16o+A0RER0IzZ0LEdYHj8ACrlMcBr7wLLjRG4bFogIP3do29qxPbdcdBwiIuqh/DINckoaoFTILBs+07Wx7DgRuVyG++LNZ3e4UJmIyP50/ts9Z1QoAr1UgtPYD5YdJ7N4QgRcXeTILdXgaEmD6DhERHSdNK0GfJZTBgBI5sLkHmHZcTJ+nq6YFxsK4OJ1XyIisn2fHC5Fm8GE4cH9MDHSV3Qcu8Ky44Q6b1XcfrQcDS16wWmIiOhaJEnCxgzzD6jJCQMhk3Fhck+w7DihcQN8MCLUG7p2E7YeLhUdh4iIriHt9AWcqWmGp6sCi8aGi45jd1h2nJBMJkNy4sWFyiaTJDgRERF1p3PZwV3jwuHlphScxv6w7DiphWPC4KVyQfGFFhw4XSs6DhERXUWVtg3fHK8CwInJN4plx0l5uLrgnvH9AQAb0rhQmYjIVn2YeR5Gk4SJkb6IDvEWHccusew4sRUJ5p1yvy2oQnlDq+A0RET0YwajCR9mmrf44VmdG8ey48SGBnkhYbAfTBKwOZP7ZRER2ZrdBVWo0urg7+mKOaNCRMexWyw7Ti45IRIA8OGhEujbTWLDEBHRZToXJi+ZGAGVi0JwGvvFsuPkkkYGI9BLhZpGHb45Xik6DhERdThd04QDpy5AJjNv+kk3jmXHySkVciybaN5MjvtlERHZjk3p5uUFM6KC0N/XQ3Aa+8ayQ1gWPwAKuQzpZ+pwsqpRdBwiIqfXom/HlsMlAIAViVyYfLNYdgihanfMjAkCwP2yiIhswRc55Whsa0eEnztuGxYoOo7dY9khAMD9iZEAgK2HS6FpNYgNQ0TkxCRJwv8OFAMAVsQPhFzOfbBuFssOAQAmD/HH8OB+aNEbsSWrRHQcIiKnlXb6AoqqGuGuVGDpRC5M7g0sOwTAvF/WyimDAADvHiyGkftlEREJ8d+Oszr3jA+H2oP7YPUGlh2yuGtsOHw9lCitb8Wujn1YiIio75y70IzdheZ/fx+cPEhwGsfBskMWbkoFlk0ynzL934GzgtMQETmfdw8WQ5KA24YHYmhQP9FxHAbLDl0mOXEgXOQyZJytw7Fyjeg4REROo7HNgC1ZpQCAh6byrE5vYtmhy4Sq3XFHbCgAWO4GICIi69uSVYomXTuGBHri1mEBouM4FJYdusLKKZEAzHMeahp1YsMQETkBo0nCe2nFAIAHpwyCTMbbzXsTyw5dYdwAX4yJ8IHeaMIHGdwNnYjI2vYUVuPchRZ4u7ngnnHhouM4HJYd6lLn2Z0N6eegazeKDUNE5OA6bwpZNmkAPFxdBKdxPCw71KW5saEI9lahtkmHHbkVouMQETmswkotDp6+AIVchvsnR4qO45BYdqhLSoXcsoXEfw+chSRxyCARkTX8b38xAGD2yGCE+7iLDeOgWHboqpZNGgCVixz5ZVpknasXHYeIyOFcaNJhW04ZAOChKbzd3FpYduiq/DxdcddY80I5DhkkIup9H2aeh77dhNhwNcYP9BUdx2Gx7FC3HuxYqJyaX4nS+haxYYiIHIjBaMKG9HMAzDeF8HZz62HZoW5Fh3hjylB/mCRgQ9o50XGIiBzGzrwKVGl1CPRSYV5cqOg4Do1lh65pZcdmdB9mnkeLvl1wGiIix9C5u/mK+IFQuSjEhnFwVi079fX1SE5OhlqthlqtRnJyMhoaGro9RqfT4fHHH0dAQAA8PT2xYMEClJaWXvG8d999F3FxcXBzc0NISAhWr15trbfh9GZEB2Ggvwe0be345EiZ6DhERHbvyPl6HC1pgKtCjvsSBoiO4/CsWnaWL1+OnJwcpKamIjU1FTk5OUhOTu72mDVr1mDbtm3YvHkz9u/fj6amJsyfPx9G48XBdi+99BKefvppPPXUUzh27Bh2796N2bNnW/OtODW5XIYHO2Y//O/AWZhMvA2diOhm/He/+aaPBWPCENBPJTiN45NJVhqgUlBQgBEjRiA9PR3x8fEAgPT0dCQmJqKwsBBRUVFXHKPRaBAYGIgNGzZgyZIlAIDy8nJERERg586dmD17Nurr6xEeHo7t27fj9ttvv6FsWq0WarUaGo0G3t7eN/4mnUiTrh2Jf92NRl073l05EdOigkRHIiKySxWaVkz9+14YTRJ2/GIqRoapRUeyGzf6+W21MztpaWlQq9WWogMACQkJUKvVOHjwYJfHHD58GAaDAUlJSZavhYWFYdSoUZZjdu3aBZPJhLKyMsTExKB///5YvHgxSkpKrppFp9NBq9Ve9qCe6adywb0TIgBwN3QiopuxIe0cjCYJ8YP8WHT6iNXKTmVlJYKCrvzpPygoCJWVlVc9xtXVFb6+l88aCA4Othxz5swZmEwm/PWvf8XLL7+MrVu3oq6uDrNmzYJer+/yddetW2dZN6RWqxEREXGT7845PTg5EjIZ8P2JGpyoahQdh4jI7jTr2rGpY4PllRwi2Gd6XHaee+45yGSybh9ZWVkA0OXMAEmSejxL4NJjTCYTDAYD/vWvf2H27NlISEjAhx9+iJMnT2Lv3r1dHr927VpoNBrLo7uzQHR1A/w9MGdkCADgzR/OCE5DRGR/Ps4qgabVgEh/D8waESw6jtPo8daqq1evxtKlS7t9TmRkJHJzc1FVVXXF92pqahAc3PUfcEhICPR6Perr6y87u1NdXY3JkycDAEJDzbMIRowYYfl+YGAgAgICcP78+S5fV6VSQaXiArDe8Mitg/FVfiU+zynDE0nDEarmPi5ERNfDYDTh7X3mhckptw6GQs4hgn2lx2d2AgICEB0d3e3Dzc0NiYmJ0Gg0yMzMtBybkZEBjUZjKS4/Nn78eCiVSuzatcvytYqKCuTn51uOmTJlCgCgqKjI8py6ujrU1tZi4MCBPX071ENjB/hi0iA/GIwS1+4QEfXAzrwKlDW0wt/TFfeM6y86jlOx2pqdmJgYzJkzBykpKUhPT0d6ejpSUlIwf/58y51YZWVliI6OthQitVqNVatW4YknnsDu3buRnZ2NFStWIDY2FjNnzgQADB8+HAsXLsQvf/lLHDx4EPn5+XjggQcQHR2N6dOnW+vt0CV+ettgAMAHGeehbTMITkNEZPskScLr35sv/z84ORJuSg4R7EtWnbOzadMmxMbGIikpCUlJSYiLi8OGDRss3zcYDCgqKkJLy8U9l9avX49FixZh8eLFmDJlCjw8PLB9+3YoFBf/Yrz//vuIj4/HvHnzcNttt0GpVCI1NRVKpdKab4c6TBsehOHB/dCka8cHGV1fOiQioov2naxFQYUW7koFkhN5FaKvWW3Oji3jnJ2bt/VwKZ7cchRBXirs++10jjphEP+iAAAgAElEQVQnIurGirczsP9ULVZOicSzd44UHcdu2dycHXJsC0aHIcTbDdWNOnyeXS46DhGRzcov02D/qVoo5DKsmsrbzUVg2aEb4uoix0NTIwEAb/xwmltIEBFdxRsdozrujAtFf18PwWmcE8sO3bBlkwbAS+WC0zXN2FNYLToOEZHNKalrwY5c89nvR24dIjiN82LZoRvm5abEfQnmhXZv/HBacBoiItvz9r4zMEnALcMCMCKMa0RFYdmhm7JySiRcFXIcKq7H4XN1ouMQEdmMumY9PsoyT+z/6W08qyMSyw7dlGBvNywaGwYAeON7biFBRNRpQ9o5tBlMGBXujclD/EXHcWosO3TTHrnVPGRwV0EVTtc0CU5DRCReq96I99KKAQCP3jqkx3tCUu9i2aGbNjTICzNjgiFJ5uvTRETObuvhEtQ16xHh5447RoWIjuP0WHaoV3RuIfHJ4TJUN7YJTkNEJI7RJOGtjg0/H546GC4KftSKxj8B6hUTIv0wboAP9EYT3uUGoUTkxFLzK3G+rgW+HkrcO4EbftoClh3qNY923G2wMf0cmnTtgtMQEfU984af5lEc9ydGwsPVRXAiAlh2qBfNignG4ABPaNvasTmTG4QSkfNJO3MBeWUauCnluJ8bftoMlh3qNXK5zHJn1jv7z8JgNAlORETUtzpHcCyeEAH/firBaagTyw71qkVjwxHQT4UKTRu+yOEGoUTkPI6Xa/H9iRrIZeaFyWQ7WHaoV7kpFZZdfV/97hSM3CCUiJzEq3tPAQDmxoZigD83/LQlLDvU65ITB0LtrsSZmmbszKsQHYeIyOpOVTdiZ77537vVM4YKTkM/xrJDva6fygUPTTGf3XllzymYeHaHiBzcK3tOQZKA2SODER3CDT9tDcsOWcWDUyLhpXJBUVUjvjleJToOEZHVnK1txhdHzWsUH58xTHAa6grLDlmF2l2JByZHAgD+veckJIlnd4jIMf1n7ymYJGBGdBBGhatFx6EusOyQ1Tw0dRA8XBU4Vq7F3qJq0XGIiHpdSV0LtmWXAeBaHVvGskNW4+fpiuQE81Ctf+0+xbM7RORwXvv+NNpNEqYODcC4Ab6i49BVsOyQVT18y2C4KeXIKWnA/lO1ouMQEfWaCk0rtmaVAgAe51kdm8ayQ1YV6KXCskkDAAD/3n1KcBoiot7zxvdnoDeaMGmQH+IH+4uOQ91g2SGre/TWIXBVyJFZXIf0MxdExyEiumnVjW34sGMPwF/wDiybx7JDVheidsO9E/oDMN+ZRURk79764Qx07SaMHeCDKUN5VsfWsexQn/jZtCFwkctw4NQFHD5XLzoOEdENq2vWY2P6xbM6MplMcCK6FpYd6hP9fT1w97hwAMArPLtDRHbsnf1n0GowIjZcjWlRgaLj0HVg2aE+8/NpQyGXAXuLapBXqhEdh4ioxzQtBrx38BwA81wdntWxDyw71GciAzyxcIz57A7X7hCRPfrfwbNo0rUjKtgLs2KCRceh68SyQ33qselDIZMB3xyvQkGFVnQcIqLr1thmwH/3nwVgPqsjl/Osjr1g2aE+NTSoH+bGhgIAXtnLuTtEZD/eTzsHbVs7Bgd6Wv4dI/vAskN9rnPS6M68CpyqbhSchojo2lr07Xin86zO9KFQ8KyOXWHZoT4XHeKNpBHBkCTglT08u0NEtm9T+nnUNesx0N8DC0aHiY5DPcSyQ0L84nbzxNHPj5bjZBXP7hCR7WrSteO1708DAH4+bQhcFPzotDf8EyMhRoWrMXuk+ezOS7tOiI5DRHRV/9t/FnXNekT6e+Ducf1Fx6EbwLJDwjyRFAWZDPgqvxL5ZZy7Q0S2R9NiwJv7zgAAfjVrOJQ8q2OX+KdGwgwP9sLCjmvf//imSHAaIqIrvbnvNBrbzHN17ozjWh17xbJDQq2ZORwKuQzfFdUgq7hOdBwiIovaJh3+d6AYAPDrpOGcq2PHWHZIqMgATyzu2BH9ha+LIEmS4ERERGb/2XsaLXojRvdXI2kEpyXbM5YdEu7xGcPgqpAj42wd9p+qFR2HiAgVmlZszDDvgWVeX8izOvaMZYeEC/Nxx30JAwAA/+DZHSKyAf/afQr6dhMmDfLDLcMCRMehm8SyQzbh59OGwl2pwNFSDXYdrxIdh4ic2LkLzdiSVQIA+M1sntVxBCw7ZBMCvVRYOSUSgHnujsnEsztEJMbL355Eu0nCbcMDMTHST3Qc6gUsO2QzHr11CLzcXFBY2YjtueWi4xCREzpR1YjPcsoAAE8mRQlOQ72FZYdshtpDiUduGQyg4ycro0lwIiJyNi99cwKSBMwZGYLY/mrRcaiXsOyQTVk5dRD8PF1xtrYZnxwpFR2HiJxIXqkGqccqIZOZ5+qQ42DZIZvST+WCn08bAsB8N4Su3Sg4ERE5i85J7ovGhGN4sJfgNNSbWHbI5qxIGIhgbxXKGlrxYcZ50XGIyAkcKq7D9ydqoJDLsGbmMNFxqJdZtezU19cjOTkZarUaarUaycnJaGho6PYYnU6Hxx9/HAEBAfD09MSCBQtQWnr55YxDhw7h9ttvh4+PD3x9fZGUlIScnBxrvhXqQ25KBVbPMP9j88re02jRtwtORESOTJIkvPC1+azO4gn9MdDfU3Ai6m1WLTvLly9HTk4OUlNTkZqaipycHCQnJ3d7zJo1a7Bt2zZs3rwZ+/fvR1NTE+bPnw+j0Xw5o7GxEbNnz8aAAQOQkZGB/fv3w9vbG7Nnz4bBYLDm26E+tGRCBPr7uqO2SYf3Dp4THYeIHNi+k7XIPFsHV4Ucj8/gWR1HJJOsNK62oKAAI0aMQHp6OuLj4wEA6enpSExMRGFhIaKirrylT6PRIDAwEBs2bMCSJUsAAOXl5YiIiMDOnTsxe/ZsZGVlYeLEiTh//jwiIiIAAHl5eYiLi8OpU6cwZMiQK15Xp9NBp9NZ/lur1SIiIgIajQbe3t7WePvUC7YeLsWTW45C7a7Evt9Oh7ebUnQkInIwkiRh0asHcLRUg5VTIvHsnSNFR6JuaLVaqNXqHn9+W+3MTlpaGtRqtaXoAEBCQgLUajUOHjzY5TGHDx+GwWBAUlKS5WthYWEYNWqU5ZioqCgEBATgnXfegV6vR2trK9555x2MHDkSAwcO7PJ1161bZ7mUplarLSWJbNtdY8MxNKgfNK0GvPbdadFxiMgB7cirwNFSDTxcFfj5tKGi45CVWK3sVFZWIigo6IqvBwUFobKy8qrHuLq6wtfX97KvBwcHW47x8vLCd999h40bN8Ld3R39+vXD119/jZ07d8LFxaXL1127di00Go3lUVJScpPvjvqCQi7DU3OiAQDv7D+L0voWwYmIyJHo2o34e2ohAPNQ00AvleBEZC09LjvPPfccZDJZt4+srCwA6HI/EUmSerzPyKXHtLa24qGHHsKUKVOQnp6OAwcOYOTIkZg7dy5aW1u7PF6lUsHb2/uyB9mH22OCkDDYD/p2E/7RsYCQiKg3vH/wHErqWhHkpULKrYNExyEr6vpUSDdWr16NpUuXdvucyMhI5Obmoqrqyg0da2pqEBwc3OVxISEh0Ov1qK+vv+zsTnV1NSZPngwA+OCDD1BcXIy0tDTI5XLL13x9ffH5559fMxvZF5lMhqfnjsCdr+zHZznleGjqIMT19xEdi4jsXH2zHv/ecxKAeVsID9cefxySHenxmZ2AgABER0d3+3Bzc0NiYiI0Gg0yMzMtx2ZkZECj0ViKy4+NHz8eSqUSu3btsnytoqIC+fn5lmNaWlogl8svOzvU+d8mE7cXcESx/dW4e2w4AOAvOwpgpTX1RORE/r3nFLRt7YgO8cI94/uLjkNWZrU1OzExMZgzZw5SUlKQnp6O9PR0pKSkYP78+ZY7scrKyhAdHW0pRGq1GqtWrcITTzyB3bt3Izs7GytWrEBsbCxmzpwJAJg1axbq6+vx2GOPoaCgAMeOHcPKlSvh4uKC6dOnW+vtkGBPzI6CykWOjLN1+LagWnQcIrJjxbXN2JBeDAB4el4MFPKeLa0g+2PVOTubNm1CbGwskpKSkJSUhLi4OGzYsMHyfYPBgKKiIrS0XFx4un79eixatAiLFy/GlClT4OHhge3bt0OhUAAAoqOjsX37duTm5iIxMRG33HILysvLkZqaitDQUGu+HRIo3Mcdq6aar6mv+6oABm4SSkQ36PmvC2EwSrhteCBuGRYoOg71AavN2bFlN3qfPonV2GbAtBe+w4VmPf5v4UgkJ0aKjkREdubwuTrc81oa5DLgq1/eiqgQ7oFlT2xuzg5Rb/NyU1r2rFn/7Ulo2zgxm4iunyRJ+POOAgDA4gkRLDpOhGWH7MrSSQMwJNATdc16Dhokoh7ZkVeB7PMN8HBV4NezhouOQ32IZYfsilIhx9o7YgCYBw2WNXQ9W4mI6FI/HiAY5O0mOBH1JZYdsjscNEhEPbUhjQMEnRnLDtmdzkGDALAtuwx5pRrBiYjIljW06PGv3Rwg6MxYdsguxfZX467OQYM7j3PQIBFdFQcIEssO2a0nOwYNpp/hoEEi6lpxbTPeTysGAPxuLgcIOiuWHbJbHDRIRNdy6QDBW4dzgKCzYtkhu/azaUPg7+mKMzXN+DDzvOg4RGRDsorrsDOvEnKZ+awOOS+WHbJrXm5KrOmYl/HiNydQ16wXnIiIbIHRJOHZL44B4ABBYtkhB7BsYgRiQr2haTXgha8LRcchIhvwQeZ5HCvXwtvNBb+ZHSU6DgnGskN2z0Uhx58WjgQAbD5UgqMlDYITEZFIdc16ywyuJ5Ki4N9PJTgRicayQw5hYqQf7h4bDkkC/vB5Pkwm3opO5KyeTy2EptWAmFBv3Bc/QHQcsgEsO+QwnpobjX4qFxwt1eCjrBLRcYhIgJySBsv///+3cCRcFPyYI5YdciBBXm6WXdGfTy1EQwsXKxM5E5NJwh8+z4ckAXePDceESD/RkchGsOyQQ3lgciSGB/dDfYsB//iG+2YROZOPskqQW6qBl8oFT82NFh2HbAjLDjkUpUKOPy4YBQDYlHEe+WXcN4vIGTS06PF8x67ma2YNR5AXdzWni1h2yOEkDvHHnaPDIEnA77lYmcgpvPB1EepbDBge3A/3Jw4UHYdsDMsOOaSn58bA01WB7PMN+ORIqeg4RGRF+WUafNAxQf1PC0dByUXJ9CP8G0EOKUTthl/cbl6s/LevzLehEpHjMZkk/L5jUfKC0WFIGOwvOhLZIJYdclgrpwzCkEBPXGjWY/2uE6LjEJEVbD1SiuzzDfB0VXD/K7oqlh1yWK4uFxcrv59WjIIKrdhARNSrNK0G/P0r86LkX9w+DCFqLkqmrrHskEObOiwAc2NDYOqYrCxJXKxM5CjW7zqBC816DAn0xMopg0THIRvGskMO7+l5I+CuVOBQcT0+yykTHYeIesHxci3eTysGAPxxwSi4uvDjjK6OfzvI4YX7uGP1jKEAgL/sKISmhYuViexZ56JkkwTMjQ3B1GEBoiORjWPZIafw8C2DMDjAE7VNOqz7qkB0HCK6CZsyzuHwuXp4uCrw9LwRouOQHWDZIaegclFg3d2xAIDNh0pw8HSt4EREdCPKG1rxt45Fyf9vdhTCfdwFJyJ7wLJDTiN+sD/uix8AAFj7aR5a9UbBiYioJyRJwjOf5aNZb8S4AT5ITowUHYnsBMsOOZWn7ohGiLcbzl1owcvfcvYOkT354mg59hRWw1Uhx9/viYNCLhMdiewEyw45FS83Jf68yDx75619Z5BXyo1CiexBXbMef9x+HADw2PShGBbsJTgR2ROWHXI6M0cEY35cKEwS8P8+yYXBaBIdiYiu4f++PI66Zj2igr3ws2lDRMchO8OyQ07puQUj4eOhREGFFm/tOyM6DhF147uiamzLLoNMBvztnljO1KEe498YckoB/VR4puOW1Ze/PYkzNU2CExFRV5p17Xh6Wz4AYOXkQRg7wFdwIrJHLDvktO4ZF45bhgVA327CU5/mwWTiVhJEtuaFr4tQ1tCK/r7ueHL2cNFxyE6x7JDTkslk+OtdsXBXKpB5tg4fHjovOhIRXeLwuXq8l1YMAPjrXbHwcHURmofsF8sOObUIPw88OTsKAPC3nYWo1LQJTkREAKBrN+KpT3IhScA94/rj1uGBoiORHWPZIaf34ORIjInwQaOuHc98xp3RiWzBf/aexsnqJgT0c8Xv58eIjkN2jmWHnJ5CLsPf74mDUiHDtwVV2JlXKToSkVM7UdWI/3x3CkDnnZOughORvWPZIQIQFeKFn00z74z+7Bf5aGjRC05E5JyMJgm//SQXBqOEmTHBmBcbKjoSOQCWHaIOj00fgqFB/VDbpOflLCJB3vzhDLLPN8BL5YI/LxoFmYxbQtDNY9kh6qByUeDFe0dDIZfhy9wKfJ5TLjoSkVPJL9PgpV1FAIDf3zkCIWo3wYnIUbDsEF1idIQPfjFjGADg95/no7S+RXAiIufQZjDiVx/lwGCUMHtkMO4d3190JHIgLDtEP/LY9CEYO8AHjW3teOLjozBy2CCR1f3tq0KcrG5CoJcK6+6O4+Ur6lUsO0Q/4qKQ4+UlY+DhqkDG2Tq8zb2ziKxq38kavHuwGADwwk/i4OfJu6+od7HsEHVhoL8n/jDfvHfWP74pwvFyreBERI6pvlmPJ7ccBQDcnzgQ06KCBCciR8SyQ3QVSyZGYGZMMAxGCWs+ykabwSg6EpFDkSQJT3+WhyqtDkMCPbH2Dg4PJOtg2SG6CplMhr/fE4uAfiqcqGrC86lFoiMROZRPj5RhZ14lXOQyvLxkLNxdFaIjkYOyatmpr69HcnIy1Go11Go1kpOT0dDQ0O0xb775JqZNmwZvb2/IZLIun38jr0t0I/z7qfD8T2IBAP89cBb7T9YKTkTkGErqWvDsF8cAAL+aNRyx/dWCE5Ejs2rZWb58OXJycpCamorU1FTk5OQgOTm522NaWlowZ84c/O53v+vV1yW6UTOig3Ff/AAAwBNbcjhdmegmGU0Sfv1xDpp07Zgw0Bc/vW2I6Ejk4GSSlcbEFhQUYMSIEUhPT0d8fDwAID09HYmJiSgsLERUVFS3x3/33XeYPn066uvr4ePj02uvCwBarRZqtRoajQbe3t438S7JWbTo2zHvX/txtrYZ8+JC8cqysbw1lugG/ee7U3g+tQiergqkrrkVEX4eoiORnbjRz2+rndlJS0uDWq22FBIASEhIgFqtxsGDB/v0dXU6HbRa7WUPop7wcHXB+iVjoJDLsCO3Ap/llImORGSX8ss0WL/rBADg2QUjWXSoT1it7FRWViIo6MpbCIOCglBZeeO7St/I665bt86yvketViMiIuKGf31yXmMifPDL283Tlf/w2TFOVybqoTaDEWs6piTPGRnCKcnUZ3pcdp577jnIZLJuH1lZWQDQ5Wl+SZJu+vR/T1937dq10Gg0lkdJSclN/frkvH4+rWO6sq4dazbnwGA0iY5EZDf+vOM4TnVMSf7r3bG8FEx9xqWnB6xevRpLly7t9jmRkZHIzc1FVVXVFd+rqalBcHBwT39Zi5CQkB6/rkqlgkqluuFfk6hT53Tlef/aj6xz9Xjh6yL8bi5ngxBdy+c5ZdiYfh4A8I97R3NKMvWpHpedgIAABAQEXPN5iYmJ0Gg0yMzMxKRJkwAAGRkZ0Gg0mDx5cs+TWvl1ia7XQH9PPP+TOPx80xG8+cMZjBvgizmjQkTHIrJZJ6sa8dQneQCAx2cMxW3DAwUnImdjtTU7MTExmDNnDlJSUpCeno709HSkpKRg/vz5ljumysrKEB0djczMTMtxlZWVyMnJwalTpwAAeXl5yMnJQV1d3XW/LpG1zY0NxaqpgwAAv9lyFOcuNAtORGSbmnXt+OnGw2g1GDFlqD/WzBwuOhI5IavO2dm0aRNiY2ORlJSEpKQkxMXFYcOGDZbvGwwGFBUVoaXl4kLP119/HWPHjkVKSgoA4NZbb8XYsWPxxRdfXPfrEvWFp+6IxviBvmjUteNnG49wOwmiH5EkCU99mofTNc0I9lbhn0vHQiHnOh3qe1abs2PLOGeHekulpg3z/rUPF5r1WDoxAn+7J050JCKbsSGtGL///Bhc5DJsfiQBEyL9REciO2dzc3aInEGI2g3/XDoWMhmw+VAJtmTxTj8iAMgpacCfvjwOwHwWlEWHRGLZIbpJU4cF4Fcd6xCe+Swfx8s5tJKcW32zHo9tOmKZp9O5vo1IFJYdol6wevpQTIsKhK7dhJ9vOgxtm0F0JCIhTCYJv/o4B2UNrYj098Dz98Zxng4Jx7JD1AvkchnWLx6DcB93FF9owW+35sIJl8MR4T/fncJ3RTVQucjx2orx8HZTio5ExLJD1Ft8PV3x6n3joFTI8FV+Jd7Zf1Z0JKI+deBULV7q2Pfqz4tGISaUN4CQbWDZIepFYyJ88Pv5IwAAf/uqEFnFdYITEfWNSk0bfvFhNkwSsGRCBO6dwD0IyXaw7BD1suSEgVgwOgztJgmPfXAE1do20ZGIrErXbsRjHxzBhWY9RoR6448LR4qORHQZlh2iXiaTybDu7lgMDeqHKq0OKe9noVXPgYPkmCRJwtpP8nD4XD283Fzw2opxcFMqRMciugzLDpEVeKpc8M4DE+DrocTRUg2e2JIDk4kLlsnxvLr3FD7NLoNCLsNr943HQH9P0ZGIrsCyQ2QlA/098UbyBCgVMuzMq8SLu4pERyLqVV/mluMf35gXJP9p4UhMHXbtTaKJRGDZIbKiSYP88Le7zVtIvLr3NLYeLhWciKh3ZJ+vxxMfHwUArJo6CPfFDxSciOjqWHaIrOye8f2xevpQAMDaT3ORceaC4EREN6e0vgUp7x+Grt2EmTFB+N3cGNGRiLrFskPUB349azjmxYbCYJTw6MbDKK5tFh2J6IY0thnw8HtZqG3SISbUmzuZk11g2SHqA3K5DC8uHo3RET5oaDHgoXcPQdPCLSXIvrQbTfjFh9korGxEkJcK7zwwAZ4qF9GxiK6JZYeoj7gpFXjr/vEIU7vhTG0zfrbpMAxGk+hYRNftzzsKsLeoBm5KOd5+YALCfNxFRyK6Liw7RH0oyMsN7zw4EZ6uChw8fQG//yyfe2iRXdiQVox3DxYDAF5aPAZx/X2E5iHqCZYdoj4WE+qNfy8fC7kM2HyoBG/tOyM6ElG3vj9Rg+e2HwcA/GZ2FObGhgpORNQzLDtEAsyIDsYz88x7aK37qhBf5VUITkTUtYIKLVZvOgKjScI94/rj59OGiI5E1GMsO0SCrJwSieSEgZAk4Jebc/DDiRrRkYguc6amCcnvZKBR145Jg/yw7u5YyGS884rsD8sOkSAymQzP3jkCc2NDoDea8MiGLBziLulkI0rrW7Di7QzUNpk393zr/glwdeFHBtkn/s0lEshFIcfLS8ZiWlQg2gwmPPS/Q8gr1YiORU6uurENK97OQLmmDUMCPbFh1SSo3ZWiYxHdMJYdIsFcXeR4fcV4xA/yQ6OuHff/NwMnqhpFxyInVd+sR/LbmSi+0IL+vu7Y+HA8/PupRMciuiksO0Q2wE2pwDsPTsToCB/Utxiw4u0MnLvAKcvUtxrbDHjwf5koqjIPDfzg4QSEqjlLh+wfyw6RjeincsF7KyciOsQL1Y06LH8rAxWaVtGxyEm06o1Y9V4WjpZq4Ofpik0Px2OAv4foWES9gmWHyIb4eLji/VWTMCjAE2UNrbjv7QzUNulExyIHp2834acbDyPzbB28VC54/6FJGBbsJToWUa9h2SGyMUFebtj4cDzCfdxxpqYZye9kch8tspp2owm/3JyN70/UwF2pwP9WTsSocLXoWES9imWHyAaF+5gXhgb0U6GgQosH381Es65ddCxyMCaThN9+koev8ivhqpDjzfvHY0Kkn+hYRL2OZYfIRg0K8MTGh823/Gafb8Cq9w6hiYWHeonJJOGZz/PxyZFSKOQy/Hv5WNwyLFB0LCKrYNkhsmHRId5476FJ6KdyQfqZOtz3dgbqm/WiY5GdMxhN+NXHOfgg4zxkMuDFe0dj9sgQ0bGIrIZlh8jGjYnwwaaH4+HjocTRkgYseTMNVdo20bHITrXqjXh0w2F8nlMOF7kM/1w6FovGhouORWRVLDtEdmB0hA+2PJqIYG8VTlQ14Z7XDnIOD/WYts2AB/6biT2F1VC5yPHW/ROwYHSY6FhEVseyQ2QnhgV7YetPJ2OgvwdK61vxk9fTUFChFR2L7ERtkw7L3kxHZrH59vINq+IxPTpIdCyiPsGyQ2RHIvw8sOWniYgO8UJNow5L3kjD4XP1omORjStraMXi19NwrFwLf09XfPhIAiYN4l1X5DxYdojsTJCXGz56JBHjBvhA29aOFW9n4IcTNaJjkY06XdOEe187iDO1zQhTu2HLTxM5R4ecDssOkR1Seyix8eF43DIsAK0GI1a9dwg78ypExyIbk1+mwb2vp6Fc04bBgZ7Y+rPJGBzYT3Qsoj7HskNkpzxcXfD2AxMwLzYUBqOE1R8cwUeHzouORTYi48wFLHszHXXNeowK98aWRxMR5sNNPck5sewQ2TGViwL/WjYWSydGwCQBv/0kDy9+UwSTSRIdjQT6LLsM9/83E426dkwa5IcPUxLg308lOhaRMCw7RHZOIZdh3d2x+OltQwAA/95zCo9syEJjG/fTcjbtRhP+/OVxrPkoB7p2E2bGBOH9hybBy00pOhqRUCw7RA5AJpPhqTui8eK9o+HqIse3BdVY9OoBnK5pEh2N+khDix4r3z2Et/efBQA8Nn0I3kieADelQnAyIvFYdogcyD3j+2PLo4kI8XbD6ZpmLHrlAPYWVouORVZWWKnFglcOYN/JWrgrFXh1+Tj8ZnY0FHKZ6GhENoFlh8jBjI7wwfbHp2LCQF806trx0HuH8OreU5AkruNxRF/lVeDu/xzE+boWRPi549OfT8a8uFDRsYhsCssOkQMK9LUiDpYAABKDSURBVFLhg5QELI8fAEkCXvi6CKs/yEaLnrumOwqTScKL3xThZ5uOoEVvxJSh/vjisamICfUWHY3I5rDsEDkoVxc5/npXLP5y1ygoFTLs6DgDUFLXIjoa3aTGNgMe2ZCFf+85BQB4eOogvLdyEnw9XQUnI7JNLDtEDu6++IH4MCUBAf1UKKxsxJ2v7OfEZTt2sqoRi149gG8LquHqIsdLi0fjmfkj4KLgP+dEV8P/O4icwIRIP2x/fApG91ejocWA+/+biT98ns/LWnbEaJLw1g9nMO/f+3G6phmhajds/Wki7h7XX3Q0Ipsnk5xw1aJWq4VarYZGo4G3N69vk/NoMxjx150FeD/tHAAg0t8DLy4ejfEDuSmkLTt/oQVPbjmKzOI6AMC0qEC88JPRCPTioEByLjf6+c2yw7JDTmj/yVr8ZutRVGjaIJcBj9w6BL+aNQwqF85ksSWSJOHDzBL8ecdxtOiN8HRV4PfzR2DJxAjIZLytnJzPjX5+W/UyVn19PZKTk6FWq6FWq5GcnIyGhoZuj3nzzTcxbdo0eHt7QyaTXfH84uJirFq1CoMGDYK7uzuGDBmCZ599Fnq93ppvhcihTB0WgNQ1t+LuceEwScDr35/GwlcO4Fi5RnQ06lClbcPKdw/hd9vy0KI3YtIgP6SuuRVLJw1g0SHqIauWneXLlyMnJwepqalITU1FTk4OkpOTuz2mpaUFc+bMwe9+97suv19YWAiTyYQ33ngDx44dw/r16/H6669f9flE1DW1uxIvLR6DN5LHw9/TFYWV5oWvr+w5iXajSXQ8pyVJEj7PKUPS+h/wXVENXF3keGZeDDanJCDCz0N0PCK7ZLXLWAUFBRgxYgTS0/9/e/cfFVXd5wH8Pb9gEIdRRIYRRn49xY/MVjERwRQ1tHq02j25nopjlpmFJ7Wz26HjmvRDsuzxaU+pJbmmW6tn07XHengSehI10Fw8UG6DpuAkNg6KIjOCAgPf/QMZJdC4IzMDl/frnDmOl++d+fBh4L7Pvd977yEkJycDAA4dOoSUlBQcO3YMcXFxt1y/qKgI6enpqKurw5AhQ245ds2aNdiwYQOqqqp6VBsPYxF1Vnu5Cct3HcWen2oAtF+YcO2cexA7fLCPKxtYLjY049++OIr8ozYAwOgIPdbOuQd/CNX5uDKivqHPHcY6ePAg9Hq9K+gAwIQJE6DX61FSUtKr71VfX4/g4JtPsGxqaoLdbu/0IKLrQgb748Mnk7B2zj3QadX4ofoSHvj3A3jrbxWw84aiHtfS2oZPik9h2p+KkH/UBrVSgWXT78TO5ycy6BD1Ao+FHZvNhtDQ0C7LQ0NDYbPZeu19Kisr8f7772PRokU3HfPWW2+55g3p9XqYTKZee38iuVAoFPjHsREoWHYf7rtzOJqdbfhoXxWmrCnC1oMWtPDQVq8TQqDgJxtm/Hk/cr40o66xBXEGHXa9kIol0++AhtfOIeoVkn+TcnJyoFAobvkoLS0FgG4n0Qkhem1yndVqxcyZM/HYY49hwYIFNx33yiuvoL6+3vWorq7ulfcnkiOjPgBb5t+LTfPGIXZ4IC42NOPVv/yEme/tx98raniPrV5y9Ew95m48hIX/eQRVtQ0IGeyHVY+Owl9fTMPdEXpfl0ckK2qpKyxevBhz58695ZioqCj8+OOPqKmp6fK18+fPw2AwSH3bLqxWK9LT05GSkoKNGzfecqy/vz/8/Xk9CqKeUigUmJZgwH13Dsf2w6fx529OoPJ8A57ZUoqJscOw/KEE3DWCG2R3WC9dwbt7juN/yn4FAPirlXgmLRrPT4mFTqvxcXVE8iQ57ISEhCAkJOR3x6WkpKC+vh6HDx/G+PHjAQDff/896uvrMXHiROmV3uDXX39Feno6kpKSsHnzZiiV3NVL5AkalRKZKVF4eEw41u+txH8Un0JJ5QX88f3v8E9jI/AvGXEI02t9XWa/cLnJiQ+LKpF3oApNzvZDgo/8wwj868x4hA8J8HF1RPLm0YsKPvDAA7Barfjoo48AAAsXLkRkZCS+/PJLAO2hZdq0adi6dasrENlsNthsNpSWluLZZ5/F/v37odPpMHLkSAQHB8NqtWLy5MkYOXIktm7dCpXq+kXQwsLCelQXz8Yick/1xUa8s+c4vvzBCgAI0Kjwz/ea8NTEKESFBPq4ur7pwuUm/Nf3p7HloAW1l9uvBzY+KhjLH0rAPaZbn2lKRJ31ySsoX7x4ES+++CJ2794NAJg9ezY++OAD16nkFosF0dHR2Lt3L6ZMmQKgfU7Qa6+91uW1Nm/ejKeeegqffPIJ5s+f3+379fRbYdghuj1lp+vw5l8rcOSXOgCAQgFMjQvF02nRmBg7jBe9A2C22rG5+BT+8oMVzdf25EQNG4TsBxIw4y4De0Tkhj4Zdvoqhh2i2yeEwIETtdhcfAp7j1+/i3qcQYenUqPw6JhwaDUD6/YTrW0C31TUYHPxKRyquuhaPjpCj6dTo/Hg3Ub4qXnYnchdDDsSMOwQ9a7K85expcSCHUfOoLG5FQAwZJAGj48ficyUSBj18p6TYr/agv/+32psOWhB9cUrAACVUoGZo8LwdGoUxo4cyj05RL2AYUcChh0iz6i/0oLPS6vxSYkFZ+qub/Snxodixl1hmBofiuBAPx9X2Tsam53Y/3MtCs01+Pr/zqLhWsjTB2jwePJIZE6IxAhOPCbqVQw7EjDsEHlWa5tAobn9cM73p64fzlEqgHujgnF/ogEZiWEYOax/3evpwuUm/L3iHArMNhw4Ues6qwoA7ggdjPmp0Xh0TDgC/AbW4Tsib2HYkYBhh8h7jtns+NtRGwrMNag42/lWLfFhOmQkGnB/YhhGhQf1yUM9ltoGFJprUGC24cgvdWi74S9mxNAAZCSGYcZdBoyPDu6T9RPJCcOOBAw7RL5RfbERheYaFJprcNhyEa03JIdQnT9GR+iRaAxC4oggJBiDYBo6CEqldwKEEALnHE0wn7XDbLXDfNaOn36th+VCY6dxo8KDcH9CGDLuMiA+TMeAQ+RFDDsSMOwQ+d6lxmZ8e+wcCn6qwf4T510Tm2802F+NBKMOCcYgJBrbA1CYXosgrQZajdKtoNHkbIX9ihMXGppw7KwD5rN2VFwLOBcamruMVysVSI4JRkZiGKYnGngBQCIfYtiRgGGHqG+52tKKH6ovoeKsHRXXAsjxGofr+jTd0agUCNJqEBSggU6rvva8/V+lUgH7lRbYrzphv9ICx9Xrz5tu8ZpKBRA7fLBrz1KiMQj3RAyBfhBv40DUF7i7/ZZ8uwgiot6m1aiQHDMMyTHDXMucrW2oqm1wHVLqCEJ1jc1obRNoaRW40NDc7d6Y36NQAEFaDe4IbQ82HXuN4sJ0A+7aQEQDAcMOEfVJapUSdxp0uNOgwyNjwl3LhRBobG6F/WoL7Fec1/5tgePq9efONgF9gAZB2mt7fa49Dwpofz7YT+21uUBE5HsMO0TUrygUCgT6qxHor4aRN14noh7gdcuJiIhI1hh2iIiISNYYdoiIiEjWGHaIiIhI1hh2iIiISNYYdoiIiEjWGHaIiIhI1hh2iIiISNYYdoiIiEjWGHaIiIhI1hh2iIiISNYYdoiIiEjWGHaIiIhI1gbkXc+FEAAAu93u40qIiIiopzq22x3b8Z4akGHH4XAAAEwmk48rISIiIqkcDgf0en2PxyuE1HgkA21tbbBardDpdFAoFD1ax263w2Qyobq6GkFBQR6ukDqw777BvvsG++4b7LtvuNN3IQQcDgdGjBgBpbLnM3EG5J4dpVKJiIgIt9YNCgriL4MPsO++wb77BvvuG+y7b0jtu5Q9Oh04QZmIiIhkjWGHiIiIZE2Vk5OT4+si+guVSoUpU6ZArR6QR/98hn33DfbdN9h332DffcNbfR+QE5SJiIho4OBhLCIiIpI1hh0iIiKSNYYdIiIikjWGHSIiIpI1hh0iIiKSNYadG6xfvx7R0dHQarVISkrCgQMHbjl+586dSExMhL+/PxITE7Fr1y4vVSovUvqel5eHSZMmYejQoRg6dCimT5+Ow4cPe7Fa+ZD6ee+wfft2KBQKPPLIIx6uUJ6k9v3SpUvIysqC0WiEVqtFQkIC8vPzvVStfEjt+3vvvYe4uDgEBATAZDJh2bJluHr1qpeqlYf9+/dj1qxZGDFiBBQKBb744ovfXWffvn1ISkqCVqtFTEwMPvzww94pRpAQQojt27cLjUYj8vLyhNlsFkuWLBGBgYHil19+6XZ8SUmJUKlUIjc3V1RUVIjc3FyhVqvFoUOHvFx5/ya1748//rhYt26dKCsrExUVFWL+/PlCr9eLM2fOeLny/k1q3ztYLBYRHh4uJk2aJB5++GEvVSsfUvve1NQkxo0bJx588EHx3XffCYvFIg4cOCDKy8u9XHn/JrXvn376qfD39xefffaZOHXqlNizZ48wGo1i6dKlXq68f8vPzxfLly8XO3fuFADErl27bjm+qqpKDBo0SCxZskSYzWaRl5cnNBqN2LFjx23XwrBzzfjx48WiRYs6LYuPjxfZ2dndjp8zZ46YOXNmp2UzZswQc+fO9ViNciS177/ldDqFTqcTW7Zs8UR5suVO351Op0hNTRUff/yxmDdvHsOOG6T2fcOGDSImJkY0Nzd7ozzZktr3rKwsMXXq1E7LXnrpJZGWluaxGuWuJ2Hn5ZdfFvHx8Z2WPffcc2LChAm3/f48jAWgubkZR44cQUZGRqflGRkZKCkp6XadgwcPdhk/Y8aMm46nrtzp+281NjaipaUFwcHBnihRltzt++uvv47hw4fjmWee8XSJsuRO33fv3o2UlBRkZWXBYDBg1KhRyM3NRWtrqzdKlgV3+p6WloYjR464DpFXVVUhPz8fDz30kMfrHchutl0tLS1FS0vLbb02r4sNoLa2Fq2trTAYDJ2WGwwG2Gy2btex2WySxlNX7vT9t7KzsxEeHo7p06d7okRZcqfvxcXF2LRpE8rLy71Roiy50/eqqip8++23eOKJJ5Cfn48TJ04gKysLTqcTr776qjfK7vfc6fvcuXNx/vx5pKWlQQgBp9OJ559/HtnZ2d4oecC62XbV6XSitrYWRqPR7ddm2LmBQqHo9H8hRJdltzOeuuduH9955x1s27YNRUVF0Gq1nipPtnrad4fDgSeffBJ5eXkICQnxVnmyJeXz3tbWhtDQUGzcuBEqlQpJSUmwWq1Ys2YNw45EUvpeVFSEVatWYf369UhOTsbJkyexZMkSGI1GrFixwhvlDljd/Zy6Wy4Vww6AkJAQqFSqLin/3LlzXVJmh7CwMEnjqSt3+t7h3XffRW5uLr755huMHj3ak2XKjtS+V1ZWwmKxYNasWa5lbW1tAAC1Wo3jx48jNjbWs0XLgDufd6PRCI1GA5VK5VqWkJAAm82G5uZm+Pn5ebRmOXCn7ytWrEBmZiYWLFgAALj77rvR0NCAhQsXYvny5VAqOQPEE262XVWr1Rg2bNhtvTZ/YgD8/PyQlJSEwsLCTssLCwsxceLEbtdJSUnpMr6goOCm46krd/oOAGvWrMEbb7yBr7/+GuPGjfN0mbIjte/x8fE4evQoysvLXY/Zs2cjPT0d5eXlMJlM3iq9X3Pn856amoqTJ0+6wiUA/PzzzzAajQw6PeRO3xsbG7sEGpVKBdF+Uo/Hah3obrZdHTduHDQaze29+G1PcZaJjlMTN23aJMxms1i6dKkIDAwUFotFCCFEZmZmp5n7xcXFQqVSidWrV4uKigqxevVqnnruBql9f/vtt4Wfn5/YsWOHOHv2rOvhcDh89S30S1L7/ls8G8s9Uvt++vRpMXjwYLF48WJx/Phx8dVXX4nQ0FDx5ptv+upb6Jek9n3lypVCp9OJbdu2iaqqKlFQUCBiY2PFnDlzfPUt9EsOh0OUlZWJsrIyAUCsXbtWlJWVuU75z87OFpmZma7xHaeeL1u2TJjNZrFp0yaeeu4J69atE5GRkcLPz0+MHTtW7Nu3z/W1yZMni3nz5nUa//nnn4u4uDih0WhEfHy82Llzp5crlgcpfY+MjBQAujxWrlzp/cL7Oamf9xsx7LhPat9LSkpEcnKy8Pf3FzExMWLVqlXC6XR6uer+T0rfW1paRE5OjoiNjRVarVaYTCbxwgsviLq6Oh9U3n/t3bu327/XHb2eN2+emDx5cqd1ioqKxJgxY4Sfn5+IiooSGzZs6JVaFEJwnxwRERHJF+fsEBERkawx7BAREZGsMewQERGRrDHsEBERkawx7BAREZGsMewQERGRrDHsEBERkawx7BAREZGsMewQERGRrDHsEBERkawx7BAREZGs/T+maxSuA3KUHAAAAABJRU5ErkJggg==",
      "text/plain": [
       "PyPlot.Figure(PyObject <matplotlib.figure.Figure object at 0x128596fd0>)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "1-element Array{PyCall.PyObject,1}:\n",
       " PyObject <matplotlib.lines.Line2D object at 0x12d5c83d0>"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "using PyPlot\n",
    "\n",
    "plot(xgrid[2:end-1], U)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Poisson equation in three dimensions\n",
    "\n",
    "$$-(u_{xx}+u_{yy}+u_{zz})=f(x,y,z),~~~0<x,y,z<1$$\n",
    "\n",
    "with zero Dirichlet boundary conditions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaQAAAGgCAYAAADl3RMjAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAHOtJREFUeJzt3W9slfX9//HXobWHP2mPaw2tZ55CTYhoQVSqRkCFqF2Q4JhRh2hlekeyKtQaBohuzIWe4TKGWwem3lAXg3JHKnMz2ikWCP4pLVXnFpHZQJU1nYs5p8A8QHv9bnxDf9aWtqdc17neV3k+khNyrnP1ul65zp8X15/zOSHHcRwBAOCzMX4HAABAopAAAEZQSAAAEygkAIAJFBIAwAQKCQBgAoUEADCBQgIAmEAhAQBMoJAAACYEupA2b96skpISjR07VjNnztTu3bv9jjSkeDyuq6++Wrm5uZo4caIWLVqkTz/91O9YaYvH4wqFQqqqqvI7ypC+/PJL3XvvvSooKND48eN1xRVXqLm52e9Ygzp16pQef/xxlZSUaNy4cbr44ov15JNPqqenx+9o/ezatUsLFy5UNBpVKBRSfX19n8cdx9G6desUjUY1btw4zZ07V5988olPaf+/wXKfPHlSq1at0vTp0zVhwgRFo1Hdd999OnLkiI+J/89Q2/vbHnzwQYVCIW3atCmDCUcusIW0bds2VVVVae3atdq/f7+uv/56zZ8/X4cPH/Y72qAaGxtVWVmp9957Tw0NDTp16pTKy8t17Ngxv6MNW1NTk+rq6nT55Zf7HWVIX3/9tWbPnq3zzjtPr7/+uv7xj3/ot7/9rc4//3y/ow1qw4YNeuaZZ1RbW6t//vOfeuqpp/Sb3/xGf/jDH/yO1s+xY8c0Y8YM1dbWDvj4U089pY0bN6q2tlZNTU0qKirSLbfcoq6urgwn7Wuw3MePH1dLS4ueeOIJtbS06JVXXtGBAwd02223+ZC0r6G292n19fV6//33FY1GM5TMBU5AXXPNNc6yZcv6TJs6daqzevVqnxKNTGdnpyPJaWxs9DvKsHR1dTlTpkxxGhoanBtvvNFZsWKF35EGtWrVKmfOnDl+x0jbggULnAceeKDPtNtvv9259957fUo0PJKc7du3997v6elxioqKnF//+te907755hsnEok4zzzzjB8RB/Td3AP54IMPHEnOoUOHMpRqaGfK/cUXXzjf//73nb///e/OpEmTnN/97nc+pEtfIPeQTpw4oebmZpWXl/eZXl5err179/qUamQSiYQkKT8/3+ckw1NZWakFCxbo5ptv9jvKsOzYsUNlZWW68847NXHiRF155ZV69tln/Y41pDlz5uitt97SgQMHJEkffvih9uzZo1tvvdXnZOlpa2tTR0dHn/dqOBzWjTfeGMj3aigUMr933dPTo4qKCq1cuVKlpaV+x0lLtt8BRuKrr75Sd3e3CgsL+0wvLCxUR0eHT6nS5ziOqqurNWfOHE2bNs3vOEN6+eWX1dLSoqamJr+jDNvnn3+uLVu2qLq6Wo899pg++OADLV++XOFwWPfdd5/f8c5o1apVSiQSmjp1qrKystTd3a3169fr7rvv9jtaWk6/Hwd6rx46dMiPSCPyzTffaPXq1VqyZIny8vL8jjOoDRs2KDs7W8uXL/c7StoCWUinhUKhPvcdx+k3zbKHHnpIH330kfbs2eN3lCG1t7drxYoVevPNNzV27Fi/4wxbT0+PysrKVFNTI0m68sor9cknn2jLli2mC2nbtm168cUXtXXrVpWWlqq1tVVVVVWKRqNaunSp3/HSFuT36smTJ7V48WL19PRo8+bNfscZVHNzs55++mm1tLQEZvt+WyAP2V1wwQXKysrqtzfU2dnZ739iVj388MPasWOHdu7cqYsuusjvOENqbm5WZ2enZs6cqezsbGVnZ6uxsVG///3vlZ2dre7ubr8jDujCCy/UZZdd1mfapZdeav7il5UrV2r16tVavHixpk+froqKCj3yyCOKx+N+R0tLUVGRJAX2vXry5EndddddamtrU0NDg/m9o927d6uzs1PFxcW979NDhw7p0Ucf1eTJk/2ON6RAFlJOTo5mzpyphoaGPtMbGho0a9Ysn1INj+M4euihh/TKK6/o7bffVklJid+RhuWmm27Sxx9/rNbW1t5bWVmZ7rnnHrW2tiorK8vviAOaPXt2v8vqDxw4oEmTJvmUaHiOHz+uMWP6vj2zsrJMXvY9mJKSEhUVFfV5r544cUKNjY3m36uny+izzz7T3/72NxUUFPgdaUgVFRX66KOP+rxPo9GoVq5cqTfeeMPveEMK7CG76upqVVRUqKysTNddd53q6up0+PBhLVu2zO9og6qsrNTWrVv16quvKjc3t/d/jpFIROPGjfM53Znl5ub2O881YcIEFRQUmD7/9cgjj2jWrFmqqanRXXfdpQ8++EB1dXWqq6vzO9qgFi5cqPXr16u4uFilpaXav3+/Nm7cqAceeMDvaP0cPXpUBw8e7L3f1tam1tZW5efnq7i4WFVVVaqpqdGUKVM0ZcoU1dTUaPz48VqyZImPqQfPHY1Gdccdd6ilpUWvvfaauru7e9+r+fn5ysnJ8Sv2kNv7u8V53nnnqaioSJdcckmmo6bP34v8zs4f//hHZ9KkSU5OTo5z1VVXBeLSaUkD3p577jm/o6UtCJd9O47j/PnPf3amTZvmhMNhZ+rUqU5dXZ3fkYaUTCadFStWOMXFxc7YsWOdiy++2Fm7dq2TSqX8jtbPzp07B3xNL1261HGc/7v0+xe/+IVTVFTkhMNh54YbbnA+/vhjf0M7g+dua2s743t1586dZnMPJEiXfYccx3Ey1H0AAJxRIM8hAQBGHwoJAGAChQQAMIFCAgCYQCEBAEygkAAAJgS+kFKplNatW6dUKuV3lLSQO7OCmlsKbnZyZ1ZQc39b4L+HlEwmFYlElEgkzI8z9W3kzqyg5paCm53cmRXU3N8W+D0kAMDoQCEBAEwwN7hqT0+Pjhw5otzc3GH9nkcymezzb1CQO7OCmlsKbnZyZ5al3I7jqKurS9FotN+o9YMxdw7piy++UCwW8zsGAOAstbe3p/V7b+b2kHJzc/2OMGKJRMLvCCMSiUQ8XX5QtwuAkUkmk4rFYml/npsrpCD+7O5pQb2yxWtsF+DclO7nORc1AABMoJAAACZQSAAAEygkAIAJFBIAwATPCmnz5s0qKSnR2LFjNXPmTO3evdurVQEARgFPCmnbtm2qqqrS2rVrtX//fl1//fWaP3++Dh8+7MXqAACjgCcjNVx77bW66qqrtGXLlt5pl156qRYtWqR4PN5n3lQq1We49NNfqAoiY4NeDJvX3/0K6nYBMDIjHXnc9T2kEydOqLm5WeXl5X2ml5eXa+/evf3mj8fjikQivbeglhEA4Oy4XkhfffWVuru7VVhY2Gd6YWGhOjo6+s2/Zs0aJRKJ3lt7e7vbkQAAAeDZ0EHfPQzkOM6Ah4bC4bDC4bBXMQAAAeH6HtIFF1ygrKysfntDnZ2d/faaAAA4zfVCysnJ0cyZM9XQ0NBnekNDg2bNmuX26gAAo4Qnh+yqq6tVUVGhsrIyXXfddaqrq9Phw4e1bNkyL1YHABgFPCmkH//4x/rvf/+rJ598Uv/+9781bdo0/fWvf9WkSZO8WB0AYBQw94uxp69fDyJjm3LY+B4SADeZ+R4SAAAjQSEBAEygkAAAJnj2xdhzkZfnYrw8D+P1OZ6gbhcAmcUeEgDABAoJAGAChQQAMIFCAgCYQCEBAEygkAAAJlBIAAATKCQAgAkUEgDABAoJAGAChQQAMIFCAgCYQCEBAEygkAAAJlBIAAATKCQAgAkUEgDABAoJAGAChQQAMIFCAgCYQCEBAEygkAAAJmT7HeBMEomE8vLyXF9uKBRyfZmZ4GVux3E8W7bXyw/ydgHQF3tIAAATKCQAgAkUEgDABAoJAGAChQQAMIFCAgCYQCEBAEygkAAAJrheSPF4XFdffbVyc3M1ceJELVq0SJ9++qnbqwEAjDKuF1JjY6MqKyv13nvvqaGhQadOnVJ5ebmOHTvm9qoAAKNIyPF4fJT//Oc/mjhxohobG3XDDTcMOX8ymVQkEmHooAwK8hA5DB0E2DPSz3HPx7JLJBKSpPz8/AEfT6VSSqVSvfeTyaTXkQAABnl6UYPjOKqurtacOXM0bdq0AeeJx+OKRCK9t1gs5mUkAIBRnh6yq6ys1F/+8hft2bNHF1100YDzDLSHFIvFOGSXQUE+NMUhO8Aec4fsHn74Ye3YsUO7du06YxlJUjgcVjgc9ioGACAgXC8kx3H08MMPa/v27XrnnXdUUlLi9ioAAKOQ64VUWVmprVu36tVXX1Vubq46OjokSZFIROPGjXN7dQCAUcL1c0hnOqb/3HPP6Sc/+cmQf89l35kX5HMlnEMC7DFzDok3MQBgJBjLDgBgAoUEADCBQgIAmOD50EHWeHmOiwsmMi+ozyfnWoH+2EMCAJhAIQEATKCQAAAmUEgAABMoJACACRQSAMAECgkAYAKFBAAwgUICAJhAIQEATKCQAAAmUEgAABMoJACACRQSAMAECgkAYAKFBAAwgUICAJhAIQEATKCQAAAmUEgAABMoJACACRQSAMCEbL8DnEkkEvFkuY7jeLJcr5ftpVAo5HeEEQvq8+nlNg/q6xBgDwkAYAKFBAAwgUICAJhAIQEATKCQAAAmUEgAABMoJACACRQSAMAEzwspHo8rFAqpqqrK61UBAALM00JqampSXV2dLr/8ci9XAwAYBTwrpKNHj+qee+7Rs88+q+9973terQYAMEp4VkiVlZVasGCBbr755kHnS6VSSiaTfW4AgHOPJ4Orvvzyy2ppaVFTU9OQ88bjcf3yl7/0IgYAIEBc30Nqb2/XihUr9OKLL2rs2LFDzr9mzRolEoneW3t7u9uRAAABEHJcHqu+vr5eP/rRj5SVldU7rbu7W6FQSGPGjFEqlerz2Hclk0nPfnpCYmj+gfDzE5nHz09gNDv9OZ5IJJSXlzfsv3P9kN1NN92kjz/+uM+0+++/X1OnTtWqVasGLSMAwLnL9ULKzc3VtGnT+kybMGGCCgoK+k0HAOA0RmoAAJiQkZ8wf+eddzKxGgBAgLGHBAAwgUICAJhAIQEATKCQAAAmZOSihpFI9wtVw8UXEkeXoD6fXi47qNsEYA8JAGAChQQAMIFCAgCYQCEBAEygkAAAJlBIAAATKCQAgAkUEgDABAoJAGAChQQAMIFCAgCYQCEBAEygkAAAJlBIAAATKCQAgAkUEgDABAoJAGAChQQAMIFCAgCYQCEBAEygkAAAJlBIAAATsv0OkGmO43i27FAo5Nmyvczt5bIlb7eLl3g++wvqNkEwsIcEADCBQgIAmEAhAQBMoJAAACZQSAAAEygkAIAJFBIAwARPCunLL7/Uvffeq4KCAo0fP15XXHGFmpubvVgVAGCUcP2LsV9//bVmz56tefPm6fXXX9fEiRP1r3/9S+eff77bqwIAjCKuF9KGDRsUi8X03HPP9U6bPHmy26sBAIwyrh+y27Fjh8rKynTnnXdq4sSJuvLKK/Xss8+ecf5UKqVkMtnnBgA497heSJ9//rm2bNmiKVOm6I033tCyZcu0fPly/elPfxpw/ng8rkgk0nuLxWJuRwIABEDIcXlEw5ycHJWVlWnv3r2905YvX66mpia9++67/eZPpVJKpVK995PJpGKxmBKJhPLy8tyM5jkGnhxYUAdX9VJQn09e4xiOZDKpSCSS9ue463tIF154oS677LI+0y699FIdPnx4wPnD4bDy8vL63AAA5x7XC2n27Nn69NNP+0w7cOCAJk2a5PaqAACjiOuF9Mgjj+i9995TTU2NDh48qK1bt6qurk6VlZVurwoAMIq4XkhXX321tm/frpdeeknTpk3Tr371K23atEn33HOP26sCAIwirl/UcLZGejLMAk74DoyLGvoL6vPJaxzDYeaiBgAARoJCAgCYQCEBAEygkAAAJrg+uOq5zMuTskE+mRzU7YL+gvpccsFEMLCHBAAwgUICAJhAIQEATKCQAAAmUEgAABMoJACACRQSAMAECgkAYAKFBAAwgUICAJhAIQEATKCQAAAmUEgAABMoJACACRQSAMAECgkAYAKFBAAwgUICAJhAIQEATKCQAAAmUEgAABMoJACACdl+B8DwOI7j2bJDoZBny5a8ze7lsr3k9Tb3SlCfyyC/xs8l7CEBAEygkAAAJlBIAAATKCQAgAkUEgDABAoJAGAChQQAMIFCAgCY4HohnTp1So8//rhKSko0btw4XXzxxXryySfV09Pj9qoAAKOI6yM1bNiwQc8884xeeOEFlZaWat++fbr//vsViUS0YsUKt1cHABglXC+kd999Vz/84Q+1YMECSdLkyZP10ksvad++fQPOn0qllEqleu8nk0m3IwEAAsD1Q3Zz5szRW2+9pQMHDkiSPvzwQ+3Zs0e33nrrgPPH43FFIpHeWywWczsSACAAQo7LowI6jqPHHntMGzZsUFZWlrq7u7V+/XqtWbNmwPkH2kOKxWJKJBLKy8tzMxrOgIEnM4/BVTOL13hmJZNJRSKRtD/HXT9kt23bNr344ovaunWrSktL1draqqqqKkWjUS1durTf/OFwWOFw2O0YAICAcb2QVq5cqdWrV2vx4sWSpOnTp+vQoUOKx+MDFhIAAJIH55COHz+uMWP6LjYrK4vLvgEAg3J9D2nhwoVav369iouLVVpaqv3792vjxo164IEH3F4VAGAUcf2ihq6uLj3xxBPavn27Ojs7FY1Gdffdd+vnP/+5cnJyhvz7kZ4Mw8hxwjfzuKghs3iNZ9ZIP8ddL6SzRSFlHm/WzKOQMovXeGaN9HOcsewAACZQSAAAEygkAIAJrl9lh+Dx+vi3l8fvOXafWUF9LnmNBwN7SAAAEygkAIAJFBIAwAQKCQBgAoUEADCBQgIAmEAhAQBMoJAAACZQSAAAEygkAIAJFBIAwAQKCQBgAoUEADCBQgIAmEAhAQBMoJAAACZQSAAAEygkAIAJFBIAwAQKCQBgAoUEADCBQgIAmJDtdwCMfo7jeLbsUCjk2bK9zB3UbeKloD6XXi8/yNslXewhAQBMoJAAACZQSAAAEygkAIAJFBIAwAQKCQBgAoUEADCBQgIAmJB2Ie3atUsLFy5UNBpVKBRSfX19n8cdx9G6desUjUY1btw4zZ07V5988olrgQEAo1PahXTs2DHNmDFDtbW1Az7+1FNPaePGjaqtrVVTU5OKiop0yy23qKur66zDAgBGr7SHDpo/f77mz58/4GOO42jTpk1au3atbr/9dknSCy+8oMLCQm3dulUPPvhgv79JpVJKpVK995PJZLqRAACjgKvnkNra2tTR0aHy8vLeaeFwWDfeeKP27t074N/E43FFIpHeWywWczMSACAgXC2kjo4OSVJhYWGf6YWFhb2PfdeaNWuUSCR6b+3t7W5GAgAEhCejfX93dFrHcc44Ym04HFY4HPYiBgAgQFzdQyoqKpKkfntDnZ2d/faaAAD4NlcLqaSkREVFRWpoaOidduLECTU2NmrWrFlurgoAMMqkfcju6NGjOnjwYO/9trY2tba2Kj8/X8XFxaqqqlJNTY2mTJmiKVOmqKamRuPHj9eSJUtcDQ4AGF3SLqR9+/Zp3rx5vferq6slSUuXLtXzzz+vn/3sZ/rf//6nn/70p/r666917bXX6s0331Rubq57qQEAo07IMfYbtslkUpFIRIlEQnl5eX7HgXHn0s87D1dQf8LcS0F9LqVgvsZH+jnOWHYAABMoJACACRQSAMAET74YC2SKl+cGgnjs3utlc34q886l55M9JACACRQSAMAECgkAYAKFBAAwgUICAJhAIQEATKCQAAAmUEgAABMoJACACRQSAMAECgkAYAKFBAAwgUICAJhAIQEATKCQAAAmUEgAABMoJACACRQSAMAECgkAYAKFBAAwgUICAJhAIQEATMj2OwBgleM4ni07FAp5tmwvc3u5bC95ub29FsTnM5lMKhKJpP137CEBAEygkAAAJlBIAAATKCQAgAkUEgDABAoJAGAChQQAMIFCAgCYkHYh7dq1SwsXLlQ0GlUoFFJ9fX3vYydPntSqVas0ffp0TZgwQdFoVPfdd5+OHDniamgAwOiTdiEdO3ZMM2bMUG1tbb/Hjh8/rpaWFj3xxBNqaWnRK6+8ogMHDui2225zJSwAYPQKOWcxdkQoFNL27du1aNGiM87T1NSka665RocOHVJxcfGQyzw95EQikVBeXt5IowGmBXXooKBi6KDMGunnuOdj2SUSCYVCIZ1//vkDPp5KpZRKpXrvJ5NJryMBAAzy9KKGb775RqtXr9aSJUvO2JLxeFyRSKT3FovFvIwEADDKs0I6efKkFi9erJ6eHm3evPmM861Zs0aJRKL31t7e7lUkAIBhnhyyO3nypO666y61tbXp7bffHvQYYjgcVjgc9iIGACBAXC+k02X02WefaefOnSooKHB7FQCAUSjtQjp69KgOHjzYe7+trU2tra3Kz89XNBrVHXfcoZaWFr322mvq7u5WR0eHJCk/P185OTnuJQcAjCppX/b9zjvvaN68ef2mL126VOvWrVNJScmAf7dz507NnTt3yOVz2TfOBVz2nVlc9p1ZGbvse+7cuYNuoCBuPACA/xjLDgBgAoUEADCBQgIAmEAhAQBM8HwsOwD9eXnxD1fwjS7n0vPJHhIAwAQKCQBgAoUEADCBQgIAmEAhAQBMoJAAACZQSAAAEygkAIAJFBIAwAQKCQBgAoUEADCBQgIAmEAhAQBMoJAAACZQSAAAEygkAIAJFBIAwAQKCQBgAoUEADCBQgIAmEAhAQBMoJAAACZk+x0AgLscx/Fs2aFQyLNle5nby2VL3m4XL1nLzR4SAMAECgkAYAKFBAAwgUICAJhAIQEATKCQAAAmUEgAABPSLqRdu3Zp4cKFikajCoVCqq+vP+O8Dz74oEKhkDZt2nRWIQEAo1/ahXTs2DHNmDFDtbW1g85XX1+v999/X9FodMThAADnjrRHapg/f77mz58/6DxffvmlHnroIb3xxhtasGDBiMMBAM4drg8d1NPTo4qKCq1cuVKlpaVDzp9KpZRKpXrvJ5NJtyMBAALA9YsaNmzYoOzsbC1fvnxY88fjcUUikd5bLBZzOxIAIABcLaTm5mY9/fTTev7554c9aN+aNWuUSCR6b+3t7W5GAgAEhKuFtHv3bnV2dqq4uFjZ2dnKzs7WoUOH9Oijj2ry5MkD/k04HFZeXl6fGwDg3OPqOaSKigrdfPPNfab94Ac/UEVFhe6//343VwUAGGXSLqSjR4/q4MGDvffb2trU2tqq/Px8FRcXq6CgoM/85513noqKinTJJZecfVoAwKiVdiHt27dP8+bN671fXV0tSVq6dKmef/5514IBAM4tIcfrn1JMUzKZVCQSUSKR4HwSYExQfzHWa9Z+edWKdD/HGcsOAGAChQQAMIFCAgCYQCEBAExwfSw7AKOXlxceBPmCiaBuF2vYQwIAmEAhAQBMoJAAACZQSAAAEygkAIAJFBIAwAQKCQBgAoUEADCBQgIAmEAhAQBMoJAAACZQSAAAEygkAIAJFBIAwAQKCQBgAoUEADCBQgIAmEAhAQBMoJAAACZQSAAAEygkAIAJFBIAwIRsvwN8l+M4kqRkMulzEgCZlEgkPFt2kD9PvNwuXkkmk4rFYr2f58NlrpC6urokSbFYzOckAICz0dXVpUgkMuz5Q066Feaxnp4eHTlyRLm5uQqFQkPOf7qJ29vblZeXl4GE7iB3ZgU1txTc7OTOLEu5HcdRV1eXotGoxowZ/pkhc3tIY8aM0UUXXZT23+Xl5fn+JIwEuTMrqLml4GYnd2ZZyZ3OntFpXNQAADCBQgIAmJC1bt26dX6HOFtZWVmaO3eusrPNHYEcFLkzK6i5peBmJ3dmBTX3aeYuagAAnJs4ZAcAMIFCAgCYQCEBAEygkAAAJlBIAAATKCQAgAn/D6JhIQH9P6k2AAAAAElFTkSuQmCC",
      "text/plain": [
       "PyPlot.Figure(PyObject <matplotlib.figure.Figure object at 0x136e6e6d0>)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "PyObject <matplotlib.image.AxesImage object at 0x13a272090>"
      ]
     },
     "execution_count": 49,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "N=5\n",
    "xgrid=linspace(0,1,N+1)\n",
    "Δx=1/N\n",
    "\n",
    "#assemling the finite different operator for the left hand side of the PDE\n",
    "T=spdiagm( (-ones(N-2), 2.0*ones(N-1),-ones(N-2)), [-1,0,1])\n",
    "I=speye(N-1,N-1)\n",
    "\n",
    "A=kron(T,I)+kron(I,T)\n",
    "\n",
    "using PyPlot\n",
    "spy(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaQAAAGgCAYAAADl3RMjAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzt3X9wFPX9x/FXIMkJMTkJP+7IEJlMzagYQEwsElSoQBxGRMdWUFHp1D9EAUmBUaN/mDptQu2I6FDTgToCWif9o8bSqQhx1CjNMIYIYwiO0oGWIElTnXgXkF4wfL5/ONzX4y7KJnvkk73nY2ZnYG+z+/nsXvLKJ/u+z6YZY4wAABhkwwa7AQAASAQSAMASBBIAwAoEEgDACgQSAMAKBBIAwAoEEgDACgQSAMAKBBIAwAoEEgDAClYH0osvvqiCggJddNFFKi4u1gcffDDYTRqQ999/X7feeqvy8vKUlpamN954I+Z1Y4wqKyuVl5enESNGaPbs2WptbR2k1vZPdXW1rr32WmVnZ2vcuHG6/fbb9emnn8ZsE4lEtHLlSo0ZM0ZZWVlauHChjh07Nkgt7p+amhpNmTJFOTk5ysnJ0YwZM7Rjx47o617o47mqq6uVlpam8vLy6Dov9LOyslJpaWkxSzAYjL7uhe/Lsz7//HPde++9Gj16tEaOHKmrr75azc3N0dcHu6/WBtKf//xnlZeX68knn9S+fft0ww03aP78+Tp69OhgN63fTp48qalTp2rjxo0JX3/mmWe0fv16bdy4UU1NTQoGg5o3b566u7svcEv7r6GhQcuXL9eePXtUX1+vb775RmVlZTp58mR0m/LyctXV1am2tla7d+/WiRMntGDBAvX29g5iy52ZMGGC1q1bp71792rv3r266aabdNttt0W/eb3Qx+9qamrSpk2bNGXKlJj1XunnVVddpfb29ujS0tISfc0L35eS1NXVpZkzZyojI0M7duzQwYMH9eyzz+qSSy6JbjPofTWW+vGPf2yWLVsWs+6KK64wjz/++CC1yF2STF1dXfT/Z86cMcFg0Kxbty667n//+5/x+/3mD3/4w2A00RWdnZ1GkmloaDDGGPPVV1+ZjIwMU1tbG93m888/N8OGDTNvvfXWYDXTFaNGjTJ//OMfPdfH7u5uU1hYaOrr682sWbPMqlWrjDHeuZZPPfWUmTp1asLXvPR9+dhjj5nrr7++z9dt6KuVI6Senh41NzerrKwsZn1ZWZkaGxsHqVXJdeTIEXV0dMT02efzadasWUO6z6FQSJKUm5srSWpubtbp06dj+pmXl6eioqIh28/e3l7V1tbq5MmTmjFjhuf6uHz5ct1yyy2aO3duzHov9fPQoUPKy8tTQUGB7rrrLh0+fFiSt74vt2/frpKSEt15550aN26cpk2bps2bN0dft6GvVgbSF198od7eXgUCgZj1gUBAHR0dg9Sq5DrbLy/12Rij1atX6/rrr1dRUZGkb/uZmZmpUaNGxWw7FPvZ0tKiiy++WD6fT8uWLVNdXZ0mTZrkqT7W1tbqo48+UnV1ddxrXunn9OnTtW3bNu3cuVObN29WR0eHSktL9eWXX3rq+/Lw4cOqqalRYWGhdu7cqWXLlumRRx7Rtm3bJNnxMyj9ghyln9LS0mL+b4yJW+c1XurzihUr9PHHH2v37t0/uO1Q7Ofll1+u/fv366uvvtJf/vIXLV26VA0NDX1uP9T62NbWplWrVmnXrl266KKLzvvrhlo/58+fH/335MmTNWPGDP3oRz/S1q1bdd1110nyxvflmTNnVFJSoqqqKknStGnT1NraqpqaGt1///3R7Qazr1aOkMaMGaPhw4fHpXJnZ2dcenvF2aoer/R55cqV2r59u959911NmDAhuj4YDKqnp0ddXV0x2w/FfmZmZuqyyy5TSUmJqqurNXXqVD3//POe6WNzc7M6OztVXFys9PR0paenq6GhQS+88ILS09MVCAQ80c9zZWVlafLkyTp06JCnvi/Hjx+vSZMmxay78soro4ViNvTVykDKzMxUcXGx6uvrY9bX19ertLR0kFqVXAUFBQoGgzF97unpUUNDw5DqszFGK1as0Ouvv6533nlHBQUFMa8XFxcrIyMjpp/t7e06cODAkOpnIsYYRSIRz/Rxzpw5amlp0f79+6NLSUmJlixZEv23F/p5rkgkok8++UTjx4/3zPelJM2cOTPuIxifffaZJk6cKMmSn0EXpHSiH2pra01GRoZ56aWXzMGDB015ebnJysoy//rXvwa7af3W3d1t9u3bZ/bt22ckmfXr15t9+/aZf//738YYY9atW2f8fr95/fXXTUtLi7n77rvN+PHjTTgcHuSWn7+HHnrI+P1+895775n29vbo8vXXX0e3WbZsmZkwYYJ5++23zUcffWRuuukmM3XqVPPNN98MYsudqaioMO+//745cuSI+fjjj80TTzxhhg0bZnbt2mWM8UYfE/lulZ0x3ujnmjVrzHvvvWcOHz5s9uzZYxYsWGCys7OjP2u88H1pjDEffvihSU9PN7/5zW/MoUOHzJ/+9CczcuRI8+qrr0a3Gey+WhtIxhjz+9//3kycONFkZmaaa665Jlo6PFS9++67RlLcsnTpUmPMt2WXTz31lAkGg8bn85kbb7zRtLS0DG6jHUrUP0nm5Zdfjm5z6tQps2LFCpObm2tGjBhhFixYYI4ePTp4je6HX/ziF9H35tixY82cOXOiYWSMN/qYyLmB5IV+Ll682IwfP95kZGSYvLw8c8cdd5jW1tbo6174vjzrb3/7mykqKjI+n89cccUVZtOmTTGvD3Zf04wx5sKMxQAA6JuV95AAAKmHQAIAWIFAAgBYgUACAFiBQAIAWIFAAgBYwepAikQiqqysVCQSGeymJBX99I5U6KNEP73Gln5a/TmkcDgsv9+vUCiknJycwW5O0tBP70iFPkr002ts6afVIyQAQOogkAAAVkja85BefPFF/e53v1N7e7uuuuoqbdiwQTfccMMPft2ZM2d0/PhxZWdnR5/jHg6Hk9VMK5ztH/0c+lKhjxL99Bq3+2mMUXd3t/Ly8jRsmINxTzImyDs7U/fmzZvNwYMHzapVq0xWVlZ0Vuvv09bW1ucEnSwsLCwsQ2dpa2tzlB1JKWqYPn26rrnmGtXU1ETXXXnllbr99tsTPgr5u0KhkC655BJHxwuFQue9rd/vH/A+vIZzAsBN4XBY+fn5+uqrr/r8+ZKI63+y6+npUXNzsx5//PGY9WVlZWpsbIzbPhKJxJQanv0znRNuVIV4uYKmvzgnAAbC6aPPXS9q+OKLL9Tb2xv3yNtAIBD3aFxJqq6ult/vjy75+fluNwkAMAQkrcru3GQ0xiRMy4qKCoVCoejS1taWrCYBACzm+p/sxowZo+HDh8eNhjo7O+NGTZLk8/nk8/kGdMxEQdfXrbG+1vc1tEzCLTbrcE4A2MD1EVJmZqaKi4tVX18fs76+vl6lpaVuHw4A4BFJ+RzS6tWrdd9996mkpEQzZszQpk2bdPToUS1btiwZhwMAeEBSAmnx4sX68ssv9fTTT6u9vV1FRUV68803NXHixGQcDgDgAdZNrnp2kr+Bctot7pfE45wA6I/+TtbKXHYAACskbS67gUqUrE4+ZOX0t3snlWapMkKg+g7AhcQICQBgBQIJAGAFAgkAYAUCCQBgBWuLGhJxepM9WcdM9Zv6FDsASAZGSAAAKxBIAAArEEgAACsQSAAAKxBIAAArWFtll2iCVafT/iTitCIv0b6pMkuM6ZcADAQjJACAFQgkAIAVCCQAgBUIJACAFQgkAIAVrK2yc/KAvmRWbDmpEKP6LjHmAwRwPhghAQCsQCABAKxAIAEArEAgAQCsQCABAKxgbZVdIm7MlebGU2edVohRfRePcwLgXIyQAABWIJAAAFYgkAAAViCQAABWGFJFDX1xY2oaN4odnOKBdvEodgBSFyMkAIAVCCQAgBUIJACAFQgkAIAVCCQAgBU8UWWXiFvVWk6quJxW5Dk5ZqpXmVF9B3gfIyQAgBUIJACAFQgkAIAVCCQAgBUIJACAFTxbZdcXm6q1nBzTpnbbhPkAAe9ghAQAsAKBBACwAoEEALACgQQAsAKBBACwguNAev/993XrrbcqLy9PaWlpeuONN2JeN8aosrJSeXl5GjFihGbPnq3W1lbXGpwsxpiES1paWtzidB9OuXHMRPtI5tNvbcM5AYYex4F08uRJTZ06VRs3bkz4+jPPPKP169dr48aNampqUjAY1Lx589Td3T3gxgIAvMvx55Dmz5+v+fPnJ3zNGKMNGzboySef1B133CFJ2rp1qwKBgF577TU9+OCDA2stAMCzXL2HdOTIEXV0dKisrCy6zufzadasWWpsbEz4NZFIROFwOGYBAKQeVwOpo6NDkhQIBGLWBwKB6Gvnqq6ult/vjy75+fluNgkAMEQkpcru3BvFZ28oJ1JRUaFQKBRd2traktEkAIDlXJ3LLhgMSvp2pDR+/Pjo+s7OzrhR01k+n08+n8/NZrjKjae3Op2Hzg3M8RaP+QABu7k6QiooKFAwGFR9fX10XU9PjxoaGlRaWurmoQAAHuN4hHTixAn985//jP7/yJEj2r9/v3Jzc3XppZeqvLxcVVVVKiwsVGFhoaqqqjRy5Ejdc889rjYcAOAtjgNp7969+slPfhL9/+rVqyVJS5cu1ZYtW/Too4/q1KlTevjhh9XV1aXp06dr165dys7Odq/VAADPSTOW/aE8HA7L7/crFAopJydnsJuTkFv3HJzcQ3LrMqXyPaS+cA8JcFd/f46n3AP63ODWzXEnP/CcFkC4ccxUQbEDYAcmVwUAWIFAAgBYgUACAFiBQAIAWIFAAgBYgSo7F9lUrUWF2MAx/RJwYTFCAgBYgUACAFiBQAIAWIFAAgBYgUACAFiBKrsLwI3qO7ce8keF2MC58dBGAPEYIQEArEAgAQCsQCABAKxAIAEArEAgAQCsQJXdIHJjrjS3qu8wMDbNYwgMVYyQAABWIJAAAFYgkAAAViCQAABWsLaowe/3x61LlRvBbvTTyT6cFkA4bV8qT1dEsQNw/hghAQCsQCABAKxAIAEArEAgAQCsQCABAKxgbZVdKBRSTk5OzDoqk+zg9DrwQLt4VN8B8RghAQCsQCABAKxAIAEArEAgAQCsQCABAKxgbZVdIm480A7x3HrIn5MKMarMEuM9jlTGCAkAYAUCCQBgBQIJAGAFAgkAYAUCCQBghSFVZdcX5kpLDreq75J5zFS5nrzHkQoYIQEArEAgAQCsQCABAKxAIAEArEAgAQCs4CiQqqurde211yo7O1vjxo3T7bffrk8//TRmm0gkopUrV2rMmDHKysrSwoULdezYMVcbfT6MMQmXtLS0hEsqc3pO+jq3Ts+5G8dM5evJOYHXOAqkhoYGLV++XHv27FF9fb2++eYblZWV6eTJk9FtysvLVVdXp9raWu3evVsnTpzQggUL1Nvb63rjAQDekWYG8KGF//73vxo3bpwaGhp04403KhQKaezYsXrllVe0ePFiSdLx48eVn5+vN998UzfffHPcPiKRiCKRSPT/4XBY+fn5CoVCysnJ6W/T+sRnN+Il+5w4+e3c6TG5nvE4Jxhs4XBYfr/f8c/xAd1DCoVCkqTc3FxJUnNzs06fPq2ysrLoNnl5eSoqKlJjY2PCfVRXV8vv90eX/Pz8gTQJADBE9TuQjDFavXq1rr/+ehUVFUmSOjo6lJmZqVGjRsVsGwgE1NHRkXA/FRUVCoVC0aWtra2/TQIADGH9njpoxYoV+vjjj7V79+4f3PbsjdZEfD6ffD5ff5vhGFPTxLPpnDg9Jg+0i2fT9QSc6NcIaeXKldq+fbveffddTZgwIbo+GAyqp6dHXV1dMdt3dnYqEAgMrKUAAE9zFEjGGK1YsUKvv/663nnnHRUUFMS8XlxcrIyMDNXX10fXtbe368CBAyotLXWnxQAAT3L0J7vly5frtdde01//+ldlZ2dH7wv5/X6NGDFCfr9fDzzwgNasWaPRo0crNzdXa9eu1eTJkzV37tykdAAA4A2OAqmmpkaSNHv27Jj1L7/8sn7+859Lkp577jmlp6dr0aJFOnXqlObMmaMtW7Zo+PDhrjQYAOBNA/ocUjL0t359oLjhG8+tc+LGLAFuHDOVr6XEexwXTn9/jnviAX1uoFornlvVWk4eLteXZB6T65m673HYhclVAQBWIJAAAFYgkAAAViCQAABWIJAAAFagyu4HUK0Vz43qO6f7cANzvCXGexy2YIQEALACgQQAsAKBBACwAoEEALACgQQAsAJVdv1AtVZibsyVlsyJW9146uz3be8lnBMMBkZIAAArEEgAACsQSAAAKxBIAAArEEgAACtQZeciKpMSs2WuNDeeOtvXflL5Wkq8x+EORkgAACsQSAAAKxBIAAArEEgAACtQ1HABcCM4nlvnxI0H/blxzFS+lhLvcbiDERIAwAoEEgDACgQSAMAKBBIAwAoEEgDAClTZDSKmpomXzOo7J5V3TlFllhjvcTjBCAkAYAUCCQBgBQIJAGAFAgkAYAUCCQBgBarsLMRcafHcqGJzeq6cVOU5nWsv1fEeRyKMkAAAViCQAABWIJAAAFYgkAAAViCQAABWoMpuiGCutMRsmSst1a+DG3iPgxESAMAKBBIAwAoEEgDACgQSAMAKjgKppqZGU6ZMUU5OjnJycjRjxgzt2LEj+nokEtHKlSs1ZswYZWVlaeHChTp27Jjrjcb/M8YkXNLS0hIuqcKNc9LXuXUi1a+DG3iPpw5HgTRhwgStW7dOe/fu1d69e3XTTTfptttuU2trqySpvLxcdXV1qq2t1e7du3XixAktWLBAvb29SWk8AMA70swAaydzc3P1u9/9Tj/72c80duxYvfLKK1q8eLEk6fjx48rPz9ebb76pm2+++bz2Fw6H5ff7FQqFlJOTM5CmpTRKZeO5dU7c+C08la+DW3iP26u/P8f7fQ+pt7dXtbW1OnnypGbMmKHm5madPn1aZWVl0W3y8vJUVFSkxsbGPvcTiUQUDodjFgBA6nEcSC0tLbr44ovl8/m0bNky1dXVadKkSero6FBmZqZGjRoVs30gEFBHR0ef+6uurpbf748u+fn5znsBABjyHAfS5Zdfrv3792vPnj166KGHtHTpUh08eLDP7c/efOxLRUWFQqFQdGlra3PaJACABzieOigzM1OXXXaZJKmkpERNTU16/vnntXjxYvX09KirqytmlNTZ2anS0tI+9+fz+eTz+frRdHwfW6bUsYlbU9M4ebgckof3uPcM+HNIxhhFIhEVFxcrIyND9fX10dfa29t14MCB7w0kAAAkhyOkJ554QvPnz1d+fr66u7tVW1ur9957T2+99Zb8fr8eeOABrVmzRqNHj1Zubq7Wrl2ryZMna+7cuclqPwDAIxwF0n/+8x/dd999am9vl9/v15QpU/TWW29p3rx5kqTnnntO6enpWrRokU6dOqU5c+Zoy5YtGj58eFIaDwDwjgF/DsltfA4pufj7ejw3Ps/i9B5Sqp/zZOI9Pvgu+OeQAABwEw/oSzFOKsRS5bdKN/qZ7Nke3BitpfL1TPVzMlQwQgIAWIFAAgBYgUACAFiBQAIAWIFAAgBYgSo7uDbHG5LHyWdruJ7xOCdDAyMkAIAVCCQAgBUIJACAFQgkAIAVCCQAgBWoskOfqExKHqfn1sm2blTfpcq15D1uF0ZIAAArEEgAACsQSAAAKxBIAAArUNQAx7gRnDxuFDu4ccxUv5a8xwcHIyQAgBUIJACAFQgkAIAVCCQAgBUIJACAFaiyg2uYmub8uTXtj5N99yXRvqkyS4z3eHIxQgIAWIFAAgBYgUACAFiBQAIAWIFAAgBYgSo7JB1zpcWzqYrNSYWYTe22Ce9xdzBCAgBYgUACAFiBQAIAWIFAAgBYgUACAFiBKjsMCqq1EnNjrjQ3njrr1lx7qXw9OSfOMUICAFiBQAIAWIFAAgBYgUACAFiBogZYhRvBibkxNY0bxQ5O8UC7eLzH+8YICQBgBQIJAGAFAgkAYAUCCQBgBQIJAGCFAQVSdXW10tLSVF5eHl0XiUS0cuVKjRkzRllZWVq4cKGOHTs24IYitRljEi5paWkJl1Tg1jnpaz9O9u3GMVP5Wkq8x6UBBFJTU5M2bdqkKVOmxKwvLy9XXV2damtrtXv3bp04cUILFixQb2/vgBsLAPCufgXSiRMntGTJEm3evFmjRo2Krg+FQnrppZf07LPPau7cuZo2bZpeffVVtbS06O2333at0QAA7+lXIC1fvly33HKL5s6dG7O+ublZp0+fVllZWXRdXl6eioqK1NjYmHBfkUhE4XA4ZgEApB7HMzXU1tbqo48+UlNTU9xrHR0dyszMjBk1SVIgEFBHR0fC/VVXV+tXv/qV02YAADzG0Qipra1Nq1at0quvvqqLLrrovL/u7I25RCoqKhQKhaJLW1ubkyYBADzCUSA1Nzers7NTxcXFSk9PV3p6uhoaGvTCCy8oPT1dgUBAPT096urqivm6zs5OBQKBhPv0+XzKycmJWYDz5aQyKVXYVK3l5Jg2tdsmqfQedxRIc+bMUUtLi/bv3x9dSkpKtGTJkui/MzIyVF9fH/2a9vZ2HThwQKWlpa43HgDgHY7uIWVnZ6uoqChmXVZWlkaPHh1d/8ADD2jNmjUaPXq0cnNztXbtWk2ePDmuAAIAgO9y/fETzz33nNLT07Vo0SKdOnVKc+bM0ZYtWzR8+HC3DwUA8JA0Y9nDNsLhsPx+v0KhEPeT0G+p/LydvrjxvB237lO4cUyup73v8f7+HGcuOwCAFXhiLDzJjSeseo0bT29166mzyTxmKl/PoX5OGCEBAKxAIAEArEAgAQCsQCABAKxAIAEArECVHVIG1VqJuVGt5Vb1nRNuVA16zVB/jzNCAgBYgUACAFiBQAIAWIFAAgBYgaIGpLyhfiM4Gdw6J8mcuNWNY6aKofIeZ4QEALACgQQAsAKBBACwAoEEALACgQQAsAJVdkAfmJomnk3VWrZViA1Ftr3HGSEBAKxAIAEArEAgAQCsQCABAKxAIAEArECVHeCQGw+08xo3qu/ceshfKldBumWw3uOMkAAAViCQAABWIJAAAFYgkAAAViCQAABWoMoOcIFNc7zZxI250tyqvsPAXIjrwAgJAGAFAgkAYAUCCQBgBQIJAGAFihqAJKLYITE3+ulkH05vvDttXypPV5Son+FwWH6/3/G+GCEBAKxAIAEArEAgAQCsQCABAKxAIAEArECVHTAIqL6zm9PrwEMb3cEICQBgBQIJAGAFAgkAYAUCCQBgBQIJAGAFR4FUWVmptLS0mCUYDEZfN8aosrJSeXl5GjFihGbPnq3W1lbXGw14lTEm4XLu9x0Pp3Omr/PqVKLr8H3VdOd7Lbme33I8QrrqqqvU3t4eXVpaWqKvPfPMM1q/fr02btyopqYmBYNBzZs3T93d3a42GgDgPY4DKT09XcFgMLqMHTtW0re/DWzYsEFPPvmk7rjjDhUVFWnr1q36+uuv9dprr7necACAtzgOpEOHDikvL08FBQW66667dPjwYUnSkSNH1NHRobKysui2Pp9Ps2bNUmNjY5/7i0QiCofDMQsAIPU4CqTp06dr27Zt2rlzpzZv3qyOjg6Vlpbqyy+/VEdHhyQpEAjEfE0gEIi+lkh1dbX8fn90yc/P70c3AABDnaNAmj9/vn76059q8uTJmjt3rv7+979LkrZu3Rrd5tybc2dv4vWloqJCoVAourS1tTlpEgDAIwZU9p2VlaXJkyfr0KFD0Wq7c0dDnZ2dcaOm7/L5fMrJyYlZAMSiWis53Kq+c+OYXM8BBlIkEtEnn3yi8ePHq6CgQMFgUPX19dHXe3p61NDQoNLS0gE3FADgbY5m+167dq1uvfVWXXrppers7NSvf/1rhcNhLV26VGlpaSovL1dVVZUKCwtVWFioqqoqjRw5Uvfcc0+y2g8A8AhHgXTs2DHdfffd+uKLLzR27Fhdd9112rNnjyZOnChJevTRR3Xq1Ck9/PDD6urq0vTp07Vr1y5lZ2cnpfEAAO9IM5Y9nCMcDsvv9ysUCnE/CfgePG8neZzeu3HjnHvpevb35zhz2QEArMATY4EhiqfOnj83ngDbn/27ccxUup6MkAAAViCQAABWIJAAAFYgkAAAVqCoAfAYbo7Hs+mcuFVgkWg/Q/1aMkICAFiBQAIAWIFAAgBYgUACAFiBQAIAWIEqOyBFpFK11vlyq/ou0XqnE7Qm85hD5XoyQgIAWIFAAgBYgUACAFiBQAIAWIFAAgBYgSo7IMV5sVproNyovnO6DzfYNGdffzBCAgBYgUACAFiBQAIAWIFAAgBYgUACAFiBKjsAcYZ6tVayuDEfoNNz5aQqz42nzn7f9snGCAkAYAUCCQBgBQIJAGAFAgkAYAUCCQBgBarsAJy3oVKtdaHZMh+gG0+d7Ws/F+JaMkICAFiBQAIAWIFAAgBYgUACAFiBogYAA0axQzy3zokbD/pz45gX4loyQgIAWIFAAgBYgUACAFiBQAIAWIFAAgBYgSo7AElj29Q0Nkhm9Z2Tyjun3Kj2+yGMkAAAViCQAABWIJAAAFYgkAAAViCQAABWcBxIn3/+ue69916NHj1aI0eO1NVXX63m5ubo68YYVVZWKi8vTyNGjNDs2bPV2trqaqMBDG3GmLglLS0t4ZIqEp0Tp+elr3043beT65Bov6FQqF/nwFEgdXV1aebMmcrIyNCOHTt08OBBPfvss7rkkkui2zzzzDNav369Nm7cqKamJgWDQc2bN0/d3d39aiAAIDWkGQfF/48//rj+8Y9/6IMPPkj4ujFGeXl5Ki8v12OPPSZJikQiCgQC+u1vf6sHH3ww7msikYgikUj0/+FwWPn5+QqFQsrJyXHaHwBDVCrPDP59knlenIxAnRwvHA7L7/c7/jnuaIS0fft2lZSU6M4779S4ceM0bdo0bd68Ofr6kSNH1NHRobKysug6n8+nWbNmqbGxMeE+q6ur5ff7o0t+fr6TJgEAPMJRIB0+fFg1NTUqLCzUzp07tWzZMj3yyCPatm2bJKmjo0N40tkcAAACrElEQVSSFAgEYr4uEAhEXztXRUWFQqFQdGlra+tPPwAAQ5yjqYPOnDmjkpISVVVVSZKmTZum1tZW1dTU6P77749ud+4w8OzNs0R8Pp98Pp/TdgMAPMbRCGn8+PGaNGlSzLorr7xSR48elSQFg0FJihsNdXZ2xo2aAOC73Kgy8yIn5yWZLsR1cBRIM2fO1Keffhqz7rPPPtPEiRMlSQUFBQoGg6qvr4++3tPTo4aGBpWWlrrQXACAVzn6k90vf/lLlZaWqqqqSosWLdKHH36oTZs2adOmTZK+TdDy8nJVVVWpsLBQhYWFqqqq0siRI3XPPfckpQMAAG9wFEjXXnut6urqVFFRoaeffloFBQXasGGDlixZEt3m0Ucf1alTp/Twww+rq6tL06dP165du5Sdne164wEA3uHoc0gXQn/r1wF4E59PSsyNZ0q5cQ8o0TH7+3OcB/QBsJpbD7TzGicP6HPykL/v289At/0hTK4KALACgQQAsAKBBACwAoEEALACgQQAsAJVdgCGJCcVYqlceSe5U313IaZrYoQEALACgQQAsAKBBACwAoEEALACgQQAsIJ1VXZnqzvC4fAgtwTAUBQKheLWpfrPk0TnRHJ2XvraR1/7zc/Pd1zdaF0gdXd3S5Ly8/MHuSUAgIHo7u6W3+8/7+2te/zEmTNndPz4cWVnZ6u7u1v5+flqa2vz9KMozv42QT+HvlToo0Q/vcbtfhpj1N3drby8PA0bdv53hqwbIQ0bNkwTJkyQ9P8fxMrJyfH0m+Es+ukdqdBHiX56jZv9dDIyOouiBgCAFQgkAIAVhldWVlYOdiO+z/DhwzV79mylp1v310VX0U/vSIU+SvTTa2zop3VFDQCA1MSf7AAAViCQAABWIJAAAFYgkAAAViCQAABWIJAAAFb4P/Si0JucN8SoAAAAAElFTkSuQmCC",
      "text/plain": [
       "PyPlot.Figure(PyObject <matplotlib.figure.Figure object at 0x11842ba90>)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "PyObject <matplotlib.image.AxesImage object at 0x136ba9390>"
      ]
     },
     "execution_count": 46,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A=kron(kron(T,I),I)+ kron(I,kron(I,T))+ kron(kron(I,T),I)\n",
    "\n",
    "size(A)\n",
    "\n",
    "using PyPlot\n",
    "spy(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "size of the matrix:(970299, 970299)\n",
      "Optimal direct solver: \n",
      "183.711107 seconds (65 allocations: 20.119 GiB, 0.14% gc time)\n",
      "\n",
      "\n",
      "\n",
      "\n",
      "\n",
      "\n"
     ]
    }
   ],
   "source": [
    "N=100\n",
    "xgrid=linspace(0,1,N+1)\n",
    "Δx=1/N\n",
    "\n",
    "f=(x,y)-> x.^2+y+z.^2\n",
    "\n",
    "#assemling the finite different operator for the left hand side of the PDE\n",
    "T=spdiagm( (-ones(N-2), 2.0*ones(N-1),-ones(N-2)), [-1,0,1])\n",
    "I=speye(N-1,N-1)\n",
    "\n",
    "#discretize the Laplace operator in 2D\n",
    "#A=kron(T,I)+kron(I,T)\n",
    "\n",
    "#discretize the Laplace operator in 3D\n",
    "A=kron(kron(T,I),I)+ kron(I,kron(I,T))+ kron(kron(I,T),I)\n",
    "\n",
    "F=ones(size(A,1))\n",
    "\n",
    "b=F.*Δx^2\n",
    "\n",
    "U=copy(F)\n",
    "Uj=copy(F)\n",
    "Ugs=copy(F)\n",
    "Usor=copy(F)\n",
    "\n",
    "println(\"size of the matrix:\", size(A))\n",
    "\n",
    "println(\"Optimal direct solver: \")\n",
    "@time U=A\\b\n",
    "\n",
    "#println(\"LU:\")\n",
    "#F=lufact(A)\n",
    "#@time U=inv(F[:U])*inv(F[:L])*b\n",
    "\n",
    "println(\"\\n\")\n",
    "#@time GaussSeidelIterations!(Ugs, A, b, 1e-12)\n",
    "\n",
    "println(\"\\n\")\n",
    "#@time JacobiIterations!(Uj, A, b, 1e-12)\n",
    "\n",
    "println(\"\\n\")\n",
    "#@time SOR!(Usor, A, b, 1.5, 1e-12)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## Convergence of iterative algorithms\n",
    "\n",
    "* linearly convergent algorithms: (**the error in the next step is a linear function of the current error**)\n",
    "$$ e_{n+1} = a e_n $$\n",
    "    \n",
    "    * number of correct digits is a linear function of the number of iterations: $$ \\textrm{num of accurate digits} = \\alpha_0 n $$\n",
    "\n",
    "\n",
    "* quadratically convergent algorithms: (**the error in the next step is a quadratic function of the current error**)\n",
    "$$ e_{n+1} = a e_n^2 $$\n",
    "\n",
    "    * number of correct digits doubles after each iterations: $$ \\textrm{num of accurate digits} = \\alpha_1 2^n$$\n",
    "\n",
    "\n",
    "* p-th order convergent algorithms:\n",
    "$$ e_{n+1} = a e_n^p $$\n",
    "\n",
    "### Possible confusion with ***order of accuracy***\n",
    "\n",
    "(For linearly convergent algorithms)\n",
    "$$e_n \\le Const\\cdot \\frac{1}{n^p}$$\n",
    "\n",
    "**p-th order accurate algorithm**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The conjugate gradient (CG) method of Hestenes and Stiefel (1952)\n",
    "\n",
    "$$Ax=b$$\n",
    "* It is optimal for symmetric positive definite systems: $A^T =A,~~~x^TAx>0,~\\forall x\\neq 0$\n",
    "* It is a type of **Krylov methods**\n",
    "* In each iteration, say iteration $k$, the algorithm chooses the next guess, $x^{(k+1)}$, by \n",
    "  * minimizes the residual, $r^{(k+1)}:=b-Ax^{(k+1)}$, \n",
    "  * in a special norm built from $A$ and $x^{(j)},j=1,2,\\cdots,k$.\n",
    "* Very simple, elegant and fast algorithm.\n",
    "\n",
    "* We have seen that solving teh model Poisson's problem in 1 and 2 dimenions, the backlash operator is faster than CG. But the situation changes when we solve the Poisson's equation in 3D, for which the matrix has more strides away from teh diagonal.\n",
    "\n",
    "* See Algorithm 11.15 in GGK or this workbook (instead of BF) for a more comprehensive view of teh CG algorithm.\n",
    "\n",
    "**Learning objectives**:\n",
    "\n",
    "* Know when CG is applicable\n",
    "* Theorem 7.31 in Burden and Faires: Under which conditions will solving the given linear system be equivalent to solving a related minimization problem.\n",
    "    * This theorem provides a hint to the \"G\" in CG.\n",
    "    * Know how to formula the method of steepest descent for the minimization problem\n",
    "* Theorem 7.32\n",
    "* Have some idea how Theorems 7.31-7.32 are related to CG."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "CG (generic function with 2 methods)"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "CG(A,b;kwargs...) = CG(x -> A*x,b;kwargs...)\n",
    "    \n",
    "function CG(A::Function, b; tol=1e-2, maxItr=1000)\n",
    "\n",
    "    n = length(b)\n",
    "\n",
    "    #initial guess\n",
    "    x=zeros(n)\n",
    "    \n",
    "    #initial residual\n",
    "    r=b-A(x)\n",
    "\n",
    "    p = copy(r)\n",
    "    \n",
    "    res0  = norm(r)\n",
    "\n",
    "    itr   = 1 \n",
    "    \n",
    "    for itr=1:maxItr\n",
    "        Ap = A(p)\n",
    "        pAp= dot(p,Ap)\n",
    "        α  = dot(r,p)/pAp\n",
    "\n",
    "        #update the approximate solution\n",
    "        x.+=α.*p\n",
    "        \n",
    "        #update the residual vector\n",
    "        r.-=α.*Ap\n",
    "            \n",
    "        #break the loop when tolerance is met\n",
    "        if norm(r)/res0 <= tol\n",
    "            break\n",
    "        end\n",
    "\n",
    "        #update search direction\n",
    "        β = -dot(Ap,r)/pAp           \n",
    "        p.=r.+β.*p\n",
    "    end\n",
    "\n",
    "    println(\"cg termines after \", itr, \" steps, with residual=\", norm(r))\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "cg termines after 371 steps, with residual=9.207057234407962e-16\n",
      "  7.314909 seconds (938 allocations: 2.711 GiB, 3.09% gc time)\n"
     ]
    }
   ],
   "source": [
    "@time uu=CG(A,b, tol=1e-14)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "using KrylovMethods"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "cg achieved desired tolerance at iteration 279. Residual norm is 9.69e-11.\n",
      "  4.136621 seconds (663 allocations: 29.637 MiB, 0.10% gc time)\n",
      "cg achieved desired tolerance at iteration 371. Residual norm is 9.21e-15.\n",
      "  5.491275 seconds (894 allocations: 29.650 MiB, 0.02% gc time)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(970299, 970299)"
      ]
     },
     "execution_count": 54,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "@time uu=KrylovMethods.cg(A, b, tol=1e-10, maxIter=1500, out=1);\n",
    "\n",
    "@time uu=KrylovMethods.cg(A, b, tol=1e-10, maxIter=1500, out=1);\n",
    "\n",
    "@time uu=KrylovMethods.cg(A, b, tol=1e-14, maxIter=1500, out=1);\n",
    "\n",
    "size(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5×5 Array{Float64,2}:\n",
       " 0.1  1.0  0.0  0.0  0.0\n",
       " 0.0  0.1  1.0  0.0  0.0\n",
       " 0.0  0.0  0.1  1.0  0.0\n",
       " 0.0  0.0  0.0  0.1  1.0\n",
       " 0.0  0.0  0.0  0.0  0.1"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "J=diagm(0.1*ones(5))+diagm(ones(4),1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5×5 Array{Float64,2}:\n",
       " 0.01  0.2   1.0   0.0   0.0 \n",
       " 0.0   0.01  0.2   1.0   0.0 \n",
       " 0.0   0.0   0.01  0.2   1.0 \n",
       " 0.0   0.0   0.0   0.01  0.2 \n",
       " 0.0   0.0   0.0   0.0   0.01"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "J^2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5×5 Array{Float64,2}:\n",
       " 0.001  0.03   0.3    1.0    0.0  \n",
       " 0.0    0.001  0.03   0.3    1.0  \n",
       " 0.0    0.0    0.001  0.03   0.3  \n",
       " 0.0    0.0    0.0    0.001  0.03 \n",
       " 0.0    0.0    0.0    0.0    0.001"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "J^3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5×5 Array{Float64,2}:\n",
       " 0.0001  0.004   0.06    0.4     1.0   \n",
       " 0.0     0.0001  0.004   0.06    0.4   \n",
       " 0.0     0.0     0.0001  0.004   0.06  \n",
       " 0.0     0.0     0.0     0.0001  0.004 \n",
       " 0.0     0.0     0.0     0.0     0.0001"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "J^4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Julia 0.6.2",
   "language": "julia",
   "name": "julia-0.6"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "0.6.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}