PROGRAM WRM3X3S1DN
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 = 2
C DIRICHLT -- NEUMANN
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=2 )
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
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='WRM2X2S1DN.FEM',STATUS='UNKNOWN' )
WRITE(1,*)'==== ONE DIMENSIONAL HELMHOLTZ EQUATION DOF=2 ===='
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)+A2*F2(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'
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
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) )
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) )
C
A(2,1) = A(1,2)
CALL INTE ( ALPHASQ, XST, XEN, NSEG, N, SAI, W, F2, F2, A(2,2) )
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,*) 'U(X)=F0(X)+',C(1),'*F1(X)+',C(2),'*F2(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 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
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)
DUDX = DERIV(F0,X)+C(1)*DERIV(F1,X)+C(2)*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, 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