Change B = 1 to B = 0.5 in both B.

Program

c..... 1D NUMERICAL INTEGRATION ELEMENT BY ELEMENT,.... c..... INTEGRATES g(x)^2f'(x) IN x..................... c..... FLEXIBL TO EXTEND TO MULI-DIMENSIONS............ c..... WRITTEN BY DR. SHAHROUZ ALIABADI................ c..... JANUARY 2020.................................... c.......COPY RIGHT RESERVED.............................

implicit none

c...... nn : number of nodes........................ c...... ne : number of elements..................... c...... nen : number of element nodes................ c...... nquad : integration option..................... c...... i,j,k,n: counter and index......................

integer*4 nn, ne, nen, nquad, i, j, k, n

c...... A : starting point for integrtion........... c...... B : ending point for integration............ c...... DX : delta x ................................ c...... FX : derivative of f w.r.t x at x ........... c...... GKISA : g function at kisa...................... c...... FXKISA: f' function at kisa..................... c.......JAC : Jacobian ............................... c.......INTEG : integral result......................... real*8 A, B, DX, FX, GKISA, FXKISA, JAC, INTEG

c...... G : value of function g at x ............... c...... F : value of function f at x ............... c...... X : x....................................... c...... EN : element connectivity ...................

real*8, allocatable :: G(:) real*8, allocatable :: F(:) real*8, allocatable :: X(:) integer*4, allocatable :: EN(:,:)

c...... GE : element level g ........................ c...... FE : element level f ........................ c...... XE : element level X ........................

real*8, allocatable :: GE(:) real*8, allocatable :: FE(:) real*8, allocatable :: XE(:) real*8, allocatable :: SH(:) real*8, allocatable :: SHX(:)

c...... KISA : integration point ...................... c...... W : weighting function .....................

real*8, allocatable :: KISA(:) real*8, allocatable :: W(:)

c...... cpu timing variables............................

real*8 tstart, tend

c....................................................... c.......start the program............................... c...... read integration option.........................

WRITE(6,'(A)')'What is the starting point in x (A) ?' READ (5,*) a WRITE(6,'(A)')'What is the ending point in x (B) ?' READ (5,*) b

WRITE(6,'(A)')'Number of nodes (odd number) ?' READ (5,*) nn IF(mod(nn-1,2).gt.0) then WRITE(6,'(A)')'Invalid number of nodes.' STOP ENDIF

WRITE(6,'(A)')'Number of element nodes(nen) ?' READ (5,*) nen

WRITE(6,'(A)')'Choose number of quadrature points' WRITE(6,'(A)')'----------------------------------' WRITE(6,'(A)')'1: One-Point Gaussian (Mid-Point)' WRITE(6,'(A)')'2: Two-Point Gaussian' WRITE(6,'(A)')'3: Three-Point Gaussian' READ (5,*) nquad

c...... initialize CPU timing...........................

CALL cpu_time(tstart)

c...... setup computational constants...................

IF (nen.eq.2) then ne = nn -1 ELSE IF (nen.eq.3) then ne = (nn-1)/2 ELSE WRITE(6,'(A)')'Invalid number for nen.' STOP ENDIF

DX = (B-A)/(nn-1)

c...... ALLOCATE MEMOORY

allocate (G(nn)) allocate (F(nn)) allocate (X(nn))

allocate (EN(nen,ne))

allocate ( GE(nen)) allocate ( FE(nen)) allocate ( XE(nen)) allocate ( SH(nen)) allocate (SHX(nen))

allocate (KISA(nquad)) allocate ( W(nquad))

c...... create data from the continous functions g and f DO i = 1,nn X(i) = A + (i-1)*dx CALL GF(G(i),F(i),X(i)) ENDDO

c...... form EN array................................... DO i = 1,ne EN(1,i) = (nen-1)*(i-1) + 1 EN(2,i) = EN(1,i) + 1 IF (nen.eq.3) THEN EN(3,i) = EN(2,i) + 1 ENDIF ENDDO

c...... call integration rule........................... CALL INTRULE (KISA,W,nquad)

c...... loop over the elements.......................... INTEG = 0.0 DO i = 1, ne

c...... localize the date............................... DO j = 1, nen XE(j) = X(EN(j,i)) GE(j) = G(EN(j,i)) FE(j) = F(EN(j,i)) ENDDO

