c **** scd22eps.for finds 2-hump sol; \Omega = [0,1]x[0,eps] *******
c ****
c **** Investigating behaviour of zero set of solution as \epsilon varies
c **** from 1.  That is, \Omega is no longer a Square.
c **** Iniital versions of this code started with "radial" initial guess,
c **** that is the internal zero cet is roughly circular.  Iterates
c **** "fall" to a sign-changing sol with zero set intersecting boundary
c **** BUT NO LONGER intersecting at the corners, as in the Square case
c **** (see scd22.for).
C **** THIS CONFIRMS numerically our conjecture that the zero set ALWAYS
c **** intersects the boundary for solutions satisfying the sign-changing
c **** variational characterization from [CCN].

c ****
C **** John M. Neuberger
C ****    * Department of Mathematics and Statistics
C ****    * Mississippi State University
C ****    * Msstate, MS 39762
C ****    
C ****    * email:  neuberge@math.msstate.edu
C ****    * web:    http://www2.msstate.edu/~neuberge
C ****    * FAX:    (601) 325-0005
C ****    * office: (601) 325-7146
C ****    * home:   (601) 338-0007
c ****
c **** SOLVES: $-\Delta u = f(u)$ on $\Omega$, zero-D conditions.
C ****    * Finds Sign Changing Solutions by Sobolev Gradient Descent
C ****      and Projections on to Submanifolds of H=H_0^{1,2}.
c ****                                                                      
C **** Based on alogithm by John M. Neuberger
C ****    * Analysis: "A sign-changing solution to a superlinear 
C ****      Dirichlet problem", with A. Castro and J. Cossio.
C ****      To appear: Rocky Mnt. J. Math., 1995. ([CCN])
C ****    * Numerical Algorithm: "A numerical method for finding 
C ****      sign-changing solutions of superlinear Dirichlet problems".
C ****      To appear: Nonlinear World, 1996. ([NUM]}
C ****    * Look for a related sequence of recent papers by Neuberger,
C ****      Castro, and/or  Cossio.
c ****
c **** LOOP: Project pos/neg parts onto $S$, 
c ****       Sobolev descent on assembled $S_1$ member,
c ****       convergence yields $w_3$ from [CCN] and [NUM].                      
c ****
c **** More current than Jan 1, 1997.
C ****
C **** Can be easily modified to:  
C ****    * Use different nonlinearites (initially f(u) = u^3 + lambda * u)
C ****      For more, see the 3 function routines f, fp, and bf and the
C ****      assignment of a value to lambda in the initial code fragment.
C ****    * Find one sign solutions...or other solutions which should exist 
C ****      when symmetry is present (we conjecture much...).
C ****    * Use more divisions (more accuracy)
C ****      and tougher convergence criteria.
C ****    * Use different initial conditions.  See the function routine g.
C ****    * Run on different regions..I have versions for disk and annulus.
C ****    * Can find semipositone and sublinear solutions with slight to
C ****      somewhat drastic modifications.  In general, this involves 
C ****      different combinations of gradient ascent/descent and projections
C ****      onto S and S_1.
C ****
c**********************************************************************

      integer n, m, num /0/, max /5000/, go /1/ 
      parameter (m = 30)
      real dx, dy, tol /.00025/
      real rstep /.25/, sstep /0.1/
      real bju, de1, de2, gnorm
      real u, u1, u2, grdj, v
      integer i, j    
      
      real lambda, eps, p1, p2
      common lambda, eps, p1, p2              
      
      dimension u(0:m,0:m), u1(0:m,0:m), u2(0:m,0:m)
      dimension grdj(0:m,0:m),v(0:m,0:m)
      
      n        = m
      eps      = 0.9
      dx       = 1.0 / n
      dy       = eps / n
      lambda   = 1.0
      p1       = 3.0
      p2       = 3.0      
      
      open(unit=1,file='d22.txt')

