c fishpk43 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 this program illustrates the use of subroutine cblktr which is c the complex version of blktri. the program solves the equation c c .5/s*(d/ds)(.5/s*du/ds)+.5/t*(d/dt)(.5/t*du/dt)-sqrt(-1)*u c (1) c = 15/4*s*t*(s**4+t**4)-sqrt(-1)*(s*t)**5 c c on the rectangle 0 .lt. s .lt. 1 and 0 .lt. t .lt. 1 c with the boundary conditions c c u(0,t) = 0 c 0 .le. t .le. 1 c u(1,t) = t**5 c c and c c u(s,0) = 0 c 0 .le. s .le. 1 c u(s,1) = s**5 c c the exact solution of this problem is u(s,t) = (s*t)**5 c c define the integers m = 50 and n = 63. then define the c grid increments deltas = 1/(m+1) and deltat = 1/(n+1). c c the grid is then given by s(i) = i*deltas for i = 1,...,m c and t(j) = j*deltat for j = 1,...,n. c c the approximate solution is given as the solution to c the following finite difference approximation of equation (1). c c .5/(s(i)*deltas)*((u(i+1,j)-u(i,j))/(2*s(i+.5)*deltas) c -(u(i,j)-u(i-1,j))/(2*s(i-.5)*deltas)) c +.5/(t(i)*deltat)*((u(i,j+1)-u(i,j))/(2*t(i+.5)*deltat) (2) c -(u(i,j)-u(i,j-1))/(2*t(i-.5)*deltat)) c -sqrt(-1)*u(i,j) c = 15/4*s(i)*t(j)*(s(i)**4+t(j)**4) c -sqrt(-1)*(s(i)*t(j))**5 c c where s(i+.5) = .5*(s(i+1)+s(i)) c s(i-.5) = .5*(s(i)+s(i-1)) c t(i+.5) = .5*(t(i+1)+t(i)) c t(i-.5) = .5*(t(i)+t(i-1)) c c the approach is to write equation (2) in the form c c am(i)*u(i-1,j)+bm(i,j)*u(i,j)+cm(i)*u(i+1,j) c +an(j)*u(i,j-1)+bn(j)*u(i,j)+cn(j)*u(i,j+1) (3) c = y(i,j) c c and then call subroutine cblktr to determine u(i,j) c c c c dimension y(75,105) ,am(75) ,bm(75) ,cm(75) , 1 an(105) ,bn(105) ,cn(105) ,s(75) , 2 t(105) ,w(1123) complex y ,am ,bm ,cm c iflg = 0 np = 1 n = 63 mp = 1 m = 50 idimy = 75 c c generate and store grid points for the purpose of computing the c coefficients and the array y. c deltas = 1./float(m+1) do 101 i=1,m s(i) = float(i)*deltas 101 continue deltat = 1./float(n+1) do 102 j=1,n t(j) = float(j)*deltat 102 continue c c compute the coefficients am,bm,cm corresponding to the s direction c hds = deltas/2. tds = deltas+deltas do 103 i=1,m temp1 = 1./(s(i)*tds) temp2 = 1./((s(i)-hds)*tds) temp3 = 1./((s(i)+hds)*tds) am(i) = cmplx(temp1*temp2,0.) cm(i) = cmplx(temp1*temp3,0.) bm(i) = -(am(i)+cm(i))-(0.,1.) 103 continue c c compute the coefficients an,bn,cn corresponding to the t direction c hdt = deltat/2. tdt = deltat+deltat do 104 j=1,n temp1 = 1./(t(j)*tdt) temp2 = 1./((t(j)-hdt)*tdt) temp3 = 1./((t(j)+hdt)*tdt) an(j) = temp1*temp2 cn(j) = temp1*temp3 bn(j) = -(an(j)+cn(j)) 104 continue c c compute right side of equation c do 106 j=1,n do 105 i=1,m y(i,j) = 3.75*s(i)*t(j)*(s(i)**4+t(j)**4)- 1 (0.,1.)*(s(i)*t(j))**5 105 continue 106 continue c c the nonzero boundary conditions enter the linear system via c the right side y(i,j). if the equations (3) given above c are evaluated at i=m and j=1,...,n then the term cm(m)*u(m+1,j) c is known from the boundary condition to be cm(m)*t(j)**5. c therefore this term can be included in the right side y(m,j). c the same analysis applies at j=n and i=1,..,m. note that the c corner at j=n,i=m includes contributions from both boundaries. c do 107 j=1,n y(m,j) = y(m,j)-cm(m)*t(j)**5 107 continue do 108 i=1,m y(i,n) = y(i,n)-cn(n)*s(i)**5 108 continue c 109 call cblktr (iflg,np,n,an,bn,cn,mp,m,am,bm,cm,idimy,y,ierror,w) iflg = iflg+1 if (iflg-1) 109,109,110 c c compute discretization error c 110 err = 0. do 112 j=1,n do 111 i=1,m z = cabs(y(i,j)-(s(i)*t(j))**5) if (z .gt. err) err = z 111 continue 112 continue iw = int(w(1)) print 1001 , ierror,err,iw stop c 1001 format (1h1,20x,25hsubroutine cblktr example/// 1 10x,46hthe output from the ncar control data 7600 was// 2 32x,10hierror = 0/ 3 18x,34hdiscretization error = 1.64572e-05/ 4 12x,33hrequired length of w array = 1123// 5 10x,32hthe output from your computer is// 6 32x,8hierror =,i2/18x,22hdiscretization error =,e12.5/ 7 12x,28hrequired length of w array =,i5) c end c