PROGRAM WRM2X2H
C=======================================================================
C     SHAPE FUNCTIONA ARE NOT EQUAL TO THE WEIGHTING FUNCTIONS.
C SOLUTION OF D2U/DXDX + ALPHASQ*U = 0 USING WEIGHTED RESIDUAL METHOD
C   WITH AN APPROXIMATING FUNCTION OF U(X)=F0(X)+A1*F1(X)+A2*F2(X)
C              AND BOUNDARY CONDITIONS OF U(0)=0. & U(1)=1.
C      ******* NON-GALERKINS WEIGHTING FUNCTIONS H1 AND H2 ********
C -------------- VARIABLE DEFNITION ----------- 6/1/1993  EIJI FUKUMORI 
C XST & XEN = INTEGRATION LIMITS. NSEG = NUMBER OF SEGMENTS IN LIMITS.
C UNKNOWN COEFFICENT (A1&A2) IN THE APPROXIMATING FUNCTION IS EVALUATED
C BY THE FOLLOWING SIMULTANEOUSEQUATION: B1 * A1 + B2 * A2 + C1 = 0.
C                                        B3 * A1 + B4 * A2 + C2 = 0. 
C=======================================================================
      IMPLICIT REAL * 8 ( A-H , O-Z )
      PARAMETER ( N = 3, NSEG=100, MULTI=10 )
      DIMENSION SAI(N) , W(N)
      COMMON / DEL /DELTAX /DOMAIN /RL  /BORDER/U0,UL
      COMMON /CONST/ ALPHASQ
      EXTERNAL F0, F1, F2, H1, H2
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=======================================================================
C    MATERIAL DATA AND BOUNDARY VALUES
      ALPHASQ=1.D0
      XST=0.D0
      XEN=1.D0
      U0 = 0.D0
      UL = 1.D0
C=======================================================================
      OPEN ( 1, FILE='WRM2X2NONGALERKIN.FEM',STATUS='UNKNOWN' )
      WRITE(1,*)' ==== DIRICHLET - DIRICHLET PROBLEM ===='
      WRITE(1,*)'  ---- NON-GALERKINS WEIGHTING FUNCTION H1 AND H2----'
      WRITE(1,*)' APPROXIMATING FUNCTION: F0(X) + A1*F1(X) + A2*F2(X)'
      WRITE(1,*)' WHERE F0(X)=X, F1(X)=X*(1-X), & F2(X)=X*X*(1-X)'
C=======================================================================
C   DELTAX = SPACIAL DEFERENTIAL LENGTH FOR DERIVATIVE EVALUATION.
      RL = XEN - XST
      DELTAX = RL / ( MULTI * NSEG )
C=======================================================================
COMPUTATION OF H(F0,F1) H(F1,F1) H(F1,F1) H(F1,F2) H(F2,F1) H(F2,F2)
      CALL INTE ( ALPHASQ, XST, XEN, NSEG, N, SAI, W, F0, H1, C1 )
      CALL INTE ( ALPHASQ, XST, XEN, NSEG, N, SAI, W, F0, H2, C2 )
      CALL INTE ( ALPHASQ, XST, XEN, NSEG, N, SAI, W, F1, H1, B1 )
      CALL INTE ( ALPHASQ, XST, XEN, NSEG, N, SAI, W, F1, H2, B3 )
      CALL INTE ( ALPHASQ, XST, XEN, NSEG, N, SAI, W, F2, H1, B2 )
      CALL INTE ( ALPHASQ, XST, XEN, NSEG, N, SAI, W, F2, H2, B4 )
C=======================================================================
C  EVALUATION OF UNKNOWN A1 AND A2 IN THE APPROXIMATING FUNCTION U(X)
      A1 = - ( C1*B4 - B2*C2 ) / ( B1*B4 - B2*B3 )
      A2 = - ( B1*C2 - C1*B3 ) / ( B1*B4 - B2*B3 )
C=======================================================================
C                     PRINTING RESULTS
      WRITE (1,100) B1, B2, C1, B3, B4, C2
      WRITE (1,110) A1, A2
  100 FORMAT( 1X,F15.10,1X,'* A1 +',F15.10,'* A2 +',F15.10,' = 0' )
  110 FORMAT(2X,'U(X) = F0(X) +',F15.10,' * F1(X) +',F15.10,'*F2(X)')
      CALL OUTPUT ( XST, XEN, A1, A2 )
      CLOSE (1)
      STOP
      END
C
C
      SUBROUTINE INTE ( ALPHASQ,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.D0
      AVE = ( X1 + X2 ) / 2.D0
      DO J = 1 , N
      X = SH * SAI(J) + AVE
      SUM = SUM + (-DERIV(G1,X)*DERIV(G2,X)+ALPHASQ*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  /BORDER/U0,UL
      F0 = U0*(1.D0-X/RL) + UL*(X/RL)
      RETURN
      END
C
C
      FUNCTION F1(X)
      IMPLICIT REAL * 8 ( A-H , O-Z )
      COMMON / DOMAIN / RL
      F1 = X/RL * ( 1.D0 - X/RL )
      RETURN
      END
C
C
      FUNCTION F2(X)
      IMPLICIT REAL * 8 ( A-H , O-Z )
      COMMON / DOMAIN / RL
      F2 = X/RL * X/RL * ( 1.D0- 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
      FUNCTION H1(X)
      IMPLICIT REAL * 8 ( A-H , O-Z )
      COMMON / DOMAIN / RL
      H1 = X/RL * ( 1.D0 - (X/RL)**2 )
      RETURN
      END
C
C
      FUNCTION H2(X)
      IMPLICIT REAL * 8 ( A-H , O-Z )
      COMMON / DOMAIN / RL
      H2 = X/RL * X/RL * ( 1.D0- X/RL )**2
      RETURN
      END
C
C
      SUBROUTINE OUTPUT ( XST,XEN,A1,A2 )
      IMPLICIT REAL * 8 ( A-H , O-Z )
      EXTERNAL F0, F1, F2
      NDIV = 10
      DX = ( XEN - XST ) / NDIV
      WRITE(1,*)'X-COORDINATE U(X) DU/DX EXACT(X) |U(X)-EXACT(X)|'
      DO I = 1 , NDIV+1
      X = DX*(I-1) + XST
      UX = F0(X) + A1*F1(X) + A2*F2(X)
      DUDX = DERIV(F0,X)+A1*DERIV(F1,X)+A2*DERIV(F2,X)
      WRITE(1,*) X, UX, DUDX, EXACT(X), DABS(UX-EXACT(X))
      END DO
      RETURN
      END
C
C
      FUNCTION EXACT(X)
      IMPLICIT REAL * 8 ( A-H , O-Z )
      COMMON / DEL /DELTAX /DOMAIN /RL  /BORDER/U0,UL
      COMMON /CONST/ ALPHASQ
      AL = DSQRT(ALPHASQ)
      A = U0
      B = (UL-U0*DCOS(AL*RL))/DSIN(AL*RL)
      EXACT = A*DCOS(AL*X) + B*DSIN(AL*X)
      RETURN
      END