Rcjp's Weblog

February 14, 1991

BSc Fortran Numerics Course – Solving Equations

Filed under: physics — rcjp @ 11:06 am

We did some very basic differential equation solving, simultaneous equations and curve fitting…


    CCCC  ******************************************************************
    CCCC  **        SOLUTION OF FIRST ORDER DIFFERENTIAL EQUATION         **
    CCCC  **                USING EULERS METHODS                          **
    CCCC  ****************************************************************** 

          EXTERNAL FUNC

          REAL XO,XOUT,Y,H,FUNC
          DATA XO,XOUT,Y/0.0,3.0,1.0/

          PRINT *,'ENTER STEP SIZE '
          READ *,H

          CALL EULER (XO,XOUT,H,Y,FUNC)

          PRINT *,'Y VALUE AT X=3.0 IS ',Y
          END

          REAL FUNCTION FUNC (X,Y)
          REAL X,Y
          FUNC=-Y
          END

          SUBROUTINE EULER (XO,XOUT,H,Y,FUNC)
          REAL I,XO,XOUT,H,Y
          DO 100 I=XO,XOUT,H
          Y=Y+H*FUNC (XO,Y)
      100 CONTINUE
          END

and


    CCCC  ******************************************************************
    CCCC  **        SOLUTION OF FIRST ORDER DIFFERENTIAL EQUATION         **
    CCCC  **                USING RUNGE-KUTTA METHODS                     **
    CCCC  ******************************************************************

          EXTERNAL FUNC

          REAL XO,XOUT,Y,H,FUNC
          DATA XO,XOUT,Y/0.0,3.0,1.0/

          PRINT *,'ENTER STEP SIZE '
          READ *,H

          CALL RK4 (XO,XOUT,H,Y,FUNC)

          PRINT *,'Y VALUE AT X=3.0 IS ',Y

          END


          REAL FUNCTION FUNC (X,Y)
          REAL X,Y
          FUNC=-Y
          END

          SUBROUTINE RK4 (XO,XOUT,H,Y,FUNC)
          REAL X,XO,XOUT,H,Y,C1,C2,C3,C4
          DO 100 X=XO,XOUT,H
          C1=H*FUNC (X,Y)
          C2=H*FUNC (X+0.5*H,Y+0.5*C1)
          C3=H*FUNC (X+0.5*H,Y+0.5*C2)
          C4=H*FUNC (X+H,Y+C3)
          Y=Y+(C1+2*C2+2*C3+C4)/6
      100 CONTINUE
          END

Simultaneous Equations


    CCCC  ******************************************************************
    CCCC  **            SOLUTION OF SIMULTANEOUS EQUATIONS                **
    CCCC  **              BY GAUSSIAN ELIMINATION METHOD                  **
    CCCC  **            WITH PARTIAL PIVOTAL CONDENSATION                 **
    CCCC  **                AND CHECK FOR SINGULARITY                     **
    CCCC  ******************************************************************

    CCCC  Initialise variables                                                                                                           
          INTEGER I,J,K,N,L
          REAL A,B,X,M,F,TOL,C,D,S
          DIMENSION A(20,20),X(20),B(20)
          COMMON A,B,X

    CCCC  Enter equations                                                                                                                

          PRINT*,'ENTER NO. OF EQUATIONS AND TOLERANCE '
          READ*,N,TOL
          DO 10 I=1,N
            PRINT114,'ENTER COEFFICIENTS OF EQU.',I
            READ*,(A(I,J),J=1,N),B(I)
       10 CONTINUE
          PRINT113,'EQUATIONS ENTERED'
          CALL PR(N)

    CCCC  PARTIAL PIVOTAL CONDENSATION                                                                                                   

          DO 5 J=1,N
            L=J
            DO 2 I=1,N
              IF (A(I,J).GT.A(L,J)) L=I
        2   CONTINUE
            DO 3 K=1,N
              C=A(J,K)
              A(J,K)=A(L,K)
              A(L,K)=C
        3   CONTINUE
          D=B(J)
          B(J)=B(L)
          B(L)=D
        5 CONTINUE
          PRINT113,'PARTIAL PIVOTAL CONDENSATION'
          CALL PR(N)

    CCCC  Gaussian elimination method

          DO 50 K=1,N-1
            DO 40 I=K+1,N
              IF (ABS(A(K,K)).LT.TOL) STOP
              M=A(I,K)/A(K,K)
              S=0
              DO 30 J=K,N
                A(I,J)=A(I,J)-M*A(K,J)
                S=S+A(I,J)**2
       30     CONTINUE
            B(I)=B(I)-M*B(K)
            IF (S.LT.TOL) THEN
              PRINT113,'SINGULARITY ERROR'
              STOP
            ENDIF
       40   CONTINUE
       50 CONTINUE
          PRINT113,'TRIANGULAR FORM BEFORE BACK SUBSTITUTION'
          CALL PR(N)

    CCCC  Back substitution & output of results                                                                                          

          PRINT113,'SOLUTIONS'
          X(N)=B(N)/A(N,N)
          PRINT115,N,X(N)
          DO 110 J=N-1,1,-1
            DO 100 L=J+1,N
              F=F+A(J,L)*X(L)
      100   CONTINUE
            X(J)=1/A(J,J)*(B(J)-F)
            PRINT115,J,X(J)
            F=0
      110 CONTINUE
      113 FORMAT (/,1X,A)
      114 FORMAT (1X,A,I1,6X)
      115 FORMAT (1X,'X(',I1,') = ',F9.5)
          END

    CCCC  Display equations                                                                                                              

          SUBROUTINE PR(N)
          INTEGER I,J,N
          COMMON A(20,20),B(20),X(20)
          DO 210 I=1,N
            DO 200 J=1,N
              PRINT215,A(I,J)
      200   CONTINUE
          PRINT216,B(I)
      210 CONTINUE
      215 FORMAT (1X,F9.4,4X,\)
      216 FORMAT (1X,'=',2X,F9.4)
          END

