Half_line_integrals.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    }
   ],
   "source": [
    "import ipywidgets\n",
    "import scipy.linalg as lin\n",
    "%pylab inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# %load cheb.py\n",
    "from numpy import *\n",
    "from numpy.fft import fft,ifft\n",
    "\n",
    "def cheb(y):\n",
    "    '''Chebyshev transform. Finds Chebyshev coefficients given y evaluated on\n",
    "    Chebyshev grid'''\n",
    "    N = len(y) - 1\n",
    "    yt = real(fft(r_[y,y[-2:0:-1]]))\n",
    "    yt = yt/(2*N)\n",
    "    yt = r_[yt[0],yt[1:N]+yt[-1:N:-1],yt[N]]\n",
    "\n",
    "    return yt\n",
    "\n",
    "def icheb(c):\n",
    "    '''Inverse Chebyshev transform. Evaluates Chebyshev series at the Chebyshev \n",
    "    grid points given Chebyshev coefficients.'''\n",
    "    N = len(c) - 1\n",
    "    y = r_[c[0],0.5*c[1:N],c[N],0.5*c[-1:N:-1]]\n",
    "    y = y*2*N\n",
    "    y = real(ifft(r_[y,y[-2:0:-1]]))[:N+1]\n",
    "\n",
    "    return y\n",
    "\n",
    "\n",
    "def Dcheb(y,interval):\n",
    "    '''Chebyshev derivative of y evaluated on Chebyshev grid in interval [a,b]'''\n",
    "    N = len(y) - 1\n",
    "    a,b = interval\n",
    "    x = 0.5*(b-a)*(cos(r_[0:N+1]*pi/N) + 1) + a\n",
    "\n",
    "    k = r_[0:N]\n",
    "\n",
    "    A = real(fft(r_[y,y[-2:0:-1]]))\n",
    "    yy = real(ifft(1j*r_[k,0,-k[-1:0:-1]]*A))\n",
    "\n",
    "    fact = 2.*(x-a)/(b-a)-1\n",
    "    fact = fact[1:-1]\n",
    "    yprime = -2./(b-a)*yy[1:N]/sqrt(1-fact**2)\n",
    "\n",
    "    A = A/(2*N)\n",
    "    A = r_[A[0],A[1:N]+A[-1:N:-1],A[N]]\n",
    "    k = r_[0:N+1]\n",
    "\n",
    "    yprime1 = sum(k**2*A)*2./(b-a)\n",
    "    yprimeN = sum((-1)**(k+1)*k**2*A)*2./(b-a)\n",
    "\n",
    "    return r_[yprime1,yprime,yprimeN]\n",
    "\n",
    "\n",
    "\n",
    "def regrid(y,M):\n",
    "    N = len(y) - 1\n",
    "    a = cheb(y)\n",
    "    if M==N:\n",
    "        return y\n",
    "    if M>N:\n",
    "        a = r_[a,zeros(M-N)]\n",
    "        return icheb(a)\n",
    "    if M<N:\n",
    "        a = a[:M+1]\n",
    "        return icheb(a)\n",
    "\n",
    "\n",
    "\n",
    "def clenshaw(x,c):\n",
    "    '''Clenshaw algorithm to evaluate Chebyshev series at x\n",
    "    assumes x is in [-1,1]'''\n",
    "    N = len(c) - 1\n",
    "    b = zeros(N+2)\n",
    "    b[-1] = 0\n",
    "    b[-2] = c[-1]\n",
    "    for r in r_[N-1:0:-1]:\n",
    "        b[r] = 2*x*b[r+1] - b[r+2] + c[r]\n",
    "    s = x*b[1] - b[2] + c[0]\n",
    "\n",
    "    return s\n",
    "\n",
    "\n",
    "def clenshaw2(x,c,change_grid = True):\n",
    "    '''Vectorized version of Clenshaw algorithm\n",
    "    Use this for Chebyshev polynomial evaluation'''\n",
    "    if change_grid:\n",
    "    \tif (min(x)!=-1) or (max(x)!=1):\n",
    "        \tx = 2*(x-min(x))/(max(x)-min(x)) - 1\n",
    "    N = len(c) - 1\n",
    "    b = zeros([N+2,len(x)])\n",
    "    b[-1,:] = 0\n",
    "    b[-2,:] = c[-1]\n",
    "    for r in r_[N-1:0:-1]:\n",
    "        b[r,:] = 2*x*b[r+1,:] - b[r+2,:] + c[r]\n",
    "    s = x*b[1,:] - b[2,:] + c[0]\n",
    "\n",
    "    return s\n",
    "\n",
    "\n",
    "\n",
    "def chebD(c,interval):\n",
    "    '''Finds derivative of Chebyshev series in spectral space\n",
    "    i.e. maps c_n--->d_n where c_n,d_n are Chebyshev coefficients\n",
    "    of f(x) and f'(x) in the interval [a,b].'''\n",
    "    N = len(c) - 1\n",
    "    a,b = interval\n",
    "    if (a!=-1.) or (b!=1.):\n",
    "        factor = 2./(b-a)\n",
    "    else: \n",
    "        factor = 1.\n",
    "\n",
    "    b = c*r_[0:N+1]\n",
    "\n",
    "\n",
    "    cp = zeros_like(b)\n",
    "\n",
    "    cp[0] = sum(b[1::2])\n",
    "\n",
    "    evens = b[2::2]\n",
    "    odds = b[1::2]\n",
    "\n",
    "    cp[1:N+1-(N%2):2] = 2*cumsum(evens[-1::-1])[-1::-1]\n",
    "    cp[2:N+1-((N+1)%2):2] = 2*cumsum(odds[-1::-1])[-2::-1]\n",
    "\n",
    "    cp = cp*factor\n",
    "\n",
    "    return cp\n",
    "\n",
    "\n",
    "def chebD_semiinf(c):\n",
    "    '''Finds the derivative of Chebyshev series in spectral space\n",
    "    i.e. maps c_n --> d_n where c_n, d_n are Chebyshev coefficients\n",
    "    of f(x) and f'(x) in the interval [0,oo)'''\n",
    "    '''To be used only for the positive half-line'''\n",
    "    \n",
    "    N = len(c) - 1\n",
    "    \n",
    "    b = c*r_[0:N+1]\n",
    "\n",
    "\n",
    "    cp = zeros_like(b)\n",
    "\n",
    "    cp[0] = sum(b[1::2])\n",
    "\n",
    "    evens = b[2::2]\n",
    "    odds = b[1::2]\n",
    "\n",
    "    cp[1:N+1-(N%2):2] = 2*cumsum(evens[-1::-1])[-1::-1]\n",
    "    cp[2:N+1-((N+1)%2):2] = 2*cumsum(odds[-1::-1])[-2::-1]\n",
    "\n",
    "    d0 = 3./4*cp[0] - cp[1]/2. + cp[2]/8.\n",
    "    d1 = -cp[0] + 7./8*cp[1] - cp[2]/2. + cp[3]/8.\n",
    "    d2 = cp[0]/4. - cp[1]/2. + 3./4*cp[2] - cp[3]/2. + cp[4]/8.\n",
    "    d3 = cp[1]/8. - cp[2]/2. + 3./4*cp[3] - cp[4]/2. + cp[5]/8.\n",
    "    \n",
    "    dn = [ cp[i-2]/8. - cp[i-1]/2. + 3./4*cp[i] - cp[i+1]/2. + cp[i+2]/8.  for i in range(4,N-1)]\n",
    "\n",
    "    dn1 = cp[N-1-2]/8. - cp[N-1-1]/2. + 3./4*cp[N-1] - cp[N-1+1]/2.\n",
    "    \n",
    "    dn2 = cp[N-2]/8. - cp[N-1]/2. + 3./4*cp[N]\n",
    "    \n",
    "    dn = r_[d0,d1,d2,d3,dn,dn1,dn2]\n",
    "\n",
    "    return dn\n",
    "\n",
    "\n",
    "def cheb2zD_semiinf(c):\n",
    "    '''Finds the Chebyshev coefficients of the operator 2z df/dz when\n",
    "    f has a series in Chebyshev rational functions Rn(z) = Tn((z-1)/(z+1)). Input\n",
    "    is the coefficients of f.'''\n",
    "    \n",
    "    N = len(c) - 1\n",
    "    \n",
    "    b = c*r_[0:N+1]\n",
    "\n",
    "\n",
    "    cp = zeros_like(b)\n",
    "\n",
    "    cp[0] = sum(b[1::2])\n",
    "\n",
    "    evens = b[2::2]\n",
    "    odds = b[1::2]\n",
    "\n",
    "    cp[1:N+1-(N%2):2] = 2*cumsum(evens[-1::-1])[-1::-1]\n",
    "    cp[2:N+1-((N+1)%2):2] = 2*cumsum(odds[-1::-1])[-2::-1]\n",
    "\n",
    "    d0 = -cp[2]/4. + cp[0]/2.\n",
    "    d1 = cp[1]/4. - cp[3]/4.\n",
    "    d2 = -cp[0]/2. + cp[2]/2. - cp[4]/4.\n",
    "    \n",
    "    dn = [ -cp[n-2]/4.  + cp[n]/2. - cp[n+2]/4 for n in range(3,N-1)] \n",
    "    \n",
    "    dn1 = -cp[N-3]/4. + cp[N-1]/2\n",
    "    dn2 = -cp[N-2]/4. + cp[N]/2\n",
    "    \n",
    "    dn = r_[d0,d1,d2,dn,dn1,dn2]\n",
    "    \n",
    "    return dn\n",
    "\n",
    "\n",
    "\n",
    "def Intcheb(y,interval):\n",
    "    '''Clenshaw-Curtis to find definite integral of function y(x) given at\n",
    "    Chebyshev grid points in some interval [a,b]'''\n",
    "    fact = 0.5*(interval[1]-interval[0])\n",
    "    b = cheb(y)\n",
    "    N = len(y) - 1\n",
    "    if N%2 == 0:\n",
    "        w = array([ 2./(-(2*k)**2+1) for k in r_[0:N/2+1]])\n",
    "    else:\n",
    "        w = array([ 2./(-(2*k)**2+1) for k in r_[0:(N-1)/2+1]])\n",
    "    return dot(b[::2],w)*fact\n",
    "\n",
    "\n",
    "\n",
    "def chebI(c,interval,x0=None,f0=None):\n",
    "    if x0==None:\n",
    "        x0=interval[0]\n",
    "\n",
    "    N = len(c) - 1\n",
    "    I = diag(1./(2*r_[0.5,r_[2:N+1]]),-1) -diag(1./(2*r_[1,r_[1:N]]),1)\n",
    "    I[0,1] = 0\n",
    "    \n",
    "    factor = (interval[1]-interval[0])/2.\n",
    "    \n",
    "    ci = dot(I,c)*factor\n",
    "    x = 2*(x0-interval[0])/(interval[1]-interval[0]) - 1\n",
    "    \n",
    "    if x==-1 and f0==None:\n",
    "        ci[0] = -sum((-1)**r_[1:N+1]*ci[1:])\n",
    "    else:\n",
    "        ci[0] = f0 - clenshaw(x,ci)\n",
    "    return ci\n",
    "\n",
    "\n",
    "\n",
    "def cheb_convolve(a,b):\n",
    "    '''Finds the product of two functions whose Chebyshev coefficients are \n",
    "    given by a and b. Output is the coefficiets of the product.'''\n",
    "\n",
    "    M = len(b)\n",
    "    N = len(a)\n",
    "    \n",
    "    if N>M:\n",
    "        b = r_[b,zeros(N-M)]\n",
    "        N = N - 1\n",
    "    elif M>N:\n",
    "        a = r_[a,zeros(M-N)]\n",
    "        N = M - 1\n",
    "    else:\n",
    "        N = N - 1\n",
    "    \n",
    "    a[0] = a[0]*2.\n",
    "    b[0] = b[0]*2.\n",
    "        \n",
    "    c0 = a[0]*b[0] + 2*dot(a[1:],b[1:])\n",
    "    \n",
    "    c1 = [ dot(a[0:k+1][::-1],b[0:k+1]) + dot(a[1:N-k+1],b[k+1:N+1]) + dot(a[k+1:N+1],b[1:N-k+1])  for k in range(1,N) ]\n",
    "\n",
    "    c2 = [ dot(a[k-N:N+1][::-1],b[k-N:N+1])  for k in range(N,2*N+1)]\n",
    "    \n",
    "    c = r_[c0/2,c1,c2]/2.\n",
    "    \n",
    "    return c[:N+1]\n",
    "    \n",
    "def cosT(d,inverse=False):\n",
    "    '''Finds the cosine transform of a given sequence'''\n",
    "    b = []\n",
    "    N = len(d)-1\n",
    "    for n in r_[0:N+1]:\n",
    "        b.append(sum(d*cos(n*r_[0:N+1]*pi/N)))\n",
    "    b = array(b)\n",
    "    if inverse:\n",
    "        return b\n",
    "    else:\n",
    "        b[0] = b[0]/(N)\n",
    "        b[1:] = b[1:]*2/(N)\n",
    "        return b"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The alternate method\n",
    "\n",
    "We assume \n",
    "$$ \\sum_n c_n R_n = \\frac{1}{(1+x)^2} \\sum_n b_n R_n $$\n",
    "\n",
    "alt_matrix and nate_matrix are used to construct the $ b_n $ s from $ R_n $ s"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def alt_matrix(N):\n",
    "    if N < 5:\n",
    "        print(\"Make N >= 5\")\n",
    "    else:\n",
    "        t0,t1,t1=3/4,1/2,1/8\n",
    "        A = np.zeros((N,N))\n",
    "        A[0,0] = 1/16\n",
    "        A[1,0:2] = [-1/4,1/16]\n",
    "        A[2,0:3] = [3/8,-1/4,1/16]\n",
    "        A[3,0:4] = [-1/4,3/8,-1/4,1/16]\n",
    "        A[4,0:5] = [1/16,-1/4,3/8,-1/4,1/16]\n",
    "        for i in range(5,N):\n",
    "            A[i,i-4:i+1] = [1/16,-1/4,3/8,-1/4,1/16]\n",
    "        A[-1,-5:] = [1/16,-1/4,3/8,-1/4,1/8]\n",
    "        return A\n",
    "\n",
    "def nate_matrix(N):\n",
    "    M = zeros((N,N))\n",
    "    K =[16,64,160,320,560]\n",
    "    d = [-1,4,-6,4]\n",
    "    for i in range(N):\n",
    "        K = K+[dot(d,K[-4:])]\n",
    "        M = M + K[i]*eye(N,N,-i)\n",
    "    M[-1,:]=0.5*M[-1,:]   \n",
    "    return M"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def integrator(c,b=[]):\n",
    "    a0 = sum(c)\n",
    "    a1 = sum(-(chebD(c,[-1,1]))*2)\n",
    "    r = c[2:]\n",
    "    rc = r[::-1]\n",
    "    #rc[rc<0.001]=0\n",
    "    A = nate_matrix(len(rc))\n",
    "    #ra = solve(A,rc)\n",
    "    ra = np.dot(A,rc)\n",
    "    a = ra[::-1]\n",
    "    a=a/2\n",
    "    Int = chebI(a,[-1,1],-1,0)\n",
    "    b = [a0,a1] +b\n",
    "    return Int,b\n",
    "\n",
    "def alligrator(c,M,b=[]):\n",
    "    a0 = sum(c)\n",
    "    a1 = sum(-(chebD(c,[-1,1]))*2)\n",
    "    #r =zeros(len(c))\n",
    "    r = c[2:]\n",
    "    rc = r[::-1]\n",
    "    ra = np.dot(M,rc)\n",
    "    a = ra[::-1]\n",
    "    a=a/2\n",
    "    Int = chebI(a,[-1,1],-1,0)\n",
    "    b = [a0,a1] +b\n",
    "    return Int,b\n",
    "\n",
    "def basis(N):\n",
    "    A=np.eye(N)\n",
    "    B=np.zeros((N-2,N))\n",
    "    M=nate_matrix(N-2)\n",
    "    i=0\n",
    "    for row in A:\n",
    "        B[:,i],k= alligrator(row,M)\n",
    "        i+=1\n",
    "    return B"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$ f_x = g $\n",
    "\n",
    "$ f(x)-f(0) = \\int_0^x g(s) ds $\n",
    "\n",
    "$ \\bar{f} = \\mathbb{I} g $"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def derivmatrix(N, interval=[-1,1]):\n",
    "    A=np.eye(N)\n",
    "    B=np.zeros((N,N))\n",
    "    i=0\n",
    "    for row in A:\n",
    "        B[i,:]= Dcheb(row,interval)\n",
    "        i+=1\n",
    "    T=B.transpose() \n",
    "    return T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1.,  2.,  3.],\n",
       "       [ 3.,  4.,  5.],\n",
       "       [ 4.,  6.,  8.]])"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = [[1,2,3],[3,4,5]]\n",
    "sum(A)\n",
    "t = ones(2)\n",
    "f =dot(t,A)\n",
    "r_[A,[f]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "N = 50\n",
    "T = basis(N)\n",
    "c1 = np.ones(N)\n",
    "D = derivmatrix(N)\n",
    "c2 = np.dot(c1,D)\n",
    "I = r_[T,[c1],[c2]]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "50"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "size(I,1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "I is the matrix we want. But, it might not be invertible. Taking the determinant,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "7.368455885374967e-40"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "lin.det(I)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That's pretty much zero so, it isn't invertible."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "9.188793387472161e-60"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "It = np.transpose(I)\n",
    "A = np.dot(I,It)\n",
    "lin.det(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "N=50\n",
    "z=np.cos(np.pi*np.linspace(0,N-1,N)/N-1)\n",
    "L = 1\n",
    "x = L*((z[1:]+1)/(1-z[1:]))\n",
    "x = r_[inf,x]\n",
    "y = x**2\n",
    "y[-1]=0\n",
    "y[-2]=0\n",
    "pI = lin.pinv(I) #pseudoinverse\n",
    "ans = np.dot(pI,y)\n",
    "#ans = lin.solve(I,y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([  3.35068530e+00,   3.90283252e+00,   4.57889611e+00,\n",
       "         5.41865015e+00,   6.47923688e+00,   7.84512961e+00,\n",
       "         9.64546143e+00,   1.20856707e+01,   1.55087683e+01,\n",
       "         2.05227394e+01,   2.82902887e+01,   4.12688875e+01,\n",
       "         6.54230219e+01,   1.18533784e+02,   2.75479426e+02,\n",
       "         1.20822943e+03,   1.41881673e+05,   8.60796321e+02,\n",
       "         2.32515116e+02,   1.05828932e+02,   6.00669458e+01,\n",
       "         3.85275243e+01,   2.67041987e+01,   1.95240972e+01,\n",
       "         1.48398774e+01,   1.16159870e+01,   9.30315337e+00,\n",
       "         7.58803659e+00,   6.28129707e+00,   5.26305273e+00,\n",
       "         4.45440466e+00,   3.80170623e+00,   3.26745122e+00,\n",
       "         2.82477939e+00,   2.45404649e+00,   2.14061770e+00,\n",
       "         1.87340999e+00,   1.64390640e+00,   1.44547537e+00,\n",
       "         1.27289183e+00,   1.12199450e+00,   9.89436838e-01,\n",
       "         8.72503518e-01,   7.68973443e-01,   6.77016304e-01,\n",
       "         5.95113600e-01,   5.21997748e-01,   4.56604694e-01,\n",
       "         3.98036727e-01,   3.45533064e-01])"
      ]
     },
     "execution_count": 43,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f1922f99f98>]"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEHCAYAAACp9y31AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAEeFJREFUeJzt3WuMXfdd7vHvg5ugQ6haWk9L5QsOyOIcc2hCNbhVCU2C\n1GBXLVahSLaiFqqUUVAicRMi8KKR4AUSlRAqTWsMWKYSSQRtAwa5SVMJTnooOfI4yknjlLQjN+Cx\nKnkal5ReUHD48WKWYTOZ8V4zs+bi/L8faWv2+t/Wbyezn1les/asVBWSpHZ820YXIElaXwa/JDXG\n4Jekxhj8ktQYg1+SGmPwS1JjNm3wJzma5HySJ3uMfUuSx5JcTPKuBX0/k+SL3eNn1q5iSboybNrg\nB44B+3qO/SfgZ4F7RxuTvAq4G3gjsBe4O8l3DVeiJF15Nm3wV9UjwIXRtiTfl+TBJKeSfCbJ/+zG\nPlNVTwD/vmCZHwcerqoLVfVV4GH6/zCRpJekl210Act0BLi9qr6Y5I3Ah4Efu8z4bcDZke3Zrk2S\nmnXFBH+S7wTeDPx5kkvN375xFUnSlemKCX7mT0v9c1Vdv4w554CbRra3A387YE2SdMXZtOf4F6qq\nrwFfSvLTAJl33ZhpDwG3JPmu7pe6t3RtktSsTRv8Se4D/h74/iSzSW4DbgVuS/L/gdPAgW7sDyeZ\nBX4a+IMkpwGq6gLwW8DJ7vGbXZskNSv+WWZJasumPeKXJK0Ng1+SGrMpr+rZunVr7dq1a6PLkKQr\nxqlTp75SVRN9xm7K4N+1axfT09MbXYYkXTGS/GPfsZ7qkaTGGPyS1BiDX5IaY/BLUmMMfklqzNjg\nT7Ijyd8keSrJ6SS/sMiYJPlgkpkkTyR5w0jfviRPd313Df0CJEnL0+eI/yLwK1W1B3gTcEeSPQvG\n7Ad2d48p4CMASbYA93T9e4BDi8yVpOZ97V//jbMXvrku+xob/FX15ap6rHv+L8DnefHNTA4AH615\njwKvTPI65m93OFNVZ6rqeeD+bqwkacT+3/sMP/o7f7Mu+1rWOf4ku4AfAv7fgq6l7nTlHbAkqYdz\n//ytddtX7+Dv7oD1ceAXu7+NP6gkU0mmk0zPzc0NvbwkqdMr+JNcxXzo/2lVfWKRIeeAHSPb27u2\npdpfpKqOVNVkVU1OTPT6cxOSpBXoc1VPgD8GPl9Vv7vEsOPAe7qre94EPFdVX2b+5ie7k1yb5Grg\nYDdWkrRB+vyRth8B3g18LsnjXdtvADsBquowcAJ4GzADfBN4b9d3McmdzN/ucAtwtKpOD/oKJEnL\nMjb4q+r/AhkzpoA7lug7wfwPBknSJuAndyWpMQa/JDXG4Jekxhj8ktQYg1+SGmPwS1JjDH5JaozB\nL0mNMfglqTEGvyQ1xuCXpMYY/JLUGINfkhpj8EtSYwx+SWqMwS9JjRl7I5YkR4G3A+er6n8v0v+r\nwK0j6/0vYKKqLiR5BvgX4AXgYlVNDlW4JGll+hzxHwP2LdVZVR+oquur6nrg14H/U1UXRobc3PUb\n+pK0CYwN/qp6BLgwblznEHDfqiqSJK2pwc7xJ/kO5v9l8PGR5gI+neRUkqkx86eSTCeZnpubG6os\nSdICQ/5y9x3A3y04zXNDdwpoP3BHkrcsNbmqjlTVZFVNTkxMDFiWJGnUkMF/kAWnearqXPf1PPAA\nsHfA/UmSVmCQ4E/yCuBG4C9H2q5J8vJLz4FbgCeH2J8kaeX6XM55H3ATsDXJLHA3cBVAVR3uhr0T\n+FRVfWNk6muBB5Jc2s+9VfXgcKVLklZibPBX1aEeY44xf9nnaNsZ4LqVFiZJWht+cleSGmPwS1Jj\nDH5JaozBL0mNMfglqTEGvyQ1xuCXpMYY/JLUGINfkhpj8EtSYwx+SWqMwS9JjTH4JakxBr8kNcbg\nl6TGjA3+JEeTnE+y6N2zktyU5Lkkj3eP94/07UvydJKZJHcNWbgkaWX6HPEfA/aNGfOZqrq+e/wm\nQJItwD3M32h9D3AoyZ7VFCtJWr2xwV9VjwAXVrD2XmCmqs5U1fPA/cCBFawjSRrQUOf435zkiSSf\nTPIDXds24OzImNmuTZK0gcbec7eHx4CdVfX1JG8D/gLYvdxFkkwBUwA7d+4coCxJ0mJWfcRfVV+r\nqq93z08AVyXZCpwDdowM3d61LbXOkaqarKrJiYmJ1ZYlSVrCqoM/yXcnSfd8b7fms8BJYHeSa5Nc\nDRwEjq92f5Kk1Rl7qifJfcBNwNYks8DdwFUAVXUYeBfw80kuAt8CDlZVAReT3Ak8BGwBjlbV6TV5\nFZKk3sYGf1UdGtP/IeBDS/SdAE6srDRJ0lrwk7uS1BiDX5IaY/BLUmMMfklqjMEvSY0x+CWpMQa/\nJDXG4Jekxhj8ktQYg1+SGmPwS1JjDH5JaozBL0mNMfglqTEGvyQ1xuCXpMaMDf4kR5OcT/LkEv23\nJnkiyeeSfDbJdSN9z3TtjyeZHrJwSdLK9DniPwbsu0z/l4Abq+oHgd8Cjizov7mqrq+qyZWVKEka\nUp9bLz6SZNdl+j87svkosH31ZUmS1srQ5/hvAz45sl3Ap5OcSjJ1uYlJppJMJ5mem5sbuCxJ0iVj\nj/j7SnIz88F/w0jzDVV1LslrgIeT/ENVPbLY/Ko6QneaaHJysoaqS5L03w1yxJ/k9cAfAQeq6tlL\n7VV1rvt6HngA2DvE/iRJK7fq4E+yE/gE8O6q+sJI+zVJXn7pOXALsOiVQZKk9TP2VE+S+4CbgK1J\nZoG7gasAquow8H7g1cCHkwBc7K7geS3wQNf2MuDeqnpwDV6DJGkZ+lzVc2hM//uA9y3Sfga47sUz\nJEkbyU/uSlJjDH5JaozBL0mNMfglqTEGvyQ1xuCXpMYY/JLUGINfkhpj8EtSYwx+SWqMwS9JjTH4\nJakxBr8kNcbgl6TGGPyS1BiDX5IaMzb4kxxNcj7JordNzLwPJplJ8kSSN4z07UvydNd315CFS5JW\nps8R/zFg32X69wO7u8cU8BGAJFuAe7r+PcChJHtWU6wkafXGBn9VPQJcuMyQA8BHa96jwCuTvA7Y\nC8xU1Zmqeh64vxsrSdpAQ5zj3wacHdme7dqWal9Ukqkk00mm5+bmBihLkrSYTfPL3ao6UlWTVTU5\nMTGx0eVI0kvWywZY4xywY2R7e9d21RLtkqQNNMQR/3HgPd3VPW8CnquqLwMngd1Jrk1yNXCwGytJ\n2kBjj/iT3AfcBGxNMgvczfzRPFV1GDgBvA2YAb4JvLfru5jkTuAhYAtwtKpOr8FrkCQtw9jgr6pD\nY/oLuGOJvhPM/2CQJG0Sm+aXu5Kk9WHwS1JjDH5JaozBL0mNMfglqTEGvyQ1xuCXpMYY/JLUGINf\nkhpj8EtSYwx+SWqMwS9JjTH4JakxBr8kNcbgl6TG9Ar+JPuSPJ1kJsldi/T/apLHu8eTSV5I8qqu\n75kkn+v6pod+AZKk5elzB64twD3AW4FZ4GSS41X11KUxVfUB4APd+HcAv1RVF0aWubmqvjJo5ZKk\nFelzxL8XmKmqM1X1PHA/cOAy4w8B9w1RnCRpeH2CfxtwdmR7tmt7kSTfAewDPj7SXMCnk5xKMrXU\nTpJMJZlOMj03N9ejLEnSSgz9y913AH+34DTPDVV1PbAfuCPJWxabWFVHqmqyqiYnJiYGLkuSdEmf\n4D8H7BjZ3t61LeYgC07zVNW57ut54AHmTx1JkjZIn+A/CexOcm2Sq5kP9+MLByV5BXAj8Jcjbdck\nefml58AtwJNDFC5JWpmxV/VU1cUkdwIPAVuAo1V1OsntXf/hbug7gU9V1TdGpr8WeCDJpX3dW1UP\nDvkCJEnLMzb4AarqBHBiQdvhBdvHgGML2s4A162qQknSoPzkriQ1xuCXpMYY/JLUGINfkhpj8EtS\nYwx+SWqMwS9JjTH4JakxBr8kNcbgl6TGGPyS1BiDX5IaY/BLUmMMfklqjMEvSY0x+CWpMb2CP8m+\nJE8nmUly1yL9NyV5Lsnj3eP9fedKktbX2DtwJdkC3AO8FZgFTiY5XlVPLRj6map6+wrnSpLWSZ8j\n/r3ATFWdqarngfuBAz3XX81cSdIa6BP824CzI9uzXdtCb07yRJJPJvmBZc4lyVSS6STTc3NzPcqS\nJK3EUL/cfQzYWVWvB34f+IvlLlBVR6pqsqomJyYmBipLkrRQn+A/B+wY2d7etf2nqvpaVX29e34C\nuCrJ1j5zJUnrq0/wnwR2J7k2ydXAQeD46IAk350k3fO93brP9pkrSVpfY6/qqaqLSe4EHgK2AEer\n6nSS27v+w8C7gJ9PchH4FnCwqgpYdO4avRZJUg9jgx/+8/TNiQVth0eefwj4UN+5kqSN4yd3Jakx\nBr8kNcbgl6TGGPyS1BiDX5IaY/BLUmMMfklqjMEvSY0x+CWpMQa/JDXG4Jekxhj8ktQYg1+SGmPw\nS1JjDH5Jakyv4E+yL8nTSWaS3LVI/63djdY/l+SzSa4b6Xuma388yfSQxUuSlm/sjViSbAHuAd4K\nzAInkxyvqqdGhn0JuLGqvppkP3AEeONI/81V9ZUB65YkrVCfI/69wExVnamq54H7gQOjA6rqs1X1\n1W7zUeZvqi5J2oT6BP824OzI9mzXtpTbgE+ObBfw6SSnkkwtv0RJ0pB63XO3ryQ3Mx/8N4w031BV\n55K8Bng4yT9U1SOLzJ0CpgB27tw5ZFmSpBF9jvjPATtGtrd3bf9NktcDfwQcqKpnL7VX1bnu63ng\nAeZPHb1IVR2pqsmqmpyYmOj/CiRJy9In+E8Cu5Ncm+Rq4CBwfHRAkp3AJ4B3V9UXRtqvSfLyS8+B\nW4AnhypekrR8Y0/1VNXFJHcCDwFbgKNVdTrJ7V3/YeD9wKuBDycBuFhVk8BrgQe6tpcB91bVg2vy\nSiRJvfQ6x19VJ4ATC9oOjzx/H/C+ReadAa5b2C5J2jh+cleSGmPwS1JjDH5JaozBL0mNMfglqTEG\nvyQ1xuCXpMYY/JLUGINfkhpj8EtSYwx+SWqMwS9JjTH4JakxBr8kNcbgl6TGGPyS1JhewZ9kX5Kn\nk8wkuWuR/iT5YNf/RJI39J0rSVpfY4M/yRbgHmA/sAc4lGTPgmH7gd3dYwr4yDLmSpLWUZ9bL+4F\nZrrbKJLkfuAA8NTImAPAR6uqgEeTvDLJ64BdPeYO5tc+9gT/9sK/r8XSkjbY335hjgvfeH6jy3hJ\n6BP824CzI9uzwBt7jNnWcy4ASaaY/9cCO3fu7FHWiz32T1/lXy++sKK5kja3FkJ/x6v+x7rsp9fN\n1tdDVR0BjgBMTk7WStZ4+JdvHLQmSXop6hP854AdI9vbu7Y+Y67qMVeStI76XNVzEtid5NokVwMH\ngeMLxhwH3tNd3fMm4Lmq+nLPuZKkdTT2iL+qLia5E3gI2AIcrarTSW7v+g8DJ4C3ATPAN4H3Xm7u\nmrwSSVIvmb8QZ3OZnJys6enpjS5Dkq4YSU5V1WSfsX5yV5IaY/BLUmMMfklqjMEvSY3ZlL/cTfIc\n8MWBlnsF8NwGrbGSecuZsxX4yjLXb9EQ3wPrZSNrXct9D732atd7Kb6nv6eqJnqNrKpN9wCObKa1\nVrrGSuYtZw4wvdH/r66Ex5DfTy/lWtdy30Ovvdr1Wn9Pb9ZTPX+1ydZa6RormTfka9e8K+m/6UbW\nupb7Hnrt1a7X9Ht6U57qUT9JpqvndbuSNr/1ek9v1iN+9XNkowuQNKh1eU97xC9JjfGIX5IaY/BL\nUmMMfklqjMH/EpXke5P8cZKPbXQtklYmyTVJ/iTJHya5dah1Df5NKMnRJOeTPLmgfV+Sp5PMJLnr\ncmtU1Zmqum1tK5W0XMt8f/8k8LGq+jngJ4aqweDfnI4B+0YbkmwB7gH2A3uAQ0n2JPnBJH+94PGa\n9S9ZUk/H6Pn+Zv52tWe7YS8MVcCmudm6/ktVPZJk14LmvcBMVZ0BSHI/cKCqfht4+/pWKGmllvP+\nBmaZD//HGfBA3SP+K8c2/usnP8x/Q2xbanCSVyc5DPxQkl9f6+IkrcpS7+9PAD+V5CMM+KcfPOJ/\niaqqZ4HbN7oOSStXVd+gu4f5kDziv3KcA3aMbG/v2iRd+db1/W3wXzlOAruTXJvkauAgcHyDa5I0\njHV9fxv8m1CS+4C/B74/yWyS26rqInAn8BDweeDPqur0RtYpafk2w/vbP9ImSY3xiF+SGmPwS1Jj\nDH5JaozBL0mNMfglqTEGvyQ1xuCXpMYY/JLUGINfkhrzH4iUipuxnUrYAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f1922f39710>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "ansr = icheb(ans)\n",
    "semilogx(z,y)"
   ]
  },
  {
   "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.6.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}