{
"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
}