PROGRAM WRM2X2
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)=0. & U(1)=1.
C -------------- VARIABLE DEFNITION ----------- 11/28/2023 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 /POWER/ POWER1, POWER2
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
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 = 1.D0
UL = 1.D0
POWER1 = 0.D0
POWER2 = 0.D0
C=======================================================================
C---------- K = POWER1, M=POWER2 -------------
DO WHILE ( POWER1+POWER2 .LE. 2.D0 )
WRITE(*,*) 'NOTE THAT K+M MUST BE GREATER THEN 2.'
WRITE(*,fmt='(a)', advance='no') 'TYPE IN K='
READ(*,*) POWER1
WRITE(*,fmt='(a)', advance='no') 'TYPE IN M='
READ(*,*) POWER2
IF ( POWER1+POWER2 .GT. 2.D0 ) THEN
WRITE(*,*) 'K AND M ARE WITHIN ACCEPTABLE RANGES.'
END IF
END DO
C=======================================================================
OPEN ( 1, FILE='WRM2X2.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)**K*(1-X/L)**M'
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
WRITE(1,*)'POWER1 =',POWER1
WRITE(1,*)'POWER2 =',POWER2
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
COMMON /POWER/ POWER1, POWER2
F2 = (X/RL)**POWER1 * ( 1.D0- X/RL )**POWER2
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.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