We were obviously using the Salford FTN77 compiler since there is a log file


    SALFORD UNIVERSITY FTN77/386 VER. 1.67      A:\GAUSSPC.FOR    15:29:42    Tuesday, 22 January  1991
    COMPILER OPTIONS:  LISTING ANSI INTL NODCLVAR NOMAP NOCHECK LOGS DYNM OFFSET LGO NOPAGETHROW NOSILENT
                       NO_OPTIMISE NOWEITEK

I can remember having to go to the machine room in the Newton building to get the listings from the old 132 column printer (everywhere else was just epson 850’s which couldn’t print that wide.)

Finite Difference


    CCCC  ******************************************************************
    CCCC  **   CALCULATES THE TEMPERATURE IN A BLOCK OF HEATED MATERIAL   **
    CCCC  **   USES INITIAL GUESSES AND RELAXATION TECHNIQUE APPLIED TO   **
    CCCC  **   THE FINITE DIFFERENCE EQUATIONS AND DRAWS A GRAPH OF THE   **
    CCCC  **   NUMBER OF ITERATIONS AGAINST THE RELAXATION COEFF. ALPHA   **
    CCCC  ******************************************************************

    CCCC  Block is heated ...                                              
    CCCC      500 500 500 500                                             
    CCCC      100 ??? ??? 100                                            
    CCCC      100 ??? ??? 100                                           
    CCCC      100 100 100 100                                          

          REAL T,RES,RESMX,DECT,W,WMAX,WMIN,R,TT
          INTEGER I,J,M,N,TMMIN,TMMAX,TNMIN,TNMAX,IT,PTS,NUM
          DIMENSION T(10,10),TT(10,10),RES(10,10),PTS(50),W(50)

          DATA  ((TT(I,J),I=1,4),J=1,4),TMMIN,TMMAX,TNMIN,TNMAX/
         1  4*500.0,100.0,2*300.0,2*100.0,2*200.0,5*100.0,1,4,1,4/

          PRINT *,'ENTER RELAXATION COEFFICIENT BOUNDARIES '
          READ *,WMIN,WMAX

          DO 250 R=WMIN,WMAX,(WMAX-WMIN)/40.0

          DO 5 I=1,4
          DO 4 J=1,4
          T(I,J)=TT(I,J)
        4 CONTINUE
        5 CONTINUE

          DO 200 IT=0,200

    CCCC  FIND RESIDUALS                                                                                                                 

            RESMX=0.0
            DO 20 M=TMMIN+1,TMMAX-1
              DO 10 N=TNMIN+1,TNMAX-1
                RES (M,N)=T(M+1,N)+T(M-1,N)+T(M,N+1)+T(M,N-1)-4.0*T(M,N)
                IF (ABS(RES(M,N)).GT.RESMX) RESMX=RES(M,N)
       10     CONTINUE
       20   CONTINUE

    CCCC  ALTER TEMPERATURE OF ALL POINTS BY RESMX/4.0*R

            DECT=RESMX/4.0*R
            DO 40 M=TMMIN+1,TMMAX-1
              DO 30 N=TNMIN+1,TNMAX-1
                T(M,N)=T(M,N)+DECT
       30     CONTINUE
       40   CONTINUE
            IF (ABS(RESMX).LT.0.1 .AND. IT.GT.0) GOTO 230
      200 CONTINUE
    C  230 PRINT *,'REL COEFF',R,'  STEPS ',IT                                                                                           
      230 NUM=NUM+1
          W(NUM)=R
          PTS(NUM)=IT
      250 CONTINUE
      300 FORMAT (1X,I3,1X,12(F6.1))

          CALL SCRSET
          CALL PPOINT (NUM,W,PTS)

          END


    CCCC  Initialise screen and clear it.                                                                                                

          SUBROUTINE SCRSET
          INCLUDE 'F77PC'
          CALL INITSCREEN
          CALL SETSCREENMODE (EGAGRAPH)
          CALL CLRVIDEO
          END


          SUBROUTINE PPOINT (NPTS,X,Y)
          INCLUDE 'F77PC'
          INTEGER NPTS,Y,N,I,XP,YP
          REAL X,XMAX,YMAX,SCALEX,SCALEY,XMIN,YMIN
          DIMENSION Y(50),X(50)
    CCCC  Draw axis                                                                                                                      

          CALL DRAW (20,330,620,330,WHITE)
          CALL DRAW (20,330,20,10,WHITE)

    CCCC  Mark points                 
          XMIN=1000
          YMIN=1000
          DO 390 N=1,NPTS
            IF (X(N).GT.XMAX) XMAX=X(N)
            IF (X(N).LT.XMIN) XMIN=X(N)
            IF (Y(N).GT.YMAX) YMAX=Y(N)
            IF (Y(N).LT.YMIN) YMIN=Y(N)
      390 CONTINUE

          SCALEX=600.0/(XMAX-XMIN)
          SCALEY=320.0/(YMAX-YMIN)
          DO 400,I=1,NPTS
          XP=INT((X(I)-XMIN)*SCALEX+20)
          YP=INT(330-(Y(I)-YMIN)*SCALEY)-2
          CALL DRAW (XP-2,YP,XP+2,YP,WHITE)
          CALL DRAW (XP,YP-2,XP,YP+2,WHITE)
      400 CONTINUE
          END

