c fishpk38 from portlib 12/30/83
c
c * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
c * *
c * f i s h p a k *
c * *
c * *
c * a package of fortran subprograms for the solution of *
c * *
c * separable elliptic partial differential equations *
c * *
c * (version 3.1 , october 1980) *
c * *
c * by *
c * *
c * john adams, paul swarztrauber and roland sweet *
c * *
c * of *
c * *
c * the national center for atmospheric research *
c * *
c * boulder, colorado (80307) u.s.a. *
c * *
c * which is sponsored by *
c * *
c * the national science foundation *
c * *
c * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
c
c program to illustrate the use of subroutine genbun to
c solve the equation
c
c (1+x)**2*(d/dx)(du/dx) - 2(1+x)(du/dx) + (d/dy)(du/dy)
c
c = 3(1+x)**4*sin(y) (1)
c
c on the rectangle 0 .lt. x .lt. 1 and -pi .lt. y .lt. pi
c with the boundary conditions
c
c (du/dx)(0,y) = 4sin(y) (2)
c -pi .le. y .le. pi
c u(1,y) = 16sin(y) (3)
c
c and with u periodic in y using finite differences on a
c grid with deltax (= dx) = 1/20 and deltay (= dy) = pi/20.
c to set up the finite difference equations we define
c the grid points
c
c x(i) = (i-1)dx i=1,2,...,21
c
c y(j) = -pi + (j-1)dy j=1,2,...,41
c
c and let v(i,j) be an approximation to u(x(i),y(j)).
c numbering the grid points in this fashion gives the set
c of unknowns as v(i,j) for i=1,2,...,20 and j=1,2,...,40.
c hence, in the program m = 20 and n = 40. at the interior
c grid point (x(i),y(j)), we replace all derivatives in
c equation (1) by second order central finite differences,
c multiply by dy**2, and collect coefficients of v(i,j) to
c get the finite difference equation
c
c a(i)v(i-1,j) + b(i)v(i,j) + c(i)v(i+1,j)
c
c + v(i,j-1) - 2v(i,j) + v(i,j+1) = f(i,j) (4)
c
c where s = (dy/dx)**2, and for i=2,3,...,19
c
c a(i) = (1+x(i))**2*s + (1+x(i))*s*dx
c
c b(i) = -2(1+x(i))**2*s
c
c c(i) = (1+x(i))**2*s - (1+x(i))*s*dx
c
c f(i,j) = 3(1+x(i))**4*dy**2*sin(y(j)) for j=1,2,...,40.
c
c to obtain equations for i = 1, we replace the
c derivative in equation (2) by a second order central
c finite difference approximation, use this equation to
c eliminate the virtual unknown v(0,j) in equation (4)
c and arrive at the equation
c
c b(1)v(1,j) + c(1)v(2,j) + v(1,j-1) - 2v(1,j) + v(1,j+1)
c
c = f(1,j)
c
c where
c
c b(1) = -2s , c(1) = 2s
c
c f(1,j) = (11+8/dx)*dy**2*sin(y(j)), j=1,2,...,40.
c
c for completeness, we set a(1) = 0.
c to obtain equations for i = 20, we incorporate
c equation (3) into equation (4) by setting
c
c v(21,j) = 16sin(y(j))
c
c and arrive at the equation
c
c a(20)v(19,j) + b(20)v(20,j)
c
c + v(20,j-1) - 2v(20,j) + v(20,j+1) = f(20,j)
c
c where
c
c a(20) = (1+x(20))**2*s + (1+x(20))*s*dx
c
c b(20) = -2*(1+x(20))**2*s
c
c f(20,j) = (3(1+x(20))**4*dy**2 - 16(1+x(20))**2*s
c + 16(1+x(20))*s*dx)*sin(y(j))
c
c for j=1,2,...,40.
c
c for completeness, we set c(20) = 0. hence, in the
c program mperod = 1.
c the periodicity condition on u gives the conditions
c
c v(i,0) = v(i,40) and v(i,41) = v(i,1) for i=1,2,...,20.
c
c hence, in the program nperod = 0.
c
c the exact solution to this problem is
c
c u(x,y) = (1+x)**4*sin(y) .
c
dimension f(22,40) ,a(20) ,b(20) ,c(20) ,
1 x(21) ,y(41) ,w(380)
c
c from the dimension statement we get that idimf = 22 and that w
c has been dimensioned
c
c 4n + (10+int(log2(n)))m = 4*20 + (10+5)*20 = 380 .
c
idimf = 22
m = 20
mp1 = m+1
mperod = 1
dx = 0.05
n = 40
nperod = 0
pi = pimach(dum)
dy = pi/20.
c
c generate grid points for later use.
c
do 101 i=1,mp1
x(i) = float(i-1)*dx
101 continue
do 102 j=1,n
y(j) = -pi+float(j-1)*dy
102 continue
c
c generate coefficients.
c
s = (dy/dx)**2
do 103 i=2,19
t = 1.+x(i)
tsq = t**2
a(i) = (tsq+t*dx)*s
b(i) = -2.*tsq*s
c(i) = (tsq-t*dx)*s
103 continue
a(1) = 0.
b(1) = -2.*s
c(1) = -b(1)
b(20) = -2.*s*(1.+x(20))**2
a(20) = -b(20)/2.+(1.+x(20))*dx*s
c(20) = 0.
c
c generate right side.
c
do 105 i=2,19
do 104 j=1,n
f(i,j) = 3.*(1.+x(i))**4*dy**2*sin(y(j))
104 continue
105 continue
t = 1.+x(20)
tsq = t**2
t4 = tsq**2
do 106 j=1,n
f(1,j) = (11.+8./dx)*dy**2*sin(y(j))
f(20,j) = (3.*t4*dy**2-16.*tsq*s+16.*t*s*dx)*sin(y(j))
106 continue
call genbun (nperod,n,mperod,m,a,b,c,idimf,f,ierror,w)
c
c compute discretiazation error. the exact solution is
c
c u(x,y) = (1+x)**4*sin(y) .
c
err = 0.
do 108 i=1,m
do 107 j=1,n
t = abs(f(i,j)-(1.+x(i))**4*sin(y(j)))
if (t .gt. err) err = t
107 continue
108 continue
print 1001 , ierror,err,w(1)
stop
c
1001 format (1h1,20x,25hsubroutine genbun example///
1 10x,46hthe output from the ncar control data 7600 was//
2 32x,10hierror = 0/
3 18x,34hdiscretization error = 9.64063e-03/
4 12x,32hrequired length of w array = 380//
5 10x,32hthe output from your computer is//
6 32x,8hierror =,i2/18x,22hdiscretization error =,e12.5/
7 12x,28hrequired length of w array =,f4.0)
c
end
c compute discretiazation error. the exact solution is