{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Chapter 28: Numerical Solution Of Partial Differential Equations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example 28.1, page no. 725"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "D=Bˆ2−4AC\n",
      "if D<0 then elliptic if D=0 then parabolic if D>0 then hyperboic\n",
      "(i) A=xˆ2, B1−yˆ2 D=4ˆ2−4∗1∗4=0 so the equation is PARABOLIC \n",
      "(ii) D=4xˆ2(yˆ2−1)\n",
      "for −inf<x<inf and −1<y<1 D<0 \n",
      "So the equation is ELLIPTIC \n",
      "(iii) A=1+xˆ2, B=5+2xˆ2, C=4+xˆ2\n",
      "D=9>0\n",
      "So the equation is HYPERBOLIC \n"
     ]
    }
   ],
   "source": [
    "print \"D=Bˆ2−4AC\"\n",
    "print \"if D<0 then elliptic if D=0 then parabolic if D>0 then hyperboic\"\n",
    "print \"(i) A=xˆ2, B1−yˆ2 D=4ˆ2−4∗1∗4=0 so the equation is PARABOLIC \"\n",
    "print \"(ii) D=4xˆ2(yˆ2−1)\"\n",
    "print \"for −inf<x<inf and −1<y<1 D<0 \"\n",
    "print \"So the equation is ELLIPTIC \"\n",
    "print \"(iii) A=1+xˆ2, B=5+2xˆ2, C=4+xˆ2\"\n",
    "print \"D=9>0\"\n",
    "print \"So the equation is HYPERBOLIC \""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example 28.2, page no. 726"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "See figure in question\n",
      "From symmetry u7=u1, u8=u2, u9=u3, u3=u1, u6=u4, u9=u7\n",
      "u5=1/4∗(2000+2000+1000+1000)=1500\n",
      "u1=1/4(0=1500+1000+2000)=1125\n",
      "u2=1/4∗(1125+1125+1000+1500)=1188\n",
      "u4=1/4(2000+1500+1125+1125)=1438\n",
      "1125 1188 1438 1500\n",
      "Iterations :\n",
      "\n",
      "1301 1414 1164 1031\n",
      "\n",
      "1250 1337 1087 1019\n",
      "\n",
      "1199 1312 1062 981\n",
      "\n",
      "1174 1287 1037 968\n",
      "\n",
      "1155 1274 1024 956\n",
      "\n",
      "1144 1265 1015 949\n"
     ]
    }
   ],
   "source": [
    "print \"See figure in question\"\n",
    "print \"From symmetry u7=u1, u8=u2, u9=u3, u3=u1, u6=u4, u9=u7\"\n",
    "print \"u5=1/4∗(2000+2000+1000+1000)=1500\"\n",
    "u5 = 1500\n",
    "print \"u1=1/4(0=1500+1000+2000)=1125\"\n",
    "u1 = 1125\n",
    "print \"u2=1/4∗(1125+1125+1000+1500)=1188\"\n",
    "u2 = 1188\n",
    "print \"u4=1/4(2000+1500+1125+1125)=1438\"\n",
    "u4 = 1438\n",
    "print u1,u2,u4,u5 \n",
    "print \"Iterations :\"\n",
    "for i in range(1,7):\n",
    "    u11 = (1000+u2+500+u4)/4\n",
    "    u22 = (u11+u1+1000+u5)/4 \n",
    "    u44 = (2000+u5+u11+u1)/4\n",
    "    u55 = (u44+u4+u22+u2)/4\n",
    "    print \"\"\n",
    "    print u55,u44,u22,u11\n",
    "    u1 = u11 \n",
    "    u2 = u22 \n",
    "    u4 = u44 \n",
    "    u5 = u55"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example 28.3, page no. 727"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "See figure in question\n",
      "To find the initial values of u1 u2 u3 u4 we assume u4=0 \n",
      "u1=1/4∗(1000+0+1000+2000)=1000\n",
      "u2=1/4(1000+500+1000+500)=625 \n",
      "u3=1/4∗(2000+0+1000+500)=875\n",
      "u4=1/4(875+0+625+0)=375\n",
      "1000 625 875 375\n",
      "Iterations:\n",
      "\n",
      "437 1000 750 1125\n",
      "\n",
      "453 1031 781 1187\n",
      "\n",
      "457 1039 789 1203\n",
      "\n",
      "458 1041 791 1207\n",
      "\n",
      "458 1041 791 1208\n",
      "\n",
      "458 1041 791 1208\n"
     ]
    }
   ],
   "source": [
    "print \"See figure in question\"\n",
    "print \"To find the initial values of u1 u2 u3 u4 we assume u4=0 \"\n",
    "print \"u1=1/4∗(1000+0+1000+2000)=1000\"\n",
    "u1 = 1000\n",
    "print \"u2=1/4(1000+500+1000+500)=625 \"\n",
    "u2 = 625\n",
    "print \"u3=1/4∗(2000+0+1000+500)=875\"\n",
    "u3 = 875\n",
    "print \"u4=1/4(875+0+625+0)=375\"\n",
    "u4 = 375\n",
    "print u1,u2,u3,u4 \n",
    "print \"Iterations:\"\n",
    "for i in range(1,7):\n",
    "    u11 = (2000+u2+1000+u3)/4 \n",
    "    u22 = (u11+500+1000+u4)/4\n",
    "    u33 = (2000+u4+u11+500)/4\n",
    "    u44 = (u33+0+u22+0)/4\n",
    "    print \"\"\n",
    "    print u44,u33,u22,u11\n",
    "    u1 = u11 \n",
    "    u2 = u22 \n",
    "    u4 = u44 \n",
    "    u3 = u33"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example 28.5, page no. 729"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Here cˆ2=4, h=1, k=1/8, therefore alpha=(cˆ2)∗k/(hˆ2)\n",
      "Using bendre−schmidits recurrence relation ie u(i)(j +1)=t∗u(i−1)(j)+t∗u(i+1)(j)+(1−2t)∗u(i,j)\n",
      "Now since u(0,t)=0=u(8,t) therefore u(0,i)=0 and u(8,j)=0 and u(x,0)=4x−1/2xˆ2 \n",
      "0.0\n",
      "-4.0\n",
      "0.0\n",
      "4.0\n",
      "8.0\n",
      "12.0\n",
      "16.0\n"
     ]
    }
   ],
   "source": [
    "import numpy\n",
    "print \"Here cˆ2=4, h=1, k=1/8, therefore alpha=(cˆ2)∗k/(hˆ2)\"\n",
    "print \"Using bendre−schmidits recurrence relation ie u(i)(j +1)=t∗u(i−1)(j)+t∗u(i+1)(j)+(1−2t)∗u(i,j)\"\n",
    "print \"Now since u(0,t)=0=u(8,t) therefore u(0,i)=0 and u(8,j)=0 and u(x,0)=4x−1/2xˆ2 \"\n",
    "c = 2\n",
    "h = 1\n",
    "k = 1/8\n",
    "t = (c**2)*k/(h**2)\n",
    "A = numpy.ones((9,9))\n",
    "for i in range(0,9):\n",
    "    for j in range (0,9):\n",
    "        A[0,i] = 0\n",
    "        A[8,i] = 0\n",
    "        A[i,0] = 4*(i-1)-1/2*(i-1)**2\n",
    "for i in range(1,8):\n",
    "    for j in range(1,7):\n",
    "        A[i,j] = t*A[i-1,j-1]+t*A[i+1,j-1]+(1-2*t)*A[i-1,j-1]\n",
    "for i in range(1,8):\n",
    "    j = 2\n",
    "    print A[i,j]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example 28.6, page no. 730"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Here cˆ2=1, h=1/3, k=1/36, therefore t=(cˆ2)∗k/(hˆ2)=1/4 \n",
      "So bendre−schmidits recurrence relation ie u(i)(j+1)=1/4(u(i−1)(j)+u(i+1)(j)+2u(i,j)\n",
      "Now since u(0,t)=0=u(1,t) therefore u(0,i)=0 and u(1,j)=0 and u(x,0)=sin(%pi)x\n",
      "-0.433012701892\n",
      "0.216506350946\n",
      "0.649519052838\n",
      "0.433012701892\n",
      "-0.216506350946\n",
      "-0.649519052838\n",
      "-0.433012701892\n"
     ]
    }
   ],
   "source": [
    "import numpy,math\n",
    "\n",
    "print \"Here cˆ2=1, h=1/3, k=1/36, therefore t=(cˆ2)∗k/(hˆ2)=1/4 \"\n",
    "print \"So bendre−schmidits recurrence relation ie u(i)(j+1)=1/4(u(i−1)(j)+u(i+1)(j)+2u(i,j)\"\n",
    "print \"Now since u(0,t)=0=u(1,t) therefore u(0,i)=0 and u(1,j)=0 and u(x,0)=sin(%pi)x\"\n",
    "c = 1.\n",
    "h = 1./3\n",
    "k = 1./36\n",
    "t = (c**2)*k/(h**2)\n",
    "A = numpy.ones((9,9))\n",
    "for i in range(0,9):\n",
    "    for j in range(0,9):\n",
    "        A[0,i] = 0\n",
    "        A[1,i] = 0\n",
    "        A[i,0] = math.sin(math.pi/3*(i-1)) \n",
    "for i in range(1,8):\n",
    "    for j in range(1,8):\n",
    "        A[i,j] = t*A[i-1,j-1]+t*A[i+1,j-1]+(1-2*t)*A[i-1,j-1]\n",
    "for i in range(1,8):\n",
    "    j = 1\n",
    "    print A[i,j]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example 28.7, page no. 732"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Here cˆ2=16, taking h=1, finding k such that cˆ2tˆ2=1 \n",
      "So bendre−schmidits recurrence relation ie u(i)(j+1)=(16tˆ2(u(i−1)(j)+u(i+1)(j))+2(1−16∗tˆ2u(i,j)−u(i)(j−1)\n",
      "Now since u(0,t)=0=u(5,t) therefore u(0,i)=0 and u(5,j)=0 and u(x,0)=xˆ2(5−x)\n",
      "Also from 1st derivative (u(i)(j+1)−u(i,j−1))/2k=g(x) and g(x)=0 in this case\n",
      "So if j=0 this gives u(i)(1)=1/2∗(u(i−1)(0)+u(i+1)(0))\n",
      "   0.0    0.0    0.0    0.0    0.0 \n",
      "   0.0    4.0    8.0    12.0    16.0 \n",
      "   0.0    12.0    24.0    36.0    48.0 \n",
      "   0.0    18.0    36.0    54.0    72.0 \n",
      "   0.0    16.0    0.0    0.0    0.0 \n"
     ]
    }
   ],
   "source": [
    "import numpy\n",
    "\n",
    "print \"Here cˆ2=16, taking h=1, finding k such that cˆ2tˆ2=1 \"\n",
    "print \"So bendre−schmidits recurrence relation ie u(i)(j+1)=(16tˆ2(u(i−1)(j)+u(i+1)(j))+2(1−16∗tˆ2u(i,j)−u(i)(j−1)\"\n",
    "print \"Now since u(0,t)=0=u(5,t) therefore u(0,i)=0 and u(5,j)=0 and u(x,0)=xˆ2(5−x)\"\n",
    "c = 4\n",
    "h =1 \n",
    "k = (h/c)\n",
    "t = k/h\n",
    "A = numpy.zeros((6,6))\n",
    "print \"Also from 1st derivative (u(i)(j+1)−u(i,j−1))/2k=g(x) and g(x)=0 in this case\"\n",
    "print \"So if j=0 this gives u(i)(1)=1/2∗(u(i−1)(0)+u(i+1)(0))\"\n",
    "for i in range(0,6):\n",
    "    for j in range(1,9):\n",
    "        A[0,i] = 0\n",
    "        A[5,i] = 0;\n",
    "        A[i,1] = (i)**2*(5-i)\n",
    "for i in range(0,4):\n",
    "    A[i+1,2] =1/2*(A[i,1]+A[i+2,1])\n",
    "for i in range(2,5):\n",
    "    for j in range(2,5):\n",
    "        A[i-1,j] = (c*t)**2*(A[i-2,j-1]+A[i,j-1])+2*(1-(c*t)**2)*A[i-1,j-1]-A[i-1,j-2]\n",
    "for i in range(0,5):\n",
    "    for j in range(0,5):\n",
    "        print  \"  \",A[i,j],\n",
    "    print \"\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example 28.8, page no. 734"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Here cˆ2=4, taking h=1, finding k such that cˆ2tˆ2=1 \n",
      "So bendre−schmidits recurrence relation ie u(i)(j+1)=(16tˆ2(u(i−1)(j)+u(i+1)(j))+2(1−16∗tˆ2u(i,j)−u(i)(j−1)\n",
      "Now since u(0,t)=0=u(4,t) therefore u(0,i)=0 and u(4,j)=0 and u(x,0)=x(4−x) \n",
      "Also from 1st derivative (u(i)(j+1)−u(i,j−1))/2 k=g(x) and g(x)=0 in this case\n",
      "So if j=0 this gives u(i)(1)=1/2∗(u(i−1)(0)+u(i+1)(0))\n",
      "   0.0    0.0    0.0    0.0    0.0 \n",
      "   3.0    0.0    -3.0    -6.0    -9.0 \n",
      "   4.0    0.0    -4.0    -8.0    -12.0 \n",
      "   3.0    0.0    -3.0    -6.0    -9.0 \n",
      "   0.0    0.0    0.0    0.0    0.0 \n"
     ]
    }
   ],
   "source": [
    "print \"Here cˆ2=4, taking h=1, finding k such that cˆ2tˆ2=1 \"\n",
    "print \"So bendre−schmidits recurrence relation ie u(i)(j+1)=(16tˆ2(u(i−1)(j)+u(i+1)(j))+2(1−16∗tˆ2u(i,j)−u(i)(j−1)\"\n",
    "print \"Now since u(0,t)=0=u(4,t) therefore u(0,i)=0 and u(4,j)=0 and u(x,0)=x(4−x) \"\n",
    "c = 2\n",
    "h = 1\n",
    "k = (h/c)\n",
    "t = k/h\n",
    "A = numpy.zeros((6,6))\n",
    "print \"Also from 1st derivative (u(i)(j+1)−u(i,j−1))/2 k=g(x) and g(x)=0 in this case\"\n",
    "print \"So if j=0 this gives u(i)(1)=1/2∗(u(i−1)(0)+u(i+1)(0))\"\n",
    "for i in range(0,6):\n",
    "    for j in range(1,9):\n",
    "        A[0,i] = 0\n",
    "        A[4,i] = 0\n",
    "        A[i,0] = (i)*(4-i)\n",
    "for i in range(0,4):\n",
    "    A[i+1,2] = 1/2*(A[i,1]+A[i+2,1])\n",
    "for i in range(2,5):\n",
    "    for j in range(2,5):\n",
    "        A[i-1,j] = (c*t)**2*(A[i-2,j-1]+A[i,j-1])+2*(1-(c*t)**2)*A[i-1,j-1]-A[i-1,j-2]\n",
    "for i in range (0,5):\n",
    "    for j in range(0,5):\n",
    "        print \"  \",A[i,j],\n",
    "    print \"\""
   ]
  }
 ],
 "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.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}