PROGRAM WRM4X4S1DN
C=======================================================================
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)
C              AND BOUNDARY CONDITIONS OF U(0)=U0. & DU/DX(L)=S;
C
C                      DEGREE OF FREEDUM = 4
C
C -------------- VARIABLE DEFNITION ------- OCT/14/2024  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:           [B] * {A} + {C} = {0.}
C--------------------      F0 = U0 + S*RL*(X/RL)
C--------------------      F1 = (X/RL)*(1.D0-X/RL) + (X/RL)
C--------------------      F2 = F1**2
C--------------------      F3 = F1**3
C--------------------      F4 = F1**4
C=======================================================================
      IMPLICIT REAL * 8 ( A-H , O-Z )
      PARAMETER ( N = 6, NSEG=10000, MULTI=10, MXN=4 )
      DIMENSION SAI(N) , W(N), A(MXN,MXN), C(MXN)
      COMMON / DEL / DELTAX     /DOMAIN/ RL      /BORDER/ U0, S
      COMMON / COFF / ALPHASQ
      EXTERNAL F0, F1, F2, F3, F4
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
C=======================================================================
      CALL GRULE ( N , SAI , W )
C=======================================================================
C------------------- MATERIAL DATA AND BOUNDARY VALUES
      ALPHASQ=1.D0
      XST=0.D0
      XEN=0.5D0
      U0 = 1.D0
      S  = 0.D0
C=======================================================================
      OPEN ( 1, FILE='WRM4X4S1DN.FEM',STATUS='UNKNOWN' )
      WRITE(1,*)'==== ONE DIMENSIONAL HELMHOLTZ EQUATION DOF=4 ===='
      WRITE(1,*)'==== DIRICHLET ------- NEUMANN PROBLEM ===='
      WRITE(1,*)'---- GALERKINS WEIGHTING FUNCTION----'
      WRITE(1,*)'# OF GL INTEGRATION SAMPLING POINTS =',N
      WRITE(1,*)'APPROXIMATING FUNCTION: U(X) = F0(X) + A1*F1(X)......'
      WRITE(1,*)'WHERE F0(X) = U0+SL(X/L)'
      WRITE(1,*)'F1(X) = (X/L)*(1-X/L) + (X/L)'
      WRITE(1,*)'F2(X) = ((X/L)**K*(1-X/L))**2'
      WRITE(1,*)'F3(X) = ((X/L)**K*(1-X/L))**3'
      WRITE(1,*)'F4(X) = ((X/L)**K*(1-X/L))**4'
C=======================================================================
C   DELTAX: SPACIAL DEFERENTIAL LENGTH FOR DERIVATIVE EVALUATION.
      RL = XEN - XST
      DELTAX = RL / ( MULTI * NSEG )
C=======================================================================
      WRITE(1,*)'LENGTH OF DOMAIN =',RL
      WRITE(1,*)'NUMBER OF SEGMENTS FOR INTEGRATION =', NSEG
      WRITE(1,*)'DX FOR DERIVATIVE EVALUATION =', DELTAX
C=======================================================================
      WRITE(1,*)'X-COORDINATE OF LEFT  END BOUNDARY =',XST
      WRITE(1,*)'X-COORDINATE OF RIGHT END BOUNDARY =',XEN
      WRITE(1,*)'ALPHASQ =', ALPHASQ
      WRITE(1,*)'NUMBER OF SEGMENTS =', NSEG
      WRITE(1,*)'DX FOR DERIVATIVE EVALUATION =', DELTAX
C----------------BOUNDARY CONDITIONS
      WRITE(1,*) 'U(X) AT X=0 =',U0
      WRITE(1,*) 'S    AT X=L =', S
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, F1, C(1) )
      CALL INTE ( ALPHASQ, XST, XEN, NSEG, N, SAI, W, F0, F2, C(2) )
      CALL INTE ( ALPHASQ, XST, XEN, NSEG, N, SAI, W, F0, F3, C(3) )
      CALL INTE ( ALPHASQ, XST, XEN, NSEG, N, SAI, W, F0, F4, C(4) )
