{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Chapter10 Elliptical Partial Differential Equations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example 10.1 page 441"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "u:\n",
      "[[   0.      0.      0.      0.      0.  ]\n",
      " [   0.      6.25   12.5    18.75   25.  ]\n",
      " [   0.     12.5    25.     37.5    50.  ]\n",
      " [   0.     18.75   37.5    56.25   75.  ]\n",
      " [   0.     25.     50.     75.    100.  ]]\n"
     ]
    }
   ],
   "source": [
    "from __future__ import division\n",
    "from numpy import arange,zeros\n",
    "from math import ceil\n",
    "\n",
    "h = 1/4\n",
    "xf = 1\n",
    "yf = 1\n",
    "x = arange(0,h+xf,h)\n",
    "y = arange(0,h+yf,h)\n",
    "m = len(y)-1\n",
    "n = len(x)-1\n",
    "\n",
    "u = zeros([m+1,n+1])\n",
    "u[m,:] = [100*xx for xx in x]\n",
    "u[:,n] = [100*yy for yy in y] # dash transpose\n",
    "u0 = u\n",
    "\n",
    "I = int(ceil(m/2)) #\n",
    "J = int(ceil(n/2)) #\n",
    "\n",
    "u[J,I] = (u0[J-2,I-2] + u0[J-2,I+2] + u0[J+2,I-2] + u0[J+2,I+2]) / 4\n",
    "\n",
    "for j in [J-1,J+1]:\n",
    "    for i in [I-1,I+1]:\n",
    "        u[j,i] = (u[j-1,i-1] + u[j-1,i+1] + u[j+1,i-1] + u[j+1,i+1]) / 4\n",
    "    \n",
    "\n",
    "# j and i may have loose 1\n",
    "j1 = [J-1,J,J,J+1]\n",
    "i1 = [I,I-1,I+1,I]\n",
    "for k in range(0,4):\n",
    "    i = i1[(k)]\n",
    "    j = j1[(k)]\n",
    "    u[j,i] = (u[j,i-1] + u[j,i+1] + u[j-1,i] + u[j+1,i]) / 4\n",
    "\n",
    "\n",
    "print \"u:\\n\",u"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example 10.3 page 442"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "   r       u11       u21       u31       u12        u22        u32        u13        u23        u33\n",
      "   ------------------------------------------------------------------------------------------------\n",
      "   1    6.2500    9.3750    6.2500    18.7500    25.0000    18.7500    43.7500    53.1250    43.7500  \n",
      "   1    7.0312    9.5703    7.0801    18.9453    25.0977    18.9819    43.0176    52.9663    42.9871  \n",
      "   1    7.1289    9.8267    7.2021    18.8110    25.1465    18.8339    42.9443    52.7695    42.9008  \n"
     ]
    }
   ],
   "source": [
    "from numpy import zeros\n",
    "m = 5\n",
    "n = 5\n",
    "u = zeros([m,n])\n",
    "u[m-1,:] = [50,100, 100, 100, 50]\n",
    "u0 = u\n",
    "I = int(ceil(m/2))-1\n",
    "J = int(ceil(n/2))-1\n",
    "\n",
    "u[J,I] = (u0[J-2,I-2] + u0[J-2,I+2] + u0[J+2,I-2] + u0[J+2,I+2]) / 4\n",
    "\n",
    "for j in [J-1,J+1]:\n",
    "    for i in [I-1,I+1]:\n",
    "        u[j,i] = (u[j-1,i-1] + u[j-1,i+1] + u[j+1,i-1] + u[j+1,i+1]) / 4\n",
    "    \n",
    "\n",
    "j1 = [J-1, J ,J ,J+1]\n",
    "i1 = [I ,I-1, I+1, I]\n",
    "for k in range(0,4):\n",
    "    i = i1[(k)]\n",
    "    j = j1[(k)]\n",
    "    u[j,i] = (u[j,i-1] + u[j,i+1] + u[j-1,i] + u[j+1,i]) / 4\n",
    "\n",
    "\n",
    "kf = 2\n",
    "tab = zeros([kf+1,(m-2)*(n-2)])\n",
    "row = []\n",
    "for j in range(1,n-1):\n",
    "    for xx in u[j,range(1,m-1)]:\n",
    "        row.append(xx)\n",
    "\n",
    "tab[0,:] = row;\n",
    "for k in range(0,kf):\n",
    "    row = [];\n",
    "    for j in range(1,n-1):\n",
    "        for i in range(1,m-1):\n",
    "            u[j,i] = (u[j,i-1] + u[j,i+1] + u[j-1,i] + u[j+1,i]) / 4;\n",
    "        for xx in u[j,range(1,m-1)]: \n",
    "            row.append(xx)\n",
    "    \n",
    "    row = [round(roww*10**4,4)/10**4 for roww in row]\n",
    "    tab[k+1,:] = row\n",
    "\n",
    "print \"%4s %9s %9s %9s %9s %10s %10s %10s %10s %10s\"%('r','u11','u21','u31','u12','u22','u32','u13','u23','u33')\n",
    "print '  ','-'*96\n",
    "for tabb in tab:\n",
    "    i=1;\n",
    "    print '  ',i,' ',\n",
    "    for tt in tabb:\n",
    "        print ' %0.4f'%(tt),' ',\n",
    "    i+=1;\n",
    "    print"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example 10.4 page 444"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "   r       u11       u21       u12       u22\n",
      "   1    1.3889    1.7670    2.0170    2.6466  \n",
      "   2    2.3349    2.6651    2.9151    3.0957  \n",
      "   3    2.7840    2.8897    3.1397    3.2079  \n",
      "   4    2.8962    2.9458    3.1958    3.2360  \n",
      "   5    2.9243    2.9598    3.2098    3.2430  \n"
     ]
    }
   ],
   "source": [
    "from numpy import arange\n",
    "h = 1/3\n",
    "x = arange(0,h+1,h)\n",
    "y = arange(0,h+1,h)\n",
    "m = len(y)\n",
    "n = len(x)\n",
    "u = zeros([m,n])\n",
    "u[m-1,range(1,n-1)] = 1\n",
    "\n",
    "kf = 5\n",
    "tab = zeros([kf,(m-2)*(n-2)])\n",
    "for k in range(0,kf):\n",
    "    row = []\n",
    "    for j in range(1,n-1):\n",
    "        for i in range(1,m-1):\n",
    "            constant = 10/9* (5 + 1/9*(i-1)**2 + 1/9*(j-1)**2)\n",
    "            u[j,i] = (u[j,i-1] + u[j,i+1] + u[j-1,i] + u[j+1,i] + constant) / 4\n",
    "        for xx in u[j,range(1,m-1)]:\n",
    "            row.append(xx)\n",
    "    \n",
    "    row = [round(roww*10**4)/10**4 for roww in row]\n",
    "    tab[k,:] = row\n",
    "\n",
    "print \"%4s %9s %9s %9s %9s\"%('r','u11','u21','u12','u22')\n",
    "i=1;\n",
    "for tabb in tab:\n",
    "    \n",
    "    print '  ',i,' ',\n",
    "    for tt in tabb:\n",
    "        print ' %0.4f'%(tt),' ',\n",
    "    i+=1;\n",
    "    print"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example 10.5 page 446"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "u:\n",
      "[[  0.    0.    0.    0.    0. ]\n",
      " [  0.    4.    8.   12.   16. ]\n",
      " [  0.    6.5  16.   22.5  32. ]\n",
      " [  0.    6.   15.   30.   48. ]\n",
      " [  0.    1.    8.   27.   64. ]]\n",
      "\n",
      "\n",
      "\t     u11       u21        u31        u12       u22        u32        u13       u23       u33\n",
      "     1    4.0000    8.0000    12.0000    6.5000    16.0000    22.5000    6.0000    15.0000    30.0000  \n",
      "     2    3.5605    7.8712    11.5228    6.3712    12.4085    21.3079    5.5228    13.8079    27.1039  \n",
      "     3    3.5606    6.7013    10.9128    5.2013    11.6421    20.2611    4.9128    12.7611    26.9886  \n",
      "     4    2.8750    6.2983    10.5929    4.7983    10.9244    20.1033    4.5929    12.6033    26.9160  \n",
      "     5    2.7568    6.0290    10.5228    4.5290    10.7975    20.0515    4.5228    12.5515    26.8981  \n",
      "     6    2.6193    5.9773    10.5045    4.4773    10.7587    20.0384    4.5045    12.5384    26.8935  \n",
      "     7    2.6127    5.9675    10.5010    4.4675    10.7520    20.0363    4.5010    12.5363    26.8931  \n",
      "     8    2.6081    5.9649    10.5002    4.4649    10.7504    20.0358    4.5002    12.5358    26.8929  \n",
      "     9    2.6073    5.9644    10.5000    4.4644    10.7501    20.0357    4.5000    12.5357    26.8929  \n",
      "     10    2.6072    5.9643    10.5000    4.4643    10.7500    20.0357    4.5000    12.5357    26.8929  \n",
      "     11    2.6072    5.9643    10.5000    4.4643    10.7500    20.0357    4.5000    12.5357    26.8929  \n"
     ]
    }
   ],
   "source": [
    "from math import ceil\n",
    "from numpy import cos,pi,sqrt,zeros\n",
    "x = range(0,5)\n",
    "y = range(0,5)\n",
    "m = len(y)\n",
    "n = len(x)\n",
    "u = zeros([m,n])\n",
    "u[m-1,:] = [xx**3 for xx in x]\n",
    "u[:,n-1] = [16*yy for yy in y]\n",
    "u0 = u\n",
    "\n",
    "I = int(ceil(m/2))-1\n",
    "J = int(ceil(n/2))-1\n",
    "\n",
    "u[J,I] = (u0[J-2,I-2] + u0[J-2,I+2] + u0[J+2,I-2] + u0[J+2,I+2]) / 4\n",
    "\n",
    "for j in [J-1,J+1]:\n",
    "    for i in [I-1,I+1]:\n",
    "        u[j,i] = (u[j-1,i-1] + u[j-1,i+1] + u[j+1,i-1] + u[j+1,i+1]) / 4\n",
    "    \n",
    "\n",
    "\n",
    "j1 = [J-1,J ,J ,J+1]\n",
    "i1 = [I, I-1, I+1, I]\n",
    "for k in range(0,4):\n",
    "    i = i1[(k)]\n",
    "    j = j1[(k)]\n",
    "    u[j,i] = (u[j,i-1] + u[j,i+1] + u[j-1,i] + u[j+1,i]) / 4\n",
    "\n",
    "print \"u:\\n\",u\n",
    "\n",
    "p = m-1\n",
    "q = n-1\n",
    "c = cos(pi/p) + cos(pi/q)\n",
    "w = 4/(2+sqrt(4-c**2))\n",
    "w = round(w*10**3)/10**3\n",
    "\n",
    "kf = 10\n",
    "tab = zeros([kf+1,(m-2)*(n-2)])\n",
    "row = []\n",
    "for j in range(1,n-1):\n",
    "    for xx in u[j,range(1,m-1)]:\n",
    "        row.append(xx)\n",
    "\n",
    "tab[0,:] = row\n",
    "for k in range(0,kf):\n",
    "    row = []\n",
    "    for j in range(1,n-1):\n",
    "        for i in range(1,m-1):\n",
    "            u[j,i] = (u[j,i-1] + u[j,i+1] + u[j-1,i] + u[j+1,i]) *w/4 + (1-w)*u[j,i]\n",
    "        for xx in u[j,range(1,m-1)]:\n",
    "            row.append(xx)\n",
    "    \n",
    "    row = [round(roww*10**4)/10**4 for roww in row]\n",
    "    tab[k+1,:] = row\n",
    "\n",
    "print \"\\n\\n\\t%8s %9s %10s %10s %9s %10s %10s %9s %9s\"%('u11','u21','u31','u12','u22','u32','u13','u23','u33')\n",
    "i=1;\n",
    "for tabb in tab:\n",
    "    \n",
    "    print '    ',i,' ',\n",
    "    for tt in tabb:\n",
    "        print ' %0.4f'%(tt),' ',\n",
    "    i+=1;\n",
    "    print\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}