Code

 C SOLUTION OF ONE DIMENSIONAL LINEAR EQUATIONS

      DIMENSION X(10),W(5),GP(5),RKE(10,4,4),FE(10,4),IECON(10,4),

     1 RK(10,10),F(10),EBC(4),IESS(4),PHI(10)

      COMMON/NODES/NE,NODE,NT

COMMON/APP/ITA

PRINT*,'GIVE THE TYPE OF APPROXIMATION REQUIRED'

PRINT*,'GIVE 1 FOR LINEAR, 2 FOR QUADRATIC'

READ*,ITA

PRINT*,'ENTER THE NUMBER OF ELEMENTS'

READ*,NE

PRINT*,'ENTER THE GRID POINTS'

READ*,(X(J),J=1,NE+1)

PRINT*,'ENTER THE TOTAL NUMBER OF NODES IN THE ENTIRE REGION'

READ*,NT

PRINT*,'ENTER THE NUMBER OF NODES FOR EACH ELEMENT'

READ*,NODE

DO 1 I=1,NE

PRINT*,'ENTER THE ELEMENT CONNECTIVITY OF ELEMENT NUMBER:',I

READ*,(IECON(I,J),J=1,NODE)

1 CONTINUE

DO 2 I=1,NT

DO 2 J=1,NT

2 RK(I,J)=0.

DO 31 I=1,NT

31 F(I)=0.

W(1)=5./9.

W(2)=8./9.

W(3)=5./9.

GP(1)=-SQRT(3./5.)

GP(2)=0

GP(3)=SQRT(3./5.)

DO 3 IE=1,NE

DO 4 I=1,NODE

DO 5 J=1,NODE

RKE(IE,I,J)=0.

5 CONTINUE

      FE(IE,I)=0

4 CONTINUE

3 CONTINUE

      DO 6 IE=1,NE

DO 7 I=1,NODE

DO 8 J=1,NODE

DO 9 L=1,3

PT=(X(IE+1)-X(IE))/2.*GP(L)+(X(IE+1)+X(IE))/2.

RKE(IE,I,J)=RKE(IE,I,J)+W(L)*(X(IE+1)-X(IE))/2.*

1(-DA(PT)*RN(IE,I,PT,X)*DRN(IE,J,PT,X)-A(PT)*DRN(IE,I,PT,X)

     2*DRN(IE,J,PT,X)+B(PT)*RN(IE,I,PT,X)*DRN(IE,J,PT,X)

     3+C(PT)*RN(IE,I,PT,X))

9 CONTINUE

8 CONTINUE

DO 10 L=1,3

PT=(X(IE+1)-X(IE))/2.*GP(L)+(X(IE+1)+X(IE))/2.

FE(IE,I)=FE(IE,I)+W(L)*(X(IE+1)-X(IE))/2.*D(PT)*RN(IE,I,PT,X)

10 CONTINUE

7 CONTINUE

6 CONTINUE

DO 11 IE=1,NE

PRINT*,'THE LOCAL STIFFNESS MATRIX OF ELEMENT NUMBER:',IE

DO 12 I=1,NODE

12 PRINT*,(RKE(IE,I,J),J=1,NODE)

PRINT*,'THE LOCAL FORCE VECTOR OF ELEMENT NUMBER:',IE

PRINT*,(FE(IE,I),I=1,NODE)

11 CONTINUE

CALL ASEMBL(RKE,FE,IECON,RK,F)

PRINT*,'THE GLOBAL STIFFNESS MATRIX IS:'

DO 15 I=1,NT

PRINT*,(RK(I,J),J=1,NT)

15 CONTINUE

PRINT*,'THE GLOBAL FORCE VECTOR IS:'

PRINT*,(F(I),I=1,NT)

PRINT*,'GIVE THE NUMBER OF ESSENTIAL BOUNDARY CONDITIONS'

READ*,M

PRINT*,'GIVE NODAL NUMBERS OF ESSENTIAL BOUNDARY CONDITIONS'

READ*,(IESS(I),I=1,M)

PRINT*,'GIVE THE VALUES OF ESSENTIAL BOUNDARY CONDITIONS'

READ*,(EBC(I),I=1,M)

DO 16 I=1,M

F(IESS(I))=EBC(I)

DO 17 J=1,NT

IF(IESS(I).EQ.J) THEN

    RK(J,J)=1.

  ELSE

    RK(IESS(I),J)=0.

      ENDIF

17 CONTINUE

16 CONTINUE

CALL GAUELM(RK,F,PHI,NT)

PRINT*,'THE VALUES OF THE NODAL PARAMETERS ARE'

PRINT*,(PHI(I),I=1,NT)

20 PRINT*,'GIVE THE POINT AT WHICH WE WANT TO FIND THE SOLUTION'

READ*,PT

DO 22 I=1,NE

IF (X(I).LE.PT .AND. PT.LE.X(I+1)) THEN

    J=1

    GOTO 25

ENDIF

22 CONTINUE

25 SOLN=0.

DO 26 I=1,NODE

26 SOLN=SOLN+PHI(IECON(J,I))*RN(J,I,PT,X)

PRINT*,'DO YOU WANT TO SOLUTION AT ANOTHER POINT?'

PRINT*,'IF YES GIVE 1, OTHERWISE GIVE ANY NUMBER OTHER THAN 1'

READ*,J

IF(J.EQ.1) GOTO 20

STOP

END

C -----------------------------------------------------------------------

      FUNCTION A(PT)

      COMMON/APP/ITA