C
      CALL INTE ( ALPHASQ, XST, XEN, NSEG, N, SAI, W, F1, F1, A(1,1) )
      CALL INTE ( ALPHASQ, XST, XEN, NSEG, N, SAI, W, F1, F2, A(1,2) )
      CALL INTE ( ALPHASQ, XST, XEN, NSEG, N, SAI, W, F1, F3, A(1,3) )
      CALL INTE ( ALPHASQ, XST, XEN, NSEG, N, SAI, W, F1, F4, A(1,4) )
C
      A(2,1) = A(1,2)
      CALL INTE ( ALPHASQ, XST, XEN, NSEG, N, SAI, W, F2, F2, A(2,2) )
      CALL INTE ( ALPHASQ, XST, XEN, NSEG, N, SAI, W, F2, F3, A(2,3) )
      CALL INTE ( ALPHASQ, XST, XEN, NSEG, N, SAI, W, F2, F4, A(2,4) )
C
      A(3,1) = A(1,3)
      A(3,2) = A(2,3)
      CALL INTE ( ALPHASQ, XST, XEN, NSEG, N, SAI, W, F3, F3, A(3,3) )
      CALL INTE ( ALPHASQ, XST, XEN, NSEG, N, SAI, W, F3, F4, A(3,4) )
C
      A(4,1) = A(1,4)
      A(4,2) = A(2,4)
      A(4,3) = A(3,4)
      CALL INTE ( ALPHASQ, XST, XEN, NSEG, N, SAI, W, F4, F4, A(4,4) )
C=======================================================================
C  EVALUATION OF UNKNOWN A1 AND A2 IN THE APPROXIMATING FUNCTION U(X)
      DO I = 1 , MXN
      C(I) = -C(I)
      END DO
      CALL SYSTEM1 ( MXN, MXN , A , C )
C=======================================================================
C                     PRINTING RESULTS
      WRITE (1,*) 'A1=',C(1)
      WRITE (1,*) 'A2=',C(2)
      WRITE (1,*) 'A3=',C(3)
      WRITE (1,*) 'A4=',C(4)
      WRITE (1,*) 'U(X)=F0(X)+',C(1),'*F1(X)+',C(2),'*F2(X)+',
     *                          C(3),'*F3(X)+',C(4),'*F4(X)'
C=======================================================================
      CALL OUTPUT ( XST, XEN, MXN, C )
      CLOSE (1)
      STOP 'NORMAL TERMINATION'
      END
C
C
      SUBROUTINE INTE ( ALPHASQ,XST,XEN,NSEG,INTEPT,SAI,W, G1,G2,TOTAL )
      IMPLICIT REAL * 8 ( A-H , O-Z )
      DIMENSION SAI(INTEPT) , W(INTEPT)
      EXTERNAL G1, G2
      TOTAL = 0.D0
      DX = ( XEN - XST ) / NSEG
      DO I = 1 , NSEG
      X1 = DX*(I-1) + XST
      X2 = X1 + DX
      SUM = 0.D0
      SH  = ( X2 - X1 ) / 2.D0
      AVE = ( X1 + X2 ) / 2.D0
      DO J = 1 , INTEPT
      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
      COMMON /BORDER/ U0, S
      F0 = U0 + S*RL*(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) + (X/RL)
      RETURN
      END
C
C
      FUNCTION F2(X)
      IMPLICIT REAL * 8 ( A-H , O-Z )
      F2 = F1(X)**2
      RETURN
      END
C
C
      FUNCTION F3(X)
      IMPLICIT REAL * 8 ( A-H , O-Z )
      F3 = F1(X)**3
      RETURN
      END
C
C
      FUNCTION F4(X)
      IMPLICIT REAL * 8 ( A-H , O-Z )
      F4 = F1(X)**4
      RETURN
      END
C
C
      FUNCTION DERIV(F,X)
      IMPLICIT REAL * 8 ( A-H , O-Z )
      COMMON / DEL / DELTAX
      EXTERNAL F
      FIP2 = F(X+2.D0*DELTAX)
      FIP1 = F(X+     DELTAX)
      FIN1 = F(X-     DELTAX)
      FIN2 = F(X-2.D0*DELTAX)
      DERIV = (-FIP2+8.D0*FIP1-8.D0*FIN1+FIN2 )/(12.D0*DELTAX)
