PROGRAM WRM2X2DD
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)+A2*F2(X)
C AND BOUNDARY CONDITIONS OF U(0)=1. & U(1)=1.
C================= DIRICHLET BC ------ DIRICHLET BC ====================
C -------------- VARIABLE DEFNITION ----------- 10/13/2024 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 2X2 SIMULTANEOUSEQUATION.
C=======================================================================
IMPLICIT REAL * 8 ( A-H , O-Z )
PARAMETER ( N = 6, NSEG=1000, MULTI=10 )
DIMENSION SAI(N) , W(N)
COMMON / DEL /DELTAX /DOMAIN /RL /BORDER/U0,UL
COMMON /CONST/ ALPHASQ
EXTERNAL F0, F1, F2
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=1.D0
U0 = 1.D0
UL = 1.D0
C=======================================================================
OPEN ( 1, FILE='WRM2X2DD.FEM',STATUS='UNKNOWN' )
WRITE(1,*)'==== DIRICHLET - DIRICHLET PROBLEM ===='
WRITE(1,*)'---- GALERKIN WEIGHTING FUNCTION ----'
WRITE(1,*)'APPROXIMATING FUNCTION: F0(X) + A1*F1(X) + A2*F2(X)'
WRITE(1,*)'WHERE F0(X) = U(0)*X/L + U(L)*(1-X/L), '
WRITE(1,*)'F1(X) = (X/L)*(1-X/L), '
WRITE(1,*)'AND F2(X) = ((X/L)*(1-X/L))**2'
WRITE(1,*)'SQUARE OF ALPHA =',ALPHASQ
WRITE(1,*)'X-COORDINATE OF LEFT END BOUNDARY =',XST
WRITE(1,*)'X-COORDINATE OF RIGHT END BOUNDARY =',XEN
WRITE(1,*)'U(0) =',U0
WRITE(1,*)'U(L) =',UL
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=======================================================================
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, C1 )
CALL INTE ( ALPHASQ, XST, XEN, NSEG, N, SAI, W, F0, F2, C2 )
CALL INTE ( ALPHASQ, XST, XEN, NSEG, N, SAI, W, F1, F1, B1 )
CALL INTE ( ALPHASQ, XST, XEN, NSEG, N, SAI, W, F1, F2, B3 )
CALL INTE ( ALPHASQ, XST, XEN, NSEG, N, SAI, W, F2, F1, B2 )
CALL INTE ( ALPHASQ, XST, XEN, NSEG, N, SAI, W, F2, F2, 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,*) 'H(F0,F1)=', C1
WRITE (1,*) 'H(F0,F2)=', C2
WRITE (1,*) 'H(F1,F1)=', B1
WRITE (1,*) 'H(F1,F2)=', B3
WRITE (1,*) 'H(F2,F1)=', B2
WRITE (1,*) 'H(F2,F2)=', B4
C
WRITE (1,*) 'A1=',A1, ' A2=', A2
WRITE (1,*) 'U(X) = F0(X)+',A1, '*F1(X)+',A2, '*F2(X)'
C
CALL OUTPUT ( XST, XEN, A1, A2 )
CLOSE (1)
STOP 'NORMAL TERMINATION'
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.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 , 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 = F1(X)**2
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,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
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