c...... loop over the number of integration points...... DO k = 1, nquad

c...... evaluate shape functions and derivatives at kisa CALL SHAPEF (SH,SHX,XE,JAC,KISA(k),nen)

c...... evaluate g and f' at integration point.......... GKISA = 0.0 FXKISA = 0.0 DO j=1,nen GKISA = GKISA + GE(j)* SH(j) FXKISA = FXKISA + FE(j)*SHX(j) ENDDO

c...... evaluate the integral of g^2 f'................. INTEG = INTEG+GKISA*GKISA*FXKISA*JAC*W(k)

ENDDO !! end number of integration point loop

ENDDO !! end number of element loop

c...... output the result............................... WRITE(6,'(A,e16.9)') 'Integration is = ', INTEG

c...... calculate cputime............................... CALL cpu_time(tend) WRITE(6,'(A,e12.5)') 'Total CPU time = ',tend-tstart

END c....................................................... c...... subroutine to assign value of g and f .......... c...... input x......................................... c...... output g and f..................................

SUBROUTINE GF(g,f,x) implicit none real*8 pi parameter (pi = 3.141592653589793) real*8 g, f, x

c...... you can define your g function here ............ g = 1.0-exp(-x) f = x**3.0

return end

c...... subroutine shape function ...................... SUBROUTINE SHAPEF (SH,SHX,XE,JAC,kisa,nen) implicit none integer*4 nen real*8 SH(nen),SHX(nen), XE(nen), kisa, JAC

IF (nen.eq.2) then SH (1) = +0.5*(1-kisa) SH (2) = +0.5*(1+kisa) SHX(1) = -0.5 SHX(2) = +0.5 JAC = XE(1)*SHX(1) + XE(2)*SHX(2) SHX(1) = SHX(1)/JAC !! derivative w.r.t x SHX(2) = SHX(2)/JAC !! derivative w.r.t x ELSE SH (1) = +0.5*kisa*(kisa-1.0) SH (2) = +1.0-kisa**2.0 SH (3) = +0.5*kisa*(kisa+1.0) SHX(1) = kisa - 0.5 SHX(2) = -2.0*kisa SHX(3) = kisa + 0.5 JAC = XE(1)*SHX(1) + XE(2)*SHX(2) + XE(3)*SHX(3) SHX(1) = SHX(1)/JAC !! derivative w.r.t x SHX(2) = SHX(2)/JAC !! derivative w.r.t x SHX(3) = SHX(3)/JAC !! derivative w.r.t x ENDIF RETURN END

c...... subroutine integration rule .................... SUBROUTINE INTRULE (kisa,w,nquad) implicit none

integer*4 nquad real*8 kisa(nquad), w(nquad)

IF (nquad.EQ.1) THEN w(1) = 2.0 kisa(1) = 0.0 ELSE IF (nquad.EQ.2) THEN w(1) = 1.0 kisa(1) = -1.0/(3.0**0.5) w(2) = W(1) kisa(2) = -kisa(1) ELSE IF (nquad.EQ.3) THEN nquad = 3 w(1) = 5.0/9.0 kisa(1) = -(3.0/5.0)**0.5 w(2) = 8.0/9.0 kisa(2) = 0.0 w(3) = w(1) kisa(3) = -kisa(1) ELSE WRITE(6,'(A)')'Invalid number for nquad.' STOP ENDIF return end

(1.27) Problem 1.2. Modify the g and f functions in program "INTELEMENT.t" as following: g(x) = 0; f() = te" Compile the program and execute it with parameters below. A=0, B=1, nn=101, nen=2, nquad=1 Record your result. Execute the program one more time with following parameters. (1.28) A=0, B=1, nn = 100,001, nen =3, nquad=3 (1.29) Compare your results are write a paragraph explaining the results. Can you find out what the exact integral is? (1.27) Problem 1.2. Modify the g and f functions in program "INTELEMENT.t" as following: g(x) = 0; f() = te" Compile the program and execute it with parameters below. A=0, B=1, nn=101, nen=2, nquad=1 Record your result. Execute the program one more time with following parameters. (1.28) A=0, B=1, nn = 100,001, nen =3, nquad=3 (1.29) Compare your results are write a paragraph explaining the results. Can you find out what the exact integral is