C      DERIV = ( F(X+DELTAX) - F(X-DELTAX) ) / ( 2.D0*DELTAX )
      RETURN
      END
C
C
      SUBROUTINE OUTPUT ( XST,XEN, MXN, C )
      IMPLICIT REAL * 8 ( A-H , O-Z )
      DIMENSION C(MXN)
      EXTERNAL F0, F1, F2, F3, F4
      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)+C(1)*F1(X)+C(2)*F2(X)+C(3)*F3(X)+C(4)*F4(X)
      DUDX = DERIV(F0,X)+C(1)*DERIV(F1,X)+C(2)*DERIV(F2,X)+
     *       C(3)*DERIV(F3,X)+C(4)*DERIV(F4,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, S
      COMMON / COFF / ALPHASQ
      A = U0
      AL = DSQRT (ALPHASQ)
      B = (AL*U0*DSIN(AL*RL)+S)/(AL*DCOS(AL*RL))
      EXACT = A*DCOS(AL*X) + B*DSIN(AL*X)
      RETURN
      END
C
C
      SUBROUTINE SYSTEM1 ( MXN, N , A , C )
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION A (MXN,MXN) , C (MXN)
      N1 = N - 1
      DO K = 1, N1
      L = K + 1
      DO I = L , N
      P = A (I,K) / A (K,K)
      IF ( P .NE. 0. ) THEN
      DO J = L , N
      A (I,J) = A (I,J) - P * A ( K , J )
      END DO
      C ( I ) = C ( I) - P * C ( K )
      END IF
      END DO
      END DO
C---- BACK SUBSTITUTION
      C (N) = C (N) / A (N,N)
      DO K = 1 , N1
      I = N - K
      L = I + 1
      P = C ( I )
      DO J = L , N
      P = P - C (J) * A (I,J)
      END DO
      C ( I ) = P / A (I,I)
      END DO
      RETURN
      END
C
C
      SUBROUTINE GRULE ( N , SAI , W )
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION SAI(N) , W(N)
      IF ( N .LT. 2 ) STOP'N<2'
      IF ( N .GT. 6 ) STOP'N>6'
      IF ( N .EQ. 2 ) THEN
      SAI(1) = DSQRT(3.D0)/3.D0
      W(1) = 1.D0
      SAI(2) = - SAI(1)
      W(2) = W(1)
      RETURN
      END IF
      IF ( N .EQ. 3 ) THEN
      SAI(1) = DSQRT(15.D0)/5.D0
      SAI(2) = 0.D0
      W(1) = 5.D0/ 9.D0
      W(2) = 8.D0/ 9.D0
      SAI(3) = - SAI(1)
      W(3) = W(1)
      RETURN
      END IF
      IF ( N .EQ. 4 ) THEN
      SAI(1) = 0.33998104358485D0
      SAI(2) = 0.86113631159405D0
        W(1) = 0.65214515486254D0
        W(2) = 0.34785484513745D0
      SAI(3) = -SAI(2)
      SAI(4) = -SAI(1)
      W(3) = W(2)
      W(4) = W(1)
      RETURN
      END IF
      IF ( N .EQ. 5 ) THEN
      SAI(1) = 0.90617984593866D0
      SAI(2) = 0.53846931010568D0
      SAI(3) = 0.D0
        W(1) = 0.23692688505619D0
        W(2) = 0.47862867049937D0
        W(3) = 5.12D0 / 9.D0
      SAI(4) = -SAI(2)
      SAI(5) = -SAI(1)
      W(4) = W(2)
      W(5) = W(1)
      RETURN
      END IF
      IF ( N .EQ. 6 ) THEN
      SAI(1) = 0.23861918608320D0
      SAI(2) = 0.66120938646626D0
      SAI(3) = 0.93246951420315D0
        W(1) = 0.46791393457269D0
        W(2) = 0.36076157304814D0
        W(3) = 0.17132449237917D0
      SAI(4) = -SAI(1)
      SAI(5) = -SAI(2)
      SAI(6) = -SAI(3)
        W(4) = W(1)
        W(5) = W(2)
        W(6) = W(3)
      END IF
      RETURN
      END