IF(ITA.EQ.1) A=-(1.+PT)

IF(ITA.EQ.2) A=-1.

R ETURN

END

C-----------------------------------------------------------------------------

      FUNCTION DA(PT)

COMMON/APP/ITA

IF(ITA.EQ.1) DA=-1.

IF(ITA.EQ.2) DA=0.

RETURN

END

C-------------------------------------------------------------------------------

      FUNCTION B(PT)

COMMON/APP/ITA

IF(ITA.EQ.1) B=-1.

IF(ITA.EQ.2) B=0.

RETURN

END

C--------------------------------------------------------------------------------

      FUNCTION C(PT)

COMMON/APP/ITA

IF(ITA.EQ.1) C=0. 

IF(ITA.EQ.2) C=0.

RETURN

END

C-------------------------------------------------------------------------------

      FUNCTION D(PT)

COMMON/APP/ITA

IF(ITA.EQ.1) D=1.+4.*PT

IF(ITA.EQ.2) D=PT

RETURN

END

C-------------------------------------------------------------------------------

      FUNCTION RN(IE,I,PT,X)

DIMENSION X(10)

COMMON/APP/ITA

IF(ITA.EQ.1) THEN 

IF (I.EQ.1) RN=(X(IE+1)-PT)/(X(IE+1)-X(IE))

IF (I.EQ.2) RN=(PT-X(IE))/(X(IE+1)-X(IE))

ENDIF

IF(ITA.EQ.2) THEN

XM=(X(IE+1)+X(IE))/2.

IF(I.EQ.1) RN=(PT-XM)*(PT-X(IE+1))/((X(IE)-XM)*(X(IE)-X(IE+1)))

IF(I.EQ.2) RN=(PT-X(IE))*(PT-X(IE+1))/((XM-X(IE))*(XM-X(IE+1)))

IF(I.EQ.3) RN=(PT-X(IE))*(PT-XM)/((X(IE+1)-X(IE))*(X(IE+1)-XM))

ENDIF

RETURN

END

C------------------------------------------------------------------------------------

FUNCTION DRN(IE,I,PT,X)

DIMENSION X(10)

COMMON/APP/ITA

IF(ITA.EQ.1) THEN

IF(I.EQ.1) DRN=-1./(X(IE+1)-X(IE))

IF(I.EQ.2) DRN=1./(X(IE+1)-X(IE))

ENDIF

IF (ITA.EQ.2) THEN

XM=(X(IE+1)+X(IE))/2.

IF(I.EQ.1) DRN=(2.*PT-(XM+X(IE+1)))/((X(IE)-XM)*(X(IE)-X(IE+1)))

IF(I.EQ.2) DRN=(2.*PT-(X(IE)+X(IE+1)))/((XM-X(IE))*(XM-X(IE+1)))

      IF(I.EQ.3) DRN=(2.*PT-(XM+X(IE)))/((X(IE+1)-X(IE))*(X(IE+1)-XM))

ENDIF

RETURN

END

C------------------------------------------------------------------------------------------

      SUBROUTINE ASEMBL(RKE,FE,IECON,RK,F)

DIMENSION RKE(10,4,4),FE(10,4),IECON(10,4),RK(10,10),F(10)

INTEGER GRN,GCN

COMMON/NODES/NE,NODE,NT

DO 1 I=1,NT

DO 2 J=1,NT

RK(I,J)=0.

2 CONTINUE

F(I)=0

1 CONTINUE

      DO 14 I=1,NE

DO 15 J=1,NODE

GRN=IECON(I,J)

DO 16 L=1,NODE

GCN=IECON(I,L)

RK(GRN,GCN)=RK(GRN,GCN)+RKE(I,J,L)

16 CONTINUE

F(GRN)=F(GRN)+FE(I,J)

15 CONTINUE

14 CONTINUE

      RETURN

END

C-----------------------------------------------------------------------------------------

      SUBROUTINE GAUELM(A,RHS,X,N)

DIMENSION A(10,10),RHS(10),X(10)

C REDUCING THE COEFFICIENT MATRIX TO UPPER TRIANGULER FORM

      DO 20 I=1,N-1

IF (A(I,I).NE.0.) GOTO 100

DO 30 K=I+1,N

IF(A(K,I).NE.0.)GOTO 40

30 CONTINUE

      PRINT*,'THE GIVEN COEFFICIENT MATRIX IS SINGULER'

STOP

40 DO 50 J=1,N

      TEMP=A(I,J)    

A(I,J)=A(K,J)

A(K,J)=TEMP

50    CONTINUE

      TEMP=RHS(I)

RHS(I)=RHS(K)

RHS(K)=TEMP

100   DO 60 K=I+1,N

      R=-A(K,I)/A(I,I)

DO 70 J=I+1,N

70    A(K,J)=A(K,J)+R*A(I,J)

      RHS(K)=RHS(K)+R*RHS(I)

60    CONTINUE

20    CONTINUE

      IF(A(N,N).NE.0.) GOTO 120

PRINT*,'THE GIVEN COEFFICIENT MATRIX IS SINGULER'

'

C     BACK SUBSTITUTION

120   X(N)=RHS(N)/A(N,N)

      DO 130 I=N-1,1,-1

SUM=0.

DO 140 J=I+1,N

140   SUM=SUM+A(I,J)*X(J)

      X(I)=(RHS(I)-SUM)/A(I,I)

130   CONTINUE

      RETURN

END