c **** initialize u
      do i = 1,n-1                                 
         do j = 1,n-1
            u(i,j) = g(float(i)/n,eps*float(j)/n)
         end do        
      end do

      call posneg(u(0,0),u1(0,0),u2(0,0),n)
      call sobray(u1(0,0),n,rstep, dx, dy)
      call sobray(u2(0,0),n,rstep, dx, dy)

      do i=1,n-1
         do j=1,n-1
            u(i,j) = u1(i,j) + u2(i,j)
         end do
      end do


c *** MAIN LOOP ***********************************************
      do while (go .eq. 1)
         num = num + 1 
         
         call gj(u(0,0), grdj(0,0),gnorm, n,dx, dy, v(0,0))
         call sobstep(u(0,0),grdj(0,0), sstep, convg, sup, n)
         call posneg(u(0,0),u1(0,0),u2(0,0),n)
         call sobray(u1(0,0),n,rstep, dx, dy)
         call sobray(u2(0,0),n,rstep, dx, dy)
         do i=1,n-1
            do j=1,n-1
               u(i,j) = u1(i,j) + u2(i,j)
            end do
         end do

c *****  DIAGNOSTICS and OUTPUT
                      
         call ju(u(0,0), bju, n, dx, dy)
         call test(u(0,0),de1, de2, n,dx, dy)
                      
         write(*,*) " num, sup, convg: ", num,sup,convg
         write(*,*) " lambda, n, eps ", lambda, n, eps
         write(*,*) "J(u)     = ",bju
         write(*,*) "max/avg  = ",de1, de2
         write(*,*) "gnorm    = ",gnorm

         WRITE (*, 110) ((u(i,j), i = n/10, 9*n/10,n/10),
     +   j = n/10, 9*n/10,n/10)
         
         write(*,*) ' '
         WRITE (*, 110) ((grdj(i,j), i = n/10, 9*n/10,n/10),
     +   j = n/10, 9*n/10,n/10)

         write (1,*) ' num, sup, convg: ', num,sup,convg
         write (1,200)
         write (1,*) "J(u) = ",bju          
         write (1,200)
         write (1,*) " lambda, n ", lambda, n
         write (1,200)
         write (1,*) "gnorm    = ",gnorm
         write (1,200)
         write (1,*) "max/avg = ",de1, de2
         write (1,200)

         WRITE (1, 110) ((u(i,j), i = n/10, 9*n/10,n/10),
     +   j = n/10, 9*n/10,n/10)
         write (1,200)
         WRITE (1, 110) ((grdj(i,j), i = n/10, 9*n/10,n/10),
     +   j = n/10, 9*n/10,n/10)
                        
         if ((num .gt. max) .or. (convg .lt. tol)) then
            go = 0                              
            close(unit=1)
         end if
      end do                            
      pause 'hit any key'

110   format (9(f8.5))
200   format (' ',/)     

      end                       
c**** END MAIN program *****************************************      
      subroutine posneg(u,u1,u2,n)
            
      real u,u1,u2
      dimension u(0:n,0:n),u1(0:n,0:n),u2(0:n,0:n)
      integer i, j, n
      
      do i = 1,n-1   
         do j = 1,n-1
            u1(i,j) = 0
            u2(i,j) = 0
            if (u(i,j) .gt. 0) then
                  u1(i,j) = u(i,j)
            else
                  u2(i,j) = u(i,j)
            end if
         end do               
      end do
           
      end                                       
c**************************  
c**** nonlinear term function (and derivative/antiderivative)
      real function f(x)
      real x
      real lambda, eps, p1, p2
      common lambda, eps, p1, p2              
      
      if (x .gt. 0) then
         f =      x**p1 + lambda*x
      else
         f =  -(-x)**p2 + lambda*x
      end if
            
      end
      real function fp(x)
      real x
      real lambda, eps, p1, p2
      common lambda, eps, p1, p2              
      
      if (x .gt. 0) then
         fp =     p1*x**(p1-1) + lambda
      else
         fp =  p2*(-x)**(p2-1) + lambda
      end if

      end
      real function bf(x)
      real x
      real lambda, eps, p1, p2
      common lambda, eps, p1, p2              

      if (x .gt. 0) then
         bf =     (1/(p1+1))*x**(p1+1) + .5*lambda*x**2
      else
         bf =  (1/(p2+2))*(-x)**(p2+1) + .5*lambda*x**2
      end if

      end
