Jeremy Knight Math 639 4/27/11 Algorithm 12.1-Poisson Equation Finite Difference
<Text-field style="Heading 1" layout="Heading 1">PoissonFiniteDiff</Text-field> To approximate the solution to the Poisson equation LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYsLUkmbWZyYWNHRiQ2KC1GIzYlLUklbXN1cEdGJDYlLUkjbWlHRiQ2JVEnJiM5NDg7RicvJSdpdGFsaWNHUSZmYWxzZUYnLyUsbWF0aHZhcmlhbnRHUSdub3JtYWxGJy1GIzYlLUkjbW5HRiQ2JFEiMkYnRjovRjhRJXRydWVGJy9GO1EnaXRhbGljRicvJTFzdXBlcnNjcmlwdHNoaWZ0R1EiMEYnLUY0NiVRInVGJ0ZDRkVGOi1GIzYlLUYxNiUtRjQ2JVEoJiM5NDg7eEYnRkNGRUY9RkdGQ0ZFLyUubGluZXRoaWNrbmVzc0dRIjFGJy8lK2Rlbm9tYWxpZ25HUSdjZW50ZXJGJy8lKW51bWFsaWduR0ZZLyUpYmV2ZWxsZWRHRjktSShtZmVuY2VkR0YkNiQtRiM2Ji1GNDYlUSJ4RidGQ0ZFLUkjbW9HRiQ2LVEiLEYnRjovJSZmZW5jZUdGOS8lKnNlcGFyYXRvckdGRC8lKXN0cmV0Y2h5R0Y5LyUqc3ltbWV0cmljR0Y5LyUobGFyZ2VvcEdGOS8lLm1vdmFibGVsaW1pdHNHRjkvJSdhY2NlbnRHRjkvJSdsc3BhY2VHUSYwLjBlbUYnLyUncnNwYWNlR1EsMC4zMzMzMzMzZW1GJy1GNDYlUSJ5RidGQ0ZFRjpGOi1GYW82LVEiK0YnRjpGZG8vRmdvRjlGaG9Gam9GXHBGXnBGYHAvRmNwUSwwLjIyMjIyMjJlbUYnL0ZmcEZgcS1GLDYoRi4tRiM2JS1GMTYlLUY0NiVRKCYjOTQ4O3lGJ0ZDRkVGPUZHRkNGRUZURldGWkZmbkZobi1GYW82LVEiPUYnRjpGZG9GXnFGaG9Gam9GXHBGXnBGYHAvRmNwUSwwLjI3Nzc3NzhlbUYnL0ZmcEZfci1GNDYlUSJmRidGQ0ZFRmhuRmBvRjo= LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYpLUkjbWlHRiQ2JVEiYUYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1JI21vR0YkNi1RJiZsZXE7RicvRjNRJ25vcm1hbEYnLyUmZmVuY2VHUSZmYWxzZUYnLyUqc2VwYXJhdG9yR0Y9LyUpc3RyZXRjaHlHRj0vJSpzeW1tZXRyaWNHRj0vJShsYXJnZW9wR0Y9LyUubW92YWJsZWxpbWl0c0dGPS8lJ2FjY2VudEdGPS8lJ2xzcGFjZUdRLDAuMjc3Nzc3OGVtRicvJSdyc3BhY2VHRkwtRiw2JVEieEYnRi9GMkY1LUYsNiVRImJGJ0YvRjItRjY2LVEiLEYnRjlGOy9GP0YxRkBGQkZERkZGSC9GS1EmMC4wZW1GJy9GTlEsMC4zMzMzMzMzZW1GJ0Y5 LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYoLUkjbWlHRiQ2JVEiY0YnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1JI21vR0YkNi1RJiZsZXE7RicvRjNRJ25vcm1hbEYnLyUmZmVuY2VHUSZmYWxzZUYnLyUqc2VwYXJhdG9yR0Y9LyUpc3RyZXRjaHlHRj0vJSpzeW1tZXRyaWNHRj0vJShsYXJnZW9wR0Y9LyUubW92YWJsZWxpbWl0c0dGPS8lJ2FjY2VudEdGPS8lJ2xzcGFjZUdRLDAuMjc3Nzc3OGVtRicvJSdyc3BhY2VHRkwtRiw2JVEieUYnRi9GMkY1LUYsNiVRImRGJ0YvRjJGOQ== Subject to the boundary conditions LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYoLUkjbWlHRiQ2JVEidUYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1JKG1mZW5jZWRHRiQ2JC1GIzYmLUYsNiVRInhGJ0YvRjItSSNtb0dGJDYtUSIsRicvRjNRJ25vcm1hbEYnLyUmZmVuY2VHUSZmYWxzZUYnLyUqc2VwYXJhdG9yR0YxLyUpc3RyZXRjaHlHRkUvJSpzeW1tZXRyaWNHRkUvJShsYXJnZW9wR0ZFLyUubW92YWJsZWxpbWl0c0dGRS8lJ2FjY2VudEdGRS8lJ2xzcGFjZUdRJjAuMGVtRicvJSdyc3BhY2VHUSwwLjMzMzMzMzNlbUYnLUYsNiVRInlGJ0YvRjJGQUZBLUY+Ni1RIj1GJ0ZBRkMvRkdGRUZIRkpGTEZORlAvRlNRLDAuMjc3Nzc3OGVtRicvRlZGam4tRiw2JVEiZ0YnRi9GMkY1RkE= if LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2JVEieEYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1JI21vR0YkNi1RIj1GJy9GM1Enbm9ybWFsRicvJSZmZW5jZUdRJmZhbHNlRicvJSpzZXBhcmF0b3JHRj0vJSlzdHJldGNoeUdGPS8lKnN5bW1ldHJpY0dGPS8lKGxhcmdlb3BHRj0vJS5tb3ZhYmxlbGltaXRzR0Y9LyUnYWNjZW50R0Y9LyUnbHNwYWNlR1EsMC4yNzc3Nzc4ZW1GJy8lJ3JzcGFjZUdGTC1GLDYlUSJhRidGL0YyRjk= or LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2JVEieEYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1JI21vR0YkNi1RIj1GJy9GM1Enbm9ybWFsRicvJSZmZW5jZUdRJmZhbHNlRicvJSpzZXBhcmF0b3JHRj0vJSlzdHJldGNoeUdGPS8lKnN5bW1ldHJpY0dGPS8lKGxhcmdlb3BHRj0vJS5tb3ZhYmxlbGltaXRzR0Y9LyUnYWNjZW50R0Y9LyUnbHNwYWNlR1EsMC4yNzc3Nzc4ZW1GJy8lJ3JzcGFjZUdGTC1GLDYlUSJiRidGL0YyRjk= and LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYoLUkjbWlHRiQ2JVEiY0YnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1JI21vR0YkNi1RJiZsZXE7RicvRjNRJ25vcm1hbEYnLyUmZmVuY2VHUSZmYWxzZUYnLyUqc2VwYXJhdG9yR0Y9LyUpc3RyZXRjaHlHRj0vJSpzeW1tZXRyaWNHRj0vJShsYXJnZW9wR0Y9LyUubW92YWJsZWxpbWl0c0dGPS8lJ2FjY2VudEdGPS8lJ2xzcGFjZUdRLDAuMjc3Nzc3OGVtRicvJSdyc3BhY2VHRkwtRiw2JVEieUYnRi9GMkY1LUYsNiVRImRGJ0YvRjJGOQ== and LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYoLUkjbWlHRiQ2JVEidUYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1JKG1mZW5jZWRHRiQ2JC1GIzYmLUYsNiVRInhGJ0YvRjItSSNtb0dGJDYtUSIsRicvRjNRJ25vcm1hbEYnLyUmZmVuY2VHUSZmYWxzZUYnLyUqc2VwYXJhdG9yR0YxLyUpc3RyZXRjaHlHRkUvJSpzeW1tZXRyaWNHRkUvJShsYXJnZW9wR0ZFLyUubW92YWJsZWxpbWl0c0dGRS8lJ2FjY2VudEdGRS8lJ2xzcGFjZUdRJjAuMGVtRicvJSdyc3BhY2VHUSwwLjMzMzMzMzNlbUYnLUYsNiVRInlGJ0YvRjJGQUZBLUY+Ni1RIj1GJ0ZBRkMvRkdGRUZIRkpGTEZORlAvRlNRLDAuMjc3Nzc3OGVtRicvRlZGam4tRiw2JVEiZ0YnRi9GMkY1RkE= if LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2JVEieUYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1JI21vR0YkNi1RIj1GJy9GM1Enbm9ybWFsRicvJSZmZW5jZUdRJmZhbHNlRicvJSpzZXBhcmF0b3JHRj0vJSlzdHJldGNoeUdGPS8lKnN5bW1ldHJpY0dGPS8lKGxhcmdlb3BHRj0vJS5tb3ZhYmxlbGltaXRzR0Y9LyUnYWNjZW50R0Y9LyUnbHNwYWNlR1EsMC4yNzc3Nzc4ZW1GJy8lJ3JzcGFjZUdGTC1GLDYlUSJjRidGL0YyRjk= or LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2JVEieUYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1JI21vR0YkNi1RIj1GJy9GM1Enbm9ybWFsRicvJSZmZW5jZUdRJmZhbHNlRicvJSpzZXBhcmF0b3JHRj0vJSlzdHJldGNoeUdGPS8lKnN5bW1ldHJpY0dGPS8lKGxhcmdlb3BHRj0vJS5tb3ZhYmxlbGltaXRzR0Y9LyUnYWNjZW50R0Y9LyUnbHNwYWNlR1EsMC4yNzc3Nzc4ZW1GJy8lJ3JzcGFjZUdGTC1GLDYlUSJkRidGL0YyRjk= and LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYoLUkjbWlHRiQ2JVEiYUYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1JI21vR0YkNi1RJiZsZXE7RicvRjNRJ25vcm1hbEYnLyUmZmVuY2VHUSZmYWxzZUYnLyUqc2VwYXJhdG9yR0Y9LyUpc3RyZXRjaHlHRj0vJSpzeW1tZXRyaWNHRj0vJShsYXJnZW9wR0Y9LyUubW92YWJsZWxpbWl0c0dGPS8lJ2FjY2VudEdGPS8lJ2xzcGFjZUdRLDAuMjc3Nzc3OGVtRicvJSdyc3BhY2VHRkwtRiw2JVEieEYnRi9GMkY1LUYsNiVRImJGJ0YvRjJGOQ== PoissonFiniteDiff:=proc(f,g,a,b,c,d,m::integer,n::integer,TOL,N) #f,g = functions defined by Poisson problem #a,b,c,d = endpoints #m>=3 , n>=3 are number of grid lines for mesh #TOL=tolerance, N= maximum iterations local i,j,h,k,L,lambda,mu,w,x,y,z,NORM: #option trace: x:=Array(1..n): y:=Array(1..n): h:=(b-a)/n: k:=(d-c)/m: for i from 1 to n-1 do x[i]:=a+i*h: end do: for j from 1 to m-1 do y[j]:=c+j*k: end do: w:=Array(1..n-1,1..m-1): lambda:=h^2/k^2: mu:=2*(1+lambda): L:=1: ##### Main While Loop ###### while L<=N do z:=(-1*h^2*f(x[1],y[m-1]) + g(a,y[m-1]) + lambda*g(x[1],d) + lambda*w[1,m-2] + w[2,m-1])/mu: NORM:=abs(z-w[1,m-1]): w[1,m-1]:=z: #8 for i from 2 to n-2 do z:=(-1*h^2*f(x[i],y[m-1]) + lambda*g(x[i],d) + w[i-1,m-1] + w[i+1,m-1] + lambda*w[i,m-2])/mu: if abs(w[i,m-1]-z)>NORM then NORM:=abs(w[i,m-1]-z): end if: w[i,m-1]:=z: end do: #9 z:=(-1*h^2*f(x[n-1],y[m-1]) + g(b,y[m-1]) + lambda*g(x[n-1],d) + w[n-2,m-1] + lambda*w[n-1,m-2])/mu: if abs(w[n-1,m-1]-z)>NORM then NORM:=abs(w[n-1,m-1]-z): end if: w[n-1,m-1]:=z: #10 for j from m-2 to 2 by -1 do #print("*********** marker 1************"): #11 z:=(-1*h^2*f(x[1],y[j]) + g(a,y[j]) + lambda*w[1,j+1] + lambda*w[1,j-1] + w[2,j])/mu: if abs(w[1,j]-z)>NORM then NORM:=abs(w[1,j]-z): end if: w[1,j]:=z: #12 for i from 2 to n-2 do z:=(-1*h^2*f(x[i],y[j]) + w[i-1,j] + lambda*w[i,j+1] + w[i+1,j] + lambda*w[i,j-1])/mu: if abs(w[i,j]-z) > NORM then NORM:= abs(w[i,j]-z): end if: w[i,j]:=z: end do: #13 z:=(-1*h^2*f(x[n-1],y[j]) + g(b,y[j]) + w[n-2,j] + lambda*w[n-1,j+1] + lambda*w[n-1,j-1])/mu: if abs(w[n-1,j]-z)>NORM then NORM:=abs(w[n-1,j]-z): end if: w[n-1,j]:=z: end do: # print("*********** marker2 ************"): #14 z:=(-1*h^2*f(x[1],y[1]) + g(a,y[1]) + lambda*g(x[1],c) + lambda*w[1,2] + w[2,1])/mu: if abs(w[1,1]-z)>NORM then NORM:=abs(w[1,1]-z): end if: w[1,1]:=z: #15 for i from 2 to n-2 do z:=(-1*h^2*f(x[i],y[1]) + lambda*g(x[i],c) + w[i-1,1] + lambda*w[i,2] + w[i+1,1])/mu: if abs(w[i,1]-z)>NORM then NORM:=abs(w[i,1]-z): end if: w[i,1]:=z: end do: #16 z:=(-1*h^2*f(x[n-1],y[1]) + g(b,y[1]) + lambda*g(x[n-1],c) + w[n-2,1] + lambda*w[n-1,2])/mu: if abs(w[n-1,1]-z)>NORM then NORM:=abs(w[n-1,1]-z): end if: w[n-1,1]:=z: #17 if NORM<=TOL then #18 for i from 1 to n-1 do for j from 1 to m-1 do print(x[i],y[j],w[i,j]): end do: end do: #19 break: end if: #20 L:=L+1: end do: ##### End main While Loop ##### if L>N then print("maximum number of iterations exceeded"): end if: print("*********** marker3 ************"): print(x,y,w); end proc:
JSFH JSFH
<Text-field style="Heading 1" layout="Heading 1">Example pg 722</Text-field> LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYjLUkjbWlHRiQ2I1EhRic= Digits:=20: f:=(x,y)->x*exp(y); g:=proc(x,y) local u: if y>=0 and y<=1 then if x=0 then u:=0: end if: if x=2 then u:=evalf(2.0*exp(y)): end if: end if: if x>=0 and x<=2 then if y=0 then u:=x: end if: if y=1 then u:=evalf(exp(1)*x): end if: end if: return u: end proc: Zio2JEkieEc2IkkieUdGJUYlNiRJKW9wZXJhdG9yR0YlSSZhcnJvd0dGJUYlKiY5JCIiIi1JJGV4cEdGJTYjOSVGLEYlRiVGJQ== interface(displayprecision=4): PoissonFiniteDiff(f,g,0.0,2.0,0.0,1.0,5,6,10^(-10),100); NiUkIjVMTExMTExMTExMISM/JCI1KysrKysrKysrP0YlJCI1aXJSIzRDPFlFMiVGJQ== NiUkIjVMTExMTExMTExMISM/JCI1KysrKysrKysrU0YlJCI1Jz0iekswZEMkWyhcRiU= NiUkIjVMTExMTExMTExMISM/JCI1KysrKysrKysrZ0YlJCI1cSVRM2gjKnpnZjInRiU= NiUkIjVMTExMTExMTExMISM/JCI1KysrKysrKysrISlGJSQiNTlOIypSIjRscStVKEYl NiUkIjVtbW1tbW1tbW1tISM/JCI1KysrKysrKysrP0YlJCI1XlknZkElKlxQXzkpRiU= NiUkIjVtbW1tbW1tbW1tISM/JCI1KysrKysrKysrU0YlJCI1c0QkekZbdnYmXCoqRiU= NiUkIjVtbW1tbW1tbW1tISM/JCI1KysrKysrKysrZ0YlJCI1a2IqXF15PCQ9OjchIz4= NiUkIjVtbW1tbW1tbW1tISM/JCI1KysrKysrKysrISlGJSQiNSgpby9yS1AmM1NbIiEjPg== NiUkIjUqKioqKioqKioqKioqKioqKioqKiEjPyQiNSsrKysrKysrKz9GJSQiNUEjR0hFLkBtPEEiISM+ NiUkIjUqKioqKioqKioqKioqKioqKioqKiEjPyQiNSsrKysrKysrK1NGJSQiNWZUR3Y9RFpTI1wiISM+ NiUkIjUqKioqKioqKioqKioqKioqKioqKiEjPyQiNSsrKysrKysrK2dGJSQiNW5pPCVmUHJVRiM9ISM+ NiUkIjUqKioqKioqKioqKioqKioqKioqKiEjPyQiNSsrKysrKysrKyEpRiUkIjUmKUg4RUUtRipmQSMhIz4= NiUkIjVMTExMTExMTEw4ISM+JCI1KysrKysrKysrPyEjPyQiNT1TJG9Cb29qKkc7RiU= NiUkIjVMTExMTExMTEw4ISM+JCI1KysrKysrKysrUyEjPyQiNSJbbHFqTEl5KCopPkYl NiUkIjVMTExMTExMTEw4ISM+JCI1KysrKysrKysrZyEjPyQiNThNXVFGaGxBSUNGJQ== NiUkIjVMTExMTExMTEw4ISM+JCI1KysrKysrKysrISkhIz8kIjVZInBaT1tER3onSEYl NiUkIjVtbW1tbW1tbW07ISM+JCI1KysrKysrKysrPyEjPyQiNUQxInBQRktVZy4jRiU= NiUkIjVtbW1tbW1tbW07ISM+JCI1KysrKysrKysrUyEjPyQiNTAvOz9mKVFlcFsjRiU= NiUkIjVtbW1tbW1tbW07ISM+JCI1KysrKysrKysrZyEjPyQiNXNTWF1kLWJdUElGJQ== NiUkIjVtbW1tbW1tbW07ISM+JCI1KysrKysrKysrISkhIz8kIjUiPVRJKnltU3M0UEYl UUEqKioqKioqKioqKn5tYXJrZXIzfioqKioqKioqKioqKjYi NiUtSSZBcnJheUclKnByb3RlY3RlZEc2Iy9JJCVpZEc2IiIqbycpUVsjLUYkNiMvRigiKlMoKVFbIy1GJDYjL0YoIipfMlJbIw==
JSFH TTdSMApJNlJUQUJMRV9TQVZFLzI0ODM4ODY2OFgqJSlhbnl0aGluZ0c2IjYiW2dsISElISEhIiciJyQiNUxMTExMTExMTEwhIz8kIjVtbW1tbW1tbW1tRikkIjUqKioqKioqKioqKioqKioqKioqKkYpJCI1TExMTExMTExMOCEjPiQiNW1tbW1tbW1tbTtGMCIiIUYmTTdSMApJNlJUQUJMRV9TQVZFLzI0ODM4ODc0MFgqJSlhbnl0aGluZ0c2IjYiW2dsISElISEhIiciJyQiNSsrKysrKysrKz8hIz8kIjUrKysrKysrKytTRikkIjUrKysrKysrKytnRikkIjUrKysrKysrKyshKUYpIiIhRjBGJg==TTdSMApJNlJUQUJMRV9TQVZFLzI0ODM5MDc1MlgsJSlhbnl0aGluZ0c2IjYiW2dsISElISEhIzUiJiIlJCI1aXJSIzRDPFlFMiUhIz8kIjVeWSdmQSUqXFBfOSlGKSQiNUEjR0hFLkBtPEEiISM+JCI1PVMkb0Jvb2oqRztGLiQiNUQxInBQRktVZy4jRi4kIjUnPSJ6SzBkQyRbKFxGKSQiNXNEJHpGW3Z2JlwqKkYpJCI1ZlRHdj1EWlMjXCJGLiQiNSJbbHFqTEl5KCopPkYuJCI1MC87P2YpUWVwWyNGLiQiNXElUTNoIyp6Z2YyJ0YpJCI1a2IqXF15PCQ9OjdGLiQiNW5pPCVmUHJVRiM9Ri4kIjU4TV1RRmhsQUlDRi4kIjVzU1hdZC1iXVBJRi4kIjU5TiMqUiI0bHErVShGKSQiNSgpby9yS1AmM1NbIkYuJCI1JilIOEVFLUYqZkEjRi4kIjVZInBaT1tER3onSEYuJCI1Ij1USSp5bVNzNFBGLkYm