c fishpk40 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 poistg to
c solve the equation
c
c (1/cos(x))(d/dx)(cos(x)(du/dx)) + (d/dy)(du/dy) =
c
c 2*y**2*(6-y**2)*sin(x)
c
c on the rectangle -pi/2 .lt. x .lt. pi/2 and
c 0 .lt. y .lt. 1 with the boundary conditions
c
c (du/dx) (-pi/2,y) = (du/dx)(pi/2,y) = 0 , 0 .le. y .le. 1 (2)
c
c u(x,0) = 0 (3)
c -pi/2 .le. x .le. pi/2
c (du/dy)(x,1) = 4sin(x) (4)
c
c using finite differences on a staggered grid with
c deltax (= dx) = pi/40 and deltay (= dy) = 1/20 .
c to set up the finite difference equations we define
c the grid points
c
c x(i) = -pi/2 + (i-0.5)dx i=1,2,...,40
c
c y(j) = (j-o.5)dy j=1,2,...,20
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,...,40 and j=1,2,...,20.
c hence, in the program m = 40 and n = 20. 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) (5)
c
c where s = (dy/dx)**2, and for i=2,3,...,39
c
c a(i) = s*cos(x(i)-dx/2)
c
c b(i) = -s*(cos(x(i)-dx/2)+cos(x(i)+dx/2))
c
c c(i) = s*cos(x(i)+dx/2)
c
c f(i,j) = 2dy**2*y(j)**2*(6-y(j)**2)*sin(x(i)) , j=1,2,...,19.
c
c to obtain equations for i = 1, we replace equation (2)
c by the second order approximation
c
c (v(1,j)-v(0,j))/dx = 0
c
c and use this equation to eliminate v(0,j) in equation (5)
c to 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) = -s*(cos(x(1)-dx/2)+cos(x(1)+dx/2))
c
c c(1) = -b(1)
c
c for completeness, we set a(1) = 0.
c to obtain equations for i = 40, we replace the derivative
c in equation (2) at x=pi/2 in a similar fashion, use this
c equation to eliminate the virtual unknown v(41,j) in equation
c (5) and arrive at the equation
c
c a(40)v(39,j) + b(40)v(40,j)
c
c + v(40,j-1) - 2v(40,j) + v(40,j+1) = f(40,j)
c
c where
c
c a(40) = -b(40) = -s*(cos(x(40)-dx/2)+cos(x(40)+dx/2))
c
c for completeness, we set c(40) = 0. hence, in the
c program mperod = 1.
c for j = 1, we replace equation (3) by the second order
c approximation
c
c (v(i,0) + v(i,1))/2 = 0
c
c to arrive at the condition
c
c v(i,0) = -v(i,1) .
c
c for j = 20, we replace equation (4) by the second order
c approximation
c
c (v(i,21) - v(i,20))/dy = 4*sin(x)
c
c and combine this equation with equation (5) to arrive at
c the equation
c
c a(i)v(i-1,20) + b(i)v(i,20) + c(i)v(i+1,20)
c
c + v(i,19) - 2v(i,20) + v(i,21) = f(i,20)
c
c where
c
c v(i,21) = v(i,20) and
c
c f(i,20) = 2*dy**2*y(j)**2*(6-y(j)**2)*sin(x(i)) - 4*dy*sin(x(i))
c
c hence, in the program nperod = 2 .
c the exact solution to this problem is
c
c u(x,y) = y**4*cos(x) .
c
dimension f(42,20) ,a(40) ,b(40) ,c(40) ,
1 w(600) ,x(40) ,y(20)
c
c from dimension statement we get value of idimf = 42. also
c note that w has been dimensioned
c 9m + 4n + m(int(log2(n))) = 360 + 80 + 160 = 600 .
c
idimf = 42
mperod = 1
m = 40
pi = pimach(dum)
dx = pi/float(m)
nperod = 2
n = 20
dy = 1./float(n)
c
c generate and store grid points for computation.
c
do 101 i=1,m
x(i) = -pi/2.+(float(i)-0.5)*dx
101 continue
do 102 j=1,n
y(j) = (float(j)-0.5)*dy
102 continue
c
c generate coefficients .
c
s = (dy/dx)**2
a(1) = 0.
b(1) = -s*cos(-pi/2.+dx)/cos(x(1))
c(1) = -b(1)
do 103 i=2,m
a(i) = s*cos(x(i)-dx/2.)/cos(x(i))
c(i) = s*cos(x(i)+dx/2.)/cos(x(i))
b(i) = -(a(i)+c(i))
103 continue
a(40) = -b(40)
c(40) = 0.
c
c generate right side of equation.
c
do 105 i=1,m
do 104 j=1,n
f(i,j) = 2.*dy**2*y(j)**2*(6.-y(j)**2)*sin(x(i))
104 continue
105 continue
do 106 i=1,m
f(i,n) = f(i,n)-4.*dy*sin(x(i))
106 continue
call poistg (nperod,n,mperod,m,a,b,c,idimf,f,ierror,w)
c
c compute discretization error. the exact solution is
c
c u(x,y) = y**4*sin(x)
c
err = 0.
do 108 i=1,m
do 107 j=1,n
t = abs(f(i,j)-y(j)**4*sin(x(i)))
if (t .gt. err) err = t
107 continue
108 continue
print 1001 , ierror,err,w(1)
stop
c
1001 format (1h1,20x,25hsubroutine poistg example///
1 10x,46hthe output from the ncar control data 7600 was//
2 32x,10hierror = 0/
3 18x,34hdiscretization error = 5.64171e-04/
4 12x,32hrequired length of w array = 560//
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 discretization error. the exact solution is