c**************************
c**** initial guess function 
      real function g(x,y)        
      real x,y, pi /3.14159/, theta, sum
      real lambda, eps, p1, p2
      common lambda, eps, p1, p2              
      
      theta = eps*pi/4.0

c      sum =       (sin(pi*x)*sin(2*pi*y/eps)) * cos(theta)
c      sum = sum + (sin(2*pi*x)*sin(pi*y/eps)) * sin(theta)
c      g   = sum
      
c **** KEEP BLOCK, goes to and falls off of "radial" solution
      d = ((x-.5)**2 + (y-.5*eps)**2) ** .5
      g = sin(pi*x)*sin(pi*y/eps)*cos(3*pi*d/2)
      
      end
c************************** 
c**** Reanka's SOR: w = 2/(1+sqr(1-m**2)), m = spect rad Jac It mat
c****               here m = 1-2*sin(pi/(2*n))**2
      subroutine ilap(b,u,n,dx, dy)
      
      real w /1.73/, total, convg, u, b, max, tol /.000001/, dx, dy
      dimension u(0:n,0:n),b(0:n,0:n)
      integer i,j,go, n
      
      go = 1
      do while (go .eq. 1)
         max = 0
         do i = 1,n-1
            do j = 1,n-1
               total = (dx**2)*(u(i-1,j)+u(i+1,j))
               total = total + (dy**2)*(u(i,j-1)+u(i,j+1))
               total = ( total-b(i,j)*(dx*dy)**2 ) / (2*(dx**2+dy**2))
               total = total - u(i,j)
               u(i,j) = u(i,j) + w*total
               convg = abs(total) / (1+abs(u(i,j)))
               if (convg .gt. max) then 
                  max = convg
               end if
            end do
         end do
         if (max .lt. tol) then
            go = 0
         end if
      end do
      end
c**************************
      subroutine gj(u,grdj,gnorm, n,dx, dy,v)
      
      real u, v, grdj, dx, dy, gnorm
      dimension u(0:n,0:n),v(0:n,0:n),grdj(0:n,0:n)
      integer i,j,n                           
      
      do i = 1,n-1
         do j = 1,n-1
            v(i,j)    =  f(u(i,j))
         end do
      end do                    
      
      call ilap(v(0,0),grdj(0,0),n,dx, dy)
      
      gnorm = 0
      do i = 1,n-1
         do j = 1,n-1
            grdj(i,j) = grdj(i,j) + u(i,j)
            gnorm = gnorm + grdj(i,j)**2 
         end do
      end do                            
      gnorm = (gnorm * dx * dy)**.5
      
      end
c**************************            
      subroutine sobray(u,n,rstep, dx, dy)
      
      real u, rstep, tol /.00005/, gamu, uu, dx, dy
      dimension u(0:n,0:n)
      integer i,j,n,max /200/, go,num

      go = 1
      num = 0
      do while (go .eq. 1)
         num = num + 1
         
         gamu = gam(u(0,0), n, dx, dy) 
         uu   = dot(u(0,0),u(0,0),n, dx, dy)

         do i = 1,n-1
            do j = 1,n-1
               u(i,j) = u(i,j) + rstep * gamu * u(i,j) / uu
            end do
         end do
         
         if ((num .gt. max) .or. (abs(gamu) .lt. tol)) then
            go = 0
         end if
      end do                    

      end
c**************************
      subroutine sobstep(u,grdj, sstep, convg, sup, n)
      
      real u, grdj, sstep, convg, uold, sup
      dimension u(0:n,0:n), grdj(0:n,0:n)
      integer i,j,n

      convg = 0
      sup   = 0                                   

      do i = 1,n-1
         do j = 1,n-1
            uold = u(i,j)
            u(i,j) = u(i,j) - sstep * grdj(i,j)
            if (abs(u(i,j)) .gt. sup) then
               sup = abs(u(i,j))
            end if
            convg = convg +  (u(i,j) - uold) ** 2
         end do
      end do                                 
      convg = (convg / (n**2)) ** .5
      
      end
