Perturbation_Analysis_Parabolic_PDE.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Perturbation Analysis for Parabolic Partial Differential Equations\n",
    "\n",
    "<h2> Cristopher Arenas</h2>\n",
    "<h3> M.Sc. Student </h3>\n",
    "<h3> Universidad Técnica Federico Santa Maria, Campus San Joaquín </h3>\n",
    "cristopher.arenas@alumnos.usm.cl"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy as sc\n",
    "from scipy.linalg import toeplitz\n",
    "import matplotlib.pyplot as plt\n",
    "from ipywidgets import *\n",
    "from support import *\n",
    "\n",
    "%matplotlib inline\n",
    "\n",
    "import matplotlib as mpl\n",
    "mpl.rcParams['lines.linewidth'] = 2\n",
    "mpl.rcParams['font.size'] = 15"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Introduction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Commonly, exist some phenomena that can be represented as mathematical models, using differential equations. In some cases, simplificationes in the moment of making the model can reduce the difficulty of the representation, but in change, the model could not consider the effect of small factors that are occurring."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Parabolic Partial Differential Equations model different process in the scope of Engineering and Science. One representative Parabolic Partial Differential Equation is the following:\n",
    "$$\n",
    "\\frac{\\partial u(x,t)}{\\partial t} = c\\,\\frac{\\partial^2 u(x,t)}{\\partial x^2} + b\\,\\frac{\\partial u(x,t)}{\\partial x} \n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Tipically, some process involved in parabolic partial differential equations are related with some quantity or a compond that travels along a fluid and changes its concentration. The process that describe the movement of the compond across the fluid is called **advection** and has a speed asociated. The process that describe the spreading of the concentration of compond in the fluid is called **diffusion** and has a coefficient of diffusion associated."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Is interesting to analize what happens when the process are acting at the same time at different rythms and the implications that carry when time advances. Next, a numerical method will be explained in order to obtain a good estimation of a general equation that models advection-difusion process in differents orders of magnitude."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Method of multiple scales"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "The method of multiple scales is a general method applicable in a wide range of problems. In this problems there are two processes acting at the same time, each one with its own scale. \n",
    "The problem of determinate $u(x,t)$ usually has an approximated value using the approximation:\n",
    "\n",
    "$$\n",
    "u(x,t, \\varepsilon) \\approx u_0(x,t) + \\varepsilon \\, u_1(x,t)\n",
    "$$\n",
    "\n",
    "This expansion works fine with a fixed value of $t$ and $\\varepsilon$ near to zero. However, there are problems if $t \\geq \\mathcal{O}(\\varepsilon ^{-1})$, since the approximation is no longer valid (secular terms can appear). There are two processes acting simultaneously that has its own scales. So, it is necessary to introduce two time variables:\n",
    "\n",
    "- $\\tau = t$ the _fast time_\n",
    "- $T = \\varepsilon t$ the _slow time_\n",
    "\n",
    "With this variables, the solution proposed by the method has the form:\n",
    "\n",
    "$$\n",
    "u(x,t,\\varepsilon) \\approx u(x,\\tau, T, \\varepsilon)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "In this maner, variations of time must be deteminated in function of the variations of the _fast time_ and _slow time_. In terms of algebra:\n",
    "\n",
    "$$\n",
    "\\frac{\\mathrm{d}}{\\mathrm{d}t} = \\left( \\frac{ \\partial }{\\partial \\tau} \\right) + \\varepsilon \\left( \\frac{\\partial}{\\partial T} \\right)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, the approximation for $u(x,t)$ that now is asymptotic for $T = \\mathcal{O}(1)$ is:\n",
    "\n",
    "$$\n",
    "u(x,t,\\varepsilon) \\approx u_0(x,\\tau,T) + \\varepsilon\\, u_1(x,\\tau,T)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## A diffusion-advection equation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Consider the following problem that shows advection around a periodic domain $0 < \\theta < 2\\pi$ with some small diffusion. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\\begin{align*}\n",
    "\\frac{\\partial u(\\theta,t)}{\\partial t} + \\frac{\\partial}{\\partial \\theta} (\\omega(\\theta)\\,u(\\theta,t)) &= \\varepsilon\\frac{\\partial^2u(\\theta,t)}{\\partial \\theta^2}\\\\\n",
    "u(\\theta,0) &= F(\\theta)\\\\\n",
    "u(0,t) &= u(2\\pi,t)\n",
    "\\end{align*}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "where $\\omega(\\theta)$ is given over $0 \\leq \\theta \\leq 2\\pi$, periodic and positive. $F(\\theta)$ is a known function and $\\varepsilon$ is a tiny value."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Application of the Method of multiple scales"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " As can be seen in previous equation, there are two process acting at simultaneously: advection on the fast time scale and difusion on the slow scale. According to the method, the next expansion is used:\n",
    " \n",
    "$$\n",
    "u(\\theta,t,\\varepsilon) \\approx u_0(\\theta,\\tau,T) + \\varepsilon\\, u_1(\\theta,\\tau,T)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Sustituing into the equation and comparing coefficients there are derivated problems to solve."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**First problem**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\n",
    "\\mathcal{O}(1):\n",
    "$\n",
    "\n",
    "$$\n",
    "\\frac{\\partial u_0}{\\partial \\tau} + \\frac{\\partial}{\\partial \\theta}(\\omega(\\theta) \\, u_0) = 0\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first problem can be rewrited as:\n",
    "\n",
    "$$\n",
    "\\frac{1}{\\omega} \\left( \\frac{\\partial}{\\partial \\tau} + \\omega \\frac{\\partial}{\\partial \\theta}\\right) (\\omega\\, u_0) = 0\n",
    "$$\n",
    "\n",
    "This is an advection problem that depends only of the fast time, no slow time are implicated. An observed moving at  a speed $\\omega$ will see the quantity $\\omega \\, u_0$ remain constant. If the value of $\\omega(\\theta)$ was $\\omega(\\theta) = \\omega$, where $\\omega$ is constant, the solution is reduced to $u_0(\\theta,\\tau,T) = F(\\theta-\\omega\\tau)$, where $F(\\theta)$ is the initial value. However, in a general scenario, is necessary to know how much advanced is done after a time $\\tau$. Consider the following initial value problem:\n",
    "\n",
    "\\begin{align*}\n",
    "\\frac{\\mathrm{d}\\theta}{\\mathrm{d}t} &= \\omega(\\theta)\\\\\n",
    "\\theta(0) &= 0\n",
    "\\end{align*}\n",
    "\n",
    "Supouse that $\\Theta(t)$ is the solution of this problem and, since $\\theta$ is a periodic variable, the function $\\Theta(t)$ has a period $P$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The next function solves the Initial Value Problem explained before using the Forward Euler Method. The return of this function is the function $\\Theta(t)$ constructed from a cubic spline. Also, the period $P$ of the function $\\Theta(t)$ was computed with the reconstructed funcion form cubic spline [2] modulo $2\\pi$. A verification of the period was employed using the Bisection method [1]. In order to determinate a value $s$ such that $\\Theta(s) = \\hat{\\theta}$, since $\\Theta(0) = 0$. Let $T(s)$ be the function:\n",
    "\n",
    "$$\n",
    "T(s) = \\Theta(s) - \\hat{\\theta}\n",
    "$$\n",
    "\n",
    "So, finding a root for $T(s)$, is possible to verify when the displacement reaches multiplies of $2\\pi$ and later, obtain the period $P$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def solve_IVP(w,t_max=10,show_plot=False,delta=1):\n",
    "    #using Explicit Euler Method to estimate de IVP\n",
    "    t = np.linspace(0,t_max,50)\n",
    "    tt = np.linspace(0,t_max,1000)\n",
    "    h = t[1] - t[0]\n",
    "    theta = np.zeros(t.shape)\n",
    "    theta[0] = 0\n",
    "    for i in range(t.shape[0]-1):\n",
    "        theta[i+1] = theta[i] + h*w(theta[i],delta)\n",
    "        \n",
    "    #compute cubic splines for estimate the funtion\n",
    "    spl = cubic_spline(t,theta,end=\"Natural\")\n",
    "    \n",
    "    if show_plot:\n",
    "        plt.figure(figsize=(8,5))\n",
    "        #plot points of the numerical solution for theta\n",
    "        plt.plot(t,theta,\"ro\",label=r\"$\\Theta(t)$\")    \n",
    "        #plot spline function using the previous points\n",
    "        plt.plot(tt,cubic_spline_eval2(tt,spl),\"-g\",label=r\"Splines for $\\Theta(t)$\")\n",
    "        \n",
    "        #get the period of the function\n",
    "        #th = [2*np.pi,4*np.pi,6*np.pi]\n",
    "        th = [2*np.pi]\n",
    "        list_s = []   \n",
    "\n",
    "        for i in th:\n",
    "            f = lambda x: cubic_spline_eval(x,spl) - i\n",
    "            #obtain the time s where theta(s) = 2pi, 4pi and 6pi\n",
    "            s = bisect(f,0,t_max,verb=False)\n",
    "            #plot a vertical line\n",
    "            plt.axvline(x=s, color='k', linestyle='--')\n",
    "            list_s.append(s)\n",
    "\n",
    "        #plot the points where theta(s) = 2pi, 4pi and 6pi\n",
    "        plt.plot(list_s,th,\"Db\")   \n",
    "        \n",
    "        #plot the spline functions mod 2pi to verify the period\n",
    "        plt.plot(tt,(cubic_spline_eval2(tt,spl)%(2*np.pi)),\"-b\",label=r\"$\\Theta(t)$ mod $2\\pi$\")\n",
    "        #other configurations of the plot\n",
    "        plt.grid(True)\n",
    "        plt.legend(loc=\"best\")\n",
    "        plt.xlabel(r\"$t$\")\n",
    "        plt.ylabel(r\"$\\Theta(t)$\")\n",
    "        plt.title(\"Period P = \"+str(list_s[0])+\" [s]\")\n",
    "        #save image in pdf\n",
    "        plt.savefig(\"Theta.pdf\")\n",
    "        plt.show()\n",
    "        \n",
    "    Theta = lambda x: cubic_spline_eval(x,spl)\n",
    "    \n",
    "    return Theta"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In order to solve the first problem, a change of variable to Lagrangian coordinates $s$ and $\\tau$ will be done. Let:\n",
    "\n",
    "$$\n",
    "\\theta(s,\\tau) = \\Theta(\\tau - s)\n",
    "$$\n",
    "\n",
    "Where $s$ is the **time delay** since $\\theta$ was zero. With that, the coordinate system consisting of one spatial variable $\\theta$ and temporal variables $\\tau$ and $T$, changes into temporal variables $\\tau$, $T$ and $s$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Also, asociated to this change of variable, exists an inverse that let to obtain $s$ knowing the time $\\tau$ and the displacement $\\theta$:\n",
    "\n",
    "$$\n",
    "s = S(\\theta,\\tau)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using the chain rule to recast the differential equation, the following equalities are valid:\n",
    "    \n",
    "\\begin{align*}\n",
    "\\left(\\frac{\\partial}{\\partial \\tau} \\right)_s &= \\left(\\frac{\\partial \\tau}{\\partial \\tau} \\right)_s \\left(\\frac{\\partial }{\\partial \\tau} \\right)_\\theta + \\left(\\frac{\\partial \\theta}{\\partial \\tau} \\right)_s \\left(\\frac{\\partial }{\\partial \\theta} \\right)_\\tau = \\left( \\frac{\\partial}{\\partial \\tau}  \\right)_\\theta + \\omega \\left( \\frac{\\partial}{\\partial \\theta} \\right)_{\\tau}\\\\\n",
    "\\left( \\frac{\\partial}{\\partial \\theta} \\right)_\\tau &= \\left(\\frac{\\partial \\tau}{\\partial \\theta} \\right)_\\tau \\left(\\frac{\\partial }{\\partial \\tau} \\right)_s + \\left(\\frac{\\partial s}{\\partial \\theta} \\right)_\\tau \\left(\\frac{\\partial }{\\partial s} \\right)_\\tau = -\\frac{1}{\\omega} \\left( \\frac{\\partial}{\\partial s} \\right)_{\\tau}\n",
    "\\end{align*}\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Replacing into the equation, the problem is reduced to:\n",
    "\n",
    "$$\n",
    "\\left(\\frac{\\partial}{\\partial \\tau} \\right)_s (\\omega\\,u_0) = 0\n",
    "$$\n",
    "\n",
    "Which has constant solution for variable $\\tau$ in this coordinate system. The value of $u_0$ can be determinated defining a function $A_0$, which depends of $s$ and $T$ and the speed $\\omega$:\n",
    "\n",
    "$$\n",
    "u_0(\\theta, \\tau, T) = \\frac{A_0(s,T)}{\\omega(\\theta)}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This shows that an observed moving at a speed $\\omega(\\theta)$ will see the quantity $\\omega\\,u_0$ remain constant. But, this constant is allowed to modify its value in the slow scale $T$. The initial value of this _constant_ is given by the initial condition of the problem:\n",
    "\n",
    "$$\n",
    "A_0(s,0) = F(\\Theta(-s))\\omega(\\Theta(-s))\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In order to deteminate $A_0$ the second problem must be aborded."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Second Problem**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\n",
    "\\mathcal{O}(\\varepsilon):\n",
    "$\n",
    "\n",
    "$$\n",
    "\\frac{\\partial u_0}{\\partial T}+ \\frac{\\partial u_1}{\\partial \\tau} + \\frac{\\partial}{\\partial \\theta} (\\omega(\\theta) \\, u_1) = \\frac{\\partial^2 u_0}{\\partial \\theta^2}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Replacing the partially determined $u_0$, the problem becomes:\n",
    "\n",
    "$$\n",
    "\\frac{1}{\\omega} \\frac{\\partial A_0}{\\partial T} + \\frac{1}{\\omega} \\left(\\frac{\\partial}{\\partial \\tau}\\right)_s (\\omega\\,u_1) = \\frac{1}{\\omega} \\left( \\frac{\\partial}{\\partial s} \\right)_{\\tau}\\frac{1}{\\omega} \\left( \\frac{\\partial}{\\partial s} \\right)_{\\tau} \\frac{A_0}{\\omega}\n",
    "$$\n",
    "\n",
    "with $\\omega = \\omega(\\Theta(\\tau-s))$ a $P$-periodic function of $s$ and $\\tau$. Thus:\n",
    "\n",
    "$$\n",
    "\\frac{\\partial A_0}{\\partial T} + \\left(\\frac{\\partial}{\\partial \\tau}\\right)_s (\\omega\\,u_1) = \\frac{1}{\\omega^2} \\frac{\\partial^2A_0}{\\partial s^2}+ \\frac{3}{2} \\frac{\\partial}{\\partial s}\\left( \\frac{1}{\\omega^2} \\right) \\frac{\\partial A_0}{\\partial s}+\\frac{1}{2}\\frac{\\partial^2}{\\partial s^2} \\left( \\frac{1}{\\omega^2}\\right) A_0\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The right side of the equation has non zero average. In order to mantain the asymptoticness of the expansion of $f$ in $T = \\mathcal{O}(1)$, the quantity $\\omega \\,u_1$ has to be bounded as $\\tau$ increases. Thus, the average with respect to $\\tau$ of the right hand side must be removed by setting equal to $\\displaystyle\\frac{\\partial A_0}{\\partial T}$. The second and third terms of the right side has zero average because they are derivatives with respect to $s$ of functions of $\\omega$ and so derivatives with respect to $\\tau$. The next equalty must be satisfied:\n",
    "\n",
    "$$\n",
    "\\frac{\\partial A_0}{\\partial T} = \\frac{1}{P} \\int_{0}^{P} \\frac{1}{\\omega^2(\\Theta(\\tau-s))} \\mathrm{d}\\tau \\frac{\\partial^2A_0}{\\partial s^2}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The previous equation corresponds to a Parabolic PDE with diffusion coeficient equals to the average value of $1/\\omega^2$, that can be solved numerically in order to detemine $A_0$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Generically**, the problem can be solve using a numerical aproximation for the integral and using the value obtained in the PDE:\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "\\frac{1}{P}\\int_{0}^{P} \\frac{1}{\\omega^2(\\Theta(\\tau-s))} \\mathrm{d}\\tau \\approx \\frac{1}{P} \\sum_{i=1}^{K} w_i \\frac{1}{\\omega^2(\\Theta(\\tau_i-s))} = c(s)\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "Where $w_i$ is the weigh used in the sum of $K$ terms and depends of the numerical method used to estimate the value of the integral. Also, this value depends of temporal variable $s$. So, let $c(s)$ the value of the integral divided by the period $P$ using a numerical method.\n",
    "\n",
    "Now, considering the PDE with the estimation of the integral as the following equation:\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "\\frac{\\partial A_0}{\\partial T} &= c(s) \\frac{\\partial^2A_0}{\\partial s^2}\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "Using an implicit finite difference approach (that is inconditionally stable in this cases) for the variable $T$ and central difference for the variable $s$:\n",
    "\n",
    "The equations to be solved for $j = 1, 2, \\ldots N$ are:\n",
    "\n",
    "$$\n",
    "- \\sigma_j\\,a_{j-1,n+1} + (1+2\\,\\sigma_j)\\,a_{j,n+1}-\\sigma_j\\,a_{j+1,n} = a_{j,n}\n",
    "$$\n",
    "\n",
    "where $\\sigma_j = \\displaystyle\\frac{c(s_j)\\Delta T}{\\Delta s^2}$. And $s$ is a periodic variable.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The next code implements the solution of this PDE, according to the previous explanation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "def solve_a0(w,F,eps=1,t_max=1,points_s=50,points_T=80,delta=1,num=1):\n",
    "    Theta = solve_IVP(w,25,False,delta)\n",
    "    th1 = lambda x: Theta(x) - 2*np.pi\n",
    "    #t1 is the period P of Theta(t)\n",
    "    t1 = bisect(th1,0,25,verb=False)\n",
    "    #t2 is two periods of Theta(t)\n",
    "    t2 = 2*t1\n",
    "\n",
    "    #take equidistant values in the s axis\n",
    "    s = np.linspace(0,t1,points_s)\n",
    "    delta_s = t1/(points_s-1)\n",
    "    \n",
    "    #tau is taked between P and 2*P for correctly evaluation of Theta(t)\n",
    "    tau = np.linspace(t1,t2,points_s)\n",
    "    #vector c for save the values of c_j\n",
    "    c = np.zeros_like(s)\n",
    "    \n",
    "    for j in range(s.shape[0]):\n",
    "        #compute the numerical quadrature of myfun with variable s\n",
    "        myfun = lambda x: 1/((w(Theta(x-s[j]),delta))**2)\n",
    "        #using simpsons rule, number points_s must be EVEN\n",
    "        val = simpsons(myfun, points_s, t1, t2)/t1\n",
    "        #saving the value of c_j\n",
    "        c[j] = val\n",
    "    \n",
    "    \n",
    "    #computing delta_T\n",
    "    delta_T = (t_max)/(points_T-1)\n",
    "    #compute the values of sigma_j\n",
    "    sigma = (delta_T/(delta_s**2))*c\n",
    "    \n",
    "    \n",
    "    #compute matrix\n",
    "    A = np.zeros((points_s,points_s))\n",
    "    for i in range(points_s):\n",
    "        A[i][i] = 1+2*sigma[i]\n",
    "        if i < points_s-1:\n",
    "            A[i+1][i] = -sigma[i]\n",
    "        if i > 0:\n",
    "            A[i-1][i] = -sigma[i]\n",
    "    A[-1][0] = -sigma[0]\n",
    "    A[0][-1] = -sigma[-1]\n",
    "    \n",
    "    \n",
    "    sols = np.zeros((points_s,points_T+1))\n",
    "    for i in range(points_s):\n",
    "        #sols[i,0] = F(Theta(t1-s[i]),num)*w(Theta(t1-s[i]),delta)\n",
    "        sols[i,0] = F(Theta(-s[i]%t1),num)*w(Theta(-s[i]%t1),delta)\n",
    "        \n",
    "    for i in range(1,points_T+1):\n",
    "        sols[:,i] = np.linalg.solve(A,sols[:,i-1])    \n",
    "    return s,sols,delta_T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def print_simulation_a0(w,F,i,eps,delta=1,num=1):\n",
    "    t_max = 1\n",
    "    s,a0,d_T = solve_a0(w,F,eps,t_max,50,501,delta,num)\n",
    "    plt.figure(figsize=(10,5))\n",
    "    plt.plot(s,a0[:,i],\"ro\")\n",
    "    plt.xlabel(\"$s$\")\n",
    "    plt.ylabel(\"$A_0(s,T)$\")\n",
    "    plt.title(r\"$A_0(s,T)$ for $T = \"+str(d_T*i)+\"$\")\n",
    "    plt.ylim(-4,4)\n",
    "    plt.grid(True)\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, in order to obtain $u_0(\\theta,\\tau,T)$, is necessary to determinate which $s$ corresponds to a value of $\\hat{\\theta}$ given. The value of $u_0(\\hat{\\theta},\\tau,T)$ is given by:\n",
    "\n",
    "$$\n",
    "u_0(\\hat{\\theta},\\tau,T) = \\frac{A_0(s,T)}{\\omega(\\hat{\\theta})}\n",
    "$$\n",
    "\n",
    "Let $R(\\alpha)$ be the function:\n",
    "\n",
    "$$\n",
    "R(\\alpha) = \\Theta(\\alpha) - \\hat{\\theta}\n",
    "$$\n",
    "\n",
    "So, finding a root for $R(\\alpha)$, the value of $\\alpha$ such $\\Theta(\\alpha) = \\hat{\\theta}$ is found. In this case, is necesary go from lagrangian coordinate system to a traditional coordinate system space-time and $\\alpha$ is interpretated as $\\tau-s$ of the lagrangian coordinates. So, $s$ is creared from the original definition of lagrangian coordinates:\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "\\Theta(\\tau-s) &= \\Theta(\\alpha) = \\hat{\\theta}\\\\\n",
    "& \\Rightarrow \\tau-s = \\alpha\\\\\n",
    "& \\Rightarrow s = \\tau - \\alpha\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "Finally, with this value of $s$ is possible to come back to traditional coordinate system to interpretate and compare results."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The next code implements the calculation of values for $u_0$ taking equdistant points in the $\\theta$ axis from $0$ to $2\\pi$. The bisection method [1] was used for obtain the zeros of $R(\\alpha)$ when was needed. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def get_u0(w,F,i,eps=1,t_max=10,points_th=50,points_T=500,delta=1,num=1):\n",
    "    Theta = solve_IVP(w,25,False,delta)\n",
    "    th1 = lambda x: Theta(x) - 2*np.pi\n",
    "    #get period of the function\n",
    "    period = bisect(th1,0,25,verb=False)\n",
    "    \n",
    "    #points to using in the s axis\n",
    "    points_s = 50\n",
    "    \n",
    "    #computing A0(s,T)\n",
    "    s, a0, delta_T = solve_a0(w,F,eps,t_max,points_s,points_T,delta,num)\n",
    "    \n",
    "    #obtain points between 0 and 2pi\n",
    "    th = np.linspace(0,2*np.pi,points_th)\n",
    "    \n",
    "    #generate structure to save values of u0\n",
    "    u_0 = np.zeros_like(th)    \n",
    "    \n",
    "    #compute cubic spline for a0 in T = i*delta_T\n",
    "    spl = cubic_spline(s,a0[:,i],end=\"Natural\")\n",
    "    a_0 = lambda x: cubic_spline_eval(x,spl)\n",
    "\n",
    "    for j in range(1,th.shape[0]):\n",
    "        fun = lambda x: Theta(x)-th[j]\n",
    "        alpha = bisect(fun,0,2*period,verb=False)\n",
    "        s2 = ((i*delta_T/eps)-alpha)%period\n",
    "        u_0[j] = a_0(s2)/w(th[j],delta)\n",
    "    u_0[0] = u_0[-1]\n",
    "    return th,u_0,delta_T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def print_simulation_u0(w,F,i,eps,delta=1,num=1):\n",
    "    t_max = 1\n",
    "    th,u0,d_T = get_u0(w,F,i,eps,t_max,50,501,delta,num)\n",
    "    plt.figure(figsize=(10,5))\n",
    "    plt.plot(th,u0,\"gD\")\n",
    "    plt.xlabel(r\"$\\theta$\")\n",
    "    plt.ylabel(r\"$u_0(\\theta,\\tau,T)$\")\n",
    "    plt.title(r\"$u_0(\\theta,\\tau,T)$ for $T = \"+str(d_T*i)+\"$, $t=\"+str(d_T*i/eps)+\"$\")\n",
    "    plt.ylim(-4,4)\n",
    "    plt.grid(True)\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Using Finite Difference Method"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, the problem is solved numerically with finite differences with an implicit method. Consider the problem in the form:\n",
    "\n",
    "$$\n",
    "\\frac{\\partial u(\\theta,t)}{\\partial t} + \\frac{\\mathrm{d}\\omega(\\theta)}{\\mathrm{d} \\theta} u(\\theta,t) + \\frac{\\partial u(\\theta,t)}{\\partial \\theta} \\omega(\\theta) = \\varepsilon \\frac{\\partial^2 u(\\theta,t)}{\\partial \\theta^2}\n",
    "$$\n",
    "\n",
    "The formulas of the **implicit** finite differences will be **backward in time** and **central in space** for the first derivative (in order to mantain the spatial derivative with a second order aproximation). Let $u_{i,n} \\approx u(\\theta_i,t_n)$, where $\\theta_i$ and $t_n$ are $N$ and $M$ equidistant points in the spatial and temporal axes, respectively.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{aligned}\n",
    "\\frac{\\partial u(\\theta,t)}{\\partial t} &= \\frac{u_{i,n}-u_{i,n-1}}{\\Delta t} + \\mathcal{O}(\\Delta t)\\\\\n",
    "\\frac{\\partial u(\\theta,t)}{\\partial \\theta} &= \\frac{u_{i+1,n}-u_{i-1,n}}{2\\,\\Delta \\theta} + \\mathcal{O}(\\Delta \\theta^2)\\\\\n",
    "\\frac{\\partial^2 u(\\theta,t)}{\\partial \\theta^2} &= \\frac{u_{i-1,n}-2\\,u_{i,n}+u_{i+1,n}}{\\Delta \\theta^2} + \\mathcal{O}(\\Delta \\theta^2)\\\\\n",
    "\\end{aligned}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "with $\\Delta \\theta = \\displaystyle\\frac{2\\pi}{N-1}$ and $\\Delta t = \\displaystyle\\frac{t_{max}}{M-1}$ ."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "Replacing into the PDE, there are a system of equation to be solved for unkowns in different values of time where the equations are given by:\n",
    "\n",
    "$$\n",
    "\\left( \\frac{\\omega(\\theta_i)\\,\\Delta t}{2\\Delta \\theta}-\\frac{\\varepsilon\\, \\Delta t}{\\Delta \\theta^2} \\right)\\,u_{i-1,n+1}+ \\left( 1+\\frac{\\mathrm{d}\\omega(\\theta_i)}{\\mathrm{d} \\theta}\\Delta t + 2\\frac{\\varepsilon\\, \\Delta t}{\\Delta \\theta^2} \\right)\\,u_{i,n+1} + \\left(- \\frac{\\omega(\\theta_i)\\,\\Delta t}{2\\Delta \\theta}-\\frac{\\varepsilon\\, \\Delta t}{\\Delta \\theta^2} \\right)\\, u_{i+1,n+1} = u_{i,n}\n",
    "$$\n",
    "\n",
    "for $0 \\leq i<N$ and $n \\in [0,M]$. For a better visualization, let be:\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "\\alpha &= \\frac{\\varepsilon\\,\\Delta t}{\\Delta \\theta^2}\\\\\n",
    "\\beta_i &= \\frac{\\omega(\\theta_i)\\,\\Delta t}{\\Delta \\theta}\\\\\n",
    "\\gamma_i &= \\omega'(\\theta_i)\\,\\Delta t\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "So, the previous equation can be written as:\n",
    "\n",
    "$$\n",
    "\\left( \\frac{1}{2}-\\frac{\\varepsilon\\, \\Delta t}{\\Delta \\theta^2} \\right)\\,u_{i-1,n+1}+ \\left( 1+\\frac{\\mathrm{d}\\omega(\\theta_i)}{\\mathrm{d} \\theta}\\Delta t + 2\\frac{\\varepsilon\\, \\Delta t}{\\Delta \\theta^2} \\right)\\,u_{i,n+1} + \\left(- \\frac{\\omega(\\theta_i)\\,\\Delta t}{2\\Delta \\theta}-\\frac{\\varepsilon\\, \\Delta t}{\\Delta \\theta^2} \\right)\\, u_{i+1,n+1} = u_{i,n}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So, the previous equation can be written as:\n",
    "\n",
    "$$\n",
    "\\left( \\frac{\\beta_i}{2}- \\alpha \\right)\\,u_{i-1,n+1}+ \\left( 1+ \\gamma_i + 2 \\alpha \\right)\\,u_{i,n+1} + \\left(-  \\frac{\\beta_i}{2} - \\alpha \\right)\\, u_{i+1,n+1} = u_{i,n}\n",
    "$$\n",
    "\n",
    "for $0 \\leq i<N$ and $n \\in [0,M]$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "An **explicit** finite difference can be employed using the forward difference for the temporal variable and mantaining the previous differences for the spatial variable:\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "\\frac{\\partial u(\\theta,t)}{\\partial t} &= \\frac{u_{i,n+1}-u_{i,n}}{\\Delta t} + \\mathcal{O}(\\Delta t)\\\\\n",
    "\\frac{\\partial u(\\theta,t)}{\\partial \\theta} &= \\frac{u_{i+1,n}-u_{i-1,n}}{2\\,\\Delta \\theta} + \\mathcal{O}(\\Delta \\theta^2)\\\\\n",
    "\\frac{\\partial^2 u(\\theta,t)}{\\partial \\theta^2} &= \\frac{u_{i-1,n}-2\\,u_{i,n}+u_{i+1,n}}{\\Delta \\theta^2} + \\mathcal{O}(\\Delta \\theta^2)\\\\\n",
    "\\end{aligned}\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So, the iterative step to be calculated will be:\n",
    "$$\n",
    "u_{i,n+1} = \\left( \\frac{\\beta_i}{2}+ \\alpha \\right)\\,u_{i-1,n}+ \\left( 1- \\gamma_i - 2 \\alpha \\right)\\,u_{i,n} + \\left(-  \\frac{\\beta_i}{2} + \\alpha \\right)\\, u_{i+1,n}\n",
    "$$\n",
    "\n",
    "for $0 \\leq i<N$ and $n \\in [0,M]$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following code implements finite difference method for this problem."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def fd_method(w,wp,F,N,M,t_max,eps,delta=1,num=1,method=\"Implicit\"):\n",
    "    Theta = np.linspace(0,2*np.pi,N)\n",
    "    dt = t_max/(M-1)\n",
    "    \n",
    "    dth = Theta[1]-Theta[0]\n",
    "\n",
    "    A = np.zeros((N-1,N-1))\n",
    "    b = np.zeros((N,M))\n",
    "\n",
    "    b[:-1,0] = F(Theta[:-1],num)\n",
    "    b[-1] = b[0]\n",
    "    \n",
    "    #defining some \"constants\" for matrix generation\n",
    "    alpha = ((eps*dt)/(dth**2))*np.ones(N)\n",
    "    beta = (dt/dth)*w(Theta,delta)\n",
    "    gamma = dt*wp(Theta)\n",
    "    #print(\"eps\",eps)\n",
    "    #print(\"dt\",dt)\n",
    "    #print(\"dth**2\",dth**2)\n",
    "    #print(\"eps*dt/dth**2 = alpha\",alpha)\n",
    "    \n",
    "    if method==\"Implicit\":\n",
    "        for i in range(N-1):\n",
    "            A[i,i] = 1+2*alpha[i]+gamma[i]\n",
    "            A[i,i-1] = -alpha[i]-0.5*beta[i]\n",
    "            if i<N-2:\n",
    "                A[i,i+1] = -alpha[i]+0.5*beta[i]\n",
    "            else:\n",
    "                A[i,0] = -alpha[i]+0.5*beta[i]\n",
    "        for j in range(M-1):\n",
    "            b[:-1,j+1] = np.linalg.solve(A,b[:-1,j])\n",
    "            b[-1,:] = b[0,:]\n",
    "    elif method==\"Explicit\":\n",
    "        for i in range(N-1):\n",
    "            A[i,i] = 1-gamma[i]-2*alpha[i]\n",
    "            A[i,i-1] = 0.5*beta[i]+alpha[i]\n",
    "\n",
    "            if i<N-2:\n",
    "                A[i,i+1] = -0.5*beta[i]+alpha[i]\n",
    "            else:\n",
    "                A[i,0] = -0.5*beta[i]+alpha[i]\n",
    "        for j in range(M-1):\n",
    "            b[:-1,j+1] = np.dot(A,b[:-1,j])\n",
    "            b[-1,:] = b[0,:]\n",
    "        \n",
    "    return b"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Using FFT"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Fast Fourier Transform will be used to determine the spatial derivative of $u(\\theta,t)$ knowing the value of $u$ in the initial time, i.e. since the value of $u(\\theta,t=0)$ is known, the value of $\\displaystyle\\frac{\\partial^n u(\\theta, t=0)}{\\partial \\theta^n}$ can be determinated. Besides, the spatial variable is periodic with period equal to $2\\pi$. Later, using a forward approach for the temporal variable, solutions of $u(\\theta,t)$ will be obtained."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Consider the problem in the form:\n",
    "\n",
    "$$\n",
    "\\frac{\\partial u(\\theta,t)}{\\partial t} + \\frac{\\mathrm{d} \\omega(\\theta)}{\\mathrm{d} \\theta} u(\\theta,t) + \\frac{\\partial u(\\theta,t)}{\\partial \\theta} \\omega(\\theta) = \\varepsilon \\frac{\\partial^2 u(\\theta,t)}{\\partial \\theta^2}\n",
    "$$\n",
    "\n",
    "Let be $w^{(k)}_{j,n}$ the value obtained using FFT for the derivative of $k$ order for $u(\\theta,t)$ evaluated in $\\theta_j$ for the time $t_n$. Let be $v_{j,n}$ the aproximation of $u(\\theta,t)$ evaluated in $\\theta_j$ for the time $t_n$. Put all the previous definitions together, the problem is reduced to an iterative method that calculates the next time step using the information of the values of $u(\\theta_j,t_n)$, the first derivative of $u(\\theta_j,t_n)$ and the second derivative of $u(\\theta_j,t_n)$ for the actual time. Consider that $\\theta_i$ and $t_n$ are $N$ and $M$ equidistant points in the spatial and temporal axes.\n",
    "\n",
    "The iterative step to be used in this case is:\n",
    "\n",
    "$$\n",
    "v_{j,n+1} = \\varepsilon\\,\\Delta t\\,w_{j,n}^{(2)} - \\omega(\\theta_j)\\,\\Delta t\\,w_{j,n}^{(1)}+(1-\\omega'(\\theta_j)\\,\\Delta t)v_{j,n}\n",
    "$$\n",
    "for $0 \\leq j<N$ and $n \\in [0,M]$.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The next code implements the spectral derivative using FFT in order to obtain the derivative of $\\nu$ order of a function $v$ valid for the spatial interval $[0,2\\pi]$. The code was extracted from [3]."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def spectralDerivativeByFFT(v,nu=1):\n",
    "    if not np.all(np.isreal(v)):\n",
    "        raise ValueError('The input vector must be real')\n",
    "    N=v.shape[0]\n",
    "    K=np.fft.fftfreq(N)*N\n",
    "    iK=(1j*K)**nu\n",
    "    v_hat=np.fft.fft(v)\n",
    "    w_hat=iK*v_hat\n",
    "    if np.mod(nu,2)!=0:\n",
    "        w_hat[int(N/2)]=0\n",
    "    return np.real(np.fft.ifft(w_hat))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The next code implements the iterative step in order to solve this kind of problem."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def fft_method(w,wp,F,N,M,t_max,eps,delta=1,num=1):\n",
    "    Theta = np.linspace(0,2*np.pi,N)\n",
    "    #T = np.linspace(0,2,M)\n",
    "    #dt = T[1]-T[0]\n",
    "    dt = t_max/(M-1)\n",
    "        \n",
    "    #defining a matrix to save all the values of the function in all the time steps\n",
    "    v = np.zeros((N,M))\n",
    "    #setting values for the initial time\n",
    "    v[:,0] = F(Theta,num)\n",
    "    \n",
    "    w1 = np.zeros(N)\n",
    "    w2 = np.zeros(N)\n",
    "    \n",
    "    #iterative step\n",
    "    for n in range(M-1):\n",
    "        #calculate the first derivative with FFT\n",
    "        w1[1:] = spectralDerivativeByFFT(v[1:,n],1)\n",
    "        w1[0] = w1[-1]\n",
    "        #calculate the second derivative with FFT\n",
    "        w2[1:] = spectralDerivativeByFFT(v[1:,n],2)\n",
    "        w2[0] = w2[-1]\n",
    "              \n",
    "        #get the values of the function in the next time\n",
    "        v[:,n+1] = eps*dt*w2 - dt*w(Theta,delta)*w1+(np.ones(N)-dt*wp(Theta))*v[:,n]\n",
    "    return v"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Numerical Examples"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For Numerical Examples, three different functions associated to the initial values will be used. These are:\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "F_1(\\theta) &= 2\\,\\sin{(\\theta)}\\\\\n",
    "F_2(\\theta) &= 2\\,\\exp{(-(\\theta-\\pi)^2)}\\\\\n",
    "F_3(\\theta) &= 2\\,\\sin^2{(0.5\\,\\theta)}\\\\\n",
    "\\end{aligned}\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def F(theta,num=1):\n",
    "    if num == 1:\n",
    "        return 2*np.sin(theta)\n",
    "    elif num==2:\n",
    "        return 2*np.exp(-(theta-np.pi)**2)\n",
    "    elif num==3:\n",
    "        return 2*np.sin(0.5*theta)**2\n",
    "    return theta*0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "According to all the theoretical part explained before, two numerical examples will be tested for the diffusion and advection equation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Again, consider the following problem that shows advection around a periodic domain $0 < \\theta < 2\\pi$ with some small diffusion. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\\begin{align*}\n",
    "\\frac{\\partial u(\\theta,t)}{\\partial t} + \\frac{\\partial}{\\partial \\theta} (\\omega(\\theta)\\,u(\\theta,t)) &= \\varepsilon\\frac{\\partial^2u(\\theta,t)}{\\partial \\theta^2}\\\\\n",
    "u(\\theta,0) &= F(\\theta)\\\\\n",
    "u(0,t) &= u(2\\pi,t)\n",
    "\\end{align*}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "where $\\omega(\\theta)$ is given over $0 \\leq \\theta \\leq 2\\pi$, periodic and positive. $F(\\theta)$ is a known function and $\\varepsilon$ is a tiny value."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Numerical Example 1**\n",
    "The first numerical example considers a constant value $\\delta$ for the speed $\\omega$. Let be:\n",
    "\\begin{align}\n",
    "\\omega(\\theta) &= \\delta\\\\\n",
    "\\omega'(\\theta) &= 0\\\\\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def w(theta, delta=1):\n",
    "    return delta*(theta**0) \n",
    "def wp(theta):\n",
    "    return theta*0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def plot_functions(w,wp,F,delta=1,num=1):\n",
    "    th = np.linspace(0,4*np.pi,200)\n",
    "    thh = np.linspace(0,2*np.pi,100)\n",
    "        \n",
    "    plt.figure(figsize=(18,5))\n",
    "\n",
    "    \n",
    "    label_fun = \"0\"\n",
    "    if num==1:\n",
    "        label_fun=r\"2\\,\\sin{(\\theta)}\"\n",
    "    elif num==2:\n",
    "        label_fun=r\"2\\,\\exp(-(\\theta-\\pi)^2)\"\n",
    "    elif num==3:\n",
    "        label_fun=r\"2\\,\\sin(0.5\\,\\theta)^2\"\n",
    "    \n",
    "    plt.subplot(121)\n",
    "    plt.plot(th,w(th,delta),\"-r\",label=r\"$w(\\theta) = \"+str(delta)+\"$\")\n",
    "    plt.plot(th,wp(th),\"-b\",label=r\"$w'(\\theta) = 0$\")\n",
    "    plt.xlim(-1,12)\n",
    "    plt.ylim(-4,4)\n",
    "    plt.axvline(x=0, color='k', linestyle='--')\n",
    "    plt.axvline(x=2*np.pi, color='k', linestyle='--')\n",
    "    plt.grid(True)\n",
    "    plt.xlabel(r\"$\\theta$\")\n",
    "    plt.legend(loc=\"best\")\n",
    "    \n",
    "    plt.subplot(122)\n",
    "    plt.plot(thh,F(thh,num),\"-g\",label = r\"$F(\\theta) = \"+label_fun+\"$\")\n",
    "    plt.xlim(-1,7)\n",
    "    plt.ylim(-4,4)\n",
    "    plt.axvline(x=0, color='k', linestyle='--')\n",
    "    plt.axvline(x=2*np.pi, color='k', linestyle='--')\n",
    "    plt.grid(True)\n",
    "    plt.xlabel(r\"$\\theta$\")\n",
    "    plt.ylabel(r\"$F(\\theta)$\")\n",
    "    plt.legend(loc=\"best\")\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Widget Javascript not detected.  It may not be installed or enabled properly.\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "6248020fa4cf4c9b9401700aef86903d"
      }
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<function __main__.plot_functions>"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "delta_widget = FloatSlider(min=0.4,max=2.0, step=0.1, value=1.0)\n",
    "fun_widget = Dropdown(options=[1,2,3],\n",
    "                      value=1, description='Function F')\n",
    "interact(plot_functions,w=fixed(w),wp=fixed(wp),F=fixed(F),delta=delta_widget,num=fun_widget)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As can be seen previously, the function $\\omega(\\theta)$ is a positive and periodic function around the interval $[0,2\\pi]$. The derivative of a constant function is zero. The function $F(\\theta)$ is shown between the same interval. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first approach for this example consists in determinate numerically the solution of the partial differential equation using Finite Differences (Explicit and Implicit) and Fast Fourier Transform."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def simulation(w,wp,F,i,N,M,t_max,eps,fd_mode,delta,num): \n",
    "    \n",
    "    #spatial axis\n",
    "    Theta = np.linspace(0,2*np.pi,N)\n",
    "    #temporal axis\n",
    "    T = np.linspace(0,t_max,M)\n",
    "    \n",
    "    y1 = fft_method(w,wp,F,N,M,t_max,eps,delta,num)\n",
    "    y2 = np.zeros_like(y1)\n",
    "    label_fd = \"No mode selected for FD\"\n",
    "    if fd_mode == \"Implicit\":\n",
    "        y2 = fd_method(w,wp,F,N,M,t_max,eps,delta,num)\n",
    "        label_fd = \"Implicit Finite Difference\"\n",
    "    elif fd_mode == \"Explicit\":\n",
    "        y2 = fd_method(w,wp,F,N,M,t_max,eps,delta,num,\"Explicit\")\n",
    "        label_fd = \"Explicit Finite Difference\"\n",
    "        \n",
    "    a1 = np.zeros(i+1)\n",
    "    a2 = np.zeros(i+1)\n",
    "    \n",
    "    for k in range(i+1):\n",
    "        a1[k] = max(abs(y1[:,k]))\n",
    "        a2[k] = max(abs(y2[:,k]))\n",
    "          \n",
    "    plt.figure(figsize=(17,13))\n",
    "\n",
    "    plt.subplot(2, 2, 1)\n",
    "    plt.title(\"Simulation at t = \"+str(T[i])+\" [s]\")\n",
    "    plt.plot(Theta[:],y1[:,i],\"-yo\", label=\"Fast Fourier Transform\")\n",
    "    plt.ylim(-5,6)\n",
    "    plt.grid(True)\n",
    "    plt.xlabel(r\"$\\theta$\")\n",
    "    plt.ylabel(r\"$u(\\theta,t)$\")\n",
    "    plt.legend(loc=\"best\")\n",
    "    \n",
    "    plt.subplot(2, 2, 2)\n",
    "    plt.title(\"Simulation at t = \"+str(T[i])+\" [s]\")\n",
    "    plt.plot(Theta[:],y2[:,i],\"-bo\", label=label_fd)\n",
    "    plt.ylim(-5,6)\n",
    "    plt.grid(True)\n",
    "    plt.xlabel(r\"$\\theta$\")\n",
    "    plt.ylabel(r\"$u(\\theta,t)$\")\n",
    "    plt.legend(loc=\"best\")\n",
    "    \n",
    "    plt.subplot(2, 2, 3)\n",
    "    plt.title(\"Max Absolute Value\")\n",
    "    plt.plot(T[:i+1],a1,\"-yo\", label=\"Fast Fourier Transform\")\n",
    "    plt.xlim(-0.1*t_max,1.1*t_max)\n",
    "    plt.grid(True)\n",
    "    plt.xlabel(r\"$t$\")\n",
    "    plt.ylabel(r\"$u(\\theta,t)$\")\n",
    "    plt.legend(loc=\"best\")\n",
    "    \n",
    "    plt.subplot(2, 2, 4)\n",
    "    plt.title(\"Max Absolute Value\")\n",
    "    plt.plot(T[:i+1],a2,\"-bo\", label=label_fd)\n",
    "    plt.xlim(-0.1*t_max,1.1*t_max)\n",
    "    plt.grid(True)\n",
    "    plt.xlabel(r\"$t$\")\n",
    "    plt.ylabel(r\"$u(\\theta,t)$\")\n",
    "    plt.legend(loc=\"best\")\n",
    "    \n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The numerical simulation can be observed on a interactive way changing some parameters:\n",
    "* $i$: solution of the equation at time $t_i = i\\,\\Delta t$. Initially $i = 0$.\n",
    "* $N$: number of points to be evaluated in the $\\theta$ axis. Initially $N = 50$.\n",
    "* $M$: number of intervals of time to be displayed. Initially $M = 80$.\n",
    "\n",
    "Aditionally, various values of $\\varepsilon$ will be used."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Widget Javascript not detected.  It may not be installed or enabled properly.\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "4836b26365a04acf833e9bd3138d1d8a"
      }
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<function __main__.simulation>"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "i_widget = IntSlider(min=0, max=79, value=0)\n",
    "N_widget = IntSlider(min=2, max=200, value=50)\n",
    "M_widget = IntSlider(min=2, max=501, value=80)\n",
    "eps_widget = Dropdown(options=[100.,10.,2.,1.,1e-1,1e-2,1e-3,1e-4,1e-5,1e-6,1e-7,\n",
    "                               1e-8,1e-9,1e-10,1e-11,1e-12,1e-13,1e-14,1e-15,1e-16,1e-17,0.],\n",
    "                      value=0, description='eps')\n",
    "delta_widget = FloatSlider(min=0.1,max=2.0, step=0.1, value=1.0)\n",
    "fun_widget = Dropdown(options=[1,2,3],\n",
    "                      value=1, description='Function F')\n",
    "fd_widget = Dropdown(options =  [\"Explicit\", \"Implicit\"], value=\"Implicit\", description = \"FD Mode\")\n",
    "time_widget = IntSlider(min=1,max=10,step=1,value=1,description=\"Max Time\")\n",
    "def update_i_values(*args):\n",
    "    i_widget.max = M_widget.value-1\n",
    "M_widget.observe(update_i_values, 'value')\n",
    "\n",
    "#def Wp(theta,delta=1):\n",
    "#    return 0*theta\n",
    "\n",
    "interact(simulation,\n",
    "         w=fixed(w),\n",
    "         wp=fixed(wp),\n",
    "         F=fixed(F),\n",
    "         i=i_widget,\n",
    "         N=N_widget,\n",
    "         M=M_widget,\n",
    "         eps=eps_widget,\n",
    "         t_max = time_widget,\n",
    "         fd_mode=fd_widget,\n",
    "         delta=delta_widget,\n",
    "         num=fun_widget)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first case $\\varepsilon = 0$ is simply an advection equation without diffusion. Since the speed $\\omega$ is constant the solution, tipically a wave is moving at this constant rate. Is expected that the maximum absolute value of the solution observed be the same in all the time elapsed.\n",
    "- The FFT has stability problems when $N$ is much greater than $M$ (apparently with a quotient $N/M$ greater to 1 and nearer to 2). Also, when the value of $M$ is tiny the simulation has stability problems more prematurely. In cases where there are no stability issues, the maximum absolute value of the solution is increasing slightly. When the value of $M$ is bigger, i.e., the time between two neighbor points is shorter, the increment of the maximun value obtained is lower that with a $M$ small. This factor maybe is associated with the numerical error of the method.\n",
    "- The Implicit FD does not have problems of stability. It can be observed that the maximum absolute value of the solutions is decresing in time slightly, and also, with bigger values of $M$ the change in this decrease is lower that the same effect with a small $M$. Again, this difference between the expected behavior and the numerical results can be interpretated as numerical error asociated.\n",
    "- The Explicit FD does not have problems of stability with in this case and its behaviour is better to comparing with FFT. The maximum value is increasing int time such as FFT. Also, with as greater is the value of $M$, the total increment of the maximum value in time $t=10$ is lower."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Cases when $\\varepsilon > 0$ are examples of a advection-diffusion equation when the solution is disipated around the axis. When the coeficient is greater, the disipation of the values is faster. The advection process, barely can be seen. The problem can be seen as the oposite as the original problem with a slow difusion, where now the advection is the slow process compared to the diffusion and if the coeficient of difusion is a big value, then the solution is practicaly a constant solution (explained for a really fast difusion, without considering what was the initial distribution of the function).\n",
    "- Again, the FFT has stability problems related to the quotient $N/M$ in the same way that the previous case. However, in this cases there is an aditional effect of the $\\varepsilon$ parameter. For $\\varepsilon = 1$, the problem of stability ocurrs with a $N/M$ quotient near to 0.10. It can be seen that at higher values of $\\varepsilon$, the problem of stability occurs sooner, so the quotient $N/M$ has to mantain as low as possible. When the stability problem arises, the maximum absolute value is disproportionately huge, and as long as the coeficient $N/M$ increases, the increment of the maximum absolute value ocurs sooner. When there are not stability problems, the behaviour of the maximum absolute value is the expected, observing a diffusion processs in the time and decresing this maximum value faster if the diffusion coefficient is bigger.\n",
    "- The Implicit FD does not have problems of stability. The behaviour of the maximun absolute value is the expected in all the cases and, as ocurrs in FFT, as bigger the diffusion coefficient is, more sooner is the time when the maximum aplitude reaches the zero value.\n",
    "- The Explicit FD has problems of stability with the quotient $N/M$ near to 0.16, but the stability issues ocurrs with minor frequency that the FFT. The behaviour in this cases, when the method is stable, is the expected behaviours for a not small diffusion coefficient."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The final cases consider $0 < \\varepsilon \\ll 1$ and is the interest for this work. Is expected to observe advection in time and apparently no diffusion process.\n",
    "\n",
    "- In the FFT the phenomenon of slow diffusion cannot be observed clearly, but the plot of maximum aplitude value shows some differences varying the paremeters. With $\\varepsilon  = 0.1$ and $\\varepsilon = 0.01$ the maximum absolute value of the solutions decrease slower as the order of magnitude of $\\varepsilon$ decrease. However, with orders of magnitude for $\\varepsilon$ smaller than $10^{-3}$ the phenomenon the maximum is increasing in time. As the order of magnitude is continue decreasing, the numerical error is reaching the results obtained for the first case.\n",
    "- In the Implicit FD during the decrease of orders of magnitude for $\\varepsilon$ it can be seen that the inclination of maximum absolute value to have lower values is the expected and for smaller values of $\\varepsilon$, the decreace is smaller too. Aproximately, for order of magnitud of $10^{-7}$ it cannot be seen differences between the maximum absolute value comparing to the first case explained before.\n",
    "- The Explicit FD has the same behaviour than FFT. When $\\varepsilon$ has a order of magnitude lower than $10^{-2}$ the decreasing maximum value in time of the solutions changes to an increasing maximum value in time. The solutions are approaching to the first case."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In general, the Fast Fourier Transform has some problems caused by certain parameter values and is a conditionally stable method. The Explicit Finite Difference is better than FFT, but has some issues on stability too. Otherwise, the Implicit Finite Difference seems to be more trustworthy, since the behaviour expected for this case is aproximately the same that the numerical results. However, all the methods introduce some numerical error, which could be seen in the first case when $\\varepsilon = 0$. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, the question is, ¿what is the effect of the small values for $\\varepsilon$ of the equation for the numerical results? First of all, using the method of multiscales, is necessary to known how is the advance of the solutions according to the speed $\\omega$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The next code shows the numerical solution of the Initial Value Problem asociated to the advance of $\\Theta$ in the time $t$.\n",
    "\n",
    "$$\n",
    "\\frac{\\mathrm{d}\\theta}{\\mathrm{d}t} = \\omega(\\theta)\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def plot_Theta(w,delta):\n",
    "    t_max = 25\n",
    "    Theta = solve_IVP(w,t_max,False,delta)\n",
    "    vTheta = np.vectorize(Theta)\n",
    "    t = np.linspace(0,t_max,50)\n",
    "    tt = np.linspace(0,t_max,1000)\n",
    "    plt.figure(figsize=(8,5))\n",
    "    #plot points of the numerical solution for theta\n",
    "    plt.plot(t,vTheta(t),\"ro\",label=r\"$\\Theta(t)$\")    \n",
    "    #plot spline function using the previous points\n",
    "    plt.plot(tt,vTheta(tt),\"-g\",label=r\"Splines for $\\Theta(t)$\")\n",
    "\n",
    "    #get the period of the function \n",
    "    f = lambda x: Theta(x) - 2*np.pi\n",
    "    #obtain the time s where theta(s) = 2pi\n",
    "    s = bisect(f,0,t_max,verb=False)\n",
    "    #plot a vertical line\n",
    "    plt.axvline(x=s, color='k', linestyle='--')\n",
    "\n",
    "    #plot the points where theta(s) = 2pi, 4pi and 6pi\n",
    "    plt.plot(s,2*np.pi,\"Db\")   \n",
    "\n",
    "    #plot the spline functions mod 2pi to verify the period\n",
    "    plt.plot(tt,vTheta(tt)%(2*np.pi),\"-b\",label=r\"$\\Theta(t)$ mod $2\\pi$\")\n",
    "    #other configurations of the plot\n",
    "    plt.grid(True)\n",
    "    plt.legend(loc=\"best\")\n",
    "    plt.xlabel(r\"$t$\")\n",
    "    plt.ylabel(r\"$\\Theta(t)$\")\n",
    "    plt.title(\"Period P = \"+str(s)+\" [s]\")\n",
    "    #save image in pdf\n",
    "    #plt.savefig(\"Theta.pdf\")\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Widget Javascript not detected.  It may not be installed or enabled properly.\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "6cc3270a9e444cb7b1caf628137b55f1"
      }
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<function __main__.plot_Theta>"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "delta_widget = FloatSlider(min=0.4,max=2.0, step=0.1, value=1.0)\n",
    "interact(plot_Theta,w=fixed(w),delta=delta_widget)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As expected, since the speed $\\omega$ in this case is $\\omega(\\theta) = 1$, a constant value. The position $\\Theta$ in time is the function $\\Theta(t) = t$, the identity function. If a different scalar value was used, the solution of the Initial Value Problem would be a multiple of the identity function by exactly this scalar value. The numerical method used in this case is appropiated since reflects the behavior of the solution. The period $P$ in this case corresponds to the value $P = 2\\pi \\approx 6.28312$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The next step is to solve the Partial Differential Equation asociated to $A_0(s,T)$:\n",
    "\n",
    "$$\n",
    "\\frac{\\partial A_0}{\\partial T}(s,T) = \\frac{1}{P} \\int_{0}^{P} \\frac{1}{\\omega^2(\\Theta(\\tau-s))} \\mathrm{d}\\tau \\frac{\\partial^2A_0}{\\partial s^2}(s,T)\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The next code solves the equation, mixing a numerical method for the integral and a finite difference method for the PDE."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Widget Javascript not detected.  It may not be installed or enabled properly.\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "3be88a68a9e44607a0d6d5eaf0ee95a0"
      }
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<function __main__.print_simulation_a0>"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "i_widget = IntSlider(min=0, max=499, value=0)\n",
    "eps_widget = Dropdown(options=[100.,10.,2.,1.,1e-1,1e-2,1e-3,1e-4,1e-5,1e-6,1e-7,\n",
    "                               1e-8,1e-9,1e-10,1e-11,1e-12,1e-13,1e-14,1e-15,1e-16,1e-17,0.],\n",
    "                      value=0, description='eps')\n",
    "delta_widget = FloatSlider(min=0.4,max=2.0, step=0.1, value=1.0)\n",
    "fun_widget = Dropdown(options=[1,2,3],\n",
    "                      value=1, description='Function F')\n",
    "interact(print_simulation_a0,w=fixed(w),F=fixed(F),i=i_widget,eps=eps_widget,delta=delta_widget,num=fun_widget)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As can be seen, the solution is the usual for a diffusion equation. The difference with traditional approaches is that the horizontal axis corresponds to a time variable that measures the delay to be applied to the function $\\Theta(t)$ such that $\\theta = 0$. Recalling that the $s$ that corresponds to a certain $\\hat{\\theta}$ is given by $s = \\tau - \\alpha$ with $\\alpha$ such that $\\Theta(\\alpha) = \\hat{\\theta}$, for some values of $\\tau$ the value of $s$ needed is negative. However, the variable is periodic on the interval $[0,P]$ and applying the modulo operator, the $s$ appears in the plot nearer of $P$ for the right side. For example, consider $T = \\tau = 0$ and the displacement $\\theta = \\pi/2$. The value of $\\alpha$ such that $\\Theta(\\alpha) = \\pi/2$ is $\\alpha = \\pi/2$ and the value of $s$ in this case is $s = 0 - \\pi/2 = -\\pi/2$. Using the periodicity of variable $s$ the value delivered to make next calculations is $s = (-\\pi/2)\\%2\\pi = 3\\pi/2$. The expected displacement of the advection process is realized because with a fixed value of displacement $\\theta$ in different moments of time, the value of $s$ obtained is different and moves around the values $[0,P]$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The next code shows the solution $u_0(\\theta,\\tau,T)$ obtained with the method of multiscales and expresed in function of the slow scale: The last part is really important, because there are not a direct implication of the $\\varepsilon$ parameter in this scale. With different values of $\\varepsilon$, the caused effect is that the diffusion ocurrs more faster or slower in the fast scale; that is the expected effect of the coefficient of diffusion."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Widget Javascript not detected.  It may not be installed or enabled properly.\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "5e196beaebab48ee8f9c8a93d8b74c27"
      }
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<function __main__.print_simulation_u0>"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "i_widget = IntSlider(min=0, max=500, value=0)\n",
    "eps_widget = Dropdown(options=[100.,10.,2.,1.,1e-1,1e-2,1e-3,1e-4,1e-5,1e-6,1e-7,\n",
    "                               1e-8,1e-9,1e-10,1e-11,1e-12,1e-13,1e-14,1e-15,1e-16,1e-17],\n",
    "                      value=1, description='eps')\n",
    "delta_widget = FloatSlider(min=0.4,max=2.0, step=0.1, value=1.0)\n",
    "fun_widget = Dropdown(options=[1,2,3],\n",
    "                      value=1, description='Function F')\n",
    "interact(print_simulation_u0,w=fixed(w),F=fixed(F),i=i_widget,eps=eps_widget,delta=delta_widget,num=fun_widget)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The final experiment for this numerical example considers the comparation of the numerical solution obtained with one pure method numerically and the solution obtained using the method of multiscales. For that, is necesary to use equivalent times for the pure methods to the results in the slow scale. This times will be changed each time that a different value of $\\varepsilon$ be selected.\n",
    "\n",
    "Consider a given $T_{max}$ in the slow time scale. For simplicity, $T_{max} = 1$ and $\\varepsilon$ will take the values $1 \\times 10^{\\,p}$ with negative values for $p$. Also, a few values of $\\varepsilon > 1$ will be taken into account, but it is not the main porpouse of this work. For example, let be $\\varepsilon = 10^{-5}$. Since $T_{max}= 1$ units of time in the slow scale, the equvalent value in the fast scale with this $\\varepsilon$ is $t_{max} = 1/\\varepsilon = 10^{5}$ units of time. With this huge number and a reasonable low value for time step a really large number of iterations (and memory if is necessary to store all the values) is necessary to do. But, with the multiscales is possible to realize few iterations and advance \"big\" steps of time without recuring to a large time step or a using too much memory. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The next code shows the comparison between the method of multiscales and a traditional approach in order to determinate the numerical values of the advection-diffusion problem where the diffusion occurs very slowly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def comparison(w,wp,F,i,eps,T_max,N,M,method=1,delta=1,num=1):\n",
    "    th,u0,d_T = get_u0(w,F,i,eps,T_max,N,M,delta,num)\n",
    "\n",
    "    t_max = T_max/eps\n",
    "    M2 = int(t_max/d_T + 1)\n",
    "       \n",
    "    label_method = \"\"\n",
    "    u1 = np.zeros((N,M))\n",
    "    if method == 1:\n",
    "        u1 = fd_method(w,wp,F,N,M2,t_max,eps,delta,num)\n",
    "        label_method = \"Implicit FD\"\n",
    "    elif method == 2:\n",
    "        u1 = fd_method(w,wp,F,N,M2,t_max,eps,delta,num,\"Explicit\")\n",
    "        label_method = \"Explicit FD\"\n",
    "    elif method == 3:\n",
    "        u1 = fft_method(w,wp,F,N,M2,t_max,eps,delta,num)\n",
    "        label_method = \"FFT\"\n",
    "\n",
    "    j = int(i/eps)\n",
    "    \n",
    "    print(\"Points Fast Time:\",M)\n",
    "    print(\"Points Slow Time:\",M2)\n",
    "    print(\"i\",i)\n",
    "    print(\"j\",j)\n",
    "    \n",
    "    plt.figure(figsize=(13,5))\n",
    "    \n",
    "    plt.subplot(121)\n",
    "    plt.plot(th,u0,\"rs\",ms=10,label=\"Multiscales\")\n",
    "    plt.plot(th,u1[:,j],\"gD\",ms=4,label=label_method)\n",
    "    plt.title(r\"Simulation at $T = \"+str(d_T*i)+\"$, $t=\"+str(d_T*i/eps)+\"$\")\n",
    "    plt.ylim(-4,4)\n",
    "    plt.grid()\n",
    "    plt.legend(loc=\"best\")\n",
    "    \n",
    "    plt.subplot(122)\n",
    "    plt.plot(th,u0,\"rs\",ms=10,label=\"Multiscales\")\n",
    "    plt.plot(th,u1[:,j],\"gD\",ms=4,label=label_method)\n",
    "    plt.title(r\"Simulation at $T = \"+str(d_T*i)+\"$, $t=\"+str(d_T*i/eps)+\"$\")\n",
    "    #plt.ylim(-4,4)\n",
    "    plt.grid()\n",
    "    plt.legend(loc=\"best\")\n",
    "    \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Widget Javascript not detected.  It may not be installed or enabled properly.\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "bfb47d64d5224a5aac71c7f22c64eed4"
      }
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<function __main__.comparison>"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAwQAAAFKCAYAAACjNpjXAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xm8HFWZ//HPkwWyQBYSuIEQiexDFEUCEgW5AWQLi4Eg\niCgBNOCMKxlXRAIybhiVERWCIDICkSUIYQ1LLiO/BBDEcRwIEVlkC0tIQi4kEMjz+6O6b/r2re6u\nXqqrqvv7fr36dbnVtZyuXOr0Oec5zzF3R0RERERE2lO/pAsgIiIiIiLJUYNARERERKSNqUEgIiIi\nItLG1CAQEREREWljahCIiIiIiLQxNQhERERERNqYGgQiIiIiIm1MDQIRERERkTamBkFMzKzTzNzM\nprdqOdLyGUVEsiotz1HVFSLtTQ2CKpnZtmY2x8yWmNkbZrbCzB41s9+a2eSky9doZvZ+M5tlZuOT\nLku9qv0sZraJmb2Tq8iivDaLsez9zOwrub+7tWb2jJnNNrOhzTxHhfOfZGZfasS5Il4v9fdE2pfq\niuxSXaG6Io5zpN2ApAuQJWY2EbgHWAdcDvwfMBjYATgQWA0szO3+37n31jW/pA31fuAsoAt4qui9\nrH3Gcp8lzADgxKJtnwM+BPw78GLB9jfd/dX6i1jST4EvAtcDs4F/yf2+m5kd4O7rm3SOcn4E3A+c\nX+d5osrCPZE2pLpCdQWqK8pRXZFG7q5XxBcwH3DgfSXeH5N0GUPK1Jkr8/Qaj5+eO74z6c/SgHtR\n92cBHgLWAAOaWO4JwHrguqLtX8h9nuObcY4K598+d55v657o1e4v1RXJf54674XqCtUVTbsnaXkl\nXoAsvYAlwCsR9+3zcC14yOwPfAd4OvfAuB/YK7fPvsC9wOvAC8CZReedlTvH+JBrPgV0RSjHpsC5\nueu+ArwJPA78ABgScq3i12Wlzp3bPhr4BfAM8Fbu5y+AUUX75e/HfgS9KP/IlWUpcGLE+9yQzxLx\nWgNz57+/yX935+bKuk/R9kG5v5NbmnGOMuf+Q4l7e2673hO92vuF6grVFaorws6tuiLFL4UMVecf\nwE5mdpS7z6vjPD8A+hMMl20EzAQWmNmngUuAOcAVwMeBc8zsSXf/XX1F72Us8BngOuBK4G2CyuVr\nwG7AQbn95gFbAjOA7wGP5rb/o9SJzWw4sIigF+BS4M+5c34O2M/M9nT31UWHfY9gOPkigofo54DL\nzOxxd/9/SX2WEBMI/r0errSjmfUDqokTfdVLDznuQdA78UDhRndfa2Z/yb1fSSPOUcocgr/nwwj+\n7bpz2xcX79jA+5L2eyLtTXVFQHVFBaorANUV6ZB0iyRLL2ASQS+GE/RMXErwR/0vIft2UrrX58/A\nRgXbj8htXwdMLNi+EUHPz+KCbbOov9dnI2BgyPHfze27Z0iZOyN+xv/IbfvXon3/Lbf9uyHnfrjo\nfowleNhfFeHfpCGfJeK//0m542dE2Hc84T0hpV59/j0LzvW/wIsl3rs6d/xGFcpT9zkqnP8G4KVm\n3Zcs3BO92veF6ooon1F1hauuiPu+ZOGepOWlEYIquPtiM9udoJfmEIL/6U8CMLM/Ejzsnohwql+5\n+1sFv/8x9/N+d3+w4HpvmdkDwIcb8gEKzpv/bzMbQDCU2h+4E/g28EGKWsJVmAq8TNATUOgigkla\nU4Ezi977ZWGZ3P05M1tKMAGvrJg/S7EP5H5W7PUBlgEfreLcy8q8N4Sg0guztmCft0rs06hzlPMB\nmntfsnBPpE2prohEdUVAdUU41RVNpgZBldz9fwl6DzCzbQiGHD8D7APcYGa7Fz3Aw/SqCNx9hZkB\nPBmy7wpgVJ3F7sPM/hU4jWBoszj97Mg6Tv1u4EF3f7two7u/nXtwfyDkmLCKcTmwTZQLxvhZiu1G\nMMz8v5V2dPe1BBVNI7wBbFHivUEF+8R9jlBmtjmwNVAxVKGB9yXV90REdUVFqitQXVGK6ormU4Og\nDu7+NHC5mf0XQc/Nh4E9CSZ6lfNOldt7XbbMe5H+Pc3sdIK0WQuA/wSeJ2jZjgUuo/nrU5T63Fbp\nwGZ9llw84/uAR3MPqkr79wc2r+ISL7t7qfvwPLCLmW3s7sW9FGMJJi9W+mLRiHOUkq+4/1xpxwbe\nl7TfE5EeqisaRnWF6opSVFfUSQ2CBnB3N7P7CR7yY2O+XD5/8WYU5Ec2s0EEE6Eej3COT+WOPcQL\nJuGY2cEh+5arVMI8QTCZbkBhz09uiHZHwnt46hHnZym0A7AJER5kOeMI78Ur5d2Uznf9J4Lc5Xuy\nIWQg/2/+foIc35U04hyl7Jb7GeXeNOq+pP2eiPShuqIX1RUB1RXhVFc0mRoEVTCzjwILi4c4zWww\nwR8LwCMxF2Np7ucB9P6f6itE7+F4h+CB19OrknsIfyNk33wWgKiz/f8AfItgaPzCgu2fJWjtXxTx\nPFHF+VkKVRMTCo2NC/09wT39MgUPI4J7OoQgywgAZjYQ2A54w93/Wcs5arBt7uc/y+4VaNR9qerz\nlLgvcd4TaWOqKyJRXRFQXRFOdUWTqUFQnZ8Co8zsRoLYwDcIWrHHE/RoXJ6LG43TncBjBCnmRhG0\noPcG9iLIrRzFtcD3gVvNbB4wjOAzhK0i+SeCdFtnmNlIgpy7T7r7/SXO/SPgGOAXZpafPLQbcEqu\n3D+KWMao4vwsharp2WhoXKi7/6+Z/QL4fO4z3sKGVRLvIUihlzeWIE3ePQSZPWo5BwBm9hSwjbtX\nGo7P9+T9p5ktJqh4r3T3Pr1sjbovNXyePvellnsiEpHqCtUVqiv6Ul2RZnGnMWqlF0HPzi+A/yF4\noL5NMKFpIXAy0K9g305Kp5LrDDm3E7LwCUFsoxdt2xG4jaCSWUmQ9mos0VPJ9Qe+STBk/CbBojc/\nIvgDd2BW0TlOJOjNyqfRu6zUuXPbNwd+CTxL8LB9NnffRhftV+5+dAFPRfg3achniXCdOwgqiGEJ\n/e31J8hY8ljucz4H/ATYpGi/8bnP1VXrOQr2fwV4LkLZhgCXE2QMceDpNN2Tcvel2nuil15RXqiu\nUF2huiJsP9UVKX5Z7oPGyszG5m7iUGBTd++ucIiIJMjMdiX4MnOyu/8m6fKIiEj6qK5oHc3KEHAe\nG2LyRCT9DiJ4yP826YKIiEhqqa5oEbGPEJjZRwgmD32PoGGgEQIRERERkZSIdVJxLo/sz4FzCOIX\nRUREREQkReIOGToN2JhgkpCIiIiIiKRMbCMEuTRn3wVOcPd1ueXWRUREREQkReIMGfoP4D53v6XS\njmY2A5gBMHjw4N3HjRtX0wXXr19Pv37NXkm9cVT+ZGW5/FkuO6j81Vq6dOkr7r550y6YUqNHj/bx\n48fXdOzrr7/O0KFDG1ugJlL5k5PlsoPKn6Qkyv7QQw9Fqy9iyvk6gSB3717AiNzrXwlyu44FBpc6\ndvfdd/daLVy4sOZj00DlT1aWy5/lsrur/NUCHvQU5K1O+qX6IruyXP4sl91d5U9SEmWPWl/ENUKw\nAzAQWBzy3rPAJQTLlYuIiIiISILiahDcC0wu2nYw8HXgUDYsXy0iIiIiIgmKpUHg7q8QLCfew8zG\n5/7zj651CEREREREUiG7s/hERERERKRuTWsQuPtl7m4aHRARaS1mdoyZ3Whmz5lZt5k9ZGafiHDc\nxmY228xeMrPXzezmgtFkERFpEo0QiIhIvU4HuoGvAEcAC4ErzewLFY77T2A68O/ANGA0cIeZDYqv\nqFLSmDFgVvk1ZkzSJRWRBotzHQIRKeG1117jpZdeYt26dQ053/Dhw3n00Ucbcq4kqPwbDBw4kC22\n2IJhw4Y15HxNcnhu7lje3Wa2FUFD4edhB5jZ1sApwMnufnlu21+BJ4ETgF/HW2Tp48UXo+8Xttho\nRwcsW9bYMonqiyJZLn+a6wo1CESa7LXXXuPFF19k7NixDB48mEas4r169Wo23XTTBpQuGSp/wN1Z\ns2YNzz33HEBmGgVFjYG8h4Gjyxx2YO7nvILzPGdm9wKHoAZBPMaMifzFf+komHkgzF4AOy6PcEDU\nBoVEpvqiryyXP811hUKGRJrspZdeYuzYsQwZMqQhD3dpHWbGkCFDGDt2LC+99FLSxanXJGBpmfd3\nBp4NmVf2aO49iUMVjYFJp8AtOwQ/l47q/d7hn+i9rYdCjBpK9YWEiaOuUINApMnWrVvH4MGDky6G\npNjgwYMbFh6QBDPbH/gYMLvMbiOBlSHbV+TekyYI+3KfbwysGATr+wU/842Ccg2FkjRyUDPVF1JO\nI+sKhQyJJEA9PVJOlv8+clmCrgRucPfLYjj/DGAGQEdHB11dXTWdp7u7u+Zj06BS+T901FFstGJF\n2XPkv9yvHASLxsHiS4LQoJkHBts812Xo/YLfT50Cfx0TNBC8oKGQP65ciFFxWbN8/5tZ9uHDh9Pd\n3djkjO+88w6rV69u6DmbKcvlj6Psa9eubcjfoxoEIiLSEGa2GXAr8DTwyQq7rwCGh2wfmXsvlLvP\nAeYATJw40Ts7O2sqa1dXF7UemwYVyx+xMRD25X72gqCBkH/P1sOItUH0T1hDId8ICGtc5BWXNcv3\nv5llf/TRRxseL5/lGHzIdvnjKPugQYPYbbfd6j6PQoZEpCazZs3CzNhhhx1C399hhx0wM2bNmhX5\nnE899RRmxk033dSz7Uc/+lFo74eZccEFF1Rb7LpcdtllmFnDe+xagZkNAW4CNgIOc/c3KhyyBBhn\nZkOLtu+ce08aqDg0qNQowMwDgy/yiy+BkWuh3/rg5+JL4MKbcg2D9cEx+YbC5x8oHWIkoroiG9Qg\nEMmakFzhmw4blshEvkGDBvHkk0/y4IMP9tr+pz/9iaeeeopBg+pPJ1/qIb948WKOOeaYus8v9TOz\nAcA1wA7Awe4eZZbbgtzPqQXn2QrYh2CUQRokLO5/9oLwL/ezc/8q+UbBoX/f0NtfqqFwwZ6lGxc9\nip5PnZMna7JxM6SkvlBdkX5qEIhkTTW5wmM2dOhQ9ttvP+bOndtr+9y5c9lvv/0YOrS487dx9tpr\nLzo6OmI7v1Tll8ChwHeBUWa2V8FrYwAzu8vM7sof4O7PApcAPzOzT5nZwQQpSJ8Gftf8j9CaSk0Q\nhvAv94VhPjsuh/lX9d1W3FCo1LgoSZON45eS+kJ1RfqpQSAidTnuuOO4+uqrcXcgyI989dVXc9xx\nx/Xar7Ozk2nTpvXa1tXVhZnxyCOPhJ57/PjxLF++nLPPPhszw8x6eoCKh4Hvvfde9tlnH4YNG8aw\nYcN4//vfzzXXXNPrfBdffDHvfe97GTRoEB0dHUybNo1Vq1YBQS/SEUccwZZbbsnQoUN5//vfzxVX\nXFHx869du5YzzzyTcePGsfHGG/O+972PW265pdc+N954I7vvvjtDhw5l5MiRfPCDH+See+6peO4M\nyfcFnw8sLnptmXuvf+5V6IvA5cBPgOuAV4ED3X1t3AVuF1FCgwq/3EdR3FAoNXKQf79smlJpG42o\nK/72t7+Fnlt1Rf3UIBCRuhx11FG8+OKL3HvvvQD88Y9/5OWXX+aoo46q+9zXX389w4cP55RTTmHx\n4sUsXryYD3zgA332e+211zjssMPYdtttue6667j22mv51Kc+xcqVG7JannvuuZx66qnsu+++/OEP\nf+BXv/pVrwweTz/9NB/+8Ie55JJLmD9/PkcffTQnnXQSV111VdkyTps2jSuuuIJvfetbzJ8/nz32\n2IMjjjiCv/zlLwD84x//YNq0aey3337Mnz+fK664gsMOO4xXX3217vuTFu4+3t2txOup3D6d7t5Z\ndNyb7n66u2/u7kPd/VB3fzKJz5BZISEhhasIRwkNKh4F6OEevCL0rpZqXNSUplRakuqKdNcVyjIk\nInUZMWIEBx98MHPnzmWfffZh7ty5HHzwwQwfHpZApjq77bYbAwYMYOutt2avvfYqud/SpUtZtWoV\nF1xwQU8GhwMP3BDAvHLlSr73ve/x5S9/mZ/85Cc92wsrosJeKnfnIx/5CM8++ywXX3wxn/jEJ0Kv\ne9ddd3HzzTdzyy23cMghh/Rcd+nSpfzHf/wH11xzDQ8//DCbbrop5513Xs9xhx56aJV3QqSECqEe\n+S/q+QxAI0JCg0IVNgKWLQvfpyg9br5xkVcuk1HU0QhpHaor0l1XaIRAROp23HHHce211/Lmm29y\n7bXX9hkCjtt2223HJptswvHHH88NN9zQq7cHgiHeNWvWcNJJJ5U8x4oVK/jiF7/INttsw8CBAxk4\ncCBz5sxh6dLSi+3eeeedjBkzhr322ou3336757X//vv3TJ5773vfy6pVqzjxxBNZsGABr7/+emM+\ntEiIsPCciqFB+ZGAwlepRkChCiMH5cKVpD2prkhvXaEGgYjU7YgjjqC7u5szzjiD119/ncMPP7yp\n1x85ciR33HEH69at4+Mf/zibb745U6ZM4YknngBg+fLgG9CWW25Z8hzTp0/n97//PV/96ldZsGAB\nf/rTnzj55JNZu7Z0OPsrr7zCsmXL2GyzzXoqhoEDBzJr1iyeeeYZAHbaaSduuOEGnnjiCQ499FBG\njx7N8ccfz8svv9zAOyBSPjynbGhQrZYtC29M5ESabBwW7qTsQy1LdUV66wqFDIlI3YYOHcphhx3G\nT3/6U4455pjQjBGDBg3irbfe6rVtRYXFk6qx1157cdttt7FmzRruvPNOTj/9dI4//njuu+8+Ro0K\nvhm98MILjB49us+xa9eu5aabbuIXv/gFp512Ws/29evXl73mZpttxtixY7niiivKZsmYMmUKU6ZM\nYdWqVdx88818+ctf5gtf+EKfjBsitao5PCfG7Cs1hysp+1DLUl2R3rpCIwQi0hCf+9znOPzww3s9\nJAttvfXWLFnSe72pBQsq5SWEjTbaqGzPS7HBgwdz+OGHc/LJJ/dkL5o0aRKDBw/mt7/9begxb775\nJuvXr2fjjTfu2bZ69WpuvPHGstfaf//9WbZsGZtssgkTJ07s8yo2fPhwjj/+eKZOnVoys5JILSKF\n59QaGlSHcuFKyj7UnlRXpLOu0AiBiDREZ2cnnZ2dJd+fOnUql1xyCV/5yleYMmUKCxcu5Lbbbqt4\n3p133pmbb76Zgw8+mE022YSddtqpz9LvN998M5deeikf+9jHeNe73sVzzz3HRRddxH777QcEk9nO\nPPNMzjjjDN566y0OPfRQ3nzzTW6++WbOOussxo4dyx577ME555zDsGHD6NevHz/4wQ8YPnw4r732\nWsmyffSjH+Wggw7iyCOP5Bvf+AYTJkzgtdde4y9/+Qtr167l+9//PhdddBGLFy/m4IMPZquttuLv\nf/8711xzDZ/+9Kej3ViRCGYvgEXjNowQRF4LoAmKJxvDhhGNlYOCcmuicftQXZHOukIjBCJZE3WI\nP2ULsUyZMoXvfe97XHvttUydOpWnn36a888/v+Jx5513HkOHDmXKlCnssccePPTQQ3322X777TEz\nvvWtb3HggQfyta99jYMPPphLL720Z59vfvOb/OpXv+LOO+/kyCOP5NRTT2XlypU9FcaVV17Jtttu\ny6c//Wm+9KUvcfTRR1d8EJsZ8+bN44QTTuBnP/sZBx10EKeeeiqLFy9m7733BmDXXXfl5Zdf5vTT\nT+fAAw/k3HPP5bOf/Sw//OEPq7l90u6K0ot2Tp7cK8tPpbUAmqrCs6fUYmkaKYhBBusL1RXJ1BXm\nBROA0mDixIlevLR1VF1dXWVbnWmn8ierWeV/9NFH+Zd/+ZeGnnP16tV9ekKyROXvq9zfiZk95O59\nx5nbTFvVF0UpPktZOioIE5q9oEQ2oSTlPsPhnwgmPq8v6JLstz4IK+oZSUi6rGU0829H9UVfWS5/\ns+sKiF5faIRAREQkY0rF38eSTajBImUfEpGmiqVBYGbTzGyRmS03s7Vm9piZfdvMNorjeiIiIu2i\nptV/UxQSkqrwJhEB4hshGAXcDXwGOAS4FDgD+Em5g0RERKS0SPH3CWQTqlbFxdJEpKliyTLk7hcV\nbVpoZsOAfzOzL3jaJi6IiIhkQLn0osWZfNIuLPtQj7A5Ex0dqWvYiLSKZs4hWA4oZEhERKRGmY6/\nrzdsSQuWicQm1gaBmfU3syFmtjfwReBXGh0QERGpTabj75ct6wlh6lq4cEM4UxEtWCbSfHGPELye\ne/0RuAf4aszXExERaWmtHH9f04RpEalb3CsVfwgYAuwJfAe4APjX4p3MbAYwA6Cjo4Ourq6aLtbd\n3V3zsWmg8ierWeUfPnw4q1evbug533nnnYafs5lU/r7Wrl2b6f+fpEZjxkQKjSkbf59RhROmvWDC\ndKs1ekTSKNYGgbv/Ofef95rZK8BvzWy2u/+jaL85wBwIFpqpdcGPzC00U0TlT1YzFyZr9MIkWV6o\nBVT+MIMGDWK33XZr6DklA+qNk09RetFqtdKEaZGsaeak4nzj4N1NvKaIxGDWrFmMHj26qdc0My64\n4IKe3zs7O5k2bVrk47u6ujAz/va3vwHw1ltvMWvWLP7yl79EPrb4NWDAhj6V6dOn99o+atQo9t57\nb37wgx+watWqKj6pyAZl4+kzkF60GpmeMC0lqb7IRn0Rd8hQoQ/nfj7ZxGuKSIv65S9/ycCBAyPv\n/4EPfIDFixez3XbbAcED/uyzz2b8+PG8//3vj3SOK664gm233bbndytKjbjzzjvzm9/8Bnfn1Vdf\n5b777uPHP/4xc+bM4e6772b8+PGRy5s1ZrY9wTyxScAE4I/u3lnhmPGE1wm/d/fjGlzEzMmH0Kwc\nBIvGtX7oTH5uRP4zj8jShGlJNdUXlcXSIDCz24A7gf8D3iFoDMwkeMj/o9yxIiJR7LLLLlXtP2zY\nMPbaa6+6rrnrrrvynve8p+T7Q4cO7XWNKVOmcOqpp7Lnnnty0kknsXDhwrqun3ITgEOB+4DoNW/g\n34H/V/D7K40qVFa1azx9vlEw88BgZKCVP6s0T1rqi8K5ZmmrL+IKGfoTMB24BrgaOBz4JvCpmK4n\nIgnKD5HeddddHHnkkQwdOpQddtiBBQsW8M477/DVr36V0aNHM3bsWH7yk94Llk+fPp19992XP/zh\nD+y8884MGjSIvffem0ceeaTsNcOGgP/6179y+OGHM2LECDbZZBP23HNP7rjjjl5lzA8B52P+Tzrp\npJ6h26eeeqpBd2SDrbfemu985zt0dXWxZMmShp8/Rea7+zh3P4agM6gaj7n7fQWvx+MoYJaUi6dv\ndfkJ06GNAbO+rzFjml5GqZ3qi9KSrC9iaRC4+5nu/h5338TdR7j7B9z95+6+Lo7ribSzpcuX8vHr\nP87S5UuTLgqnnnoqe++9N9dffz3bbLMN06ZN4/Of/zyrV6/myiuvZNq0acycOZP777+/13H//Oc/\nOf300znzzDO58sorWbVqFQcddBBr166NfO0lS5bw4Q9/mBdeeIELL7yQ66+/nqlTp/LMM8+E7n/3\n3XcD8O1vf5vFixezePFittxyy7LXeOedd3j77bd7XuvXr49Uto9+9KMA3HfffZE/T9a4e7SbIZG0\nVTy9FixrGtUXAdUXfTVzDoGINNjS5UuZdMkkVq5dyaRLJrH4lMXsOGrHxMrzqU99iq9+NVhuZOut\nt2bChAk89thjPQ/TAw44gN///vfMmzePD37wgz3HLV++nBtvvJEPfehDAOy+++5st912XHbZZZx2\n2mmRrn322WczfPhw/vjHPzJ48GBgw4M1zB577AHAdtttF3louDh29IwzzuDcc8+teNzWW28NwIv6\n4lLKb8xsM+Al4CrgDHdfk3CZEtVW8fSlJkIXxVwvHaVQonqovtggqfri61//esXjkqov1CAQyaj8\nw33FmhU4zoo1KxJ/yO+///49/7399tsDsN9++/Vs69evH9tuuy3PPfdcr+M233zznoc7wDbbbMPu\nu+/OAw88EPkBf/fdd3PCCSf0PNzjMHfu3J5JZgBbbbVVpOO0QHtJbwK/ABYArwGdwNeB7YAjkytW\nOiiefoN2m2DdaKovelN90ZcaBCIZNfP2maxcuxIneHg4zsq1K5l5+0zmHz8/kTKNGDGi57832mij\nPtvy24uHdjfffPM+59piiy144YUXIl97+fLlFYdw6zVhwoSyk4pLyVdoHRnOER8Hd38B+HzBpi4z\nexH4pZm9z93/p/iYVl/IsrPo93ILkKWx/FFFuf+duZ9RJlg3814082+nUQtZfunmL4XWF1+6+Utc\nPfXqus9fzptvvom793yON954A4CBAwf2+WyDBg3qta1///50d3f3bFu3bh2jR4/uc9xmm23GM888\n02v72rVre37Ph+/kf1++fDmbbbZZyXubL+Prr7/O6tWr6e7u7nPOUvLHjh8/np122qnXe/lFLNet\nW8f69etDz/XEE08A0f/tG7WIpRoEIhk1+6DZLHp2UU+Pj2GMGDSC2QfNTrpoVXv55Zf7bHvppZeY\nMGFC5HOMGjWqqgZEMy1YEAR+T5o0KeGSZMK1wC+B3YE+DYKWX8iyoyNSTPxbI0ems/wRVXP/oyxY\n1sx70cy/nUYtZHn+lPN7jRDk64vzp5wf+6KQG2+8MWbWc50hQ4YAQZad4msPGjSo17b+/fszYMCA\nnm0DBw7klVde6XPcq6++yoQJE3ptLzxX8XlGjRrFq6++WvKzF5cxnzK0uHxRji2UX8Ry4MCB9OvX\nL/RcixYtAmDy5MmR/m0atYhlMxcmE5EG2nHUjiw+ZTEjB4+kn/Vj5OCRiceE1urll1/ueQhCMGns\nz3/+M3vuuWfkc+y///5cffXVkSeW5UcwqpmIVotnn32W7373u0yePLlPb5GE8qKfrWnMmPCMOcWN\ngY6O0AXIFs2bl0y5E9BWE6xjovqiN9UXfWmEQCTD8g/5L938Jc6fcn4mH+4Q9NaccMIJnHvuuQwe\nPJizzjqLLbbYgunTp0c+x1lnncUee+zBRz7yEWbOnMmoUaN4+OGHGTVqFCeffHKf/TfaaCPe/e53\nc/XVV/Oe97yHQYMGseuuu/Y8+Gvx+uuvc9999+HurFy5kkWLFnHhhRey6aab8pvf/Kbm87aZfG7A\nhxItRdz6Sc9VAAAgAElEQVSiThjURPT2mmAdI9UXG6i+6EsjBCIZt+OoHbl66tWZfbgDvOtd7+LH\nP/4xs2bN4rjjjmPTTTfl9ttvZ9CgQZHPsdNOO3HvvfcyevRoPvOZzzB16lSuvfZattlmm5LHXHjh\nhbzyyisccMAB7LHHHjz//PN1fY4lS5YwadIk9t57bz75yU9y9913M3PmTB5++OGy5WgFZjbEzKaZ\n2TRgLLB5/nczG5Lb53Ezu6TgmFlmNtvMjjKzA8zsHOCnwDx3/2synyQZS0fB4Z8Ifkpf+UbBoX9X\nY6Aeqi8Cqi9CuHuqXrvvvrvXauHChTUfmwYqf7KaVf5HHnmk4ed87bXXGn7OZjnxxBN9t912S7oY\ndYnj/pf7OwEe9BQ8rwtfwHiCMJ+w1/jcPk8BlxUccxzwILAKeAt4HDgH2DjKNTNdXxQEAD02Ct/s\na3i/7wQ/HxtVFCAUIvHy1ylS+UODpUq8mqiZ9171RW9Zry+aXVe4R68vFDIkIiJ1c/enAKuwz/ii\n3+cCc+MrVfpFyaDTtiJOsK57YTMRUciQiIhIUspl0Gl7y5ZFGx8otbCZiESmEQIRSdRll13WkDzb\nIlk0e0Gw0FZ+hEAZdCIYMyb6yIEaCy1F9UV8NEIgIiKSkPxk2ZFrod/64KfChSpQhiaRhlODQERE\nJEHKoFMfZWgSqZ9ChkQS4O49Kx+KFAsSQ0g72XH5hlV3Jbr8pOyVg4LQq1ZsUKm+kFIaWVdohECk\nyQYOHMiaNWuSLoak2Jo1axg4cGDSxZA4Rc2Moww6JRVmaFpfkKGplUYKVF9IOY2sK9QgEGmyLbbY\ngueee4433nhDPcHSi7vzxhtv8Nxzz7HFFlskXRyJkzLo1K0dMjSpvpAwcdQVChkSabJhw4YB8Pzz\nz7Nu3bqGnHPt2rVVrdKYNir/BgMHDqSjo6Pn70QyThlxYtMOGZpUX/SV5fKnua5Qg0AkAcOGDWvo\nF76uri522223hp2v2VR+aVnKiBOb/GTs/ByCES2aoUn1RW9ZLn+ay66QIRERkSZRRpzGUoYmkcbQ\nCIGIiEgTtENGnCQoQ5NI/TRCICIiErN2yIjTNMrQJNJwsTQIzOwYM7vRzJ4zs24ze8jMPhHHtURE\nRNKuHTLiNI0yNIk0XFwjBKcD3cBXgCOAhcCVZvaFmK4npYwZA2a1v8aMSfoTiIhk3uwFwaRXWx/8\n3ooZcRIXtb5TvSbSR1xzCA5391cKfr/bzLYiaCj8PKZrtreoqe2qpcwXIiJ1a5eMOIlSRieRmsUy\nQlDUGMh7GNgqjusJkR9wpTJclM18oR4WEZG6KSNOcymjk0h0zZxUPAlY2sTrtaZSQ6JFwh6E+Ult\nt+zQezJbqe1lqYdFRKRq+Yw4agzEq6Z6TaSNNaVBYGb7Ax8DZjfjei0twhfxsAdhqQwXt29XPvNF\nNSMHnZMna+RARKSQMuI0nTI6iVQv9nUIzGw8cCVwg7tfVmKfGcAMgI6ODrq6umq6Vnd3d83HpkGU\n8ncW/b50VJClYvaCoMep8EHoBQ/C970YnuHi01NLZ76YvaCGnNkvvpjZf4Ms//1kueyg8rcCM9se\n+CrBaPAE4I/u3hnhuOHAzwg6jfoBNwFfdPfW6ENXppumK5fRSesViISLtUFgZpsBtwJPA58stZ+7\nzwHmAEycONE7Oztrul5XVxe1HpsGPeWPOEE4bJGbUg9C92ASW76hkM9wcfn1cPzRfbd//oHwhkW+\nUVDcECmU1X+DLP/9ZLnsoPK3iAnAocB9wMAqjrsa2BH4DLAe+CHwB2CfRhdQ2sPsBUGdWFyvKaOT\nSGmxhQyZ2RCCnp6NgMPc/Y24rtVyqmgMFA+Jfv6B8NR2F90cfJkfuRb6rQ9+Lr4EDvpH+PYL9izd\nw6LYTBEJMd/dx7n7McD/RTnAzCYBBwInuvt17n49cAKwt5kdEGNZG08pL1MjP3m7uF7TvA2R0uJa\nmGwAcA2wA3Cwu78Ux3XaRVgcf6mRgAv2LP0gLJXhImx7qZzZhSMHJWMzVQmKtB13X1/DYYcAL7r7\nfxec5wHgydx72aGUl6mijE4i1YlrhOCXBEPH3wVGmdleBa+NY7pmSyrVG19ukZtyD8JSGS6Kt5fq\nYSk3clCWKkER6WtnYEnI9kdz72WWUl4mTxmdRKKLq0GQ/3p4PrC46LVlTNfMnqIh5s7Jk3ulEC2X\nKaHSkGgjHoTVjBzkYzNVCYpIFUYCK0O2r8i9l0kKq0yIMjqJ1CyWScXuPj6O87acCr3mlTIl5L+w\nl5rcG6qjo3zWi6I1DfINi8LfS622GTbJWT0zItJIac1K15n7WSrTW+HzsBHXzXpmq1jKP3du9H3r\nuLbufbKyXP40lz32tKMSXXHmniiZEoq/sPfiXn0hOjoqNlTCGiJRKsGwBdQqNlBEpJWtADYP2T4y\n914fac9KFyXlZSOum/XMVlkuf5bLDip/ktJc9mauVCxlhA0x15UpodYh0WXLgoZE8atIcUhSuUqw\nLM0tEGlnSwifK1BqbkHqVQqrFBFJIzUIUiDKXIGSmRLCvry7N77XvUIDQ3MLRKQGtwJjzGzv/AYz\nmwhsm3svc5TyMqWUFlakLDUImqHUgyinUu96KjIlVBg5KFcJaoKdSOszsyFmNs3MpgFjgc3zv+fW\npcHMHjezS/LHuPtiYAFwuZkdZWYfA64A7nX3O5P4HI2glJcppLSwImWpQdAMFR4wNQ8xpyxTQlgl\nWG70Q0RayhYE689cA+wF7FLw+xa5fQYA/YuOOxa4B7gUuBx4CJjahPLGKhUdOVKSRq1FelODIAHF\nD6JIQ8zNCAtqgIbNLRCRTHH3p9zdSryeyu0z3t2nFx230t1PcvcR7j7M3Y9391eS+Ax1UcrLzNCo\ntUhfahA0WakHUWaHmOucWxBsVCyniGRcqbDKDHTktBONWouEU4OgiSo9iDI5xFxQCXYtXFjV3IKK\nFMspIiINpFFrkXBqEDRRzQ+ijA8xVxr9UCyniIg0g9LCioRTg6CRKmQTivQgatEh5lKjH4rlFBGR\nZlFaWJFwahA0UsQVftviQRRhVEOxnCKSKcpl3xIyO2dPJEZqEMQoLBSmbR5E5SbY5SiWU0QyRbns\nW0Ym5+yJxEgNgpiUC4XRgyhQcwYi9b6JSApo/lOGKC2sSFlqEMSg1lCYt0aObE4BU6LmECr1volI\nwjT/KWOUFlakLDUIYhApFCbkQbRo3rxEypukciFU6n0TkTTS/CcRaTVqENSiEdmEpEdYCJV630Qk\nrTT/SURajRoEtVA2ofpUiNFU75uIpJk6fUSk1ahB0ABtnU2oFqViOXPU+yYiaaZOnxaiVLIigBoE\ndVM2ocZT75uIpJ06fVqEUsmKAGoQ1KXm0BalNSsrUu+benBEJGHq9Gk9SmYh7UoNgjrUmk1Iac0q\nq6n3TT04IhIn5bJvaUpmIe0stgaBmW1vZheZ2V/N7B0z64rrWklRaEu8SvW+qQdHRBKhXPYtS8ks\npN3FOUIwATgUeAxYGuN1EqOJZQ0WoVdNPTgiItJoSmYh7S7OBsF8dx/n7scA/xfjdeJTYb0B0MSy\nhqqQfUg9OCIiEgeN+Eu7i61B4O7r4zp300SMSdfEsuZQD45IepnZLmZ2l5m9YWbPm9k5Zta/wjHj\nzcxDXnObVW4R0Ii/iCYVV6Hq2HVNLGso9eCIpJOZjQTuBBw4EjgHmAmcHfEU/w5MKnh9O4ZiRqfc\n9G1JI/7SzgYkXYCsyIerrBwEi8aFPCwKQlskHvmHdf7fYUSpdKTFOjo0yU8kXqcBg4Gj3P014A4z\nGwbMMrMf5baV85i73xd7KaNSbvq2lR/xF2k3qWgQmNkMYAZAR0cHXV1dNZ2nu7u75mPDdOZ+Fsau\ne0HseuGX0UZct9Hlb7a4yt9Z8N/5RsHMA4ORgajpSKOUK8v3P8tlB5W/BRwC3F70xX8u8ENgX2B+\nIqVqkKWjqnzmSHZ0dERr2GnEX1pcKhoE7j4HmAMwceJE7+zsrOk8XV1d1HpsOeVi1/M9CY24blzl\nb5Zmlb9UD065SjtKubJ8/7NcdlD5W8DOwN2FG9z9n2b2Ru69Sg2C35jZZsBLwFXAGe6+JpaSVqni\n6LBkm0aPRQDNIYhEsesJUzpSkbQbCawM2b4i914pbwK/AE4B9gcuAj5HMLqQOGU2E5F2kYoRgsSN\nGVN2yDBS7LrEp1QPTm6+QJSQLhFJH3d/Afh8waYuM3sR+KWZvc/d/6f4mGaEmHbmfkYZHU4qVCzr\nYWpZLn+Wyw4qf5LSXPbYGgRmNoRgYTKAscAwM5uW+/0Wd38jrmtXLUL8YE2x69IUUSptEYnVCmB4\nyPaRufeqcS3wS2B3oE+DoJkhprMXBGFC+c6GsNHhpELFsh6mluXyZ7nsoPInKc1ljzNkaAvgmtxr\nL2CXgt+3iPG6dSuVXlTrDaSTQrpEEreEYK5ADzMbBwzJvVcNL/qZGOWmb2NlUs92Tp6s1LPScuJc\nmOwpd7cSr6fium69aopFV/aBRKnSFkncrcBBZrZpwbZjgTXAPVWeKz+S/FAjClYv5aZvU0o9K21G\nk4oLRJpA5t73pSwFiatYaWtRIZE4XUgwQXiemR2Qi/OfBfykMBWpmT1uZpcU/D7LzGab2VG5484B\nfgrMc/e/NvkzlKTRYal6YVKRjFGDoEC5WHRJv6orbfXsiDSEu68gyBLUnyDF6NkEX+zPKtp1QG6f\nvCUE6xT8BrgFOB44L/czOVFHfTU63BaUxU7agbIMFYgygUxSJOKCMlpUSCR+7v4IsF+FfcYX/T6X\nlKQY7UWjvpKjLHbSLjRCUECx6BmzbFl4CFcB9eyIiEitFDkg7UINgiKaQNY6tKiQiIjUQ1nspF20\nV4OgVBqxIppA1hrUsyMiIvVQ5IC0i/ZqENQ7iVQTyDJFPTsiUlaZXPPKSCZ5ihyQdtBeDYIiZdOI\nKb1o5qlnR0TKUq55iUiRA9Lq2rZBoMmm7UE9OyISlXLNSw+lnpU205YNAk02bS9le3bClqNXeIBI\n21EnkfRSKoudO10LFypyQFpOWzYINNm0xdXbY6PwAJG2ok4iEWl3bdkg0GTTFhdhfQJQeICIBNRJ\nJCLtri0bBJpsKgoPEJE8dRKJSLtrzQZBhPUGNNm0fSk8QEQKqZNIRNpdazYIIsaAK41Ye1J4gIgU\nUyeRRDJmzIYEFFq7QlpIazYIilQdK640Yi1N4QEiEkadRFKR1q6QFtXyDYKKseJagKztKDxARADl\nmpe6KTmFtIqWbhAoVlxKUXiAiJTLNa9OIqlEySmklbR0g0Cx4lJO1AXLFBMqIiKF1OEoraalGwSK\nFZdetGCZiIg0gDocpdW0dINAseLSixYsExGRBlCHo7Sa2BoEZraLmd1lZm+Y2fNmdo6Z9Y/reqUo\nVlyqoZhQkdrU+sw3s+Fm9hszW2Fmq8zsCjPT/3mSaupwlFYTS4PAzEYCdwIOHAmcA8wEzm7YRUIW\nH+vJDVxEqeQkCsWESmxKLZbYIvNU6nzmXw10Ap8BpgN7AH9oaAHL1RctcP8lGepwlFgkVF8MaOjZ\nNjgNGAwc5e6vAXeY2TBglpn9KLetPo2I51YqOSlQLiZ0/lXJlk0yrvVzl9f0zDezScCBwL7u/t+5\nbc8B95vZAe5+Z0NK1/r3XxKS73AUaZiEnldxhQwdAtxeVAnMJagw9o3jgmXjvpVKTiJQTKg0SwvO\nU6n1mX8I8GK+MQDg7g8AT+bei0UL3n9pFq1dIU3WrOdVXA2CnYElhRvc/Z/AG7n3Gkpx39IIigmV\nZmjR51Wtz/w+x+U8WuG4mrXo/ZdmWbaMroULtXaFNEUzn1dxNQhGAitDtq/IvdcwivuWRqoYE9pC\ncd/SfC38vKr1md+0ugJa+v6LSItp9vMqrjkEVTGzGcAMgI6ODrq6uioe05n7GSXuO8r50qC7uzsz\nZQ2T1fJ3Fv1edUzoiy8m/rmzeu/zWr38nbmfrfS8Sorqi0Cr/z+TZlkuO6j8SYpS9s7cz2Y/r8xD\n8rDXfVKzl4BfuPvZRdtfB2a5+3mljp04caI/+OCDUS4C9G5Beb8g7rtPqEcMnzEOXV1ddHZ2Jl2M\nmmW2/GPGRJqcs3RU8D/i7AUhIwcJ/41l9t7ntHz5G/y8MrOH3H1i/SVvjFqf+WZ2NbC5u08u2n4z\ngLtPKXdd1RedSRejZlkuf5bLDip/kiKVPaH6Iq6QoSUUxX+a2ThgCOHxojVT3LfULbdgWZ+40AKK\nO5ZGaOHnVa3P/D7H5ZSaW1CXFr7/ItJimv28iqtBcCtwkJltWrDtWGANcE+jL6ZcwBInxR1LI7Xo\n86rWZ/6twBgz2zu/wcwmAtvm3mu4Fr3/krQWX2tEktHM51VcDYILgTeBeWZ2QC7mcxbwk4asQRBC\ni49JXMrF8Yn0UupLQZEWfF5Feuab2eNmdkn+d3dfDCwALjezo8zsY8AVwL0NW4MgRAvef0ma1rqQ\nmDTreRVLg8DdVwD7A/2B+QSrVf4UOKthF1EuYGkSrU8gkdVb2Wf0eVXFM39Abp9CxxKMIlwKXA48\nBExtaAFVX0iTaa0LqVlCz6vYsgy5+yPAfnGdPyzHb5Ynmkh65YfsJp0SjAyMUNyxRJTmieiNFuWZ\n7+7jQ7atBE7KveKh+kKaKB9munIQLBqn+kKqlNAaFnGFDIm0FMUdS7U0EV2k/WjOmWSVGgQiEZWN\n49PkMSmgLwUi7UlzziSr1CAQKaXe+DxNHmtb+lIg0p4050yySg0CkVJy6xP0eRXR5DEppi8FIu1J\na11IVqlBIFIHxYlLGH0pEGlfmnMmkaRs7Qo1CERqpDjxNlb0IO+cPLnPegP6UiDSvrTWhVSUsrUr\n1CAQqZHixNtYxAe0vhSItAmtdSF1Sjr8WA0CkRopTlzyqn6Q60uBSGspNees+JVQjnlJtzSEH6tB\nIFIjxYkLRHiQ60uBiIiUkJbwYzUIROpQMU484UlCEq+0PMhFRCSb0hJ+rAaBSJ1qihPXGgUtIS0P\nchERyaa0hB+rQSBSrSriv5OeJCTxSsuDXEREsikt4cdqEIhUq9zksQJpmCQk8UrLg1xEUiplueYl\nndKQploNApEYKLa8RZSqzAuk4UEuIimVslzzkl5Jp6lWg0AkBootbxFab0BEGkhhpNIjZWtXqEEg\nEgPFlrcmrTcgIrVSGKn0krK1K9QgEImBYstbj9YbEJFaKYxU0k4NApGY1LRGgSaWpZIqcxGph8JI\nJe3UIBCJUdWx5ZpYlkqqzCszs8+a2d/NbK2ZPWRm+0c4ZpaZecjr4GaUWaRZFEYqaacGgUgjRYwZ\n18SybFFlXp6ZfQK4ELgcOAT4P+AmM3tPhMNXAZOKXotjKqpIIhRG2sYKstV1Tp6c2tSzahCINFKp\nSUIFNLEse1SZVzQL+K27f9fdFwLTgceBb0Q49m13v6/otSrGsookQimK21RGUs/G0iAws2PNbJ6Z\nvZAb/p0ex3VEskax6Cml9QZqZmbbAjsCV+e3uft64BqC0QIRyVGKYklrhEBcIwTTgPHATTGdXyST\nFIueUlpvoB47534uKdr+KLCZmW1e4fgRZvaKma0zs4fN7KjGF1EkISnLNS/JSnOEQFwNgmPd/QPA\nV2I6v0gmKRY9G6rtwXlr5Mh4C5Ru+Q+/smj7iqL3wzwOfA04BjgaeB64To0CaRkpyzUvyUl7hMCA\nOE6aGy4WkSL5sJNJpwQjAyMUi546+Yf2ykGwaFzIv0/RnBCARV1ddDathPEzs+HAlpX2c/fiUYGq\nuPvviq47H1gEfAeYV6JsM4AZAB0dHXR1ddV07e7u7pqPTQOVPzlZLjuo/M3WmftZLkJg/lXBtiQ/\nVywNAhEpLd8omHlgMDIQuj5BsY4O9SA1QWEPjhf04LRho+0Y4OII+xkbRgKG03uUID8ysIKI3N3N\nbB7wQzPr7+7vhOwzB5gDMHHiRO/s7Ix6+l66urqo9dg0UPmTk+Wyg8qflNkLgk6mfP0SFiGQ5OdS\nliGRBGh9gnTSHI+Au//a3a3SK7d7fpRg56LT7Ay86u4vV3v53EtEpGWkPVtdpBGCuIePNQQcUPmT\nFWf5PzRyJButqNxRunRU6ZGDcmXTva9PZ+5nlB6csHImXf4kufsTZraUYFThdgAz65f7/dZqzmVm\nRjCX4H/CRgdEWsaYMdE6ejQ63FIqRggkKGrIUDXDx1XTEHBA5U9WrOV/9dXw7QXhQZVi18uVTfe+\nMaLM8QgrZ1rKn6BZwO/M7Cng/wEnAjsAx+d3MLN9gbuA/d39nty2e4DrCEYZhgKfBT4IfKyJZRdp\nvozkppfGy0cIpE2kkKEqh49FpEppzz7QMrTeQCzc/SrgNIIFyW4DdgUOc/e/FexmQH96dxw9DnwZ\nuAH4HbApMMXdb2xCsUVSI6256aUBMpJ6VpOKRVIgSvYBaYAq1xuQ6Nz9YsqMJLt7F0WjyO5+SszF\nEkm9ipnNJNsKQr7SPJoc10rFu5jZNDYM+040s2m5IWMRKaL1CZJRda+cFg8SkQbS6LCkRVxZhj5O\nsGz9f+V+/7fc72fHdD2RTIuUfSAs1GXMmMTKnHUVV4zU4kEiEjNlNpO0iKVB4O6zSswz6IzjeiKt\noKbYdU04q4l65UQkDTQ6LGmhdQhEUqTU+gSacNZY6pUTkTRIe256qUKppBUFr87Jk1M7sq8GgUiS\nIsSkVwxtkd7KPZRz1CsnImmhzGYtIuOpZNUgEEnSsmXhseo5Cm2pQYSHrXrlRCRNql69XlIvayP7\nahCIpJhCW+pX6qGsXjkRSUxGctNLbbI4sq8GgUiKRQptycclKvtQH5UeyuqVE5FElBodVmazzMvq\nyL4aBCIpVnNoS0pjFJuproeyeuVERKQGWR3ZV4NAJOXKhbZkLUaxmSI9lNUrJyIiDZTVpBVqEIhk\nQFhoSxZjFBuqVDahnKw+lEVEJLuymrRCDQKRNKoQspLVGMWGqhAWldWHsoi0sQi57DVPLP2ymLRC\nDQKRNKqQjjSrMYpxCgufyuJDWUTaWMZz2csGWUtaoQaBSAZFzT7ULr1K5cKnsvZQFhHJ0zyxDMl4\nKlk1CEQySNmHNqg5fCqlD2UREdA8scyJkEq2a+HC1CatUINAJKPaJvtQhcnDNWcTSulDWURE88Sk\n2dQgEMmwtsg+VGFUQ9mERKTVaJ6YNJsaBCJZ0ojsQxmfW1A8+qFsQiLSatTRIc2mBoFIlsSVfSgj\ncwtKjX4om1CyzOxYM5tnZi+YmZvZ9CqO/bCZ3W9ma83sSTP7YoxFFckEdXSkVAunhVWDQKSFVOpV\nSvXcggpzBSqNfiibUKKmAeOBm6o5yMy2B24HngQOBS4CfmJmn2l0AUWyRh0dKdTCaWHVIBBpIeV6\nlVI/t6DCA7Tm0Q9lE2qGY939A8BXqjzuq8DzwAnufre7/wCYA5xlVtAaFGlT6uhIt1R3slVJDQKR\nFhPWq5TFuQXFD9pIMbXKJpQId19f46GHAPPc/e2CbXOBrYH31F0wkazJeC77dpL6TrYqqUEg0oKK\ne5VSNbcgJDSoc/LkXruEPWgVU9tazGwoMA5YUvTWo7mfOze3RCIpECGXvTo6kteKaWHVIBBpAW+N\nHFn2/brmFjR60lSFRka5B61ialvKiNzPlUXbV+R+lv+jFhFJSCumhR3Q6BOa2TBgJsFQ8E7AGmAx\n8HV3X9ro64kILJo3j87Ozr5v5MKw81+kJ50SPLRGhMwtWDkIFo2r4ot2pdGDMWMijTAsHRU8RGcv\nCK5b7kGbH/WYf1WE8knNzGw4sGWl/dy9uHc/dmY2A5gB0NHRQVdXV03n6e7urvnYNFD5k5NU2T90\n1FFstGJFxf3eGjmSRfPmlXw/y/ceki1/Z+7n7AVBfbkiV1eFhbCGlTHN977hDQLgXcBngUuAM4Ah\nwDeB+81sV3d/JoZrikgF+UZB4Zfvwt54L+iNL2wUFH9h76XOeZ9hjZEoD9pQiqltpGOAiyPsV88f\nQH5kYHjR9vzIQOg3H3efQzDxmIkTJ3poQziCrq6u8EZ0Rqj8yUms7BEaAwAbrVhRtnxZvveQjvKX\n62TLCytjGspeShwhQ08C27n7me5+h7vfQJBObiBwcgzXE5GIqp1b0MhJU8VhSaVCgyDCXAHF1MbK\n3X/t7lbpVec1Xgeeoe9cgfzvTR99EMmSVspwk0WtFsLa8AaBu7/u7muKtr0KPA1s1ejriUgZFXrN\ny80tqDRpqlRlFLY9rGFRrjHSag9aKelWYKqZ9S/YdixBQ+FvyRRJJP1aLcNNVrVSWtimTCo2s82B\n7QHNIRBppgorG5fL3FPuC3upyihse6mGxecfKD/RuZUetK3OzHYxs2nAx3KbJprZNDPbt2Cffc3s\n7cJtwHkEKUb/y8wmm9nXgFOBc9wL/lBFpEcrZrjJjBZOC9usLEOzgW7gsiZdT0QiKtUbX2r04PMP\nhFdGt28Xvv20w8IbFhfsWWMa0Qw+aNvAx4FrgP/K/f5vud/PLtjHgP4UzDtw98eBgwk6jG4F/hWY\n6e6/bkKZRTKpFTPcZEYLp4WNNKm4nowTZvY54ATgaHcPreqVNSKg8icry+WvtuydRb+HZe4pNWmq\nVGX06anh292DY8MmCYdNdC7UtXBh+AdI2b9Tlv92GsHdZwGzKuzTRcgkZHe/F9gzjnKJtKKaEy+I\nlBE1y1BNGSfM7Ajg5wQpR68vdZCyRgRU/mRlufxVl72jI1JK0LAv7KUqo8uvh+OP7rv9opuDc5XK\nxlAujWhW/j2y/LcjItkSJcON1Cli2mw6OjI5GhAmUshQLRknzOzDBEvQX+ju58VReBGpUblhzyLF\nsfyl5h0c9I/SIUA1TRJWaJCISCglXohZlMZANftlQCxzCMxsAjAfuA34YhzXEJGYRPgiXqoyKldJ\nlZsk3LVwYUvEYIqINIsSLzRPO6R4jWOl4i0IGgLdwH8Ce9qGxYtec/dHGn1NEWmgUl/EixYhKxXq\nU5vTK48AABDfSURBVPVKwhoJEBGpLGKop56pjRW2gGYrNsLiGCHYhSCN3DhgIbC44PXLGK4nIs1Q\nbyXT0dEy2RhERJquhTPcpFU7pXht+AhBqUwSIpJxqmRERKSNlEvxWtVIeAY0ax0CEREREWm0MWOC\nkM5KrzFjki5p5pRaj6cVU7yqQSAiIiKSVW2YEadZSmXV0xwCEREREUmtdsiI00ztkuK14XMIRERE\nRKT52iUjTrNVnT0vgzRCICIiIpJx7ZQRJ3ZRs+q1UIpXNQhEREREMq5cRhypUhumeFWDQERERCTj\n2ikjTsMoQ1MPNQhEREREMq6dMuI0jDI09VCDQERERKQFtEtGnLi0c4YmZRkSERERaRHtkBEnDu2e\noUkjBCIiIiJZ1YYZcRpNGZrUIBARERHJrlIZcYobAC++CGZ0Tp7cdhNmK1GGJjUIRESkTmZ2rJnN\nM7MXzMzNbHrE42bl9i9+HRxzkUVanybMRqYMTWoQiIhI/aYB44Gbajh2FTCp6LW4YSUTEaC9J8xW\nogxNmlQsIiL1O9bd15vZJsBnqjz2bXe/L45CiUig3SfMRpFvFMw8MBgZaLf7oxECERGpi7uvT7oM\nIhJOE2ajy2doarfGAKhBICIiyRphZq+Y2Toze9jMjkq6QCKtRBNmy1CGph5qEIiISFIeB74GHAMc\nDTwPXKdGgUjjaMIsQSalwsxK+VfxhOqOjvCMTcuWJVPuJtIcAhER6cXMhgNbVtrP3ZfUcx13/13R\ndecDi4DvAPNKlG0GMAOgo6ODrq6umq7d3d1d87FpoPInJytl78z9zMfG5+cQjAiZMJuFz5NXy/3v\nrCLjUpz3Is1/O2oQiIhIsWOAiyPsZ428qLu7mc0Dfmhm/d39nZB95gBzACZOnOidnZ01Xaurq4ta\nj00DlT85WSx7pQmzWfo8jbj/S0clcy/S/LejkCEREenF3X/t7lbpFdflcy8RaaB2njBbKD/J+pYd\nNLm6UCwNAjO70MyWmFm3ma0ws/82swPiuJaIiLQGMzOCuQT/EzY6ICJVqGYibFh8fQuuYKyMS6XF\nFTI0GLgAeAzYCDgFuNXM9lG+aRGR1mJmuwC7AINymyaaWTfwsrvfk9tnX+AuYP+CbfcA1wFLgKHA\nZ4EPAh9r7icQaUGlJsJaxMG9FlzBuFzGpflXJVu2pMXSIHD3Ewt/N7PbgCeB4wE1CEREWsvHgbMK\nfv+33OseNsxtNKA/vecdPA58mWAC83rgz8AUd7815vKKSE65ePpWM3tBsDDbilyjoC0zLpXQlEnF\n7v6Oma0kGC0QEZEW4u6zgFkV9umiaBKyu58SW6FEpKJ2W8E4SsaldhXbpGILDDCzUWb2FWAH4NK4\nriciIiIi0bRrPH2+UXDo39UYKBRnlqFjgXXAK8B3gWPd/YEYryciIiIiEbTkCsalFiArooxLfUUK\nGapxkZrbgT2A0cAngblmdmhu2Lj4/FpoBpU/aVkuf5bLDiq/iEiztWQ8fb0ToavJzNRios4hqHqR\nGndfATyY+/U2M9sKOAf4SPFBWmgmoPInK8vlz3LZQeUXEWm2doinLzth2rXcSaFIIUMNWqTmYWDb\n+ossIiIiIvWqGE+f4fUJtABZdZqyUnFusZlJBKlHRURERCQFqo6nz8D6BO06YboeDW8QmNk+Zna9\nmX3azDrN7GjgBmAv4PuNvp6IiIiIRPPWyJGR9ls6Cg7/RDa/RLfkhOmYxTFC8AzwDvA9gonF5+eu\ns7e73xLD9UREREQkgkXz5gXx88WvAlkPt5m9IJgTYeuD31tiwnTMGt4gcPen3H2au2/t7hvnfh7m\n7osbfS0RERERaZxMhNtUSC+anxsxci30Wx/8bLUJ043WlDkEIiIiIpJ+mQi3iTCPQQuQVUcNAhER\nEREBshluU2q+gxYgi04NAhEREREBIobbpCgdaU3zHdp4AbJS1CAQERERkR41hdskkI400nyHsAnU\ny5Y1vaxppwaBiIiIiPRSKtymqelIQyYPd06e3PN2JuY7ZIQaBCIiIiLtLkIYTdPTkVYYdcjifIe0\nUoNAREREpN0tW1Z2fYJI4Tkxzy0oHp1QetHGUYNARERERMqqOTynQXMLSo1OKL1oY6hBICIiIiJl\nVQrPqWtuQYWFxiqNTii9aP3UIBARkZqZ2TAzO9vMHjCzVWa2zMyuN7MdIx6/i5ndZWZvmNnzZnaO\nmfWPu9wiUp1y4TkV5xaEfdkvfFUYRah5dELpRSNTg0BEROrxLuCzwO3ANOBUYEvgfjMbV+5AMxsJ\n3Ak4cCRwDjATODvOAotIbcLCcyLNLahS8WhDpMnDSi9aFzUIRESkHk8C27n7me5+h7vfABwKDARO\nrnDsacBg4KjcsRcSNAZON7NhsZZaRGpSHJ5Tqfe+VChRue3Fow2aPBw/NQhERKRm7v66u68p2vYq\n8DSwVYXDDwFud/fXCrbNJWgk7NvQgopIbSqE3ZTrvS8VSlRpe9hogyYPx0sNAhERaSgz2xzYHlha\nYdedgSWFG9z9n8AbufdEJGkV0pGW6r2H8C/3t29X+kt/pdEGTR6OjxoEIiLSaLOBbuCyCvuNBFaG\nbF+Re09EMiCs977Ul/tPTy39pb/mhcY0ebhuA5IugIiIpIuZDSeYGFyWuy8p3mZmnwNOAI5294b3\n45nZDGAGQEdHB11dXTWdp7u7u+Zj00DlT06Wyw6NK39n0e/53vu82Qtg0bhgBMD7bfhyf/n1cPzR\nfbfPXrChYTHplKCRMCJkrsBN8+ezySab9C1QBv5N0vy3owaBiIgUOwa4OMJ+1usXsyOAnwNfd/fr\nIxy/Ahgesn1k7r0+3H0OMAdg4sSJ3tnZGeEyfXV1dVHrsWmg8icny2WHBpa/o6NsutByX+7LfenP\nv58fMSgOD9pkk00ye//T/LejkCEREenF3X/t7lbpVXiMmX2YYELwhe5+XsRLLaForkAuVekQiuYW\niEjKVJhbAKUnAleaIKy5As2nEQIREamLmU0A5gO3AV+s4tBbga+a2abuvjq37VhgDXBPY0spIkko\nDiWqtL0szRWIjUYIRESkZma2BUFDoBv4T2BPM9sr99qlYL9tzOxtM/t0weEXAm8C88zsgNz8gFnA\nT4pSkYpIVtT7pb2jI3zkQQuNxUojBCIiUo9dgK1z/72w6L172DD30ID+FHREufsKM9sfuIBghGEl\n8FOCRoGIZJG+tGdS7A0CM/sS8DPgOnefFvf1RESkedy9i6LJxSX2eypsP3d/BNiv4QUTEZHIYg0Z\nyg0lzwJejvM6IiIiIiJSm7hHCL4P3ASMi/k6IiIiIiJSg9hGCMxsT+DjwDfiuoaIiIiIiNQnlgaB\nmRnB4jQ/cvfn4riGiIiIiIjUL64RgpOADuDHMZ1fREREREQaINIcAjMbDmxZaT93X5Lb9/vAF9x9\nTcTzzwBm5H7tNrPHohwXYjTwSo3HpoHKn6wslz/LZQeVv1rbNPFaqfXQQw+9YmZP13i4/uaSleXy\nZ7nsoPInKYmyR6ovzIuWmQ7dyewzwMWV9nN3M7MfEqSQ+2jBWzcBy4ETgdXu/k6UwlXLzB5094lx\nnLsZVP5kZbn8WS47qPzSfFn/N1P5k5PlsoPKn6Q0lz3SCIG7/xr4dcRz7gRMBFaEvLcC2Ae4N+K5\nREREREQkRnGkHf02wUJkhX4GrALOAv43hmuKiIiIiEgNGt4gcPe/FW8zs5XAK7kVLeM0J+bzx03l\nT1aWy5/lsoPKL82X9X8zlT85WS47qPxJSm3ZI80hqPsiZl0EDYJpsV9MREREREQia0qDQERERERE\n0im2lYqbycx2MbO7zOwNM3vezM4xs/5JlysKM9vezC4ys7+a2Tu50ZRMMLNjzOxGM3vOzLrN7CEz\n+0TS5YrKzKaZ2SIzW25ma83sMTP7tpltlHTZqmVmY3P/Bm5mmyRdnkrMbHqurMWv05IuW1RmNsDM\nvmFmfzezN83sWTP7adLlktKyXFeA6oskqb5IjuqL5ohjUnFTmdlI4E7gEeBIYDtgNkFj59sJFi2q\nCcChwH3AwITLUq3TgSeBrxDk1T0UuNLMRrv7zxMtWTSjgLuB84CVwJ7ALGAM8PnkilWT84BuYGjS\nBanSfkDheiVPJFWQGlxGUP6zgSXAOGCXJAskpbVAXQGqL5Kk+iJ5qi9ilPmQITP7/+3db4hUVRzG\n8e+TQbaw2WZEFJIGQdCLXpYW2z8Ig6BeSAZFSARFgb0IDaGCCgWDCstKgugfSJFZorgtVCQtFAWJ\nQWaClREIoWxsZWTarxfnjtg0rjPuOueemecDy505l4Hfi53zzLn3nHNXAMuBiyJiompbTvVFbbTV\nlaTTIuKf6vUG4NyIuDZvVe2pOvL9TW3rgfkRMS9TWVMiaSXwADAUhXw5JA0D7wOrSB39YET8nreq\nyUlaArxKAbW2ImkhsBm4PCJ25q7HTqz0rADnRd04L7rDedEdvTBl6CZgtKkzfws4E7gmT0nta3Tu\nJWru3CvbgQu6Xcs0OgAUcwu4mu7wPPAE5T65sUR3Ax/XuXO3/yk6K8B5UUPOC2tHEXnRCwOCS0m3\nX46KiJ+Ag9U56675wO7cRXRC0gxJA5KuBpYCL5VytQe4DzgDeCF3ISdpj6TD1Xzce3MX04ErgN2S\n1kqaqOakb5RU8o+bXuesqB/nRXc5L/IoIi+KX0MADJHm8zUbr85Zl0i6AbiVNBouyR+kThLgDWBZ\nxlraJmk28CRwZ0T8LSl3SZ3YBzwKfAHMAG4H1kkaiIhaLbQ6jvOBJcAOUu2DwFPAe5KuLOgHQj9x\nVtSI86K7nBdZFZEXvTAgsBqQNBdYD2yKiNeyFtO5BcAAaZHYY8Ba4P6sFbVnJfB5RGzNXUinImIU\nGD2maUTSTOARSWsKmBqh6u+WiDgAIGkfsI20cOyjjLWZ1ZrzIgvnRT5F5EUvDAjGgVkt2oeqc3aK\nSToHGAH2AndkLqdjEfFV9XJM0n7gdUlPR8SenHVNRtJlpCtrw5LOrpoHquMsSUci4s/Wn66tDcBt\nwFzqv3vEOPB9o3OvjAGHSDtH1KKDt/9wVtSA86L7nBfZFZEXvTAg2EXT/E9Jc0j/7LtafsKmjaQB\nYAtpYdXNEXEwc0lT1ejs5wG17eCBS0jbDn7W4tzPwCvAPV2taOqi6Vhn3wIzW7QLqPvVqn7lrMjM\neZGN8yKvIvKiFwYEI8AySYMR8VvVtpi0V+22fGX1PkmnA++QOpsFEfFL5pKmw1XV8YesVZzYGHBd\nU9tC4GHS/t51v2LSyiLSzhd7cxfShi3A401bKQ6TQndHvrJsEs6KjJwXWTkv8ioiL3phQLCOtNJ/\no6TVwMWkfaWfKWRf6QHSFxLgQuAsSYuq91trfgXlRVLtDwKzq0VLDdsj4q88ZbVH0gekBxV9Axwh\nde4PAW/X+fYvHN3C75Nj26p5uQCf1n2vZknvkhaIfU1aJLa4+ltawHxQgJdJ/c5mSatIi8RWAx9G\nxFjWyux4is4KcF7k5LzIx3nRHcUPCCJivNqtYC3pwQ+/As+SOvoSnEe6anKsxvt5wI9draYzN1bH\nNS3O1b12gC9JK//nAodJV0lWkH442Kn1HWlO6xzSbdOdwF0R8WbWqtoUEROSrgeeI+1lfwjYRHoK\nq9VQD2QFOC9ycl7k47zoguKfVGxmZmZmZievFx5MZmZmZmZmJ8kDAjMzMzOzPuYBgZmZmZlZH/OA\nwMzMzMysj3lAYGZmZmbWxzwgMDMzMzPrYx4QmJmZmZn1MQ8IzMzMzMz6mAcEZmZmZmZ97F9DdH0R\no8Wi6wAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f1872515080>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "i_widget = IntSlider(min=0, max=79, value=0)\n",
    "N_widget = IntSlider(min=2, max=200, value=50)\n",
    "M_widget = IntSlider(min=2, max=400, value=80)\n",
    "eps_widget = Dropdown(options=[10,2,1.,1e-1,1e-2,1e-3,1e-4,1e-5,1e-6,1e-7,\n",
    "                               1e-8,1e-9,1e-10,1e-11,1e-12,1e-13,1e-14,1e-15,1e-16,1e-17],\n",
    "                      value=1, description='eps')\n",
    "delta_widget = FloatSlider(min=0.4,max=2.0, step=0.1, value=1.0)\n",
    "fun_widget = Dropdown(options=[1,2,3],value=1, description='Function F')\n",
    "time_widget = FloatSlider(min=0.1,max=10,step=0.1,value=1,description=\"Max Time\")\n",
    "method_widget=Dropdown(options={\"Implicit FD\":1, \"Explicit FD\":2, \"FFT\":3},value=1)\n",
    "def update_i_values(*args):\n",
    "    i_widget.max = M_widget.value-1\n",
    "M_widget.observe(update_i_values, 'value')\n",
    "interact(comparison,\n",
    "         w=fixed(w),\n",
    "         wp=fixed(wp),\n",
    "         F=fixed(F),\n",
    "         i=i_widget,\n",
    "         eps=eps_widget,\n",
    "         T_max=time_widget,\n",
    "         N=N_widget,\n",
    "         M=M_widget,\n",
    "         method=method_widget,\n",
    "         delta=delta_widget,\n",
    "         num=fun_widget)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "It can be seen that the result are very near when there are not a big difference between scales. However, when the time starts to be too much large to the traditional computation, the method of multiscales starts to show differences in the numerical simulatons obtained."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Numerical Example 2**\n",
    "The second numerical example consider a variable speed $\\omega$ that depends of the current position $\\theta$. A positive parameter $\\delta$ will be empoyed in order to avoid having positions of $\\theta$ without speed.\n",
    "\\begin{align}\n",
    "\\omega(\\theta) &= \\sin^2{\\left(\\frac{\\theta}{2}\\right)} +\\delta\\\\\n",
    "\\omega'(\\theta) &= \\frac{1}{2}\\sin{(\\theta)}\\\\\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def w2(theta,delta=1):\n",
    "    return (np.sin(0.5*theta))**2+delta\n",
    "\n",
    "def wp2(theta):\n",
    "    return 0.5*np.sin(theta)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Widget Javascript not detected.  It may not be installed or enabled properly.\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "5f3c515dd1434a3ab17f4d1efa20c134"
      }
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<function __main__.plot_functions>"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "delta_widget = FloatSlider(min=0.4,max=2.0, step=0.1, value=1.0)\n",
    "fun_widget = Dropdown(options=[1,2,3],\n",
    "                      value=1, description='Function F')\n",
    "interact(plot_functions,w=fixed(w2),wp=fixed(wp2),F=fixed(F),delta=delta_widget,num=fun_widget)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "The delta parameter only desplaces the speed $\\omega$ in the vertical axis. This causes that with a greater delta, the speed of the solutions will be greater too. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "Now, with the same numerical methods used in previous example, some analysis of the solutions in cases where the speed is not constant will be realized. Consider the next code that uses the functions $\\omega(\\theta)$ and $w'(\\theta)$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Widget Javascript not detected.  It may not be installed or enabled properly.\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "c2851117b44f4b2e8dd3fcf7ce277d67"
      }
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<function __main__.simulation>"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "i_widget = IntSlider(min=0, max=79, value=0)\n",
    "N_widget = IntSlider(min=2, max=200, value=50)\n",
    "M_widget = IntSlider(min=2, max=400, value=80)\n",
    "eps_widget = Dropdown(options=[100.,10.,2.,1.,1e-1,1e-2,1e-3,1e-4,1e-5,1e-6,1e-7,\n",
    "                               1e-8,1e-9,1e-10,1e-11,1e-12,1e-13,1e-14,1e-15,1e-16,1e-17,0.],\n",
    "                      value=0, description='eps')\n",
    "delta_widget = FloatSlider(min=0.1,max=2.0, step=0.1, value=1.0)\n",
    "fun_widget = Dropdown(options=[1,2,3],\n",
    "                      value=1, description='Function F')\n",
    "fd_widget = Dropdown(options =  [\"Explicit\", \"Implicit\"], value=\"Implicit\", description = \"FD Mode\")\n",
    "time_widget = IntSlider(min=1,max=10,step=1,value=1,description=\"Max Time\")\n",
    "def update_i_values(*args):\n",
    "    i_widget.max = M_widget.value-1\n",
    "M_widget.observe(update_i_values, 'value')\n",
    "\n",
    "interact(simulation,\n",
    "         w=fixed(w2),\n",
    "         wp=fixed(wp2),\n",
    "         F=fixed(F),\n",
    "         i=i_widget,\n",
    "         N=N_widget,\n",
    "         M=M_widget,\n",
    "         eps=eps_widget,\n",
    "         t_max = time_widget,\n",
    "         fd_mode=fd_widget,\n",
    "         delta=delta_widget,\n",
    "         num=fun_widget)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "The first case with $\\varepsilon = 0$, where only there is a advection process is different when the speed is not a constant value. The maximum value is increasing in time where the solutions are pushed. When time is more than 1 second the maximum value of the solutions decreases in time. \n",
    "- FFT apparently has no problems of stability, but with the increment of the max time to 10 seconds the problem arises. \n",
    "- Implicit FD does not have problems of stability.\n",
    "- Explicit FD has the same problems of stability that FFT, but the problems starts later in this case."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "The second case $\\varepsilon > 0$, shows the expected diffusion with some advection that is not been well observed as the order of magnitude of $\\varepsilon$ increases. \n",
    "- FFT has problems of stability since the $\\varepsilon$ parameter has an efect in the moment of time when this problem arises. As greater this value is, sooner is the moment of inestability.\n",
    "- Explicit FD does not have problem of stability. \n",
    "- Implicit FD has similar problems that FFT and also, the $\\varepsilon$ parameter has an influence in the stability of the method, but FFT seems to be more inestable that Implicit FD."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The third case, and interest of this work is $0 < \\varepsilon \\ll 1$.\n",
    "- FFT does not have stability problems with the default parameters, but with an increment of max time for the simulation great oscilations can be seen at 3 seconds. In order to observate the different phenomena it is necessary to increase as much as possible the number of points $M$ and keep $N$ low. With a considerable value of max time is necesarry to keep as low as possible the time steps. With an order of magnitude lower than $10^{-2}$ the behaviour of the solution starts to seems similar to the first case when $\\varepsilon = 0$.\n",
    "- Implicit FD does not have stability problems and in this cases it is possible to observe the whole view of the different functions with different parameters. It can be seen that the variable speed causes thatthe maximum values in time oscillates, but with a inclination to decrease in time. With an order of magnitude lower than $10^{-2}$ the behaviour of the solution starts to seems similar to the first case when $\\varepsilon = 0$.\n",
    "- Explicit FD has problems of stability but it seems to be more stable that FFT. The behaviour is in general the same that the explained en previous cases. With an order of magnitude lower than $10^{-2}$ the behaviour of the solution starts to seems similar to the first case when $\\varepsilon = 0$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "In general, the methods does not let to visualize what is happening with low orders of magnitude for the diffusion coefficient (lower to $10^{-2}$). So, again, with using the multiscales of time the behaviour of this situation could be observed."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "The next plot shows the displacement of a certain value in order of time. In a different form to the example with constant value for speed, in this case is expected that the dispacement in time be non-linear."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Widget Javascript not detected.  It may not be installed or enabled properly.\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "0dbcfae7a1a240318805565eb6be5ac0"
      }
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<function __main__.plot_Theta>"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "delta_widget = FloatSlider(min=0.1,max=2.0, step=0.1, value=1.0)\n",
    "interact(plot_Theta,w=fixed(w2),delta=delta_widget)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this case, it can be seen that the role of $\\delta$ parameter is important. If $\\delta$ was zero, then for certain values of positions could not be movement, since $\\omega(\\theta) = 0$ for some $\\theta$. This could causes an infinity value for the period. It can be seen that with $\\delta = 0.1$ the period is equal to 19 seconds approx, and with increments of this parameter the period quickly reduce its value. \n",
    "\n",
    "The displacement is not a constant, such as in the first numerical example, and seems to be like a curves that increases its value in time. A important thing to control with $\\delta$ parameter is that if this value is too high, then the displacement is like a stright line and start to seems like the first numerical example. So is important to mantaint it bounded. \n",
    "\n",
    "The period obtain in this case is not simple to calculate. In first numerical case, period was a multiple of $2\\pi$ because, changes in $\\delta$ only changes the inclination of the straight line that defines the displacement. Using the same value of delta, the period obtain in this cases is lower that the linear case. The variable speed helps to move faster that with a constant value."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, the function $A_0$ calculates from the diffusion equation is shown below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Widget Javascript not detected.  It may not be installed or enabled properly.\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "233e2ca42d2a4ae9ab952a840e64695e"
      }
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<function __main__.print_simulation_a0>"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "i_widget = IntSlider(min=0, max=499, value=0)\n",
    "eps_widget = Dropdown(options=[100.,10.,2.,1.,1e-1,1e-2,1e-3,1e-4,1e-5,1e-6,1e-7,\n",
    "                               1e-8,1e-9,1e-10,1e-11,1e-12,1e-13,1e-14,1e-15,1e-16,1e-17,0.],\n",
    "                      value=0, description='eps')\n",
    "delta_widget = FloatSlider(min=0.4,max=2.0, step=0.1, value=1.0)\n",
    "fun_widget = Dropdown(options=[1,2,3],\n",
    "                      value=1, description='Function F')\n",
    "interact(print_simulation_a0,w=fixed(w2),F=fixed(F),i=i_widget,eps=eps_widget,delta=delta_widget,num=fun_widget)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It can be seen that in this case the shape of the function is a little bit different for function 1. Functions 2 and 3 appears to be similar to the linear speed case, but the amplitude in the initial value is different. In order of advances of time in slow scale it can be seen a regular behaviour that expreses the diffusion process. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, with that result, the function $u_0$ can be reconstructed using the multiscale. The following code shows the result of computing $u_0$ using the previous values for $A_0$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Widget Javascript not detected.  It may not be installed or enabled properly.\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "f585d556977b4c049dc43d5c9593443e"
      }
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<function __main__.print_simulation_u0>"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "i_widget = IntSlider(min=0, max=499, value=0)\n",
    "eps_widget = Dropdown(options=[100.,10.,2.,1.,1e-1,1e-2,1e-3,1e-4,1e-5,1e-6,1e-7,\n",
    "                               1e-8,1e-9,1e-10,1e-11,1e-12,1e-13,1e-14,1e-15,1e-16,1e-17],\n",
    "                      value=1, description='eps')\n",
    "delta_widget = FloatSlider(min=0.4,max=2.0, step=0.1, value=1.0)\n",
    "fun_widget = Dropdown(options=[1,2,3],\n",
    "                      value=1, description='Function F')\n",
    "interact(print_simulation_u0,w=fixed(w2),F=fixed(F),i=i_widget,eps=eps_widget,delta=delta_widget,num=fun_widget)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It can be seen that the functions are reconstructed in a good way and shows similar behaviour that the Implicit Finite Difference (more stable). It can be seen clearly that there is two process involved, the advection that produces movement of the values around the spatial axis and diffusion that shows that the values tends to spread."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, the comparison between a traditional approach with one scale and the approximations with multiscales are shown next."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Widget Javascript not detected.  It may not be installed or enabled properly.\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "add2746ddfca4b83b965c231da2fb414"
      }
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<function __main__.comparison>"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAwQAAAFKCAYAAACjNpjXAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xm8HFWZ//HPkwWyQBYSuIEQiexDFEUCEgW5AWQLi4Eg\niCgBNOCMKxlXRAIybhiVERWCIDICkSUIYQ1LLiO/BBDEcRwIEVlkC0tIQi4kEMjz+6O6b/r2re6u\nXqqrqvv7fr36dbnVtZyuXOr0Oec5zzF3R0RERERE2lO/pAsgIiIiIiLJUYNARERERKSNqUEgIiIi\nItLG1CAQEREREWljahCIiIiIiLQxNQhERERERNqYGgQiIiIiIm1MDQIRERERkTamBkFMzKzTzNzM\nprdqOdLyGUVEsiotz1HVFSLtTQ2CKpnZtmY2x8yWmNkbZrbCzB41s9+a2eSky9doZvZ+M5tlZuOT\nLku9qv0sZraJmb2Tq8iivDaLsez9zOwrub+7tWb2jJnNNrOhzTxHhfOfZGZfasS5Il4v9fdE2pfq\niuxSXaG6Io5zpN2ApAuQJWY2EbgHWAdcDvwfMBjYATgQWA0szO3+37n31jW/pA31fuAsoAt4qui9\nrH3Gcp8lzADgxKJtnwM+BPw78GLB9jfd/dX6i1jST4EvAtcDs4F/yf2+m5kd4O7rm3SOcn4E3A+c\nX+d5osrCPZE2pLpCdQWqK8pRXZFG7q5XxBcwH3DgfSXeH5N0GUPK1Jkr8/Qaj5+eO74z6c/SgHtR\n92cBHgLWAAOaWO4JwHrguqLtX8h9nuObcY4K598+d55v657o1e4v1RXJf54674XqCtUVTbsnaXkl\nXoAsvYAlwCsR9+3zcC14yOwPfAd4OvfAuB/YK7fPvsC9wOvAC8CZReedlTvH+JBrPgV0RSjHpsC5\nueu+ArwJPA78ABgScq3i12Wlzp3bPhr4BfAM8Fbu5y+AUUX75e/HfgS9KP/IlWUpcGLE+9yQzxLx\nWgNz57+/yX935+bKuk/R9kG5v5NbmnGOMuf+Q4l7e2673hO92vuF6grVFaorws6tuiLFL4UMVecf\nwE5mdpS7z6vjPD8A+hMMl20EzAQWmNmngUuAOcAVwMeBc8zsSXf/XX1F72Us8BngOuBK4G2CyuVr\nwG7AQbn95gFbAjOA7wGP5rb/o9SJzWw4sIigF+BS4M+5c34O2M/M9nT31UWHfY9gOPkigofo54DL\nzOxxd/9/SX2WEBMI/r0errSjmfUDqokTfdVLDznuQdA78UDhRndfa2Z/yb1fSSPOUcocgr/nwwj+\n7bpz2xcX79jA+5L2eyLtTXVFQHVFBaorANUV6ZB0iyRLL2ASQS+GE/RMXErwR/0vIft2UrrX58/A\nRgXbj8htXwdMLNi+EUHPz+KCbbOov9dnI2BgyPHfze27Z0iZOyN+xv/IbfvXon3/Lbf9uyHnfrjo\nfowleNhfFeHfpCGfJeK//0m542dE2Hc84T0hpV59/j0LzvW/wIsl3rs6d/xGFcpT9zkqnP8G4KVm\n3Zcs3BO92veF6ooon1F1hauuiPu+ZOGepOWlEYIquPtiM9udoJfmEIL/6U8CMLM/Ejzsnohwql+5\n+1sFv/8x9/N+d3+w4HpvmdkDwIcb8gEKzpv/bzMbQDCU2h+4E/g28EGKWsJVmAq8TNATUOgigkla\nU4Ezi977ZWGZ3P05M1tKMAGvrJg/S7EP5H5W7PUBlgEfreLcy8q8N4Sg0guztmCft0rs06hzlPMB\nmntfsnBPpE2prohEdUVAdUU41RVNpgZBldz9fwl6DzCzbQiGHD8D7APcYGa7Fz3Aw/SqCNx9hZkB\nPBmy7wpgVJ3F7sPM/hU4jWBoszj97Mg6Tv1u4EF3f7two7u/nXtwfyDkmLCKcTmwTZQLxvhZiu1G\nMMz8v5V2dPe1BBVNI7wBbFHivUEF+8R9jlBmtjmwNVAxVKGB9yXV90REdUVFqitQXVGK6ormU4Og\nDu7+NHC5mf0XQc/Nh4E9CSZ6lfNOldt7XbbMe5H+Pc3sdIK0WQuA/wSeJ2jZjgUuo/nrU5T63Fbp\nwGZ9llw84/uAR3MPqkr79wc2r+ISL7t7qfvwPLCLmW3s7sW9FGMJJi9W+mLRiHOUkq+4/1xpxwbe\nl7TfE5EeqisaRnWF6opSVFfUSQ2CBnB3N7P7CR7yY2O+XD5/8WYU5Ec2s0EEE6Eej3COT+WOPcQL\nJuGY2cEh+5arVMI8QTCZbkBhz09uiHZHwnt46hHnZym0A7AJER5kOeMI78Ur5d2Uznf9J4Lc5Xuy\nIWQg/2/+foIc35U04hyl7Jb7GeXeNOq+pP2eiPShuqIX1RUB1RXhVFc0mRoEVTCzjwILi4c4zWww\nwR8LwCMxF2Np7ucB9P6f6itE7+F4h+CB19OrknsIfyNk33wWgKiz/f8AfItgaPzCgu2fJWjtXxTx\nPFHF+VkKVRMTCo2NC/09wT39MgUPI4J7OoQgywgAZjYQ2A54w93/Wcs5arBt7uc/y+4VaNR9qerz\nlLgvcd4TaWOqKyJRXRFQXRFOdUWTqUFQnZ8Co8zsRoLYwDcIWrHHE/RoXJ6LG43TncBjBCnmRhG0\noPcG9iLIrRzFtcD3gVvNbB4wjOAzhK0i+SeCdFtnmNlIgpy7T7r7/SXO/SPgGOAXZpafPLQbcEqu\n3D+KWMao4vwsharp2WhoXKi7/6+Z/QL4fO4z3sKGVRLvIUihlzeWIE3ePQSZPWo5BwBm9hSwjbtX\nGo7P9+T9p5ktJqh4r3T3Pr1sjbovNXyePvellnsiEpHqCtUVqiv6Ul2RZnGnMWqlF0HPzi+A/yF4\noL5NMKFpIXAy0K9g305Kp5LrDDm3E7LwCUFsoxdt2xG4jaCSWUmQ9mos0VPJ9Qe+STBk/CbBojc/\nIvgDd2BW0TlOJOjNyqfRu6zUuXPbNwd+CTxL8LB9NnffRhftV+5+dAFPRfg3achniXCdOwgqiGEJ\n/e31J8hY8ljucz4H/ATYpGi/8bnP1VXrOQr2fwV4LkLZhgCXE2QMceDpNN2Tcvel2nuil15RXqiu\nUF2huiJsP9UVKX5Z7oPGyszG5m7iUGBTd++ucIiIJMjMdiX4MnOyu/8m6fKIiEj6qK5oHc3KEHAe\nG2LyRCT9DiJ4yP826YKIiEhqqa5oEbGPEJjZRwgmD32PoGGgEQIRERERkZSIdVJxLo/sz4FzCOIX\nRUREREQkReIOGToN2JhgkpCIiIiIiKRMbCMEuTRn3wVOcPd1ueXWRUREREQkReIMGfoP4D53v6XS\njmY2A5gBMHjw4N3HjRtX0wXXr19Pv37NXkm9cVT+ZGW5/FkuO6j81Vq6dOkr7r550y6YUqNHj/bx\n48fXdOzrr7/O0KFDG1ugJlL5k5PlsoPKn6Qkyv7QQw9Fqy9iyvk6gSB3717AiNzrXwlyu44FBpc6\ndvfdd/daLVy4sOZj00DlT1aWy5/lsrur/NUCHvQU5K1O+qX6IruyXP4sl91d5U9SEmWPWl/ENUKw\nAzAQWBzy3rPAJQTLlYuIiIiISILiahDcC0wu2nYw8HXgUDYsXy0iIiIiIgmKpUHg7q8QLCfew8zG\n5/7zj651CEREREREUiG7s/hERERERKRuTWsQuPtl7m4aHRARaS1mdoyZ3Whmz5lZt5k9ZGafiHDc\nxmY228xeMrPXzezmgtFkERFpEo0QiIhIvU4HuoGvAEcAC4ErzewLFY77T2A68O/ANGA0cIeZDYqv\nqFLSmDFgVvk1ZkzSJRWRBotzHQIRKeG1117jpZdeYt26dQ053/Dhw3n00Ucbcq4kqPwbDBw4kC22\n2IJhw4Y15HxNcnhu7lje3Wa2FUFD4edhB5jZ1sApwMnufnlu21+BJ4ETgF/HW2Tp48UXo+8Xttho\nRwcsW9bYMonqiyJZLn+a6wo1CESa7LXXXuPFF19k7NixDB48mEas4r169Wo23XTTBpQuGSp/wN1Z\ns2YNzz33HEBmGgVFjYG8h4Gjyxx2YO7nvILzPGdm9wKHoAZBPMaMifzFf+komHkgzF4AOy6PcEDU\nBoVEpvqiryyXP811hUKGRJrspZdeYuzYsQwZMqQhD3dpHWbGkCFDGDt2LC+99FLSxanXJGBpmfd3\nBp4NmVf2aO49iUMVjYFJp8AtOwQ/l47q/d7hn+i9rYdCjBpK9YWEiaOuUINApMnWrVvH4MGDky6G\npNjgwYMbFh6QBDPbH/gYMLvMbiOBlSHbV+TekyYI+3KfbwysGATr+wU/842Ccg2FkjRyUDPVF1JO\nI+sKhQyJJEA9PVJOlv8+clmCrgRucPfLYjj/DGAGQEdHB11dXTWdp7u7u+Zj06BS+T901FFstGJF\n2XPkv9yvHASLxsHiS4LQoJkHBts812Xo/YLfT50Cfx0TNBC8oKGQP65ciFFxWbN8/5tZ9uHDh9Pd\n3djkjO+88w6rV69u6DmbKcvlj6Psa9eubcjfoxoEIiLSEGa2GXAr8DTwyQq7rwCGh2wfmXsvlLvP\nAeYATJw40Ts7O2sqa1dXF7UemwYVyx+xMRD25X72gqCBkH/P1sOItUH0T1hDId8ICGtc5BWXNcv3\nv5llf/TRRxseL5/lGHzIdvnjKPugQYPYbbfd6j6PQoZEpCazZs3CzNhhhx1C399hhx0wM2bNmhX5\nnE899RRmxk033dSz7Uc/+lFo74eZccEFF1Rb7LpcdtllmFnDe+xagZkNAW4CNgIOc/c3KhyyBBhn\nZkOLtu+ce08aqDg0qNQowMwDgy/yiy+BkWuh3/rg5+JL4MKbcg2D9cEx+YbC5x8oHWIkoroiG9Qg\nEMmakFzhmw4blshEvkGDBvHkk0/y4IMP9tr+pz/9iaeeeopBg+pPJ1/qIb948WKOOeaYus8v9TOz\nAcA1wA7Awe4eZZbbgtzPqQXn2QrYh2CUQRokLO5/9oLwL/ezc/8q+UbBoX/f0NtfqqFwwZ6lGxc9\nip5PnZMna7JxM6SkvlBdkX5qEIhkTTW5wmM2dOhQ9ttvP+bOndtr+9y5c9lvv/0YOrS487dx9tpr\nLzo6OmI7v1Tll8ChwHeBUWa2V8FrYwAzu8vM7sof4O7PApcAPzOzT5nZwQQpSJ8Gftf8j9CaSk0Q\nhvAv94VhPjsuh/lX9d1W3FCo1LgoSZON45eS+kJ1RfqpQSAidTnuuOO4+uqrcXcgyI989dVXc9xx\nx/Xar7Ozk2nTpvXa1tXVhZnxyCOPhJ57/PjxLF++nLPPPhszw8x6eoCKh4Hvvfde9tlnH4YNG8aw\nYcN4//vfzzXXXNPrfBdffDHvfe97GTRoEB0dHUybNo1Vq1YBQS/SEUccwZZbbsnQoUN5//vfzxVX\nXFHx869du5YzzzyTcePGsfHGG/O+972PW265pdc+N954I7vvvjtDhw5l5MiRfPCDH+See+6peO4M\nyfcFnw8sLnptmXuvf+5V6IvA5cBPgOuAV4ED3X1t3AVuF1FCgwq/3EdR3FAoNXKQf79smlJpG42o\nK/72t7+Fnlt1Rf3UIBCRuhx11FG8+OKL3HvvvQD88Y9/5OWXX+aoo46q+9zXX389w4cP55RTTmHx\n4sUsXryYD3zgA332e+211zjssMPYdtttue6667j22mv51Kc+xcqVG7JannvuuZx66qnsu+++/OEP\nf+BXv/pVrwweTz/9NB/+8Ie55JJLmD9/PkcffTQnnXQSV111VdkyTps2jSuuuIJvfetbzJ8/nz32\n2IMjjjiCv/zlLwD84x//YNq0aey3337Mnz+fK664gsMOO4xXX3217vuTFu4+3t2txOup3D6d7t5Z\ndNyb7n66u2/u7kPd/VB3fzKJz5BZISEhhasIRwkNKh4F6OEevCL0rpZqXNSUplRakuqKdNcVyjIk\nInUZMWIEBx98MHPnzmWfffZh7ty5HHzwwQwfHpZApjq77bYbAwYMYOutt2avvfYqud/SpUtZtWoV\nF1xwQU8GhwMP3BDAvHLlSr73ve/x5S9/mZ/85Cc92wsrosJeKnfnIx/5CM8++ywXX3wxn/jEJ0Kv\ne9ddd3HzzTdzyy23cMghh/Rcd+nSpfzHf/wH11xzDQ8//DCbbrop5513Xs9xhx56aJV3QqSECqEe\n+S/q+QxAI0JCg0IVNgKWLQvfpyg9br5xkVcuk1HU0QhpHaor0l1XaIRAROp23HHHce211/Lmm29y\n7bXX9hkCjtt2223HJptswvHHH88NN9zQq7cHgiHeNWvWcNJJJ5U8x4oVK/jiF7/INttsw8CBAxk4\ncCBz5sxh6dLSi+3eeeedjBkzhr322ou3336757X//vv3TJ5773vfy6pVqzjxxBNZsGABr7/+emM+\ntEiIsPCciqFB+ZGAwlepRkChCiMH5cKVpD2prkhvXaEGgYjU7YgjjqC7u5szzjiD119/ncMPP7yp\n1x85ciR33HEH69at4+Mf/zibb745U6ZM4YknngBg+fLgG9CWW25Z8hzTp0/n97//PV/96ldZsGAB\nf/rTnzj55JNZu7Z0OPsrr7zCsmXL2GyzzXoqhoEDBzJr1iyeeeYZAHbaaSduuOEGnnjiCQ499FBG\njx7N8ccfz8svv9zAOyBSPjynbGhQrZYtC29M5ESabBwW7qTsQy1LdUV66wqFDIlI3YYOHcphhx3G\nT3/6U4455pjQjBGDBg3irbfe6rVtRYXFk6qx1157cdttt7FmzRruvPNOTj/9dI4//njuu+8+Ro0K\nvhm98MILjB49us+xa9eu5aabbuIXv/gFp512Ws/29evXl73mZpttxtixY7niiivKZsmYMmUKU6ZM\nYdWqVdx88818+ctf5gtf+EKfjBsitao5PCfG7Cs1hysp+1DLUl2R3rpCIwQi0hCf+9znOPzww3s9\nJAttvfXWLFnSe72pBQsq5SWEjTbaqGzPS7HBgwdz+OGHc/LJJ/dkL5o0aRKDBw/mt7/9begxb775\nJuvXr2fjjTfu2bZ69WpuvPHGstfaf//9WbZsGZtssgkTJ07s8yo2fPhwjj/+eKZOnVoys5JILSKF\n59QaGlSHcuFKyj7UnlRXpLOu0AiBiDREZ2cnnZ2dJd+fOnUql1xyCV/5yleYMmUKCxcu5Lbbbqt4\n3p133pmbb76Zgw8+mE022YSddtqpz9LvN998M5deeikf+9jHeNe73sVzzz3HRRddxH777QcEk9nO\nPPNMzjjjDN566y0OPfRQ3nzzTW6++WbOOussxo4dyx577ME555zDsGHD6NevHz/4wQ8YPnw4r732\nWsmyffSjH+Wggw7iyCOP5Bvf+AYTJkzgtdde4y9/+Qtr167l+9//PhdddBGLFy/m4IMPZquttuLv\nf/8711xzDZ/+9Kej3ViRCGYvgEXjNowQRF4LoAmKJxvDhhGNlYOCcmuicftQXZHOukIjBCJZE3WI\nP2ULsUyZMoXvfe97XHvttUydOpWnn36a888/v+Jx5513HkOHDmXKlCnssccePPTQQ3322X777TEz\nvvWtb3HggQfyta99jYMPPphLL720Z59vfvOb/OpXv+LOO+/kyCOP5NRTT2XlypU9FcaVV17Jtttu\ny6c//Wm+9KUvcfTRR1d8EJsZ8+bN44QTTuBnP/sZBx10EKeeeiqLFy9m7733BmDXXXfl5Zdf5vTT\nT+fAAw/k3HPP5bOf/Sw//OEPq7l90u6K0ot2Tp7cK8tPpbUAmqrCs6fUYmkaKYhBBusL1RXJ1BXm\nBROA0mDixIlevLR1VF1dXWVbnWmn8ierWeV/9NFH+Zd/+ZeGnnP16tV9ekKyROXvq9zfiZk95O59\nx5nbTFvVF0UpPktZOioIE5q9oEQ2oSTlPsPhnwgmPq8v6JLstz4IK+oZSUi6rGU0829H9UVfWS5/\ns+sKiF5faIRAREQkY0rF38eSTajBImUfEpGmiqVBYGbTzGyRmS03s7Vm9piZfdvMNorjeiIiIu2i\nptV/UxQSkqrwJhEB4hshGAXcDXwGOAS4FDgD+Em5g0RERKS0SPH3CWQTqlbFxdJEpKliyTLk7hcV\nbVpoZsOAfzOzL3jaJi6IiIhkQLn0osWZfNIuLPtQj7A5Ex0dqWvYiLSKZs4hWA4oZEhERKRGmY6/\nrzdsSQuWicQm1gaBmfU3syFmtjfwReBXGh0QERGpTabj75ct6wlh6lq4cEM4UxEtWCbSfHGPELye\ne/0RuAf4aszXExERaWmtHH9f04RpEalb3CsVfwgYAuwJfAe4APjX4p3MbAYwA6Cjo4Ourq6aLtbd\n3V3zsWmg8ierWeUfPnw4q1evbug533nnnYafs5lU/r7Wrl2b6f+fpEZjxkQKjSkbf59RhROmvWDC\ndKs1ekTSKNYGgbv/Ofef95rZK8BvzWy2u/+jaL85wBwIFpqpdcGPzC00U0TlT1YzFyZr9MIkWV6o\nBVT+MIMGDWK33XZr6DklA+qNk09RetFqtdKEaZGsaeak4nzj4N1NvKaIxGDWrFmMHj26qdc0My64\n4IKe3zs7O5k2bVrk47u6ujAz/va3vwHw1ltvMWvWLP7yl79EPrb4NWDAhj6V6dOn99o+atQo9t57\nb37wgx+watWqKj6pyAZl4+kzkF60GpmeMC0lqb7IRn0Rd8hQoQ/nfj7ZxGuKSIv65S9/ycCBAyPv\n/4EPfIDFixez3XbbAcED/uyzz2b8+PG8//3vj3SOK664gm233bbndytKjbjzzjvzm9/8Bnfn1Vdf\n5b777uPHP/4xc+bM4e6772b8+PGRy5s1ZrY9wTyxScAE4I/u3lnhmPGE1wm/d/fjGlzEzMmH0Kwc\nBIvGtX7oTH5uRP4zj8jShGlJNdUXlcXSIDCz24A7gf8D3iFoDMwkeMj/o9yxIiJR7LLLLlXtP2zY\nMPbaa6+6rrnrrrvynve8p+T7Q4cO7XWNKVOmcOqpp7Lnnnty0kknsXDhwrqun3ITgEOB+4DoNW/g\n34H/V/D7K40qVFa1azx9vlEw88BgZKCVP6s0T1rqi8K5ZmmrL+IKGfoTMB24BrgaOBz4JvCpmK4n\nIgnKD5HeddddHHnkkQwdOpQddtiBBQsW8M477/DVr36V0aNHM3bsWH7yk94Llk+fPp19992XP/zh\nD+y8884MGjSIvffem0ceeaTsNcOGgP/6179y+OGHM2LECDbZZBP23HNP7rjjjl5lzA8B52P+Tzrp\npJ6h26eeeqpBd2SDrbfemu985zt0dXWxZMmShp8/Rea7+zh3P4agM6gaj7n7fQWvx+MoYJaUi6dv\ndfkJ06GNAbO+rzFjml5GqZ3qi9KSrC9iaRC4+5nu/h5338TdR7j7B9z95+6+Lo7ribSzpcuX8vHr\nP87S5UuTLgqnnnoqe++9N9dffz3bbLMN06ZN4/Of/zyrV6/myiuvZNq0acycOZP777+/13H//Oc/\nOf300znzzDO58sorWbVqFQcddBBr166NfO0lS5bw4Q9/mBdeeIELL7yQ66+/nqlTp/LMM8+E7n/3\n3XcD8O1vf5vFixezePFittxyy7LXeOedd3j77bd7XuvXr49Uto9+9KMA3HfffZE/T9a4e7SbIZG0\nVTy9FixrGtUXAdUXfTVzDoGINNjS5UuZdMkkVq5dyaRLJrH4lMXsOGrHxMrzqU99iq9+NVhuZOut\nt2bChAk89thjPQ/TAw44gN///vfMmzePD37wgz3HLV++nBtvvJEPfehDAOy+++5st912XHbZZZx2\n2mmRrn322WczfPhw/vjHPzJ48GBgw4M1zB577AHAdtttF3louDh29IwzzuDcc8+teNzWW28NwIv6\n4lLKb8xsM+Al4CrgDHdfk3CZEtVW8fSlJkIXxVwvHaVQonqovtggqfri61//esXjkqov1CAQyaj8\nw33FmhU4zoo1KxJ/yO+///49/7399tsDsN9++/Vs69evH9tuuy3PPfdcr+M233zznoc7wDbbbMPu\nu+/OAw88EPkBf/fdd3PCCSf0PNzjMHfu3J5JZgBbbbVVpOO0QHtJbwK/ABYArwGdwNeB7YAjkytW\nOiiefoN2m2DdaKovelN90ZcaBCIZNfP2maxcuxIneHg4zsq1K5l5+0zmHz8/kTKNGDGi57832mij\nPtvy24uHdjfffPM+59piiy144YUXIl97+fLlFYdw6zVhwoSyk4pLyVdoHRnOER8Hd38B+HzBpi4z\nexH4pZm9z93/p/iYVl/IsrPo93ILkKWx/FFFuf+duZ9RJlg3814082+nUQtZfunmL4XWF1+6+Utc\nPfXqus9fzptvvom793yON954A4CBAwf2+WyDBg3qta1///50d3f3bFu3bh2jR4/uc9xmm23GM888\n02v72rVre37Ph+/kf1++fDmbbbZZyXubL+Prr7/O6tWr6e7u7nPOUvLHjh8/np122qnXe/lFLNet\nW8f69etDz/XEE08A0f/tG7WIpRoEIhk1+6DZLHp2UU+Pj2GMGDSC2QfNTrpoVXv55Zf7bHvppZeY\nMGFC5HOMGjWqqgZEMy1YEAR+T5o0KeGSZMK1wC+B3YE+DYKWX8iyoyNSTPxbI0ems/wRVXP/oyxY\n1sx70cy/nUYtZHn+lPN7jRDk64vzp5wf+6KQG2+8MWbWc50hQ4YAQZad4msPGjSo17b+/fszYMCA\nnm0DBw7klVde6XPcq6++yoQJE3ptLzxX8XlGjRrFq6++WvKzF5cxnzK0uHxRji2UX8Ry4MCB9OvX\nL/RcixYtAmDy5MmR/m0atYhlMxcmE5EG2nHUjiw+ZTEjB4+kn/Vj5OCRiceE1urll1/ueQhCMGns\nz3/+M3vuuWfkc+y///5cffXVkSeW5UcwqpmIVotnn32W7373u0yePLlPb5GE8qKfrWnMmPCMOcWN\ngY6O0AXIFs2bl0y5E9BWE6xjovqiN9UXfWmEQCTD8g/5L938Jc6fcn4mH+4Q9NaccMIJnHvuuQwe\nPJizzjqLLbbYgunTp0c+x1lnncUee+zBRz7yEWbOnMmoUaN4+OGHGTVqFCeffHKf/TfaaCPe/e53\nc/XVV/Oe97yHQYMGseuuu/Y8+Gvx+uuvc9999+HurFy5kkWLFnHhhRey6aab8pvf/Kbm87aZfG7A\nhxItRdz6Sc9VAAAgAElEQVSiThjURPT2mmAdI9UXG6i+6EsjBCIZt+OoHbl66tWZfbgDvOtd7+LH\nP/4xs2bN4rjjjmPTTTfl9ttvZ9CgQZHPsdNOO3HvvfcyevRoPvOZzzB16lSuvfZattlmm5LHXHjh\nhbzyyisccMAB7LHHHjz//PN1fY4lS5YwadIk9t57bz75yU9y9913M3PmTB5++OGy5WgFZjbEzKaZ\n2TRgLLB5/nczG5Lb53Ezu6TgmFlmNtvMjjKzA8zsHOCnwDx3/2synyQZS0fB4Z8Ifkpf+UbBoX9X\nY6Aeqi8Cqi9CuHuqXrvvvrvXauHChTUfmwYqf7KaVf5HHnmk4ed87bXXGn7OZjnxxBN9t912S7oY\ndYnj/pf7OwEe9BQ8rwtfwHiCMJ+w1/jcPk8BlxUccxzwILAKeAt4HDgH2DjKNTNdXxQEAD02Ct/s\na3i/7wQ/HxtVFCAUIvHy1ylS+UODpUq8mqiZ9171RW9Zry+aXVe4R68vFDIkIiJ1c/enAKuwz/ii\n3+cCc+MrVfpFyaDTtiJOsK57YTMRUciQiIhIUspl0Gl7y5ZFGx8otbCZiESmEQIRSdRll13WkDzb\nIlk0e0Gw0FZ+hEAZdCIYMyb6yIEaCy1F9UV8NEIgIiKSkPxk2ZFrod/64KfChSpQhiaRhlODQERE\nJEHKoFMfZWgSqZ9ChkQS4O49Kx+KFAsSQ0g72XH5hlV3Jbr8pOyVg4LQq1ZsUKm+kFIaWVdohECk\nyQYOHMiaNWuSLoak2Jo1axg4cGDSxZA4Rc2Moww6JRVmaFpfkKGplUYKVF9IOY2sK9QgEGmyLbbY\ngueee4433nhDPcHSi7vzxhtv8Nxzz7HFFlskXRyJkzLo1K0dMjSpvpAwcdQVChkSabJhw4YB8Pzz\nz7Nu3bqGnHPt2rVVrdKYNir/BgMHDqSjo6Pn70QyThlxYtMOGZpUX/SV5fKnua5Qg0AkAcOGDWvo\nF76uri522223hp2v2VR+aVnKiBOb/GTs/ByCES2aoUn1RW9ZLn+ay66QIRERkSZRRpzGUoYmkcbQ\nCIGIiEgTtENGnCQoQ5NI/TRCICIiErN2yIjTNMrQJNJwsTQIzOwYM7vRzJ4zs24ze8jMPhHHtURE\nRNKuHTLiNI0yNIk0XFwjBKcD3cBXgCOAhcCVZvaFmK4npYwZA2a1v8aMSfoTiIhk3uwFwaRXWx/8\n3ooZcRIXtb5TvSbSR1xzCA5391cKfr/bzLYiaCj8PKZrtreoqe2qpcwXIiJ1a5eMOIlSRieRmsUy\nQlDUGMh7GNgqjusJkR9wpTJclM18oR4WEZG6KSNOcymjk0h0zZxUPAlY2sTrtaZSQ6JFwh6E+Ult\nt+zQezJbqe1lqYdFRKRq+Yw4agzEq6Z6TaSNNaVBYGb7Ax8DZjfjei0twhfxsAdhqQwXt29XPvNF\nNSMHnZMna+RARKSQMuI0nTI6iVQv9nUIzGw8cCVwg7tfVmKfGcAMgI6ODrq6umq6Vnd3d83HpkGU\n8ncW/b50VJClYvaCoMep8EHoBQ/C970YnuHi01NLZ76YvaCGnNkvvpjZf4Ms//1kueyg8rcCM9se\n+CrBaPAE4I/u3hnhuOHAzwg6jfoBNwFfdPfW6ENXppumK5fRSesViISLtUFgZpsBtwJPA58stZ+7\nzwHmAEycONE7Oztrul5XVxe1HpsGPeWPOEE4bJGbUg9C92ASW76hkM9wcfn1cPzRfbd//oHwhkW+\nUVDcECmU1X+DLP/9ZLnsoPK3iAnAocB9wMAqjrsa2BH4DLAe+CHwB2CfRhdQ2sPsBUGdWFyvKaOT\nSGmxhQyZ2RCCnp6NgMPc/Y24rtVyqmgMFA+Jfv6B8NR2F90cfJkfuRb6rQ9+Lr4EDvpH+PYL9izd\nw6LYTBEJMd/dx7n7McD/RTnAzCYBBwInuvt17n49cAKwt5kdEGNZG08pL1MjP3m7uF7TvA2R0uJa\nmGwAcA2wA3Cwu78Ux3XaRVgcf6mRgAv2LP0gLJXhImx7qZzZhSMHJWMzVQmKtB13X1/DYYcAL7r7\nfxec5wHgydx72aGUl6mijE4i1YlrhOCXBEPH3wVGmdleBa+NY7pmSyrVG19ukZtyD8JSGS6Kt5fq\nYSk3clCWKkER6WtnYEnI9kdz72WWUl4mTxmdRKKLq0GQ/3p4PrC46LVlTNfMnqIh5s7Jk3ulEC2X\nKaHSkGgjHoTVjBzkYzNVCYpIFUYCK0O2r8i9l0kKq0yIMjqJ1CyWScXuPj6O87acCr3mlTIl5L+w\nl5rcG6qjo3zWi6I1DfINi8LfS622GTbJWT0zItJIac1K15n7WSrTW+HzsBHXzXpmq1jKP3du9H3r\nuLbufbKyXP40lz32tKMSXXHmniiZEoq/sPfiXn0hOjoqNlTCGiJRKsGwBdQqNlBEpJWtADYP2T4y\n914fac9KFyXlZSOum/XMVlkuf5bLDip/ktJc9mauVCxlhA0x15UpodYh0WXLgoZE8atIcUhSuUqw\nLM0tEGlnSwifK1BqbkHqVQqrFBFJIzUIUiDKXIGSmRLCvry7N77XvUIDQ3MLRKQGtwJjzGzv/AYz\nmwhsm3svc5TyMqWUFlakLDUImqHUgyinUu96KjIlVBg5KFcJaoKdSOszsyFmNs3MpgFjgc3zv+fW\npcHMHjezS/LHuPtiYAFwuZkdZWYfA64A7nX3O5P4HI2glJcppLSwImWpQdAMFR4wNQ8xpyxTQlgl\nWG70Q0RayhYE689cA+wF7FLw+xa5fQYA/YuOOxa4B7gUuBx4CJjahPLGKhUdOVKSRq1FelODIAHF\nD6JIQ8zNCAtqgIbNLRCRTHH3p9zdSryeyu0z3t2nFx230t1PcvcR7j7M3Y9391eS+Ax1UcrLzNCo\ntUhfahA0WakHUWaHmOucWxBsVCyniGRcqbDKDHTktBONWouEU4OgiSo9iDI5xFxQCXYtXFjV3IKK\nFMspIiINpFFrkXBqEDRRzQ+ijA8xVxr9UCyniIg0g9LCioRTg6CRKmQTivQgatEh5lKjH4rlFBGR\nZlFaWJFwahA0UsQVftviQRRhVEOxnCKSKcpl3xIyO2dPJEZqEMQoLBSmbR5E5SbY5SiWU0QyRbns\nW0Ym5+yJxEgNgpiUC4XRgyhQcwYi9b6JSApo/lOGKC2sSFlqEMSg1lCYt0aObE4BU6LmECr1volI\nwjT/KWOUFlakLDUIYhApFCbkQbRo3rxEypukciFU6n0TkTTS/CcRaTVqENSiEdmEpEdYCJV630Qk\nrTT/SURajRoEtVA2ofpUiNFU75uIpJk6fUSk1ahB0ABtnU2oFqViOXPU+yYiaaZOnxaiVLIigBoE\ndVM2ocZT75uIpJ06fVqEUsmKAGoQ1KXm0BalNSsrUu+benBEJGHq9Gk9SmYh7UoNgjrUmk1Iac0q\nq6n3TT04IhIn5bJvaUpmIe0stgaBmW1vZheZ2V/N7B0z64rrWklRaEu8SvW+qQdHRBKhXPYtS8ks\npN3FOUIwATgUeAxYGuN1EqOJZQ0WoVdNPTgiItJoSmYh7S7OBsF8dx/n7scA/xfjdeJTYb0B0MSy\nhqqQfUg9OCIiEgeN+Eu7i61B4O7r4zp300SMSdfEsuZQD45IepnZLmZ2l5m9YWbPm9k5Zta/wjHj\nzcxDXnObVW4R0Ii/iCYVV6Hq2HVNLGso9eCIpJOZjQTuBBw4EjgHmAmcHfEU/w5MKnh9O4ZiRqfc\n9G1JI/7SzgYkXYCsyIerrBwEi8aFPCwKQlskHvmHdf7fYUSpdKTFOjo0yU8kXqcBg4Gj3P014A4z\nGwbMMrMf5baV85i73xd7KaNSbvq2lR/xF2k3qWgQmNkMYAZAR0cHXV1dNZ2nu7u75mPDdOZ+Fsau\ne0HseuGX0UZct9Hlb7a4yt9Z8N/5RsHMA4ORgajpSKOUK8v3P8tlB5W/BRwC3F70xX8u8ENgX2B+\nIqVqkKWjqnzmSHZ0dERr2GnEX1pcKhoE7j4HmAMwceJE7+zsrOk8XV1d1HpsOeVi1/M9CY24blzl\nb5Zmlb9UD065SjtKubJ8/7NcdlD5W8DOwN2FG9z9n2b2Ru69Sg2C35jZZsBLwFXAGe6+JpaSVqni\n6LBkm0aPRQDNIYhEsesJUzpSkbQbCawM2b4i914pbwK/AE4B9gcuAj5HMLqQOGU2E5F2kYoRgsSN\nGVN2yDBS7LrEp1QPTm6+QJSQLhFJH3d/Afh8waYuM3sR+KWZvc/d/6f4mGaEmHbmfkYZHU4qVCzr\nYWpZLn+Wyw4qf5LSXPbYGgRmNoRgYTKAscAwM5uW+/0Wd38jrmtXLUL8YE2x69IUUSptEYnVCmB4\nyPaRufeqcS3wS2B3oE+DoJkhprMXBGFC+c6GsNHhpELFsh6mluXyZ7nsoPInKc1ljzNkaAvgmtxr\nL2CXgt+3iPG6dSuVXlTrDaSTQrpEEreEYK5ADzMbBwzJvVcNL/qZGOWmb2NlUs92Tp6s1LPScuJc\nmOwpd7cSr6fium69aopFV/aBRKnSFkncrcBBZrZpwbZjgTXAPVWeKz+S/FAjClYv5aZvU0o9K21G\nk4oLRJpA5t73pSwFiatYaWtRIZE4XUgwQXiemR2Qi/OfBfykMBWpmT1uZpcU/D7LzGab2VG5484B\nfgrMc/e/NvkzlKTRYal6YVKRjFGDoEC5WHRJv6orbfXsiDSEu68gyBLUnyDF6NkEX+zPKtp1QG6f\nvCUE6xT8BrgFOB44L/czOVFHfTU63BaUxU7agbIMFYgygUxSJOKCMlpUSCR+7v4IsF+FfcYX/T6X\nlKQY7UWjvpKjLHbSLjRCUECx6BmzbFl4CFcB9eyIiEitFDkg7UINgiKaQNY6tKiQiIjUQ1nspF20\nV4OgVBqxIppA1hrUsyMiIvVQ5IC0i/ZqENQ7iVQTyDJFPTsiUlaZXPPKSCZ5ihyQdtBeDYIiZdOI\nKb1o5qlnR0TKUq55iUiRA9Lq2rZBoMmm7UE9OyISlXLNSw+lnpU205YNAk02bS9le3bClqNXeIBI\n21EnkfRSKoudO10LFypyQFpOWzYINNm0xdXbY6PwAJG2ok4iEWl3bdkg0GTTFhdhfQJQeICIBNRJ\nJCLtri0bBJpsKgoPEJE8dRKJSLtrzQZBhPUGNNm0fSk8QEQKqZNIRNpdazYIIsaAK41Ye1J4gIgU\nUyeRRDJmzIYEFFq7QlpIazYIilQdK640Yi1N4QEiEkadRFKR1q6QFtXyDYKKseJagKztKDxARADl\nmpe6KTmFtIqWbhAoVlxKUXiAiJTLNa9OIqlEySmklbR0g0Cx4lJO1AXLFBMqIiKF1OEoraalGwSK\nFZdetGCZiIg0gDocpdW0dINAseLSixYsExGRBlCHo7Sa2BoEZraLmd1lZm+Y2fNmdo6Z9Y/reqUo\nVlyqoZhQkdrU+sw3s+Fm9hszW2Fmq8zsCjPT/3mSaupwlFYTS4PAzEYCdwIOHAmcA8wEzm7YRUIW\nH+vJDVxEqeQkCsWESmxKLZbYIvNU6nzmXw10Ap8BpgN7AH9oaAHL1RctcP8lGepwlFgkVF8MaOjZ\nNjgNGAwc5e6vAXeY2TBglpn9KLetPo2I51YqOSlQLiZ0/lXJlk0yrvVzl9f0zDezScCBwL7u/t+5\nbc8B95vZAe5+Z0NK1/r3XxKS73AUaZiEnldxhQwdAtxeVAnMJagw9o3jgmXjvpVKTiJQTKg0SwvO\nU6n1mX8I8GK+MQDg7g8AT+bei0UL3n9pFq1dIU3WrOdVXA2CnYElhRvc/Z/AG7n3Gkpx39IIigmV\nZmjR51Wtz/w+x+U8WuG4mrXo/ZdmWbaMroULtXaFNEUzn1dxNQhGAitDtq/IvdcwivuWRqoYE9pC\ncd/SfC38vKr1md+0ugJa+v6LSItp9vMqrjkEVTGzGcAMgI6ODrq6uioe05n7GSXuO8r50qC7uzsz\nZQ2T1fJ3Fv1edUzoiy8m/rmzeu/zWr38nbmfrfS8Sorqi0Cr/z+TZlkuO6j8SYpS9s7cz2Y/r8xD\n8rDXfVKzl4BfuPvZRdtfB2a5+3mljp04caI/+OCDUS4C9G5Beb8g7rtPqEcMnzEOXV1ddHZ2Jl2M\nmmW2/GPGRJqcs3RU8D/i7AUhIwcJ/41l9t7ntHz5G/y8MrOH3H1i/SVvjFqf+WZ2NbC5u08u2n4z\ngLtPKXdd1RedSRejZlkuf5bLDip/kiKVPaH6Iq6QoSUUxX+a2ThgCOHxojVT3LfULbdgWZ+40AKK\nO5ZGaOHnVa3P/D7H5ZSaW1CXFr7/ItJimv28iqtBcCtwkJltWrDtWGANcE+jL6ZcwBInxR1LI7Xo\n86rWZ/6twBgz2zu/wcwmAtvm3mu4Fr3/krQWX2tEktHM51VcDYILgTeBeWZ2QC7mcxbwk4asQRBC\ni49JXMrF8Yn0UupLQZEWfF5Feuab2eNmdkn+d3dfDCwALjezo8zsY8AVwL0NW4MgRAvef0ma1rqQ\nmDTreRVLg8DdVwD7A/2B+QSrVf4UOKthF1EuYGkSrU8gkdVb2Wf0eVXFM39Abp9CxxKMIlwKXA48\nBExtaAFVX0iTaa0LqVlCz6vYsgy5+yPAfnGdPyzHb5Ynmkh65YfsJp0SjAyMUNyxRJTmieiNFuWZ\n7+7jQ7atBE7KveKh+kKaKB9munIQLBqn+kKqlNAaFnGFDIm0FMUdS7U0EV2k/WjOmWSVGgQiEZWN\n49PkMSmgLwUi7UlzziSr1CAQKaXe+DxNHmtb+lIg0p4050yySg0CkVJy6xP0eRXR5DEppi8FIu1J\na11IVqlBIFIHxYlLGH0pEGlfmnMmkaRs7Qo1CERqpDjxNlb0IO+cPLnPegP6UiDSvrTWhVSUsrUr\n1CAQqZHixNtYxAe0vhSItAmtdSF1Sjr8WA0CkRopTlzyqn6Q60uBSGspNees+JVQjnlJtzSEH6tB\nIFIjxYkLRHiQ60uBiIiUkJbwYzUIROpQMU484UlCEq+0PMhFRCSb0hJ+rAaBSJ1qihPXGgUtIS0P\nchERyaa0hB+rQSBSrSriv5OeJCTxSsuDXEREsikt4cdqEIhUq9zksQJpmCQk8UrLg1xEUiplueYl\nndKQploNApEYKLa8RZSqzAuk4UEuIimVslzzkl5Jp6lWg0AkBootbxFab0BEGkhhpNIjZWtXqEEg\nEgPFlrcmrTcgIrVSGKn0krK1K9QgEImBYstbj9YbEJFaKYxU0k4NApGY1LRGgSaWpZIqcxGph8JI\nJe3UIBCJUdWx5ZpYlkqqzCszs8+a2d/NbK2ZPWRm+0c4ZpaZecjr4GaUWaRZFEYqaacGgUgjRYwZ\n18SybFFlXp6ZfQK4ELgcOAT4P+AmM3tPhMNXAZOKXotjKqpIIhRG2sYKstV1Tp6c2tSzahCINFKp\nSUIFNLEse1SZVzQL+K27f9fdFwLTgceBb0Q49m13v6/otSrGsookQimK21RGUs/G0iAws2PNbJ6Z\nvZAb/p0ex3VEskax6Cml9QZqZmbbAjsCV+e3uft64BqC0QIRyVGKYklrhEBcIwTTgPHATTGdXyST\nFIueUlpvoB47534uKdr+KLCZmW1e4fgRZvaKma0zs4fN7KjGF1EkISnLNS/JSnOEQFwNgmPd/QPA\nV2I6v0gmKRY9G6rtwXlr5Mh4C5Ru+Q+/smj7iqL3wzwOfA04BjgaeB64To0CaRkpyzUvyUl7hMCA\nOE6aGy4WkSL5sJNJpwQjAyMUi546+Yf2ykGwaFzIv0/RnBCARV1ddDathPEzs+HAlpX2c/fiUYGq\nuPvviq47H1gEfAeYV6JsM4AZAB0dHXR1ddV07e7u7pqPTQOVPzlZLjuo/M3WmftZLkJg/lXBtiQ/\nVywNAhEpLd8omHlgMDIQuj5BsY4O9SA1QWEPjhf04LRho+0Y4OII+xkbRgKG03uUID8ysIKI3N3N\nbB7wQzPr7+7vhOwzB5gDMHHiRO/s7Ix6+l66urqo9dg0UPmTk+Wyg8qflNkLgk6mfP0SFiGQ5OdS\nliGRBGh9gnTSHI+Au//a3a3SK7d7fpRg56LT7Ay86u4vV3v53EtEpGWkPVtdpBGCuIePNQQcUPmT\nFWf5PzRyJButqNxRunRU6ZGDcmXTva9PZ+5nlB6csHImXf4kufsTZraUYFThdgAz65f7/dZqzmVm\nRjCX4H/CRgdEWsaYMdE6ejQ63FIqRggkKGrIUDXDx1XTEHBA5U9WrOV/9dXw7QXhQZVi18uVTfe+\nMaLM8QgrZ1rKn6BZwO/M7Cng/wEnAjsAx+d3MLN9gbuA/d39nty2e4DrCEYZhgKfBT4IfKyJZRdp\nvozkppfGy0cIpE2kkKEqh49FpEppzz7QMrTeQCzc/SrgNIIFyW4DdgUOc/e/FexmQH96dxw9DnwZ\nuAH4HbApMMXdb2xCsUVSI6256aUBMpJ6VpOKRVIgSvYBaYAq1xuQ6Nz9YsqMJLt7F0WjyO5+SszF\nEkm9ipnNJNsKQr7SPJoc10rFu5jZNDYM+040s2m5IWMRKaL1CZJRda+cFg8SkQbS6LCkRVxZhj5O\nsGz9f+V+/7fc72fHdD2RTIuUfSAs1GXMmMTKnHUVV4zU4kEiEjNlNpO0iKVB4O6zSswz6IzjeiKt\noKbYdU04q4l65UQkDTQ6LGmhdQhEUqTU+gSacNZY6pUTkTRIe256qUKppBUFr87Jk1M7sq8GgUiS\nIsSkVwxtkd7KPZRz1CsnImmhzGYtIuOpZNUgEEnSsmXhseo5Cm2pQYSHrXrlRCRNql69XlIvayP7\nahCIpJhCW+pX6qGsXjkRSUxGctNLbbI4sq8GgUiKRQptycclKvtQH5UeyuqVE5FElBodVmazzMvq\nyL4aBCIpVnNoS0pjFJuproeyeuVERKQGWR3ZV4NAJOXKhbZkLUaxmSI9lNUrJyIiDZTVpBVqEIhk\nQFhoSxZjFBuqVDahnKw+lEVEJLuymrRCDQKRNKoQspLVGMWGqhAWldWHsoi0sQi57DVPLP2ymLRC\nDQKRNKqQjjSrMYpxCgufyuJDWUTaWMZz2csGWUtaoQaBSAZFzT7ULr1K5cKnsvZQFhHJ0zyxDMl4\nKlk1CEQySNmHNqg5fCqlD2UREdA8scyJkEq2a+HC1CatUINAJKPaJvtQhcnDNWcTSulDWURE88Sk\n2dQgEMmwtsg+VGFUQ9mERKTVaJ6YNJsaBCJZ0ojsQxmfW1A8+qFsQiLSatTRIc2mBoFIlsSVfSgj\ncwtKjX4om1CyzOxYM5tnZi+YmZvZ9CqO/bCZ3W9ma83sSTP7YoxFFckEdXSkVAunhVWDQKSFVOpV\nSvXcggpzBSqNfiibUKKmAeOBm6o5yMy2B24HngQOBS4CfmJmn2l0AUWyRh0dKdTCaWHVIBBpIeV6\nlVI/t6DCA7Tm0Q9lE2qGY939A8BXqjzuq8DzwAnufre7/wCYA5xlVtAaFGlT6uhIt1R3slVJDQKR\nFhPWq5TFuQXFD9pIMbXKJpQId19f46GHAPPc/e2CbXOBrYH31F0wkazJeC77dpL6TrYqqUEg0oKK\ne5VSNbcgJDSoc/LkXruEPWgVU9tazGwoMA5YUvTWo7mfOze3RCIpECGXvTo6kteKaWHVIBBpAW+N\nHFn2/brmFjR60lSFRka5B61ialvKiNzPlUXbV+R+lv+jFhFJSCumhR3Q6BOa2TBgJsFQ8E7AGmAx\n8HV3X9ro64kILJo3j87Ozr5v5MKw81+kJ50SPLRGhMwtWDkIFo2r4ot2pdGDMWMijTAsHRU8RGcv\nCK5b7kGbH/WYf1WE8knNzGw4sGWl/dy9uHc/dmY2A5gB0NHRQVdXV03n6e7urvnYNFD5k5NU2T90\n1FFstGJFxf3eGjmSRfPmlXw/y/ceki1/Z+7n7AVBfbkiV1eFhbCGlTHN977hDQLgXcBngUuAM4Ah\nwDeB+81sV3d/JoZrikgF+UZB4Zfvwt54L+iNL2wUFH9h76XOeZ9hjZEoD9pQiqltpGOAiyPsV88f\nQH5kYHjR9vzIQOg3H3efQzDxmIkTJ3poQziCrq6u8EZ0Rqj8yUms7BEaAwAbrVhRtnxZvveQjvKX\n62TLCytjGspeShwhQ08C27n7me5+h7vfQJBObiBwcgzXE5GIqp1b0MhJU8VhSaVCgyDCXAHF1MbK\n3X/t7lbpVec1Xgeeoe9cgfzvTR99EMmSVspwk0WtFsLa8AaBu7/u7muKtr0KPA1s1ejriUgZFXrN\ny80tqDRpqlRlFLY9rGFRrjHSag9aKelWYKqZ9S/YdixBQ+FvyRRJJP1aLcNNVrVSWtimTCo2s82B\n7QHNIRBppgorG5fL3FPuC3upyihse6mGxecfKD/RuZUetK3OzHYxs2nAx3KbJprZNDPbt2Cffc3s\n7cJtwHkEKUb/y8wmm9nXgFOBc9wL/lBFpEcrZrjJjBZOC9usLEOzgW7gsiZdT0QiKtUbX2r04PMP\nhFdGt28Xvv20w8IbFhfsWWMa0Qw+aNvAx4FrgP/K/f5vud/PLtjHgP4UzDtw98eBgwk6jG4F/hWY\n6e6/bkKZRTKpFTPcZEYLp4WNNKm4nowTZvY54ATgaHcPreqVNSKg8icry+WvtuydRb+HZe4pNWmq\nVGX06anh292DY8MmCYdNdC7UtXBh+AdI2b9Tlv92GsHdZwGzKuzTRcgkZHe/F9gzjnKJtKKaEy+I\nlBE1y1BNGSfM7Ajg5wQpR68vdZCyRgRU/mRlufxVl72jI1JK0LAv7KUqo8uvh+OP7rv9opuDc5XK\nxlAujWhW/j2y/LcjItkSJcON1Cli2mw6OjI5GhAmUshQLRknzOzDBEvQX+ju58VReBGpUblhzyLF\nsfyl5h0c9I/SIUA1TRJWaJCISCglXohZlMZANftlQCxzCMxsAjAfuA34YhzXEJGYRPgiXqoyKldJ\nlZsk3LVwYUvEYIqINIsSLzRPO6R4jWOl4i0IGgLdwH8Ce9qGxYtec/dHGn1NEWmgUl/EixYhKxXq\nU5vTK48AABDfSURBVPVKwhoJEBGpLGKop56pjRW2gGYrNsLiGCHYhSCN3DhgIbC44PXLGK4nIs1Q\nbyXT0dEy2RhERJquhTPcpFU7pXht+AhBqUwSIpJxqmRERKSNlEvxWtVIeAY0ax0CEREREWm0MWOC\nkM5KrzFjki5p5pRaj6cVU7yqQSAiIiKSVW2YEadZSmXV0xwCEREREUmtdsiI00ztkuK14XMIRERE\nRKT52iUjTrNVnT0vgzRCICIiIpJx7ZQRJ3ZRs+q1UIpXNQhEREREMq5cRhypUhumeFWDQERERCTj\n2ikjTsMoQ1MPNQhEREREMq6dMuI0jDI09VCDQERERKQFtEtGnLi0c4YmZRkSERERaRHtkBEnDu2e\noUkjBCIiIiJZ1YYZcRpNGZrUIBARERHJrlIZcYobAC++CGZ0Tp7cdhNmK1GGJjUIRESkTmZ2rJnN\nM7MXzMzNbHrE42bl9i9+HRxzkUVanybMRqYMTWoQiIhI/aYB44Gbajh2FTCp6LW4YSUTEaC9J8xW\nogxNmlQsIiL1O9bd15vZJsBnqjz2bXe/L45CiUig3SfMRpFvFMw8MBgZaLf7oxECERGpi7uvT7oM\nIhJOE2ajy2doarfGAKhBICIiyRphZq+Y2Toze9jMjkq6QCKtRBNmy1CGph5qEIiISFIeB74GHAMc\nDTwPXKdGgUjjaMIsQSalwsxK+VfxhOqOjvCMTcuWJVPuJtIcAhER6cXMhgNbVtrP3ZfUcx13/13R\ndecDi4DvAPNKlG0GMAOgo6ODrq6umq7d3d1d87FpoPInJytl78z9zMfG5+cQjAiZMJuFz5NXy/3v\nrCLjUpz3Is1/O2oQiIhIsWOAiyPsZ428qLu7mc0Dfmhm/d39nZB95gBzACZOnOidnZ01Xaurq4ta\nj00DlT85WSx7pQmzWfo8jbj/S0clcy/S/LejkCEREenF3X/t7lbpFdflcy8RaaB2njBbKD/J+pYd\nNLm6UCwNAjO70MyWmFm3ma0ws/82swPiuJaIiLQGMzOCuQT/EzY6ICJVqGYibFh8fQuuYKyMS6XF\nFTI0GLgAeAzYCDgFuNXM9lG+aRGR1mJmuwC7AINymyaaWTfwsrvfk9tnX+AuYP+CbfcA1wFLgKHA\nZ4EPAh9r7icQaUGlJsJaxMG9FlzBuFzGpflXJVu2pMXSIHD3Ewt/N7PbgCeB4wE1CEREWsvHgbMK\nfv+33OseNsxtNKA/vecdPA58mWAC83rgz8AUd7815vKKSE65ePpWM3tBsDDbilyjoC0zLpXQlEnF\n7v6Oma0kGC0QEZEW4u6zgFkV9umiaBKyu58SW6FEpKJ2W8E4SsaldhXbpGILDDCzUWb2FWAH4NK4\nriciIiIi0bRrPH2+UXDo39UYKBRnlqFjgXXAK8B3gWPd/YEYryciIiIiEbTkCsalFiArooxLfUUK\nGapxkZrbgT2A0cAngblmdmhu2Lj4/FpoBpU/aVkuf5bLDiq/iEiztWQ8fb0ToavJzNRios4hqHqR\nGndfATyY+/U2M9sKOAf4SPFBWmgmoPInK8vlz3LZQeUXEWm2doinLzth2rXcSaFIIUMNWqTmYWDb\n+ossIiIiIvWqGE+f4fUJtABZdZqyUnFusZlJBKlHRURERCQFqo6nz8D6BO06YboeDW8QmNk+Zna9\nmX3azDrN7GjgBmAv4PuNvp6IiIiIRPPWyJGR9ls6Cg7/RDa/RLfkhOmYxTFC8AzwDvA9gonF5+eu\ns7e73xLD9UREREQkgkXz5gXx88WvAlkPt5m9IJgTYeuD31tiwnTMGt4gcPen3H2au2/t7hvnfh7m\n7osbfS0RERERaZxMhNtUSC+anxsxci30Wx/8bLUJ043WlDkEIiIiIpJ+mQi3iTCPQQuQVUcNAhER\nEREBshluU2q+gxYgi04NAhEREREBIobbpCgdaU3zHdp4AbJS1CAQERERkR41hdskkI400nyHsAnU\ny5Y1vaxppwaBiIiIiPRSKtymqelIQyYPd06e3PN2JuY7ZIQaBCIiIiLtLkIYTdPTkVYYdcjifIe0\nUoNAREREpN0tW1Z2fYJI4Tkxzy0oHp1QetHGUYNARERERMqqOTynQXMLSo1OKL1oY6hBICIiIiJl\nVQrPqWtuQYWFxiqNTii9aP3UIBARkZqZ2TAzO9vMHjCzVWa2zMyuN7MdIx6/i5ndZWZvmNnzZnaO\nmfWPu9wiUp1y4TkV5xaEfdkvfFUYRah5dELpRSNTg0BEROrxLuCzwO3ANOBUYEvgfjMbV+5AMxsJ\n3Ak4cCRwDjATODvOAotIbcLCcyLNLahS8WhDpMnDSi9aFzUIRESkHk8C27n7me5+h7vfABwKDARO\nrnDsacBg4KjcsRcSNAZON7NhsZZaRGpSHJ5Tqfe+VChRue3Fow2aPBw/NQhERKRm7v66u68p2vYq\n8DSwVYXDDwFud/fXCrbNJWgk7NvQgopIbSqE3ZTrvS8VSlRpe9hogyYPx0sNAhERaSgz2xzYHlha\nYdedgSWFG9z9n8AbufdEJGkV0pGW6r2H8C/3t29X+kt/pdEGTR6OjxoEIiLSaLOBbuCyCvuNBFaG\nbF+Re09EMiCs977Ul/tPTy39pb/mhcY0ebhuA5IugIiIpIuZDSeYGFyWuy8p3mZmnwNOAI5294b3\n45nZDGAGQEdHB11dXTWdp7u7u+Zj00DlT06Wyw6NK39n0e/53vu82Qtg0bhgBMD7bfhyf/n1cPzR\nfbfPXrChYTHplKCRMCJkrsBN8+ezySab9C1QBv5N0vy3owaBiIgUOwa4OMJ+1usXsyOAnwNfd/fr\nIxy/Ahgesn1k7r0+3H0OMAdg4sSJ3tnZGeEyfXV1dVHrsWmg8icny2WHBpa/o6NsutByX+7LfenP\nv58fMSgOD9pkk00ye//T/LejkCEREenF3X/t7lbpVXiMmX2YYELwhe5+XsRLLaForkAuVekQiuYW\niEjKVJhbAKUnAleaIKy5As2nEQIREamLmU0A5gO3AV+s4tBbga+a2abuvjq37VhgDXBPY0spIkko\nDiWqtL0szRWIjUYIRESkZma2BUFDoBv4T2BPM9sr99qlYL9tzOxtM/t0weEXAm8C88zsgNz8gFnA\nT4pSkYpIVtT7pb2jI3zkQQuNxUojBCIiUo9dgK1z/72w6L172DD30ID+FHREufsKM9sfuIBghGEl\n8FOCRoGIZJG+tGdS7A0CM/sS8DPgOnefFvf1RESkedy9i6LJxSX2eypsP3d/BNiv4QUTEZHIYg0Z\nyg0lzwJejvM6IiIiIiJSm7hHCL4P3ASMi/k6IiIiIiJSg9hGCMxsT+DjwDfiuoaIiIiIiNQnlgaB\nmRnB4jQ/cvfn4riGiIiIiIjUL64RgpOADuDHMZ1fREREREQaINIcAjMbDmxZaT93X5Lb9/vAF9x9\nTcTzzwBm5H7tNrPHohwXYjTwSo3HpoHKn6wslz/LZQeVv1rbNPFaqfXQQw+9YmZP13i4/uaSleXy\nZ7nsoPInKYmyR6ovzIuWmQ7dyewzwMWV9nN3M7MfEqSQ+2jBWzcBy4ETgdXu/k6UwlXLzB5094lx\nnLsZVP5kZbn8WS47qPzSfFn/N1P5k5PlsoPKn6Q0lz3SCIG7/xr4dcRz7gRMBFaEvLcC2Ae4N+K5\nREREREQkRnGkHf02wUJkhX4GrALOAv43hmuKiIiIiEgNGt4gcPe/FW8zs5XAK7kVLeM0J+bzx03l\nT1aWy5/lsoPKL82X9X8zlT85WS47qPxJSm3ZI80hqPsiZl0EDYJpsV9MREREREQia0qDQERERERE\n0im2lYqbycx2MbO7zOwNM3vezM4xs/5JlysKM9vezC4ys7+a2Tu50ZRMMLNjzOxGM3vOzLrN7CEz\n+0TS5YrKzKaZ2SIzW25ma83sMTP7tpltlHTZqmVmY3P/Bm5mmyRdnkrMbHqurMWv05IuW1RmNsDM\nvmFmfzezN83sWTP7adLlktKyXFeA6oskqb5IjuqL5ohjUnFTmdlI4E7gEeBIYDtgNkFj59sJFi2q\nCcChwH3AwITLUq3TgSeBrxDk1T0UuNLMRrv7zxMtWTSjgLuB84CVwJ7ALGAM8PnkilWT84BuYGjS\nBanSfkDheiVPJFWQGlxGUP6zgSXAOGCXJAskpbVAXQGqL5Kk+iJ5qi9ilPmQITP7/+3db4hUVRzG\n8e+TQbaw2WZEFJIGQdCLXpYW2z8Ig6BeSAZFSARFgb0IDaGCCgWDCstKgugfSJFZorgtVCQtFAWJ\nQWaClREIoWxsZWTarxfnjtg0rjPuOueemecDy505l4Hfi53zzLn3nHNXAMuBiyJiompbTvVFbbTV\nlaTTIuKf6vUG4NyIuDZvVe2pOvL9TW3rgfkRMS9TWVMiaSXwADAUhXw5JA0D7wOrSB39YET8nreq\nyUlaArxKAbW2ImkhsBm4PCJ25q7HTqz0rADnRd04L7rDedEdvTBl6CZgtKkzfws4E7gmT0nta3Tu\nJWru3CvbgQu6Xcs0OgAUcwu4mu7wPPAE5T65sUR3Ax/XuXO3/yk6K8B5UUPOC2tHEXnRCwOCS0m3\nX46KiJ+Ag9U56675wO7cRXRC0gxJA5KuBpYCL5VytQe4DzgDeCF3ISdpj6TD1Xzce3MX04ErgN2S\n1kqaqOakb5RU8o+bXuesqB/nRXc5L/IoIi+KX0MADJHm8zUbr85Zl0i6AbiVNBouyR+kThLgDWBZ\nxlraJmk28CRwZ0T8LSl3SZ3YBzwKfAHMAG4H1kkaiIhaLbQ6jvOBJcAOUu2DwFPAe5KuLOgHQj9x\nVtSI86K7nBdZFZEXvTAgsBqQNBdYD2yKiNeyFtO5BcAAaZHYY8Ba4P6sFbVnJfB5RGzNXUinImIU\nGD2maUTSTOARSWsKmBqh6u+WiDgAIGkfsI20cOyjjLWZ1ZrzIgvnRT5F5EUvDAjGgVkt2oeqc3aK\nSToHGAH2AndkLqdjEfFV9XJM0n7gdUlPR8SenHVNRtJlpCtrw5LOrpoHquMsSUci4s/Wn66tDcBt\nwFzqv3vEOPB9o3OvjAGHSDtH1KKDt/9wVtSA86L7nBfZFZEXvTAg2EXT/E9Jc0j/7LtafsKmjaQB\nYAtpYdXNEXEwc0lT1ejs5wG17eCBS0jbDn7W4tzPwCvAPV2taOqi6Vhn3wIzW7QLqPvVqn7lrMjM\neZGN8yKvIvKiFwYEI8AySYMR8VvVtpi0V+22fGX1PkmnA++QOpsFEfFL5pKmw1XV8YesVZzYGHBd\nU9tC4GHS/t51v2LSyiLSzhd7cxfShi3A401bKQ6TQndHvrJsEs6KjJwXWTkv8ioiL3phQLCOtNJ/\no6TVwMWkfaWfKWRf6QHSFxLgQuAsSYuq91trfgXlRVLtDwKzq0VLDdsj4q88ZbVH0gekBxV9Axwh\nde4PAW/X+fYvHN3C75Nj26p5uQCf1n2vZknvkhaIfU1aJLa4+ltawHxQgJdJ/c5mSatIi8RWAx9G\nxFjWyux4is4KcF7k5LzIx3nRHcUPCCJivNqtYC3pwQ+/As+SOvoSnEe6anKsxvt5wI9draYzN1bH\nNS3O1b12gC9JK//nAodJV0lWkH442Kn1HWlO6xzSbdOdwF0R8WbWqtoUEROSrgeeI+1lfwjYRHoK\nq9VQD2QFOC9ycl7k47zoguKfVGxmZmZmZievFx5MZmZmZmZmJ8kDAjMzMzOzPuYBgZmZmZlZH/OA\nwMzMzMysj3lAYGZmZmbWxzwgMDMzMzPrYx4QmJmZmZn1MQ8IzMzMzMz6mAcEZmZmZmZ97F9DdH0R\no8Wi6wAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f1872ed2630>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "i_widget = IntSlider(min=0, max=79, value=0)\n",
    "N_widget = IntSlider(min=2, max=200, value=50)\n",
    "M_widget = IntSlider(min=2, max=400, value=80)\n",
    "eps_widget = Dropdown(options=[1.,1e-1,1e-2,1e-3,1e-4,1e-5,1e-6,1e-7,\n",
    "                               1e-8,1e-9,1e-10,1e-11,1e-12,1e-13,1e-14,1e-15,1e-16,1e-17],\n",
    "                      value=1, description='eps')\n",
    "delta_widget = FloatSlider(min=0.4,max=2.0, step=0.1, value=1.0)\n",
    "fun_widget = Dropdown(options=[1,2,3],value=1, description='Function F')\n",
    "time_widget = FloatSlider(min=0.1,max=10,step=0.1,value=1,description=\"Max Time\")\n",
    "method_widget=Dropdown(options={\"Implicit FD\":1, \"Explicit FD\":2, \"FFT\":3},value=1)\n",
    "def update_i_values(*args):\n",
    "    i_widget.max = M_widget.value-1\n",
    "M_widget.observe(update_i_values, 'value')\n",
    "interact(comparison,\n",
    "         w=fixed(w2),\n",
    "         wp=fixed(wp2),\n",
    "         F=fixed(F),\n",
    "         i=i_widget,\n",
    "         eps=eps_widget,\n",
    "         T_max=time_widget,\n",
    "         N=N_widget,\n",
    "         M=M_widget,\n",
    "         method=method_widget,\n",
    "         delta=delta_widget,\n",
    "         num=fun_widget)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "In this case, it seems to be more differences between both approximations. The presence of a non-constant speed seems to affect the results of the approximations with high orders of magintude for $\\varepsilon$. The first function $F(\\theta)$ has the major differences between methods. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Conclusions\n",
    "\n",
    "In this work, a method that uses different scales in order to obtain an approximation asymtotic at $T = \\mathcal{O}(1)$ was studied and developed with numerical examples. \n",
    "\n",
    "It is important to consider the effect of defining the _slow scale_ $T$, since with this scale is possible to avoid secular terms that ocurs in the _fast scale_.\n",
    "\n",
    "When a numerical method using a unique _fast scale_ was compared with the multiscale, several differences appear as time is increasing. Also, using a slow scale is less expensive that fast scale when $\\varepsilon$ is too small and the effect of diffusion in time is too long.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## References\n",
    "[1] Jupyter Notebooks of SC1. Claudio Torres. https://github.com/tclaudioe/Scientific-Computing/tree/master/SC1.\n",
    "\n",
    "[2] Jupyter Notebooks of SC2. Claudio Torres. https://github.com/tclaudioe/Scientific-Computing/tree/master/SC2.\n",
    "\n",
    "[3] Jupyter Notebooks of SC5. Claudio Torres. https://github.com/tclaudioe/Scientific-Computing/tree/master/SC5.\n",
    "\n",
    "[4] Using observe to update values in a widget. http://ipywidgets.readthedocs.io/en/latest/examples/Using%20Interact.html#Arguments-that-are-dependent-on-each-other"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Pendiente**\n",
    "- poner datos informativos \n",
    "- algunos 3D\n",
    "- comparison en el mismo tiempo al cambiar el M"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}