URL www.oocities.org/reibaretti2004 www.oocities.org/serienumerica
c Poisson equation del ^2 U = f(x,y)= x *exp(y) / solution is
c u= x *exp(y) see Numerical Analysis Burden & Faires Chapter 11
dimension u(0:100,0:100),uiter(0:100,0:100)
f(x,y)= x*exp(y)
x1=0.
x2=2.
y1=0.
y2=1.
nx=5
ny=nx
dx=(x2-x1)/float(nx)
dy=(y2-y1)/float(ny)
iter=20
c boundary conditions are set (are given)
do 10 j=0,ny
y=dy*float(j)
u(0,j)=0.
u(nx,j)=x2*exp(y)
uiter(0,j)=0.
uiter(nx,j)=x2*exp(y)
10 continue
do 20 i=0,nx
x=dx*float(i)
u(i,0)=x
u(i,ny)=x*exp(y2)
uiter(i,0)=x
uiter(i,ny)=x*exp(y2)
20 continue
c intial values are given for interior points
do 30 i=1,nx-1
do 30 j=1,ny-1
x=x1+dx*float(i)
y=y1+dy*float(j)
u(i,j)=.25*(u(i,0)+u(i,ny)+u(0,j)+u(nx,j))
30 continue
c start iterations
denom=2./dx**2+2./dy**2
do 40 it=1,iter
do 50 i=1,nx-1
do 50 j=1,ny-1
x=x1+dx*float(i)
y=y1+dy*float(j)
c**************corregir
uiter(i,j)=(u(i+1,j)+u(i-1,j))/dx**2+ (u(i,j+1)+u(i,j-1))/dy**2
& -f(x,y)
uiter(i,j)=uiter(i,j)/denom
50 continue
c relabel uiter and u
do 60 i=1,nx-1
do 60 j=1,ny-1
u(i,j)=uiter(i,j)
60 continue
40 continue
c print part
do 70 j=1,ny,1
do 70 i=1,nx,1
x=dx*float(i)
y=dy*float(j)
print 100,x,y,u(i,j),f(x,y)
70 continue
100 format(1x,'x,y,u, uexact=',4(2x,e12.5))
stop
end
RUN with 20 iterations
x,y,u, uexact= 0.40000E+00 0.20000E+00 0.48961E+00 0.48856E+00
x,y,u, uexact= 0.80000E+00 0.20000E+00 0.97879E+00 0.97712E+00
x,y,u, uexact= 0.12000E+01 0.20000E+00 0.14669E+01 0.14657E+01
x,y,u, uexact= 0.16000E+01 0.20000E+00 0.19551E+01 0.19542E+01
x,y,u, uexact= 0.20000E+01 0.20000E+00 0.24428E+01 0.24428E+01
x,y,u, uexact= 0.40000E+00 0.40000E+00 0.59863E+00 0.59673E+00
x,y,u, uexact= 0.80000E+00 0.40000E+00 0.11958E+01 0.11935E+01
x,y,u, uexact= 0.12000E+01 0.40000E+00 0.17925E+01 0.17902E+01
x,y,u, uexact= 0.16000E+01 0.40000E+00 0.23880E+01 0.23869E+01
x,y,u, uexact= 0.20000E+01 0.40000E+00 0.29836E+01 0.29836E+01
x,y,u, uexact= 0.40000E+00 0.60000E+00 0.73056E+00 0.72885E+00
x,y,u, uexact= 0.80000E+00 0.60000E+00 0.14604E+01 0.14577E+01
x,y,u, uexact= 0.12000E+01 0.60000E+00 0.21885E+01 0.21865E+01
x,y,u, uexact= 0.16000E+01 0.60000E+00 0.29168E+01 0.29154E+01
x,y,u, uexact= 0.20000E+01 0.60000E+00 0.36442E+01 0.36442E+01
x,y,u, uexact= 0.40000E+00 0.80000E+00 0.89143E+00 0.89022E+00
x,y,u, uexact= 0.80000E+00 0.80000E+00 0.17820E+01 0.17804E+01
x,y,u, uexact= 0.12000E+01 0.80000E+00 0.26722E+01 0.26706E+01
x,y,u, uexact= 0.16000E+01 0.80000E+00 0.35617E+01 0.35609E+01
x,y,u, uexact= 0.20000E+01 0.80000E+00 0.44511E+01 0.44511E+01
x,y,u, uexact= 0.40000E+00 0.10000E+01 0.10873E+01 0.10873E+01
x,y,u, uexact= 0.80000E+00 0.10000E+01 0.21746E+01 0.21746E+01
x,y,u, uexact= 0.12000E+01 0.10000E+01 0.32619E+01 0.32619E+01
x,y,u, uexact= 0.16000E+01 0.10000E+01 0.43493E+01 0.43493E+01
x,y,u, uexact= 0.20000E+01 0.10000E+01 0.54366E+01 0.54366E+01