c ********************************************************************
c ******** H = H_0^{1,2} inner product routine ***********************
c COMPUTE $rangle u, v \langle$ = $\int \nabla u \cdot \nabla v dx$
      real function dot(u,v,n, dx, dy)
            
      real u,v,sum, term, dx, dy
      dimension u(0:n,0:n), v(0:n,0:n)
      integer i,j,n
           
      sum = 0
      do i=1,n
         do j =1,n
                        
c           average adjacent grid corners to do differencing at centers
            
c           top center            
            u1 = (u(i-1,j)+u(i-1,j-1))/2.
            u2 = (u(i,j)+u(i,j-1))/2.
            
c           bottom center
            u3 = (u(i-1,j-1)+u(i,j-1))/2.
            u4 = (u(i-1,j)+u(i,j))/2.
            
c           left center
            v1 = (v(i-1,j)+v(i-1,j-1))/2.
            v2 = (v(i,j)+v(i,j-1))/2.

c           right center
            v3 = (v(i-1,j-1)+v(i,j-1))/2.
            v4 = (v(i-1,j)+v(i,j))/2.

c           (u_x v_x + u_y v_y) dx dy
            term = (u2-u1)*(v2-v1)/(dx**2) + (u4-u3)*(v4-v3)/(dy**2)
c            term = (u2-u1)*(v2-v1) + (u4-u3)*(v4-v3)
         
            sum = sum + term
            
         end do  
      end do  
      
      dot = sum * dx * dy
c      dot = sum
      end
c ********************************************************************
c**** compute gamma = $||u||^2 - \int u f(u) dx$
      real function gam(u,n, dx, dy)
     
      real u,sum, dx, dy
      dimension u(0:n,0:n)
      integer i,j,n
      
      sum = 0
      do i=1,n
         do j=1,n                
            
c           average 4 grid corners to do calculation at centers
            ui    =  (u(i,j)+u(i-1,j)+u(i-1,j-1)+u(i,j-1))/4
            sum   =  sum - ui * f(ui)

         end do
      end do  
      gam = (sum * dx * dy) + dot(u(0,0),u(0,0),n, dx, dy)
      
      end                                        
c ********************************************************************
c**** compute J(u), dx2/dx2 = 1 in term
c**** $J(u) = ||u||^2/2 - \int F(u) dx$
      subroutine ju(u,bju,n,dx, dy)
     
      real u, dx, dy, bju, ui
      integer i,j,n
      dimension u(0:n,0:n)
      
      bju = (dot(u(0,0),u(0,0),n, dx, dy) / 2)
      do i=1,n
         do j=1,n

c           average 4 grid corners to do calculation at centers
            ui    =  (u(i,j)+u(i-1,j)+u(i-1,j-1)+u(i,j-1))/4
            bju   =  bju - bf(ui)*dx*dy

          end do
      end do       

      end      
c ********************************************************************
c****************************************
c**** compute Convergence Criteria
c *** max = sup norm of PDE(u), avg = $L^2$ norm of PDE(u)
      subroutine test(u,max, avg, n,dx, dy)
     
      real u, term, dx, dy, max , avg
      integer i,j,n
      dimension u(0:n,0:n)        
      avg =    0
      max =    0
            
      sum = 0
      do i=1,n-1
         do j=1,n-1                                       
         
c           SECOND DIFFerencing, use grid points.  Improvable?
            term = (u(i+1,j) + u(i-1,j) - 2*u(i,j) ) / dx**2
            term = term + ( u(i,j+1) + u(i,j-1) - 2*u(i,j) ) / dy**2
            term = (term + f(u(i,j)))**2
            avg = avg + term
            
            if (term .gt. max) then
               max = term
            end if
         end do
      end do                           
      max = max ** .5
      avg = (avg * dx*dy)**.5

      end
c ********** END OF FILE **********************************************