Least Squares Fitting


    CCCC  ******************************************************************
    CCCC  **      DRAWS A GRAPH OF DATA HELD IN FILE 'GDATA.DAT'          **
    CCCC  **               AND DRAW LEAST SQUARES FIT                     **
    CCCC  ******************************************************************

    CCCC  Data file includes the number of points                                                                                        
    CCCC  and X,Y coordinates                                                                                                            

          INCLUDE 'F77PC'
          INTEGER NPTS,I
          REAL X,Y,GRAD,INTCPT
          DIMENSION X(30),Y(30)

          OPEN (5,FILE='A:\LSQ\GDATA.DAT',STATUS='OLD')
          READ (5,*) NPTS
          DO 10 I=1,NPTS
            READ (5,*) X(I),Y(I)
       10 CONTINUE

          CALL LSFIT (X,Y,NPTS,GRAD,INTCPT)
          CALL SCRSET
          CALL GRAFIT (X,Y,NPTS,GRAD,INTCPT)
          PRINT *,'GRADIENT=',GRAD
          PRINT *,'INTERCEPT=',INTCPT
          CLOSE (5)
          END


    CCCC  Calculation of least squares fit to data                         
    CCCC  with coordinates held in array X,Y.                             
    CCCC  Returns GRAD-Gradient and INTCPT the intercept                 
    CCCC  on the y -axis                                                

          SUBROUTINE LSFIT (X,Y,NPTS,GRAD,INTCPT)
          INTEGER NPTS,I
          REAL X,Y,SX,SXX,SXY,SY,DEL,GRAD,INTCPT
          DIMENSION X(*),Y(*)
          DO 250 I=1,NPTS
            SX=SX+X(I)
            SXX=SXX+X(I)*X(I)
            SXY=SXY+X(I)*Y(I)
            SY=SY+Y(I)
      250 CONTINUE
          DEL=NPTS*SXX-(SX*SX)
          GRAD=(NPTS*SXY-SX*SY)/DEL
          INTCPT=(SXX*SY-SXY*SX)/DEL
          END


    CCCC  Initialise screen and clear it.

          SUBROUTINE SCRSET
          INCLUDE 'F77PC'
          CALL INITSCREEN
          CALL SETSCREENMODE (EGAGRAPH)
          CALL CLRVIDEO
          END


    CCCC  Plots the data, draws the axis                               
    CCCC  and the line of best fit                                    

          SUBROUTINE GRAFIT (X,Y,NPTS,GRAD,INTCPT)
          INCLUDE 'F77PC'
          INTEGER I,NPTS,XP,YP,Y1,Y2
          REAL X,Y,SCALEX,SCALEY,GRAD,INTCPT
          DIMENSION X(*),Y(*)

    CCCC  Draw axis                                                                                                                      
          CALL DRAW (20,330,620,330,WHITE)
          CALL DRAW (20,330,20,10,WHITE)

    CCCC  Mark points                                                                                                                    
          SCALEX=600/X(NPTS)
          SCALEY=320/Y(NPTS)
          DO 300,I=1,NPTS
          XP=INT(X(I)*SCALEX+20)
          YP=INT(330-Y(I)*SCALEY)
          CALL DRAW (XP-2,YP,XP+2,YP,WHITE)
          CALL DRAW (XP,YP-2,XP,YP+2,WHITE)
      300 CONTINUE

    CCCC  Draw best fit line                                                                                                             
          Y1=INT(330-INTCPT*SCALEY)
          Y2=INT(330-(GRAD*X(NPTS)+INTCPT)*SCALEY)
          CALL DRAW (20,Y1,620,Y2,WHITE)
          END
Advertisements

Leave a Comment »

No comments yet.

RSS feed for comments on this post. TrackBack URI

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

Create a free website or blog at WordPress.com.

%d bloggers like this: