Poisson equation
by Reinaldo Baretti Machín
www.oocities.org/serienumerica

c FORTRAN CODE
c   POISSON equation in TWO DIMENSIONS ,see Burden &
c     Fairies p 575 Chap 11 , Schaum Numerical Analysis Chapter 29
c  f(x,y) is the source or density/ phi(x,y) is the exact solution
c   Del**2 (Phi(x,y)) =4
      dimension u(0:100,0:100), u1(0:100,0:100),delta(0:100,0:100)
      f(x,y)=4.
      phi(x,y)=x**2 + y**2
      uleft(x1,y)=phi(x1,y)
      ubot(x,y1)=phi(x,y1)
      uright(x2,y)=phi(x2,y)
      utop(x,y2)=phi(x,y2)
      x1=0.
      x2=1.
      y1=0.
      y2=1.
      ny=4
      nx=4
      dx=(x2-x1)/float(nx)
      dy=(y2-y1)/float(ny)
      ainv=1./(2./dx**2+2./dy**2)
      niter=20
c     boundary conditions are imposed
      do 10 i=0,nx
      x=dx*float(i)
      u(i,0)=ubot(x,y1)
      u1(i,0)= ubot(x,y1)
      u(i,ny)=utop(x,y2)
      u1(i,ny)=utop(x,y2)
10    continue
      do 15 j=0,ny
      y=dy*float(j)
      u(0,j)=uleft(x1,y)
      u1(0,j)=uleft(x1,y)
      u(nx,j)=uright(x2,y)
      u1(nx,j)= uright(x2,y)
15    continue
c initial solution for u is taken as some average
      do 20 i=1,nx-1
      do 20 j=1,ny-1
      x=dx*float(i)
      y=dy*float(j)
      u(i,j)=.25*( uleft(x1,y)+ uright(x2,y)+ ubot(x,y2)+
     $ utop(x,y2))
20    continue
c   iteration procedure is used for the midpoint u1 (x,y)
      do 100 iter=1,niter
      do 30 i=1,nx-1
      do 30 j=1,ny-1
      x=dx*float(i)
      y=dy*float(j)
      u1(i,j)=-f(x,y)+(1./dx**2)*(u(i+1,j)+u(i-1,j))+(1./dy**2)*
     $ (u(i,j+1) +u(i,j-1))
      u1(i,j)=ainv*u1(i,j)
30    continue
c introduce tolerance criteria here e.g. abs (u1(i,j)-u(i,j)).le.1.0e-3
c   u(i,j) is redefined
      do 40 i=1,nx-1
      do 40 j=1,ny-1
      u(i,j)=u1(i,j)
40    continue
c      print 120,u(nx/2,ny/2),phi(x2/2.,y2/2.)
      if(iter.eq.niter)then
c  prints results for interior points
      do 50 j=1,ny-1
      do 50 i=1,nx-1
      x=dx*float(i)
      y=dy*float(j)
      print 110,x,y,u(i,j),phi(x,y)
c      print*,' '
50    continue
      endif
100   continue
110   format('x,y,u,phi=',4(2x,e11.4))
120   format('u(nx/2,ny/2),phi(x2/2./y2/2.)=',2(4x,e10.3))
      stop
      end
      
      RUN
      x,y,u,phi=   0.2500E+00   0.2500E+00   0.1253E+00   0.1250E+00
x,y,u,phi=   0.5000E+00   0.2500E+00   0.3130E+00   0.3125E+00
x,y,u,phi=   0.7500E+00   0.2500E+00   0.6253E+00   0.6250E+00
x,y,u,phi=   0.2500E+00   0.5000E+00   0.3130E+00   0.3125E+00
x,y,u,phi=   0.5000E+00   0.5000E+00   0.5007E+00   0.5000E+00
x,y,u,phi=   0.7500E+00   0.5000E+00   0.8130E+00   0.8125E+00
x,y,u,phi=   0.2500E+00   0.7500E+00   0.6253E+00   0.6250E+00
x,y,u,phi=   0.5000E+00   0.7500E+00   0.8130E+00   0.8125E+00
x,y,u,phi=   0.7500E+00   0.7500E+00   0.1125E+01   0.1125E+01