Pade_approx.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Using matplotlib backend: Qt5Agg\n",
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    }
   ],
   "source": [
    "%pylab"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "21.0\n"
     ]
    }
   ],
   "source": [
    "def C_z(z_values, f_values, z): #z_values are the points at which we know the value of our function, f_values the values of our function at those points, z the point where we want to evaluate the continued fraction C_z;\n",
    "\n",
    "    N = z_values.size\n",
    "    \n",
    "    g_mat = np.zeros([N,N]) #this is the matrix whose elements are: g_mat[i,j] = g_i(z_j), cfr eqn. (A2). Note that only diagonal and upper diagonal elements matter. The coefficients a_i will be found on the diagonal: a_i = g_mat[i,i];\n",
    "    \"\"\"\n",
    "    for i in range(N):     #filling the first row with the initial data;\n",
    "        g_mat[0,i] = f_values[i]\"\"\"\n",
    "    #can be done diretly:\n",
    "    g_mat[0]=f_values\n",
    "    \n",
    "    for i in range(N-1): #filling the rest of the matrix using the recursion relation;\n",
    "        for j in range(N-(i+1)):\n",
    "            g_mat[i+1,i+1+j] = (g_mat[i,i]-g_mat[i,j+i+1])/((z_values[j+1+i]-z_values[i])*g_mat[i,j+i+1])\n",
    "\n",
    "    \"\"\"\n",
    "    a_coeff = np.zeros(N)\n",
    "    \n",
    "    for i in range (N):\n",
    "        a_coeff[i] = g_mat[i,i]\n",
    "    \"\"\"\n",
    "    #also more effenciently:\n",
    "    a_coeff = np.diag(g_mat)\n",
    "\n",
    "    #now we evaluate the contiued fraction at point z with previously computed coeeficients a_i;\n",
    "\n",
    "    A = np.zeros(N+1)\n",
    "    B = np.zeros(N+1)\n",
    "    A[1] = a_coeff[0]\n",
    "    B[0] = 1\n",
    "    B[1] = 1\n",
    "    \n",
    "    for i in range(N-1):\n",
    "        A[i+2] = A[i+1] + (z-z_values[i])*a_coeff[i+1]*A[i]\n",
    "        B[i+2] = B[i+1] + (z-z_values[i])*a_coeff[i+1]*B[i]\n",
    "\n",
    "    return A[N]/B[N]\n",
    "\n",
    "#try: if C_z is computed for one of the given points, it returns the expected result: the code could be actually working properly!\n",
    "\n",
    "\n",
    "points = np.array([11,2,3,12,22,123,9,15,55])\n",
    "values = np.array([11,21,3,7,1,12,43,54,88])\n",
    "\n",
    "print(C_z(points,values,2))\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "f = lambda z: C_z(points,values,z)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "fvectorized = np.vectorize(f)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 11.,  21.,   3.,   7.,   1.,  12.,  43.,  54.,  88.])"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fvectorized(points)"
   ]
  }
 ],
 "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.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}