PROGRAM ELEVEN
C=======================================================================
C SOLUTION OF D2U/DXDX + ALPHA*U = 0 USING WEIGHTED RESIDUAL METHOD
C         WITH AN APPROXIMATING FUNCTION OF U(X)=F0(X)+A1*F1(X)
C              AND BOUNDARY CONDITIONS OF U(0)=1. & DU/DX(L)=S;
C -------------- VARIABLE DEFNITION ----------- 12/2/2004  EIJI FUKUMORI 
C   XST & XEN: INTEGRATION LIMITS; NSEG: NUMBER OF SEGMENTS IN LIMITS;
C   UNKNOWN COEFFICENT (A1) IN THE APPROXIMATING FUNCTION IS EVALUATED
C    BY THE FOLLOWING EQUATION:           B1 * A1 + C1 = -S
C=======================================================================
      IMPLICIT REAL * 8 ( A-H , O-Z )
      PARAMETER ( N = 3, ALPHA=1., XST=0., XEN=0.5, NSEG=100, MULTI=10 )
      DIMENSION SAI(N) , W(N)
      COMMON / DEL / DELTAX /DOMAIN / RL, S, U0
      EXTERNAL F0, F1
C=======================================================================
C            THREE-SAMPLING-POINT GAUSS INTEGRATION METHOD
C            N: NUMBER OF SAMPLING POINTS IN EACH SEGMENET
C   SAI(I) & W(I): NON-DIMENSIONALIZED COORDINATE & WEIGHTING FACTOR
      DATA SAI/-0.7745966692415D0,0.0000000000000D0, 0.7745966692415D0/
      DATA W  / 0.5555555555555D0,0.8888888888888D0, 0.5555555555555D0/
C=======================================================================
      RL = XEN - XST
      U0 = 1.D0
      WRITE (*,240)
  240 FORMAT( 'Type in S = '  $ )
      READ(*,*) S
C=======================================================================
      OPEN ( 1, FILE='ELEVEN.FEM',STATUS='UNKNOWN' )
      WRITE(1,*)' APPROXIMATING FUNCTION: F0(X) + A1*F1(X)'
C=======================================================================
C   DELTAX: SPACIAL DEFERENTIAL LENGTH FOR DERIVATIVE EVALUATION.
      DELTAX = ( XEN - XST ) / ( MULTI * NSEG )
C=======================================================================
C                COMPUTATION OF H(F0,F1) AND H(F1,F1)
      CALL INTE ( ALPHA, XST, XEN, NSEG, N, SAI, W, F0, F1, C1 )
      CALL INTE ( ALPHA, XST, XEN, NSEG, N, SAI, W, F1, F1, B1 )
C=======================================================================
C      EVALUATION OF UNKNOWN A1 IN THE APPROXIMATING FUNCTION U(X)
      A1 = (-S - C1) / B1
C=======================================================================
C                     PRINTING RESULTS
      WRITE (1,100) B1, C1, -S
      WRITE (1,110) A1
  100 FORMAT ( 1X, F20.10, 1X, '* A1 +', F20.10, ' =', F20.10 )
  110 FORMAT ( 2X, 'U(X) = F0(X) + ',F15.10, ' * F1(X)' )
      X = 0.5
      UOFX = F0(X) + A1*F1(X)
      WRITE (1,*) ' X=',X, '      U(X)=',UOFX
      CALL INTE1 ( ALPHA,XST, XEN, NSEG, N, SAI, W, A1,SQERROR )
      WRITE (1,*) ' SQUARE ERROR =', SQERROR
      CLOSE (1)
      STOP
      END
C
C
      SUBROUTINE INTE ( ALPHA,XST,XEN,NSEG, N,SAI,W, G1,G2, TOTAL )
      IMPLICIT REAL * 8 ( A-H , O-Z )
      DIMENSION SAI(N) , W(N)
      EXTERNAL G1, G2
      TOTAL = 0.
      DX = ( XEN - XST ) / NSEG
      DO I = 1 , NSEG
      X1 = DX*(I-1)
      X2 = X1 + DX
      SUM = 0.
      SH  = ( X2 - X1 ) / 2.
      AVE = ( X1 + X2 ) / 2.
      DO J = 1 , N
      X = SH * SAI(J) + AVE
      SUM = SUM + (-DERIV(G1,X)*DERIV(G2,X)+ALPHA*G1(X)*G2(X)) * W(J)
      END DO
      TOTAL = TOTAL + SH * SUM
      END DO
      RETURN
      END
C
C
      FUNCTION F0(X)
      IMPLICIT REAL * 8 ( A-H , O-Z )
      COMMON / DOMAIN / RL, S, U0
      F0 = U0 - S*RL* ( 1 - X/RL )* ( X/RL )
      RETURN
      END
C
C
      FUNCTION F1(X)
      IMPLICIT REAL * 8 ( A-H , O-Z )
      COMMON / DOMAIN / RL, S, U0
      F1 = (1-X/RL )* (X/RL ) + (X/RL )
      RETURN
      END
C
C
      FUNCTION DERIV(F,X)
      IMPLICIT REAL * 8 ( A-H , O-Z )
      COMMON / DEL / DELTAX
      EXTERNAL F
      DERIV = ( F(X+DELTAX) - F(X-DELTAX) ) / ( 2.*DELTAX )
      RETURN
      END
C
C
      SUBROUTINE INTE1 ( ALPHA,XST,XEN,NSEG, N,SAI,W,A1,TOTAL )
      IMPLICIT REAL * 8 ( A-H , O-Z )
      COMMON  /DOMAIN / RL, S, U0
      DIMENSION SAI(N) , W(N)
      TOTAL = 0.
      DX = ( XEN - XST ) / NSEG
      C=DSQRT(ALPHA)
      A = ( S+U0*C*DSIN(C*RL) ) / (C*DCOS(C*RL))
      B = U0
      DO I = 1 , NSEG
      X1 = DX*(I-1)
      X2 = X1 + DX
      SUM = 0.
      SH  = ( X2 - X1 ) / 2.
      AVE = ( X1 + X2 ) / 2.
      DO J = 1 , N
      X = SH * SAI(J) + AVE
      EXACT = A*DSIN(C*X) + B*DCOS(C*X)
      APPRO = F0(X) + A1*F1(X)
      SUM = SUM + W(J)*(APPRO-EXACT)**2
      END DO
      TOTAL = TOTAL + SH * SUM
      END DO
      X = 0.
      N = 10
      DX =  ( XEN - XST ) / N
      WRITE(1,*) ' X  APPROXIMATION EXACT'
      WRITE(1,*) X, U0, U0
      DO I = 1 , N
      X = X + DX
      EXACT = A*DSIN(C*X) + B*DCOS(C*X)
      APPRO = F0(X) + A1*F1(X)
      WRITE (1,*) X, APPRO , EXACT
      END DO
      RETURN
      END